Basic FEM

Table of contents

The finite element method for linear elasticity

The partial differential equation for the static linear elasticity problem

$$\def\arraystretch{1.5} \begin{array}{lc} -\nabla \cdot \bm \sigma =f & \text{in }\Omega \\ \bm \sigma = 2\mu \bm \varepsilon + \lambda \bm I \nabla \cdot \bm u & \text{in }\Omega \\ \bm u = \bm g & \text{on }\partial \Omega_D \\ \bm \sigma \cdot \bm n = \bm t & \text{on }\partial \Omega_N \end{array}$$

By multiplying (1) with a test function $\bm v$ and integrating by parts, see here, we end up with the weak form of the linear elasticity equation, which states: Find the displacement field $\bm u$, with $\bm u = \bm g$ on the boundary $\partial \Omega_D$, such that

$$\boxed{ \int_{\Omega } \bm \sigma (\bm u) : \bm \varepsilon \left( \bm v\right) d\Omega =\int_{\Omega } \bm f\cdot \bm vd\Omega +\int_{\partial \Omega } \bm{t} \cdot \bm v ds \quad \forall \bm v \text{ that are } 0 \text{ on } \partial \Omega_D}$$

The weak form of the linear elasticity problem is posed on a domain $\Omega$ which needs to be discretized (meshed). This means dividing the domain into a finite amount of elements such that the problem can be tackled numerically, see Figure 1. In this section we will discuss this discretization and the choice of the weight function $v$ using the Galerkin method [2] for the 2D linear elasticity problem, but the reader is also referred to the Galerkin Method section for more information. For completeness we shall give a brief introduction to the method here.

The weighted residual method and the Galerkin method

The Galerkin method belongs to a family of weighted residual methods. Consider a differential equation of the form

$$- \dfrac{d^2}{dx^2} k u(x) = f(x) \quad \text{ for } x \in [a,b]$$

where enough boundary conditions are given such that a unique solution can exist, e.g., $u(a)=0$ and $du/dx(b)=0$. Now, an approximation method aims at solving (3). This is done by multiplying (3) by an arbitrary function $v(x)$ (the weight function) to obtain

$$v (f + k \dfrac{d^2 u}{dx^2} )$$

This can be integrated over the domain of interest

$$\boxed{- \int_a^b \dfrac{d^2 u}{dx^2} v dx = \int_a^b f v dx}$$

See here for an in-depth explanation of the reasons for integrating the equation. Integrating the term with the higher order derivative by parts leads to the weak form:

$$\int_a^b k \dfrac{d u}{d x} \dfrac{d v}{d x} dx = \int_a^b f v dx$$

In what follows we will not (yet) use the weak form, but just the integral form of the differential equation multiplied by the weight function to keep things general.

The solution to the differential equation can be approximated by a linear combination of some known functions $\varphi(x)$ and unknown parameters $a_i$:

$$u(x) \approx U(x) = \sum_{i=1}^n a_i \varphi_i(x)$$

where $a_1, ..., a_n$ are unknown parameters and $\varphi_1(x), ..., \varphi_n(x)$ are pre-specified functions, called trial or basis functions. The problem becomes to find these parameters $a_1, ..., a_n$ which is done by substituting $u(x)$ by $U(x)$ in (5), which leads to the weighted residual method,

$$\int_a^b \left( \dfrac{d^2}{dx^2} U(x) + f \right) v dx = 0$$

with $U(x)$ not satisfying the equation exactly in general, we get a residual

$$r(x) = \dfrac{d^2}{dx^2}U(x) + f \neq 0$$

which can be written as

$$\int_a^b r(x) v dx = 0$$

meaning that for a certain choice of $v$ determines the unknowns $a_1, ..., a_n$ and thus the approximating solution $U(x)$ by minimizing the signed integral. Let $v$ be a general function and constructed as a series of $n$ linear combinations of known functions $V_j(x)$ and arbitrary parameters $c_j$ such that

$$v(x) = V_1 c_1 + V_2 c_2 + ... + V_n c_n = \sum_j V_j c_j$$

Since the weight functions are arbitrary and $V_j(x)$ are known, it must be concluded that $c_j$ are arbitrary and moreover, since $c_j$ are not depended on $x$ we can pull them out of the integral and get

$$\sum_j c_j \int_a^b r V_j dx = 0$$

since this should hold for arbitrary $c_j$ we conclude that the following must hold

