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

\[\nabla \cdot \sigma+f=0 \quad \text { in } \Omega\]
Strain tensor.

Let us note \(\varepsilon(u): \Omega \rightarrow \mathbb{R}^{3,3}\) the linearized deformation field defined by

\[\varepsilon(u)=\frac{1}{2}\left(\nabla u+\nabla u^T\right) .\]
\(\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\}\).
Stress-strain relation.

In the framework of linear elasticity, the stress field is expressed as a function of the linearized strain field as

\[\sigma(u)=\lambda \operatorname{tr}(\varepsilon(u)) \mathcal{I}_3+2 \mu \varepsilon(u)\]

where \(\lambda\) and \(\mu\) are the Lamé coefficients and \(\mathcal{I}_3\) the identity matrix in \(\mathbb{R}^{3,3}\). Using the definition of the deformation tensor, we get

\[\sigma(u)=\lambda(\nabla \cdot u) \mathcal{I}_3+\mu\left(\nabla u+\nabla u^T\right)\]

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 \(\nu\) such that

Relation between Lamé coefficients and Young’s modulus.
\[\left\{\begin{aligned} E=\mu \frac{3 \lambda+2 \mu}{\lambda+\mu} \quad \text { and } \quad v = \frac{\lambda} {2 (\lambda + \mu)}\\ \lambda = \frac{E \ \nu} {(1 + \nu) \ (1 - 2\nu)},\quad \mu = \frac{E} {2 (1 + \nu)} \end{aligned} \right.\]

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.

Homogeneous Dirichlet problem

we impose the boundary condition

\[u=0 \quad \text { on } \partial \Omega\]
Pure traction problem (or Neumann problem)

we impose the boundary condition

\[\sigma(u) \cdot n=g \quad \text { on } \partial \Omega\]

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 Dirichlet-Neumann conditions.

2. The mathematical framework

Weak formulation.

On \(\left[H^1(\Omega)\right]^3 \times\left[H^1(\Omega)\right]^3\), we introduce the bilinear form

\[a(v, w)=\int_{\Omega} \sigma(v): \varepsilon(w)=\int_{\Omega} \lambda(\nabla \cdot v)(\nabla \cdot w)+\int_{\Omega} 2 \mu \varepsilon(v): \varepsilon(w).\]
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

\[\mathcal{R}=\left\{z \in\left[H^1(\Omega)\right]^3 ; \exists(\alpha, \beta) \in \mathbb{R}^3 \times \mathbb{R}^3, \forall x \in \Omega, z(x)=\alpha+\beta \times x\right\}\]

we have the equivalence

\[(z \in \mathcal{R}) \Longleftrightarrow\left(\forall v \in \left[H^1(\Omega)\right]^3, a(z, v)=0\right) .\]

We note that a rigid displacement consists in the compound of a translation and a rotation and that

\[\mathcal{R} \cap\left[H_0^1(\Omega)\right]^3=\{0\}\]

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

\[V_{\mathrm{D}}=\left[H_0^1(\Omega)\right]^3\]

and the linear form \(f_{\mathrm{D}} \in V_{\mathrm{D}}^{\prime}\) such that for any \(w \in V_{\mathrm{D}}\),

\[f_{\mathrm{D}}(w)=\int_{\Omega} f \cdot w\]
Problem Formulation

We obtain the following problem:

\[\text { Find } u \in V_{\mathrm{D}} \text { such that } a(u, w)=f_{\mathrm{D}}(w), \quad \forall w \in V_{\mathrm{D}} .\]

The well-posedness of the previous problem follows from the following inequality

Lemma First Korn inequality.

Let \(\Omega\) be a domain of \(\mathbb{R}^3\). Let \(\|\varepsilon(v)\|_{0, \Omega}=\left(\int_{\Omega} \varepsilon(v): \varepsilon(v)\right)^{\frac{1}{2}}\). There exists a constant \(\boldsymbol{\kappa}_{\Omega}\) such that [[1st-korn-ineq]]]

\[\forall v \in\left[H_0^1(\Omega)\right]^3, \quad \kappa_{\Omega}|v\|_{1, \Omega}.\]

The inequality implies the coercivity of the bilinear form \(a\) on \(V_{\mathrm{D}}=\left[H_0^1(\Omega)\right]^3\) since

\[\forall v \in\left[H_0^1(\Omega)\right]^3, \quad a(v, v) \geq 2 \mu|\varepsilon(v)\|_{0, \Omega}^2 \geq 2 \mu \kappa_{\Omega}^2|v||_{1, \Omega}^2 .\]

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\).

Necessary condition for the existence of a solution.

A necessary condition for the existence of a solution is that

\[\forall z \in \mathcal{R}, \quad \int_{\Omega} f \cdot z+\int_{\partial \Omega} g \cdot z=0 .\]

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

\[V_{\mathrm{N}}=\left\{v \in\left[H^1(\Omega)\right]^3; \int_{\Omega} v=0; \int_{\Omega} \nabla \times v=0\right\},\]

and the linear form \(f_{\mathrm{N}} \in V_{\mathrm{N}}^{\prime}\) such that for any \(w \in V_{\mathrm{N}}\),

\[f_{\mathrm{N}}(w)=\int_{\Omega} f \cdot w+\int_{\partial \Omega} g \cdot w .\]

We obtain the following problem:

\[\left\{\begin{array}{l} \text { Find } u \in V_{\mathrm{N}} \text { such that } \\ a(u, w)=f_{\mathrm{N}}(w), \quad \forall w \in V_{\mathrm{N}} \end{array}\right.\]

The well-posedness of this problem follows from the following inequality

Lemma Second inequality of Korn.

Let \(\Omega\) be a domain of \(\mathbb{R}^3\). There exists a constant \(\kappa_{\Omega}^{\prime}\) such that [[2nd-korn-ineq]]

\[\forall v \in\left[H^1(\Omega)\right]^3, \quad \boldsymbol{\kappa}_{\Omega}^{\prime}\|v\|_{1, \Omega} \leq \|\varepsilon(v)\|_{0, \Omega}+\|v\|_{0, \Omega}\]

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 \(\left\{\widehat{K}, \widehat{P}, \widehat{Sigma}\right\}\) a Lagrangian finite element of degree \(k \geq 1\).

3.1. Homogeneous Dirichlet problem

In order to construct a \(V_{\mathrm{D}}\)-conformal approximation space, we pose

\[V_{\mathrm{D} h}=\left[P_{\mathrm{c}, h}^k\right]^3 \cap\left[H_0^1(\Omega)\right]^3,\]

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:

\[\left\{\begin{array}{l} \text { Search } u_h \in V_{\mathrm{D} h} \text { such as } \\ a\left(u_h, w_h\right)=f_{\mathrm{D}}\left(w_h\right), \quad \forall w_h \in V_{\mathrm{D} h}, \end{array}\right.\]

which is clearly well-posed since \(a\) is coercive on \(V_{\mathrm{D}}\) and \(V_{\mathrm{D}h} \subset V_{\mathrm{D}}\).

Theorem Convergence.

With the above assumptions, we suppose that the unique solution \(u\) is in \(\left[H^{k+1}(\Omega) \cap H_0^1(\Omega)\right]^3\). Then, there exists a constant \(c\) such that for all \(h\),

\[\left\|u-u_h\right\|_{1, \Omega} \leq c b^k|u|_{k+1, \Omega}\]

Moreover, if the problem is regularizing (i.e., if there exists a constant \(c_s\) such that for any \(f \in\left[L^2(\Omega)\right]^3\), the unique solution satisfies \(\|u|_{2, \Omega} \leq c_s\|f|_{0, \Omega}\)), there exists a constant \(c\) such that for any \(h\),

\[\left\|u-u_h\right\|_{0, \Omega} \leq c h^{k+1}|u|_{k+1, \Omega} .\]
Elements of proof

The estimation results from the lemma of Céa and the interpolation theorem which we apply component by component. The estimate results from the Aubin-Nitsche lemma.

3.2. Pure traction problem

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

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

  2. 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 \(\left\{\left(a_i-a_0\right) \times \tau_i\right\}_{1 \leq i \leq 3}\) forms a basis of \(\mathbb{R}^3\);

  3. impose that the displacement at the node \(a_i\) along the direction \(\boldsymbol{\tau}_i\) is zero.

Pure Traction approximation space

This leads to the approximation space

\[\begin{aligned} V_{\mathrm{N} h}=\left\{v_h \in\left[\mathcal{C}^0(\bar{\Omega})\right]^3 ;\right. & \forall K \in \mathcal{T}_h, v_h \circ T_K \in[\widehat{P}]^3 ; \\ & \left.v_h\left(a_0\right)=0 ; v_h\left(a_i\right) \cdot \tau_i=0, i \in\{1,2,3\}\right\} . \end{aligned}\]

We consider the following approximated problem:

\[\left\{\begin{array}{l} \text { Find } u_h \in V_{\mathrm{N} h} \text { such as } \\ a\left(u_h, w_h\right)=f_{\mathrm{N}}\left(w_h\right), \quad \forall w_h \in V_{\mathrm{N} h} . \end{array}\right.\]

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.

Theorem Convergence.

With the above assumptions, we suppose that the unique solution u of continuous problem is in \(\left[H^{k+1}(\Omega)\right]^3 \cap V_{\mathrm{N}}\). Then there exists a constant \(c\) such that for all \(h\),

\[\left\|u-u_h\right\|_{1, \Omega} \leq c h^k|u|_{k+1, \Omega} .\]

Moreover, if \(g=0\) and if the problem is regularizing, (i.e. if there exists a constant \(c_s\) such that for any \(f \in\left[L^2(\Omega)\right]^3\), the unique solution of problem with \(g=0\) satisfies \(\left.\|u\|_{2, \Omega} \leq c_s\|f\|_{0, \Omega}\right)\)), there exists a constant \(c\) such that for any \(h\),

\[\left\|u-u_h\right\|_{0, \Omega} \leq c h^{k+1}|u|_{k+1, \Omega}\]

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.

3.3. Example

# To be done

4. Virtually incompressible materials: loss of coercivity

We are now interested in materials whose ratio of Lamé coefficients is such that ]

\[\frac{\lambda}{\mu} \gg 1\]

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.

Definition Coercivity and Continuity constants.

We pose

\[\begin{aligned} & \alpha_a=\inf _{v \in V} \frac{a(v, v)}{\|v\|_{1, \Omega}^2} \\ & \omega_a=\sup _{v \in V} \sup _{w \in V} \frac{a(v, w)}{\|v\|_{1, \Omega}\|w\|_{1, \Omega}} \end{aligned}\]

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

\[\left\|u-u_h\right\|_{1, \Omega} \leq c \frac{\lambda}{\mu} h^k|u|_{k+1, \Omega}\]

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.

4.1. Example of coercivity loss

# TODO

Bibliography

  • [grisvard] GRISVARD (P.), Singularities in Boundary Value Problems. Masson, Paris, France, 1992.