Introduction Finite Element Method (FEM)

We will again use the 1-D steady state heat diffusion equation as an example and solve it using FEM. Formally speaking the steady-state diffusion equation is an elliptic PDE. Elliptic equations describe a large variety of steady state problems in earth sciences including diffusion, plate flexure, and incompressible flows.

Types of PDEs

Besides elliptic PDEs, there are, for example, also parabolic PDEs, which include the transient term, and hyperbolic PDEs, which describe e.g. wave propagation. We will learn about those a bit later in the course.

Governing equation - strong form

Let’s go on and solve the elliptic 1-D heat diffusion equation with the Finite Element Method (FEM):

(90)\[\frac{\partial}{\partial x}k\frac{\partial T_{ex}}{\partial x}=0.\]

This is called the strong form of the governing equation and the function \(T_{ex}(x)\) by a (simpler) approximation:

Approximate solution

(91)\[T_{ex} \cong \tilde{T}= \sum_{j=1}^{n} N_j T_j = N_j T_j.\]

In FEM we divide our computational domain into simple geometeric entities, the finite elements, and those finite elements have grid points. For example, in 2-D we can divide the domain in triangles with three nodes and in 1-D in simple lines with two nodes each. Think of our simple FD descritization, if we were to use the same descritization in FEM, the intervals between grid points would be our finite elements.

../_images/shapeFunction_Linear_1D.svg

Fig. 4 Example of 1D linear finite element shape functions

These linear independent functions \(N\), called shape functions in FEM, are associated with the chosen discretization (Fig. 4). Each node is associated with one shape function and that shape function has the value \(1\) at that node and is zero at every other node. In our simple 1-D case, the shape functions vary linearly from one to zero over the adjacent elements of each node. With these properties, the sum over all all shape functions is always one, and the coefficients of those shape functions happen to be the unknown variable values at the nodes, the nodal temperatures in our case.

Let’s plug the approximated solution into the strong form of our equation:

(92)\[R = \frac{\partial}{\partial x}k\frac{\partial N_j T_j }{\partial x}\neq 0.\]

Note how we use the Einstein convention and sum over repeated indices. We can now write the MWR form of the equation:

(93)\[\int_XW_i\frac{\partial}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx=0\ \ \ \ \ \ \ i=1,2,...,n\]

Galerkin method

In finite elements, the most frequently used weighting method is the Galerkin method that we already know from the previous session. We therefore use the interpolation functions also as weighting functions:

(94)\[\int_XN_i\frac{\partial}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx=0\ \ \ \ \ \ \ i=1,2,...,n\]

The next step is to reduce the order of differentiation by using partial integration. Remember, integration by parts looks like this:

(95)\[\int_{\Omega} fg' d\Omega = -\int_{\Omega} f'g d\Omega + \oint_{\Gamma} fg d\Gamma\]

Weak form

Applying this to equation (95) results in the weak form of our governing equation:

(96)\[\begin{split}\begin{align} \begin{split} \int_XN_i\frac{\partial}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx=0\ \ \ \ \ \ \ i=1,2,...,n\\ \Rightarrow \\ -\int_X \frac{\partial N_i}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx + \oint_{\Gamma} N_ik\frac{\partial N_j T_j }{\partial x}d\Gamma=0\ \ \ \ \ \ \ i=1,2,...,n\\ \Rightarrow \\ \int_X \frac{\partial N_i}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx + \oint_{\Gamma}\vec{q}\vec{n}d\Gamma=0\ \ \ \ \ \ \ i=1,2,...,n\\ \end{split} \end{align}\end{split}\]

Close inspection of the boundary integral reveals that it is the heat flow through the boundaries of the modeling domain. In 1-D these are simply numbers of heat flow at the sides and in 2-D/3-D heat fluxes through the sides of the modeling domain. For the moment we will neglect it and only continue with the solution inside our domain:

(97)\[\int_X \frac{\partial N_i}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx =0\ \ \ \ \ \ \ i=1,2,...,n\]

How can we solve this using the FEM? The basic idea is to split the integrals into subdomains and the solution will be the sum of the subdomains. In FEM, those subdomains are the finite elements:

(98)\[\begin{split}\begin{align} \begin{split} \int_X \frac{\partial N_i}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx =0\ \ \ \ \ \ \ i=1,2,...,n\\ \Rightarrow \\ \int_X \frac{\partial N_i}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx = \sum_{Elements} \int_{X_e} \frac{\partial N_i}{\partial x}k\frac{\partial N_j T_j }{\partial x}dx = 0\ \ \ \ \ \ \ i=1,2,...,n\\ \end{split} \end{align}\end{split}\]

While mathematically we are always allowed to split an integral into a sum of integrals, in FEM this is particularly useful. The catch is that the shape function associated with a node is only non-zero in the elements connected to that node. This in turn means that in equation (98), we can do the integration of each element completely independent from all other elements!

In the following sessions, we will go through all the steps involved in solving equation (98). We will first do this in a simple 1-D case before moving on to the general 2-D.

