19.1. Fluid-structure interaction with Taylor-Hood elements

In this section we use Taylor-Hood elements for the Navier-Stokes equations and Lagrangian elements for the elastic wave equation. We recall from the introduction of FSI that the fluid unknowns are the ALE velocity \(\hat u_f\) and pressure \(p_f\), while the solid unknowns are the displacement \(d_s\) and velocity \(u_s=\dot d_s\). The mesh displacement in the fluid domain is denoted by \(d_m\) with mesh velocity \(w_m=\dot d_m\).

The ALE Navier-Stokes equations read

\[\begin{split} \begin{align*} &\int_{\hat\Omega_f}J_m\Big(\rho_f\frac{\partial \hat{u}_f}{\partial t}\cdot\hat{v} +\rho_f(\nabla_{\hat x}\hat{u}_fF_m^{-1})(\hat{u}_f-w_m)\cdot\hat{v} +2\rho_f\nu_f \mathrm{sym}(\nabla_{\hat x}\hat{u}_fF_m^{-1}):\mathrm{sym}(\nabla_{\hat x}\hat{v}F_m^{-1}) -\mathrm{tr}(\nabla_{\hat x}\hat{v}F_m^{-1})p_f\Big)\,d\hat x = 0&&\qquad \forall \hat{v},\\ &\int_{\hat\Omega_f}J_m\,\mathrm{tr}(\nabla_{\hat x}\hat{u}_fF_m^{-1})q\,d\hat x = 0&&\qquad \forall q. \end{align*} \end{split}\]

The elastic wave equation is written as first order system

\[\begin{split} \begin{align*} & \int_{\Omega_s}\frac{\partial d_s}{\partial t}\cdot v \,dX = \int_{\Omega_s}u_s\cdot v\,dX &&\qquad \forall v,\\ &\int_{\Omega_s}\rho_s\frac{\partial u_s}{\partial t}\cdot w+(F_s\Sigma_s):\nabla w\,dX = 0&&\qquad \forall w. \end{align*} \end{split}\]

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 to decrease with the distance.

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

tau = 0.004
tend = 5
order = 3

velocity_dirichlet = "inlet|wall|circ|circ_inner"
mesh_dirichlet = "inlet|wall|outlet|circ|circ_inner"
fixed_stokes_boundary = "wall|inlet|circ|interface|circ_inner"

rho_s, poisson_s, mu_s, rho_f, nu_f, U = 1e4, 0.4, 0.5 * 1e6, 1e3, 1e-3, 1
lambda_s = 2 * mu_s * poisson_s / (1 - 2 * poisson_s)

# Parabolic inflow profile at the inlet
par = Parameter(0)
u_inflow = par * CoefficientFunction((4 * U * 1.5 * y * (0.41 - y) / (0.41 * 0.41), 0))


# magnitude of inflow velocity over time
def Force(t):
    if t < 0:
        return 0
    elif t < 2:
        return (1 - cos(pi / 2.0 * t)) / 2.0
    else:
        return 1


# directly start with Stokes solution instead of increasing inflow
start_with_stokes = True


def GenerateMesh(order, maxh=0.2):
    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"
    solid.edges.Min(X).name = "circ_inner"

    domain_fluid = (fluid - circle) - solid
    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 use Taylor-Hood elements for the Navier-Stokes equations and \(H^1\)-conforming elements for the elastic wave equation. Since the velocity is continuous across the interface, we can use one global velocity space whose restriction represents \(\hat u_f\) in the fluid and \(u_s\) in the solid. Similarly, one global displacement space represents the solid displacement \(d_s\) and the mesh displacement \(d_m\). With the definedon flag, the pressure space is restricted to the fluid domain.

V = VectorH1(mesh, order=order, dirichlet=velocity_dirichlet)
Q = H1(mesh, order=order - 1, definedon="fluid")
D = VectorH1(mesh, order=order, dirichlet=mesh_dirichlet)

X = V * Q * D
Y = V * Q
(u, p, d), (v, q, w) = X.TnT()

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

velocity, pressure, mesh_deformation = gf_solution.components
velocity_old, pressure_old, mesh_deformation_old = gf_solution_old.components

gradu_old = Grad(velocity_old)
gradd_old = Grad(mesh_deformation_old)

I = Id(mesh.dim)


def CalcStresses(A):
    F = A + I
    C = F.trans * F
    E = 0.5 * (C - I)
    J = Det(F)
    Finv = Inv(F)
    return (F, C, E, J, Finv)


F, C, E, J, Finv = CalcStresses(Grad(d))
F_old, C_old, E_old, J_old, Finv_old = CalcStresses(gradd_old)


def Stress(mat):
    return mu_s * mat + lambda_s / 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, check_unused=False to avoid warning not solving on the solid
stokes = BilinearForm(Y, symmetric=True, check_unused=False)
stokes += (
    nu_f * rho_f * 2 * InnerProduct(Sym(Grad(u)), Sym(Grad(v)))
    - div(u) * q
    - div(v) * p
    - 1e-8 * p * q
) * dx("fluid")
stokes.Assemble()

compile_integrands = False

