Basic FEA Theory – Part 5 – Solve For Displacements

Introduction

This article is the last article of my series on basic FEA theory. The series started with two Python code examples (Full Finite Element Solver in 100 Lines of Python and Full Finite Element Solver in 200 Lines of Python), and articles on (1) virtual work, (2) isoparametric shape functions, (3) strain calculations, and (4) element stiffness matrix.

The focus of this article is on how to put everything together and solve for the displacements. I will also go through the Python code for the whole FE solver.

Final Formulation of the Principle of Virtual Work (for One Element)

From my previous article, we know that the principal of virtual work can be written: \[\displaystyle \int_{V^{(m)}} \delta\boldsymbol{\varepsilon}^\top \boldsymbol{\sigma} dV =\int_{V^{(m)}} \delta\mathbf{u}^\top \mathbf{b} dV + \int_{A^{(m)}} \delta \mathbf{u}^\top \mathbf{t} dA + \delta \mathbf{u}^\top \mathbf{f}.\] Earlier I also showed that \(\mathbf{u}=\mathbf{H}\mathbf{U}\), from which we get \(\delta\mathbf{u}^\top = \delta\mathbf{U}^\top \mathbf{H}^\top\).

From this we can write the principal of virtual work in matrix form: \[\delta\mathbf{U}^\top \left[ \displaystyle\sum_{p}^{np} \mathbf{B}^\top \mathbf{C} \mathbf{B} J w_p \right]\mathbf{U} = \delta\mathbf{U}^\top \left[ \sum_p^{np} \mathbf{H}^\top \mathbf{b} J w_p \right] + \delta\mathbf{U}^\top \left[ \sum_p^{np} \mathbf{H}^\top \mathbf{t} J w_p \right] + \delta\mathbf{U}^\top \mathbf{f}.\]

Note that \(\delta \mathbf{U}^\top\) can be cancelled out from this equation! Super cool!

The final form therefore becomes: \[\left[ \displaystyle\sum_{p}^{np} \mathbf{B}^\top \mathbf{C} \mathbf{B} J w_p \right]\mathbf{U} = \left[ \sum_p^{np} \mathbf{H}^\top \mathbf{b} J w_p \right] + \left[ \sum_p^{np} \mathbf{H}^\top \mathbf{t} J w_p \right] + \mathbf{f}.\] This equation can be written in short form as: \( \mathbf{K} \mathbf{U} = \mathbf{R}\). The varibles \(\mathbf{K}\) and \(\mathbf{R}\) are known, and we can solve for the displacements \(\mathbf{U}\). That is the essence of how the linear FE method works!

Python Code Review

In this section I will review the original Python code that I started this series with.

				
					def shape(xi):
	x,y = tuple(xi)
	N = [(1.0-x)*(1.0-y), 
	(1.0+x)*(1.0-y), (1.0+x)*(1.0+y), 
	(1.0-x)*(1.0+y)]
	return 0.25*np.array(N)
def gradshape(xi):
	x,y = tuple(xi)
	dN = [[-(1.0-y),  (1.0-y), 
	(1.0+y), -(1.0+y)],
        [-(1.0-x), -(1.0+x), 
        (1.0+x),  (1.0-x)]]
	return 0.25*np.array(dN)

The first part of the code is the definition of the shape functions, and the gradient of the shape functions. \[N_1=0.25(1-\xi)(1-\eta)\] \[N_2=0.25(1+\xi)(1-\eta)\] \[N_3=0.25(1+\xi)(1+\eta)\] \[N_4=0.25(1-\xi)(1+\eta)\]

The Python code then defines the FE mesh (nodes, and elements). The next step is to calculate the \(\mathbf{C}\) matrix.

				
					E = 100.0
v = 0.48
C = E/(1.0+v)/(1.0-2.0*v) * np.array([
    [1.0-v,v,0.0], 
    [v,1.0-v,0.0],
    [0.0,0.0,0.5-v]])

\[\mathbf{C}^{(m)}_{ij} = \displaystyle\frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0\\ 0 & 0 & 0.5-\nu \end{bmatrix}\]

The main portion of the code is the formulation of the stiffness matrix. This is perhaps also the most challenging part to understand. Based on the review in the previous parts of this series we are now in a position to understand this in detail. The first step is formulate the \(\mathbf{B}\) matrix.

				
					dN = gradshape(q)
J  = np.dot(dN, xIe).T
dN = np.dot(np.linalg.inv(J), dN)
B[0,0::2] = dN[0,:]
B[1,1::2] = dN[1,:]
B[2,0::2] = dN[1,:]
B[2,1::2] = dN[0,:]

\[\mathbf{B}^{(m)}=\begin{bmatrix} M_{11} & 0 & M_{12} & 0 & M_{13} & 0 & M_{14} & 0 \\ 0 & M_{21} & 0 & M_{22} & 0 & M_{23} & 0 & M_{24} \\ M_{11} & M_{21} & M_{12} & M_{22} & M_{13} & M_{23} & M_{14} & M_{24} \end{bmatrix}\]

The element stiffness matrix can now be calculated in one line of Python code.

				
					Ke += np.dot(np.dot(B.T,C),B) * np.linalg.det(J)

\[\displaystyle\sum_{p}^{np} \mathbf{B}^\top \mathbf{C} \mathbf{B} J w_p\]

The global stiffness matrix is formulated to adding the individual element stiffness matrices.

				
					for i,I in enumerate(c):
	for j,J in enumerate(c):
		K[2*I,2*J]     += Ke[2*i,2*j]
		K[2*I+1,2*J]   += Ke[2*i+1,2*j]
		K[2*I+1,2*J+1] += Ke[2*i+1,2*j+1]
		K[2*I,2*J+1]   += Ke[2*i,2*j+1]
...

The final step is to simply solve the matrix equation. In Python this can be done on one line:

				
					u = np.linalg.solve(K, f)

Summary

In this series on the linear FE method I have shown how to write a complete linear FE solver in about 100 lines of Python code.
One reason for doing this is to illustrate the power of a high level computer language (like Python). I definitely recommend that you learn how to code in Python (or other similar languages, e.g. Julia).
Another reason for this series was to demonstrate that it is possible to understand the linear FE method from first principles without too much effort.
It is easy to extend the code presented there to cover other elements types. I will leave that as an exercise.
I will probably not do this, but it is rather easy to extend the linear FE method to a non-linear FE solver.