14.2. Fluid-structure interaction with H(div)-conforming HDG

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).]

14.2.1. H(div)-conforming HDG for Navier-Stokes with up-winding

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, Notebook Stokes HDG, [Lehrenfeld, Schöberl. High order exactly divergence-free Hybrid Discontinuous Galerkin Methods for unsteady incompressible flows. Computer Methods in Applied Mechanics and Engineering, (2016).].

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

\[\begin{align*} & 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),\\ & M^f_h(u,v)=\int_{\Omega_f} u_T\cdot v_T\,dx,\\ & D^f_h(v,p)=-\int_{\Omega_f} p\,\mathrm{div}(v_T)\,dx, \end{align*}\]

For the convection we use upwinding

\[\begin{align*} 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), \end{align*}\]

with the upwind variable \({u}^{\mathrm{up}}\) defined as

\[\begin{align*} & {u}^{\mathrm{up}}= ({u}_T\cdot {n}){n}+\begin{cases} {u}_T^{{\tau}} & \text{ if }{u}\cdot {n} \geq 0\\ {u}_F^{{\tau}} & \text{ if } {u}\cdot {n} < 0 \end{cases}. \end{align*}\]

14.2.2. Navier Stokes in ALE (Arbitrary Lagrangian Eulerian) formulation

Let \(\Phi(x,t)=\mathrm{id} + d\) 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

\[\begin{align*} u\circ\Phi =P_{\Phi}[\hat{u}]= \frac{1}{J} F \hat{u},\qquad F=\nabla \Phi,\quad J=\det(F). \end{align*}\]

As a result the chain rule now reads

\[\begin{align*} \nabla_xu\circ\Phi = \nabla_{\hat{x}} (P_{\Phi}[\hat{u}])F^{-1}, \end{align*}\]

The Piola transformation also depends on the time. Thus we obtain additional terms from the chain rule of the time derivative

\[\begin{align*} \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} \end{align*}\]

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.

14.2.3. Interface conditions

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

\[\begin{align*} &\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,\\ &\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, \end{align*}\]

which have to be transformed again in the ALE setting.

14.2.4. Benchmark example

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.

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.

from ngsolve import *
from netgen.occ import *
import netgen.geom2d as csg
from ngsolve.webgui import Draw
import ipywidgets as widgets


tau = 0.004
tend = 5
order = 3

# HDG stabilization parameter
alpha = Parameter(5)

# solid density, Poisson ratio, first Lame parameter, fluid density,
# dynamic viscosity, and maximal inflow velocity
rhos, nus, mus, rhof, nuf, U = 1e4, 0.4, 0.5 * 1e6, 1e3, 1e-3, 1
# second Lame parameter for solid
ls = 2 * mus * nus / (1 - 2 * nus)

# inflow profile
u_inflow = CF((4 * U * 1.5 * y * (0.41 - y) / (0.41 * 0.41), 0))


def GenerateMesh(order, maxh=0.203):
    circle = Circle((0.2, 0.2), r=0.05).Face()
    circle.edges.name = "circ"
    fluid = Rectangle(2.5, 0.41).Face()
    fluid.faces.name = "fluid"
    fluid.edges.Min(X).name = "inlet"
    fluid.edges.Max(X).name = "outlet"
    fluid.edges.Min(Y).name = "wall"
    fluid.edges.Max(Y).name = "wall"
    solid = (
        MoveTo(0.248989794855664, 0.19).Rectangle(0.6 - 0.248989794855664, 0.02).Face()
    )
    solid.faces.name = "solid"
    solid.edges.name = "interface"

    domain_fluid = (fluid - circle) - solid
    solid.edges["circ"].name = "circ_inner"
    solid.edges.Min(X).name = "circ_inner"
    domain = Glue([domain_fluid, solid])

    mesh = Mesh(OCCGeometry(domain, dim=2).GenerateMesh(maxh=maxh))
    mesh.Curve(order)

    return mesh


mesh = GenerateMesh(order=order)
Draw(mesh);

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.

