Introduction
The goal of this tutorial is to build step by step a rather simple model using MoReFEM: elastic deformation of a bar upon which a force is applied; time step is constant and volumic mass is the same throughout the bar.
The library was first built with a use in macOS and XCode in mind, but it has been extended since and the procedure should guide you to make the project work in both macOS and Linux environments.
Description of the mechanical problem
In this tutorial we are solving a simple 2D problem: the deformation of a bar within the theory of linear elasticity using a finite element method. Here is a quick description of the equations depicting the problem and the resolution method that we will implement step by step; our model will in fact be generic enough to work with 3D cases.
The strong formulation of the equilibrium equations of a dynamic mechanical problem on a fixed domain $ \Omega_t $ reads:
\begin{cases}
\displaystyle \rho(\underline{f}(\underline{x}) - \underline{\ddot{y}}(\underline{x})) + \underline{\nabla}_{\underline{x}} \cdot \underline{\underline{\sigma}}(\underline{x}) = 0 \quad \text{in $\Omega_t$} \\
\displaystyle \underline{\underline{\sigma}}(\underline{x}) \cdot \underline{n} = \underline{g}(\underline{x}) \quad \text{on } \Gamma^N \\
\displaystyle \underline{y}(\underline{x}) = \underline{0} \quad \text{on } \, \Gamma^D \\
\displaystyle \underline{\underline{\sigma}}^T(\underline{x}) = \underline{\underline{\sigma}}(\underline{x})
\end{cases}
The domain $\Omega_t$ over which the solid is defined has a fixed boundary $\Gamma^D$ where Dirichlet conditions are applied (here we impose a zero displacement condition) and another boundary $\Gamma^N$ where Neumann conditions are applied (for instance external surfacic loadings in 3D). $\underline{\underline{\sigma}}$ is the Cauchy stress tensor, $\rho$ the volumic mass in the current configuration, $\underline{f}$ represents the applied external forces, $\underline{\ddot{y}}$ is the acceleration, $\underline{n}$ is the outward-pointing normal of the boundary $\Gamma^N $ and $\underline{g}$ the external surfacic load.
The weak formulation of the fundamental law of dynamics in the deformed configuration reads:
\forall \underline{y}^* \in \mathcal{V}, \quad \int_{\Omega_t} \rho(\underline{f}-\underline{\ddot{y})} \cdot \underline{y}^* \textrm{d}\Omega_t + \int_{\Omega_t} \underline{\underline{\sigma}}:\underline{\nabla}_{\underline{x}}\underline{y}^* \textrm{d}\Omega_t = \int_{\Gamma^N} \underline{g}\cdot \underline{y}^* \text{d}S
In order to solve this equation without having to deal with the changes of positions due to the fact that each variable is expressed in the reference configuration, we can rewrite this equation within a Total Lagrangian formalism. Here we will only consider a surfacic loading $ \underline{g}(\underline{x})$ without a volumic contribution (namely $ \underline{f}(\underline{x}) = \underline{0} $). Without going into the details of the chain rules involved, the variational formulation of the fundamental law of dynamics in total Lagrangian formalism within elasticity theory reads:
\forall \underline{y}^* \in \mathcal{V}, \quad \displaystyle \int_{\Omega _{0}}^{} \rho_0 \underline{\ddot{y}} \cdot \underline{y}^{*}\textrm{d}\Omega_0 + \int_{\Omega_{0}}^{} \left( \underline{\underline{\underline{\underline{\text{A}}}}} : \underline{\underline{\varepsilon}}(\underline{y}) \right) : \underline{\underline{\varepsilon}}(\underline{y}^*) \textrm{d}\Omega_0
= \int_{\Gamma_0^N} \underline{g}_0 \cdot \underline{y}^* \text{d}S_0
Where $ \underline{\underline{\varepsilon}}(\underline{y}) = \underline{\underline{\nabla}} \, \underline{y} + (\underline{\underline{\nabla}} \, \underline{y})^T $ is the linearized strain tensor and the consitutive behaviour law $ \underline{\underline{\sigma}} = \underline{\underline{\underline{\underline{\text{A}}}}} : \underline{\underline{\varepsilon}}(\underline{y}) $.
Here we are only solving for one unknown: the displacement $ \underline{y}$. In order to do so, we will place ourselves within the plane strain hypothesis, namely $ \varepsilon_{33} = 0$ which implies:
\underline{\underline{\sigma}} = \begin{pmatrix}
\sigma_{xx} & \sigma_{xy} & 0 & \\
\sigma_{xy} & \sigma_{yy} & 0 & \\
0 & 0 & \sigma_{zz} &
\end{pmatrix} \quad \text{with} \quad \displaystyle \sigma_{zz} = \frac{\lambda}{2(\mu+\lambda)}(\sigma_{xx}+\sigma_{yy})
where $\lambda$ and $\mu$ are the Lamé coefficients of the elastic solid.
Thus the constitutive behaviour law in 2D reads, using engineering notation:
\begin{pmatrix}
\sigma_{xx}(\underline{y}) & \\
\sigma_{yy}(\underline{y}) & \\
\sigma_{xy}(\underline{y}) &
\end{pmatrix} = \begin{pmatrix}
\lambda + 2\mu & \lambda & 0 & \\
\lambda & \lambda + 2\mu & 0 & \\
0 & 0 & \mu &
\end{pmatrix}\begin{pmatrix}
\varepsilon_{xx}({\underline{y}}) &\\
\varepsilon_{yy}({\underline{y}}) &\\
2\varepsilon_{xy}({\underline{y}}) &
\end{pmatrix} = \underline{\underline{\hat{\text{A}}}} \cdot \underline{\hat{\varepsilon}} (\underline{y})
Finally, we are left with the following set of equations and boundary conditions:
\begin{cases}
\displaystyle \forall \underline{y}^* \in \mathcal{V}, \quad \displaystyle \int_{\Omega _{0}}^{} \rho_0 \, \underline{y}^{*} \cdot \underline{\ddot{y}} \, \textrm{d}\Omega_0 + \int_{\Omega_{0}}^{} \underline{\hat{\varepsilon}}(\underline{y}^*)^T \cdot \underline{\underline{\hat{\text{A}}}} \cdot \underline{\hat{\varepsilon}}(\underline{y}) \, \textrm{d}\Omega_0
= \int_{\Gamma_0^N} \underline{y}^* \cdot \underline{g}_0 \, \text{d}S_0 \\
\displaystyle \underline{y}(\underline{x}) = \underline{0} \quad \text{on } \, \Gamma^D
\end{cases}
This bilinear system (with respect to the virtual displacement field $ \underline{y}^* $ and $\underline{y}$) will be solved using the finite element method.
Resolution of the linear system
Spatial discretization
In order to solve our previously defined bilinear system (with respect to $\underline{y}$ and $\underline{y}^*$), we will be using the standard Galerkin method. It consists of approximating the function of interest (the displacement field in our case) by a finite sum of known shape functions (polynomials usually) $\phi_k(\underline{\xi})$ weighted by unknown coefficients $y_{jk} $ where $k \in [1, \, N + 1] $, $ N $ being the order of the shape functions used. In 2D, the discretization of the displacement field gives:
y_j = \sum_{k=1}^{N+1} y_{jk} \phi_k (\underline{\xi}) \quad j \in [x,y]
\underline{y}(\underline{x})_{2D} = \begin{pmatrix}
y_{x} & \\
y_{y} &
\end{pmatrix} = \begin{pmatrix}
\phi_1 & \dots & \phi_{N+1} & 0 & \dots & 0 & \\
0 & \dots & 0 & \phi_1 & \dots & \phi_{N+1} &
\end{pmatrix} \begin{pmatrix} y_{x1} & \\
\vdots & \\
y_{xN+1} & \\
y_{y1} & \\
\vdots & \\
y_{yN+1} &
\end{pmatrix} = \underline{\underline{\mathbb{N}}} \cdot \underline{\mathbb{U}}_h
To discretize the linearized strain tensor $\underline{\underline{\varepsilon}} $ we can note that:
\hat{\underline{\varepsilon}}(\underline{y}) = \begin{pmatrix} \varepsilon_{xx} \\
\varepsilon_{yy} & \\
2\varepsilon_{xy} &
\end{pmatrix} = \begin{pmatrix}
\partial_x y_x & \\
\partial_y y_y & \\
\partial_y y_x + \partial_x y_y &
\end{pmatrix} =
\begin{pmatrix}
1 & 0 & 0 & 0 & \\
0 & 0 & 0 & 1 & \\
0 & 1 & 1 & 0 &
\end{pmatrix} \begin{pmatrix}
\partial_x y_x & \\
\partial_y y_x & \\
\partial_x y_y & \\
\partial_y y_y & \\
\end{pmatrix}
We can then apply the spatial discretization in a similar fashion using the derivatives of the shape functions:
\hat{\underline{\varepsilon}}(\underline{y}) = \begin{pmatrix} \varepsilon_{xx} \\
\varepsilon_{yy} & \\
2\varepsilon_{xy} &
\end{pmatrix} =
\begin{pmatrix}
1 & 0 & 0 & 0 & \\
0 & 0 & 0 & 1 & \\
0 & 1 & 1 & 0 &
\end{pmatrix} \begin{pmatrix}
\partial_x \phi_1 & \dots & \partial_x \phi_{N+1} & 0 & \dots & 0 & \\
\partial_y \phi_1 & \dots & \partial_y \phi_{N+1} & 0 & \dots & 0 & \\
0 & \dots & 0 & \partial_x \phi_1 & \dots & \partial_x \phi_{N+1} & \\
0 & \dots & 0 & \partial_y \phi_1 & \dots & \partial_y \phi_{N+1} &
\end{pmatrix} \begin{pmatrix} y_{x1} & \\
\vdots & \\
y_{xN+1} & \\
y_{y1} & \\
\vdots & \\
y_{yN+1} &
\end{pmatrix} = \underline{\underline{\mathbb{B}}} \cdot \underline{\mathbb{U}}_h
Plugging these discretized forms into our equilibrium equation gives:
\forall \, \underline{\mathbb{U}}^*_h \in \mathcal{V}_h, \quad \displaystyle \int_{\Omega_0} \rho_0 \underline{\mathbb{U}}^{*T}_h \underline{\underline{\mathbb{N}}} ^T \cdot \underline{\underline{\mathbb{N}}} \cdot \underline{\dot{\mathbb{V}}}_h \, \text{d}\Omega_0 + \int_{\Omega_{0}}^{} \underline{\mathbb{U}}^{*T}_h \cdot \underline{\underline{\mathbb{B}}}^T \cdot \underline{\underline{\hat{\text{A}}}} \cdot \underline{\underline{\mathbb{B}}} \cdot \underline{\mathbb{U}}_h \, \textrm{d}\Omega_0
= \int_{\Gamma_0^N} \underline{\mathbb{U}}^{*T} \cdot \underline{\underline{\mathbb{N}}}^T \cdot \underline{g}_0 \, \text{d}S_0
where $\underline{
{\mathbb{V}}}_h $ is time derivative of the unknown coefficients $@dot{y} _{jk}$ relative to the velocity field (which is itself the time derivative of the unknown weighting coefficients of the displacement field).
