# Issue with Von Mises stress calculation

Hi all,

I am facing a weird issue while trying to calculate the Von Mises stress from a displacement field.
I obtain a negative result, which should not be as this Von mises stress is the defined as a square root of some kind of tensor L2 norm.

I am using the following code:

``````from fenics import *
from dolfin import *
from ufl import nabla_div

E = 62000
nu = 0.24
Sig_y = 7.0
mu = E/(2*(1+nu))
lambda_ = nu*E/((1+nu)*(1-2*nu))
Utop = -0.1

N = 16
h = 1/N
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), N, N, N)
tol=1E-8

class BoundaryBottom(SubDomain):
def inside(self, x, on_boundary):
return x < h+tol

class BoundaryTop(SubDomain):
def inside(self, x, on_boundary):
return  x > 1.0 - (h+tol)

sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryTop = BoundaryTop()
boundaryBottom.mark(sub_domains, 1)
boundaryTop.mark(sub_domains, 2)

# mesh
vtk_file = File("VTK/Subdomains.pvd")
vtk_file << sub_domains

# Define Functional spaces V for variational formulation and projecting quantities
# sequential version
V = VectorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "CG", 1)
Ws = TensorFunctionSpace(mesh, "CG", 1)
# functions and BC)
u = TrialFunction(V)
v = TestFunction(V)
bc1= DirichletBC(V, Constant((0,0,0)), boundaryBottom)
bc2= DirichletBC(V.sub(1), Utop, boundaryTop)
bc=[bc1,bc2]
d=u.geometric_dimension()

def epsilon(u):

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

a = inner(sigma(u), epsilon(v))*dx
f= Constant((0,0,0))
L = dot(f,v)*dx
u=Function(V)
solve(a==L, u, bc, solver_parameters={"linear_solver": "petsc","preconditioner": "ilu"})

# Post treatment computations
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
von_Mises = project(sqrt (3/2.* inner (s, s) ), W, solver_type="cg")

vtkfile = File('VTK/u.pvd')
vtkfile << u

vtkfile = File('VTK/VM.pvd')
vtkfile << von_Mises
``````

Any idea of where the problem could come from ?

I use a cg type projection as other type quickly lead to memory issue …

Quentin

You are projecting a first order continuous space, and is therefore observing Gibbs phenomenon,
As the stresses are cellwise discontinuous.
See for instance

http://jsdokken.com/FEniCS23-tutorial/src/approximations.html
that describes this in detail.

2 Likes

Good morning Dokken,

thanks for the link an projection on a 0 order DG space solves the problem.

Best regards,

Quentin

Just a detail that came to my mind:

In the above example I had on purpose used a quadratic functional space for the displacement in order to obtain a stress which is linear per element.

I would have naturally rely on finite element shape function to calculate displacement gradient, then strain and stress to compute the Von mises.

Is there an easy way to do that using FEniCS ?

Thanks,

Quentin

FEniCS and FEniCSx does not rely on isoparametric elements, i.e. you can use a higher order function space on a linear mesh.

Selecting an appropriate space for the stresses becomes slightly more complicated, as the product of two linear functions (the gradient) becomes a quadratic function, and the square root of a quadratic function is not in a function space.
Choosing a higher order DG space (degree>=2) for representing the stresses are recommended

1 Like