Linear Elasticity Problem

Linear Elasticity Problem#

Analysis of deformable structures is essential in engineering, with the equations of linear elasticity being fundamental in this field. In this section, we present the variational formulation of linear elasticity equations using the principle of virtual work.

The PDE problem#

The governing equations for small elastic deformations of a body \( \Omega \) can be expressed as:

\[\begin{split} \begin{align} -\nabla \cdot \sigma(u) &= f & \text{in } \Omega \\ \sigma(u) &= \kappa \text{tr}(\epsilon(u))I + 2 \mu \epsilon(u) \end{align} \end{split}\]

where \( \sigma \) is the stress tensor, \( f \) represents the body force per unit volume, \( \kappa \) and \( \mu \) are Lamé’s elasticity parameters for the material, \( I \) denotes the identity tensor, and \( \epsilon \) is the symmetric strain tensor. The displacement vector field is denoted by \( u \).

By substituting \( \epsilon(u) \) into \( \sigma \), we obtain:

\[ \sigma(u) = \kappa (\nabla \cdot u)I + \mu(\nabla u + (\nabla u)^T) \]

The Variational Formulation#

The variational formulation of the linear elasticity equations involves forming the inner product of the PDE with a vector test function \( v \in \hat{V} \) and integrating over the domain \( \Omega \). This yields:

\[ \int_{\Omega} \sigma : \nabla v \, \mathrm{d} x = \int_{\Omega} f \cdot v \, \mathrm{d} x \]

Integrating the term \( \nabla \cdot \sigma \cdot v \) by parts, considering boundary conditions, we obtain:

\[ \int_{\Omega} \sigma : \nabla v \, \mathrm{d} x = \int_{\Omega} f \cdot v \, \mathrm{d} x + \int_{\partial \Omega_T} T \cdot v \, \mathrm{d} s \]

where \( T \) represents the traction vector on the part \( \partial \Omega_T \) of the boundary where it’s prescribed.

This leads to the variational formulation: Find \( u \in V \) such that

\[ a(u, v) = L(v) \quad \forall v \in \hat{V} \]

where

\[\begin{split} \begin{align} a(u, v) &= \int_{\Omega} \sigma(u) : \nabla v \, \mathrm{d} x \\ L(v) &= \int_{\Omega} f \cdot v \, \mathrm{d} x + \int_{\partial \Omega_T} T \cdot v \, \mathrm{d} s \end{align} \end{split}\]

This formulation can be alternatively expressed as:

\[ a(u, v) = \int_{\Omega} \sigma(u) : \epsilon(v) \, \mathrm{d} x \]

where \( \epsilon(v) = \frac{1}{2} (\nabla v + (\nabla v)^T) \) is the symmetric strain tensor.

This variational formulation is essential for solving linear elasticity problems numerically using methods like the finite element method (FEM).

Formal Model#

from sympde.expr import BilinearForm, LinearForm, integral
from sympde.expr     import find, EssentialBC, Norm, SemiNorm
from sympde.topology import VectorFunctionSpace, Cube, element_of
from sympde.calculus import grad, dot, inner, outer, cross, div
from sympde.core import Constant
from sympde.core import Matrix, Vector, Transpose

domain = Cube()

I = Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]], name='I')

kappa = Constant('kappa', is_real=True)
mu    = Constant('mu',    is_real=True)
rho   = Constant('rho',   is_real=True)

epsilon = lambda w: 0.5*(grad(w) + Transpose(grad(w)))
sigma   = lambda w: kappa * div(w) * I + 2 * mu * epsilon(w)

V = VectorFunctionSpace('V', domain)

x,y,z = domain.coordinates

u,v = [element_of(V, name=i) for i in ['u', 'v']]

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 6
      4 from sympde.calculus import grad, dot, inner, outer, cross, div
      5 from sympde.core import Constant
----> 6 from sympde.core import Matrix, Vector, Transpose
      8 domain = Cube()
     10 I = Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]], name='I')

ImportError: cannot import name 'Vector' from 'sympde.core' (/opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages/sympde/core/__init__.py)

# bilinear form
a = BilinearForm((u,v), integral(domain , inner(sigma(u), epsilon(v))))

L = 1
W = 1 #0.2
delta = W / L

#mu = 1
#rho = 1
#kappa = 1.25
g = 0.4 * delta**2

# linear form
f = Vector([0, 0, -rho*g], name='f')
l = LinearForm(v, integral(domain, dot(f,v)))

# Dirichlet boundary conditions
bc = [EssentialBC(u, 0, domain.boundary)]

# Variational problem
equation   = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)

Discretization#

We shall need the discretize function from PsyDAC.

from psydac.api.discretization import discretize

degree = [2,2,2]
ncells = [8,8,8]

# Create computational domain from topological domain
domain_h = discretize(domain, ncells=ncells, comm=None)

# Create discrete spline space
Vh = discretize(V, domain_h, degree=degree)

# Discretize equation
equation_h = discretize(equation, domain_h, [Vh, Vh])

Solving the PDE#

uh = equation_h.solve(mu=1, rho=1, kappa=1.25)
Contents