FEM implementation

Element stiffness matrix \(A_{el}\)

Let’s look at a single 1-D element. If we use linear shape functions, one 1-D element has two nodes – and only the shape functions of those two nodes will be non-zero. The so-called connectivity, which connects elements and nodes (which nodes belong to which element) is easy in 1-D: element number k will have the nodes i=k and i=k+1 (Element 1 has nodes 1 and 2). This means that for each element we will get two equations:

(99)\[\begin{split}\begin{bmatrix} \int_{X_e} \frac{\partial N_1}{\partial x}k\frac{\partial N_1}{\partial x}dx & \int_{X_e} \frac{\partial N_1}{\partial x}k\frac{\partial N_2}{\partial x}dx \\ \int_{X_e} \frac{\partial N_2}{\partial x}k\frac{\partial N_1}{\partial x}dx & \int_{X_e} \frac{\partial N_2}{\partial x}k\frac{\partial N_2}{\partial x}dx \end{bmatrix} \begin{bmatrix} T_1 \\ T_2 \end{bmatrix}=0\end{split}\]

We get a \(2x2\) so-called stiffness matrix for every element. In the end we will have to sum the contributions from every element into a global stiffness matrix, which will again be \([n x n]\). Note: The node numbering in equation (99) is local! I.e. node number 1 is the first node of an element k and has the global node number k, while node number 2 is the second node of the element and has the global node number k+1!

We have two linear shape functions for the points i and i+1 of each element (see Fig. 4):

(100)\[\begin{split}\begin{align} \begin{split} N_i(x) = 1 - \frac{x-x_i}{x_{i+1} - x_i} &= 1 - \frac{x-x_i}{\Delta x}\\ N_{i+1}(x) = \frac{x-x_i}{x_{i+1} - x_i} &= \frac{x-x_i}{\Delta x}\\ \end{split} \end{align}\end{split}\]

Matrix assembly

The above exercise results in the integrated element stiffness matrix of steady-state diffusion:

(101)\[\begin{split}Ael = \begin{bmatrix} \frac{k}{\Delta x} & -\frac{k}{\Delta x} \\ -\frac{k}{\Delta x} & \frac{k}{\Delta x} \end{bmatrix}\end{split}\]

The global stiffness matrix \(A\) is the sum of all element stiffness matrices; we simply have to add them all together using the connectivity information, i.e. which nodes belong to which elements (something that we call EL2NOD in the codes).

If we assume three linear elements, this looks like this

(102)\[\begin{split}A = \begin{bmatrix} \frac{k}{\Delta x} & -\frac{k}{\Delta x} & 0 & 0\\ -\frac{k}{\Delta x} & \frac{k}{\Delta x} & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & \frac{k}{\Delta x} & -\frac{k}{\Delta x} & 0 \\ 0 & -\frac{k}{\Delta x} & \frac{k}{\Delta x} & 0\\ 0 & 0 & 0 & 0\\ \end{bmatrix}+ \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & \frac{k}{\Delta x} & -\frac{k}{\Delta x} \\ 0 & 0 & -\frac{k}{\Delta x} & \frac{k}{\Delta x} \\ \end{bmatrix}\end{split}\]

Implementation

Let’s put everything together and write our first finite element code! We will do this by programming a FEM solution to the steady-state heat diffusion equation but this time with a source term.

(103)\[\frac{\partial}{\partial x}k\frac{\partial T_{ex}}{\partial x}+Q=0.\]

The weak form will then look like this (signs get flipped during partial integration):

(104)\[\begin{split}\begin{align} \begin{split} \int_X \frac{\partial N_i}{\partial x}k\frac{\partial N_j }{\partial x}T_j - N_iQdx =0\ \ \ \ \ \ \ i=1,2,...,n\\ \Rightarrow \\ \int_X \frac{\partial N_i}{\partial x}k\frac{\partial N_j }{\partial x} T_j - N_i Qd x = \sum_{Elements} \int_{X_e} \frac{\partial N_i}{\partial x}k\frac{\partial N_j }{\partial x} T_j - N_i Q dx = 0\ \ \ \ \ \ \ i=1,2,...,n\\ \end{split} \end{align}\end{split}\]

After another symbolic integration (try it) of the source term,we get the following element stiffness matrix and element system of equations:

(105)\[\begin{split}Ael = \begin{bmatrix} \frac{k}{\Delta x} & -\frac{k}{\Delta x} \\ -\frac{k}{\Delta x} & \frac{k}{\Delta x} \end{bmatrix} \begin{bmatrix} T_i \\ T_{i+1} \end{bmatrix} = \begin{bmatrix} \frac{\Delta x}{2}Q_{el} \\ \frac{\Delta x}{2}Q_{el} \end{bmatrix}\end{split}\]

Where the first matrix is the element stiffness matrix, the second vector are the unknown temperatures of the element, and the right-hand side has the integration source term \(Q_{el}\) which here is constant value (per element).

Let’s put this into python!