mesh_velocity = (d - mesh_deformation_old) / tau
grad_u = Grad(u) * Finv
grad_v = Grad(v) * Finv
grad_u_old = gradu_old * Finv_old
grad_v_at_old_geometry = Grad(v) * Finv_old

bfa = BilinearForm(X, symmetric=False, check_unused=False, condense=True)
########################### Fluid: Navier-Stokes ##########################
# M du/dt
bfa += (rho_f / tau * (InnerProduct(0.5 * (J + J_old) * (u - velocity_old), v))).Compile(
    compile_integrands, wait=True
) * dx("fluid")
# symmetric stress div (eps u)
bfa += (
    0.5
    * rho_f
    * nu_f
    * (
        InnerProduct(J * 2 * Sym(grad_u), Sym(grad_v))
        + InnerProduct(J_old * 2 * Sym(grad_u_old), grad_v_at_old_geometry)
    )
).Compile(compile_integrands, wait=True) * dx("fluid")
# Convection and mesh-velocity
bfa += (
    0.5
    * rho_f
    * (
        InnerProduct(J * grad_u * (u - mesh_velocity), v)
        + InnerProduct(
            J_old
            * grad_u_old
            * (velocity_old - mesh_velocity),
            v,
        )
    )
).Compile(compile_integrands, wait=True) * dx("fluid")
# Pressure/Constraint implicit
bfa += (-J * (Trace(grad_v) * p + Trace(grad_u) * q)).Compile(
    compile_integrands, wait=True
) * dx("fluid")

########################### Solid: elastic wave ##########################
# M du/dt
bfa += (rho_s / tau * InnerProduct(u - velocity_old, v)).Compile(
    compile_integrands, wait=True
) * dx("solid")
# Material law
bfa += (InnerProduct(F * Stress(E) + F_old * Stress(E_old), Grad(v))).Compile(
    compile_integrands, wait=True
) * dx("solid")
# dd/dt = u
bfa += (InnerProduct(u + velocity_old - 2.0 * mesh_velocity, w)).Compile(
    compile_integrands, wait=True
) * dx("solid")


########################## Deformation extension ##########################
def minCF(a, b):
    return IfPos(a - b, b, a)


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


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 * mu_s * (1 / sqrt(gf_dist * gf_dist + 1e-12)) * NeoHookExt(C)).Compile(
        compile_integrands, wait=True
    )
    * dx("fluid")
)

Draw(1 / sqrt(gf_dist * gf_dist + 1e-12), mesh, "h(x)", min=0.1, max=100, order=3);

proxy not matching space, checking if it is working anyway

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 corresponding 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(fixed_stokes_boundary))
# set all pressure dofs as active for Stokes
bt_stokes[V.ndof :] = True

gf_stokes = GridFunction(Y)
stokes_velocity, stokes_pressure = gf_stokes.components
res_stokes = gf_stokes.vec.CreateVector()
inv_stokes = stokes.mat.Inverse(bt_stokes, inverse="sparsecholesky")

Reset and draw.

t = 0
i = 0


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


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

if start_with_stokes:
    par.Set(1)
    stokes_velocity.Set(u_inflow, BND, definedon=mesh.Boundaries("inlet"))
    res_stokes.data = stokes.mat * gf_stokes.vec
    gf_stokes.vec.data -= inv_stokes * res_stokes
    velocity.vec.data = stokes_velocity.vec
    pressure.vec.data = stokes_pressure.vec
else:
    gf_solution.vec[:] = 0
gf_solution_old.vec.data = gf_solution.vec

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


gf_history = GridFunction(X, multidim=0)

Time loop.

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

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

        # update inflow by extending inflow profile as Stokes solution
        if t < 2 + tau / 2.0 and not start_with_stokes:
            par.Set(Force(t) - Force(t - tau))
            stokes_velocity.Set(
                u_inflow, BND, definedon=mesh.Boundaries("inlet")
            )
            res_stokes.data = stokes.mat * gf_stokes.vec
            gf_stokes.vec.data -= inv_stokes * res_stokes
            velocity.vec.data += stokes_velocity.vec
            pressure.vec.data += stokes_pressure.vec

        solvers.Newton(bfa, gf_solution, maxit=10, maxerr=1e-2, printing=False)

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

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

        times.append(t)
        CalcForces(disp_x, disp_y)

        gf_solution_old.vec.data = gf_solution.vec

        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[2],
    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()

Using a Taylor-Hood discretization for the Navier-Stokes equations has the drawback that velocity fields are not exactly divergence-free pointwise:

\[ \begin{align*} \int_{\Omega_f}\mathrm{div}(u_f)\,q\,dx = 0 \quad\forall q\quad\nRightarrow \quad\mathrm{div}(u_f)\equiv 0. \end{align*} \]

Using \(H(\mathrm{div})\)-conforming Brezzi-Douglas-Marini or Raviart-Thomas finite elements yields exactly divergence-free solutions. The ALE transformation and interface coupling become a bit more involved: [Neunteufel, Schöberl. Fluid-structure interaction with H(div)-conforming finite elements. Computers and Structures, (2021).], see Notebook FSI with \(H(\mathrm{div})\)-conforming HDG.