Linear Elasticity Equations

Summary of the Equations of Linear Elasticity

The stiffness matrix of linear elastic isotropic material contains two parameters:

E, Young's modulus (elastic modulus)
ν, Poisson’s ratio

Define the following quantities.

$\begin{array}{l} σ = stress \\ f = body force \\ ε = strain \\ u = displacement \end{array}$

The equilibrium equation is

$- \nabla \cdot σ = f$

The linearized, small-displacement strain-displacement relationship is

$ε = \frac{1}{2} (\nabla u + \nabla u^{T})$

The balance of angular momentum states that stress is symmetric:

$σ_{i j} = σ_{j i}$

The Voigt notation for the constitutive equation of the linear isotropic model is

$[\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{33} \\ σ_{23} \\ σ_{13} \\ σ_{12} \end{matrix}] = \frac{E}{(1 + ν) (1 - 2 ν)} [\begin{matrix} 1 - ν & ν & ν & 0 & 0 & 0 \\ ν & 1 - ν & ν & 0 & 0 & 0 \\ ν & ν & 1 - ν & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 - 2 ν & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 - 2 ν & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 - 2 ν \end{matrix}] [\begin{matrix} ε_{11} \\ ε_{22} \\ ε_{33} \\ ε_{23} \\ ε_{13} \\ ε_{12} \end{matrix}]$

The expanded form uses all the entries in σ and ε takes symmetry into account.

[\begin{matrix} σ_{11} \\ σ_{12} \\ σ_{13} \\ σ_{21} \\ σ_{22} \\ σ_{23} \\ σ_{31} \\ σ_{32} \\ σ_{33} \end{matrix}] = \frac{E}{(1 + ν) (1 - 2 ν)} [\begin{matrix} 1 - ν & 0 & 0 & 0 & ν & 0 & 0 & 0 & ν \\ • & 1 - 2 ν & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ • & • & 1 - 2 ν & 0 & 0 & 0 & 0 & 0 & 0 \\ • & • & • & 1 - 2 ν & 0 & 0 & 0 & 0 & 0 \\ • & • & • & • & 1 - ν & 0 & 0 & 0 & ν \\ • & • & • & • & • & 1 - 2 ν & 0 & 0 & 0 \\ • & • & • & • & • & • & 1 - 2 ν & 0 & 0 \\ • & • & • & • & • & • & • & 1 - 2 ν & 0 \\ • & • & • & • & • & • & • & • & 1 - ν \end{matrix}] [\begin{matrix} ε_{11} \\ ε_{12} \\ ε_{13} \\ ε_{21} \\ ε_{22} \\ ε_{23} \\ ε_{31} \\ ε_{32} \\ ε_{33} \end{matrix}]

(1)

In the preceding diagram, • means the entry is symmetric.

3D Linear Elasticity Problem

The toolbox form for the equation is

$- \nabla \cdot (c \otimes \nabla u) = f$

But the equations in the summary do not have ∇u alone, it appears together with its transpose:

$ε = \frac{1}{2} (\nabla u + \nabla u^{T})$

It is a straightforward exercise to convert this equation for strain ε to ∇u. In column vector form,

$\nabla u = [\begin{matrix} \partial u_{x} / \partial x \\ \partial u_{x} / \partial y \\ \partial u_{x} / \partial z \\ \partial u_{y} / \partial x \\ \partial u_{y} / \partial y \\ \partial u_{y} / \partial z \\ \partial u_{z} / \partial x \\ \partial u_{z} / \partial y \\ \partial u_{z} / \partial z \end{matrix}]$

Therefore, you can write the strain-displacement equation as

$ε = [\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 \\ 0 & \frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}] \nabla u \equiv A \nabla u$

where A stands for the displayed matrix. So rewriting Equation 1, and recalling that • means an entry is symmetric, you can write the stiffness tensor as

