{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fluid-structure interaction with Taylor-Hood elements\n", "\n", "In this section we will use Taylor-Hood elements for discretizing the Navier-Stokes equations and Lagrangian elements for the elastic wave equations. We recall that, [see introduction of FSI](fsi_intro.ipynb), the Navier-Stokes equations in ALE form read\n", "\\begin{align*}\n", "&\\int_{\\Omega^f}J(\\rho\\frac{\\partial \\hat{u}}{\\partial t}\\cdot\\hat{v}+\\rho((\\hat{u}-\\dot{d})\\cdot\\nabla)\\hat{u}F^{-1}\\cdot\\hat{v}+2\\rho\\nu \\mathrm{sym}(\\nabla\\hat{u}F^{-1}):\\mathrm{sym}(\\nabla\\hat{v}F^{-1})-\\mathrm{tr}(\\nabla\\hat{v}F^{-1})p)\\,dx = 0&&\\qquad \\forall \\hat{v},\\\\\n", "&\\int_{\\Omega^f}J\\,\\mathrm{tr}(\\nabla\\hat{u}F^{-1})q\\,dx = 0&&\\qquad \\forall q.\n", "\\end{align*}\n", "and the elastic wave equation as first order system\n", "\\begin{align*}\n", "& \\int_{\\Omega^s}\\frac{\\partial d}{\\partial t}\\cdot v \\,dx = \\int_{\\Omega^s}u\\cdot v\\,dx &&\\qquad \\forall v,\\\\\n", "&\\int_{\\Omega^s}\\rho\\frac{\\partial u}{\\partial t}\\cdot w+(F\\Sigma):\\nabla w\\,dx = 0&&\\qquad \\forall w.\n", "\\end{align*}\n", "\n", "We will consider the following benchmark proposed in [Turek, Hron. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow. In: Fluid-Structure Interaction: Modelling, Simulation, Optimisation, 2006]. It is based on the classical flow around cylinder benchmark [Schäfer, Turek, Durst, Krause, Rannacher. Benchmark computations of laminar flow around a cylinder. In: Flow simulation with high-performance computers II, 1996], where additionally an elastic flag is \"glued\" behind the obstacle.\n", "\n", "\n", "\n", "We choose the space dependent function $h(x)$ in the deformation extension problem to be large close to the elastic solid's tip and decreases with the distance." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ngsolve import *\n", "from netgen.occ import *\n", "from ngsolve.webgui import Draw\n", "import ipywidgets as widgets\n", "\n", "tau = 0.004\n", "tend = 5\n", "order = 3\n", "\n", "rhos, nus, mus, rhof, nuf, U = 1e4, 0.4, 0.5 * 1e6, 1e3, 1e-3, 1\n", "ls = 2 * mus * nus / (1 - 2 * nus)\n", "\n", "# Parabolic inflow profile at the inlet\n", "par = Parameter(0)\n", "u_inflow = par * CoefficientFunction((4 * U * 1.5 * y * (0.41 - y) / (0.41 * 0.41), 0))\n", "\n", "\n", "# magnitude of inflow velocity over time\n", "def Force(t):\n", " if t < 0:\n", " return 0\n", " elif t < 2:\n", " return (1 - cos(pi / 2.0 * t)) / 2.0\n", " else:\n", " return 1\n", "\n", "\n", "# directly start with Stokes solution instead of increasing inflow\n", "start_with_stokes = True\n", "\n", "\n", "def GenerateMesh(order, maxh=0.2):\n", " circle = Circle((0.2, 0.2), r=0.05).Face()\n", " circle.edges.name = \"circ\"\n", " fluid = Rectangle(2.5, 0.41).Face()\n", " fluid.faces.name = \"fluid\"\n", " fluid.edges.Min(X).name = \"inlet\"\n", " fluid.edges.Max(X).name = \"outlet\"\n", " fluid.edges.Min(Y).name = \"wall\"\n", " fluid.edges.Max(Y).name = \"wall\"\n", " solid = (\n", " MoveTo(0.248989794855664, 0.19).Rectangle(0.6 - 0.248989794855664, 0.02).Face()\n", " )\n", " solid.faces.name = \"solid\"\n", " solid.edges.name = \"interface\"\n", " solid.edges.Min(X).name = \"circ_inner\"\n", "\n", " domain_fluid = (fluid - circle) - solid\n", " domain = Glue([domain_fluid, solid])\n", "\n", " mesh = Mesh(OCCGeometry(domain, dim=2).GenerateMesh(maxh=maxh))\n", " mesh.Curve(order)\n", "\n", " return mesh\n", "\n", "\n", "mesh = GenerateMesh(order=order)\n", "Draw(mesh);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the spatial discretization we will use the Taylor-Hood elements for the Navier-Stokes equations and also H1-conforming elements for the elastic wave equation. Thus, we can use one global space for the velocity and the displacement. With the definedon flag, we can tell the pressure space, that it lives only on the fluid domain." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "V = VectorH1(mesh, order=order, dirichlet=\"inlet|wall|circ|circ_inner\")\n", "Q = H1(mesh, order=order - 1, definedon=\"fluid\")\n", "D = VectorH1(mesh, order=order, dirichlet=\"inlet|wall|outlet|circ|circ_inner\")\n", "\n", "X = V * Q * D\n", "Y = V * Q\n", "(u, p, d), (v, q, w) = X.TnT()\n", "\n", "gf_solution = GridFunction(X)\n", "gf_solution_old = GridFunction(X)\n", "\n", "velocity, pressure, deformation = gf_solution.components\n", "velocity_old, pressure_old, deformation_old = gf_solution_old.components\n", "\n", "gradu_old = Grad(velocity_old)\n", "gradd_old = Grad(deformation_old)\n", "\n", "I = Id(mesh.dim)\n", "\n", "\n", "def CalcStresses(A):\n", " F = A + I\n", " C = F.trans * F\n", " E = 0.5 * (C - I)\n", " J = Det(F)\n", " Finv = Inv(F)\n", " return (F, C, E, J, Finv)\n", "\n", "\n", "F, C, E, J, Finv = CalcStresses(Grad(d))\n", "F_old, C_old, E_old, J_old, Finv_old = CalcStresses(gradd_old)\n", "\n", "\n", "def Stress(mat):\n", " return mus * mat + ls / 2 * Trace(mat) * I" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the time discretization we will use the Crank-Nicolson method\n", "\n", "\\begin{align*}\n", "\\int_{t^n}^{t^{n+1}} f(s)\\,ds \\approx \\frac{\\tau}{2}(f(t^{n+1})+f(t^n)).\n", "\\end{align*}\n", "\n", "Only the pressure constraint is handled with implicit Euler." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# For Stokes problem, check_unused=False to avoid warning not solving on the solid\n", "stokes = BilinearForm(Y, symmetric=True, check_unused=False)\n", "stokes += (\n", " nuf * rhof * 2 * InnerProduct(Sym(Grad(u)), Sym(Grad(v)))\n", " - div(u) * q\n", " - div(v) * p\n", " - 1e-8 * p * q\n", ") * dx(\"fluid\")\n", "stokes.Assemble()\n", "\n", "true_compile = False\n", "\n", "bfa = BilinearForm(X, symmetric=False, check_unused=False, condense=True)\n", "########################### Fluid: Navier-Stokes ##########################\n", "# M du/dt\n", "bfa += (rhof / tau * (InnerProduct(0.5 * (J + J_old) * (u - velocity_old), v))).Compile(\n", " true_compile, wait=True\n", ") * dx(\"fluid\")\n", "# symmetric stress div (eps u)\n", "bfa += (\n", " 0.5\n", " * rhof\n", " * nuf\n", " * (\n", " InnerProduct(J * 2 * Sym(Grad(u) * Finv), Sym(Grad(v) * Finv))\n", " + InnerProduct(J_old * 2 * Sym(gradu_old * Finv_old), (Grad(v) * Finv_old))\n", " )\n", ").Compile(true_compile, wait=True) * dx(\"fluid\")\n", "# Convection and mesh-velocity\n", "bfa += (\n", " 0.5\n", " * rhof\n", " * (\n", " InnerProduct(J * (Grad(u) * Finv) * (u - (d - deformation_old) / tau), v)\n", " + InnerProduct(\n", " J_old\n", " * (gradu_old * Finv_old)\n", " * (velocity_old - (d - deformation_old) / tau),\n", " v,\n", " )\n", " )\n", ").Compile(true_compile, wait=True) * dx(\"fluid\")\n", "# Pressure/Constraint implicit\n", "bfa += (-J * (Trace(Grad(v) * Finv) * p + Trace(Grad(u) * Finv) * q)).Compile(\n", " true_compile, wait=True\n", ") * dx(\"fluid\")\n", "\n", "########################### Solid: elastic wave ##########################\n", "# M du/dt\n", "bfa += (rhos / tau * InnerProduct(u - velocity_old, v)).Compile(\n", " true_compile, wait=True\n", ") * dx(\"solid\")\n", "# Material law\n", "bfa += (InnerProduct(F * Stress(E) + F_old * Stress(E_old), Grad(v))).Compile(\n", " true_compile, wait=True\n", ") * dx(\"solid\")\n", "# dd/dt = u\n", "bfa += (InnerProduct(u + velocity_old - 2.0 / tau * (d - deformation_old), w)).Compile(\n", " true_compile, wait=True\n", ") * dx(\"solid\")\n", "\n", "\n", "########################## Deformation extension ##########################\n", "def minCF(a, b):\n", " return IfPos(a - b, b, a)\n", "\n", "\n", "gf_dist = GridFunction(H1(mesh, order=2))\n", "gf_dist.Set(\n", " minCF(\n", " (x - 0.6) * (x - 0.6) + (y - 0.19) * (y - 0.19),\n", " (x - 0.6) * (x - 0.6) + (y - 0.21) * (y - 0.21),\n", " )\n", ")\n", "\n", "\n", "def NeoHookExt(C, mu=1, lam=1):\n", " return 0.5 * mu * (Trace(C - I) + 2 * mu / lam * Det(C) ** (-lam / 2 / mu) - 1)\n", "\n", "\n", "bfa += Variation(\n", " (1e-20 * mus * (1 / sqrt(gf_dist * gf_dist + 1e-12)) * NeoHookExt(C)).Compile(\n", " true_compile, wait=True\n", " )\n", " * dx(\"fluid\")\n", ")\n", "\n", "Draw(1 / sqrt(gf_dist * gf_dist + 1e-12), mesh, \"h(x)\", min=0.1, max=100, order=3);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To increase the inflow velocity depending on time, we have to extend the new velocity into the domain, before solving the system. This can be done by solving a Stokes problem with the new velocity as dirichlet data only on the fluid domain. To tell the solver on which domain it should work, we have to define the according degrees of freedom, which is done in terms of bitarrays." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bt_stokes = Y.FreeDofs() & ~Y.GetDofs(mesh.Materials(\"solid\"))\n", "bt_stokes &= ~Y.GetDofs(mesh.Boundaries(\"wall|inlet|circ|interface|circ_inner\"))\n", "# set all pressure dofs as active for Stokes\n", "bt_stokes[V.ndof :] = True\n", "\n", "gf_stokes = GridFunction(Y)\n", "res_stokes = gf_stokes.vec.CreateVector()\n", "inv_stokes = stokes.mat.Inverse(bt_stokes, inverse=\"sparsecholesky\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reset and draw." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = 0\n", "i = 0\n", "\n", "\n", "# Calculate quantities of interest\n", "def CalcForces(disp_x, disp_y):\n", " dmidx, dmidy = deformation(0.6, 0.2)\n", " disp_x.append(dmidx)\n", " disp_y.append(dmidy)\n", " return\n", "\n", "\n", "disp_x = [0]\n", "disp_y = [0]\n", "times = [0]\n", "\n", "if start_with_stokes:\n", " par.Set(1)\n", " gf_stokes.components[0].Set(u_inflow, definedon=mesh.Boundaries(\"inlet\"))\n", " res_stokes.data = stokes.mat * gf_stokes.vec\n", " gf_stokes.vec.data -= inv_stokes * res_stokes\n", " velocity.vec.data = gf_stokes.components[0].vec\n", " pressure.vec.data = gf_stokes.components[1].vec\n", "else:\n", " gf_solution.vec[:] = 0\n", "gf_solution_old.vec.data = gf_solution.vec\n", "\n", "scene_u = Draw(\n", " velocity,\n", " mesh.Materials(\"fluid\"),\n", " \"velocity\",\n", " deformation=deformation,\n", " order=3,\n", " max=2,\n", ")\n", "scene_p = Draw(pressure, mesh.Materials(\"fluid\"), \"pressure\", deformation=deformation)\n", "scene_d = Draw(deformation, mesh, \"deformation\")\n", "\n", "\n", "gf_history = GridFunction(X, multidim=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Time loop." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tw = widgets.Text(value=\"t = 0\")\n", "display(tw)\n", "\n", "with TaskManager():\n", " while t < tend - tau / 2.0:\n", " t += tau\n", "\n", " # update inflow by extending inflow profile as Stokes solution\n", " if t < 2 + tau / 2.0 and not start_with_stokes:\n", " par.Set(Force(t) - Force(t - tau))\n", " gf_stokes.components[0].Set(\n", " u_inflow, BND, definedon=mesh.Boundaries(\"inlet\")\n", " )\n", " res_stokes.data = stokes.mat * gf_stokes.vec\n", " gf_stokes.vec.data -= inv_stokes * res_stokes\n", " velocity.vec.data += gf_stokes.components[0].vec\n", " pressure.vec.data += gf_stokes.components[1].vec\n", "\n", " solvers.Newton(bfa, gf_solution, maxit=10, maxerr=1e-2, printing=False)\n", "\n", " if i % 10 == 0:\n", " scene_u.Redraw()\n", " scene_d.Redraw()\n", " scene_p.Redraw()\n", "\n", " if i % 20 == 0:\n", " gf_history.AddMultiDimComponent(gf_solution.vec)\n", "\n", " times.append(t)\n", " CalcForces(disp_x, disp_y)\n", "\n", " gf_solution_old.vec.data = gf_solution.vec\n", "\n", " i += 1\n", " tw.value = f\"t = {round(t,5)}\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Draw(\n", " gf_history.components[0],\n", " mesh.Materials(\"fluid\"),\n", " animate=True,\n", " min=0,\n", " max=2,\n", " autoscale=True,\n", " deformation=gf_history.components[2],\n", " order=3,\n", ");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Draw results over time. They become periodic after some time." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "plt.plot(times, disp_x)\n", "plt.xlabel(\"time\")\n", "plt.ylabel(\"disp_x\")\n", "plt.grid(True)\n", "plt.show()\n", "\n", "plt.plot(times, disp_y)\n", "plt.xlabel(\"time\")\n", "plt.ylabel(\"disp_y\")\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using a Taylor-Hood discretization for the Navier-Stokes equations has the draw back of not exactly divergence-free velocity fields\n", "\n", "\\begin{align*}\n", "\\int_{\\Omega}\\mathrm{div}(u)\\,q\\,dx = 0 \\quad\\forall q\\quad\\nRightarrow \\quad\\mathrm{div}(u)\\equiv 0.\n", "\\end{align*}\n", "\n", "Using H(div)-conforming Brezzi-Douglas-Marini or Raviart-Thomas finite elements yield to exactly divergence-free solutions. The ALE transformation and interface coupling becomes a bit more involved:\n", "[Neunteufel, Schöberl. Fluid-structure interaction with H(div)-conforming finite elements. Computers \\& Structures, (2021).], see [Notebook FSI with H(div)-conforming HDG](fsi_HDiv.ipynb)." ] }, { "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": 4 }