Mechanical Finite Element Simulations Explained


Finite element (FE) simulations are powerful predictive tools, and are essential to engineering design, simplifying design processes, and reducing prototyping costs. However, unfortunately and quite frequently, these get performed in a pure "black-box" fashion, by which the whole process of obtaining a solution is not fully understood by the user. While this is sometimes acceptable, as the corresponding FE package usually takes care of most details, it is nonetheless useful to understand what is going on behind the scenes. In this blog we explore the basic equations solved in FE applied to deformable continuum solid mechanics, and how these are treated from a computational perspective.

  • Strong form

    Consider a generic body in three-dimensional space, \Omega , and its boundary \partial\Omega . In traditional continuum solid mechanics, the balance of momentum equations are solved, which, under the assumption of no inertia (i.e. the system can then be considered as static or quasi-static) and small strains, are 123

    \begin{align}\text{div} \pmb{\sigma} &= \mathbf{b} \text{ in } \Omega, \\ \pmb{\sigma} \mathbf{n} &= \mathbf{t} \, \text{on }\partial\Omega,\end{align}

  • where \text{div} is the divergence operator, \pmb{\sigma} is the stress, \mathbf{b} are the volumetric body forces, \mathbf{n} is the normal to \partial\Omega , and \mathbf{t} is the surface traction. Equations 1 and 2 can be easily understood by integrating the first over the volume of the body and then applying the divergence theorem, leading to
  • \begin{align*}\int_\Omega \text{div} \, \pmb{\sigma} \, \text{d}V \stackrel{\substack{\text{divergence} \\ \text{theorem}}}{=} \int_{\partial\Omega} \pmb{\sigma} \mathbf{n} \, \text{d}A = \int_{\partial\Omega} \textbf{t} \, \text{d}A = \int_{\Omega} \mathbf{b} \,\text{d}V,\end{align*}


    which basically states that Eq. 1 is equivalent to a force balance over the considered body \Omega . Eq. 2 can be further decomposed into

    \begin{align}\pmb{\sigma} \mathbf{n}= \overline{\mathbf{t}} \text{ on } \partial\Omega_N, \\ \pmb{\sigma} \mathbf{n} = \tilde{\mathbf{t}} \text{ on } \partial\Omega_D, \end{align}

    where \bar{\mathbf{t}} is the prescribed traction at the portion of the boundary where Neumann boundary conditions are applied, \partial\Omega_N ; \tilde{\mathbf{t}} is the traction that arises as a reaction at the portion of the boundary where Dirichlet boundary conditions are applied \partial\Omega_D , and also \partial\Omega = \partial\Omega_N \cup \partial\Omega_D and \partial\Omega_N \cap \partial\Omega_D = \varnothing . As a brief reminder, Neumann boundary conditions prescribe derivatives of the solution field, whilst Dirichlet boundary conditions prescribes the solution field itself.

    The term "strong form" stated in the section title refers to the fact that Eq. 1 is enforced point-wise, locally, at each point of the considered body \Omega.

    Weak form

    Eqs. 1 and 2 are not usable from a practical perspective because their solution on complicated geometries and/or with sophisticated material models cannot, in general, be found analytically, and thus we resort to numerical solutions instead, such as the finite element method.

  • In order to apply a finite element formulation, the weak form of Eqs. 1 and 2 must be used, which in the context of solid mechanics results in the principle of virtual work. The term "weak form" here refers to the fact that the equations are not enforced locally, but rather in an average sense, i.e. in an integral sense. Therefore, Eq. 1 is first integrated over the volume of the considered body \Omega and then contracted with an arbitrary virtual displacement \pmb{\eta}, which is zero at boundaries where Dirichlet boundary conditions are applied \partial\Omega_D. This leads to

    \begin{align}\int_\Omega \left( \text{div} \, \pmb{\sigma} - \mathbf{b} \right) \cdot \pmb{\eta} \, \text{d}V = 0.\end{align}

    By using integration by parts, Eq. 5 becomes

    \begin{align*}\int_\Omega \left[ \pmb{\sigma} : \nabla\pmb{\eta} - \text{div} \left( \pmb{\sigma} \pmb{\eta} \right) -\mathbf{b} \right] \text{d}V = 0,\end{align*}

    and finally, applying the divergence theorem leads to the final expression of the principle of virtual work,

    \begin{align*}\int_\Omega \pmb{\sigma} : \nabla\pmb{\eta} \, \text{d}V - \int_\Omega \mathbf{b} \cdot \pmb{\eta} \, \text{d}V - \int_{\partial\Omega} \mathbf{t} \cdot \pmb{\eta} \, \text{d}A = 0,\end{align*}

    where \cdot and : are, respectively, single and double dot products. Finally, by applying Eqs. 3 and 4, and considering that \pmb{\eta} = 0 at \partial\Omega_D, the final expression of the principle of virtual work is obtained as

    \begin{align*}\int_\Omega \pmb{\sigma} : \nabla\pmb{\eta} \, \text{d}V - \int_\Omega \mathbf{b} \cdot \pmb{\eta} \,\text{d}V - \int_{\partial\Omega_N} \bar{\mathbf{t}} \cdot \pmb{\eta} \, \text{d}A = 0,\end{align*}

    Finite element discretization

  • Figure 1 shows the discretization process, in which a complex geometry is converted into a mesh of n^e finite elements and n_n nodes, where each element e is composed of n^e_n nodes.

    Fig. 1. (Left) Original Geometry, (Right) Discretized geometry with finite elements, each triangle here represents an element e occupying a domain \Omega^e


    Using this finite element discretization, the discretized geometry \Omega_h is defined as

    \begin{align*}\Omega_h = \bigcup^{n^e}_{e = 1} \Omega^e. \end{align*}

    By using the appropriate finite element interpolation functions,

    \begin{align*}\mathbf{u} \left( \mathbf{X}, t \right) = \sum^{n^e_n}_{a = 1} N_a \left( \mathbf{X} \right) \mathbf{u}_a \left( t \right), \\ \pmb{\eta} \left( \mathbf{X}, t \right) = \sum^{n^e_n}_{a = 1} N_a \left( \mathbf{X} \right) \pmb{\eta}_a \left( t \right),\end{align*}

    where \mathbf{X} are the undeformed coordinates of the body \Omega, t is time, N_a is the finite element interpolation function corresponding to node a of the element e, and \mathbf{u}_a is the displacement at node a ; therefore, the contribution of a single node a of element e to the discretized principle of virtual work is obtained as

    \begin{align*}\int_{\Omega^{e}} \left( \pmb{\eta}_a \cdot \pmb{\sigma} \nabla N_a \right) \text{d}V - \int_{\Omega^{e}} \left( \mathbf{b} \cdot N_a \pmb{\eta}_a \right) \text{d}V - \int_{\partial\Omega^{e}_N} \left( \bar{\mathbf{t}} \cdot N_a \pmb{\eta}_a \right) \text{d}A,\end{align*}

    where \partial\Omega^e is the boundary of element e. The contribution from all elements e containing the node a ( a \in e) to the discretized principle of virtual work is

    \begin{align*} \sum^{n_e}_{\substack{e = 1 \\ a \in e}} & \left[ \int_{\Omega^{e}} \left( \pmb{\eta}_a \cdot \pmb{\sigma} \nabla N_a \right) \text{d}V - \int_{\Omega^{e}} \left( \mathbf{b} \cdot N_a \pmb{\eta}_a \right) \text{d}V - \int_{\partial\Omega^{e}_N} \left( \bar{\mathbf{t}} \cdot N_a \pmb{\eta}_a \right) \text{d}A \right] \\ &= \pmb{\eta}_a \cdot \sum_{\substack{e = 1 \\ a \in e}}^{n_e} \left[ \int_{\Omega^{e}} \pmb{\sigma} \nabla N_a \, \text{d}V -\int_{\Omega^{e}} N_a \mathbf{b} \, \text{d}V - \int_{\partial\Omega^{e}_N} N_a \bar{\mathbf{t}} \, \text{d}A \right].\end{align*}

  • Finally, by summing the contribution from all the nodes in the mesh, the final expression for the discretized principle of virtual work can be obtained as

    \begin{align}\sum^{n_n}_{a = 1} \pmb{\eta}_a \cdot \sum^{n_e}_{\substack{e = 1 \\ a \in e}} \left[ \int_{\Omega^{e}} \pmb{\sigma} \nabla N_a \, \text{d}V - \int_{\Omega^{e}} N_a \mathbf{b} \, \text{d}V -\int_{\partial\Omega^{e}_N} N_a \bar{\mathbf{t}} \, \text{d}A \right] = 0.\end{align}

  • Equation 6 is equivalent to

    \begin{align}\mathbf{H}^\text{T} \left(\mathbf{F}_\text{int} - \mathbf{F}_\text{ext} \right) = 0,\end{align}

    by rearranging the terms as

    \begin{align*} \mathbf{H} &=\begin{Bmatrix} \pmb{\eta}_1 & \pmb{\eta}_2 & \dots & \pmb{\eta}_{n_n} \end{Bmatrix}^\text{T}, \\ \mathbf{F}_\text{int} &= \begin{Bmatrix} \mathbf{F}_{\text{int},1} & \mathbf{F}_{\text{int},2} & \dots & \mathbf{F}_{\text{int},n_n} \end{Bmatrix}^\text{T}, \\ \mathbf{F}_\text{ext} &=\begin{Bmatrix}\mathbf{F}_{\text{ext},1} &\mathbf{F}_{\text{ext},2} &\dots &\mathbf{F}_{\text{ext},n_n}\end{Bmatrix}^\text{T},\end{align*}

    where, respectively, the nodal internal and external forces are

    \begin{align*} & \mathbf{F}_{\text{int},a} = \sum^{n_e}_{\substack{e = 1 \\ a \in e}} \int_{\Omega^e} \pmb{\sigma} \nabla N_a \, \text{d}V, \\ & \mathbf{F}_{\text{ext},a} = \sum^{n_e}_{\substack{e = 1 \\ a \in e}} \left[ \int_{\Omega^e} N_a \mathbf{b} \, \text{d}V + \int_{\partial\Omega^{e}_N} N_a \bar{\mathbf{t}} \, \text{d}A \right]. \end{align*}

    The nodal internal forces can be recast in the familiar matrix-vector notation, such that

    \begin{align*}\mathbf{F}_{\text{int},a} = \sum^{n_e}_{\substack{e = 1 \\ a \in e }} \int_{\Omega^e} (\mathbf{B}^\text{T}_a \{{\pmb{\sigma}}\}) \text{d}V,\end{align*}


    \begin{align*} \mathbf{B}_a = \begin{bmatrix}\frac{\partial N_a}{\partial X_1} & 0 & 0 \\ 0 & \frac{\partial N_a}{\partial X_2} & 0 \\ 0 & 0 & \frac{\partial N_a}{\partial X_3} \\ \frac{\partial N_a}{\partial X_2} & \frac{\partial N_a}{\partial X_1} & 0 \\ 0 & \frac{\partial N_a}{\partial X_3} & \frac{\partial N_a}{\partial X_2} \\ \frac{\partial N_a}{\partial X_3} & 0 & \frac{\partial N_a}{\partial X_1} \ \end{bmatrix}\end{align*}


    \begin{align*}{\pmb{\sigma}} = \begin{Bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{23} & \sigma_{13} \end{Bmatrix}^\text{T}.\end{align*}

    \mathbf{B}_a is the symmetric gradient operator such that the vector projection of the infinitesimal strain \pmb{\varepsilon} can be obtained at an integration point as

    \begin{align*}\pmb{\varepsilon} = \nabla^\text{sym} \mathbf{u} = \frac{1}{2} \left[ \frac{\partial \mathbf{u}}{\partial \mathbf{X}} + \left( \frac{\partial \mathbf{u}}{\partial \mathbf{X}} \right)^\text{T} \right] = \frac{1}{2} \sum^{n^e_n}_{a = 1} \left[ \mathbf{u}_a \otimes \frac{\partial N_a}{\partial \mathbf{X}} + \frac{\partial N_a}{\partial \mathbf{X}} \otimes \mathbf{u}_a \right],\end{align*}

    and therefore

    \begin{align*}\{{\pmb{\varepsilon}}\} &= \begin{Bmatrix} \varepsilon_{11} & \varepsilon_{22} & \varepsilon_{33} & 2\varepsilon_{12} & 2\varepsilon_{23} & 2\varepsilon_{13} \end{Bmatrix}^\text{T} \\ &= \sum^{n^e_n}_{a = 1} \Bigg\lbrace u_1 \frac{\partial N_a}{X_1} \ u_2 \frac{\partial N_a}{X_2} \ u_3 \frac{\partial N_a}{X_3} \ \cdots \\ &\quad \quad \quad \cdots \ u_1 \frac{\partial N_a}{\partial X_2} + u_2 \frac{\partial N_a}{\partial X_1} \ u_2 \frac{\partial N_a}{\partial X_3} + u_3 \frac{\partial N_a}{\partial X_2} \ u_1 \frac{\partial N_a}{\partial X_3} + u_3 \frac{\partial N_a}{\partial X_1} \Bigg\rbrace ^\text{T} \\ &= \sum^{n^e_n}_{a = 1} \mathbf{B}_a \begin{Bmatrix} u_{a,1} & u_{a,2} & u_{a,3} \end{Bmatrix}^\text{T}. \end{align*}

  • Equation 7 must hold for all arbitrary virtual displacements and thus it can be recast in the residual form

    \begin{align} \mathbf{r} = \sum^{n_n}_{a = 1} \sum^{n_e}_{\substack{e = 1 \\ a \in e}} \left[ {\int_{\Omega^e} \mathbf{B}^\text{T}_a \{{\pmb{\sigma}}\} \, \text{d}V} - \left( \int_{\Omega^e} N_a \mathbf{b} \ \text{d}V + \int_{\partial\Omega^{e}_N} N_a \bar{\mathbf{t}} \, \text{d}A \right)\right] = \mathbf{0} \end{align}

    which is well-recognized. Equation 8 is generally nonlinear on the displacements \mathbf{u}, which is the sought function; from here onwards and for the sake of simplicity, this equation is written without the summations by rearranging the appropriate terms in matrix-vector form as

    \begin{align*} \mathbf{R} &= \mathbf{F}_\text{int}- \mathbf{F}_\text{ext} \\ &= \bigg\lbrace \int_{\Omega_h} \mathbf{B}^\text{T} \{{\pmb{\sigma}}\} \, \text{d}V \bigg\rbrace - \underbrace{\bigg\lbrace \int_{\Omega_h} \mathbf{N}^\text{T} \mathbf{b} \, \text{d}V \bigg\rbrace}_{\text{internal virtual work}} + \underbrace{ \bigg\lbrace \int_{\partial\Omega_h,N } \mathbf{N}^\text{T} \overline{\mathbf{t}} \, \text{d}A \bigg\rbrace}_{\text{external virtual work}} = \mathbf{0}, \end{align*}


    \begin{align*} \mathbf{R} &= \begin{Bmatrix}\mathbf{F}_{\text{int}, 1} - \mathbf{F}_{\text{ext}, 1} & \mathbf{F}_{\text{int}, 2} - \mathbf{F}_{\text{ext}, 2} & \dots & \mathbf{F}_{\text{int}, n_n} - \mathbf{F}_{\text{ext}, n_n}\end{Bmatrix}^\text{T}, \\ \mathbf{B} &=\begin{bmatrix}\mathbf{B}_1 & \mathbf{B}_2 & \dots & \mathbf{B}_{n^e_n}\end{bmatrix}, \\ \mathbf{N} &=\begin{bmatrix}N_1 & 0 & 0 & N_2 & 0 & 0 & \dots & N_{n^e_n} & 0 & 0 \\0 & N_1 & 0 & 0 & N_2 & 0 & \dots & 0 & N_{n^e_n} & 0 \\0 & 0 & N_1 & 0 & 0 & N_2 & \dots & 0 & 0 & N_{n^e_n}\end{bmatrix}.\end{align*}

  • Coming up next!

    The next blog will continue from these equations and ellaborate on how they become either linear or nonlinear, and also develop one of the most common types of nonlinearity, i.e. material nonlinearity and in particular, plasticity.


    • T. Belytschko et al. Nonlinear Finite Elements for Continua and Structures. Wiley, 2013. ISBN: 9781118700082.
    • J. Bonet and R.D. Wood. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press, 2008. ISBN: 9780511755446.
    • E.A. de Souza Neta, D. Peric, and D.R.J. Owen. Computational Methods for Plasticity: Theory and Applications. Wiley, 2011. ISBN: 9781119964544.
  • Francesc Levrero-Florencio
    Francesc Levrero-Florencio

    I am a mathematical modeler with expertise in high performance computing, who is able to design and improve current prediction models through the implementation of novel and efficient numerical methods.

    Discover how customers like you found success by
    leaving traditional engineering simulation behind

    Try OnScale following
    our simulation guides

    Simulate Now

    Discuss your engineering
    applications with us

    Request a Demo