Plane Stress Analysis Comparison between FEniCSx and COMSOL

dolfinx
1

Hello, I am using the following tutorial (Isotropic and orthotropic plane stress elasticity — Computational Mechanics Numerical Tours with FEniCSx) to solve a simple Plane Stress Problem. The left edge is restrained and the right edge is having a force.

But when I solve the same problem in COMSOL, I get difference in displacement solution.

What could be the reason for this deviation?

Here is my mesh file:

cl__1 = 1e+22;
Point(1) = {0, 0, 0, cl__1};
Point(2) = {2, 0, 0, cl__1};
Point(3) = {2, 0.5, 0, cl__1};
Point(4) = {0, 0.5, 0, cl__1};
Physical Line("1") = {1};
Physical Line("2") = {2};
Physical Line("3") = {3};
Physical Line("4") = {4};
Physical Surface("mysurface") = {1};

And here is my code:

import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import sys
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
import dolfinx.cpp as _cpp
from dolfinx import mesh, fem, io, common, default_scalar_type, la
import dolfinx.fem.petsc
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities, locate_entities
from ufl import (FacetNormal, FiniteElement, as_matrix, Identity, Measure, TestFunction, TrialFunction, VectorElement, as_vector, sqrt, div, tr, dot, ds, dx, inner, lhs, grad, nabla_grad, nabla_div, rhs, sym, as_tensor)

# Mesh
mesh_path = "BM111.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(mesh_path, MPI.COMM_WORLD, gdim=2)

import ufl

gdim = 2

E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
z = E/(1-nu**2)

def strain(u, repr ="vectorial"):
    eps_t = sym(grad(u))
    if repr =="vectorial":
        return as_vector([eps_t[0,0], eps_t[1,1], (2*eps_t[0,1])])
    elif repr =="tensorial":
        return eps_t

C = as_matrix([[z, z*nu, 0.0],[z*nu, z,  0.0],[0.0, 0.0, 0.5*z*(1-nu)]])

def stress(u, repr ="vectorial"):
    sigv = dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))

#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
uh = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

T = fem.Constant(domain, 100000000.0)

n = FacetNormal(domain)

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)

L_form = dot(T*n,u_) * ds(6)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()

left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(8))

ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)
bcs = [fem.dirichletbc(ux0, left_dofs, V.sub(0)), fem.dirichletbc(uy0, left_dofs, V.sub(1))]

WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=uh, bcs=bcs)
uh = WPsolve.solve()
uh.name = "u"

Here is the displacement plot from COMSOL:

BM1
BM1920×705 22.6 KB

and FEniCSx:

Screenshot from 2024-10-31 14-49-27
Screenshot from 2024-10-31 14-49-271301×520 10.4 KB

2

In your problem, you applied T = 1e8 Pa, so the x-directional stress should be nearly 1e8 Pa in the most of the domain.
Given that the elastic modulus is 2e11 Pa, the x-directional strain should be around 0.5e-3 in the most of the domain.
Therefore, the x-directional elongation should be roughly 1e-3 m, which seems to agree with your fenics result.
(Could anyone verify if my understanding of solid mechanics is correct?)

Please review your COMSOL setup, especially the natural boundary condition.
Did you applied pressure (traction) or total force?

1 Like
3

Thank you @jpark. I was making this mistake of applied loading. I was applying total force but now when I applied pressure, the result turned out to be correct. Thanks once again.