The linear elasticity equation with Coefficient Form PDEs toolbox
1. Introduction
This notebook deals with the finite element approximation of deformable continuum mechanics problems in three space dimensions. The domain \(\Omega \subset \mathbb{R}^3\) represents a deformable medium, initially at rest and to which an external loading \(f: \Omega \rightarrow \mathbb{R}^3\) is applied. The objective is to determine the displacement field \(u: \Omega \rightarrow \mathbb{R}^3\) induced by this loading once the medium has reached equilibrium. It is assumed that the deformations are small enough to be modeled in the framework of linear elasticity. For simplicity, it is also assumed that the medium is isotropic.
We note \(\sigma: \Omega \rightarrow \mathbb{R}^{3,3}\) the stress field. The equilibrium condition is written as follows
Let us note \(\varepsilon(u): \Omega \rightarrow \mathbb{R}^{3,3}\) the linearized deformation field defined by
\(\nabla u\) is the field with values in \(\mathbb{R}^{3,3}\) of components \((\nabla u)_{i j}=\partial_j u_i, i, j \in\{1,2,3}\). 
In the framework of linear elasticity, the stress field is expressed as a function of the linearized strain field as
where \(\lambda\) and \(\mu\) are the Lamé coefficients and \(\mathcal{I}_3\) the identity matrix in \(\mathb{R}^{3,3}\). Using the definition of the deformation tensor, we get
The Lamé coefficients are phenomenological coefficients which, for thermodynamic reasons, are constrained by the relations \(\mu>0\) and \(\lambda+\frac{2}{3} \mu \geq 0\). For simplicity, we assume that \(\lambda \geq 0\) and that the medium is homogeneous so that the coefficients \(\lambda\) and \(\mu\) are constants. In some applications, it is more convenient to introduce the Young’s modulus \(E\) and the Poisson’s ratio \(v\) such that
The Poisson coefficient is such that \(1 \leq \nu<frac{1}{2}\); moreover, \(\nu \geq 0\) if \(\lambda \geq 0\). A Poisson’s ratio very close to \(\frac{1}{2}\) (i.e. Lamé coefficients such that the ratio \(\frac{1}{2}\) is very large) characterizes a practically incompressible material.
The model problem must be completed by boundary conditions. In the following, two cases are considered. The homogeneous Dirichlet problem:: we impose the boundary condition
 The pure traction problem (or Neumann problem)

we impose the boundary condition
where \(g: \partial \Omega \rightarrow \mathbb{R}^3\) represents a force field applied on the edge of \(\Omega\) and \(n\) the exterior normal to \(\Omega\).
Other boundary conditions are possible such as, for example, mixed DirichletNeumann conditions.
2. The mathematical framework
On \(\left[H^1(\Omega)\right]^3 \times\left[H^1(\Omega)\right]^3\), we introduce the bilinear form
For two matrices \(\sigma\) and \(\varepsilon\) in \(\mathbb{R}^{3,3}, \sigma: \varepsilon\) denotes their maximal contraction, which is equal to \(\sum_{i, j=1}^3 \sigma_{i j} \varepsilon_{i j}\). 
It is clear that \(a \in \mathcal{L}\left(\left[H^1(\Omega)\right]^3 \times\left[H^1(\Omega)\right]^3 ; \mathbb{R}\right)\). Moreover, the bilinear form \(a\) is symmetric and positive. On the other hand, the bilinear form \(a\) is singular. Indeed, by introducing the vector space of rigid displacements
we have the equivalence
We note that a rigid displacement consists in the compound of a translation and a rotation and that
since the only rigid displacement that preserves the boundary is the zero displacement.
2.1. Homogeneous Dirichlet problem.
In order to write the homogeneous Dirichlet problem in weak form, we introduce the functional space
and the linear form \(f_{\mathrm{D}} \in V_{\mathrm{D}}^{\prime}\) such that for any \(w \in V_{\mathrm{D}}\),
We obtain the following problem:
The wellposedness of the previous problem follows from the following inequality .Lemma First Korn inequality.
The inequality implies the coercivity of the bilinear form \(a\) on stem: [V_{\mathrm{D}}=left[H_0^{1(\Omega)\right\]}3] since
2.2. Pure traction problem
The mathematical study of the pure traction problem requires some precautions. Indeed, we cannot look for the solution \(u\) in \(\left[H^1(\Omega)\right\)^3] and ask that for any \(w \in\left[H^1(\Omega)\right\)^3], \(a(u, w)=\int_{\Omega} f \cdot w+\int_{\partial \Omega} g \cdot w\) because the bilinear form \(a\) is singular on \(\left[H^1(\Omega)\right\)^3].
A necessary condition for the existence of a solution is that
This equation expresses the fact that the resultant of all external forces and their moments are zero. Moreover, the solution \(u\), if it exists, is determined only to the nearest rigid displacement. We therefore consider the functional space
and the linear form \(f_{\mathrm{N}} \in V_{\mathrm{N}}^{\prime}\) such that for any \(w \in V_{\mathrm{N}}\),
We obtain the following problem:
The wellposedness of this problem follows from the following inequality
We show that the inequality implies the coercivity of the bilinear form \(a\) on \(V_{\mathrm{N}}\).
In continuum mechanics, the test function \(w\) involved in the weak formulations is interpreted as an admissible virtual displacement field and the weak formulations express the principle of virtual work. Moreover, the bilinear form \(a\) being symmetric and coercive on \(V_{\mathrm{D}}\) and \(V_{\mathrm{N}}\), the unique solution, respectively, minimizes on \(V_{\mathrm{D}}\) and \(V_{\mathrm{N}}\) the energy functional
\[\mathcal{E}_{\mathrm{D}}(v)=\frac{1}{2} \lambda \int_{\Omega}(\nabla \cdot v)^2+\frac{1}{2} \mu \int_{\Omega} \varepsilon(v): \varepsilon(v)f_{\mathrm{D}}(v)\]
and
\[\mathcal{E}_{\mathrm{N}}(v)=\frac{1}{2} \lambda \int_{\Omega}(\nabla \cdot v)^2+\frac{1}{2} \mu \int_{\Omega} \varepsilon(v): \varepsilon(v)f_{\mathbb{N}}(v)\]
We find the principle of least energy. The quadratic terms in \(v\) represent the elastic energy of deformation and the linear terms the potential energy under the external force field. 
3. Conformal approximation
We consider a conformal approximation of problems by Lagrangian finite elements. We suppose that \(\Omega\) is a polyhedron of \(\mathbb{R}^3\) and we consider a regular and conformal family of affine meshes of \(\Omega\) that we note \(\left{\mathcal{T}_h\right}_{h>0}\). We choose as reference finite element \(\widehat{K}, \widehat{P}, \widehat{Sigma}}\) a Lagrangian finite element of degree \(k \geq 1\).
 Homogeneous Dirichlet problem

