penalty-checkpoint.ipynb 8.14 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212
{
 "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
}