2-D FEM - elliptic problems
Let’s move on to 2D! We will again use our steady-state heat diffusion equation as an example and will follow basically the same steps as during the 1-D example.
Governing equation - strong form
The strong form looks like this, we have just added the second dimension.
2-D shape functions
We now use 2-D shape functions in our approximate solution:
where the shape functions \(N_j(x,y)\) are now function of both spatial dimensions. An example of a 2-D FEM mesh and associated shape functions are shown in Fig. 5. Note that the structured quad mesh and the bi-linear shape functions are just one possibility of using the FEM in 2-D. There are many other possibility, like unstructured triangle meshes and higher order shape functions.
Tip
Never forget that the shape functions \(N(x,y)\) are defined over the entire modeling domain in global coordinates. We will just evaluate them per element using local coordinates.
Note that the 2-D shape functions still sum to \(1\) everywhere:
Weak form
We proceed by substituting equation (116) into equation (115) and by using the Galerkin method we get:
We do the integration by parts and obtain:
Note how the signs flipped during integration by parts and by substituting \(\vec{q}=-k\nabla T\). The line integral has a clear physical meaning now, it’s the heat flow in and out of the modeling domain.
We go on and split the modeling domain into finite elements and the first integral turns into a sum over those elements. The line integral is again obmitted for the moment.
Element connectivity
The basic concept of the finite element method is to solve/assemble the system of equations, for e.g. heat diffusion, on each element and add all element contributions together to obtain the global matrix equation. Every element has an element number and a certain number of nodes. We will initially use quadratic elements with four nodes.
For a FEM mesh with bi-linear quad elements, we will need two matrices: a matrix we will call CGCOORD that contains the global coordinates of each grid point (it therefore has the size [nnod,2], where 2 is the number of dimensions (x,y) and nnod is the total number of nodes) and a matris so we call ELEM2NODE, which contains for each element the indices of the four nodes that belong to it (it therefore has the size [nel, nnodel], where nel is the total number of elements and nnodel is the number of nodes per element). See Fig. 10 for more details. We had already used these matrices in 1-D but now their meanings become clear.
We will use a connectivity as shown in Fig. 13. Element 0 has, for example, the global nodes 0, 1, 6, 5. Note that the local element numbering is always counterclockwise (in this class). It therefore contributes to the calculations of those four temperatures. The element stiffness matrix will now be [4,4]. The contribution of all elements is added/assembled into the glob
Numerical integration
The next step is to learn how we can integrate the element stiffness matrix. We have seen that the general finite element form, called weak form, of the heat diffusion equation can be written as (with boundary terms omitted):
The main problem with solving equation (121) is that the integral may be very difficult to evaluate. In finite elements we can use elements of very different geometries to mesh complex objects, which makes the 2-D integrals very hard to solve. We will therefore use numerical integration techniques like Gauss-Legendre quadrature. The problem is (or good thing as we will see later) that basically all quadrature rules (for quads) apply to rectangular domains that are 2x2 in size and have a local coordinate system \((−1 \leq (\xi,\eta) \leq 1)\). The integration of a function over a 2d domain can than be expressed as:
where M and N denote the number of Gauss points, \((\xi,\eta)\) denote the Gauss point coordinates, and \(W_I\) and \(W_J\) denote the corresponding Gauss weights. The selection of the number of Gauss points is based on the formula \(N=int[(p+1)/2]+1\), where p the polynomial degree to which the integrand is approximated. We will use linear shape functions so that we have 2 Gauss points (in every direction). In this case the local coordinates are \(\pm0.5773502692\) and the weights are all \(1\).
Gauss points
Order |
Integration point |
\(\xi\) |
\(\eta\) |
weight |
---|---|---|---|---|
\(1\) |
\(1\) |
\(0\) |
\(0\) |
\(2\) |
\(2\) |
\(1\) |
\(-0.577\) |
\(-0.577\) |
\(1\) |
\(2\) |
\(2\) |
\(+0.577\) |
\(-0.577\) |
\(1\) |
\(3\) |
\(1\) |
\(-0.775\) |
\(-0.775\) |
\(0.556\) |
\(3\) |
\(2\) |
\(0\) |
\(-0.775\) |
\(0.889\) |
\(3\) |
\(3\) |
\(+0.775\) |
\(-0.775\) |
\(0.556\) |
Isoparametric elements
So far so good, but how do we solve the integral in equation (121) using local coordinates? We will have to do a coordinate transformation for the purpose of numerical integration that relates the element \(\Omega_e\) to the master element area \(\hat{\Omega}\) - or, equivalently, \((x,y)\) to \((ξ,η)\). This can be done by using interpolation functions in terms of local coordinates. We will use so-called isoparametric elements in which we use the same number of interpolation functions for coordinate mapping and primary variable interpolation. This implies for our four node element that the local variables vary bi-linearly over the element.
where \(x_j\) are the nodal coordinates. Our primary variable (temperature) interpolation functions are still in global coordinates:
We can express the shape functions N in terms of local coordinates using equation (123) and the chain rule of differentiation:
Or in matrix form:
Jacobian matrix
equation (126) gives the relationship between the derivatives of \(N\) with respect to the global and local coordinates. The matrix in equation (126) is called the Jacobian matrix, J, of the coordinate transformation. Unfortunately this relationship is in the wrong direction. If we have a look at our starting equation equation (121), we realize that we need to relate the derivative with respect to global coordinates to the derivative with respect to local coordinates (which is the opposite direction). Fortunately this can be done quite easily:
Using equation (123) we can spell out the different components oft the Jacobian matrix:
Finally the element area dA=dxdy of the original element in global coordinates is transformed to:
We have now all the information we need to solve the integral in equation (121) numerically:
It should be noted that the transformation of a quadrilateral element of a mesh to a master element \(\hat{\Omega}\) is solely for the purpose of numerically evaluating the integrals in equation (121) . No transformation of the physical domain or elements is involved in the finite element analysis.
Excercises
FEM mesh
The first excercise is about building the FEM mesh and element connectivity matrix that contains the information on which nodes belong to which element. Have a look at Fig. 10 and follow the logic of that Figure.
Create the global coordinate vector GCOORD. Assume that the length of the box in x direction (lx) is 1 ([0,1]) and the length of box is y direction (ly) is also 1 ([0, 1]). The number of nodes in x is 5 (nx) and the number of nodes in y is 4 (ny).
Create the connectivity matrix ELEM2NODE. The element numbering should look like this:
Tip
Hint loop over all nodes and create their respective x and y coordinates. The node numbering should be like in Fig. 10 (see jupyter notebook for a discussion).
Functions like ceil or floor might help you with the connectivity matrix; it might also help to first compute the row we are in and then compute the index of the lower left node and move on from there.
Shape functions
The second excercise is about spelling out the shape functions and their derivatives. Remember that the shape functions are defined in local coordinates. For bi-linear quads, they have this functional form:
where \(T_{j=1..4}\) are the four nodal temperatures, \(\tilde{T}\) is the temperature at an (integration) point inside the element, \(N\) are the four shape functions, and \((\xi,\eta)\) are the two local coordinates between -1 and 1.
Implement the shape functions in a jupyter notebook, also compute the eight derivatives with local coordinates \((\xi,\eta)\), and try them out on a single element with coordinates -1,-1 to 1,1.
Implementation
Now we have assembled all the pieces that we need to make a finite element model and solve it. Listing 1 shows the complete python code for our first solver.
Python code
1import numpy as np
2from tabulate import tabulate
3from scipy.sparse.linalg import spsolve
4from scipy.sparse import csr_matrix
5import matplotlib as mpl
6import matplotlib.pyplot as plt
7from shapes import shapes
8
9#geometry
10lx = 1
11ly = 1
12nx = 51
13ny = 51
14nnodel = 4
15dx = lx/(nx-1)
16dy = ly/(ny-1)
17w_i = 0.2 # width inclusion
18h_i = 0.2 # heigths inclusion
19
20# model parameters
21k1 = 1
22k2 = 0.001
23Ttop = 0
24Tbot = 1
25
26nex = nx-1
27ney = ny-1
28nnod = nx*ny
29nel = nex*ney
30GCOORD = np.zeros((nnod,2))
31T = np.zeros(nnod) #initial T, not strictly needed
32
33# global coordinates
34id = 0
35for i in range(0,ny):
36 for j in range(0,nx):
37 GCOORD[id,0] = -lx/2 + j*dx
38 GCOORD[id,1] = -ly/2 + i*dy
39 id = id + 1
40
41# FEM connectivity
42EL2NOD = np.zeros((nel,nnodel), dtype=int)
43
44for iel in range(0,nel):
45 row = iel//nex
46 ind = iel + row
47 EL2NOD[iel,:] = [ind, ind+1, ind+nx+1, ind+nx]
48
49# Gauss integration points
50nip = 4
51gauss = np.array([[ -np.sqrt(1/3), np.sqrt(1/3), np.sqrt(1/3), -np.sqrt(1/3)], [-np.sqrt(1/3), -np.sqrt(1/3), np.sqrt(1/3), np.sqrt(1/3)]]).T.copy()
52
53# Storage
54Rhs_all = np.zeros(nnod)
55I = np.zeros((nel,nnodel*nnodel))
56J = np.zeros((nel,nnodel*nnodel))
57K = np.zeros((nel,nnodel*nnodel))
58
59for iel in range(0,nel):
60 ECOORD = np.take(GCOORD, EL2NOD[iel,:], axis=0 )
61 Ael = np.zeros((nnodel,nnodel))
62 Rhs_el = np.zeros(nnodel)
63
64 for ip in range(0,nip):
65 # 1. update shape functions
66 xi = gauss[ip,0]
67 eta = gauss[ip,1]
68 N, dNds = shapes(xi, eta)
69
70 # 2. set up Jacobian, inverse of Jacobian, and determinant
71 Jac = np.matmul(dNds,ECOORD) #[2,nnodel]*[nnodel,2]
72 invJ = np.linalg.inv(Jac)
73 detJ = np.linalg.det(Jac)
74
75 # 3. get global derivatives
76 dNdx = np.matmul(invJ, dNds) # [2,2]*[2,nnodel]
77
78 # 4. compute element stiffness matrix
79 kel = k1
80 if abs(np.mean(ECOORD[:,0]))<w_i and abs(np.mean(ECOORD[:,1]))<h_i:
81 kel=k2
82 Ael = Ael + np.matmul(dNdx.T, dNdx)*detJ*kel # [nnodel,1]*[1,nnodel] / weights are missing, they are 1
83
84 # 5. assemble right-hand side, no source terms, just here for completeness
85 Rhs_el = Rhs_el + np.zeros(nnodel)
86
87 # assemble coefficients
88 I[iel, :] = np.tile(EL2NOD[iel, :], (nnodel, 1)).T.ravel()
89 J[iel, :] = np.tile(EL2NOD[iel, :], (nnodel, 1)).ravel()
90 K[iel, :] = Ael.ravel()
91
92 Rhs_all[EL2NOD[iel,:]] += Rhs_el
93
94
95# Create the global stiffness matrix using a sparse representation
96A_all = csr_matrix((K.ravel(), (I.ravel(), J.ravel())), shape=(nnod, nnod))
97
98# indices and values at top and bottom
99i_bot = np.arange(0,nx, dtype=int)
100i_top = np.arange(nx*(ny-1),nx*ny, dtype=int)
101Ind_bc = np.concatenate((i_bot, i_top))
102Val_bc = np.concatenate((np.ones(i_bot.shape)*Tbot, np.ones(i_top.shape)*Ttop ))
103
104# smart way of boundary conditions that keeps matrix symmetry
105Free = np.arange(0,nnod)
106Free = np.delete(Free, Ind_bc)
107TMP = A_all[:,Ind_bc]
108Rhs_all = Rhs_all - TMP.dot(Val_bc)
109
110# solve reduced system
111T[Free] = spsolve(A_all[np.ix_(Free, Free)],Rhs_all[Free])
112T[Ind_bc] = Val_bc
113
114# postprocessing - heat flow
115Q_x = np.zeros(nel)
116Q_y = np.zeros(nel)
117Ec_x = np.zeros(nel)
118Ec_y = np.zeros(nel)
119
120for iel in range(0,nel):
121 # 0. get element coordinates
122 ECOORD = np.take(GCOORD, EL2NOD[iel,:], axis=0 )
123
124 # 1. update shape functions
125 xi = 0
126 eta = 0
127 N, dNds = shapes(xi, eta)
128
129 # 2. set up Jacobian, inverse of Jacobian, and determinant
130 Jac = np.matmul(dNds,ECOORD) #[2,nnodel]*[nnodel,2]
131 invJ = np.linalg.inv(Jac)
132 detJ = np.linalg.det(Jac)
133
134 # 3. get global derivatives
135 dNdx = np.matmul(invJ, dNds) # [2,2]*[2,nnodel]
136
137 # 4. heat flux per element
138 kel = k1
139 if abs(np.mean(ECOORD[:,0]))<0.25 and abs(np.mean(ECOORD[:,1]))<0.25:
140 kel=k2
141 Q_x[iel] = -kel*np.matmul(dNdx[0,:], np.take(T, EL2NOD[iel,:]))
142 Q_y[iel] = -kel*np.matmul(dNdx[1,:], np.take(T, EL2NOD[iel,:]))
143 Ec_x[iel] = np.mean(ECOORD[:,0])
144 Ec_y[iel] = np.mean(ECOORD[:,1])
145
146# plotting
147fig = plt.figure()
148left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
149ax = fig.add_axes([left, bottom, width, height])
150
151X = np.reshape(GCOORD[:,0], (ny,nx))
152Y = np.reshape(GCOORD[:,1], (ny,nx))
153T = T.reshape((ny,nx))
154
155cp = plt.contourf(X, Y, T, cmap='gist_yarg')
156plt.colorbar(cp)
157plt.quiver(Ec_x, Ec_y, Q_x, Q_y, np.sqrt(np.square(Q_x) + np.square(Q_y)), cmap='hot')
158
159ax.set_title('Temperature with heat flow vectors')
160ax.set_xlabel('x')
161ax.set_ylabel('y')
162plt.show()
The key components of the code are an element loop (which corresponds to the sum over the integrals in equation (120) and an integration loop (which corresponds to the loop over integration points in equation (120)). Most of the work is done inside the integration loop:
shape functions are calculated at the local coordinates of the integration point
The Jacobian matrix, its determinant, and the inverse Jacobian are is calculated according to equation (127)
Global derivatives of shape functions at local coordinates are calculated according to equation (128)
The element stiffness matrix and the element force vector are integrated.
After the integration
the local stiffness matrices and their global coefficients are stored
and the element force vector is added to the global force vector
Let’s have a closer look at some of the different steps.
Jacobian:
The calculation of the Jacobian matrix is implemented via equation (127) as a matrix-vector multiplication:
which is implemented as:
ECOORD = np.take(GCOORD, EL2NOD[iel,:], axis=0 )
Jac = np.matmul(dNds,ECOORD) #[2,nnodel]*[nnodel,2]
Global derivatives of shape functions
The global derivatives are calculated according to equation (128), which in matrix form is:
which is implemented as:
dNdx = np.matmul(invJ, dNds) # [2,2]*[2,nnodel]
Element stiffness matrix
The final step is to put the equation (130) into finite element form. Also this now quite easy:
where the first matrix is the element stiffness matrix, which is implemented as:
Ael = Ael + np.matmul(dNdx.T, dNdx)*detJ*kel # [nnodel,1]*[1,nnodel] / weights are missing, they are 1
Boundary conditions
Before we can solve anything, we need to set boundary conditions. In the 1-D case, we had set the boundary conditions in the usual way by setting the entire row to zero, putting a 1 on the main diagonal, and setting the Rhs to the specified temperature. We had noted before that this procedure breaks matrix symmetry. We will therefore use a more elaborate method here, which was described in [Dabrowski et al., 2008] . We add the columns that refer to nodes with boundary conditions to the righ-hand side and afterwards solve the smaller system of equations! This is implemented as
# indices and values at top and bottom
i_bot = np.arange(0,nx, dtype=int)
i_top = np.arange(nx*(ny-1),nx*ny, dtype=int)
Ind_bc = np.concatenate((i_bot, i_top))
Val_bc = np.concatenate((np.ones(i_bot.shape)*Tbot, np.ones(i_top.shape)*Ttop ))
# smart way of boundary conditions that keeps matrix symmetry
Free = np.arange(0,nnod)
Free = np.delete(Free, Ind_bc)
TMP = A_all[:,Ind_bc]
Rhs_all = Rhs_all - TMP.dot(Val_bc)
# solve reduced system
T[Free] = spsolve(A_all[np.ix_(Free, Free)],Rhs_all[Free])
T[Ind_bc] = Val_bc
First we find all the global indices of the nodes with boundary conditions and build a vector of equal length with all the fixed values. Afterwards we move the columns to the Rhs and solve the smaller system of equations. Note how the Rhs is a vector full of zeros at the moment. This will change when we look at the transient problem!
Excercise
You can download the python code here (2d_fem_steady_state.py
and shapes.py
). Try to get it to work and make sure that you understand the key parts. I good way is to use the Visual Studio Code debugger and to step through the code. Make sure that you understand the various matrix shapes and the connection between the equations presented above and the matrix multiplications in the code.