This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Bogus output for uniaxial stress (Neumann) boundary conditions

+1 vote

I would like to set constant stress (Neumann) boundary conditions. As a simple 2D test case, I solve plane stress force balance (divergence of sigma = 0), with tractions of P on the right and left and 0 on the top and bottom of the unit square.

In weak form (the displacement term makes the solution unique):

Vv = VectorFunctionSpace(mesh, "Lagrange", 1)
tau = TestFunction(Vv) 
u = TrialFunction(Vv)
epsilon = sym(nabla_grad(u))
sigma = E/(1+nu)* (epsilon) + E*nu/((1-nu**2)) * tr(epsilon)*Identity(d)
a = inner(sigma, sym(nabla_grad(tau)))*dx + inner(u,tau)*dx 
L = dot(P,tau)*ds

So the expected result is that sigma_xx = P everywhere and sigma_yy= sigma_xy = 0 everywhere. However, the output has strange behavior near the boundaries:
enter image description here

What am I doing wrong?

Here is the minimal working code:

from dolfin import *
import numpy as np

'''Check method for sheet with uniaxial stress boundary.'''

##############################
# System Params ##############
############################## 
E , nu = 1.0, 0.4          #J/m^3, unitless

############################## 
# Mesh and Function Space ####
##############################
mesh =UnitSquareMesh(100, 100)
Vv = VectorFunctionSpace(mesh, "Lagrange", 1)
Vf = FunctionSpace(mesh, "Lagrange", 1)
Vt = TensorFunctionSpace(mesh, "Lagrange", 1)

##############################
# Define variational problem #
##############################
print 'Defining test and trial functions...'
tau = TestFunction(Vv)
u = TrialFunction(Vv)

# Neumann BC
P = Expression(('E*0.05*(x[0]-xc)','0.0'),E=E,xc=0.5)

# functions for variational problem
d = u.geometric_dimension()
epsilon = sym(nabla_grad(u))
# PLANE STRESS
sigma = E/(1+nu)* (epsilon) + E*nu/((1-nu**2)) * tr(epsilon)*Identity(d)

## Solve force balance
a = inner(sigma, sym(nabla_grad(tau)))*dx + inner(u,tau)*dx 
L = dot(P,tau)*ds
u = Function(Vv)
solve(a == L, u) 

## View the result
sigma = E/(1+nu)* (sym(nabla_grad(u))) + E*nu/((1-nu**2)) * tr(sym(nabla_grad(u)))*Identity(d)
s_V = project( sigma, Vt)
sxx = project(s_V[0,0], Vf);           sxxv = sxx.vector().array()
sxy = project(s_V[0,1], Vf);           sxyv = sxy.vector().array()
syx = project(s_V[1,0], Vf);           syxv = syx.vector().array()
syy = project(s_V[1,1], Vf);           syyv = syy.vector().array()

plot(sxx, interactive=True)
plot(syy, interactive=True)
asked Oct 7, 2015 by npmitchell FEniCS Novice (600 points)

I think that:

a = inner(sigma, sym(nabla_grad(tau)))*dx + inner(u,tau)*dx

would be:

a = inner(sigma, sym(nabla_grad(tau)))*dx

?

commented Oct 7, 2015 by newuser FEniCS Novice (650 points)

newuser, a term like inner(u,tau)*dx is necessary to make the displacement solution unique -- otherwise u is defined only up to translations in x and y. For instance, see this example:
http://fenicsproject.org/documentation/dolfin/dev/python/demo/documented/neumann-poisson/python/documentation.html
However, it does seem this term somehow messes things up.

commented Oct 7, 2015 by npmitchell FEniCS Novice (600 points)
edited Oct 7, 2015 by npmitchell

Ok, so it's a penalty term to remove rigid body motion if I'm right. Thanks!

commented Oct 7, 2015 by newuser FEniCS Novice (650 points)

Yes, but evidently the term I wrote is somehow incorrect, since when I use your answer to write the problem correctly without the term, I get the right stress field and a displacement field that is somewhat arbitrary, but when I add this term to account for the translation, I get the wrong stress field...

commented Oct 7, 2015 by npmitchell FEniCS Novice (600 points)

2 Answers

0 votes
 
Best answer

npmitchell, the following code seems to work:

from mshr import *
from dolfin import *

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = DOLFIN_EPS
        return on_boundary and abs(x[0]) < tol

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = DOLFIN_EPS
        return on_boundary and abs(x[0] - 1) < tol

E , nu = 1.0, 0.4          #J/m^3, unitless

mesh = UnitSquareMesh(100, 100)

Vv = VectorFunctionSpace(mesh, "Lagrange", 1)
Vf = FunctionSpace(mesh, "Lagrange", 1)
Vt = TensorFunctionSpace(mesh, "Lagrange", 1)
Q = VectorFunctionSpace(mesh, 'R', 0, 3)
W = MixedFunctionSpace([Vv, Q])

u, p = TrialFunctions(W)
tau, q = TestFunctions(W)

left_boundary = LeftBoundary()        
right_boundary = RightBoundary()        

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
left_boundary.mark(boundary_parts, 1)
right_boundary.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]  

P1 = Constant((-1, 0))
P2 = Constant((1, 0))

