The equations of linear elasticity

Analysis of structures is one of the major activities of modern engineering, which likely makes the PDE modeling the deformation of elastic bodies the most popular PDE in the world. It takes just one page of code to solve the equations of 2D or 3D elasticity in FEniCS, and the details follow below.

(click and drag to rotate the 3D scenes)

PDE problem

The equations governing small elastic deformations of a body $\Omega$ can be written as

\begin{aligned} (3.20) & - \nabla \cdot σ & = f in Ω, \\ (3.21) & σ & = λ tr (ε) I + 2 μ ε, \\ (3.22) & ε & = \frac{1}{2} (\nabla u + (\nabla u)^{⊤}), \end{aligned}

$\begin{align} -\nabla\cdot\sigma &= f\hbox{ in }\Omega, \tag{3.20}\\ \sigma &= \lambda\,\hbox{tr}\,(\varepsilon) I + 2\mu\varepsilon, \tag{3.21}\\ \varepsilon &= \frac{1}{2}\left(\nabla u + (\nabla u)^{\top}\right), \tag{3.22} \end{align}$ where

σ

$\sigma$ is the stress tensor,

f

$f$ is the body force per unit volume,

λ

$\lambda$ and

μ

$\mu$ are

Lam\'e's

$\text{Lam\'e's}$ elasticity parameters for the material in

Ω

$\Omega$ ,

I

$I$ is the identity tensor,

t r

$\mathrm{tr}$ is the trace operator on a tensor,

ε

$\varepsilon$ is the symmetric strain-rate tensor (symmetric gradient), and

u

$u$ is the displacement vector field. We have here assumed isotropic elastic conditions.

We combine (3.21) and (3.22) to obtain

\begin{matrix} (3.23) & σ = λ (\nabla \cdot u) I + μ (\nabla u + (\nabla u)^{⊤}) . \end{matrix}

$\begin{equation} \sigma = \lambda(\nabla\cdot u)I + \mu(\nabla u + (\nabla u)^{\top})\tp \tag{3.23} \end{equation}$ Note that (3.20)--(3.22) can easily be transformed to a single vector PDE for

u

$u$ , which is the governing PDE for the unknown

u

