Unexpected results in the application of Hertzian contact tutorial

General
1

Hi, I went through this amazing tutorial on Hertzian contact in 3D spatial space and decided to implement it in 2D but the results of Contact Pressure plot and Displacement plot are quite unexpected when I try to indent on point (0.5, 1) .
Plots -

image
image1602×862 92.7 KB

Expectations - I want the contact pressure to be more in the area of contact like this -
image

I want displacement plot like this -

WhatsApp Image 2023-06-26 at 16.40.25
WhatsApp Image 2023-06-26 at 16.40.25899×1599 105 KB

How can I achieve both by tweaking the boundary conditions ? I would love to hear on the same.

Minimum Working Example -

from dolfin import *
import numpy as np
from mshr import *
import matplotlib.pyplot as plt

domain = Rectangle(Point(0., 0.), Point(1., 1.))
mesh = generate_mesh(domain, 10)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.) and on_boundary

def bottom(x, on_boundary):
        return near(x[1], 0.) and on_boundary

facets = MeshFunction("size_t", mesh, 2)
facets.set_all(0)
bc_top = Top()
bc_top.mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)

R = 0.5
d = 0.02
obstacle = Expression("-d+(pow(x[0] - 0.5,2) + 1)/2/R", d=d, R=R, degree=2)

V = VectorFunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 1)
V0 = FunctionSpace(mesh, "DG", 0)

u = Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)
gap = Function(V2, name="Gap")
p = Function(V0, name="Contact pressure")

bc =[DirichletBC(V, Constant((0., 0.)), bottom), DirichletBC(V.sub(1), obstacle, bc_top)]

E = Constant(10.)
nu = Constant(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+abs(x))/2.

pen = Constant(1e4)
form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[1], ppos(u[1]-obstacle))*ds(1)
J = derivative(form, u, du)

problem = NonlinearVariationalProblem(form, u, bc, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["linear_solver"] = "cg"
solver.parameters["newton_solver"]["preconditioner"] = "ilu"

solver.solve()

p.assign(-project(sigma(u)[1, 1], V0))
gap.assign(project(obstacle-u[1], V2))

a = sqrt(R*d)
F = 4/3.*float(E)/(1-float(nu)**2)*a*d
p0 = 3*F/(2*pi*a**2)
print("Maximum pressure FE: {0:8.5f} Hertz: {1:8.5f}".format(max(np.abs(p.vector().get_local())), p0))
print("Applied force    FE: {0:8.5f} Hertz: {1:8.5f}".format(4*assemble(p*ds(1)), F))

file_results = XDMFFile("2d_contact_penalty_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.write(gap, 0.)
file_results.write(p, 0.)

Thank You.