Primer on Finite Element Method

The Finite Element Method (FEM) is a numerical method for solving partial differential equations (PDEs) and integral equations. It is a powerful tool for solving problems in solid mechanics.

Here we present a short primer on FEM to bring the reader up to speed with the basics of FEM. The reader is assumed to have a basic understanding of calculus and linear algebra. For a detailed treatment of FEM, the reader is referred to the following books and courses:

Books:

Courses:

Finding Equilibrium

In solid mechanics, we are interested in finding the equilibrium configuration of a solid body. The equilibrium configuration of a body is described by the linear momentum balance equation:

\[ \nabla \cdot{} \boldsymbol{\sigma} + \boldsymbol{b} = \rho \ddot{\boldsymbol{u}} \]

where \(\boldsymbol{\sigma}\) is the stress tensor, \(\boldsymbol{b}\) is the body force per unit volume, \(\rho\) is the density, and \(\ddot{\boldsymbol{u}}\) is the acceleration.

The above equation is referred to as the strong form of the problem as it must be satisfied point-wise in the entire domain. For quasi-static problems, the acceleration term is zero and for simplicity we will assume that the body force is zero and the equation simplifies to:

\[ \nabla \cdot{} \boldsymbol{\sigma} = 0 \]

The stress tensor is related to the strain tensor by the constitutive law. For a linear elastic material, the stress tensor is given by:

\[ \boldsymbol{\sigma} = \lambda \text{tr}(\boldsymbol{\varepsilon})\mathbb{I} + 2 \mu \boldsymbol{\varepsilon}, \]

where \(\lambda\) and \(\mu\) are the first Lame’s constants and the shear modulus. \(\mathbb{I}\) is the identity tensor. \(\boldsymbol{\varepsilon}\) is the strain tensor and \(\text{tr}(\boldsymbol{\varepsilon})\) is the trace of the strain tensor. The strain tensor is given by:

\[ \boldsymbol{\varepsilon} = \dfrac{1}{2} \left( \nabla \boldsymbol{u} + \nabla \boldsymbol{u}^T \right) \]

where \(\nabla \boldsymbol{u}\) is the gradient of the displacement field.

The above strong form of the problem is often difficult to solve analytically and therefore, we rely on weak form or variational form of the problem.

The Principle of Minimum Energy

The equilibrium configuration of a body can also be described from a principle of minimum energy which states that the configuration (i.e. displacements \(\boldsymbol{u}\)) that makes the total potential energy functional stationary is the equilibrium configuration. The generic form of total potential energy functional of a body is given by:

\[ \Psi = \Psi_\text{int}(\boldsymbol{u}) + \Psi_\text{ext}(\boldsymbol{u}) \]

where \(\Psi_\text{int}(\boldsymbol{u})\) is the internal energy due to process inherent to the body and \(\Psi_\text{ext}(\boldsymbol{u})\) is the external energy due to the external forces. The internal energy of a body is often sum of energies due to different physical processes inherent to the body which is given by: \[ \Psi_\text{int}(\boldsymbol{u}) = \Psi_\text{e}(\boldsymbol{u}) + \ldots \]

where \(\Psi_\text{e}(\boldsymbol{u})\) is the strain energy due to the deformation of the material, \(\ldots\) represents other energy terms. During this course, we will come across different physical phenomena where this \(\ldots\) will be replaced by different energy terms. For example, in case of fracture problems, the \(\ldots\) will be replaced by the fracture energy which is dissipated during the fracture process. Similarly, in case of contact problems, the \(\ldots\) will be replaced by the contact energy.

Key Idea

Express the total potential energy functional \(\Psi\) as a sum of internal energy \(\Psi_\text{int}(\boldsymbol{u})\) and external energy \(\Psi_\text{ext}(\boldsymbol{u})\) as given below: \[ \Psi = \Psi_\text{int}(\boldsymbol{u}) ~ + ~ \Psi_\text{ext}(\boldsymbol{u}) \]

For now, we will assume that only the strain energy contributes to the internal energy. Therefore, the total potential energy functional is given by:

\[ \Psi = \int_{\Omega} \psi (\boldsymbol{u})~ \text{d}\Omega ~ \underbrace{- \int_{\Gamma} \boldsymbol{t} \cdot{} \boldsymbol{u} ~\text{d}\Gamma}_{\Psi_\text{ext}} \]

