Intended learning outcomes

• Derive the weak form of the equilibrium from the strong form of the equilibrium equations

• Derive the finite element form from the weak form

• Calculate strains, stresses and element nodal force contributions for simple elements.

The Finite Element Method (FEM)

We use the finite element method to solve partial differential equations describing a boundary value problem. Here, we discuss the solution of the mechanical boundary value problem. The goal is only to briefly introduce the method, and not to discuss the many interesting aspects related to its implementation and use in various settings. This is left for specialized courses.

Strong form

The strong form of our boundary value problem is described here, and summarized by the following figure and equations.

$$\begin{aligned} \boldsymbol{ \sigma}(\boldsymbol{ \epsilon}(\underline{\boldsymbol{ u}}(\underline{\boldsymbol{ x}})))\cdot\underline{\boldsymbol{ \nabla}} + \underline{\boldsymbol{ b}} &= \underline{\boldsymbol{ 0}},\phantom{(\underline{\boldsymbol{ x}})} \quad \underline{\boldsymbol{ x}} \in \Omega \\ \underline{\boldsymbol{ u}}(x) &= \underline{\boldsymbol{ u}}_{ \mathrm{ D} }(\underline{\boldsymbol{ x}}), \quad \underline{\boldsymbol{ x}} \in \Gamma_{ \mathrm{ D} } \\ \underline{\boldsymbol{ n}} \cdot \boldsymbol{ \sigma}(\boldsymbol{ \epsilon}(\underline{\boldsymbol{ u}}(\underline{\boldsymbol{ x}}))) &= \underline{\boldsymbol{ t}}_{ \mathrm{ N} }(\underline{\boldsymbol{ x}}), \quad \underline{\boldsymbol{ x}} \in \Gamma_{ \mathrm{ N} } \end{aligned}$$

Weak form

The weak form is an integrated form of the strong form, which makes the problem possible to approximate using the Finite Element Method. To obtain this form, we apply two "tricks":

1. Scalar multiply the strong form with an arbitrary test displacement field $\underline{\boldsymbol{ \delta u}}(\underline{\boldsymbol{ x}})$

2. Integrate over the domain $\Omega$

Because the test displacement field $\underline{\boldsymbol{ \delta u}}(\underline{\boldsymbol{ x}})$ is arbitrary, the following weak form of the boundary value problem is in practice equivalent to Equation (1). (There are some restrictions wrt. to the smoothness of the functions $\underline{\boldsymbol{ u}}(\underline{\boldsymbol{ x}})$ and $\underline{\boldsymbol{ \delta u}}(\underline{\boldsymbol{ x}})$ from mathematics.) For brevity, from hereon we don't write out the dependence on $\underline{\boldsymbol{ x}}$ and $\boldsymbol{ \epsilon}$ for $\underline{\boldsymbol{ u}}$, $\underline{\boldsymbol{ \delta u}}$, $\boldsymbol{ \sigma}$, $\underline{\boldsymbol{ u}}_{ \mathrm{ D} }$, and $\underline{\boldsymbol{ t}}_{ \mathrm{ N} }$.

$$\begin{aligned} \int_\Omega \underline{\boldsymbol{ \delta u}} \cdot \left[ \boldsymbol{ \sigma}\cdot\underline{\boldsymbol{ \nabla}}\right]\, \mathrm{d} \Omega &= - \int_\Omega \underline{\boldsymbol{ \delta u}} \cdot \underline{\boldsymbol{ b}}\, \mathrm{d} \Omega \\ \underline{\boldsymbol{ u}} &= \underline{\boldsymbol{ u}}_{ \mathrm{ D} } \; \mathrm{on}\, \Gamma_{ \mathrm{ D} } \\ \underline{\boldsymbol{ n}} \cdot \boldsymbol{ \sigma}(\boldsymbol{ \epsilon}(\underline{\boldsymbol{ u}})) &= \underline{\boldsymbol{ t}}_{ \mathrm{ N} } \; \mathrm{on}\, \Gamma_{ \mathrm{ N} } \end{aligned}$$

If we then apply the Green-Gauss theorem to the first equation, we obtain

$$\begin{aligned} \int_\Gamma \underline{\boldsymbol{ \delta u}}\cdot \boldsymbol{ \sigma} \cdot \underline{\boldsymbol{ n}}\, \mathrm{d} \Gamma - \int_\Omega \left[ \underline{\boldsymbol{ \delta u}}\otimes\underline{\boldsymbol{ \nabla}}\right] : \boldsymbol{ \sigma}\, \mathrm{d} \Omega &= - \int_\Omega \underline{\boldsymbol{ \delta u}} \cdot \underline{\boldsymbol{ b}}\, \mathrm{d} \Omega \\ \end{aligned}$$

