{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 6.1.3 Reissner-Mindlin plate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To avoid shear locking when the thickness $t$ becomes small several methods and elements have been proposed. We discuss two approaches:\n", "\n", "1: Mixed Interpolation of Tensorial Components (MITC)\n", "\n", "2: Rotations in Nèdèlec space using the Tangential-Displacement Normal-Normal-Stress method (TDNNS)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mixed Interpolation of Tensorial Components (MITC)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ngsolve import *\n", "from netgen.geom2d import SplineGeometry\n", "from ngsolve.webgui import Draw" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up benchmark example with exact solution. Young modulus $E$, Poisson ratio $\\nu$, shear correction factor $\\kappa$ and corresponding Lamè parameters $\\mu$ and $\\lambda$. Geometry parameters are given by thickness $t$ and radius $R$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "E, nu, k = 10.92, 0.3, 5/6\n", "mu = E/(2*(1+nu))\n", "lam = E*nu/(1-nu**2)\n", "#shearing modulus\n", "G = E/(2*(1+nu))\n", "#thickness, shear locking with t=0.1\n", "t = 0.1\n", "R = 5\n", "#force\n", "fz = 1\n", "\n", "order = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Due to symmetry we only need to mesh one quarter of the circle." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sg = SplineGeometry()\n", "pnts = [ (0,0), (R,0), (R,R), (0,R) ]\n", "pind = [ sg.AppendPoint(*pnt) for pnt in pnts ]\n", "sg.Append(['line',pind[0],pind[1]], leftdomain=1, rightdomain=0, bc=\"bottom\")\n", "sg.Append(['spline3',pind[1],pind[2],pind[3]], leftdomain=1, rightdomain=0, bc=\"circ\")\n", "sg.Append(['line',pind[3],pind[0]], leftdomain=1, rightdomain=0, bc=\"left\")\n", "mesh = Mesh(sg.GenerateMesh(maxh=R/3))\n", "mesh.Curve(order)\n", "Draw(mesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Depending on the boundary conditions (simply supported or clamped) we have different exact solutions for the vertical displacement. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "r = sqrt(x**2+y**2)\n", "xi = r/R\n", "Db = E*t**3/(12*(1-nu**2))\n", "\n", "clamped = True\n", "\n", "#exact solution for simply supported bc\n", "w_s_ex = -fz*R**4/(64*Db)*(1-xi**2)*( (6+2*nu)/(1+nu) - (1+xi**2) + 8*(t/R)**2/(3*k*(1-nu)))\n", "#Draw(w_s_ex, mesh, \"w_s_ex\")\n", "#exact solution for clamped bc\n", "w_c_ex = -fz*R**4/(64*Db)*(1-xi**2)*( (1-xi**2) + 8*(t/R)**2/(3*k*(1-nu)))\n", "Draw(w_c_ex, mesh, \"w_c_ex\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use (lowest order) Lagrangian finite elements for rotations and the vertical deflection together with additional internal bubbles for order >1." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if clamped:\n", " fesB = VectorH1(mesh, order=order, orderinner=order+1, dirichletx=\"circ|left\", dirichlety=\"circ|bottom\", autoupdate=True)\n", "else:\n", " fesB = VectorH1(mesh, order=order, orderinner=order+1, dirichletx=\"left\", dirichlety=\"bottom\", autoupdate=True)\n", "fesW = H1(mesh, order=order, orderinner=order+1, dirichlet=\"circ\", autoupdate=True)\n", "fes = FESpace( [fesB, fesW], autoupdate=True ) \n", "(beta,u), (dbeta,du) = fes.TnT()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Direct approach where shear locking may occur\n", "$$\\frac{t^3}{12}\\int_{\\Omega} 2\\mu\\, \\varepsilon(\\beta):\\varepsilon(\\delta\\beta) + \\lambda\\, \\text{div}(\\beta)\\text{div}(\\delta\\beta)\\,dx + t\\kappa\\,G\\int_{\\Omega}(\\nabla u-\\beta)\\cdot(\\nabla\\delta u-\\delta\\beta)\\,dx = \\int_{\\Omega} f\\,\\delta u\\,dx,\\qquad \\forall (\\delta u,\\delta\\beta). $$\n", "\n", "Adding interpolation (reduction) operator $\\boldsymbol{R}:[H^1_0(\\Omega)]^2\\to H(\\text{curl})$. Spaces are chosen according to [Brezzi, Fortin and Stenberg. Error analysis of mixed-interpolated elements for Reissner-Mindlin plates. Mathematical Models and Methods in Applied Sciences 1, 2\n", " (1991), 125-151.]\n", "$$\\frac{t^3}{12}\\int_{\\Omega} 2\\mu\\, \\varepsilon(\\beta):\\varepsilon(\\delta\\beta) + \\lambda\\, \\text{div}(\\beta)\\text{div}(\\delta\\beta)\\,dx + t\\kappa\\,G\\int_{\\Omega}\\boldsymbol{R}(\\nabla u-\\beta)\\cdot\\boldsymbol{R}(\\nabla\\delta u-\\delta\\beta)\\,dx = \\int_{\\Omega} f\\,\\delta u\\,dx,\\qquad \\forall (\\delta u,\\delta\\beta). $$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Gamma = HCurl(mesh,order=order,orderedge=order-1, autoupdate=True)\n", "\n", "a = BilinearForm(fes)\n", "a += t**3/12*(2*mu*InnerProduct(Sym(grad(beta)),Sym(grad(dbeta))) + lam*div(beta)*div(dbeta))*dx\n", "#a += t*k*G*InnerProduct( grad(u)-beta, grad(du)-dbeta )*dx\n", "a += t*k*G*InnerProduct(Interpolate(grad(u)-beta,Gamma), Interpolate(grad(du)-dbeta,Gamma))*dx\n", "\n", "f = LinearForm(fes)\n", "f += -fz*du*dx\n", "\n", "gfsol = GridFunction(fes, autoupdate=True)\n", "gfbeta, gfu = gfsol.components" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define function for solving the problem." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def SolveBVP():\n", " fes.Update()\n", " gfsol.Update()\n", " with TaskManager():\n", " a.Assemble()\n", " f.Assemble()\n", " inv = a.mat.Inverse(fes.FreeDofs(), inverse=\"sparsecholesky\")\n", " gfsol.vec.data = inv * f.vec" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve, compute error, refine, ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l = []\n", "for i in range(5):\n", " print(\"i = \", i)\n", " SolveBVP()\n", " if clamped:\n", " norm_w = sqrt(Integrate(w_c_ex*w_c_ex, mesh))\n", " err = sqrt(Integrate((gfu-w_c_ex)*(gfu-w_c_ex), mesh))/norm_w\n", " else:\n", " norm_w = sqrt(Integrate(w_s_ex*w_s_ex, mesh))\n", " err = sqrt(Integrate((gfu-w_s_ex)*(gfu-w_s_ex), mesh))/norm_w\n", " print(\"err = \", err)\n", " l.append ( (fes.ndof, err ))\n", " mesh.Refine()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convergence plot with matplotlib." ] }, { "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", "ndof,err = zip(*l)\n", "plt.plot(ndof,err, \"-*\")\n", "plt.ion()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## TDNNS method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Invert material law of Hooke and reset mesh." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def CMatInv(mat, E, nu):\n", " return (1+nu)/E*(mat-nu/(nu+1)*Trace(mat)*Id(2))\n", "\n", "mesh = Mesh(sg.GenerateMesh(maxh=R/3))\n", "mesh.Curve(order)\n", "Draw(mesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead of using Lagrangian elements together with an interpolation operator we can directly use H(curl) for the rotation $\\beta$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "order=1\n", "if clamped:\n", " fesB = HCurl(mesh, order=order-1, dirichlet=\"circ\", autoupdate=True)\n", " fesS = HDivDiv(mesh, order=order-1, dirichlet=\"\", autoupdate=True) \n", "else:\n", " fesB = HCurl(mesh, order=order-1)\n", " fesS = HDivDiv(mesh, order=order-1, dirichlet=\"circ\", autoupdate=True)\n", "fesW = H1(mesh, order=order, dirichlet=\"circ\", autoupdate=True)\n", "\n", "fes = FESpace( [fesW, fesB, fesS], autoupdate=True ) \n", "(u,beta,sigma), (du,dbeta,dsigma) = fes.TnT()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use the TDNNS method with the stress (bending) tensor $\\boldsymbol{\\sigma}$ as proposed in [Pechstein and Schöberl. The TDNNS method for Reissner-Mindlin plates. Numerische Mathematik 137, 3 (2017), 713-740]\n", "\\begin{align*}\n", "&\\mathcal{L}(u,\\beta,\\boldsymbol{\\sigma})=-\\frac{6}{t^3}\\|\\boldsymbol{\\sigma}\\|^2_{\\mathcal{C}^{-1}} + \\langle \\boldsymbol{\\sigma}, \\nabla \\beta\\rangle + \\frac{tkG}{2}\\|\\nabla u-\\beta\\|^2 -\\int_{\\Omega} f\\cdot u\\,dx\\\\\n", "&\\langle \\boldsymbol{\\sigma}, \\nabla \\beta\\rangle:= \\sum_{T\\in\\mathcal{T}}\\int_T\\boldsymbol{\\sigma}:\\nabla \\beta\\,dx -\\int_{\\partial T}\\boldsymbol{\\sigma}_{nn}\\beta_n\\,ds\n", "\\end{align*}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = specialcf.normal(2)\n", " \n", "a = BilinearForm(fes)\n", "a += (-12/t**3*InnerProduct(CMatInv(sigma, E, nu),dsigma) + InnerProduct(dsigma,grad(beta)) + InnerProduct(sigma,grad(dbeta)))*dx\n", "a += ( -((sigma*n)*n)*(dbeta*n) - ((dsigma*n)*n)*(beta*n) )*dx(element_boundary=True)\n", "a += t*k*G*InnerProduct( grad(u)-beta, grad(du)-dbeta )*dx\n", "\n", "f = LinearForm(fes)\n", "f += -fz*du*dx\n", "\n", "gfsol = GridFunction(fes, autoupdate=True)\n", "gfu, gfbeta, gfsigma = gfsol.components" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def SolveBVP():\n", " fes.Update()\n", " gfsol.Update()\n", " with TaskManager():\n", " a.Assemble()\n", " f.Assemble()\n", " inv = a.mat.Inverse(fes.FreeDofs())\n", " gfsol.vec.data = inv * f.vec" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve, compute error, and refine." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l = []\n", "for i in range(5):\n", " print(\"i = \", i)\n", " SolveBVP()\n", " if clamped:\n", " norm_w = sqrt(Integrate(w_c_ex*w_c_ex, mesh))\n", " err = sqrt(Integrate((gfu-w_c_ex)*(gfu-w_c_ex), mesh))/norm_w\n", " else:\n", " norm_w = sqrt(Integrate(w_s_ex*w_s_ex, mesh))\n", " err = sqrt(Integrate((gfu-w_s_ex)*(gfu-w_s_ex), mesh))/norm_w\n", " print(\"err = \", err)\n", " l.append ( (fes.ndof, err ))\n", " mesh.Refine()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convergence plot with matplotlib." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.yscale('log')\n", "plt.xscale('log')\n", "plt.xlabel(\"ndof\")\n", "plt.ylabel(\"relative error\")\n", "ndof,err = zip(*l)\n", "plt.plot(ndof,err, \"-*\")\n", "\n", "plt.ion()\n", "plt.show()" ] }, { "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.8.10" } }, "nbformat": 4, "nbformat_minor": 2 }