$\begin{matrix} σ = \frac{E}{(1 + ν) (1 - 2 ν)} [\begin{matrix} 1 - ν & 0 & 0 & 0 & ν & 0 & 0 & 0 & ν \\ • & 1 - 2 ν & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ • & • & 1 - 2 ν & 0 & 0 & 0 & 0 & 0 & 0 \\ • & • & • & 1 - 2 ν & 0 & 0 & 0 & 0 & 0 \\ • & • & • & • & 1 - ν & 0 & 0 & 0 & ν \\ • & • & • & • & • & 1 - 2 ν & 0 & 0 & 0 \\ • & • & • & • & • & • & 1 - 2 ν & 0 & 0 \\ • & • & • & • & • & • & • & 1 - 2 ν & 0 \\ • & • & • & • & • & • & • & • & 1 - ν \end{matrix}] A \nabla u \\ = \frac{E}{(1 + ν) (1 - 2 ν)} [\begin{matrix} 1 - ν & 0 & 0 & 0 & ν & 0 & 0 & 0 & ν \\ 0 & 1 / 2 - ν & 0 & 1 / 2 - ν & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 / 2 - ν & 0 & 0 & 0 & 1 / 2 - ν & 0 & 0 \\ 0 & 1 / 2 - ν & 0 & 1 / 2 - ν & 0 & 0 & 0 & 0 & 0 \\ ν & 0 & 0 & 0 & 1 - ν & 0 & 0 & 0 & ν \\ 0 & 0 & 0 & 0 & 0 & 1 / 2 - ν & 0 & 1 / 2 - ν & 0 \\ 0 & 0 & 1 / 2 - ν & 0 & 0 & 0 & 1 / 2 - ν & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 / 2 - ν & 0 & 1 / 2 - ν & 0 \\ ν & 0 & 0 & 0 & ν & 0 & 0 & 0 & 1 - ν \end{matrix}] \nabla u \end{matrix}$

Make the definitions

$\begin{array}{l} μ = \frac{E}{2 (1 + ν)} \\ λ = \frac{E ν}{(1 + ν) (1 - 2 ν)} \\ \frac{E (1 - ν)}{(1 + ν) (1 - 2 ν)} = 2 μ + λ \end{array}$

and the equation becomes

$σ = [\begin{matrix} 2 μ + λ & 0 & 0 & 0 & λ & 0 & 0 & 0 & λ \\ 0 & μ & 0 & μ & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & μ & 0 & 0 & 0 & μ & 0 & 0 \\ 0 & μ & 0 & μ & 0 & 0 & 0 & 0 & 0 \\ λ & 0 & 0 & 0 & 2 μ + λ & 0 & 0 & 0 & λ \\ 0 & 0 & 0 & 0 & 0 & μ & 0 & μ & 0 \\ 0 & 0 & μ & 0 & 0 & 0 & μ & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & μ & 0 & μ & 0 \\ λ & 0 & 0 & 0 & λ & 0 & 0 & 0 & 2 μ + λ \end{matrix}] \nabla u \equiv c \nabla u$

If you are solving a 3-D linear elasticity problem by using PDEModel instead of femodel, use the elasticityC3D(E,nu) function (included in your software) to obtain the c coefficient. This function uses the linearized, small-displacement assumption for an isotropic material. For examples that use this function, see StationaryResults.

By default, the toolbox uses the zero Neumann boundary condition and assumes that the boundary is stress-free.

Plane Stress

Plane stress is a condition that prevails in a flat plate in the x-y plane, loaded only in its own plane and without z-direction restraint. For plane stress, σ₁₃ = σ₂₃ = σ₃₁ = σ₃₂ = σ₃₃ = 0. Assuming isotropic conditions, the Hooke's law for plane stress gives the following strain-stress relation:

$[\begin{matrix} ε_{11} \\ ε_{22} \\ 2 ε_{12} \end{matrix}] = \frac{1}{E} [\begin{matrix} 1 & - ν & 0 \\ - ν & 1 & 0 \\ 0 & 0 & 2 + 2 ν \end{matrix}] [\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{12} \end{matrix}]$

Inverting this equation, obtain the stress-strain relation:

$(\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{12} \end{matrix}) = \frac{E}{1 - ν^{2}} (\begin{matrix} 1 & ν & 0 \\ ν & 1 & 0 \\ 0 & 0 & \frac{1 - ν}{2} \end{matrix}) (\begin{matrix} ε_{11} \\ ε_{22} \\ 2 ε_{12} \end{matrix})$

Convert the equation for strain ε to ∇u.

$ε = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] \nabla u \equiv A \nabla u$

Now you can rewrite the stiffness matrix as

