{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Reissner-Mindlin plate\n", "In this section we discretize and solve the Reissner-Mindlin plate (see [Reissner-Mindlin and Kirchhoff-Love plates](plates_derivation.ipynb) for a derivation) with two different methods and discuss the appearance of shear locking. The total energy consisting of bending and shearing reads with the vertical deflection $w$ and the linearized rotation vector $\\beta$ \n", "\n", "\\begin{align*}\n", "\\mathcal{W}_{\\mathrm{Reissner-Mindlin}}(w,\\beta) = \\frac{1}{2}\\int_{\\Omega}\\| \\varepsilon(\\beta) \\|_{\\mathbb{D}}^2 + \\frac{\\kappa\\,G}{t^2} \\| \\nabla w - \\beta \\|^2\\,dx - \\int_{\\Omega}f\\cdot w\\,dx,\n", "\\end{align*}\n", "\n", "where $t$ denotes the thickness of the plate, $\\varepsilon(\\cdot)$ is the symmetric part of the gradient and the elasticity tensor is $\\mathbb{D}\\varepsilon=\\frac{E}{12(1-\\nu^2)}((1-\\nu)\\varepsilon+\\nu\\,\\mathrm{tr}(\\varepsilon)I_{2\\times 2})$ (with $E$ Young's modulus, $\\nu$ Poisson ratio). Further, $G= \\frac{E}{2(1+\\nu)}$ and $\\kappa=5/6$ denote the shearing modulus and shear correction factor, respectively.\n", "\n", "We are especially interested in positive, but small thicknesses $t>0$. Therefore, the shear energy acts as a penalty enforcing the equality $\\beta = \\nabla w$. If there are not enough finite element functions such that the discrete so-called Kirchhoff constraint can be satisfied, we observe bad results known as shear locking." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider a benchmark example [Katili. A new discrete Kirchhoff-Mindlin element based on Mindlin-Reissner plate theory and assumed shear strain fields-Part I: An extended DKT element for thick-plate bending analysis. International Journal for Numerical Methods in Engineering 36, 11 (1993), 1859-1883] with exact solution. It consists of a circle domain, with clamped boundary and a uniform body force (e.g. gravity) is acting on it. Young modulus $E$, Poisson ratio $\\nu$, and shear correction factor $\\kappa$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ngsolve import *\n", "from netgen.occ import *\n", "from ngsolve.webgui import Draw\n", "\n", "E, nu, k = 10.92, 0.3, 5 / 6\n", "G = E / (2 * (1 + nu))\n", "fz = -1\n", "\n", "\n", "def DMat(mat, E, nu):\n", " return E / (12 * (1 - nu**2)) * ((1 - nu) * mat + nu * Trace(mat) * Id(2))\n", "\n", "\n", "def DMatInv(mat, E, nu):\n", " return (\n", " 12\n", " * (1 - nu**2)\n", " / E\n", " * (1 / (1 - nu) * mat - nu / (1 - nu**2) * Trace(mat) * Id(2))\n", " )\n", "\n", "\n", "R = 5\n", "\n", "\n", "def GenerateMesh(order, maxh=1):\n", " circ = Circle((0, 0), R).Face()\n", " circ.edges.name = \"circ\"\n", " return Mesh(OCCGeometry(circ, dim=2).GenerateMesh(maxh=maxh * R / 3)).Curve(order)\n", "\n", "\n", "def ExactSolution(thickness=0.1):\n", " r = sqrt(x**2 + y**2)\n", " xi = r / R\n", " Db = E * thickness**3 / (12 * (1 - nu**2))\n", "\n", " w_ex = (\n", " fz\n", " * thickness**3\n", " * R**4\n", " / (64 * Db)\n", " * (1 - xi**2)\n", " * ((1 - xi**2) + 8 * (thickness / R) ** 2 / (3 * k * (1 - nu)))\n", " )\n", " return w_ex\n", "\n", "\n", "mesh = GenerateMesh(order=2, maxh=1)\n", "w_ex = ExactSolution(thickness=0.1)\n", "Draw(w_ex, mesh, deformation=True, euler_angles=[-60, 5, 30]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simple discretization with Lagrange elements\n", "First we consider a simple discretization by using Lagrange elements for the displacements $w$ and rotations $\\beta$ of same polynomial degree $k$. Then we can directly discretize the Reissner-Mindlin plate equation: Find $(w,\\beta)\\in \\mathrm{Lag}^k\\times [\\mathrm{Lag}^k]^2$ with $w=\\beta=0$ on $\\partial\\Omega$ such that for all $(v,\\delta)\\in \\mathrm{Lag}^k\\times [\\mathrm{Lag}^k]^2$\n", "\n", "\\begin{align*}\n", "\\frac{1}{12}\\int_{\\Omega}\\mathbb{D}\\varepsilon(\\beta):\\varepsilon(\\delta)\\,dx + \\frac{\\kappa G}{t^2}\\int_{\\Omega}(\\nabla w-\\beta)\\cdot(\\nabla v-\\delta)\\,dx = \\int_{\\Omega}f\\,v\\,dx.\n", "\\end{align*}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def SolveReissnerMindlinLagrange(mesh, order=1, thickness=0.1, draw=False):\n", " fesB = VectorH1(mesh, order=order, dirichlet=\"circ\")\n", " fesW = H1(mesh, order=order, dirichlet=\"circ\")\n", " fes = fesB * fesW\n", " (beta, w), (delta, v) = fes.TnT()\n", "\n", " a = BilinearForm(fes, symmetric=True)\n", " # bending part\n", " a += InnerProduct(DMat(Sym(Grad(beta)), E, nu), Sym(Grad(delta))) * dx\n", " # shearing part\n", " a += k * G / thickness**2 * InnerProduct(Grad(w) - beta, Grad(v) - delta) * dx\n", " a.Assemble()\n", "\n", " f = LinearForm(fz * v * dx).Assemble()\n", "\n", " gf_solution = GridFunction(fes)\n", " _, gf_w = gf_solution.components\n", "\n", " inv = a.mat.Inverse(fes.FreeDofs(), inverse=\"sparsecholesky\")\n", " gf_solution.vec.data = inv * f.vec\n", "\n", " if draw:\n", " Draw(gf_w, mesh, \"w\")\n", "\n", " return gf_w, fes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve for different thicknesses for a sequence of refined meshes. Different polynomial order can be tested." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results = []\n", "thicknesses = [1, 0.1, 0.01, 0.001]\n", "# try out polynomial order 1 and 2\n", "order = 2\n", "\n", "with TaskManager():\n", " for t in thicknesses:\n", " results.append([])\n", " w_ex = ExactSolution(thickness=t)\n", "\n", " for i in range(5):\n", " mesh = GenerateMesh(order=order, maxh=0.5**i)\n", " gfw, fes = SolveReissnerMindlinLagrange(\n", " mesh, order=order, thickness=t, draw=False\n", " )\n", "\n", " # relative L2-error of displacement\n", " err = sqrt(Integrate((gfw - w_ex) * (gfw - w_ex), mesh)) / sqrt(\n", " Integrate(w_ex * w_ex, mesh)\n", " )\n", " results[-1].append((fes.ndof, err))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot error. We observe immense shear locking for polynomial order $k=1$ as the pre-asymptotic regime increases rapidly. Increasing the polynomial order to $k=2$ (or $k=3$) mitigates the problem, but it is still present." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "plt.yscale(\"log\")\n", "plt.xscale(\"log\")\n", "plt.xlabel(\"ndof\")\n", "plt.ylabel(\"relative error\")\n", "for i, result in enumerate(results):\n", " ndof, err = zip(*result)\n", " plt.plot(ndof, err, \"-*\", label=\"thickness=\" + str(thicknesses[i]))\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## TDNNS method\n", "\n", "Next, we apply the [TDNNS method](TDNNS_intro.ipynb) to discretize the bending term $\\int_{\\Omega}\\| \\varepsilon(\\beta) \\|_{\\mathbb{D}}^2$ as in [Pechstein and Schöberl. The TDNNS method for Reissner-Mindlin plates. Numerische Mathematik 137, 3 (2017), 713-740]. \n", "\n", " \n", " \n", "\n", "Instead of using Lagrangian elements, we discretize the rotations $\\beta$ in the Nedelec space. Thus, all discrete gradients $\\nabla w$ are in the space of rotations, and we don't have the locking problem from before. By introducing the linearized bending moment tensor $\\sigma=\\mathbb{D}\\varepsilon(\\beta)$ we obtain the following formulation.\n", "\n", "Find $(\\sigma,\\beta, w) \\in H(\\mathrm{div} \\mathrm{div})\\times H(\\mathrm{curl})\\times H^1(\\Omega)$ such that for all $(\\tau,\\delta, v) \\in H(\\mathrm{div} \\mathrm{div})\\times H(\\mathrm{curl})\\times H^1(\\Omega)$\n", "\n", "\\begin{align*}\n", "\\begin{array}{ccccll}\n", "-\\int_{\\Omega} \\mathbb{D}^{-1} \\sigma : \\tau\\,dx & - & \\left< \\mathrm{div} \\tau, \\beta \\right> & = & 0 \\\\\n", "-\\left< \\mathrm{div} \\sigma, \\delta \\right> & + & \\frac{\\kappa\\, G}{t^2} \\int_{\\Omega} (\\nabla w - \\beta) (\\nabla v - \\delta)\\,dx & = & \\int_{\\Omega} f v\\,dx,\n", "\\end{array}\n", "\\end{align*}\n", "\n", "with the TDNNS duality pairing\n", "\n", "\\begin{align*}\n", "\\langle \\mathrm{div} \\boldsymbol{\\sigma}, \\beta\\rangle &= \\sum_{T\\in\\mathcal{T}}\\int_T\\mathrm{div}\\boldsymbol{\\sigma}\\cdot\\beta\\,dx -\\int_{\\partial T}\\boldsymbol{\\sigma}_{nt}\\beta_t\\,ds\\\\\n", "&=-\\sum_{T\\in\\mathcal{T}}\\int_T\\boldsymbol{\\sigma}:\\nabla \\beta\\,dx +\\int_{\\partial T}\\boldsymbol{\\sigma}_{nn}\\beta_n\\,ds =-\\langle \\boldsymbol{\\sigma}, \\nabla \\beta\\rangle\n", "\\end{align*}\n", "\n", "and the inverted elasticity tensor $\\mathbb{D}^{-1}\\varepsilon = \\frac{12(1-\\nu^2)}{Et^3}\\big(\\frac{1}{1-\\nu}\\varepsilon-\\frac{\\nu}{1-\\nu^2}\\mathrm{tr}(\\varepsilon)I_{2\\times2}\\big)$. Here $\\partial T$ denotes the element-boundary (edges) of $T$, $n$ the normal vector and $t$ the tangential vector on $\\partial T$. With e.g. $\\sigma_{nt}:= t^\\top\\sigma n$ we denote the normal tangential component of $\\sigma$.\n", "\n", "We use Lagrangian elements of order $k$ for the displacement $w$, and one order less for the rotations $\\beta$ and the bending moment tensor $\\sigma$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def SolveReissnerMindlinTDNNS(mesh, order=1, thickness=0.1, draw=False):\n", " fesB = HCurl(mesh, order=order - 1, dirichlet=\"circ\")\n", " fesS = HDivDiv(mesh, order=order - 1, dirichlet=\"\")\n", " fesW = H1(mesh, order=order, dirichlet=\"circ\")\n", "\n", " fes = fesW * fesB * fesS\n", " (w, beta, sigma), (v, delta, tau) = fes.TnT()\n", "\n", " n = specialcf.normal(2)\n", "\n", " a = BilinearForm(fes, symmetric=True)\n", " # bending part\n", " a += (\n", " -InnerProduct(DMatInv(sigma, E, nu), tau)\n", " + InnerProduct(tau, Grad(beta))\n", " + InnerProduct(sigma, Grad(delta))\n", " ) * dx\n", " a += -(sigma[n, n] * delta * n + tau[n, n] * beta * n) * dx(element_boundary=True)\n", " # shearing part\n", " a += k * G / thickness**2 * InnerProduct(Grad(w) - beta, Grad(v) - delta) * dx\n", " a.Assemble()\n", "\n", " f = LinearForm(fz * v * dx).Assemble()\n", "\n", " gf_solution = GridFunction(fes)\n", " gf_w, _, _ = gf_solution.components\n", "\n", " inv = a.mat.Inverse(fes.FreeDofs(), inverse=\"\")\n", " gf_solution.vec.data = inv * f.vec\n", " if draw:\n", " Draw(gf_w, mesh, \"w\")\n", "\n", " return gf_w, fes" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results = []\n", "\n", "# try out polynomial order 1 and 2\n", "order = 1\n", "\n", "thicknesses = [1, 0.1, 0.01, 0.001]\n", "with TaskManager():\n", " for t in thicknesses:\n", " results.append([])\n", " w_ex = ExactSolution(thickness=t)\n", "\n", " for i in range(5):\n", " mesh = GenerateMesh(order=order, maxh=0.5**i)\n", " gfw, fes = SolveReissnerMindlinTDNNS(\n", " mesh, order=order, thickness=t, draw=False\n", " )\n", "\n", " # relative L2-error of displacement\n", " err = sqrt(Integrate((gfw - w_ex) * (gfw - w_ex), mesh)) / sqrt(\n", " Integrate(w_ex * w_ex, mesh)\n", " )\n", " results[-1].append((fes.ndof, err))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the TDNNS method we obtain a locking free discretization already for polynomial degree $k=1$ for the displacements." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "plt.yscale(\"log\")\n", "plt.xscale(\"log\")\n", "plt.xlabel(\"ndof\")\n", "plt.ylabel(\"relative error\")\n", "for i, result in enumerate(results):\n", " ndof, err = zip(*result)\n", " plt.plot(ndof, err, \"-*\", label=\"t=\" + str(thicknesses[i]))\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the limit $t \\to 0$ the shear energy term can be understood as a penalty formulation enforcing the Kirchhoff-Love assumption\n", "\n", "\\begin{align*}\n", "\\beta = \\nabla w.\n", "\\end{align*}\n", "\n", "Thus, in the limit case, the total energy simplifies by eliminating the rotation $\\beta$ to\n", "\n", "\\begin{align*}\n", "\\mathcal{W}_{\\mathrm{Kirchhoff-Love}}(w)=\\frac{1}{2}\\int_{\\Omega}\\| \\varepsilon(\\nabla w) \\|_{\\mathbb{D}}^2-\\int_{\\Omega}f\\,w,\n", "\\end{align*}\n", "\n", "which is the Kirchhoff-Love plate model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Drawback of the current implementation of the TDNNS formulation for Reissner-Mindlin plates: We have a mixed saddle point problem leading to an indefinite system matrix. This can be easily overcome with hybridization techniques [NGSolve docu - Static condensation](https://docu.ngsolve.org/latest/i-tutorials/unit-1.4-staticcond/staticcond.html#1.4-Static-Condensation) by breaking the normal-normal continuity of $\\sigma$ and reinforcing it by a Lagrange multiplier $\\alpha$. After static condensation a symmetric positive definite system is regained.\n", "\n", "Exercise:\n", "* Implement hybridization for TDNNS for Reissner-Mindlin plates" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 2 }