The verification of thermoelastic von Mises stress

dolfinx
1

Hello everyone, I have calculated the von Mises stress of a rectangular prism at 100°C (with x = 0 fixed), and validated the results using Ansys under the same conditions. The displacement results are quite accurate, but there is a significant difference in the von Mises stress results. Could anyone explain the reason for this? Here is my code：

from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

L = 1
W = 0.2
E = 2e5
nu = 0.3
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1.2e-5

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],[20, 6, 6], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))

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

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

def eps(v):
    return 0.5 * (ufl.grad(v) + ufl.grad(v).T)
def sigma(v, dT):
    return (lmbda * ufl.tr(eps(v)) - alpha * (3 * lmbda + 2 * mu) * dT) * ufl.Identity(3) + 2.0 * mu * eps(v)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
Wint = ufl.inner(sigma(u, 100), eps(v)) * ufl.dx
a = ufl.lhs(Wint)
L = ufl.rhs(Wint)
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
#Mises stress
s = sigma(uh, 100) - 1. / 3 * ufl.tr(sigma(uh, 100)) * ufl.Identity(3)
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))

V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)
#save results
with io.XDMFFile(domain.comm, "out_elasticity/u.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(uh)
with io.XDMFFile(domain.comm, "out_elasticity/von_mises_stress.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(stresses)

Followed by the successful displacement validation results:

0c0d53fe63a25e8110a3290afef823f
0c0d53fe63a25e8110a3290afef823f761×603 25.2 KB

and then the von Mises stress results:
1fea5b40ed2dd2a98341c67d72eb2e9
1fea5b40ed2dd2a98341c67d72eb2e9747×619 19.3 KB

2

In general, the results depend on the mesh and the order of the elements.

Additionally, as far as I know, ANSYS uses averaging to compute nodal stresses from the element stresses surrounding a node. In contrast, in your FEniCS results, the stress is constant across each element because you used a DG element of order 0. Clearly, the two results cannot be the same.

3

Thank you for your reply. I am using exactly the same mesh, and the reason for the discrepancy is likely due to the different calculation methods you mentioned later. Is there any way to obtain consistent results, similar to the displacement results?

4

Using P1 elements for displacement, any stressed does clearly not live in P1 (which is what Andys is using.

By refining your grid, you should get a result that looks more and more equivalent.

5

Thank you, I modified the method for calculating the Mises stress in ANSYS and obtained results that are quite close.

6

Thank you. After modifying the method for calculating the Mises stress in ANSYS, I further refined the mesh and obtained more accurate results.