In order to construct a \(V_{\mathrm{D}}\)conformal approximation space, we pose
The elements of \(V_{\mathrm{D} h}\) are the vector fields of which each component is in \(P_{\mathrm{c}, h}^k\) and which cancel on the boundary of \(\Omega\). We consider the following approximate problem:
which is clearly wellposed since \(a\) is coercive on \(V_{\mathrm{D}}\) and \(V_{\mathrm{D}} h} \subset V_{\mathrm{D}}\).
The estimation results from the lemma of Céa and the interpolation theorem which we apply component by component. The estimate results from the AubinNitsche lemma.
 Pure traction problem

For the pure traction problem, one way to eliminate the arbitrary rigid displacement at the discrete level is to:

impose that the displacement at one mesh node, \(a_0\), is zero;

choose three other mesh nodes, \(a_1, a_2, a_3\), and three unit vectors, \(\tau_1, \tau_2, \tau_3\), such that the set stem: [\left{\left(a_ia_0\right) \times \tau_i\right}_{1 \leq i \leq 3}] forms a basis of \(\mathbb{R}^3\);

impose that the displacement at the node \(a_i\) along the direction \(\boldsymbol{\tau}_i\) is zero.
This leads to the approximation space
We consider the following approximated problem:
Using the second Korn inequality, we show that the bilinear form \(a\) is coercive on \(V_{\mathrm{N} h}\) so that the discrete problem is well posed.
A sufficient condition for the problems to be regularizing is that the polyhedron \(\Omega\) is convex and that \(g=0\); see, for example, the book by Grisvard. 
4. Virtually incompressible materials: loss of coercivity
This situation occurs when the Poisson’s ratio is very close to \(\frac{1}{2}\), i.e. for practically incompressible materials. For such materials, it is observed that if we use a mesh that is not fine enough, the discrete solution is polluted by spurious oscillations. This phenomenon, called coercivity loss. The ratio \(\frac{lambda}{\mu}\) being very large, it is not reasonable to absorb it in the generic constants \(c\) appearing in the error estimates. We consider the bilinear form \(a\) defined above. We pose
where \(V\) is the functional space on which the continuous problem is posed. Under the hypothesis, we show that the ratio \(\frac{\omega_a}{\alpha_a}\) is of order \(\frac{\lambda}{\mu}\).
Using the convergence analysis presented in the previous section, we obtain the error estimate
with a constant \(c\) independent of \(h\) and the ratio \(\frac{\lambda}{\mu}\). This estimate shows that the mesh must be fine enough for the error to be controlled.