{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fluid-structure interaction with H(div)-conforming HDG\n", "\n", "In this section we will use $H(\\mathrm{div})$-conforming HDG elements for discretizing the Navier-Stokes equations and Lagrangian elements for the elastic wave equations following [Neunteufel, Schöberl. Fluid-structure interaction with H(div)-conforming finite elements. Computers \\& Structures, (2021).]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## H(div)-conforming HDG for Navier-Stokes with up-winding\n", "We quickly recap the H(div)-conforming HDG method for Navier-Stokes with up-winding for the convection term, as we will need to reformulate it in terms of an ALE formulation, for details see e.g. [NGSolve docu DG/HDG splitting](https://docu.ngsolve.org/latest/i-tutorials/unit-3.7-opsplit/opsplit.html), [Notebook Stokes HDG](../CFD/stokes_hdg.ipynb), [Lehrenfeld, Schöberl. High order exactly divergence-free Hybrid Discontinuous Galerkin Methods for unsteady incompressible flows. Computer Methods in Applied Mechanics and Engineering, (2016).].\n", "\n", "To this end, let $u=(u_T,u_F)$, $v=(v_T,u_F)$ with $u_T,v_T$ in the H(div)-conforming $BDM^k$ space and $u_F$, $v_F$ in the vector-valued tangential facet space $F_h^k$. We define the tangential jump $[\\![ v^\\tau]\\!] = v_T^\\tau-v_F^\\tau$ with $u^\\tau= u - (u\\cdot n)n$ the tangential projection. Further, $p$ is in the space of piece-wise polynomials $Q_h^{k-1}$. Then the diffusion, mass, and divergence terms of the Navier-Stokes equations are given by\n", "\n", "\\begin{align*}\n", "& A_h^f(u,v)=\\sum_{T\\in \\mathcal{T}^f}\\Big(\\int_T2\\nu\\varepsilon(u_T):\\nabla v_T\\,dx+\\int_{\\partial T}\\big(-2\\nu\\varepsilon(u_T)n\\cdot[\\![ v^\\tau]\\!]- 2\\nu\\varepsilon(v_T)n\\cdot[\\![ u^\\tau]\\!]+\\frac{\\nu\\alpha k^2}{h}[\\![ u^{\\tau}]\\!]\\cdot[\\![ v^{\\tau}]\\!]\\big)\\,ds\\Big),\\\\\n", "& M^f_h(u,v)=\\int_{\\Omega_f} u_T\\cdot v_T\\,dx,\\\\\n", "& D^f_h(v,p)=-\\int_{\\Omega_f} p\\,\\mathrm{div}(v_T)\\,dx,\n", "\\end{align*}\n", "\n", "For the convection we use upwinding \n", "\n", "\\begin{align*}\n", "C^f_h({u},v)= \\sum_{T\\in\\mathcal{T}^f}\\Big(-\\int_T\\nabla v_T{u}_T\\cdot{u}_T\\,dx+\\int_{\\partial T} {u}_T\\cdot {n}\\,{u}^{\\mathrm{up}}\\cdot v_T\\,ds+\\int_{\\partial T_{\\mathrm{out}}} {u}_T\\cdot {n}\\,({u}_F-{u}_T)^{{\\tau}}\\cdot v_F\\,ds\\Big),\n", "\\end{align*}\n", "\n", "with the upwind variable ${u}^{\\mathrm{up}}$ defined as\n", "\n", "\\begin{align*}\n", "& {u}^{\\mathrm{up}}= ({u}_T\\cdot {n}){n}+\\begin{cases}\n", "{u}_T^{{\\tau}} & \\text{ if }{u}\\cdot {n} \\geq 0\\\\\n", "{u}_F^{{\\tau}} & \\text{ if } {u}\\cdot {n} < 0\n", "\\end{cases}.\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Navier Stokes in ALE (Arbitrary Lagrangian Eulerian) formulation\n", "\n", "Let $\\Phi(x,t)=\\mathrm{id} + d$\n", "describe the movement of the mesh, where $d$ is the displacement. For Taylor-Hood elements we choose the relation of the velocity on the deformed and initial configuration $u\\circ\\Phi=\\hat{u}$, i.e., the continuity of $u$ is preserved if $\\hat{u}$ is continuous. For a velocity field in $H(\\mathrm{div})$, however, we want to preserve the normal-continuity to obtain an exactly divergence-free solution on both configurations. Therefore, we use the Piola transformation\n", "\n", "\\begin{align*}\n", "u\\circ\\Phi =P_{\\Phi}[\\hat{u}]= \\frac{1}{J} F \\hat{u},\\qquad F=\\nabla \\Phi,\\quad J=\\det(F).\n", "\\end{align*}\n", "\n", "As a result the chain rule now reads\n", "\n", "\\begin{align*}\n", "\\nabla_xu\\circ\\Phi = \\nabla_{\\hat{x}} (P_{\\Phi}[\\hat{u}])F^{-1},\n", "\\end{align*}\n", "\n", "\n", "\n", "\n", "\n", "The Piola transformation also depends on the time. Thus we obtain additional terms from the chain rule of the time derivative\n", "\n", "\\begin{align*}\n", "\\frac{\\partial u}{\\partial t}\\circ\\Phi = \\frac{1}{J}(\\nabla_{\\hat{x}}\\dot{d}-\\mathrm{tr}(\\nabla_{\\hat{x}}\\dot{d}F^{-1})F)\\hat{u}+P_{\\Phi}[\\frac{\\partial \\hat{u}}{\\partial t}]- \\nabla_{\\hat{x}} P_{\\Phi}[\\hat{u}]F^{-1}\\dot{d}\n", "\\end{align*}\n", "\n", "\n", "To avoid implementing all the chain rules involving the deformation gradient and Jacobian determinants, see [Neunteufel, Schöberl. Fluid-structure interaction with H(div)-conforming finite elements. Computers \\& Structures, (2021).] for details, we use the ``dx(deformation=gf_d)`` flag where ``gf_d`` is the GridFunction of the deformation. Then NGSolve automatically adds e.g. the deformation gradients for the chain rules. Only the additional terms coming from the chain rule of the time derivative has to be added by hand." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interface conditions\n", "To guarantee that the stresses from the fluid and solid are in equilibrium at the interface we need to force the continuity of the velocity at the interface as we use Lagrangian elements for the solid and $H(\\mathrm{div})$-conforming HDG for the fluid (we still use one global Lagrange space for the displacement and deformation extension). Therefore, we use a Langrange multiplier $\\mu = (\\mu_1,\\mu_2)$ at the interface to enforce the continuity of the velocity\n", "\n", "\\begin{align*}\n", "&\\int_{\\Gamma_I} ({u}_T^f-u^s)^{n}\\,\\lambda_1\\,ds+\\int_{\\Gamma_I} (v_T^{f}-v^{s})^{n}\\,{\\mu}_1\\,ds=0,\\\\\n", "&\\int_{\\Gamma_I}({u}_F^f-u^s)^{t}\\,\\lambda_2\\,ds +\\int_{\\Gamma_I} (v_F^{f}-v^{s})^{t}\\,\\mu_2\\, ds =0,\n", "\\end{align*}\n", "which have to be transformed again in the ALE setting." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Benchmark example\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", "import netgen.geom2d as csg\n", "from ngsolve.webgui import Draw\n", "import ipywidgets as widgets\n", "\n", "\n", "tau = 0.004\n", "tend = 5\n", "order = 3\n", "\n", "# HDG stabilization parameter\n", "alpha = Parameter(5)\n", "\n", "# solid density, Poisson ratio, first Lame parameter, fluid density,\n", "# dynamic viscosity, and maximal inflow velocity\n", "rhos, nus, mus, rhof, nuf, U = 1e4, 0.4, 0.5 * 1e6, 1e3, 1e-3, 1\n", "# second Lame parameter for solid\n", "ls = 2 * mus * nus / (1 - 2 * nus)\n", "\n", "# inflow profile\n", "u_inflow = CF((4 * U * 1.5 * y * (0.41 - y) / (0.41 * 0.41), 0))\n", "\n", "\n", "def GenerateMesh(order, maxh=0.203):\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", "\n", " domain_fluid = (fluid - circle) - solid\n", " solid.edges[\"circ\"].name = \"circ_inner\"\n", " solid.edges.Min(X).name = \"circ_inner\"\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": [ "V1 = VectorH1(mesh, order=order, dirichlet=\"wall|outlet|inlet|circ|circ_inner\")\n", "V2 = HDiv(\n", " mesh,\n", " order=order,\n", " dirichlet=\"wall|inlet|circ|circ_inner\",\n", " definedon=\"fluid\",\n", " hodivfree=True,\n", ")\n", "V3 = TangentialFacetFESpace(\n", " mesh, order=order, dirichlet=\"wall|inlet|circ|circ_inner|outlet\", definedon=\"fluid\"\n", ")\n", "V4 = VectorH1(mesh, order=order, dirichlet=\"circ|circ_inner\", definedon=\"solid\")\n", "Q = L2(mesh, order=0, definedon=\"fluid\", lowest_order_wb=True)\n", "# Lagrange multiplier forcing continuity of velocities across the interface\n", "V5 = SurfaceL2(mesh, order=order, dirichlet=\"wall|outlet|inlet|circ|circ_inner\") ** 2\n", "\n", "X = V2 * V3 * Q * V1 * V4 * V5\n", "Y = V2 * V3 * Q\n", "\n", "(u, uhat, p, d, v, f), (ut, uhatt, pt, dt, vt, ft) = X.TnT()\n", "\n", "gf_solution = GridFunction(X)\n", "gf_solution_old = GridFunction(X)\n", "\n", "gfu, gfvelhat, pressure, gfd, gfv, gff = gf_solution.components\n", "uold, uhatold, pold, dold, vold, fold = gf_solution_old.components\n", "\n", "I = Id(mesh.dim)\n", "F = Grad(d) + I\n", "C = F.trans * F\n", "Fold = Grad(dold) + I\n", "Cold = Fold.trans * Fold\n", "\n", "\n", "n, h = specialcf.normal(2), specialcf.mesh_size\n", "\n", "\n", "def tang(vec):\n", " return vec - (vec * n) * n\n", "\n", "\n", "def norma(vec):\n", " return (vec * n) * n\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\n", "stokes = BilinearForm(Y, symmetric=True, check_unused=False)\n", "stokes += (\n", " nuf * rhof * InnerProduct(2 * Sym(Grad(u)), Sym(Grad(ut)))\n", " - div(u) * pt\n", " - div(ut) * p\n", " - 1e-8 * p * pt\n", ") * dx(\"fluid\")\n", "stokes += (\n", " nuf\n", " * rhof\n", " * (\n", " InnerProduct(2 * Sym(Grad(u)) * n, tang(uhatt - ut))\n", " + InnerProduct(2 * Sym(Grad(ut)) * n, tang(uhat - u))\n", " + alpha * order * order / h * InnerProduct(tang(uhatt - ut), tang(uhat - u))\n", " )\n", ") * dx(\"fluid\", element_boundary=True)\n", "stokes.Assemble()\n", "\n", "true_compile = False\n", "\n", "######################################### SOLID #############################################\n", "bfa = BilinearForm(X, symmetric=False, check_unused=False, condense=True)\n", "bfa += (\n", " (1 / tau * (d - dold) - 0.5 * (v + vold)) * dt\n", " + rhos / tau * (v - vold) * vt\n", " + 2 * InnerProduct(F * Stress(0.5 * (0.5 * (C + Cold) - I)), Grad(vt))\n", ").Compile(true_compile, True) * dx(\"solid\")\n", "\n", "\n", "def minCF(a, b):\n", " return IfPos(a - b, b, a)\n", "\n", "\n", "gfdist = GridFunction(H1(mesh, order=2))\n", "gfdist.Set(minCF((x - 0.6) ** 2 + (y - 0.19) ** 2, (x - 0.6) ** 2 + (y - 0.21) ** 2))\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 / sqrt(gfdist * gfdist + 1e-12) * NeoHookExt(C) * dx(\"fluid\")\n", ").Compile(true_compile, True)\n", "\n", "################################# FLUID ##############################################\n", "graddp = 1 / (tau) * (Grad(d) - Grad(dold))\n", "relv = u - 1 / tau * (d - dold)\n", "relvn = relv * n\n", "upwind = (u * n) * n + (IfPos(relvn, tang(u), tang(uhat)))\n", "\n", "cfinner = (\n", " nuf * rhof * InnerProduct(Sym(Grad(u)) + Sym(Grad(uold)), Sym(Grad(ut)))\n", " - div(u) * pt\n", " - div(ut) * p\n", " - rhof * InnerProduct(0.5 * (u + uold), Grad(ut) * relv)\n", " + rhof * InnerProduct(0.5 * graddp * (u + uold), ut)\n", " + rhof / tau * (u - uold) * ut\n", ")\n", "cfouter = nuf * rhof * (\n", " InnerProduct(Sym(Grad(u) + Grad(uold)) * n, tang(uhatt - ut))\n", " + InnerProduct(2 * Sym(Grad(ut)) * n, 0.5 * tang(uhat - u + uhatold - uold))\n", " + alpha\n", " * order\n", " * order\n", " / h\n", " * InnerProduct(tang(uhatt - ut), 0.5 * tang(uhat - u + uhatold - uold))\n", ") + rhof * (\n", " relvn * upwind * ut\n", " + IfPos(relvn, 1, 0) * 0.5 * relvn * tang(uhat - u + uhatold - uold) * uhatt\n", ")\n", "bfa += cfinner.Compile(true_compile, True) * dx(\"fluid\", deformation=gfd)\n", "bfa += cfouter.Compile(true_compile, True) * dx(\n", " \"fluid\", element_boundary=True, deformation=gfd\n", ")\n", "# interface\n", "bfa += (\n", " (norma(v - u.Trace()) + tang(v - uhat.Trace())) * ft\n", " + (norma(vt - ut.Trace()) + tang(vt - uhatt.Trace())) * f\n", ").Compile(true_compile, True) * ds(\"interface\", deformation=gfd)" ] }, { "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", "inv_stokes = stokes.mat.Inverse(bt_stokes, inverse=\"\")\n", "\n", "gf_stokes = GridFunction(Y)\n", "res_stokes = gf_stokes.vec.CreateVector()\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", "gf_solution.components[0].vec.data = gf_stokes.components[0].vec\n", "gf_solution.components[1].vec.data = gf_stokes.components[1].vec\n", "gf_solution.components[2].vec.data = gf_stokes.components[2].vec\n", "\n", "scene_u = Draw(\n", " gfu,\n", " mesh.Materials(\"fluid\"),\n", " \"velocity\",\n", " deformation=gfd,\n", " order=3,\n", " max=2,\n", ")\n", "scene_p = Draw(pressure, mesh.Materials(\"fluid\"), \"pressure\", deformation=gfd)\n", "scene_d = Draw(gfd, mesh, \"deformation\")\n", "\n", "gf_history = GridFunction(X, multidim=0)\n", "\n", "w = gf_solution.vec.CreateVector()\n", "res = gf_solution.vec.CreateVector()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Time loop." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculate quantities of interest\n", "def CalcForces(disp_x, disp_y):\n", " dmidx, dmidy = gfd(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", "t = 0\n", "i = 0\n", "\n", "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", " gf_solution_old.vec.data = gf_solution.vec\n", "\n", " # Newton\n", " for step in range(2):\n", " bfa.AssembleLinearization(gf_solution.vec)\n", " inv = bfa.mat.Inverse(bfa.space.FreeDofs(bfa.condense), inverse=\"\")\n", " bfa.Apply(gf_solution.vec, res)\n", "\n", " if bfa.condense:\n", " res.data += bfa.harmonic_extension_trans * res\n", " w.data = inv * res\n", " if bfa.condense:\n", " w.data += bfa.harmonic_extension * w\n", " w.data += bfa.inner_solve * res\n", "\n", " gf_solution.vec.data -= w\n", "\n", " if i % 4 == 0:\n", " scene_u.Redraw()\n", " scene_p.Redraw()\n", " scene_d.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", " 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[3],\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": "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 }