Partial differential equations » Solving the elasticity equation

Showcases the resolution of the elasticity system.

Image

In this example we will seek to solve the linearized elasticity system. Most of the reference material in the introduction is taken from [1].

Introduction

In the elasticity setting, we are most typically interested in characterizing the elastic displacement $ u_\Omega : \Omega \rightarrow \mathbb{R}^d $ of a mechanical structure $ \Omega \subset \mathbb{R}^d $ , which in the mathematical setting is most usually defined as a bounded, Lipschitz domain, whose boundary is decomposed into three disjoint pieces: $ \Gamma_D \cup \Gamma_N \cup \Gamma_0$ . Each piece of the boundary has the following semantic meanings:

  • on the segment $ \Gamma_D \neq \emptyset $ is where the structure is "fixed" or "attached", which is to say that the structure cannot be displaced or moved
  • on the segment $ \Gamma_N $ is where the traction loads $ g : \Gamma_N \rightarrow \mathbb{R}^d $ are applied
  • the remaining segment $ \Gamma_0 $ corresponds to a traction-free boundary; this is the region where the shape is free to move about

Additionally, $ \Omega $ can be subjected to body forces $ f : \Omega \rightarrow \mathbb{R}^d $ . In general, this framework gives rise to the linearized elasticity system:

\[ \left\{ \begin{aligned} - \nabla \cdot \sigma (u_\Omega) &= 0 && \mathrm{in} \ \Omega \\ u_\Omega &= 0 && \mathrm{on} \ \Gamma_D \\ \sigma(u_\Omega) \cdot n &= 0 && \mathrm{on} \ \Gamma_0 \\ \sigma(u_\Omega) \cdot n &= f && \mathrm{on} \ \Gamma_N \end{aligned} \right. , \]

where $ \sigma(u_\Omega) $ is known as the stress tensor of the displacement field $ u_\Omega $ . Furthermore, this stress tensor is related to the strain tensor

\[ e(u_\Omega) := \dfrac{1}{2}\left( \nabla u + \nabla u^T \right) \]

via Hooke's law

\[ \sigma(u_\Omega) = Ae(u_\Omega) \ , \]

where for any $ e $ in the set $ \mathbb{R}^{d \times d}_s $ of symmetric $ d \times d $ matrices,

\[ Ae = 2 \mu e + \lambda \mathrm{tr}(e) I \ . \]

Here $ I $ is the identity $ d \times d $ matrix, and $ \lambda $ , $ \mu $ are the Lamé parameters of the constituent material, satisfying $ \mu > 0 $ and $ \lambda + 2\mu / d > 0 $ .

Weak Formulation

As with the Poisson example, we work with the variational formulation. We seek to solve:

\[ \text{Find } u_\Omega \in H^1_{\Gamma_D}(\Omega)^d \text{ s.t. } \forall v \in H^1_{\Gamma_D}(\Omega)^d, \quad \int_{\Omega} Ae(u_\Omega) : e(v) \ dx = \int_{\Gamma_N} g \cdot v \ dx \ , \]

where $ H^1_{\Gamma_D}(\Omega)^d := \{ u \in H^1(\Omega)^d, \ u = 0 \text{ on } \Gamma_D \} $ . This problem is well-posed by the Lax-Milgram theorem (Korn's inequality provides coercivity).

For implementation, the bilinear form can be expanded as:

\[ \int_{\Omega} \lambda (\mathrm{div} \ u) \ (\mathrm{div} \ v) \ dx + \int_{\Omega} \mu \, e(u) : e(v) \ dx = \int_{\Gamma_N} g \cdot v \ dx \ . \]

Implementation

Setup

We load the mesh using Mesh::load(), define a vector-valued P1 finite element space with $ d $ components (one per spatial dimension), and set the material parameters:

#include <Rodin/Solver.h>
#include <Rodin/Geometry.h>
#include <Rodin/Variational.h>
#include <Rodin/Solid.h>

using namespace Rodin;
using namespace Rodin::Geometry;
using namespace Rodin::Variational;

// Define boundary attributes
Attribute Gamma = 1, GammaD = 2, GammaN = 3, Gamma0 = 4;

// Load mesh
Mesh mesh;
mesh.load(meshFile, IO::FileFormat::MFEM);
mesh.getConnectivity().compute(1, 2);

// Vector-valued P1 space (d components)
size_t d = mesh.getSpaceDimension();
P1 fes(mesh, d);

// Lamé coefficients
const Real lambda = 0.5769, mu = 0.3846;

// Applied traction force
VectorFunction f{0, -1};

Problem Formulation

The variational problem is expressed through the Problem class. The LinearElasticityIntegral provides a convenient shorthand for the full elasticity bilinear form. BoundaryIntegral applies the traction load on boundary segment $ \Gamma_N $ , and DirichletBC fixes the displacement on $ \Gamma_D $ :

TrialFunction u(fes);
TestFunction  v(fes);

Problem elasticity(u, v);
elasticity = LinearElasticityIntegral(u, v)(lambda, mu)
           - BoundaryIntegral(f, v).over(GammaN)
           + DirichletBC(u, VectorFunction{0, 0}).on(GammaD);

Solving and Output

We solve with the CG solver and save the results using IO::FileFormat::MFEM:

Solver::CG(elasticity).solve();

// Save in MFEM format
u.getSolution().save("Elasticity.gf", IO::FileFormat::MFEM);
mesh.save("Elasticity.mesh", IO::FileFormat::MFEM);

For visualization in ParaView, export using IO::XDMF:

#include <Rodin/IO/XDMF.h>

IO::XDMF xdmf("Elasticity");
xdmf.grid().setMesh(mesh).add("displacement", u.getSolution());
xdmf.write();

The full source code may be found below:

« Solving the Poisson equation | Partial differential equations | Solving the Helmholtz equation »