Compression of a box

[Salim]

Hi,
I am a beginner in fenics, i would like to simulate the compression of a box, i did the simulation but resulats (using paraview) was illogical especially for the Von Mises. Can you find the mistakes made in my code:

from fenics import *
from ufl import nabla_div
rho = 1500e-9
E = 3000.2
nu = 0.38
mu = E/(2*(1+nu))
lambda_ = (nu*E)/((1+nu)*(1-2*nu))
tol=1E-9
mesh = BoxMesh(Point(0,0,0),Point(100,100,100),20,20,20)
class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]<0+tol
class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]>100-tol
sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryBottom.mark(sub_domains, 1)
boundaryTop = BoundaryTop()
boundaryTop.mark(sub_domains, 2)
V = VectorFunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)
bc1= DirichletBC(V, Constant((0,0,0)), boundaryBottom)
bc2= DirichletBC(V, Constant((0,0,-30)), 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)*ds
u=Function(V)
solve(a==L, u, bc)
vtkfile = File('u.pvd')
vtkfile << u
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
W = FunctionSpace(mesh, 'CG', 0)
von_Mises = project(sqrt (3/2.* inner (s, s) ), W)
vtkfile = File('s.pvd')
vtkfile << von_Mises

Thanks in advance for your help!

[dokken]

Please format your code using 3x` such that it is possible to copy-paste the code to an editor and run it without modifications.
Please also note that your VonMises stresses should be projected into a DG 0 space, as the displacement lives in CG 1. This has also been noted in other questions such as:

Results obtained by FEniCS mismatches Abaqus

The code you supplied is not complete, and therefore there is not much one can do to investigate this. Note that the von Mises stresses should be projected into W (the discontinuous space).

[Salim]

I think that now the code can be copied, Thanks in advance for your help! But beside if I project the VonMises stresses into a DG 0 space, i have the same distributions and shape

[dokken]

As long as you save it to a PVD, the data will be stored at vertices, and not be properly discontinuous, you should use XDMFile.write_checkpoint, as shown in Loading xdmf data back in - #4 by dokken

[bleyerj]

Hi,
if you restrain the lateral displacements on the top and bottom surfaces, you are not simulating simple compression. This will result in high stress concentrations near the border/corners of you domain as illustrated in your figure.
Consider something like:

bc1= DirichletBC(V.sub(2), Constant(0.), boundaryBottom)
bc2= DirichletBC(V.sub(2), Constant(-30), boundaryTop)

You also need to remove rigid body motions by fixing horizontal displacement and rotations about the Z axis : either by enforcing u_x=0 (resp. u_y=0) on a symmetry plane such as x=0 (resp. y=0) or by fixing appropriate displacements at some specific points if you are not satisfied with symmetry conditions

[Salim]

I saved the Von mises in XDMF file using:

save = XDMFFile('u_xdmf.xdmf')
save.write_checkpoint(von_Mises, "mises", 0)

But I get the same result as before, the distribution of von mises if different than the resulat that i get using comsol for example

[dokken]

see @kamensky’s comment above

[Salim]

Thank you very much for the help. Now, the von mises distribution is logical.

but there is a difference between the comsol results and fenics results, so i decided to use the same mesh in both, I exported the mesh created in comsol in stl format then i converted into xml format using meshio to use it in fenics. But, using the new mesh leads to illogical results.

Have you an idea about the cause and how can I resolve it?
Thanks in advance for your help!

[Ekrem_Ekici]

It would be nice if you rescale the colormaps for sake of comparison.

[Salim]

It’s done, I changed the color and data range

[dokken]

Please add the code that you are currently executing as to reflect the changes that you have made since the initial post.

[Salim]

from fenics import *
from ufl import nabla_div
E = 2500
nu = 0.4
mu = E/(2*(1+nu))
lambda_ = (nu*E)/((1+nu)*(1-2*nu))
tol=1E-6
mesh = Mesh('Box.xml')
class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]<0+tol
class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]>100-tol
sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryBottom.mark(sub_domains, 1)
boundaryTop = BoundaryTop()
boundaryTop.mark(sub_domains, 2)
V = VectorFunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)
bc1= DirichletBC(V, Constant((0,0,0)), boundaryBottom)
bc2= DirichletBC(V.sub(2), Constant(-30), 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)*ds
u=Function(V)
solve(a==L, u, bc)
vtkfile = File('u.pvd')
vtkfile << u
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
W = FunctionSpace(mesh, 'P', 1)
von_Mises = project(sqrt (3/2.* inner (s, s) ), W)
vtkfile = File('s.pvd')
vtkfile << von_Mises
save = XDMFFile('Von_mises.xdmf') 
save.write_checkpoint(von_Mises, "mises", 0) 

I has the same illogical distribution when I I project the VonMises stresses into a DG 0 space.
In this code, as I said before, i converted the box from STL format using meshio into xml format to use it in fenics. When I used the built-in mesh i has logical results, but it’s not the case with an imported mesh

[dokken]

You get quite different results if you use CG 1 or DG 0 visualization. Here is your code with a built in mesh, using DG-0:

**strong text**from fenics import *
from ufl import nabla_div
E = 2500
nu = 0.4
mu = E/(2*(1+nu))
lambda_ = (nu*E)/((1+nu)*(1-2*nu))
tol=1E-6
N = 20
mesh = BoxMesh(Point(0,0,0),Point(100,100,100), N, N, N)
class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]<0+tol
class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]>100-tol
sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryBottom.mark(sub_domains, 1)
boundaryTop = BoundaryTop()
boundaryTop.mark(sub_domains, 2)
V = VectorFunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)
bc1= DirichletBC(V, Constant((0,0,0)), boundaryBottom)
bc2= DirichletBC(V.sub(2), Constant(-30), 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)*ds
u=Function(V)
solve(a==L, u, bc)
vtkfile = File('u.pvd')
vtkfile << u
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
W = FunctionSpace(mesh, 'DG', 0)
von_Mises = project(sqrt (3/2.* inner (s, s) ), W)
save = XDMFFile('Von_mises.xdmf') 
save.write_checkpoint(von_Mises, "mises", 0) 

and you get:

as compared to

which I claim is radically different solutions, as the magnitude of the stress is different
DG 0 ranges from (630.882 to 1477.44) on the surface
CG 1 ranges from (633.903 to 1403.21) on the surface

@Salim

When I used the built-in mesh i has logical results, but it’s not the case with an imported mesh

You need to supply your mesh if you want any further guidance on this subject. Could you please provide your original mesh file?

[Salim]

Mesh (Stl format)
Mesh (Xml format after the conversion)

[dokken]

As you can see from the xml data in your file, your mesh in only a surface mesh (it does not contain tetrahedron, only triangles). This is because your step file only contains the triangular elements.