and we can insert the known traction on the Neumann boundary $\Gamma_{ \mathrm{ N} }$ as

$$\begin{aligned} \int_{\Gamma_{ \mathrm{ N} }} \underline{\boldsymbol{ \delta u}}\cdot \underline{\boldsymbol{ t}}_{ \mathrm{ N} }\, \mathrm{d} \Gamma_{ \mathrm{ N} } + \int_{\Gamma_{ \mathrm{ D} }} \underline{\boldsymbol{ \delta u}}\cdot \boldsymbol{ \sigma} \cdot \underline{\boldsymbol{ n}}\, \mathrm{d} \Gamma_{ \mathrm{ D} } - \int_\Omega \left[ \underline{\boldsymbol{ \delta u}}\otimes\underline{\boldsymbol{ \nabla}}\right] : \boldsymbol{ \sigma}\, \mathrm{d} \Omega &= - \int_\Omega \underline{\boldsymbol{ \delta u}} \cdot \underline{\boldsymbol{ b}}\, \mathrm{d} \Omega \\ \underline{\boldsymbol{ u}} &= \underline{\boldsymbol{ u}}_{ \mathrm{ D} } \; \mathrm{on}\, \Gamma_{ \mathrm{ D} } \end{aligned}$$

while keeping the known displacements on $\Gamma_{ \mathrm{ D} }$.

Finite element form

The weak form in Equation (4) is also very difficult to solve, perhaps arguably more difficult than the strong form in Equation (1) since we need to deal with the arbitrary function $\underline{\boldsymbol{ \delta u}}$. The main issue for both of these forms, is that the solution is a continuous function $\underline{\boldsymbol{ u}}$. It would be much easier to solve if we could somehow parameterize the function, e.g. if we could write that (in 1D) $u(x) = a_1 x + a_2 x^2 + a_3 x^3$ and then just determine the coefficients $a_i$. This is what we will do, namely that we will approximate both $\underline{\boldsymbol{ u}}$ and $\underline{\boldsymbol{ \delta u}}$ using so-called base functions:

$$\begin{aligned} \underline{\boldsymbol{ u}}(\underline{\boldsymbol{ x}}) &\approx \underline{\boldsymbol{ N}}_i(\underline{\boldsymbol{ x}}) a_i \\ \underline{\boldsymbol{ \delta u}}(\underline{\boldsymbol{ x}}) &\approx \underline{\boldsymbol{ N}}^\delta_i(\underline{\boldsymbol{ x}}) \delta a_i \end{aligned}$$

where $a_i$ and $\delta a_i$ are the coefficients that we need to determine and $\underline{\boldsymbol{ N}}_i$ and $\underline{\boldsymbol{ N}}^\delta_i$ are shape functions that we postulate can be used to describe the solution and test space. So we restrict our solution to solutions that can be described by this discretization. Here, we will use the standard Galerkin method of having $\underline{\boldsymbol{ N}}_i = \underline{\boldsymbol{ N}}^\delta_i$, which is described in detail in specific finite element courses. Inserting these approximations into our weak form, we obtained the discretized FE-problem

$$\begin{aligned} \delta a_i \left[\int_{\Gamma_{ \mathrm{ N} }} \underline{\boldsymbol{ N}}_i \cdot \underline{\boldsymbol{ t}}_{ \mathrm{ N} }\, \mathrm{d} \Gamma_{ \mathrm{ N} } + \int_{\Gamma_{ \mathrm{ D} }}\underline{\boldsymbol{ N}}_i \cdot \boldsymbol{ \sigma} \cdot \underline{\boldsymbol{ n}}\, \mathrm{d} \Gamma_{ \mathrm{ D} } - \int_\Omega \left[ \underline{\boldsymbol{ N}}_i\otimes\underline{\boldsymbol{ \nabla}}\right] : \boldsymbol{ \sigma}\, \mathrm{d} \Omega\right] &= - \delta a_i \int_\Omega \underline{\boldsymbol{ N}}_i \cdot \underline{\boldsymbol{ b}}\, \mathrm{d} \Omega \\ a_i \underline{\boldsymbol{ N}}_i = \underline{\boldsymbol{ u}}_{ \mathrm{ D} } \; \mathrm{on}\, \Gamma_{ \mathrm{ D} } \end{aligned}$$

