Comparative Analysis: Ansys vs. FEniCSx [Displacement in Linear Elasticity]

[jaribeiro]

Hi everyone,

I’m currently conducting a validation study of my code against Ansys results. My project involves a 2D plane stress case in linear elasticity, and my code is based on this tutorial by @dokken.

However, I’m encountering some discrepancies in the displacement results:

• displacement in x:
  

  

• displacement in y:
  

  

I suspect that the differences may be due to the precision settings used in Ansys compared to those in FEniCSx. Does anyone have any suggestions on how to address this? Also, how can I modify the precision settings in FEniCSx?

I’ve tried altering the petsc_options, such as setting the absolute tolerance ("ksp_atol") and relative tolerance ("ksp_rtol"), but it doesn’t seem to have any effect. Any advice or insights would be greatly appreciated.

Thanks!

[bleyerj]

Are you sure that you are indeed doing plane stress computations in FEniCS ? The tutorial by @dokken considers the general linear isotropic constitutive equation for 3D elasticity. IN 2D this reduces to plane strain conditions. For plane stress, you need to change the constitutive equation very slightly, see here:
https://bleyerj.github.io/comet-fenicsx/tours/linear_problems/isotropic_orthotropic_elasticity/isotropic_orthotropic_elasticity.html#equation-isotropic-plane-stress
I suggest you use a direct solver to avoid any issue with solver tolerances but based on the large differences I suggest your problem is more likely due to plane stress/plane strain conditions.

[jaribeiro]

Thank you, @bleyerj! I’ve also used your old tutorial for guidance, and I already have the stress and strain calculations for plane stress conditions sorted out, so that aspect seems to be fine.

I plotted the error squared for the displacement in both the x and y directions, and observed that each has an error squared of 1:

Therefore, I think that the problem needs to be something minor

[bleyerj]

Then I suggest, test against a much simpler example in terms of geometry and boundary conditions. Upon mesh refinement, do you see a reduction in the error ?

[dokken]

I would also strongly suggest that you share the code that you are running.
Even if @bleyerj is addressing very good points, it is hard to tell what is wrong without the code at hand.

[jaribeiro]

I’m going to run some tests, increasing the refinement and the element degree, and I’ll share the results here soon.

These new results should be more accurate and closely aligned with reality.

However, I’m intrigued as to why the two software programs yield different results under the same study conditions. Do you have any insight into why there might be a discrepancy?

[jaribeiro]

I have a modular code, which is complex and time-consuming for you to assist with. I will create a simpler version and post it here soon, thanks @dokken

[jaribeiro]

@dokken and @bleyerj, here’s my code, simplified as much as possible (you’ll need to import the mesh). Thank you for the help!

The boundary conditions are:

You can find the mesh at this link
(The link also includes a Python script to facilitate the process and the results)

code:

import dolfinx
import json
from mpi4py import MPI
import numpy as np
from petsc4py.PETSc import ScalarType
import shutil
import time
import ufl


# Material
def eps(v):
    return ufl.sym(ufl.grad(v))


E = 200000.0
nu = 0.3
model = "plane_stress"

shear_modulus = E / 2 / (1 + nu)
lame_modulus = E * nu / (1 + nu) / (1 - 2 * nu)


def get_strain_tensor(displacement):
    strain = 0.5 * (ufl.nabla_grad(displacement) +
                    ufl.nabla_grad(displacement).T)
    strain_zz = -lame_modulus / (lame_modulus + 2 * shear_modulus) * (
        strain[0, 0] + strain[1, 1])

    strain = ufl.as_tensor([[strain[0, 0], strain[0, 1], 0],
                            [strain[1, 0], strain[1, 1], 0], [0, 0, strain_zz]])
    return strain


def get_stress_tensor(displacement):
    lame_modulus_2 = 2 * lame_modulus * shear_modulus / (lame_modulus +
                                                       2 * shear_modulus)

    dim = len(displacement)
    strain = 0.5 * (ufl.nabla_grad(displacement) +
                    ufl.nabla_grad(displacement).T)
    stress = lame_modulus_2 * ufl.nabla_div(displacement) * ufl.Identity(
        dim) + 2 * shear_modulus * strain

    stress = ufl.as_tensor([[stress[0, 0], stress[0, 1], 0],
                            [stress[1, 0], stress[1, 1], 0], [0, 0, 0]])
    return stress