$$\int_a^b r V_j dx = 0$$

which in fact is a series of $n$ equations of the form

$$\begin{align*} \int_{a}^{b} r V_{1}dx & =0\\ \int_{a}^{b} r V_{2}dx & =0\\ \cdots\\ \int_{a}^{b} r V_{n}dx & =0 \end{align*}$$

substituting the approximation into the integral we get

$$\int_{a}^{b}\left(\sum_{i=1}^{n}\dfrac{d^{2}\varphi_{i}}{dx^{2}}a_{i}+f\right)V_{j}dx=0$$

where $a_i$ is not dependent on $x$ and thus can be pulled out of the integral, leading to the $n$-by-$n$ system

$$- \left( \int_{a}^{b}{\bf V}^{T}\dfrac{d^{2}\bm{\varphi}}{dx^{2}}dx \right) {\bf a } = \int_a^b {\bf V}^T f dx$$

where

$${\bf V}^T=\begin{bmatrix}V_{1}\\ V_{2}\\ \vdots\\ V_{n} \end{bmatrix}_{n \times 1},\quad\bm{\varphi}=[\varphi_{1},\varphi_{2},...,\varphi_{n}]_{1 \times n},\quad{\bf a}=\begin{bmatrix}a_{1}\\ a_{2}\\ \vdots\\ a_{n} \end{bmatrix}_{n \times 1}$$

and

$${\bf S} = - \int_{a}^{b}{\bf V}^{T}\dfrac{d^{2}\bm{\varphi}}{dx^{2}}dx,\quad{\bf f}=\int_{a}^{b}{\bf V}^{T}fdx$$

thus

$${\bf Sa}={\bf f}$$

The procedure leading from the differential equation multiplied by the weight function to the system of equations is known as the weighted residual method and specific methods are obtained by making a choice of the weight function $v$, e.g., the least squares method where $v = \frac{\partial r}{\partial a_i}$, but the most popular one for the FEM is the Galerkin method. See Chapter 8 in [1] for more information.

The Galerkin method

Galerkin suggested that given that the trial function or approximation has the form $U(x)=\sum_{i=1}^n a_i \varphi_i(x)$, the weight function $v$ should have the same form as the trial function, i.e., $v(x) = \sum^n_{i=1} c_i \varphi_i(x)$, meaning we should simply choose the basis functions as the weight functions, thus the Galerkin method takes the form:

$$\boxed{\int_{a}^{b}r\varphi_{j}dx=0}$$

and loosely speaking, setting $V_j = \varphi_j$, and in the case of the weak form the derivative of the basis weight function become $\frac{d V_j}{d x} = \frac{d \varphi_j}{d x}$ such that the weak form in (5) becomes:

$$\int_{a}^{b}k\dfrac{dU}{dx}\dfrac{d\varphi_{j}}{dx}dx=\int_{a}^{b}f\varphi_{j}dx$$

and substituting the approximation gives

$$\int_{a}^{b}k\left(\sum_{i=1}^{n}\dfrac{d\varphi_{i}}{dx}a_{i}\right)\dfrac{d\varphi_{j}}{dx}dx=\int_{a}^{b}f\varphi_{j}dx$$

which after pulling out $a_i$ of the integral, we get $j=1,...,n$ equations:

$$\sum_{i=1}^{n}\left(\int_{a}^{b}k\dfrac{d\varphi_{i}}{dx}\dfrac{d\varphi_{j}}{dx}dx\right)a_{i}=\int_{a}^{b}f\varphi_{j}dx$$

or using tensor notation

$$\left(\int_{a}^{b}k\dfrac{d\bm{\varphi}^{T}}{dx}\dfrac{d\bm{\varphi}}{dx}dx\right){\bf a}=\int_{a}^{b}\bm{\varphi}^{T}fdx$$

with

$$\boxed{ {\bf S}=\int_{a}^{b}k\dfrac{d\bm{\varphi}^{T}}{dx}\dfrac{d\bm{\varphi}}{dx}dx,\quad{\bf f}=\int_{a}^{b}\bm{\varphi}^{T}fdx }$$

The Galerkin method for the elasticity problem in 2D

In the previous section, Galerkin's method has only been discussed for approximations that are defined on the whole interval and also continuous on the whole interval. Since the first use was at a time where computers were not digital yet, the problems were kept small and the approximations simple. However, with the advent of the first digital computers came the use of Galerkin's method with piecewise polynomials which constitutes the finite element method.