$u$ (Navier's equation). In the derivation of the variational formulation, however, it is convenient to keep the equations split as above.

Variational formulation

The variational formulation of (3.20)--(3.22) consists of forming the inner product of (3.20) and a vector test function $v\in \hat{V}$ , where $\hat{V}$ is a vector-valued test function space, and integrating over the domain $\Omega$ :

- \int_{Ω} (\nabla \cdot σ) \cdot v d x = \int_{Ω} f \cdot v d x .

$-\int_\Omega (\nabla\cdot\sigma) \cdot v \dx = \int_\Omega f\cdot v\dx\tp$ Since

\nabla \cdot σ

$\nabla\cdot\sigma$ contains second-order derivatives of the primary unknown

u

$u$ , we integrate this term by parts:

- \int_{Ω} (\nabla \cdot σ) \cdot v d x = \int_{Ω} σ : \nabla v d x - \int_{\partial Ω} (σ \cdot n) \cdot v d s,

$-\int_\Omega (\nabla\cdot\sigma) \cdot v \dx = \int_\Omega \sigma : \nabla v\dx - \int_{\partial\Omega} (\sigma\cdot n)\cdot v \ds,$ where the colon operator is the inner product between tensors (summed pairwise product of all elements), and

n

$n$ is the outward unit normal at the boundary. The quantity

σ \cdot n

$\sigma\cdot n$ is known as the traction or stress vector at the boundary, and is often prescribed as a boundary condition. We here assume that it is prescribed on a part

\partial Ω_{T}

$\partial\Omega_T$ of the boundary as

σ \cdot n = T

$\sigma\cdot n = T$ . On the remaining part of the boundary, we assume that the value of the displacement is given as a Dirichlet condition. We thus obtain

\int_{Ω} σ : \nabla v d x = \int_{Ω} f \cdot v d x + \int_{\partial Ω_{T}} T \cdot v d s .

$\int_\Omega \sigma : \nabla v \dx = \int_\Omega f\cdot v \dx + \int_{\partial\Omega_T} T\cdot v\ds\tp$ Inserting the expression (3.23) for

σ

$\sigma$ gives the variational form with

u

$u$ as unknown. Note that the boundary integral on the remaining part

\partial Ω ∖ \partial Ω_{T}

$\partial\Omega\setminus\partial\Omega_T$ vanishes due to the Dirichlet condition.

We can now summarize the variational formulation as: find $u\in V$ such that

\begin{matrix} (3.24) & a (u, v) = L (v) \forall v \in \hat{V}, \end{matrix}

$\begin{equation} a(u,v) = L(v)\quad\forall v\in\hat{V}, \tag{3.24} \end{equation}$ where

\begin{aligned} (3.25) & a (u, v) & = \int_{Ω} σ (u) : \nabla v d x, \\ (3.26) & σ (u) & = λ (\nabla \cdot u) I + μ (\nabla u + (\nabla u)^{⊤}), \\ (3.27) & L (v) & = \int_{Ω} f \cdot v d x + \int_{\partial Ω_{T}} T \cdot v d s . \end{aligned}

$\begin{align} a(u,v) &= \int_\Omega\sigma(u) :\nabla v \dx, \tag{3.25} \\ \sigma(u) &= \lambda(\nabla\cdot u)I + \mu(\nabla u + (\nabla u)^{\top}), \tag{3.26}\\ L(v) &= \int_\Omega f\cdot v\dx + \int_{\partial\Omega_T} T\cdot v\ds\tp \tag{3.27} \end{align}$

One can show that the inner product of a symmetric tensor $A$ and an anti-symmetric tensor $B$ vanishes. If we express $\nabla v$ as a sum of its symmetric and anti-symmetric parts, only the symmetric part will survive in the product $\sigma :\nabla v$ since $\sigma$ is a symmetric tensor. Thus replacing $\nabla u$ by the symmetric gradient $\epsilon(u)$ gives rise to the slightly different variational form

\begin{matrix} (3.28) & a (u, v) = \int_{Ω} σ (u) : ε (v) d x, \end{matrix}

$\begin{equation} a(u,v) = \int_\Omega\sigma(u) :\varepsilon(v) \dx, \tag{3.28} \end{equation}$ where

ε (v)

$\varepsilon(v)$ is the symmetric part of

\nabla v

$\nabla v$ :

ε (v) = \frac{1}{2} (\nabla v + (\nabla v)^{⊤}) .

$\varepsilon(v) = \frac{1}{2}\left(\nabla v + (\nabla v)^{\top}\right)\tp$ The formulation (3.28) is what naturally arises from minimization of elastic potential energy and is a more popular formulation than (3.25).

FEniCS implementation

Test problem

As a test example, we will model a clamped beam deformed under its own weight in 3D. This can be modeled by setting the right-hand side body force per unit volume to $f=(0,0,-\varrho g)$ with $\varrho$ the density of the beam and $g$ the acceleration of gravity. The beam is box-shaped with length $L$ and has a square cross section of width $W$ . We set $u=\ub = (0,0,0)$ at the clamped end, $x=0$ . The rest of the boundary is traction free; that is, we set $T = 0$ .

FEniCS implementation

We first list the code and then comment upon the new constructions compared to the previous examples we have seen.

from fenics import *

# Scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

This example program can be found in the file elasticbeam.py.

Vector function spaces

The primary unknown is now a vector field $u$ and not a scalar field, so we need to work with a vector function space:

V = VectorFunctionSpace(mesh, 'P', 1)

With u = Function(V) we get u as a vector-valued finite element function with three components for this 3D problem.

Constant vectors

For the boundary condition $u=(0, 0, 0)$ , we must set a vector value to zero, not just a scalar. Such a vector constant is specified as Constant((0, 0, 0)) in FEniCS. The corresponding 2D code would use Constant((0, 0)). Later in the code, we also need f as a vector and specify it as Constant((0, 0, rho*g)).

`nabla_grad`

The gradient and divergence operators now have a prefix nabla_. This is strictly not necessary in the present problem, but recommended in general for vector PDEs arising from continuum mechanics, if you interpret $\nabla$ as a vector in the PDE notation; see the box about nabla_grad in the section Variational formulation.

Stress computation

As soon as the displacement u is computed, we can compute various stress measures. We will compute the von Mises stress defined as $\sigma_M = \sqrt{\frac{3}{2}s:s}$ where $s$ is the deviatoric stress tensor

s = σ - \frac{1}{3} t r (σ) I .

$s = \sigma - \frac{1}{3}\mathrm{tr}\,(\sigma)\,I\tp$ There is a one-to-one mapping between these formulas and the FEniCS code:

s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
von_Mises = sqrt(3./2*inner(s, s))

The von_Mises variable is now an expression that must be projected to a finite element space before we can visualize it:

V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)