$[\begin{matrix} σ_{11} \\ σ_{12} \\ σ_{21} \\ σ_{22} \end{matrix}] = [\begin{matrix} \frac{E}{1 - ν^{2}} & 0 & 0 & \frac{E ν}{1 - ν^{2}} \\ 0 & \frac{E}{2 (1 + ν)} & \frac{E}{2 (1 + ν)} & 0 \\ 0 & \frac{E}{2 (1 + ν)} & \frac{E}{2 (1 + ν)} & 0 \\ \frac{E ν}{1 - ν^{2}} & 0 & 0 & \frac{E}{1 - ν^{2}} \end{matrix}] \nabla u = [\begin{matrix} \frac{2 μ (μ + λ)}{2 μ + λ} & 0 & 0 & \frac{2 λ μ}{2 μ + λ} \\ 0 & μ & μ & 0 \\ 0 & μ & μ & 0 \\ \frac{2 λ μ}{2 μ + λ} & 0 & 0 & \frac{2 μ (μ + λ)}{2 μ + λ} \end{matrix}] \nabla u$

Plane Strain

Plane strain is a deformation state where there are no displacements in the z-direction, and the displacements in the x- and y-directions are functions of x and y but not z. The stress-strain relation is only slightly different from the plane stress case, and the same set of material parameters is used.

For plane strain, ε₁₃ = ε₂₃ = ε₃₁ = ε₃₂ = ε₃₃ = 0. Assuming isotropic conditions, the stress-strain relation can be written as follows:

$(\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{12} \end{matrix}) = \frac{E}{(1 + ν) (1 - 2 ν)} (\begin{matrix} 1 - ν & ν & 0 \\ ν & 1 - ν & 0 \\ 0 & 0 & \frac{1 - 2 ν}{2} \end{matrix}) (\begin{matrix} ε_{11} \\ ε_{22} \\ 2 ε_{12} \end{matrix})$

Convert the equation for strain ε to ∇u.

$ε = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] \nabla u \equiv A \nabla u$

Now you can rewrite the stiffness matrix as

$[\begin{matrix} σ_{11} \\ σ_{12} \\ σ_{21} \\ σ_{22} \end{matrix}] = [\begin{matrix} \frac{E (1 - ν)}{(1 + ν) (1 - 2 ν)} & 0 & 0 & \frac{E ν}{(1 + ν) (1 - 2 ν)} \\ 0 & \frac{E}{2 (1 + ν)} & \frac{E}{2 (1 + ν)} & 0 \\ 0 & \frac{E}{2 (1 + ν)} & \frac{E}{2 (1 + ν)} & 0 \\ \frac{E ν}{(1 + ν) (1 - 2 ν)} & 0 & 0 & \frac{E (1 - ν)}{(1 + ν) (1 - 2 ν)} \end{matrix}] \nabla u = [\begin{matrix} 2 μ + λ & 0 & 0 & λ \\ 0 & μ & μ & 0 \\ 0 & μ & μ & 0 \\ λ & 0 & 0 & 2 μ + λ \end{matrix}] \nabla u$

Axisymmetric Analysis

Axisymmetric analysis speeds up simulations by simplifying 3-D solids using their symmetry around the axis of rotation and analyzing only the 2-D axisymmetric section. Use polar coordinates r,θ,z for radial, circumferential, and axial directions. If z is the axis of rotation, the symmetry around the z-axis means that the stress components are independent of the θ coordinate. The stress equilibrium equations for an axisymmetric structural mechanics are as follows:

$\begin{array}{l} \frac{\partial σ_{r r}}{\partial r} + \frac{\partial σ_{r z}}{\partial z} + \frac{1}{r} (σ_{r r} - σ_{θ θ}) + f_{r} = 0 \\ \frac{\partial τ_{r z}}{\partial r} + \frac{\partial σ_{z z}}{\partial z} + \frac{τ_{r z}}{r} + f_{z} = 0 \end{array}$

τ is the shear stress, and γ is the shear strain. Assuming isotropic conditions, the stress-strain relation can be written as follows:

$(\begin{matrix} σ_{r} \\ σ_{z} \\ \begin{array}{l} σ_{θ} \\ τ_{r z} \end{array} \end{matrix}) = \frac{E}{(1 + ν) (1 - 2 ν)} (\begin{matrix} 1 - ν & ν & 0 & 0 \\ ν & 1 - ν & 0 & 0 \\ 0 & 0 & 1 - ν & 0 \\ 0 & 0 & 0 & \frac{1 - 2 ν}{2} \end{matrix}) (\begin{matrix} ε_{r} \\ ε_{z} \\ \begin{array}{l} ε_{θ} \\ γ_{r z} \end{array} \end{matrix})$