V1 = VectorH1(mesh, order=order, dirichlet="wall|outlet|inlet|circ|circ_inner")
V2 = HDiv(
    mesh,
    order=order,
    dirichlet="wall|inlet|circ|circ_inner",
    definedon="fluid",
    hodivfree=True,
)
V3 = TangentialFacetFESpace(
    mesh, order=order, dirichlet="wall|inlet|circ|circ_inner|outlet", definedon="fluid"
)
V4 = VectorH1(mesh, order=order, dirichlet="circ|circ_inner", definedon="solid")
Q = L2(mesh, order=0, definedon="fluid", lowest_order_wb=True)
# Lagrange multiplier forcing continuity of velocities across the interface
V5 = SurfaceL2(mesh, order=order, dirichlet="wall|outlet|inlet|circ|circ_inner") ** 2

X = V2 * V3 * Q * V1 * V4 * V5
Y = V2 * V3 * Q

(u, uhat, p, d, v, f), (ut, uhatt, pt, dt, vt, ft) = X.TnT()

gf_solution = GridFunction(X)
gf_solution_old = GridFunction(X)

gfu, gfvelhat, pressure, gfd, gfv, gff = gf_solution.components
uold, uhatold, pold, dold, vold, fold = gf_solution_old.components

I = Id(mesh.dim)
F = Grad(d) + I
C = F.trans * F
Fold = Grad(dold) + I
Cold = Fold.trans * Fold


n, h = specialcf.normal(2), specialcf.mesh_size


def tang(vec):
    return vec - (vec * n) * n


def norma(vec):
    return (vec * n) * n


def Stress(mat):
    return mus * mat + ls / 2 * Trace(mat) * I

For the time discretization we will use the Crank-Nicolson method

\[\begin{align*} \int_{t^n}^{t^{n+1}} f(s)\,ds \approx \frac{\tau}{2}(f(t^{n+1})+f(t^n)). \end{align*}\]

Only the pressure constraint is handled with implicit Euler.

# For Stokes problem
stokes = BilinearForm(Y, symmetric=True, check_unused=False)
stokes += (
    nuf * rhof * InnerProduct(2 * Sym(Grad(u)), Sym(Grad(ut)))
    - div(u) * pt
    - div(ut) * p
    - 1e-8 * p * pt
) * dx("fluid")
stokes += (
    nuf
    * rhof
    * (
        InnerProduct(2 * Sym(Grad(u)) * n, tang(uhatt - ut))
        + InnerProduct(2 * Sym(Grad(ut)) * n, tang(uhat - u))
        + alpha * order * order / h * InnerProduct(tang(uhatt - ut), tang(uhat - u))
    )
) * dx("fluid", element_boundary=True)
stokes.Assemble()

true_compile = False

######################################### SOLID #############################################
bfa = BilinearForm(X, symmetric=False, check_unused=False, condense=True)
bfa += (
    (1 / tau * (d - dold) - 0.5 * (v + vold)) * dt
    + rhos / tau * (v - vold) * vt
    + 2 * InnerProduct(F * Stress(0.5 * (0.5 * (C + Cold) - I)), Grad(vt))
).Compile(true_compile, True) * dx("solid")


def minCF(a, b):
    return IfPos(a - b, b, a)


gfdist = GridFunction(H1(mesh, order=2))
gfdist.Set(minCF((x - 0.6) ** 2 + (y - 0.19) ** 2, (x - 0.6) ** 2 + (y - 0.21) ** 2))


def NeoHookExt(C, mu=1, lam=1):
    return 0.5 * mu * (Trace(C - I) + 2 * mu / lam * Det(C) ** (-lam / 2 / mu) - 1)


bfa += Variation(
    1e-20 * mus / sqrt(gfdist * gfdist + 1e-12) * NeoHookExt(C) * dx("fluid")
).Compile(true_compile, True)

################################# FLUID ##############################################
graddp = 1 / (tau) * (Grad(d) - Grad(dold))
relv = u - 1 / tau * (d - dold)
relvn = relv * n
upwind = (u * n) * n + (IfPos(relvn, tang(u), tang(uhat)))