where \(\boldsymbol{t}\) is the traction vector on the boundary of the domain, \(\Gamma\) is the boundary of the domain, \(\Omega\) is the domain, \(\psi\) is the strain energy density, and \(\Psi(\boldsymbol{u})\) is the free energy or total potential energy functional. The strain energy density is given by:

\[ \psi (\boldsymbol{u}) = \dfrac{1}{2} \boldsymbol{\sigma}(\boldsymbol{\varepsilon}) : \boldsymbol{\varepsilon}(\boldsymbol{u}) \]

The stationary point for the energy functional is defined as the displacements \(\boldsymbol{u}\) that makes the first variation of the energy functional zero. Therefore, \(\boldsymbol{u}\) is an equilibrium configuration if and only if:

\[ \delta \Psi = \int_{\Omega} \delta \psi (\boldsymbol{u}) ~\text{d}\Omega - \int_{\Gamma} \delta (\boldsymbol{t} \cdot{} \boldsymbol{u}) ~\text{d}\Gamma = 0 \]

where \(\delta\) is the small variation in quantity caused by an arbitrary small change in the displacement field \(\delta\boldsymbol{u}\). Applying variational calculus, we can write the first variation of the total potential energy functional as:

\[ \delta \Psi = \int_{\Omega} \delta (\boldsymbol{\sigma}:\boldsymbol{\varepsilon}) \text{d}\Omega - \int_{\Gamma} \boldsymbol{t} \cdot{} \delta \boldsymbol{u} ~\text{d}\Gamma = 0 \]

\[ \delta \Psi = \int_{\Omega} \boldsymbol{\sigma}:\delta\boldsymbol{\varepsilon} \text{d}\Omega - \int_{\Gamma} \boldsymbol{t} \cdot{} \delta \boldsymbol{u} ~\text{d}\Gamma = 0 \]

Using the property \(\delta \boldsymbol{\varepsilon} = \delta \nabla \boldsymbol{u}=\nabla \delta \boldsymbol{u}\), we can express the above equation as:

\[ \delta \Psi = \int_{\Omega} \boldsymbol{\sigma} : \nabla \delta \boldsymbol{u} \text{d}\Omega - \int_{\Gamma} \boldsymbol{t} \cdot{} \delta \boldsymbol{u} ~\text{d}\Gamma = 0 \]

Note: The above equation is also known as the principle of virtual work.

Applying the identity \(\nabla (\boldsymbol{\sigma} :\delta \boldsymbol{u}) = \nabla \boldsymbol{\sigma} \cdot{} \delta \boldsymbol{u} + \boldsymbol{\sigma}: \nabla \delta \boldsymbol{u}\) and the divergence theorem, we can express the above equation as:

\[ \int_{\Omega} \nabla \cdot{} \boldsymbol{\sigma} \cdot{} \delta \boldsymbol{u} \text{d}\Omega - \int_{\Gamma} \boldsymbol{t} \cdot{} \delta \boldsymbol{u} ~\text{d}\Gamma = 0 \] This is the variational form or the weak form of the equilibrium equation. We can rewrite it as:

\[ \big(\int_{\Omega} \nabla \cdot{} \boldsymbol{\sigma}~\text{d}\Omega - \int_{\Gamma} \boldsymbol{t} ~\text{d}\Gamma \big) \cdot{} \delta \boldsymbol{u} = 0 \] which is the equilibrium equation of the body. And since \(\delta \boldsymbol{u}\) is arbitrary, we can conclude that:

\[ \underbrace{\int_{\Omega} \nabla \cdot{} \boldsymbol{\sigma}~\text{d}\Omega}_{\boldsymbol{f}_\text{int}} - \underbrace{\int_{\Gamma} \boldsymbol{t} ~\text{d}\Gamma}_{\boldsymbol{f}_\text{ext}} = \mathbf{0} \] where \(\boldsymbol{f}_\text{int}\) is the internal force and \(\boldsymbol{f}_\text{ext}\) is the external force. The above equation is the force balance equation between the internal forces and the external forces and is also the equilibrium equation of the body.

Key Idea

Therefore, we can conclude that the equilibrium configuration of a body is the configuration that makes the free energy stationary.

