I am trying to simulate a linear elastic problem with a small amount of contact. To get the best results I was hoping to implement a Lagrange multiplier approach to solve for the displacements with constraints on the contact gap. When I read literature I typically find it is already in the matrix form and I am unsure of the best way to convert it into the variational form to solve in fenicsx.
Given a displacement u, a contact pressure K, and the normal to each node on the contact surface n the signorini contact conditions are as follows:
which also correspond to the KKT conditions where P is the lagrange multiplier and u \cdot n - g \leq 0 is the inequality constraint. The energy minimization is linear elastic and minimizes
with u = u_{applied} on \Gamma_{u} with \sigma being the stress and \epsilon being the strain. I tried to follow the tutorial for constrained minimization of the poisson equation, but since the gap function is not a PDE it doesnt seem to form a linear problem when solving the adjoint equation, therefore I tried to derive the variational form as follows:
with G(u) being the gap function seen as
and <x>_{+} being x when x is positive and 0 when x is negative, automatically satisfying the inequality constraint. To satisfy the equality constraint that G(u) = 0, the previous and following equations are solved using a block solver for both u and k being a CG 2d vector and CG scalar values respectively.
The minimum working example code is included below written on dolfinx 0.10.0 installed on docker.
from mpi4py import MPI
import gmsh
import numpy as np
import dolfinx.fem.petsc
import ufl
from dolfinx import mesh
from BoltGroupComputer import BoltGroup
import pyvista
plotstress = True
# physical variables
N = 25
L = 1.0
D = 2.0*L
# modifying the applied load changes the max error (tested at finite difference step of 1e-6)
u_disp = 0.1
# material properties
E_val = 200.0
nu_val = 0.3
iniK = 50.0
# values for loading and contact
R = .5
d = 0.0
# create mesh
comm = MPI.COMM_WORLD
domain = mesh.create_rectangle(
comm,
[np.array([-L, -D]), np.array([L, 0])],
[N, int(round(N*D/L))],
cell_type=mesh.CellType.quadrilateral,
)
# locate all exterior faces (for later boolean operations)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
gdim = domain.geometry.dim
boundary_facets_all = dolfinx.mesh.exterior_facet_indices(domain.topology)
# define the mesh marker for each facets on the mesh
num_facets = domain.topology.index_map(fdim).size_local
markers = np.zeros(num_facets, dtype=np.int32)
top_mark,bottom_mark,free_mark = 1,2,3
def topBnds(x):
return np.isclose(x[1],0)
def bottomBnds(x):
return np.isclose(x[1],-D)
# locate different regions of facets where the boundary conditions will be applied
top_facets, markers[top_facets] = mesh.locate_entities_boundary(domain, fdim, topBnds), top_mark
bottom_facets, markers[bottom_facets] = mesh.locate_entities_boundary(domain, fdim, bottomBnds), bottom_mark
free_facets, markers[free_facets] = np.setdiff1d(boundary_facets_all,np.union1d(top_facets,bottom_facets)), free_mark
facet_tags = mesh.meshtags(domain, fdim, np.arange(num_facets, dtype=np.int32), markers)
subdomain_tags = np.hstack(top_facets)
gamma, gamma_to_domain = dolfinx.mesh.create_submesh(domain, fdim, subdomain_tags)[0:2]
V = dolfinx.fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
Q = dolfinx.fem.functionspace(gamma, ("Lagrange", 1))
W = ufl.MixedFunctionSpace(V, Q)
dx = ufl.Measure("dx", domain=domain)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
E = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(E_val))
nu = dolfinx.fem.Constant(domain, nu_val)
x = ufl.SpatialCoordinate(domain)
f = dolfinx.fem.Constant(domain, np.zeros(gdim, dtype=dolfinx.default_scalar_type))
v, w = ufl.TestFunctions(W)
u = dolfinx.fem.Function(V, name="Displacement")
k = dolfinx.fem.Function(Q, name="lambda")
k.x.array[:] = iniK
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
def epsilon(w):
return ufl.sym(ufl.grad(w))
def sigma(w, mu, lmbda):
ew = epsilon(w)
gdim = ew.ufl_shape[0]
return 2.0 * mu * epsilon(w) + lmbda * ufl.tr(ufl.grad(w)) * ufl.Identity(gdim)
F = ufl.inner(sigma(u, mu, lmbda), epsilon(v)) * dx
F -= ufl.inner(f, v) * dx # body force
# define the obstacle norm and gap functions
centroid = ufl.as_vector([0.0,R-d])
x = ufl.SpatialCoordinate(domain)
dist = ufl.sqrt(ufl.dot(centroid-x,centroid-x))
gap = dist-R
norm = (centroid-x)/dist
deformedGap = ufl.dot(norm,u)-gap
# derived equations
G = ufl.conditional(ufl.ge(deformedGap,0),deformedGap,0)
dGdu = ufl.derivative(G,u,v)
F += k*dGdu*ds(top_mark)
F += ufl.inner(G,w)*ds(top_mark)
residual = ufl.extract_blocks(F)
du, dpsi = ufl.TrialFunctions(W)
jac = ufl.derivative(F, u, du) + ufl.derivative(F, k, dpsi)
J = ufl.extract_blocks(jac)
bc_bottom = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, dolfinx.default_scalar_type((u_disp))),
dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, facet_tags.find(bottom_mark)), V.sub(1))
bc_fixed = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, dolfinx.default_scalar_type((0))),
dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, facet_tags.find(bottom_mark)), V.sub(0))
bcs = [bc_bottom, bc_fixed]
V_DG = dolfinx.fem.functionspace(domain, ("DG", 1, (domain.geometry.dim,)))
stress_space, stress_to_disp = V_DG.sub(0).collapse()
von_mises = dolfinx.fem.Function(stress_space, name="von_Mises")
u_dg = dolfinx.fem.Function(V_DG, name="u")
s = sigma(u, mu, lmbda) - 1.0 / 3 * ufl.tr(sigma(u, mu, lmbda)) * ufl.Identity(len(u))
von_Mises = ufl.sqrt(3.0 / 2 * ufl.inner(s, s))
stress_expr = dolfinx.fem.Expression(von_Mises, stress_space.element.interpolation_points)
entity_maps = [gamma_to_domain]
solver = dolfinx.fem.petsc.NonlinearProblem(
residual,
u=[u, k],
J=J,
bcs=bcs,
entity_maps=entity_maps,
petsc_options={
"snes_monitor": None,
"ksp_type": "preonly",
"pc_type": "hypre",
"ksp_error_if_not_converged": True,
# "snes_error_if_not_converged": True,
"snes_max_it": 1000,
},
petsc_options_prefix="signorini_",
)
solver.solve()
von_mises.interpolate(stress_expr)
u_dg.interpolate(u)
print("max deformation: " + str(np.max(np.abs(u.x.array))))
print("avg pressure value: " + str(np.mean(k.x.array)) + " and std: " + str(np.std(k.x.array)))
if (plotstress):
grid = dolfinx.plot.vtk_mesh(u_dg.function_space)
pyvista_grid = pyvista.UnstructuredGrid(*grid)
values = u_dg.x.array.reshape(-1, domain.geometry.dim)
values_padded = np.zeros((values.shape[0], 3))
values_padded[:, : domain.geometry.dim] = values
pyvista_grid.point_data["u"] = values_padded
warped = pyvista_grid.warp_by_vector("u")
stresses = np.zeros_like(u_dg.x.array)
stresses[stress_to_disp] = von_mises.x.array
stresses = stresses.reshape(-1, domain.geometry.dim)[:, 0]
warped.point_data["von_Mises"] = stresses
plotter = pyvista.Plotter()
plotter.add_mesh(pyvista_grid, style="wireframe", color="b")
plotter.add_mesh(warped, scalars="von_Mises", lighting=True, show_edges=True)
plotter.show_bounds()
plotter.view_xy()
plotter.show()
After running the simulation, the displacement is incorrect, but has a correct trend. One of the main concerns is during solving the value of K is not modified from the initial guess, and only u changes, making it behave like a worse performing penalty method but I’m not sure why since 2 equations are solved for. I think it has to do with how the inequality constraint is modified, but Im not sure how else to do it. In the problem, the top surface is the contact surface and has a circular indenter fixed above the domain. the bottom of the domain has a displacement applied of 0.1 upwards to make it contact the circular indenter (hertzian contact). The correct solution is included below.
Before trying the lagrange multiplier method I implemented a penalty approach and the latent variable approach from the fenicsx examples, and while they both work well at larger displacements, it seems like the smaller contact distances in the problem I am simulating cause convergence issues and have less accurate results than desired. This should be solved by the lagrange multiplier method since that directly solves for the contact pressure, but I just dont know how to formulate it, any help or tips would be greatly appreciated!
Thank you!