d = u.geometric_dimension()
epsilon = sym(grad(u))
sigma = E/(1+nu)* (epsilon) + E*nu/((1-nu**2)) * tr(epsilon)*Identity(d)

trans = [Constant((1, 0)), Constant((0, 1))]
rota = [Expression(('-x[1]', 'x[0]'))]
n = trans+rota

## Solve force balance
a = inner(sigma, sym(grad(tau)))*dx - sum(p[i]*inner(tau, n[i])*dx for i in range(len(n)))\
   -sum(q[i]*inner(u, n[i])*dx for i in range(len(n)))
L = dot(P1,tau)*ds(1) + dot(P2,tau)*ds(2)

A, b = assemble_system(a, L)

u = Function(W)
solve(A, u.vector(), b)

uh, nh = u.split()

## View the result
sigma = E/(1+nu)* (sym(grad(uh))) + E*nu/((1-nu**2)) * tr(sym(grad(uh)))*Identity(d)
s_V = project( sigma, Vt)
sxx = project(s_V[0,0], Vf);           sxxv = sxx.vector().array()
sxy = project(s_V[0,1], Vf);           sxyv = sxy.vector().array()
syx = project(s_V[1,0], Vf);           syxv = syx.vector().array()
syy = project(s_V[1,1], Vf);           syyv = syy.vector().array()

plot(sxx, interactive=True)
plot(syy, interactive=True)

I introduced the Langrange multiplier to assure the uniqueness of $u$ (block both translation and rotation of rigid body motion).

  1. I found this interesting discussion. You can build a null space and pass it to PETSc solver in Fenics!

  2. Thanks to MiroK

answered Oct 7, 2015 by newuser FEniCS Novice (650 points)
selected Oct 8, 2015 by npmitchell
+1 vote

Hi, I have the same problem but for another geometry (a sphere). Personally, I would write everything in details to figure out the error...

Here is my solution based on your code, I modified it a bit (use the fact that this problem is symmetry, so I solved on the half of square):

from mshr import *
from dolfin import *

'''Check method for sheet with uniaxial stress boundary.'''

##############################
# System Params ##############
############################## 
E , nu = 1.0, 0.4          #J/m^3, unitless

############################## 
# Mesh and Function Space ####
##############################
mesh = UnitSquareMesh(100, 100)
Vv = VectorFunctionSpace(mesh, 'Lagrange', 1)
Vf = FunctionSpace(mesh, 'Lagrange', 1)
Vt = TensorFunctionSpace(mesh, 'Lagrange', 1)

VX, VY = Vv.split()

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-7
        return on_boundary and abs(x[0]) < tol

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-7
        return on_boundary and abs(x[0] - 1) < tol

left_boundary = LeftBoundary()        
right_boundary = RightBoundary()        

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
right_boundary.mark(boundary_parts, 1)
ds = Measure("ds")[boundary_parts]   

# The displacement $u_x$ is blocked
bc = DirichletBC(VX, Constant(0.), left_boundary)
##############################
# Define variational problem #
##############################
print 'Defining test and trial functions...'
tau = TestFunction(Vv)
u = TrialFunction(Vv)

# Neumann BC
P = Constant((1., 0.))

# functions for variational problem
d = u.geometric_dimension()
epsilon = sym(grad(u))
# PLANE STRESS
sigma = E/(1+nu)*(epsilon) + E*nu/((1-nu**2)) * tr(epsilon)*Identity(d)

## Solve force balance
a = inner(sigma, sym(nabla_grad(tau)))*dx
L = dot(P,tau)*ds(1)

u_h = Function(Vv)

solve(a==L, u_h, bc)

## View the result
sigma = E/(1+nu)* (sym(nabla_grad(u_h))) + E*nu/((1-nu**2)) * tr(sym(nabla_grad(u_h)))*Identity(d)

s_V = project( sigma, Vt)
sxx = project(s_V[0,0], Vf);           sxxv = sxx.vector().array()
sxy = project(s_V[0,1], Vf);           sxyv = sxy.vector().array()
syx = project(s_V[1,0], Vf);           syxv = syx.vector().array()
syy = project(s_V[1,1], Vf);           syyv = syy.vector().array()

plot(sxx, interactive=True, rescale=False,scalarbar=True)
plot(syy, interactive=True, rescale=False,scalarbar=True)
plot(sxy, interactive=True, rescale=False,scalarbar=True)
answered Oct 7, 2015 by newuser FEniCS Novice (650 points)

newuser: This is helpful, but I can't have any Dirichlet conditions in my actual application (for which this is a test). A trivial change to your suggestion overcomes this. Removing the Dirichlet condition,

right_boundary.mark(boundary_parts, 1)
left_boundary.mark(boundary_parts, 2)
...
P = Expression(('2.0*(x[0]-xc)','0.0'),xc=0.5)
...
a = inner(sigma, sym(nabla_grad(tau)))*dx
L = dot(P,tau)*ds(1) + dot(P,tau)*ds(2)
u = Function(Vv)
solve(a == L, u) 
...

The remaining issue is ensuring the uniqueness of u.

commented Oct 7, 2015 by npmitchell FEniCS Novice (600 points)
edited Oct 7, 2015 by npmitchell
...