Issue with applying Dirichlet boundary condition (Linear Elasticity)

1

Hi
I have a unit cube and I want to have only 1 element along each edge. I want to implement a bi-axial tension simulation. With that, I want to fix the left and bottom faces and apply equal normal stretch (Non-Zero Dirichlet boundary condition) on the right and top faces. However it seems like that the boundary conditions are not applied on the left and right boundaries as shown here:
222
However I am expecting something like this:

333

Does anybody know the right way to apply such Dirichlet boundary conditions? Here is the code:

from fenics import *

mu = 1
lambda_ = 1.25
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)
V = VectorFunctionSpace(mesh, 'CG', 1)

# Define boundary condition
tol = 1E-14

def LEFT(x, on_boundary):
    return on_boundary and x[0] < tol

def RIGHT(x, on_boundary):
    return on_boundary and abs(x[0] - 1.0) < tol

def BOTTOM(x, on_boundary):
    return on_boundary and x[1] < tol

def TOP(x, on_boundary):
    return on_boundary and abs(x[1] - 1.0) < tol

bc_left = DirichletBC(V, Constant((0, 0, 0)), LEFT)
bc_right = DirichletBC(V, Constant((0.1, 0, 0)), RIGHT)
bc_bottom = DirichletBC(V, Constant((0, 0, 0)), BOTTOM)
bc_top = DirichletBC(V, Constant((0, 0.1, 0)), TOP)

bc_all = [bc_left,bc_right, bc_bottom,  bc_top]
# Define strain and stress
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc_all)

File("displacement.pvd") << u
2

You cannot prescribe both \boldsymbol{U}=0.1\boldsymbol{e}_y and \boldsymbol{U}=0.1\boldsymbol{e}_x both on right and top, since the nodes belonging both to right and top cannot fulfill simultaneously both conditions. What you probably want is to impose \boldsymbol{U}=0.1x\boldsymbol{e}_x+0.1y\boldsymbol{e}_y everywhere on the boundary which will yield a uniform strain E_{xx}=0.1, E_{yy}=0.1

3

Thanks for looking into this issue. I found out this problem is not related to existence of some nodes on the right and top faces simultaneously. It seems like the order of the Dirichlet boundary conditions is important. For example if I use:

bc_all = [bc_right,bc_top,bc_left,bc_bottom]

instead of :

bc_all = [bc_left,bc_right,bc_bottom,bc_top]

The results are reasonable. Here is the displacement:

444

It is strange why the order of the boundary conditions has influence on the results. Anyways, that was a fix to this issue.

4

Hi,
I think that the issue is still there since what you show is only a Y displacement, there is no X displacement. Depending on the order of the boundary conditions, the conflicting bc will be satisfied for some node but not the other. If you try with a more refined mesh, say 5x5x5, you will see that you don’t have a proper strain state. You should use instead something like:

def border(x, on_boundary):
    return on_boundary
bc_all = DirichletBC(V, Expression(("0.1*x[0]","0.1*x[1]","0"), degree=1), border)
1 Like
5

Hi
You are absolutely right. I did not check the X displacement. After checking that I noticed there is not displacement in X direction. Defining the boundary condition based on what you suggested makes sense.
However after looking into the results, it seems like the bottom and left boundaries are not fixed:
55
Which is different from the second figure I provided above.

6

What you want to do will not yield a uniform strain state. What type of mathematical conditions do you want to precisely write, we cannot help you based on a picture only.

7

What I am trying to do is simulation of a bi-axial stretch test. From the experimental point of view, the left and bottom faces should be fixed and displacement should be applied on the right and top faces. I am trying to find a way to apply these boundary conditions. It is possible to do it using a commercial software (e.g. ANSYS) that I showed in the second figure. In this case you can see that the bottom and left faces are fixed (Zero displacement). However I was not been able to reproduce such displacement contour in FEniCS.

8

Again, this will not enforce a uniform stretch state if that’s what you are trying to achieve.
However, if you really want to do that how about imposing

bc_all = [DirichletBC(V, Expression(("0.1*x[0]","0.1*x[1]","0"), degree=1), TOP),
          DirichletBC(V, Expression(("0.1*x[0]","0.1*x[1]","0"), degree=1), RIGHT),
          DirichletBC(V, Constant((0, 0, 0)), BOTTOM),
          DirichletBC(V, Constant((0, 0, 0)), LEFT)]