What is the most accurate way to recover the stress tensor?

Question

Hi,

I am trying to solve a simple solid mechanics problem where I need the value of the stress tensor. I currently am trying to recover it using project(), but the results show that the boundary conditions are not being followed. This indicates to me that I'm doing something wrong, either in the way that I set up the problem, or the way that I recover the stress tensor. Read on for the details:

The model is of an isotropic linearly elastic solid, and I'm solving it in static equilibrium. The domain is a 3d rectangular prism which is fixed on one surface (the left face) with a Dirichlet boundary condition, and has stress-free Neumann boundary conditions on all other surfaces.

Currently, after solving for the displacements u, I recover the value of the stress vector on the top surface (surface traction) as follows:

def sigma(v):
    gdim = v.geometric_dimension()
    strain = sym(grad(v))
    I = Identity(gdim)

    return 2.0*mu*strain + lam*tr(strain)*I

# the surface normal of the top should be (0, 0, 1)
surface_traction_top = project(sigma(u) * Constant((0, 0, 1)))

If I sample points from surface_traction_top along the top of my domain, it shows substantial non-zero values near the left edge. I've plotted the normal stress on the top surface from surface_traction_top to illustrate this:

As I understand it, my Neumann boundary condition demands that the surface traction equal a vector of zeros on all boundaries except the left face, which is clamped in place by the Dirichlet boundary condition. This is clearly obeyed for the most part in the above plot, but near that edge, which it shares with the clamped left face, there is substantial deviation.

I have read that in solid mechanics, it is not straightforward to recover the stress tensor in the finite element method, which is what led me to this question. Is there a better way to calculate the stress tensor in FEniCS?

Edit:
For reference, here's my Python code. I tried to trim it down as much as possible, but it's still a bit lengthy:

# Should run on FENICS 1.6.0

# Based on:
    # FENICS elasticity demo by Garth N. Wells
    # Hyperelasticity example by Anders Logg from Oxford University's course, Python in Scientific Computing

from dolfin import *
import numpy as np
import os

import matplotlib
if 'DISPLAY' not in os.environ.keys():
    # allow making figures without X server -- must be run before pyplot is imported
    matplotlib.use('Agg')
import matplotlib.pyplot as plt

# make nicer figures if seaborn is available
try:
    import seaborn as sns
    sns.set_style('white')
    sns.set_context('talk')
except ImportError:
    pass



# specify parameters
####################

# Adjust log level
set_log_level(PROGRESS)

# Turn on optimization
parameters["form_compiler"]["cpp_optimize"] = True

# dimensions of cantilever
length = 1.25e-2
width = 0.5e-2
thickness = 0.1e-2

# Material parameters
rho   = 1450.0       # kg / m^3 density of PVC plastic
mu    = 0.0023 * 1e9 # N  / m^2 Lame parameter mu for PVC plastic
lam = 0.0105 * 1e9 # N  / m^2 Lame parameter lambda for PVC plastic

# number of mesh elements along each axis
nx = 40
ny = 15
nz = 3




# make mesh and define problem
##############################

mesh = BoxMesh(Point(0,0,0), Point(length, width, thickness), nx, ny, nz)

# Define function space for displacement vector
V = VectorFunctionSpace(mesh, "Lagrange", 2)

# gravity
body_force  = Constant((0.0, 0.0, -9.8*rho)) 

# Stress computation
def sigma(v):
    gdim = v.geometric_dimension()
    strain = sym(grad(v))
    I = Identity(gdim)

    return 2.0*mu*strain + lam*tr(strain)*I

bnd_left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
bnd_top = CompiledSubDomain("near(x[2], side) && on_boundary", side = thickness)

# Mark boundaries
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)

bnd_left.mark(boundaries, 1)
bnd_top.mark(boundaries, 2)

# Redefine boundary measure
ds = ds[boundaries]

# Define boundary conditions
# clamp left face
bcs = [DirichletBC(V, (0, 0, 0), boundaries, 1)]




# solve for displacement
########################

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), grad(v))*dx
L = inner(body_force, v)*dx

u = Function(V)

solve(a == L, u, bcs)

uOut = File("displacement.pvd")
uOut << u




# postprocessing: plot normal stress on top surface
###################################################

# normal vector in z direction
nz = Constant((0, 0, 1))

# get solution for normal stress component
surface_traction_top_normal = project(dot(dot(sigma(u),nz),nz))

# get coordinates spanning top face
# (x, y, z=thickness)
x_sweep = np.linspace(0, length, 10 * nx)
y_sweep = np.linspace(0, width, 10 * ny)

XY_X, XY_Y = np.meshgrid(x_sweep, y_sweep)
XY_X = XY_X.flatten()
XY_Y = XY_Y.flatten()
XY_Z_hi = thickness * np.ones(XY_Y.size)


# store interpolated solution at discrete points
surface_traction_top_normal_discrete = np.zeros(XY_Y.size)

for i, (x,y,z) in enumerate(zip(XY_X, XY_Y, XY_Z_hi)):
    surface_traction_top_normal_discrete[i] = surface_traction_top_normal(Point(x,y,z))

# find magnitude of stress data, so that we can place 0 in the middle of the plotting scale
magnitude = max(abs(surface_traction_top_normal_discrete))

# unflatten stress data so it can be plotted coherently
surface_traction_top_normal_discrete = surface_traction_top_normal_discrete.reshape((len(y_sweep), len(x_sweep)))

# plot
plt.figure()
ax = plt.axes()
ax.set_aspect('equal')

plt.pcolor(x_sweep, y_sweep, surface_traction_top_normal_discrete, cmap='RdBu_r', vmin=-magnitude, vmax=magnitude)
plt.colorbar(label='$t_{top} \\cdot e_3$')
plt.xlim([min(x_sweep), max(x_sweep)])
plt.savefig('top_stress_e3.png')

Answer 1

You could project to a DG1 field

 nz = Constant((0, 0, 1))
 surface_traction_top = project(dot(sigma(u),nz), VectorFunctionSpace(mesh,'DG',1))
 surface_traction_top_normal = project(dot(dot(sigma(u),nz),nz) FunctionSpace(mesh,'DG',1))

Comments

When I tried this, it still gave me the same problem:

I edited my question to include my full code, if that helps. I applied your suggestion around line 119.

As you can see in the figure, it's a rather coarse mesh. But I found that changing the number of elements doesn't help that much. In particular, increasing the number of elements in the x-direction merely puts a higher surface stress in a smaller area.

---

I see that you are using a CG2 field, so I guess you should be using a DG1 stress representation. I have edited my answer accordingly.

---

Thank you, this mostly solves the problem. If you have time, I have one more quick question about numerical methods.

I started computing the surface integral of normal stress to get a quick statistic. It should equal zero if the Neumann BCs are being obeyed.

surf_integral_top_t3 = assemble(surface_traction_top_normal * ds(2))

Before applying your advice, this had a value of about 1e-3.

With your advice, it becomes about 3e-4, which is better, though by no means perfect.

I tried also refining the mesh near that interface, which led to a progressive decrease in the value. As I refined the mesh more, it reached negative 2e-4.

For comparison, I tried the same thing, but projected onto CG1 instead of DG1. For low mesh densities, it didn't work as well as DG1, but as I increased the density, it converged to -5e-5, which is closer to the right answer.

Why did CG1 perform better? I feel as though DG1 should be better able to represent that sharp edge where the two boundaries meet, but that was only for low mesh densities. Is DG1 less stable in general?