Issue with Von Mises stress calculation

General
1

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[1] < h+tol
    
    
class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return  x[1] > 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):
    return sym(grad(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 …

thanks in advance,

Quentin

Screenshot 2023-11-14 162231
Screenshot 2023-11-14 1622311469×1517 230 KB

2

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
3

Good morning Dokken,

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

Best regards,

Quentin

4

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

5

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
6

Thanks for the quick reply.

just need to get used to the philosophy of FEniCS.