Based on calculus of variation, we can express the first variation of a functional in terms of its functional derivative (see Appendix C — The First Variation and the Functional Derivative for more details). The first variation of the total potential energy functional \(\delta \Psi = 0\) thus can mathematically be expressed as:

\[ \delta \Psi = \dfrac{\partial \Psi}{\partial \boldsymbol{u}} \delta \boldsymbol{u} = 0, \] Since the variation \(\delta \boldsymbol{u}\) is arbitrary, we can conclude that: \[ \forall \delta \boldsymbol{u} \implies \dfrac{\partial \Psi}{\partial \boldsymbol{u}} = \mathbf{0}. \] where the condition \(\frac{\partial \Psi}{\partial \boldsymbol{u}}\) represents the minima of the energy functional (given \(\frac{\partial^2 \Psi}{\partial \boldsymbol{u}^2} \geq 0\)). Thus, above mathematical identity proves that finding the stationary point for the energy functional is equivalent to minimizing the energy functional and finding the roots of the equation \(\frac{\partial \Psi}{\partial \boldsymbol{u}} = \mathbf{0}\). We can expand the equation by taking the derivatives of the strain energy and the external energy with respect to the displacements \(\boldsymbol{u}\) \[ \begin{align*} \frac{\partial \Psi}{\partial \boldsymbol{u}} &= \boldsymbol{0}\\ \frac{\partial \Psi_{e}}{\partial \boldsymbol{u}} + \frac{\partial \Psi_\text{ext}}{\partial \boldsymbol{u}} &= \boldsymbol{0}\\ \boldsymbol{f}_\text{int} - \boldsymbol{f}_\text{ext} &= \boldsymbol{0}\\ \end{align*} \]

Key Idea

The internal forces \(\boldsymbol{f}_\text{int}\) is the derivative of the strain energy with respect to the displacements \(\boldsymbol{u}\).

\[ \boldsymbol{f}_\text{int} = \dfrac{\partial \Psi_{e}}{\partial \boldsymbol{u}} \]

To solve the equation \(\boldsymbol{f}_\text{int} - \boldsymbol{f}_\text{ext} = \boldsymbol{0}\), one can use any derivative based method such as Newton’s method which states that,

\[ \Delta \boldsymbol{u} = \boldsymbol{u}_0 -(\dfrac{\partial f}{\partial \boldsymbol{u}})^{-1} f(\boldsymbol{u_0}) \] where \(\boldsymbol{u}_0\) is the initial guess for the displacements and \(\Delta \boldsymbol{u}\) is the increment in the displacements. Replacing \(f\) with \(\boldsymbol{f}_\text{int} - \boldsymbol{f}_\text{ext}\), and assuming \(\boldsymbol{u}_0\) as zeros, we can write the equation as:

\[ \begin{align*} \Delta \boldsymbol{u} &= -(\dfrac{\partial \boldsymbol{f}_\text{int}}{\partial \boldsymbol{u}})^{-1} (\boldsymbol{f}_\text{int} - \boldsymbol{f}_\text{ext})\\ &= \mathbf{K}^{-1} (\boldsymbol{f}_\text{ext} - \boldsymbol{f}_\text{int}) \end{align*} \]

where \(\frac{\partial \boldsymbol{f}_\text{int}}{\partial \boldsymbol{u}}\) is the stiffness matrix \(\mathbf{K}\). So, the stiffness matrix is the first derivative of the internal force vector with respect to the displacements.

Key Idea

The stiffness matrix is the first derivative of the internal force vector with respect to the displacements or the second derivative of the strain energy with respect to the displacements.

\[ \mathbf{K} = \dfrac{\partial \boldsymbol{f}_\text{int}}{\partial \boldsymbol{u}} = \dfrac{\partial^2 \Psi_{e}}{\partial \boldsymbol{u}^2} \]

So if we can compute the total strain energy \(\Psi_\text{e}(\boldsymbol{u})\), we can compute the internal forces \(\boldsymbol{f}_\text{int}\) and the stiffness matrix \(\mathbf{K}\) by differentiating the total strain energy with respect to the displacements \(\boldsymbol{u}\).

Note