Since $\underline{\boldsymbol{ \delta u}}$ should be arbitrary, so should $\delta a_i$. So by letting $\delta a_i=1$ and $\delta a_j=0$ for all $j\neq i$ for all $i$ (e.g. $\delta a_1=1$, then $\delta a_2=\delta a_3=\cdots=\delta a_n=0$ and $\delta a_2=1$, then $\delta a_1=\delta a_3=\cdots=\delta a_n=0$ and so on), we have $n$ equations ($i=1,2,\cdots,n$) that we can write as

$$\begin{aligned} f^{ \mathrm{ int} }_i &= f^{ \mathrm{ ext} }_i \\ f^{ \mathrm{ int} }_i &= \int_\Omega \left[ \underline{\boldsymbol{ N}}_i\otimes\underline{\boldsymbol{ \nabla}}\right] : \boldsymbol{ \sigma}\, \mathrm{d} \Omega \\ f^{ \mathrm{ ext} }_i &= \int_{\Gamma_{ \mathrm{ N} }} \underline{\boldsymbol{ N}}_i \cdot \underline{\boldsymbol{ t}}_{ \mathrm{ N} }\, \mathrm{d} \Gamma_{ \mathrm{ N} } + \int_{\Gamma_{ \mathrm{ D} }}\underline{\boldsymbol{ N}}_i \cdot \boldsymbol{ \sigma} \cdot \underline{\boldsymbol{ n}}\, \mathrm{d} \Gamma_{ \mathrm{ D} } + \int_\Omega \underline{\boldsymbol{ N}}_i \cdot \underline{\boldsymbol{ b}}\, \mathrm{d} \Omega \end{aligned}$$

where $f^{ \mathrm{ int} }$ and $f^{ \mathrm{ ext} }$ are the so-called internal (due to stresses) and external (due to applied loads) forces, respectively.

Base functions

The base functions $\underline{\boldsymbol{ N}}$ have many names, e.g. basis, shape, test, trial, etc. functions. These are generally used interchangably, potentially leading to some confusion. In general, test often denote the functions approximating the test function, $\underline{\boldsymbol{ \delta u}}$ (yes, the same name), i.e. $\underline{\boldsymbol{ N}}^\delta$. Trial functions often describe the functions approximating the solution $\underline{\boldsymbol{ u}}$, i.e. $\underline{\boldsymbol{ N}}$. Shape function is often used to describe the geometry. Finally, base or basis functions are often used to describe all of the above. For the $n$ equations $f^{ \mathrm{ int} }_i = f^{ \mathrm{ ext} }_i$ above to form a solvable equation system, we require (1) that the base functions are independent. That is, for each $j\in[1,n]$, there doesn't exist coefficients $a_i$ where $a_j=0$ such that $\underline{\boldsymbol{ N}}_j = a_i \underline{\boldsymbol{ N}}_i$. Secondly, we require that the solution is constrained by Dirichlet boundary conditions such that rigid body motions are prevented.

In practice, we often choose so-called nodal base functions. Consider the nodal positions $\underline{\boldsymbol{ x}}_i$. We then have a set of scalar base functions $N_i$, which are defined such that

$$\begin{aligned} N_i(\underline{\boldsymbol{ x}}_j) = \left\lbrace \begin{matrix} 1 & i=j \\ 0 & i\neq j \end{matrix}\right. \end{aligned}$$

In one dimension for linear base functions, this choice looks like

for 3 elements. If we consider the definition of the base function inside one element, between nodes $x_i$ and $x_{i+1}$, we have

$$\begin{aligned} N_i(x) = 1 - \frac{x-x_i}{x_{i+1}-x_{i}}, \quad N_{i+1}(x) = \frac{x-x_i}{x_{i+1}-x_{i}} \end{aligned}$$

As we increase to two dimensions, we have a more complicated situation. The left figure below shows only the triangular element, and then each of the base function in the other three figures.

The linear base functions inside an element with nodal coordinates $\underline{\boldsymbol{ x}}_1$, $\underline{\boldsymbol{ x}}_2$, and $\underline{\boldsymbol{ x}}_3$, ordered counter-clockwise around the element as in the figure above, are

$$\begin{aligned} N_1(x,y) &= \frac{x_2 y_3 - x_3 y_2 + [y_2-y_3]x + [x_3-x_2]y}{2A} \\ N_2(x,y) &= \frac{x_3 y_1 - x_1 y_3 + [y_3-y_1]x + [x_1-x_3]y}{2A} \\ N_3(x,y) &= \frac{x_1 y_2 - x_2 y_1 + [y_1-y_2]x + [x_2-x_1]y}{2A} \end{aligned}$$

