Introduction
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}
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.
\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
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*}
\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}
\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*}
where
\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*}
and
\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*}
\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*}
where
\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.
References
- 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.