To define a Galerkin Finite Element method for the elasticity problem in (2), we must introduce a finite dimensional subspace $V_h$, i.e., a set of piece-wise integrable functions, typically polynomials of low order (which translate to linear or quadratic triangle elements in 2D)

$$u\left(x,y\right)\approx U\left(x,y\right)=\sum_i \varphi_i \left(x,y\right)u_i$$

so that the functions $\varphi_i$ constitute a base for $V_h$, i.e., any $v\in V_h$ can be written as a linear combination of the basis functions $\varphi_i$ as discussed in the previous section. The interpolation points at $(x_i, y_i)$, $u_i \overset{2D}{=} {\left(\begin{array}{cc} u_x^i & u_y^i \end{array}\right)}^T$ are the nodal displacement vectors and make up the interpolant of $u$. The discrete displacement field for an $n$-noded mesh is represented by the displacement vector $\mathbf{u}={\left\lbrack \begin{array}{ccccccc} u_x^1 & u_y^1 & u_x^2 & u_y^2 & \cdots & u_x^n & u_y^n \end{array}\right\rbrack }^T$. This defines our discrete solution field on our mesh. We will in later sections take a look at the construction of elements and basis functions. Suffice to say for now we choose the basis functions, $\varphi_i$ to have a local support, meaning they have the value $\varphi_i =1$ in the node $i$ and $\varphi_i =0$ in all other nodes of the mesh. Thus $u_i$ has the value $u$ in node $i$. The value of $u$ between the nodes is interpolated using the basis function and we shall discuss the choice of polynomial order and its consequence on the approximation later.

As a first example: For a linear triangle, $K$, in 2D the displacement field takes the form ${\mathbf{u} }_K ={\left\lbrack \begin{array}{cccccc} u_x^1 & u_y^1 & u_x^2 & u_y^2 & u_x^3 & u_y^3 \end{array}\right\rbrack }^T$.

The displacement field on element $K$ is approximated using the interpolant, expressed on Voigt form using basis functions.

$$\def\arraystretch{1.5} u_K \left(\bm x\right)\approx U_K \left(\bm x \right)={\underset{\Phi_K \left(\bm x \right)}{\underbrace{\left\lbrack \begin{array}{cccccc} \varphi_1 \left(\bm x\right) & 0 & \varphi_2 \left(\bm x\right) & 0 & \varphi_3 \left(\bm x\right) & 0\\ 0 & \varphi_1 \left(\bm x\right) & 0 & \varphi_2 \left(\bm x\right) & 0 & \varphi_3 \left(\bm x\right) \end{array}\right\rbrack } } }_{2\times6 } {\left\lbrack \begin{array}{c} u_x^1 \\ u_y^1 \\ u_x^2 \\ u_y^2 \\ u_x^3 \\ u_y^3 \end{array}\right\rbrack }_{6 \times 1}$$ $$\def\arraystretch{1.5} =\left\lbrack \begin{array}{c} \varphi_1 u_x^1 +\varphi_2 u_x^2 +\varphi_3 u_x^3 \\ \varphi_1 u_y^1 +\varphi_2 u_y^2 +\varphi_3 u_y^3 \end{array}\right\rbrack =\left\lbrack \begin{array}{c} U_x \left(x,y\right)\\ U_y \left(x,y\right) \end{array}\right\rbrack \approx \left\lbrack \begin{array}{c} u_{K,x} \left(x,y\right)\\ u_{K,y} \left(x,y\right) \end{array}\right\rbrack$$

Looking at the virtual strain energy (or the double contraction between the stress and virtual strain)

$$\bm \sigma \left(\bm u\right): \bm \varepsilon \left( \bm v\right)=2\mu \bm \varepsilon(\bm u) : \bm \varepsilon(\bm v) +\lambda \mathrm{tr}\bm \varepsilon(\bm u) \bm I : \bm \varepsilon(\bm v)$$

We need to rewrite the weak form to facilitate the formulation of a FEM. Introducing the Voigt notation, in 2D we have

