Error on the deformation tensor

I/O
1

Hello,

while trying to solve a linear elasticity problem where pure bending condition is applied to the beam as follows

g1024

which should produce values for eps_xy and eps_yy equal to zero.
However, while plotting the deformation tensor i obtain the following results for the various components:

Here is the code used for the solution of the linear elasticity problem:

from dolfin import *
import ufl as uf

# --- Mesh
L, H = 1, 0.1
nel = 75
mesh = RectangleMesh(Point(0, 0), Point(L, H), 100, 10, 'crossed')
ndim = mesh.topology().dim()

# --- Material Parameter
E = Constant(70e3)
nu = Constant(0.)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
k = lmbda + mu*2/3.

# --- Function Spaces
V = VectorElement("Lagrange", mesh.ufl_cell(), 1)
V_u = FunctionSpace(mesh, V)

u = Function(V_u)
du = TrialFunction(V_u)
v = TestFunction(V_u)
u.rename('u', 'u')

S3 = TensorFunctionSpace(mesh, 'DG', 0)
Eps = Function(S3)
Eps.rename('E', 'E')

# ---Boundaries
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
left = Left()

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], L)
right = Right()

# ---Measures
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 3)

# ---Displacement BC
u_R = Expression("t*(0.5-x[1]/H)", H=H, t=0.05, degree=1)
bcu_l = DirichletBC(V_u, Constant((0., 0.)), boundaries, 1)
bcu_r = DirichletBC(V_u.sub(0), u_R, boundaries, 3)
bcu = [bcu_r, bcu_l]

# ---Definitions
def eps(u):
    return sym(grad(u))

def sigma(u):
    return 2. * mu * (eps(u)) + lmbda * tr(eps(u)) * Identity(ndim)

# --- Energy functional
psi = 0.5 * inner(sigma(u), eps(u))*dx
E_u = derivative(psi, u, v)
E_du = uf.replace(E_u, {u: du})
problem_u = LinearVariationalProblem(lhs(E_du), rhs(E_du), u, bcu)
solver_u = LinearVariationalSolver(problem_u)

# --- Solver
solver_u.solve()
File('u.pvd') << u
Eps.assign(project(eps(u), S3))
File("E.pvd") << Eps

Am i making mistakes in the solution of the elasticity problem or is this problem tied to the projection of the deformation tensor for the output?

Thanks in advance

2

These strains are a lot smaller than \epsilon_{xx}, and the error decreases with refinement. In general, you expect some discretization error when your displacement function space does not contain the exact solution. In the case of pure bending, your exact solution involves quadratic variation of the vertical displacement u_y. You can see that the error in \epsilon_{xy} and \epsilon_{yy} becomes negligible when you use quadratic elements for the displacement:

#V = VectorElement("Lagrange", mesh.ufl_cell(), 1)
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
2 Likes