In this course, to solve any problem, we will always follow these steps:

  • Compute the total internal energy \(\Psi_\text{int}(\boldsymbol{u})\)
  • Differentiate it with respect to the displacements \(\boldsymbol{u}\) to compute the internal forces \(\boldsymbol{f}_\text{int}\) and the stiffness matrix \(\mathbf{K}\).
  • Solve the equilibrium equation \(\boldsymbol{f}_\text{int} - \boldsymbol{f}_\text{ext} = \boldsymbol{0}\) to find the displacements \(\boldsymbol{u}\).

Using FEM to compute the total strain energy

We use Finite Element discretization to compute the total strain energy. In FEM, we discretize the domain into a finite number of elements. The collection of all elements and how each element is connected with its neighbors (referred as the connectivity) is called the mesh. Each element is a small subdomain of the domain. The domain is then approximated by a set of basis functions or shape functions. What discretization does is that instead of evaluating the displacement field at every point in the domain, we evaluate it at the nodes of the mesh. The values residing at the nodes of the mesh are called the nodal values. We will represent the nodal values as \(\hat{\boldsymbol{u}}\). The nodal values for an element \(e\) is denoted as \(\hat{\boldsymbol{u}}_e\).

Discretization of the domain into elements. Each finite element in physical space (shown in blue) is mapped to a reference element (shown in green) in the reference space.

Integrating the strain energy over the domain

The strain energy of a domain which is given as:

\[ \Psi_\text{e}(\boldsymbol{u}) = \int_{\Omega} \frac{1}{2} \boldsymbol{\sigma} : \boldsymbol{\varepsilon} \, \text{d}\Omega \]

where \(\boldsymbol{\sigma}\) is the stress tensor and \(\boldsymbol{\varepsilon}\) is the strain tensor evaluated at the quadrature points. The finite element method way of integrating this over a discretized domain is to approximate the integral as a sum of integrals over the elements in the domain.

\[ \Psi_\text{e}(\boldsymbol{u}) \approx \sum_{e \in \mathcal{E}} \sum_{\xi \in \mathcal{Q}} \frac{1}{2} \boldsymbol{\sigma}(\xi) : \boldsymbol{\varepsilon}(\xi) \, \mathrm{det}\boldsymbol{J}(\xi) \, w(\xi) \]

where \(\mathcal{E}\) is the set of all elements in the domain, \(\mathcal{Q}\) is the set of all quadrature points in the domain, \(J(\xi)\) is the Jacobian matrix that maps the gradient in physical space to the reference element, and \(\boldsymbol{\sigma}(\xi)\) and \(\boldsymbol{\varepsilon}(\xi)\) are the stress and strain tensors evaluated at the quadrature point \(\xi\). The number and the values of quadrature points ( along with the weights \(w(\xi)\) associated with them) depends on the type of gauss-quadrature rule used, check Gaussian quadrature for more details.

Note

In this course, we will use linear triangular elements for the discretization with one quadrature point (\(\xi\)) located at the centroid of the element.

In order to compute the above integration of the strain energy over the domain, we need to perform 2 mathematical operations:

  • Evaluate strain tensors from the nodal displacements at the quadrature points.
  • Evaluate stress tensors from the strain tensors.

We will now discuss how to perform these operations in a finite element method.

Evaluation of nodal values at quadrature points

The evaluation of a nodal values at a quadrature point is given by:

\[ \boldsymbol{u}(\xi) = \boldsymbol{N}(\xi) \hat{\boldsymbol{u}}_e \]

where \(\boldsymbol{N}(\xi)\) is the shape function matrix and \(\hat{\boldsymbol{u}}_e\) is the vector of nodal displacements for the element \(e\).

For a linear triangle element, the shape function matrix is given as:

\[ \boldsymbol{N}(\xi) = \begin{bmatrix} \underbrace{1 - \xi_1 - \xi_2}_{N_1}, & \underbrace{\xi_1}_{N_2}, & \underbrace{\xi_2}_{N_3} \end{bmatrix} \]

where \(\xi_1\) and \(\xi_2\) are the natural coordinates of the quadrature point. The shape functions take value 1 at the node corresponding to the shape function and 0 at the other nodes. The figure below shows the shape functions for a linear triangle element.

The displacement values at a quadrature point is given as:

\[ \boldsymbol{u}(\xi_1, \xi_2) = N_1(\xi_1, \xi_2) \hat{\boldsymbol{u}}_e + N_2(\xi_1, \xi_2) \hat{\boldsymbol{u}}_e + N_3(\xi_1, \xi_2) \hat{\boldsymbol{u}}_e. \]

Gradient of nodal values and strain tensor

The strain tensor is given as

\[ \boldsymbol{\epsilon}(x) = \dfrac{1}{2} \left( \nabla \boldsymbol{u}(x) + \nabla \boldsymbol{u}(x)^T \right) \]

where \(\nabla \boldsymbol{u}(x)=\left( \frac{\partial \boldsymbol{u}(x)}{\partial x_1}, \frac{\partial \boldsymbol{u}(x)}{\partial x_2}, \frac{\partial \boldsymbol{u}(x)}{\partial x_3} \right)\) is the gradient of the nodal displacements evaluated at point \(\boldsymbol{x}\). As shown in previous section, the displacement values can be evaluated at quadrature points (\(\boldsymbol{\xi}\)) using the shape function matrix. Using that relation, we can write the gradient of the nodal displacements at the quadrature point as,

\[ \nabla \boldsymbol{u}(x) = \nabla (\boldsymbol{N}(\boldsymbol{\xi}) \hat{\boldsymbol{u}}_e) = \frac{\partial}{\partial \boldsymbol{x}} (\boldsymbol{N}(\boldsymbol{\xi}) \hat{\boldsymbol{u}}_e), \]

where \(\nabla \boldsymbol{N}(\boldsymbol{\xi})\) is the gradient of the shape function matrix evaluated at the quadrature point \(\boldsymbol{\xi}\) and \(\nabla \boldsymbol{u}_e\) is the gradient of the nodal displacements for the element \(e\). We can use chain rule to express the gradient (\(\nabla = \frac{\partial}{\partial \boldsymbol{x}}\)) in terms of the gradient in the reference coordinate system (\(\boldsymbol{\xi}\)) as:

\[ \nabla \boldsymbol{u}(\boldsymbol{x}) =\frac{\partial}{\partial \boldsymbol{x}} (\boldsymbol{N}(\boldsymbol{\xi}) \hat{\boldsymbol{u}}_e) = \frac{\partial \boldsymbol{N}(\boldsymbol{\xi}) \hat{\boldsymbol{u}}_e}{\partial \boldsymbol{\xi}}\frac{\partial \boldsymbol{\xi}}{\partial \boldsymbol{x}}, \]

where \(\frac{\partial \boldsymbol{\xi}}{\partial \boldsymbol{x}}\) is the gradient of the reference coordinate system in the physical coordinate system. We know that a Jacobian matrix is the gradient of the the mapping function \(\boldsymbol{\xi} \to \boldsymbol{x}\) in reference coordinate system. i.e. \[ \boldsymbol{J} = \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{\xi}} = \begin{bmatrix} \frac{\partial x_1}{\partial \xi_1} & \frac{\partial x_1}{\partial \xi_2} \\ \frac{\partial x_2}{\partial \xi_1} & \frac{\partial x_2}{\partial \xi_2} \end{bmatrix} \]

Therefore, we can use this relation to express inverse gradient \(\frac{\partial \boldsymbol{\xi}}{\partial \boldsymbol{x}}\) simply as the inverse of the Jacobian matrix,

\[ \frac{\partial \boldsymbol{\xi}}{\partial \boldsymbol{x}} = \boldsymbol{J}^{-1}. \]

Using this relations, we can express the gradient of the nodal displacements in physical space as,

\[ \nabla \boldsymbol{u}(\boldsymbol{x}) =\frac{\partial \boldsymbol{N}(\boldsymbol{\xi}) \hat{\boldsymbol{u}}_e}{\partial \boldsymbol{\xi}}\boldsymbol{J}^{-1} \]

Since \(\hat{\boldsymbol{u}}_e\) is a vector of nodal displacements in the physical space and doesnot depends on the reference coordinate system, we can take it outside the gradient operator as a constant vector. The gradient of displacement values is thus given as,

\[ \nabla \boldsymbol{u}(\boldsymbol{x}) =\frac{\partial \boldsymbol{N}(\boldsymbol{\xi})}{\partial \boldsymbol{\xi}}\boldsymbol{J}^{-1} \cdot{} \hat{\boldsymbol{u}}_e \]

