Hi,
I made a simple setup with a rectangular box, one side has a zero displacement boundary condition within a circle, and the other side has a traction boundary condition that integrates to 1 N. The material parameters are E=1e4 and nu=0.2.
When I run the code, reducing the mesh size from 2 mm, dividing by 2^0.25, the displacement at the center of the traction circle, and the maximum principal stress seem to diverge instead of converging to a fixed value. A self-contained code with the loop computing these values, together with the mesh sizes, is shown below.
Why is this happening? I thought this only happens with a nodal loading, is it because the traction is a step function at the edge of the circle? Could this be fixed by having something like a cosine function going from 0 at the circle’s edge to 1 (or some other value to compensate so the total traction integrates to 1 N) at a smaller circle? How is this done in practice? The contact used in a corresponding setup was with a 1 mm thick and 4 mm diameter rubber pad. Is the zero displacement condition also causing this? How could that be fixed?
Thank you,
Alex
import numpy as np
import time
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, io, default_scalar_type
import dolfinx.fem.petsc
r_loading = 2.0e-3
def in_sphere(x, point, r_range):
return ((x[0]-point[0])**2 + (x[1]-point[1])**2 + (x[2]-point[2])**2) < r_range**2
def mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal):
dim = domain.topology.dim
fdim = dim - 1
domain.topology.create_connectivity(fdim, dim)
def left_facets(xq):
return in_sphere(xq, facet_node, r_loading)
def right_facets(xq):
return in_sphere(xq, grasp_node, r_loading)
facets_left = mesh.locate_entities(domain, fdim, left_facets)
facets_right = mesh.locate_entities(domain, fdim, right_facets)
facet_indices = np.concatenate([facets_left, facets_right])
facet_markers = np.concatenate([np.full_like(facets_left, 1, dtype=np.int32),
np.full_like(facets_right, 2, dtype=np.int32)])
# Sort by facet index
sort_perm = np.argsort(facet_indices)
facet_indices = facet_indices[sort_perm]
facet_markers = facet_markers[sort_perm]
mt = mesh.meshtags(domain, fdim, facet_indices, facet_markers)
return mt
def elasticity_parameters(E=10.0e4, nu=0.20):
"""Lamé parameters from E, nu."""
lambda_ = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
mu = E * 0.5 / (1.0 + nu)
return lambda_, mu
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u, lambda_, mu):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2.0 * mu * epsilon(u)
def build_primal_problem(domain, facet_tags, lambda_, mu, node_normal,
p_traction=1.0/(np.pi*r_loading**2), body_force=(0.0, 0.0, 0.0)):
dim = domain.topology.dim
fdim = dim - 1
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f_vec = fem.Constant(domain, default_scalar_type(body_force))
# Traction direction: push along +x (change sign if needed)
T = fem.Constant(domain, -p_traction * default_scalar_type((node_normal[0], node_normal[1], node_normal[2])))
a = ufl.inner(sigma(u, lambda_, mu), epsilon(v)) * ufl.dx
L = ufl.dot(f_vec, v) * ufl.dx + ufl.dot(T, v) * ds(2)
# Dirichlet BC on left jaw (tag 1)
left_facets = facet_tags.find(1)
domain.topology.create_connectivity(fdim, dim)
left_dofs = fem.locate_dofs_topological(V, fdim, left_facets)
u_D = np.array((0.0, 0.0, 0.0), dtype=default_scalar_type)
bc_left = fem.dirichletbc(u_D, left_dofs, V)
return V, a, L, [bc_left]
def principal_stress(domain, uh, lambda_, mu):
#Compute principal stresses via invariants (analytical eigenvalue solution)
sigma_uh = sigma(uh, lambda_, mu)
q = ufl.tr(sigma_uh) / 3 # mean stress (I1/3)
B = sigma_uh - q * ufl.Identity(3) # deviatoric stress tensor
j = ufl.tr(B * B) / 2 # J2 = 1/2 tr(B^2)
b = ufl.tr(B * B * B) / 3 # J3 = 1/3 tr(B^3)
p = 2 / ufl.sqrt(3) * ufl.sqrt(j + 1e-14) # p = 2/sqrt(3)*sqrt(J2) (small eps**2 added)
r = 4 * b / (p**3) # r = 4 J3 / p^3
r = ufl.max_value(ufl.min_value(r,+1-1e-7),-1+1e-7) # clamp r to [-1+eps, 1-eps] to avoid acos domain error
phi = ufl.acos(r) / 3
lambda2 = q + p * ufl.cos(phi) # largest principal stress (max)
V_dg = fem.functionspace(domain, ("DG", 0))
pmax_expr = fem.Expression(lambda2, V_dg.element.interpolation_points)
pmax = fem.Function(V_dg, name="PrincipalStress")
pmax.interpolate(pmax_expr)
return pmax, V_dg
def array_dist_squared(x, p):
return (x[:,0]-p[0])**2 + (x[:,1]-p[1])**2 + (x[:,2]-p[2])**2
def fea_solve(domain, facet_node, grasp_node, grasp_node_normal):
#lambda_, mu = elasticity_parameters2(domain, grasp_node, 4.0e-3)
lambda_, mu = elasticity_parameters()
facet_tags = mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal)
Vp, a_p, L_p, bcs_p = build_primal_problem(domain, facet_tags, lambda_, mu, grasp_node_normal)
problem = dolfinx.fem.petsc.LinearProblem(a_p, L_p, bcs=bcs_p,
petsc_options_prefix="linear_problem",
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
uh = problem.solve()
# Principal stress field
pmax, V_dg = principal_stress(domain, uh, lambda_, mu)
# MPS maximum value not considering low E area
dof_coords = V_dg.tabulate_dof_coordinates()
indices = np.where(np.logical_and(array_dist_squared(dof_coords, facet_node) > 5.0e-3**2, \
array_dist_squared(dof_coords, grasp_node) > 5.0e-3**2))
max_pmax = np.sort(pmax.x.array[indices])[-1]
# Get displacement at node where load is applied
dof_coords = Vp.tabulate_dof_coordinates()
distances = np.linalg.norm(dof_coords - grasp_node, axis=1)
dof_index = np.argmin(distances)
u_node = np.sqrt(uh.x.array[3*dof_index]**2 + uh.x.array[3*dof_index+1]**2 + uh.x.array[3*dof_index+2]**2)
return domain, uh, pmax, max_pmax, u_node
# Part FEA
##################################################################
length = 20.0e-3
width = 10.0e-3
grasp_facet_node = np.array([0.0, width/2, width/2])
grasp_node2 = np.array([length, width/2, width/2])
grasp_node_normals = np.array([1.0, 0.0, 0.0 ])
cell_size = 2.0e-3
num_loops = 9
cell_sizes = np.empty((num_loops,1))
u_disps = np.empty((num_loops,1))
stresses = np.empty((num_loops,1))
for n in range(num_loops):
n_l = int(length/cell_size)
n_w = int(width/cell_size )
domain = mesh.create_box(MPI.COMM_WORLD, ((0.0, 0.0, 0.0), (length,width,width)), (n_l,n_w,n_w), mesh.CellType.tetrahedron)
t_start = time.perf_counter()
domain_f, uh, pmax, max_stress, u_node = fea_solve(domain, grasp_facet_node, grasp_node2, grasp_node_normals)
t_end = time.perf_counter()
print(f"FEA elapsed time: {t_end - t_start:.6f} seconds")
cell_sizes[n] = cell_size
u_disps[n] = u_node
stresses[n] = max_stress
cell_size = cell_size/2.0**0.25
with io.XDMFFile(domain_f.comm, "max_p_stress_part.xdmf", "w") as xdmf:
xdmf.write_mesh(domain_f)
uh.name = "u"
xdmf.write_function(uh)
pmax.name = "PrincipalStress"
xdmf.write_function(pmax)
print('\n')
print(cell_sizes)
print('\n')
print(u_disps)
print('\n')
print(stresses)