cfinner = (
    nuf * rhof * InnerProduct(Sym(Grad(u)) + Sym(Grad(uold)), Sym(Grad(ut)))
    - div(u) * pt
    - div(ut) * p
    - rhof * InnerProduct(0.5 * (u + uold), Grad(ut) * relv)
    + rhof * InnerProduct(0.5 * graddp * (u + uold), ut)
    + rhof / tau * (u - uold) * ut
)
cfouter = nuf * rhof * (
    InnerProduct(Sym(Grad(u) + Grad(uold)) * n, tang(uhatt - ut))
    + InnerProduct(2 * Sym(Grad(ut)) * n, 0.5 * tang(uhat - u + uhatold - uold))
    + alpha
    * order
    * order
    / h
    * InnerProduct(tang(uhatt - ut), 0.5 * tang(uhat - u + uhatold - uold))
) + rhof * (
    relvn * upwind * ut
    + IfPos(relvn, 1, 0) * 0.5 * relvn * tang(uhat - u + uhatold - uold) * uhatt
)
bfa += cfinner.Compile(true_compile, True) * dx("fluid", deformation=gfd)
bfa += cfouter.Compile(true_compile, True) * dx(
    "fluid", element_boundary=True, deformation=gfd
)
# interface
bfa += (
    (norma(v - u.Trace()) + tang(v - uhat.Trace())) * ft
    + (norma(vt - ut.Trace()) + tang(vt - uhatt.Trace())) * f
).Compile(true_compile, True) * ds("interface", deformation=gfd)

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.

bt_stokes = Y.FreeDofs() & ~Y.GetDofs(mesh.Materials("solid"))
bt_stokes &= ~Y.GetDofs(mesh.Boundaries("wall|inlet|circ|interface|circ_inner"))
inv_stokes = stokes.mat.Inverse(bt_stokes, inverse="")

gf_stokes = GridFunction(Y)
res_stokes = gf_stokes.vec.CreateVector()
gf_stokes.components[0].Set(u_inflow, definedon=mesh.Boundaries("inlet"))
res_stokes.data = stokes.mat * gf_stokes.vec
gf_stokes.vec.data -= inv_stokes * res_stokes
gf_solution.components[0].vec.data = gf_stokes.components[0].vec
gf_solution.components[1].vec.data = gf_stokes.components[1].vec
gf_solution.components[2].vec.data = gf_stokes.components[2].vec

scene_u = Draw(
    gfu,
    mesh.Materials("fluid"),
    "velocity",
    deformation=gfd,
    order=3,
    max=2,
)
scene_p = Draw(pressure, mesh.Materials("fluid"), "pressure", deformation=gfd)
scene_d = Draw(gfd, mesh, "deformation")

gf_history = GridFunction(X, multidim=0)

w = gf_solution.vec.CreateVector()
res = gf_solution.vec.CreateVector()

Time loop.

# Calculate quantities of interest
def CalcForces(disp_x, disp_y):
    dmidx, dmidy = gfd(0.6, 0.2)
    disp_x.append(dmidx)
    disp_y.append(dmidy)
    return


disp_x = [0]
disp_y = [0]
times = [0]

t = 0
i = 0

tw = widgets.Text(value="t = 0")
display(tw)

with TaskManager():
    while t < tend - tau / 2.0:
        t += tau

        gf_solution_old.vec.data = gf_solution.vec

        # Newton
        for step in range(2):
            bfa.AssembleLinearization(gf_solution.vec)
            inv = bfa.mat.Inverse(bfa.space.FreeDofs(bfa.condense), inverse="")
            bfa.Apply(gf_solution.vec, res)

            if bfa.condense:
                res.data += bfa.harmonic_extension_trans * res
            w.data = inv * res
            if bfa.condense:
                w.data += bfa.harmonic_extension * w
                w.data += bfa.inner_solve * res

            gf_solution.vec.data -= w

        if i % 4 == 0:
            scene_u.Redraw()
            scene_p.Redraw()
            scene_d.Redraw()

        if i % 20 == 0:
            gf_history.AddMultiDimComponent(gf_solution.vec)

        times.append(t)
        CalcForces(disp_x, disp_y)
        i += 1
        tw.value = f"t = {round(t,5)}"

Draw(
    gf_history.components[0],
    mesh.Materials("fluid"),
    animate=True,
    min=0,
    max=2,
    autoscale=True,
    deformation=gf_history.components[3],
    order=3,
);

Draw results over time. They become periodic after some time.

import matplotlib.pyplot as plt

plt.plot(times, disp_x)
plt.xlabel("time")
plt.ylabel("disp_x")
plt.grid(True)
plt.show()

plt.plot(times, disp_y)
plt.xlabel("time")
plt.ylabel("disp_y")
plt.grid(True)
plt.show()