def get_von_mises_stress(displacement):
    d = len(displacement) + 1
    deviatoric_stress = get_stress_tensor(displacement) - (1. / 3) * ufl.tr(
        get_stress_tensor(displacement)) * ufl.Identity(d)

    von_mises_stress = ufl.sqrt(3. / 2 *
                                ufl.inner(deviatoric_stress, deviatoric_stress))

    return von_mises_stress


# Mesh
mesh_path = "mesh.msh"
dolfinx_mesh, _, _ = dolfinx.io.gmshio.read_from_msh(mesh_path,
                                                     MPI.COMM_WORLD,
                                                     gdim=2)


# Boundary
def bottom(x):
    return np.isclose(x[1], min(x[1]))


def top(x):
    return np.isclose(x[1], max(x[1]))


def left(x):
    return np.isclose(x[0], min(x[0]))


def right(x):
    return np.isclose(x[0], max(x[0]))


# Solve
v_func_space = dolfinx.fem.VectorFunctionSpace(dolfinx_mesh, ("S", 1))
u_trial_func = ufl.TrialFunction(v_func_space)
v_test_func = ufl.TestFunction(v_func_space)

# Bilinear form
stress_u = get_stress_tensor(u_trial_func)
strain_v = get_strain_tensor(v_test_func)
bilinear_form_a = ufl.inner(stress_u, strain_v) * ufl.dx

# Lienar form
facet_indices = []
facet_markers = []
facet_dim = dolfinx_mesh.topology.dim - 1

facet_indices_boundary_left = dolfinx.mesh.locate_entities_boundary(
    dolfinx_mesh, facet_dim, left)
facet_indices.append(facet_indices_boundary_left)
facet_markers.append(np.full_like(facet_indices_boundary_left, 0))

facet_indices_boundary_top = dolfinx.mesh.locate_entities_boundary(
    dolfinx_mesh, facet_dim, top)
facet_indices.append(facet_indices_boundary_top)
facet_markers.append(np.full_like(facet_indices_boundary_top, 1))

facet_indices_boundary_right = dolfinx.mesh.locate_entities_boundary(
    dolfinx_mesh, facet_dim, right)
facet_indices.append(facet_indices_boundary_right)
facet_markers.append(np.full_like(facet_indices_boundary_right, 2))

facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)

facet_tag = dolfinx.mesh.meshtags(dolfinx_mesh, facet_dim,
                                  facet_indices[sorted_facets],
                                  facet_markers[sorted_facets])

ds = ufl.Measure("ds", domain=dolfinx_mesh, subdomain_data=facet_tag)

tension_left = dolfinx.fem.Constant(dolfinx_mesh, ScalarType((1, 0)))
tension_top = dolfinx.fem.Constant(dolfinx_mesh, ScalarType((0, 1)))
tension_right = dolfinx.fem.Constant(dolfinx_mesh, ScalarType((-1, 0)))
tensions = [tension_left, tension_top, tension_right]

linear_form_l = 0
for i, tension in enumerate(tensions):
    linear_form_l += ufl.dot(tension, v_test_func) * ds(i)

# BCs
sub_0, _ = v_func_space.sub(0).collapse()
edges = dolfinx.mesh.locate_entities_boundary(dolfinx_mesh,
                                              dolfinx_mesh.topology.dim - 1,
                                              bottom)
boundary_dofs = dolfinx.fem.locate_dofs_topological(
    (v_func_space.sub(0), sub_0), dolfinx_mesh.topology.dim - 1, edges)
disp_func = dolfinx.fem.Function(sub_0)
x_disp_const = dolfinx.fem.Constant(dolfinx_mesh, ScalarType(0.0))
x_disp_expr = dolfinx.fem.Expression(x_disp_const,
                                     sub_0.element.interpolation_points())
disp_func.interpolate(x_disp_expr)
bc_x = dolfinx.fem.dirichletbc(disp_func, boundary_dofs, v_func_space.sub(0))

sub_1, _ = v_func_space.sub(1).collapse()
edges = dolfinx.mesh.locate_entities_boundary(dolfinx_mesh,
                                              dolfinx_mesh.topology.dim - 1,
                                              bottom)
boundary_dofs = dolfinx.fem.locate_dofs_topological(
    (v_func_space.sub(1), sub_1), dolfinx_mesh.topology.dim - 1, edges)
