Retrieving imposed traction boundary in a linear elastic problem

1

Dear FEniCS users,

I define a simple linear elasticity problem with Neumann boundary conditions. I would like to retrieve the imposed traction on the boundary, but the one I get is hardly similar to the initial one. Here is a minimal working example.

import matplotlib.pyplot as plt
from dolfin import *
from mshr import *

# Creating mesh

domain = Circle(Point(0.0,0.0),1)
mesh = generate_mesh(domain,10)
V = VectorFunctionSpace(mesh,'CG',1)

# Assigning constants
E = 30
nu = 0.3
m = E/(2*(1+nu))
l = (E*nu)/((1+nu)*(1-2*nu))

# Define boundary as subdomain
class boundary(SubDomain):
    def inside(self,x,on_boundary):
        return(on_boundary)
Boundary = boundary()
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
Boundary.mark(boundaries,1)
ds = ds(domain=mesh, subdomain_data=boundaries)


def epsilon(u):
    return(0.5*(nabla_grad(u) + nabla_grad(u).T))
def sigma(u,lambda_=1.,mu=1.):
    return(lambda_*div(u)*Identity(d) + 2*mu*epsilon(u))

u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0.0,0.0))
Tt = Constant((1.0,1.0))
a = inner(sigma(u,l,m), epsilon(v))*dx
L = dot(f, v)*dx + dot(Tt,v)*ds(1)
u = Function(V)
solve(a == L, u)

Tt_trial, Tt_test = TrialFunction(V), TestFunction(V)
a = inner(Tt_trial, Tt_test)*ds(1)
L = inner(Tt,       Tt_test)*ds(1)
A = assemble(a, keep_diagonal=True)
A.ident_zeros() 
b = assemble(L)
T_true = Function(V)
solve(A, T_true.vector(), b, 'cg', 'amg') 

n = FacetNormal(mesh)
u_trial, u_test = TrialFunction(V), TestFunction(V)
a = inner(u_trial, u_test)*ds(1)
L = inner(dot(sigma(u,l,m),n),       u_test)*ds(1)
A = assemble(a, keep_diagonal=True)
A.ident_zeros()
b = assemble(L)
T = Function(V)

solve(A, T.vector(), b, 'cg', 'amg') 
plot(T_true, title="Imposed traction")
plot(T, title="Retrieved traction")

Running this I get the following plot :
Retrieved_traction

Could anyone please show me where my error lies ?
Thanks in advance.

2

Please note that the traction, i.e. sigma(u) does not lie in the same function space as the deformation.
You should project your traction into a function space of lower degree and continuity, i.e,


n = FacetNormal(mesh)
VDG = VectorFunctionSpace(mesh, "DG", 0)
u_trial, u_test = TrialFunction(VDG), TestFunction(VDG)
a = inner(u_trial, u_test)*ds(1)
L = inner(dot(sigma(u,l,m),n),       u_test)*ds(1)
A = assemble(a, keep_diagonal=True)
A.ident_zeros()
b = assemble(L)
T = Function(VDG)

Secondly, you cannot solve a linear elasticity problem without Dirichlet boundary condition’s as it is not well defined. Then you would have to do a null space removal of the problem, see for instance: Bitbucket

3

Thank you Dokken for your fast response. As for the Dirichlet boundary condition it is no problem for me to add some. Here is the code I would use :

domain = Circle(Point(0.0,0.0),1)
mesh = generate_mesh(domain,10)

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

E = 30
nu = 0.3

m = E/(2*(1+nu))
l= (E*nu)/((1+nu)*(1-2*nu))

class boundary(SubDomain):
    def inside(self,x,on_boundary):
        return(on_boundary)
    
class top(SubDomain):
    def inside(self,x,on_boundary):
        return(on_boundary and (x[1]>-0.6))

class bottom(SubDomain):
    def inside(self,x,on_boundary):
        return(on_boundary and (x[1]<=-0.6))

Boundary = boundary()
Top = top()
Bottom = bottom()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
Top.mark(boundaries,1)
Bottom.mark(boundaries,2)

Tt = Constant((0.0, 1.0))
ds = ds(domain=mesh, subdomain_data=boundaries)

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

def sigma(u,lambda_=1.,mu=1.):
    return(lambda_*div(u)*Identity(d) + 2*mu*epsilon(u))

bc = DirichletBC(V,(0.0,0.0),boundaries,2)

u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0.0,0.0))

a = inner(sigma(u,l,m), epsilon(v))*dx
L = dot(f, v)*dx + dot(Tt,v)*ds(1)

u = Function(V)
solve(a == L, u, bc)

When I follow your suggestion and project the traction T on a discontinuous galerkin space, I get two issues : 1) the traction is not defined only on the boundary nodes, 2) the tractions get spurious at the bottom of the sphere, while they seem rather continuous at the top, as shown in the following picture.

T_DG

Do you think there might be a way to get better approximations of the tractions with the computed displacements ? Thanks again for your amazing response.

4

I would suggest using Paraview for the visualization, as it preserves the piecewise constants of the DG 0 space.
If you do the following:

File("traction.pvd") << T
File("u.pvd") <<u

you can visualize it as

traction
traction1280×720 185 KB

where the left figure is the traction, the right the undeformed and deformed mesh.

5

Thank you Dokken for your responses ! Best.