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 strucure 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 $ 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 $ . In general, as we have said in the Poisson example, we are mainly interested in the variational formulation of this equation. That is we seek to solve the problem:

\[ \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_{\Omega} f \cdot v \ dx + \int_{\Gamma_N} g \cdot v \ dx \ , \]

As far as the mathematical setting is concerned $ f \in L^2(\Omega)^d $ , $ g \in L^2(\Gamma_N)^d $ and the space $ H^1_{\Gamma_D}(\Omega)^d $ is defined by:

\[ H^1_{\Gamma_D}(\Omega)^d := \{ u \in H^1(\Omega)^d, \ u = 0 \text{ on } \Gamma_D = 0 \} \ . \]

In particular, this problem is well posed via the classical Lax-Milgram theory. In this example we will be implementing the variational formulation using Rodin's Variational module.

Implementation

For simplicity, we will disregard any body forces by taking $ f = 0 $ . In this vein, we also note that the integral in the variational formulation may be expanded and rewritten into

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

which should allow us to more easily translate it into code. In fact, the implementation follows exactly the same philosophy as the Poisson case:

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

// Load mesh
Mesh Omega = Mesh::load(meshFile);

// Functions
int d = 2;
FiniteElementSpace<H1> Vh(Omega, d);

// Lamé coefficients
auto mu     = ScalarFunction(0.3846),
     lambda = ScalarFunction(0.5769);

// Pull force
auto f = VectorFunction{0, -1};

We define the relevant boundary attributes (which we should know beforehand), we load the mesh, we build the (vectorial) finite element space, and define our parameters. The only real difference will be the variational formulation! But even this is easy to implement if you know what the weak formulation looks like:

// Define problem
TrialFunction u(Vh);
TestFunction  v(Vh);
Problem elasticity(u, v);
elasticity = Integral(lambda * Div(u), Div(v))
           + Integral(
               mu * (Jacobian(u) + Jacobian(u).T()), 0.5 * (Jacobian(v) + Jacobian(v).T()))
           - BoundaryIntegral(f, v).over(GammaN)
           + DirichletBC(u, VectorFunction{0, 0}).on(GammaD);

In this case we have just written out the definition of the terms $ e(u_\Omega) $ and $ e(v) $ . Rodin is smart enough to figure out how to build the underlying linear system to solve and assemble it when we need to solve it.

This brings us to the last point which is solving the variational problem. In exactly the same fashion as the previous example, we need only choose a Solver which should allows us to solve the problem.

// Solve problem
Solver::CG().setMaxIterations(200)
            .setRelativeTolerance(1e-12)
            .printIterations(true)
            .solve(elasticity);

// Save solution
u.getGridFunction().save("u.gf");
Omega.save("Omega.mesh");

The full source code may be found below:

« Solving the Poisson equation | Partial differential equations | Solving on surfaces with flat elements »