where $A=\left[[x_2 y_3 - x_3 y_2] - [x_1 y_3 - x_3 y_1] + [x_1 y_2 - x_2 y_1]\right]/2$ is the element's area. It is left as an exercise for the reader to verify that $N_i(x_j,y_j)=\delta_{ij}$. For a detailed derivation, please see Niels Ottosen and Hans Petersson, 1992: "Introduction to the Finite Element Method", Section 7.3.1.

However, the base functions discussed so far are scalar base functions. But since we want to approximate vector functions, we need vector base functions. These are easily constructed for each node as

$$\begin{aligned} \underline{\boldsymbol{ N}}_{[3i-2]}(\underline{\boldsymbol{ x}}) &= [N_i(\underline{\boldsymbol{ x}}), 0, 0]\\ \underline{\boldsymbol{ N}}_{[3i-1]}(\underline{\boldsymbol{ x}}) &= [0, N_i(\underline{\boldsymbol{ x}}), 0] \\ \underline{\boldsymbol{ N}}_{[3i]}(\underline{\boldsymbol{ x}}) &= [0, 0, N_i(\underline{\boldsymbol{ x}})] \end{aligned}$$

for three dimensions. More formally and short, we can write $\underline{\boldsymbol{ N}}_i(\underline{\boldsymbol{ x}}) = N_j(\underline{\boldsymbol{ x}})\underline{\boldsymbol{ e}}_{ k}$ where $j$ is the node number, $m$ is the dimension, $1\leq k\leq m$, and $i=m(j-1)+k$. If we have $n_{ \mathrm{ nod} }$ nodes, we then have $n=m\,n_{ \mathrm{ nod} }$ base functions $\underline{\boldsymbol{ N}}_i$. And as mentioned above, there are equally many coefficients $a_i$ and equations $\underline{\boldsymbol{ f}}^{ \mathrm{ int} }_i=\underline{\boldsymbol{ f}}^{ \mathrm{ ext} }_i$ as base functions.

Evaluating the integrals

The force vectors $f^{ \mathrm{ int} }_i$ and $f^{ \mathrm{ ext} }_i$ in Equation (7) are calculated by solving the integral over the entire domain. But the construction of base functions above ensures that one base function is only non-zero in a few elements, as illustrated below.

This makes it possible to calculate contributions elementwise. If we evaluate $\underline{\boldsymbol{ N}}_i\otimes\underline{\boldsymbol{ \nabla}}$, for the base functions above, this becomes a constant tensor as $N_j(\underline{\boldsymbol{ x}})$ is linear inside the element. Consequently, the strain $\underline{\boldsymbol{ \epsilon}} = \left[ \underline{\boldsymbol{ u}}\otimes\underline{\boldsymbol{ \nabla}}\right]^{ \mathrm{sym} } = a_i \left[ \underline{\boldsymbol{ N}}_i\otimes\underline{\boldsymbol{ \nabla}}\right]^{ \mathrm{sym} }$ is also constant. And since the stress is a function of the strain, then $\boldsymbol{ \sigma}(\boldsymbol{ \epsilon})$ is constant inside the element, implying that the entire integrand $\left[ \underline{\boldsymbol{ N}}_i\otimes\underline{\boldsymbol{ \nabla}}\right] : \boldsymbol{ \sigma}$ is constant. So for triangular elements with linear base functions, the integrals are trivial to solve. For cases when the base functions are not linear, e.g. quadrilateral elements (4 nodes) or quadratic (6-noded) triangular elements, we need to use numerical integration techniques, see e.g. Niels Ottosen and Hans Petersson, 1992: "Introduction to the Finite Element Method", Chapter 20.

Solving a linear finite element problem

For the finite element problem to be linear, we require that the material behavior is linear. For this example, we will use linear elasticity, such that $\boldsymbol{ \sigma}=\textbf{\textsf{ E}}:\boldsymbol{ \epsilon}$. For this case, our $f^{ \mathrm{ int} }_i$ simplifies to

$$\begin{aligned} f^{ \mathrm{ int} }_i &= \int_\Omega \left[ \underline{\boldsymbol{ N}}_i\otimes\underline{\boldsymbol{ \nabla}}\right]^{ \mathrm{sym} } : \textbf{\textsf{ E}}: \left[ \underline{\boldsymbol{ N}}_j\otimes\underline{\boldsymbol{ \nabla}}\right]^{ \mathrm{sym} } \, \mathrm{d} \Omega a_j\\ \end{aligned}$$