disp_func = dolfinx.fem.Function(sub_1)
y_disp_const = dolfinx.fem.Constant(dolfinx_mesh, ScalarType(0.0))
y_disp_expr = dolfinx.fem.Expression(y_disp_const,
                                     sub_1.element.interpolation_points())
disp_func.interpolate(y_disp_expr)
bc_y = dolfinx.fem.dirichletbc(disp_func, boundary_dofs, v_func_space.sub(1))

bc = [bc_x, bc_y]

linear_variational_problem = dolfinx.fem.petsc.LinearProblem(
    a=bilinear_form_a,
    L=linear_form_l,
    bcs=bc,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu"
    })
displacement_func = linear_variational_problem.solve()

# von mises
v_func_space = dolfinx.fem.FunctionSpace(dolfinx_mesh, ("S", 1))
tensor = get_von_mises_stress(displacement_func)
field_expr = dolfinx.fem.Expression(tensor,
                                    v_func_space.element.interpolation_points())
von_mises_func = dolfinx.fem.Function(v_func_space)
von_mises_func.interpolate(field_expr)

# Save to xdmf
fields = {displacement_func, von_mises_func}

with dolfinx.io.XDMFFile(dolfinx_mesh.comm, "test.xdmf", "w") as xdmf:
    xdmf.write_mesh(dolfinx_mesh)
    displacement_func.name = "displacement"
    xdmf.write_function(displacement_func)
    von_mises_func.name="von_mises"
    xdmf.write_function(von_mises_func)

[dokken]

@jaribeiro

v_func_space = dolfinx.fem.FunctionSpace(dolfinx_mesh, ("S", 1))
tensor = get_von_mises_stress(displacement_func)
field_expr = dolfinx.fem.Expression(tensor,
                                    v_func_space.element.interpolation_points())
von_mises_func = dolfinx.fem.Function(v_func_space)
von_mises_func.interpolate(field_expr)

You should not put Von-Mises stresses in “S”. They should be in DG-0. See for instance: Projection and interpolation — FEniCS Tutorial @ Sorbonne
and the text above.

[jaribeiro]

Thanks, @dokken! I attempted your approach, but it appears that the von Mises stress results differ significantly from the Ansys ones. I’m also not quite clear on why I need to use DG-0. Could you provide more detailed information or resources to help me understand better? I would really appreciate it.

Also, I’m wondering about the varying displacements. Could these variations be influencing the von Mises stress calculations?

[driley]

Hi @jaribeiro, from a quick look it seems that ANSYS uses an alternative definition of their Von Mises stress. Notably, they show at least on their website the use of the total stress tensor.

As for DG0, it has to do with the fact that to compute strain and therefore stress you take the derivative of displacement so your stress/strain states are now piecewise discontinuous.

[jaribeiro]

Thank you, @driley!

Your point is insightful. It seems they first calculate the stress and then the von Mises stress, similar to what’s discussed in this article: ANSYS blog on equivalent stress.

From my understanding, in relation to DG-0, to compute the stress tensor, the calculation would be as follows, correct?

v_func_space = dolfinx.fem.TensorFunctionSpace(
            dolfinx_mesh,
            ("DG", 0),
            shape=(3, 3))
tensor = get_stress_tensor(displacement_func)
field_expr = dolfinx.fem.Expression(tensor,
                                    v_func_space.element.interpolation_points())
stress_func = dolfinx.fem.Function(v_func_space)
stress_func.interpolate(field_expr)

However, I’ve noticed some discrepancies in the stress results using FEniCSx.

Stress XX:

Stress XY:

Stress YY:

Do you have any insight on why this might be occurring? Is it normal to see minor differences in the displacements?

[dokken]

Do you know what quadrature rules ANSYS use to integrate contributions for each non-affine quadrilateral?

As the mapping is non-affine, the jacobians are not constant per element, and a suitable quadrature rule has to be chosen.

Also, looking at your stress plots, it does not seem to be visualized as cellwise cosntants.
This is how your code visualizes on my system:

when using DG-0 for the von MIses stresses.

[jaribeiro]

This is an insightful observation, thanks @dokken. At the moment, I’ve been assuming they are affine quadrilaterals, which could affect the displacement results differently.

ANSYS utilizes Gauss quadrature. Do you have any idea how this could be implemented in Fenicsx? I’ve never attempted anything similar before.

Also, regarding the stress tensor, can I use the code provided above? Or will introducing the quadrature require changes to it? Thanks