$$\def\arraystretch{1.5} \bm \varepsilon_V \left(\bm u\right)=\left\lbrack \begin{array}{c} \frac{\partial u_x }{\partial x}\\ \frac{\partial u_y }{\partial y}\\ \frac{\partial u_x }{\partial x}+\frac{\partial u_y }{\partial y} \end{array}\right\rbrack \text{ and } \mathbf D=\left\lbrack \begin{array}{ccc} 2\mu +\lambda & 0 & 0\\ 0 & 2\mu +\lambda & 0\\ 0 & 0 & \mu \end{array}\right\rbrack$$

Then we may write

$$\int_{\Omega } \bm \sigma \left(\bm u\right) : \bm\varepsilon \left(\bm v\right) d\Omega =\int_{\Omega } 2\mu \bm \varepsilon(\bm u) : \bm \varepsilon(\bm v) + \lambda \mathrm{tr} \bm \varepsilon(\bm u) \bm I : \bm \varepsilon(\bm v) d\Omega =\int_{\Omega } \bm \varepsilon_V {\left( \bm u \right)}^T \mathbf D \bm \varepsilon_V \left(\bm v\right) d\Omega$$

which is more suitable for numerical implementation.

$\bm \varepsilon_V {\left( \bm u \right)}$ can also be written

$$\bm \varepsilon_V {\left(\bm u \right)} = \left\lbrack \begin{array}{cc} \dfrac{\partial }{\partial x} & 0\\ 0 & \dfrac{\partial }{\partial y}\\ \dfrac{\partial }{\partial y} & \dfrac{\partial }{\partial x} \end{array}\right\rbrack \left\lbrack \begin{array}{c} u_x \\ u_y \end{array}\right\rbrack$$

which allows for breaking out the displacement field out of the integral. Introducing the discrete version of this on matrix form

$$\def\arraystretch{2.5} \mathbf B=\left\lbrack \begin{array}{cc} \dfrac{\partial }{\partial x} & 0\\ 0 & \dfrac{\partial }{\partial y}\\ \dfrac{\partial }{\partial y} & \dfrac{\partial }{\partial x} \end{array}\right\rbrack \left\lbrack \begin{array}{ccccc} \varphi_1 \left(\mathit{\mathbf{x} }\right) & 0 & \varphi_2 \left(\mathit{\mathbf{x} }\right) & 0 & \cdots \\ 0 & \varphi_1 \left(\mathit{\mathbf{x} }\right) & 0 & \varphi_2 \left(\mathit{\mathbf{x} }\right) & \cdots \end{array}\right\rbrack =\left\lbrack \begin{array}{ccccc} \dfrac{\partial \varphi_1 }{\partial x} & 0 & \dfrac{\partial \varphi_2 }{\partial x} & 0 & \cdots \\ 0 & \dfrac{\partial \varphi_1 }{\partial y} & 0 & \dfrac{\partial \varphi_2 }{\partial y} & \cdots \\ \dfrac{\partial \varphi_1 }{\partial y} & \dfrac{\partial \varphi_1 }{\partial x} & \dfrac{\partial \varphi_2 }{\partial y} & \dfrac{\partial \varphi_2 }{\partial x} & \cdots \end{array}\right\rbrack$$

The $\mathbf B$ matrix contains spatial derivatives of the basis function of our element. We shall see in the next section how the basis functions can be derived for any element and in the section after that how the derivatives are derived and computed.

With these definitions, the matrix formulation corresponding to the FEM formulation of (1) becomes

$$\boxed{ \left(\int_{\Omega } {\mathbf{B} }^T \mathit{\mathbf{D} }\;\mathbf{B}\;d\Omega \right)\mathbf{u}=\int_{\Omega } \Phi^T \mathit{\mathbf{f} }d\Omega +\int_{\partial \Omega } \Phi^T \mathit{\mathbf{t} }\;d\;s }$$

or

$$\boxed{ \mathbf{S}\;\mathbf{u}=\mathbf{f}+\mathbf{g} }$$

Mandel notation

An alternative to the Voigt notation is the Mandel notation

$$\bm \varepsilon_M {\left(\bm u \right)} =\left\lbrack \begin{array}{cc} \dfrac{\partial }{\partial x} & 0\\ 0 & \dfrac{\partial }{\partial y}\\ \dfrac{1}{\sqrt{2} }\dfrac{\partial }{\partial y} & \dfrac{1}{\sqrt{2} }\dfrac{\partial }{\partial x} \end{array}\right\rbrack \left\lbrack \begin{array}{c} u_x \\ u_y \end{array}\right\rbrack$$

and the corresponding discrete form

