Poisson Equation
Toppic
The Poisson equation is the simplest example of the PDE’s considerd in Solid Mechanics. It is an eliptical PDE, and is simplified compared to linear elasticity in the sense that its solution is a scalar field, instead fo the vector field found in elasticity problems. This makes Poisson’s equation a good start to explore numerical solving strategies for Solid Mechanics problems.
Bram Lagerweij
11 Feb 2020
Table of Contents
1 Laplace Equation
The most basic description of the Laplace equation is given by:
Where the entirety of the boundary \(\partial\Omega\) is the union of these to boundary conditions that do not intersect.
The following images summarizes this.
2 Poisson equation
In case of non-homogeneous formulations the Laplace equations is called the Poisson equation.
The boundary condition can still be defined in the same way as in the Laplace equation. An example of such a Poisson problem in 1D is a statically determinate Euler-Bernoulli beam problem. Solving a these linear beam problem can be done with finite differences.
The PDE described by
Where \(M\) is the internal bending moment of the beam. This beam has a length \(L\) and a stiffness EI. In general these kinds of problems can not be solved directly in this way, as it is not always possible to describe the moment explicitly, but because our cantilever beam is statically determinate it can be done. Now we’ll be exploring two examples to introduce the different types of boundary conditions.
Example 1: Dirichlet
In this example we consider a beam with a length of 1000mm which is simply supported at \(x=0\) and \(x=250\). Simply supported means that the displacement \(u\) at those points is fixed and equals 0. That is our ODE becomes:
where I did compute the moment equation explicitly already. To derive \(u''\) a central difference scheme is used,
We’ll be evaluating this derivative an \(N\) regularly distributed points in our domain. And if we note \(x_n\) as the location of one of these points than we can note the derivative as:
This is implemented into a matrix format by finitedifference.Dxx()
, such that:
where \(u\) is a vector with the field at all the discrete points and \(u''\) the derivative that was calculated. This does however not yet specify the way to analyze the derivative at the first and last points. After all that would require the calculation of \(u\) outside the domain. As a result the matrix will have an empty first and last row.
This and the right hand side (\(f\)) of the Poisson equation are available through finitedifference.poisson()
.
You would expect that we can solve the system of equations:
but that is not true, as we’ll have to deal with the boundary conditions as well, without those the problem is singular. To be specific we know that \(u(0)=0\) and \(u(L/4)=0\), this can be used to make the problem determinate. Lets say that \(x_0 = 0\) and \(x_n = L/4\) then we can add the following to equations to our system of equations:
these two equations can be placed in the still empty first and last row of our stiffness matrix and right hand side. That is in the first row we make the first element equal to 1 and the rest all equal to 0. Similarly the right hand side of the first degree of freedom is set to 0. In the last row we set the degree of freedom that corresponds to \(x_n\) to 1 and the rest to 0, here we do also set the right hand side of the last row equal to zero (see lines 53 to 61 in the code below).
1# Importing External modules
2import sys
3import numpy as np
4from scipy.sparse.linalg import spsolve
5
6# Importing my own scripts
7sys.path.insert(1, '../src')
8from finitedifference import poisson
9
10
11def moment(P, L):
12 """
13 Moment as a function of :math:`x` of the double simply supported beam.
14
15 Parameters
16 ----------
17 P : float
18 Applied load.
19 L : float
20 Length of the beam
21
22 Returns
23 -------
24 callable
25 The moment :math:`M(x)` of the beam.
26 """
27
28 def fun(x):
29 shape = np.shape(x)
30
31 if len(shape) == 0:
32 if x < L / 4:
33 m = -3 * P * x
34 else:
35 m = P * (x - L)
36 else:
37 m = np.zeros_like(x)
38 ind = np.where(x < L / 4) # where x < L/4
39 m[ind] = -3 * P * x[ind]
40
41 ind = np.where(L / 4 <= x) # where L/4 < x
42 m[ind] = P * (x[ind] - L)
43 return m
44
45 return fun
46
47
48if __name__ == '__main__':
49 # Define properties of the problem.
50 L = 1000 # 1000 mm length
51 P = 1 # 1 N load
52 EI = 187500000 # Beam bending stiffness Nmm^4
53
54 # Discretion of the space.
55 dof = 101 # Number of nodes
56 x, dx = np.linspace(0, L, dof, retstep=True)
57
58 # Calculate the internal Moment.
59 M = moment(P, L) # Create a callable for the moment in Nmm
60
61 # Create linear problem.
62 K, f = poisson(dof, dx, M, c=EI)
63
64 # Boundary condition u(0) = 0
65 K[0, 0] = 1
66 f[0] = 0
67
68 # Boundary condition u(L/4) = 0 For this purpose we use
69 # the last row of the matrix, this row is not yet used.
70 index = int(dof / 4)
71 K[-1, index] = 1
72 f[-1] = 0
73
74 # Solve the problem.
75 u = spsolve(K, f)
Example 2: Dirichlet and Neumann
The approach follows exactly what was described in example 1, except of course the constraints. Our problem is formulate following:
where the moment did change as well because the loading conditions changed. That is after discritization our system of equations is represented by:
Now as for the boundary conditions, for the first row we again fill the first element with a 1 and leave the rest 0. In the right hand side we set the value of the forcing term equal to zero. As a result the first row reads:
Now for the Neumann boundary it is a bit more tricky. The derivative \(u'(x_0)\) can be approximated with a backwards finite difference:
I’ll put this in the last row as that one is not yet populated. That means that we have to populate the first element of the last row with a -1, the second element of that row with a 1 and set the last element of the right hand side to zero as well. (see lines 64 to 72 below)
1# Importing External modules
2import sys
3import numpy as np
4from scipy.sparse.linalg import spsolve
5
6# Importing my own scripts
7sys.path.insert(1, '../src')
8from finitedifference import poisson
9
10
11def moment(P, L):
12 """
13 Moment as a function of :math:`x` of the cantilever beam.
14
15 Parameters
16 ----------
17 P : float
18 Applied load.
19 L : float
20 Length of the beam
21
22 Returns
23 -------
24 callable
25 The moment :math:`M(x)` of the beam.
26 """
27 def fun(x):
28 return P*(x-L)
29 return fun
30
31
32if __name__ == '__main__':
33 # Define properties of the problem.
34 L = 1000 # 1000 mm length
35 P = 1 # 1 N load
36 EI = 187500000 # Beam bending stiffness Nmm^4
37
38 # Discretion of the space.
39 dof = 101 # Number of nodes
40 x, dx = np.linspace(0, L, dof, retstep=True)
41
42 # Calculate the internal Moment.
43 M = moment(P, L) # Create a callable for the moment in Nmm
44
45 # Create linear problem.
46 K, f = poisson(dof, dx, M, c=EI)
47
48 # Boundary condition u(0) = 0
49 K[0, 0] = 1
50 f[0] = 0
51
52 # Boundary condition u'(0) = 0 with a finite difference.
53 # For this purpose we use the last row of the matrix
54 # This row is not yet used
55 K[-1, 0] = -1 / dx
56 K[-1, 1] = 1 / dx
57 f[-1] = 0
58
59 # Solve the problem.
60 u = spsolve(K, f)