v_func_space = dolfinx.fem.TensorFunctionSpace(
            dolfinx_mesh,
            ("DG", 0),
            shape=(3, 3))
tensor = get_stress_tensor(displacement_func)
field_expr = dolfinx.fem.Expression(tensor,
                                    v_func_space.element.interpolation_points())
stress_func = dolfinx.fem.Function(v_func_space)
stress_func.interpolate(field_expr)

Referring to the plot, it seems I initially opened the incorrect one, but now I have the correct one, thanks for catching the error:

[dokken]

@jaribeiro

This is an insightful observation, thanks @dokken. At the moment, I’ve been assuming they are affine quadrilaterals, which could affect the displacement results differently.

The quadrilaterals in your grid does not look affine (not all of them are parallelograms as far as I can tell), which changes the quadrature degree required.

You can use any of the quadrature rules defined in basix, or use your own,
See:How to choose a specific quadrature rule? - #6 by dokken
or Total number of gauss point - #2 by dokken

[jaribeiro]

Thanks, @dokken! I followed your links and it seems that I only need to change where I have the dx and ds, and add the quadrature properties there, right?

bilinear_form_a = ufl.inner(stress_u, strain_v) * ufl.dx(metadata={"quadrature_degree": quad_degree, "quadrature_rule":quad_scheme})

linear_form_l += ufl.dot(tension, v_test_func) * ds(i, metadata={"quadrature_degree": quad_degree, "quadrature_rule": quad_scheme})

And, do I also need to use DG-0 for the stress tensor, similar to how I used it for the von Mises stress?

v_func_space = dolfinx.fem.TensorFunctionSpace(
            dolfinx_mesh,
            ("DG", 0),
            shape=(3, 3))
tensor = get_stress_tensor(displacement_func)
field_expr = dolfinx.fem.Expression(tensor,
                                    v_func_space.element.interpolation_points())
stress_func = dolfinx.fem.Function(v_func_space)
stress_func.interpolate(field_expr)

[dokken]

@jaribeiro

Thanks, @dokken! I followed your links and it seems that I only need to change where I have the dx and ds, and add the quadrature properties there, right?

bilinear_form_a = ufl.inner(stress_u, strain_v) * ufl.dx(metadata={"quadrature_degree": quad_degree, "quadrature_rule":quad_scheme})

linear_form_l += ufl.dot(tension, v_test_func) * ds(i, metadata={"quadrature_degree": quad_degree, "quadrature_rule": quad_scheme})

Yes

@jaribeiro

And, do I also need to use DG-0 for the stress tensor, similar to how I used it for the von Mises stress?

Yes.
The stress tensor contains the derivative of a piecewise linear function, whose derivative in any direction is a element wise constant function

[jaribeiro]

Many thanks to @dokken, I am beginning to make significant progress in resolving my problem.

@dokken

The stress tensor contains the derivative of a piecewise linear function, whose derivative in any direction is a element wise constant function

Despite utilizing DG-0, I observe that each stress component is assigned per node rather than per cell, which is different from the approach used in von Mises stress. Is this a reasonable outcome? Could this be attributed to the use of TensorFunctionSpace instead of the FunctionSpace (von Mises stress case)?

Considering this, can I find a way to get the von Mises stress values at the nodes?

[dokken]

@jaribeiro

Despite utilizing DG-0, I observe that each stress component is assigned per node rather than per cell, which is different from the approach used in von Mises stress. Is this a reasonable outcome? Could this be attributed to the use of TensorFunctionSpace instead of the FunctionSpace (von Mises stress case)?

This is probably due to limitations in the XDMF-file format.
What i would suggest to do is to use VTXWriter, and instead of interpolating into DG 0 interpolate into a DG-1 space (as this is supported by VTX, and will be the same as DG 0 if the stress is DG-0).

[dokken]

@jaribeiro

Considering this, can I find a way to get the von Mises stress values at the nodes?

The values are not well defined at nodes when using a first order serendipity (equivalent to Lagrange) element, as the derivatives are discontinuous at both facets and vertices between cells.

You should in such a case create a special operator to average this if you want to.

Note that projections introduces Gibbs phenomenon, as shown in

@dokken

You should not put Von-Mises stresses in “S”. They should be in DG-0. See for instance: Projection and interpolation — FEniCS Tutorial @ Sorbonne
and the text above.