as $\textbf{\textsf{ E}}$ is minor symmetric and $\boldsymbol{ \epsilon}=\left[ \underline{\boldsymbol{ N}}_j\otimes\underline{\boldsymbol{ \nabla}}\right]^{ \mathrm{sym} } a_j$. Note that everything inside the integral is now independent of the solution $a_j$, and we can write this as

$$\begin{aligned} K_{ij} &= \int_\Omega \left[ \underline{\boldsymbol{ N}}_i\otimes\underline{\boldsymbol{ \nabla}}\right]^{ \mathrm{sym} } : \textbf{\textsf{ E}}: \left[ \underline{\boldsymbol{ N}}_j\otimes\underline{\boldsymbol{ \nabla}}\right]^{ \mathrm{sym} } \, \mathrm{d} \Omega \\ f^{ \mathrm{ int} }_i &= K_{ij} a_j \end{aligned}$$

We can then remind ourselves that $f^{ \mathrm{ ext} }_i$ is given as

$$\begin{aligned} f^{ \mathrm{ ext} }_i &= \int_{\Gamma_{ \mathrm{ N} }} \underline{\boldsymbol{ N}}_i \cdot \underline{\boldsymbol{ t}}_{ \mathrm{ N} }\, \mathrm{d} \Gamma_{ \mathrm{ N} } + \int_{\Gamma_{ \mathrm{ D} }}\underline{\boldsymbol{ N}}_i \cdot \boldsymbol{ \sigma} \cdot \underline{\boldsymbol{ n}}\, \mathrm{d} \Gamma_{ \mathrm{ D} } + \int_\Omega \underline{\boldsymbol{ N}}_i \cdot \underline{\boldsymbol{ b}}\, \mathrm{d} \Omega \end{aligned}$$

which we will simply denote $f_i$ from now on. Our complete equation system is then

$$\begin{aligned} K_{ij} a_j &= f_i \\ \underline{\underline{ K}} \underline{ a} &= \underline{ f} \end{aligned}$$

where the latter equation is just in matrix-vector form.

However, for degrees of freedom $i$ on $\Gamma_{ \mathrm{ D} }$, $f_i$ is unknown. On the other hand, that implies that the corresponding $a_i$ is known. We can partition our equation system into "free" degrees of freedom, $\underline{a}_{ \mathrm{ F} }$, where the solution $a_i$ is unknown (but $f_i$ is known) and prescribed degrees of freedom, $\underline{a}_{ \mathrm{ P} }$, where the solution $a_i$ is known (but $f_i$ is unknown). This gives

$$\begin{aligned} \begin{bmatrix} \underline{\underline{ K}}_{ \mathrm{ FF} } & \underline{\underline{ K}}_{ \mathrm{ FP} } \\ \underline{\underline{ K}}_{ \mathrm{ PF} } & \underline{\underline{ K}}_{ \mathrm{ PP} } \end{bmatrix} \begin{bmatrix} \underline{ a}_{ \mathrm{ F} } \\ \underline{ a}_{ \mathrm{ P} } \end{bmatrix}= \begin{bmatrix} \underline{ f}_{ \mathrm{ F} } \\ \underline{ f}_{ \mathrm{ P} } \end{bmatrix} \end{aligned}$$

where $\underline{ a}_{ \mathrm{ P} }$ and $\underline{ f}_{ \mathrm{ F} }$ are known, while $\underline{ a}_{ \mathrm{ F} }$ and $\underline{ f}_{ \mathrm{ P} }$ are unknown. We can then solve this by first solving the linear equation system

$$\begin{aligned} \underline{\underline{ K}}_{ \mathrm{ FF} } \underline{ a}_{ \mathrm{ F} } = \underline{ f}_{ \mathrm{ F} } - \underline{\underline{ K}}_{ \mathrm{ FP} } \underline{ a}_{ \mathrm{ P} } \end{aligned}$$

to find $\underline{ a}_{ \mathrm{ F} }$. $\underline{ f}_{ \mathrm{ P} }$ is then given as

$$\begin{aligned} \underline{ f}_{ \mathrm{ P} } = \underline{\underline{ K}}_{ \mathrm{ PF} }\underline{ a}_{ \mathrm{ F} } + \underline{\underline{ K}}_{ \mathrm{ PP} }\underline{ a}_{ \mathrm{ P} } \end{aligned}$$