{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Hertzian contact with a rigid indenter using a penalty approach\n",
"\n",
"In this numerical tour, we explore the formulation of frictionless contact between a rigid surface (the indenter) and an elastic domain, representing an infinite half-space for the present case. Contact will be solved using a penalty formulation, that is small interpenetration between the solid and the identer will be authorized. Prerequisites for this tour are the formulation of an linear elasticity problem and the resolution of a general nonlinear variational problem."
]
},
{
"cell_type": "raw",
"metadata": {
"raw_mimetype": "text/restructuredtext"
},
"source": [
".. figure:: contact_problem.png\n",
" :scale: 20%\n",
" :align: center\n",
" \n",
" Hertzian contact problem of a rigid spherical indenter on a semi-infinite domain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The elastic (isotropic $E,\\nu$) solid will be represented by a 3D cubic domain of unit dimension and contact will take place on its top surface $z=0$, centered at $x=y=0$. Symmetry conditions will be applied on the $x=0$ and $y=0$ surfaces whereas the bottom surface $z=-1$ will be fully fixed. If contact appears on a small region of extent $a\\ll 1$, the problem can then be considered to be a good approximation of contact on a semi-infinite domain.\n",
"\n",
"The rigid indenter will not be explictly modeled (in terms of mesh in particular) but instead its distance with respect to the solid top surface will be given as an ``Expression``. We will also consider that the indenter radius $R$ is sufficiently large with respect to the contact region characteristic size $a$ so that the spherical surface can be approximated by a parabola. In this case, the distance between such an indenter and the top surface can be written as:\n",
"\n",
"\\begin{equation}\n",
"h(x,y) = h_0 + \\dfrac{1}{2R}(x^2+y^2)\n",
"\\end{equation}\n",
"\n",
"where $h_0$ is the initial gap between both surfaces at $x=y=0$. Obviously, if $h_0>0$ there is no contact between both surfaces, contact appears only if $h_0=-d<0$ where $d$ will be the indenter depth inside the surface. This classical problem admits the following known analytical solution **[REF?]**:\n",
"\n",
"* the contact area is of circular shape and radius $a=\\sqrt{Rd}$\n",
"* the force exerted by the indenter onto the surface is $F=\\dfrac{4}{3}\\dfrac{E}{1-\\nu^2}ad$\n",
"* the pressure distribution on the contact region is given by $p(r) = p_0\\sqrt{1-(r/a)^2}$ where $p_0=3F/(2\\pi a^2)$ is the maximal pressure\n"
]
},
{
"cell_type": "raw",
"metadata": {
"raw_mimetype": "text/restructuredtext"
},
"source": [
".. figure:: contact_solution.png\n",
" :scale: 20%\n",
" :align: center\n",
" \n",
" Indent at depth :math:`d` (left) and contact solution (right)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Contact problem formulation and penalty approach\n",
"\n",
"The unilateral contact (Signorini) condition on the top surface $\\Gamma$ writes as:\n",
"\n",
"\\begin{equation}\n",
"g \\geq 0, p\\geq 0, g\\cdot p =0 \\text{ on }\\Gamma\n",
"\\end{equation}\n",
"\n",
"where $g=h-u$ is the gap between the obstacle surface and the solid surface and $p$ is the pressure. One of the most simple way to solve approximately this contact condition consists in replacing the previous complementary conditions by the following penalized condition:\n",
"\n",
"\\begin{equation}\n",
"p = k\\ppos(-g)\n",
"\\end{equation}\n",
"\n",
"where $k$ is a large penalizing stiffness coefficient. With the previous relation, the pressure will be positive but a small negative gap will be authorized. "
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
"Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
"Maximum pressure FE: 1.42763078681 Hertz: 1.39916433487:\n",
"Applied force FE: 0.0320642723007 Hertz: 0.029304029304:\n"
]
}
],
"source": [
"from dolfin import *\n",
"import numpy as np\n",
"\n",
"N = 30\n",
"mesh = UnitCubeMesh.create(N, N, N/2, CellType.Type_hexahedron)\n",
"mesh.coordinates()[:, :2] = mesh.coordinates()[:, :2]**2\n",
"mesh.coordinates()[:, 2] = -mesh.coordinates()[:, 2]**2\n",
"\n",
"class Top(SubDomain):\n",
" def inside(self, x, on_boundary):\n",
" return near(x[2], 0.) and on_boundary\n",
"def symmetry_x(x, on_boundary):\n",
" return near(x[0], 0) and on_boundary\n",
"def symmetry_y(x, on_boundary):\n",
" return near(x[1], 0) and on_boundary\n",
"def bottom(x, on_boundary):\n",
" return near(x[2], -1) and on_boundary\n",
" \n",
"# exterior facets MeshFunction\n",
"facets = MeshFunction(\"size_t\", mesh, 2)\n",
"facets.set_all(0)\n",
"Top().mark(facets, 1)\n",
"ds = Measure('ds', subdomain_data=facets)\n",
"\n",
"R = 0.5\n",
"d = 0.02\n",
"obstacle = Expression(\"-d+(pow(x[0],2)+pow(x[1], 2))/2/R\", d=d, R=R, degree=2)\n",
"\n",
"E = Constant(10.)\n",
"nu = Constant(0.3)\n",
"mu = E/2/(1+nu)\n",
"lmbda = E*nu/(1+nu)/(1-2*nu)\n",
"\n",
"\n",
"V = VectorFunctionSpace(mesh, \"CG\", 1)\n",
"V2 = FunctionSpace(mesh, \"CG\", 1)\n",
"V0 = FunctionSpace(mesh, \"DG\", 0)\n",
"\n",
"u = Function(V, name=\"Displacement\")\n",
"du = TrialFunction(V)\n",
"u_ = TestFunction(V)\n",
"h = Function(V2)\n",
"gap = Function(V2, name=\"Gap\")\n",
"p = Function(V0, name=\"Contact pressure\")\n",
"\n",
"\n",
"def border(x, on_boundary):\n",
" return on_boundary\n",
"bc =[DirichletBC(V, Constant((0., 0., 0.)), bottom),\n",
" DirichletBC(V.sub(0), Constant(0.), symmetry_x),\n",
" DirichletBC(V.sub(1), Constant(0.), symmetry_y)]\n",
"\n",
"def eps(v):\n",
" return sym(grad(v))\n",
"\n",
"def sigma(v):\n",
" return lmbda*tr(eps(v))*Identity(3) + 2.0*mu*eps(v)\n",
"\n",
"def ppos(x):\n",
" return (x+abs(x))/2.\n",
"\n",
"pen = Constant(1e4)\n",
"form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[2], ppos(u[2]-obstacle))*ds(1)\n",
"J = derivative(form, u, du)\n",
"problem = NonlinearVariationalProblem(form, u, bc, J=J)\n",
"solver = NonlinearVariationalSolver(problem)\n",
"solver.parameters[\"newton_solver\"][\"linear_solver\"] = \"cg\"\n",
"solver.parameters[\"newton_solver\"][\"preconditioner\"] = \"ilu\"\n",
"\n",
"\n",
"h.interpolate(obstacle)\n",
"solver.solve()\n",
"\n",
"p.assign(-project(sigma(u)[2, 2], V0))\n",
"gap.assign(project(h-u[2], V2))\n",
"\n",
"#pl = plot(p)\n",
"#plt.colorbar(pl)\n",
"#plt.show()\n",
"\n",
"a = sqrt(R*d)\n",
"F = 4/3.*float(E)/(1-float(nu)**2)*a*d\n",
"p0 = 3*F/(2*pi*a**2)\n",
"print(\"Maximum pressure FE: {} Hertz: {}:\".format(max(np.abs(p.vector().get_local())), p0))\n",
"print(\"Applied force FE: {} Hertz: {}:\".format(4*assemble(p*ds(1)), F))\n",
"\n",
"file_results = XDMFFile(\"contact_penalty_results.xdmf\")\n",
"file_results.parameters[\"flush_output\"] = True\n",
"file_results.parameters[\"functions_share_mesh\"] = True\n",
"file_results.write(u, 0.)\n",
"file_results.write(gap, 0.)\n",
"file_results.write(p, 0.)"
]
}
],
"metadata": {
"celltoolbar": "Raw Cell Format",
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.12"
}
},
"nbformat": 4,
"nbformat_minor": 2
}