$$\def\arraystretch{2.5} \mathbf B_{\varepsilon } :=\left\lbrack \begin{array}{ccccc} \dfrac{\partial \varphi_1 }{\partial x} & 0 & \dfrac{\partial \varphi_2 }{\partial x} & 0 & \cdots \\ 0 & \dfrac{\partial \varphi_1 }{\partial y} & 0 & \dfrac{\partial \varphi_2 }{\partial y} & \cdots \\ \dfrac{1}{\sqrt{2} }\dfrac{\partial \varphi_1 }{\partial y} & \dfrac{1}{\sqrt{2} }\dfrac{\partial \varphi_1 }{\partial x} & \dfrac{1}{\sqrt{2} }\dfrac{\partial \varphi_2 }{\partial y} & \dfrac{1}{\sqrt{2} }\dfrac{\partial \varphi_2 }{\partial x} & \cdots \end{array}\right\rbrack$$

we have $\bm \varepsilon :\bm \varepsilon =\bm \varepsilon_M \cdot \bm \varepsilon_M$ and $\lambda \mathrm{tr} \bm \varepsilon \bm I : \bm \varepsilon$ can be written $\lambda \nabla \cdot \bm u \bm I : \bm \varepsilon$

$$\nabla \cdot \bm u=\left\lbrack \begin{array}{cc} \dfrac{\partial }{\partial x} & \dfrac{\partial }{\partial y} \end{array}\right\rbrack \left\lbrack \begin{array}{c} u_x \\ u_y \end{array}\right\rbrack$$

for which the corresponding discrete system is given by

$$\mathbf B_{\mathrm{div} } :=\left\lbrack \begin{array}{ccccc} \dfrac{\partial \varphi_{\;1} }{\partial x} & \dfrac{\partial \varphi_{\;1} }{\partial y} & \dfrac{\partial \varphi_{\;2} }{\partial x} & \dfrac{\partial \varphi_{\;2} }{\partial y} & \cdots \end{array}\right\rbrack$$

Using the Mandel notation, the discrete system becomes

$$\left(\int_{\Omega } 2\mu {\mathit{\mathbf{B} } }_{\varepsilon }^T {\mathit{\mathbf{B} } }_{\varepsilon } +\lambda {\mathit{\mathbf{B} } }_{\mathrm{div} }^T {\mathit{\mathbf{B} } }_{\mathrm{div} } d\Omega \right)\mathbf{u}=\int_{\Omega } \Phi^T \mathit{\mathbf{f} }d\Omega +\int_{\partial \Omega } \Phi^T \mathit{\mathbf{t} }\;d\;s$$

Note that the linear system is the same both in the Voigt and Mandel form, however there are some benefits with the Mandel form. The expression is kept in the general form of the Hooke's Law. The expression is separated into the volumetric and deviatoric terms, which facilitates special treatment of integration known as hourglass control.

Element formulations

For an element $K$, e.g., linear 2D element

We have (using the Mandel notation), the element stiffness matrix

$$\mathbf S_K =\int_K \left(2\mu {\mathit{\mathbf{B} } }_{\varepsilon }^T {\mathit{\mathbf{B} } }_{\varepsilon } +\lambda {\mathit{\mathbf{B} } }_D {\mathit{\mathbf{B} } }_D \right)d\;K$$

the element load vector

$$\mathbf f_K =\int_K \Phi^T \mathit{\mathbf{f} }\left(x\right)\;d\;K$$

the traction (external force) vector

$$\mathbf g_E =\int_E \Phi_E^T \mathit{\mathbf{t} }\left(x\right)\;d\;E$$

where $E$ denotes edge (in 2D, where as in 3D it would be a surface). The edge element is one spatial dimension lower, we are formulating the equations on an edge instead of a triangle. For the edge element we have

$$\bm \Phi_E \overset{2D\;\mathrm{lin} }{=} \left\lbrack \begin{array}{cccc} \varphi_1^E & 0 & \varphi_2^E & 0\\ 0 & \varphi_1^E & 0 & \varphi_2^E \end{array}\right\rbrack ,\; \bm \varphi^E \overset{2D\;\mathrm{lin} }{=} \left\lbrack \begin{array}{cc} 1+\xi & \xi \end{array}\right\rbrack$$

Assembling these element matrices into the global system leads to the system

$$\mathbf{S}\;\mathbf{u}=\mathbf{f}+\mathbf{g}$$

References