NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4fa62fbd-0920-4990-892e-29ef20062ee1",
   "metadata": {},
   "source": [
    "# Elastic Pendulum\n",
    "\n",
    "$$\n",
    "\\textit{Here is presented a simple example of an elastic pendulum}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a9e13aaa-6c05-4759-b4f0-b4f25a685439",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6408ec0",
   "metadata": {},
   "source": [
    "\n",
    "The equation of Elastodynamics is:\n",
    "\\begin{align*}\n",
    "\\rho \\ddot u &= \\nabla \\cdot \\sigma + f^{\\text{ext}} \\\\\n",
    "\\sigma &= C : E \\\\\n",
    "E &= \\frac{1}{2} (\\nabla u + \\nabla u^T + \\nabla u ^T\\nabla u)\n",
    "\\end{align*}\n",
    "\n",
    "where $\\rho$ is the density, $u$ is the displacement, $\\sigma$ is the stress, $C$ is the elasticity tensor, $\\varepsilon$ is the strain, and $f^{\\text{ext}}$ is the external force.\n",
    "\n",
    "Suppose that at $t=0$ the displacement and the velocity are zero, moreover the external force is the gravity force $f^{\\text{ext}} = - \\vec e _y$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e35649b7-f3a2-4641-bf33-ab3d3aeb6e2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f53036120d7e406e83cf772964300ca2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from netgen.occ import *\n",
    "bar = Rectangle(0.1,1).Face()\n",
    "bar.edges.Max(Y).name=\"top\"\n",
    "bar.faces.name=\"bar\"\n",
    "bar.faces.maxh=0.05\n",
    "\n",
    "hole = Circle((0.05, 0), 0.1).Face()\n",
    "hole.faces.name=\"hole\"\n",
    "\n",
    "circ = Circle((0.05, 0), 0.2).Face()\n",
    "circ.faces.name=\"circ\"\n",
    "circ.faces.maxh=0.1\n",
    "\n",
    "shape = Glue([bar - (circ), circ - hole])\n",
    "\n",
    "shape = shape.Rotate(Axis((0.05, 1, 0), (0, 0, 1)), 180)\n",
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.1))\n",
    "mesh.Curve(4)\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2c0bddb",
   "metadata": {},
   "source": [
    "Given the disolacement, we define the Right Cauchy-Green tensor as:\n",
    "\n",
    "\\begin{align*}\n",
    "C = (I + \\nabla u)^T (I + \\nabla u)\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a225848b-237b-47a5-ac94-a3d4f05f24e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def C(u):\n",
    "    F = Id(2) + Grad(u)\n",
    "    return F.trans * F\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d06ad841",
   "metadata": {},
   "source": [
    "The Lamé parameters $\\lambda$ and $\\mu$ are derived from the Young's modulus $E$ and the Poisson's ratio $\\nu$ as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7433049d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# E, nu = 210, 0.2\n",
    "E, nu = 80, 0.4  # rubber\n",
    "\n",
    "mu = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcb762ab",
   "metadata": {},
   "source": [
    "We can define the hyperelastic material strain energy density function as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5c7f9009",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def NeoHooke (C):\n",
    "    return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86283920",
   "metadata": {},
   "source": [
    "The upper boundary is free to swing, to do so we control the $x$ and $y$ mean value of the upper boundary using a Number space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "17ef5205-2520-4739-8814-69b744f5ee17",
   "metadata": {},
   "outputs": [],
   "source": [
    "V = VectorH1(mesh, order=2)\n",
    "Q = NumberSpace(mesh, definedon=mesh.Boundaries('top'))\n",
    "fes = V * Q**2\n",
    "(u,q), (v, p) = fes.TnT()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1bada08",
   "metadata": {},
   "source": [
    "$u$ ... displacement, $v = \\dot u $ ... velocity, $a = \\dot v$ .. acceleration\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1c8bbfc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfut = GridFunction(V, multidim=0)\n",
    "\n",
    "# define the needed GridFunctiond required in the scheme\n",
    "gfu = GridFunction(fes)\n",
    "gfv = GridFunction(fes)\n",
    "gfa = GridFunction(fes)\n",
    "\n",
    "gfuold = GridFunction(fes)\n",
    "gfvold = GridFunction(fes)\n",
    "gfaold = GridFunction(fes)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c423eb06",
   "metadata": {},
   "source": [
    "with new acceleration, with elasticity operator $K$:\n",
    "\n",
    "$$\n",
    "a^{n+1} = f - K(u^{n+1})\n",
    "$$\n",
    "\n",
    "Here the unknowns are the displacement $u^{n+1}$ and the acceleration $a^{n+1}$, the velocity are to be determined via the time-stepping scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0b983054",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "bfa = BilinearForm(fes)\n",
    "bfa += Variation(NeoHooke(C(u))*dx) \n",
    "bfa += (InnerProduct(u, p) + InnerProduct(v, q)) * ds('top') # we add the constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93cd3fe7",
   "metadata": {},
   "source": [
    "\n",
    "Newmark scheme:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{u^{n+1}-u^n}{\\tau} = \\frac{v^n+v^{n+1}}{2} \\\\\n",
    "\\frac{v^{n+1}-v^n}{\\tau} = \\frac{a^n+a^{n+1}}{2} \\\\\n",
    "\\end{align*}\n",
    "\n",
    "That becomes:\n",
    "\\begin{align*}\n",
    "v^{n+1} = \\frac{2}{\\tau} (u^{n+1}-u^n) - v^n \\\\\n",
    "a^{n+1} = \\frac{2}{\\tau} (v^{n+1}-v^n) - a^n\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0125b77f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "tau = 0.025\n",
    "tend = 10\n",
    "force = CF((0, -1))\n",
    "\n",
    "\n",
    "vel_new = 2/tau * (u-gfuold.components[0]) - gfvold.components[0]\n",
    "acc_new = 2/tau * (vel_new-gfvold.components[0]) - gfaold.components[0]\n",
    "\n",
    "# need to add to the bilinear form since it depends on the current valurs of the GridFunctions\n",
    "\n",
    "bfa += acc_new*v*dx\n",
    "bfa += -force*v*dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "171c5690",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d2e1ad0daade47c4875dcad8d3267061",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ngsolve.solvers import Newton\n",
    "gfu.vec[:] = 0\n",
    "gfut.AddMultiDimComponent(gfu.components[0].vec)\n",
    "#scene = Draw(gfu.components[0], mesh, \"deformation\", deformation=True)  \n",
    "t = 0\n",
    "i =1\n",
    "with TaskManager():\n",
    "    while t < tend:\n",
    "        i += 1\n",
    "        t += tau\n",
    "        Newton(a=bfa, u=gfu, printing=False, inverse=\"sparsecholesky\")\n",
    "\n",
    "        if i % 5 == 0:\n",
    "            #scene.Redraw()\n",
    "            gfut.AddMultiDimComponent(gfu.components[0].vec)\n",
    "        gfv.vec[:] = 2/tau * (gfu.vec-gfuold.vec) - gfvold.vec\n",
    "        gfa.vec[:] = 2/tau * (gfv.vec-gfvold.vec) - gfaold.vec\n",
    "\n",
    "        gfuold.vec[:] = gfu.vec\n",
    "        gfvold.vec[:] = gfv.vec\n",
    "        gfaold.vec[:] = gfa.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2b46d75d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "986029d7fdc54531896d9d64bd6d28cd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "settings = {\"Multidim\": {\n",
    "    \"speed\" : 100\n",
    "}}\n",
    "\n",
    "#Draw(gfut, mesh, animate=True, min=0, \n",
    "#     max=0.5, autoscale=True, deformation=True)#, settings = settings);\n",
    "\n",
    "Draw(gfut, mesh, interpolate_multidim=True, animate=True, autoscale = False, min = 0, max = 2, deformation = True, speed = 100, **settings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9315bc4",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}