\[\require{physics} \renewcommand{\vec}[1]{\underline{#1}}\]

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

1 Laplace Equation

The most basic description of the Laplace equation is given by:

\[\begin{split}\grad^2 u(\vec{m}) &= \pdv[2]{u}{x} + \pdv[2]{u}{y} = 0 \qquad& \forall \vec{m} \in \Omega \\ \quad\text{s.t.:}& \quad u(\vec{m}) = \vec{\tilde{u}}(\vec{m}) & \forall \vec{m} \in \mathcal{S}_u\\ & \quad \grad {u}(\vec{m}) = \tilde{\vec{t}}(\vec{m}) & \forall \vec{m} \in \mathcal{S}_t\end{split}\]

Where the entirety of the boundary \(\partial\Omega\) is the union of these to boundary conditions that do not intersect.

\[\begin{split}\partial\Omega = \mathcal{S}_u \cup \mathcal{S}_t \\ 0 = \mathcal{S}_u \cap \mathcal{S}_t\end{split}\]

The following images summarizes this.

_images/Domain.svg

Fig. 39 A domain \(\Omega\) subjected to the Laplace equation with combined boundary conditions.

2 Poisson equation

In case of non-homogeneous formulations the Laplace equations is called the Poisson equation.

\[\begin{split}\grad^2 u(\vec{m}) &= \pdv[2]{u}{x} + \pdv[2]{u}{y} = \vec{b}(\vec{m}) \qquad& \forall \vec{m} \in \Omega \\ \quad\text{s.t.:}& \quad u(\vec{m}) = \vec{\tilde{u}}(\vec{m}) & \forall \vec{m} \in \mathcal{S}_u\\ & \quad \grad {u}(\vec{m}) = \tilde{\vec{t}}(\vec{m}) & \forall \vec{m} \in \mathcal{S}_t\end{split}\]

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

\[\begin{split}EI u''(x) = M(x) & \qquad \forall x \in\Omega = [0, L]\\\end{split}\]

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

_images/Simple_Beam_Drawing.svg

Fig. 40 A beam that is simply supported at \(x=0\) and \(250\) mm and subjected to a point load.

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:

\[\begin{split}EI \, u''(x) = & M(x) \quad \forall \quad 0 \leq x \leq 1000\\ \text{where:} \quad & M(x) = \begin{cases} -3Px & 0 \leq x \leq L/4\\ P(x - L) & L/4 \leq x \leq L \end{cases}\\ \text{s.t.:} \quad &u(0) = 0 \\ & u(L/4) = 0\end{split}\]

where I did compute the moment equation explicitly already. To derive \(u''\) a central difference scheme is used,

\[u''(x) = \frac{u(x-dx) - 2 u(x) + u(x+dx)}{dx^2}\]

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:

\[u''(x_n) = \frac{u(x_{n-1}) - 2 u(x_n) + u(x_{n+1})}{dx^2}\]

This is implemented into a matrix format by finitedifference.Dxx(), such that:

\[u'' = D_{xx} u\]

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:

\[EI\,D_{xx}\, u = f\]

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:

\[u_0 = 0 \qq{and} u_n = 0\]

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).

_images/Simply_Solution.svg

Fig. 41 The finite difference solution of the beam problem seems to be in good agreement with the exact result. This simulation was run with 101 degrees of freedom.

 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

_images/Cantilever_Drawing.svg

Fig. 42 A cantilever beam is fixed in the wall of the left and subjected to a point load at the right. This type of constraint, called an endcast, limits both the displacement and rotation, that is \(u(0)=0\) and \(u'(0)=0\).

The approach follows exactly what was described in example 1, except of course the constraints. Our problem is formulate following:

\[\begin{split}EI \, u''(x) = M(x) & \quad \forall \quad 0 \leq x \leq 1000\\ \text{where:} \quad & M(x) = P(x-L)\\ \text{s.t.:} \quad &u(0) = 0 \\ & u'(0) = 0\end{split}\]

where the moment did change as well because the loading conditions changed. That is after discritization our system of equations is represented by:

\[EI\,D_{xx}\, u = f\]

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:

\[u(x_0) = 0\]

Now for the Neumann boundary it is a bit more tricky. The derivative \(u'(x_0)\) can be approximated with a backwards finite difference:

\[u'(x_0) = \frac{-u(x_0) + u(x_1}{dx} =0\]

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)

_images/Cantilever_Solution.svg

Fig. 43 The finite difference solution of the beam problem seems to be in good agreement with the exact result. This simulation was run with 101 degrees of freedom.

 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)