Elasticity problem pure Neumann boundary condition.

Question

Hi, I am solving Elasticity problem with pure Neumann boundary condition.
Problem definition: Square plate subjected to pressure (2 unit) in y-direction at top and bottom.

But my FEniCS solution is not matching with Abaqus solution.

Abaqus solution: at (0.0, 1.0) u1 = 0.073 and u2 = 0.00275
FEniCS solution at (0.0,1.0) u1 = -0.4302 and u2 = -0.3206

Can anyone point out the mistake to fix this issue.

Thanks & Regards

# test elasticity 
# purpose:
#=======================================
from dolfin import *
mesh = UnitSquareMesh(2,2)
E, nu = Constant(10.0), Constant(0.3)
plot(mesh, title = 'mesh', interactive = True)
ndim = mesh.geometry().dim()
V = VectorFunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
du = Function(V)

mu    = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu )*(1.0-2.0*nu))

def epsilon(v):
      return 0.5*(grad(v)  + grad(v).T)

def sigma(u):
      return 2.0*mu*epsilon(u) + lmbda*tr(epsilon(u))*Identity(2)

class top(SubDomain):
        def inside(self,x,on_boundary):
               tol = 1e-10
               return abs(x[1]-1.0) < tol and on_boundary

class bottom(SubDomain):
          def inside(self,x,on_boundary):
                tol = 1e-10
                return abs(x[1]) < tol and on_boundary
Top = top()
Bottom = bottom()


boundaries = FacetFunction("size_t",mesh)
boundaries.set_all(0)
Top.mark(boundaries, 1)
Bottom.mark(boundaries, 2)
gR = Expression(('0.0','2'))
gF = Expression(('0.0', '-2'))
ds = ds(domain=mesh , subdomain_data=boundaries)

E_du = inner(grad(v),sigma(u))*dx - dot(gR,v)*ds(1) - dot(gF,v)*ds(2)
solve(lhs(E_du)==rhs(E_du), du)
plot(du,interactive = True)
u1, u2 = split(du)
plot(u1, interactive = True)
plot(u2, interactive = True)

Answer 1

Hi, comparing point values for this problem should be done with care. In 2d the Neumann problem in linear elasticity has 3 singular modes, see e.g. the discussion here. That is, for compatible right hand side (the one which vanishes on singular modes) the solution of your weak form is unique up to a linear combination of these modes. The direct solver invoked by solve(lhs(E_du)==rhs(E_du), du) finds some solution but you don't know how the modes were fixed. Consult the singular Poisson demo to see how to get more control over the obtained solution.