Key Idea

The gradient of any quantity \(\phi(\boldsymbol{x})\) in physical space is thus given as, \[ \nabla \phi(\boldsymbol{x}) =\frac{\partial \boldsymbol{N}(\boldsymbol{\xi})}{\partial \boldsymbol{\xi}}\boldsymbol{J}^{-1} \cdot{} \phi_{e} \]

Once we have the gradient of the displacement values, we can evaluate the strain tensor at the quadrature point by taking the transpose of the gradient of the displacement values and adding it to the gradient of the displacement values, as given below:

\[ \boldsymbol{\epsilon}(\boldsymbol{x}) = \dfrac{1}{2} \left( \nabla \boldsymbol{u}(\boldsymbol{x}) + \nabla \boldsymbol{u}(\boldsymbol{x})^T \right) \]

We can now use the above operations to integrate the strain energy over the domain.

Note

In this course, we will use tatva library to perform these 3 operations. The basic usage of how to perform these operations is given Appendix A — Basics of tatva. Please check that section for more details.

Using differentiation to compute the internal forces and stiffness matrix

Once we have the total internal energy \(\Psi_\text{int}(\boldsymbol{u})\) (or the total strain energy \(\Psi_\text{e}(\boldsymbol{u})\) as shown above), we can differentiate it with respect to the nodal displacements \(\hat{\boldsymbol{u}}\) to compute the internal forces \(\hat{\boldsymbol{f}}_\text{int}\) and the stiffness matrix \(\hat{\mathbf{K}}\). Therefore,

\[ \hat{\boldsymbol{f}}_\text{int} = \dfrac{\partial \Psi_\text{e}}{\partial \hat{\boldsymbol{u}}},\quad \text{and} \quad \hat{\mathbf{K}} = \dfrac{\partial \hat{\boldsymbol{f}}_\text{int}}{\partial \hat{\boldsymbol{u}}} \]

Since we define computer code to compute the total internal energy \(\Psi_\text{int}(\boldsymbol{u})\) in tatva (see Appendix A — Basics of tatva), we cannot differentiate these python functions analytically. We will use the numerical technique called automatic differentiation to differentiate these python functions. For this course, we will use the JAX library to perform the automatic differentiation. We dont need to know the details of how automatic differentiation works. What we need to know is how to use it compute the derivatives of the python functions.

Note

In this course, we will use JAX library to compute the internal forces and the stiffness matrix from the total internal energy \(\Psi_\text{int}(\boldsymbol{u})\). The basic usage of how to perform these operations is given Appendix B — Basics of JAX. Please check that section for more details.

This way of directly differentiating the internal energy functional to compute the internal forces and the stiffness matrix is different from the standard Finite Element Method. In FEM, we compute the local stiffness matrix at element level (using shape functions and its derivaties) and then assemble this local stiffness matrix into a global stiffness matrix based on the connectivity of the mesh. Whereas here we differetiate the scalar internal energy with respect to the nodal displacement values and this way we directly get out internal force vector. Differentiating this nodal force vector again with nodal displacement gives us our global stiffness matrix. So we do not have to go through the manual assemble process.

Core Concepts

To solve a problem, in this course, we will follow the following steps:

  • Discretize the domain into elements and nodes.
  • Compute the total internal energy \(\Psi_\text{int}(\boldsymbol{u})\) by integrating the strain energy over the domain.
  • Differentiate the total internal energy with respect to the nodal displacements \(\hat{\boldsymbol{u}}\) to compute the internal forces \(\hat{\boldsymbol{f}}_\text{int}\) and the stiffness matrix \(\hat{\mathbf{K}}\).
  • Solve the equilibrium equation \(\hat{\boldsymbol{f}}_\text{int} - \hat{\boldsymbol{f}}_\text{ext} = \hat{\boldsymbol{0}}\) to find the nodal displacements \(\hat{\boldsymbol{u}}\).

References

Bathe, K. J. 2006. Finite Element Procedures. Prentice Hall.
Hughes, Thomas J R. 2000. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover.
Zienkiewicz, O. C., R. L. Taylor, and J. Z. Zhu. 2005. The Finite Element Method: Its Basis and Fundamentals. Butterworth-Heinemann.