Help on solving 2D plane strain problem

Question

Hi there,

Is it possible to solve 2D plane strain problem in a 3D mesh?

Based on the Elasticity demo, I first calculate the stress manually by setting each strain in Y direction to be 0. The code is as follows:

def strain(v):
    strain = as_tensor([\
        [v[0].dx(0), v[0].dx(1), v[0].dx(2)],\
        [v[1].dx(0), v[1].dx(1), v[1].dx(2)],\
        [v[2].dx(0), v[2].dx(1), v[2].dx(2)]\
        ])
    return strain

def sigma(strain):
    sigma = as_tensor([ \
        [2.0*mu*strain[0,0] + lmbda*(strain[0,0] + strain[2,2]), 0.0,                               mu*(strain[0,2] + strain[2,0])],\
        [0.0,                                                    lmbda*(strain[0,0] + strain[2,2]), 0.0],\
        [mu*(strain[2,0] + strain[0,2]),                         0.0,                               2.0*mu*strain[2,2] + lmbda*(strain[0,0] + strain[2,2])] \
        ])
    return sigma

In order to solve the problem, I guess we also need to constrain the displace in Y direction anywhere to be 0. However, I'm not sure how to do this? Could you please give me some clue? It'll be great if you've some experience in 2D plane strain problem.

The complete code is also attached. Be aware that it's not working yet.

from dolfin import *

# Mesh and function space
mesh = UnitCubeMesh(10,10,10)
V = VectorFunctionSpace(mesh, "CG", 1)

# Constants
E  = 10.0
nu = 0.3
mu    = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

# Weakform
u = Function(V)
v = TestFunction(V)
f = Constant((0.0, 0.0, 0.0))

def strain(v):
    strain = as_tensor([\
        [v[0].dx(0), v[0].dx(1), v[0].dx(2)],\
        [v[1].dx(0), v[1].dx(1), v[1].dx(2)],\
        [v[2].dx(0), v[2].dx(1), v[2].dx(2)]\
        ])
    return strain

def sigma(strain):
    sigma = as_tensor([ \
        [2.0*mu*strain[0,0] + lmbda*(strain[0,0] + strain[2,2]), 0.0,                               mu*(strain[0,2] + strain[2,0])],\
        [0.0,                                                    lmbda*(strain[0,0] + strain[2,2]), 0.0],\
        [mu*(strain[2,0] + strain[0,2]),                         0.0,                               2.0*mu*strain[2,2] + lmbda*(strain[0,0] + strain[2,2])] \
        ])
    return sigma

F = inner(sigma(strain(u)), grad(v))*dx - inner(f, v)*dx

# Boundary condition
# Bottom boundary
class Bottom_surf(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
         # tolerance for coordinate comparisons
        return  abs(x[2]) < tol

bottom_posi = Bottom_surf()

class Bottom_disp(Expression):
    def eval(self, values, x):
        values[0] = .0
        values[1] = .0
        values[2] = .0

    def value_shape(self):
        return (3,)

bottom_expr = Bottom_disp()

# Top boundary
class Top_surf(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
         # tolerance for coordinate comparisons
        return  abs(x[2] - 1.0) < tol

top_posi = Top_surf()

class Top_disp(Expression):
    def eval(self, values, x):
        values[0] = .0
        values[1] = .0
        values[2] = -0.5

    def value_shape(self):
        return (3,)

top_expr = Top_disp()

bc_b = DirichletBC(V, bottom_expr, bottom_posi)
bc_t = DirichletBC(V, top_expr, top_posi)
bcs = [bc_t, bc_b]

# Compute solution
solve(F == 0, u, bcs)

# Plot solution
plot(u, mode = "displacement", interactive=True)

Answer 1

In fact I think you can use the usual definition of the strain/stress

V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)
eps = sym(grad(u))
Id = Identity(u.geometric_dimension())
sigma = 2.0*mu*eps + lmbda*tr(eps)*Id

and block the Y component of the displacement V.sub(1) by using the following BC

def domain(x, on_boundary):
   return True

bc = DirichletBC(V.sub(1), Constant(0.0), domain)

Comments

The alternative is to use a mixed function space, but it is probably more messy. I use that strategy for plane stress.