Penalty Contact with circular shaped boundary Condition

variational formulation
1

Hello everyone,

i am currently working on a contact model for a multibody contact of a rigid circular shaped indenter and an elastic domain. The indenter should be realized as a boundary condition because of a considerable difference in stiffness (steel compared to rubber).

This topic refers to:

However, I have difficulties creating a penalty term which use both x and y displacements of the rigid indenter as an additional energy term in the formulation. I tried lots of different ways to implement a circle equation but nothing worked. My idea was to use something like this:

def in_circle(x):
    return np.less_equal( ((x[0]-center_x)**2 + (x[1]-center_y)**2), ((D/2)**2))

The MWE should realize the sketch below, where a circular boundary should cause a static deformation realized with linear elasticity and the penalty approach:

image
image630×628 34.4 KB

The code for the dolfinx MWE + mesh generation using gmsh:

import gmsh
import numpy as np
from dolfinx import fem
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
from petsc4py.PETSc import ScalarType, Options
from ufl import dx, inner, derivative, tr, Identity, sym, grad, TrialFunction, TestFunction, Measure

gmsh.initialize()
gmsh.model.add("rectangle")
lc=1e-1

# generate points for inner contour
gmsh.model.geo.addPoint(1, 1, 0, lc, 1)
gmsh.model.geo.addPoint(-1, 1, 0, lc, 2)
gmsh.model.geo.addPoint(-1, -1, 0, lc, 3)
gmsh.model.geo.addPoint(1, -1, 0, lc, 4)
# generate inner contour
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
# generate points for outer contour
gmsh.model.geo.addPoint(2, 2, 0, lc, 5)
gmsh.model.geo.addPoint(-2, 2, 0, lc, 6)
gmsh.model.geo.addPoint(-2, -2, 0, lc, 7)
gmsh.model.geo.addPoint(2, -2, 0, lc, 8)
# generate outer contour
gmsh.model.geo.addLine(5, 6, 5)
gmsh.model.geo.addLine(6, 7, 6)
gmsh.model.geo.addLine(7, 8, 7)
gmsh.model.geo.addLine(8, 5, 8)
# generate curveloops
gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)
gmsh.model.geo.addCurveLoop([5, 6, 7, 8], 2)
# generate surface
gmsh.model.geo.addPlaneSurface([2, 1], 1)
gmsh.model.geo.synchronize()
# generate physical groups for boundary conditions
inside, outside, surface = 1, 2, 3
gmsh.model.addPhysicalGroup(1, [1, 2, 3, 4], inside, name="Inner Contour")
gmsh.model.addPhysicalGroup(1, [5, 6, 7, 8], outside, name="Outer Contour")
gmsh.model.addPhysicalGroup(2, [1], surface, name="Surface")
# finish gmsh
gmsh.model.mesh.generate(2)
gmsh.write("rectangle.msh") 
# create a DOLFINx mesh 
dim = 2
gmsh_model_rank = 0
domain_com = MPI.COMM_SELF
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, domain_com, gmsh_model_rank, dim)
domain.name = "Stator"
cell_markers.name = f"{domain.name}_cells"
facet_markers.name = f"{domain.name}_facets"

with XDMFFile(domain.comm, f"rectangle_mesh_{MPI.COMM_SELF.rank}.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_meshtags(cell_markers)
    domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
    file.write_meshtags(facet_markers)

gmsh.finalize()

tdim = domain.topology.dim
fdim = tdim - 1

# define inner boundry
metadata = {"quadrature_degree": 4}
ds = Measure('ds', domain=domain, subdomain_data=facet_markers, metadata=metadata)

# define function space
V = fem.VectorFunctionSpace(domain, ("CG",1))

# definition of circular shape
########## Insert 'circle' here ##########

# test and trial functions
u = fem.Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)

# define outer contour as homogenous boundary condition
u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_markers.find(outside)), V)

# constituive relation
E = fem.Constant(domain, ScalarType(10.))
nu = fem.Constant(domain, ScalarType(0.3))
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
def ppos(x):
    return (x+np.abs(x))/2.

pen = fem.Constant(domain, ScalarType(1e4))
form = inner(sigma(u), eps(u_))*dx
#form = inner(sigma(u), eps(u_))*dx +pen*dot(u_, ppos(u-circle))*ds(1)
J = derivative(form, u, du)

problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=J)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)

with XDMFFile(domain.comm, "rec_output.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(u)

I use FEniCSx 0.6.0 inside a Docker container.
Many thanks in advance!