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.
The equations governing small elastic deformations of a body can be written as
The variational formulation of (3.20)--(3.22) consists of forming the inner product of (3.20) and a vector test function , where is a vector-valued test function space, and integrating over the domain :
We can now summarize the variational formulation as: find such that
One can show that the inner product of a symmetric tensor and an anti-symmetric tensor vanishes. If we express as a sum of its symmetric and anti-symmetric parts, only the symmetric part will survive in the product since is a symmetric tensor. Thus replacing by the symmetric gradient gives rise to the slightly different variational form
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 with the density of the beam and the acceleration of gravity. The beam is box-shaped with length and has a square cross section of width . We set at the clamped end, . The rest of the boundary is traction free; that is, we set .
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 < 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.
The primary unknown is now a vector field and not a scalar field, so we need to work with a vector function space:
V = VectorFunctionSpace(mesh, 'P', 1)
u = Function(V) we get
u as a vector-valued finite element
function with three components for this 3D problem.
For the boundary condition , 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)).
The gradient and divergence operators now have a prefix
This is strictly not necessary in the present problem, but
recommended in general for vector PDEs arising from continuum mechanics,
if you interpret as a vector in the PDE notation;
see the box about
nabla_grad in the section Variational formulation.
As soon as the displacement
u is computed, we can compute various
stress measures. We will compute the von Mises stress defined as
where is the deviatoric stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d) von_Mises = sqrt(3./2*inner(s, s))
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)