Solve the linear elasticity problem by Fenics

Question

To better understand fenics, I try to solve the linear elasticity problem: a sphere with a spherical cavity inside, subjected to uniform pressure on the outer boundary of the sphere. The variational formulation of this problem is:

$$\int_V \boldsymbol{\sigma(u)} : \boldsymbol{\varepsilon(v)} \mathrm{d}V = \int_{\partial V} \mathbf{t}\cdot\boldsymbol{v} \mathrm{d}S$$

The minimal working sample code:

from mshr import Sphere, generate_mesh
from dolfin import *

# Radius of outer and inner sphere
oradius, iradius = 5., 1.
tol = 1e-3

#Material parameters
nu = 0.3
mu = 1.0
Young = 2.*mu*(1.+nu)
lmbda = 2.*mu*nu/(1.-2.*nu)

# Geometry
outer_sphere = Sphere(Point(0., 0., 0.), oradius)
inner_sphere = Sphere(Point(0., 0., 0.), iradius)
g3d = outer_sphere - inner_sphere
mesh = generate_mesh(g3d, 8)

#Define the outer sphere boundary
class OuterBoundary(SubDomain):
    def inside(self, x, on_boundary):
        r2 = sqrt(x[0]**2 + x[1]**2 + x[2]**2)
        return on_boundary and abs(r2-oradius) < tol

#Define the inner sphere boundary
class InnerBoundary(SubDomain):
    def inside(self, x, on_boundary):
        r2 = sqrt(x[0]**2 + x[1]**2 + x[2]**2)
        return on_boundary and abs(r2-iradius) < tol

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
outer_boundary = OuterBoundary()
inner_boundary = InnerBoundary()

outer_boundary.mark(boundary_parts, 1)
outer_boundary.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]

#Function Space
VV = VectorFunctionSpace(mesh, "Lagrange", 2, 3) # displacement space

u = TrialFunction(VV)
v = TestFunction(VV)

eps_u = sym(grad(u))
eps_v = sym(grad(v))
sig_u = lmbda*tr(eps_u)*delta + 2.*mu*eps_u

# Traction vector [1., 1., 1.] applied on the outer boundary
t = Constant((1., 1., 1.))

F = inner(sig_u, eps_v)*dx - inner(t, v)*ds(1) 
a, L = lhs(F), rhs(F)

u_h = Function(VV)

A, b = assemble_system(a, L)

solve(A, u_h.vector(), b, 'lu')

ux, uy, uz = u_h.split()

plot(ux, mode = "displacement", interactive=True, wire_frame=True, axes=True,rescale=False,scalarbar=True)

However, the above code show a null displacement field ($u_x$). Everything seems ok to me. Could you please show me where I'm wrong?

Thanks!

Answer 1

Hi, the problem is that you are marking outer boundary on boundary_parts twice, overwriting the previous step.

Modify:

outer_boundary.mark(boundary_parts, 1)
outer_boundary.mark(boundary_parts, 2)

by this:

outer_boundary.mark(boundary_parts, 1)
inner_boundary.mark(boundary_parts, 2)

Comments

Thanks for pointing out this typo. But it's not the problem.

Answer 2

Hi, the first problem with your code is that the vector b is zero (check print b.norm("l2") once the vector is assembled). This is probably because the conditions for facets to be marked as one or two are never met. You can see that even after you have marked the boundary_parts, all the facet values are still zero. I will let you figure this out and move to the more serious problem of your code. The system as you assemble it is singular. The matrix has 6 zero eigenvalues and any solution is only determined up to the 3 translations and 3 rotations that are the eigenvectors of the zero eigenvalues. To fix this lack of
uniqueness consider

from mshr import Sphere, generate_mesh
from dolfin import *

# Radius of outer and inner sphere
oradius, iradius = 5., 1.
tol = 1e-3

#Material parameters
nu = 0.3
mu = 1.0
Young = 2.*mu*(1.+nu)
lmbda = 2.*mu*nu/(1.-2.*nu)

# Geometry
outer_sphere = Sphere(Point(0., 0., 0.), oradius)
inner_sphere = Sphere(Point(0., 0., 0.), iradius)
g3d = outer_sphere - inner_sphere
mesh = generate_mesh(g3d, 8)

#Function Space
VV = VectorFunctionSpace(mesh, "Lagrange", 1, 3) # displacement space
Q = VectorFunctionSpace(mesh, 'R', 0, 6)
W = MixedFunctionSpace([VV, Q])

u, p = TrialFunctions(W)
v, q = TestFunctions(W)

eps_u = sym(grad(u))
eps_v = sym(grad(v))
sig_u = lmbda*tr(eps_u)*Identity(3) + 2.*mu*eps_u

t = Constant((1, 1, 1))

# Nullspace of rigid motions
# Translation
Z_transl = [Constant((1, 0, 0)), Constant((0, 1, 0)), Constant((0, 0, 1))]
# Rotations
Z_rot = [Expression(('0', 'x[2]', '-x[1]')),
         Expression(('-x[2]', '0', 'x[0]')),
         Expression(('x[1]', '-x[0]', '0'))]
# All
Z = Z_transl + Z_rot

a = inner(sig_u, eps_v)*dx\
   -sum(p[i]*inner(v, Z[i])*dx for i in range(len(Z)))\
   -sum(q[i]*inner(u, Z[i])*dx for i in range(len(Z))) 

L = inner(t, v)*ds

w_h = Function(W)

A, b = assemble_system(a, L)
print b.norm('l2')

solve(A, w_h.vector(), b, 'lu')

uh, ph = w_h.split(deepcopy=True) 

plot(uh, mode = "displacement", interactive=True, wire_frame=True, axes=True,rescale=False,scalarbar=True)

Here 6 Lagrange multipliers are introduced to enforce orthogonality of the unique solution with respect to the given nullspace (Z) formed by the rigid motions. If you do not like the multiplier approach you can pass the nullspace as 6 $l^2$-orthogonal vectors to the Krylov method, see here for how this is done for singular Poisson problem. Let me know if you have troubles making it work. Also, please tag the question with singular elasticity tag.

Comments

Thank you for the detailed response! I will let you know if it works.

---

It worked like a charm! So I have to block the rotation to make the problem well-posed. But still not find the issue of boundary.mark().

---

Both translations and rotations must be removed. As far as marking goes, consider the following

from mshr import Sphere, generate_mesh
from dolfin import *
import numpy as np

# Radius of outer and inner sphere
oradius, iradius = 5., 1.

# Geometry
outer_sphere = Sphere(Point(0., 0., 0.), oradius)
inner_sphere = Sphere(Point(0., 0., 0.), iradius)
g3d = outer_sphere - inner_sphere
mesh = generate_mesh(g3d, 8)

facet_f = FacetFunction('size_t', mesh, 0)
# Mark all external facets
DomainBoundary().mark(facet_f, 1)
# Now mark the inner boundary by checking that radius is smaller than radius
# of mid sphere
center = Point(0., 0., 0.)
mid = 0.5*(oradius - iradius)
for facet in SubsetIterator(facet_f, 1):
    if facet.midpoint().distance(center) < mid:
        facet_f[facet] = 2

# See that both values are now present
print len(np.where(facet_f.array() == 1)[0])
print len(np.where(facet_f.array() == 2)[0])