Higher-order time-stepping schemes for fluid-structure interaction problems

TL;DR

This paper develops and analyzes second-order time-stepping schemes for fluid-structure interaction problems using a distributed Lagrange multiplier approach, demonstrating their stability and effectiveness through numerical tests.

Contribution

It introduces and studies second-order time integration schemes for fluid-structure interaction using a distributed Lagrange multiplier framework, with proven stability and validated numerical results.

Findings

The schemes are unconditionally stable under certain conditions.

Numerical tests confirm the theoretical stability results.

Second-order accuracy is achieved in fluid-structure interaction simulations.

Abstract

We consider a recently introduced formulation for fluid-structure interaction problems which makes use of a distributed Lagrange multiplier in the spirit of the fictitious domain method. In this paper we focus on time integration methods of second order based on backward differentiation formulae and on the Crank-Nicolson method. We show the stability properties of the resulting method; numerical tests confirm the theoretical results.

Figures

Tables

Table 1. Table 1. Mesh parameters

	DOFs ${𝐮}_h$	DOFs $p_h$	DOFs ${𝐗}_h$	DOFs ${\lambda }_h$
coarse mesh (M = $8$)	$578$	$209$	$306$	$306$
fine mesh (M = $16$)	$2,178$	$801$	$1,122$	$1,122$

Table 2. Table 2. Convergence results for the fully implicit scheme on the coarse mesh

Fluid velocity
	BDF1		BDF2		CNm		CNt
$\Delta t$	$L^2$ error	rate	$L^2$ error	rate	$L^2$ error	rate	$L^2$ error	rate
$0.05$	$7.63\cdot {10}^{-2}$		$2.97\cdot {10}^{-2}$		$2.37\cdot {10}^{-1}$		$2.42\cdot {10}^{-1}$
$0.025$	$4.11\cdot {10}^{-2}$	$0.89$	$4.90\cdot {10}^{-3}$	$2.60$	$6.24\cdot {10}^{-2}$	$1.92$	$6.02\cdot {10}^{-2}$	$2.00$
$0.0125$	$2.13\cdot {10}^{-2}$	$0.95$	$1.13\cdot {10}^{-3}$	$2.11$	$1.21\cdot {10}^{-2}$	$2.36$	$1.10\cdot {10}^{-2}$	$2.45$
$0.00625$	$1.08\cdot {10}^{-2}$	$0.97$	$2.86\cdot {10}^{-4}$	$1.98$	$2.03\cdot {10}^{-3}$	$2.58$	$9.95\cdot {10}^{-4}$	$3.47$

Table 3. Table 3. Convergence results for the fully implicit scheme on the fine mesh

Fluid velocity
	BDF1		BDF2		CNm		CNt
$\Delta t$	$L^2$ error	rate	$L^2$ error	rate	$L^2$ error	rate	$L^2$ error	rate
$0.05$	$9.05\cdot {10}^{-2}$		$3.62\cdot {10}^{-2}$		$2.28\cdot {10}^{-1}$		$2.26\cdot {10}^{-1}$
$0.025$	$4.87\cdot {10}^{-2}$	$0.89$	$5.05\cdot {10}^{-3}$	$2.84$	$6.23\cdot {10}^{-2}$	$1.87$	$6.04\cdot {10}^{-2}$	$1.91$
$0.0125$	$2.54\cdot {10}^{-2}$	$0.94$	$1.20\cdot {10}^{-3}$	$2.07$	$2.28\cdot {10}^{-2}$	$1.45$	$2.07\cdot {10}^{-2}$	$1.54$
$0.00625$	$1.29\cdot {10}^{-2}$	$0.98$	$3.53\cdot {10}^{-4}$	$1.77$	$5.27\cdot {10}^{-3}$	$2.11$	$4.03\cdot {10}^{-3}$	$2.36$

Table 4. Table 4. Maximum iterates of the nonlinear solver on the coarse mesh

$\Delta t$	BDF1	BDF2	CNm	CNt
$0.05$	$5$	$5$	$4$	$7$
$0.025$	$4$	$4$	$3$	$4$
$0.0125$	$3$	$3$	$3$	$3$
$0.00625$	$3$	$3$	$2$	$3$

Table 5. Table 5. Maximum iterates of the nonlinear solver on the fine mesh

$\Delta t$	BDF1	BDF2	CNm	CNt
$0.05$	$10$	$5$	$6$	$6$
$0.025$	$6$	$5$	$5$	$4$
$0.0125$	$6$	$4$	$4$	$4$
$0.00625$	$4$	$4$	$3$	$3$

Table 6. Table 6. Convergence results for the semi-implicit scheme on the coarse mesh

Fluid velocity
	BDF1		BDF2		CNm		CNt
$\Delta t$	$L^2$ error	rate	$L^2$ error	rate	$L^2$ error	rate	$L^2$ error	rate
$0.05$	$7.78\cdot {10}^{-2}$		$3.05\cdot {10}^{-2}$		$2.49\cdot {10}^{-1}$		$2.58\cdot {10}^{-1}$
$0.025$	$4.17\cdot {10}^{-2}$	$0.90$	$7.89\cdot {10}^{-3}$	$1.95$	$6.24\cdot {10}^{-2}$	$2.00$	$6.74\cdot {10}^{-2}$	$1.94$
$0.0125$	$2.17\cdot {10}^{-2}$	$0.95$	$3.14\cdot {10}^{-3}$	$1.33$	$1.24\cdot {10}^{-2}$	$2.33$	$2.64\cdot {10}^{-2}$	$1.35$
$0.00625$	$1.10\cdot {10}^{-2}$	$0.97$	$1.29\cdot {10}^{-3}$	$1.29$	$2.25\cdot {10}^{-3}$	$2.47$	$3.18\cdot {10}^{-3}$	$3.06$

Table 7. Table 7. Convergence results for the semi-implicit scheme on the fine mesh

Fluid velocity
	BDF1		BDF2		CNm		CNt
$\Delta t$	$L^2$ error	rate	$L^2$ error	rate	$L^2$ error	rate	$L^2$ error	rate
$0.05$	$9.18\cdot {10}^{-2}$		$3.89\cdot {10}^{-2}$		$2.36\cdot {10}^{-1}$		$2.39\cdot {10}^{-1}$
$0.025$	$5.05\cdot {10}^{-2}$	$0.86$	$8.59\cdot {10}^{-3}$	$2.18$	$7.54\cdot {10}^{-2}$	$1.64$	$7.06\cdot {10}^{-2}$	$1.76$
$0.0125$	$2.63\cdot {10}^{-2}$	$0.94$	$3.32\cdot {10}^{-3}$	$1.37$	$4.24\cdot {10}^{-2}$	$0.83$	$2.22\cdot {10}^{-2}$	$1.67$
$0.00625$	$1.33\cdot {10}^{-2}$	$0.98$	$1.40\cdot {10}^{-3}$	$1.24$	$2.19\cdot {10}^{-2}$	$0.96$	$4.19\cdot {10}^{-3}$	$2.40$

Table 8. Table 8. Maximum residual in the semi-implicit scheme on the coarse mesh

$\Delta t$	BDF1	BDF2	CNm	CNt
$0.05$	$3.83\cdot {10}^{-3}$	$3.83\cdot {10}^{-3}$	$9.87\cdot {10}^{-3}$	$9.64\cdot {10}^{-2}$
$0.025$	$2.09\cdot {10}^{-3}$	$2.30\cdot {10}^{-3}$	$1.24\cdot {10}^{-3}$	$2.17\cdot {10}^{-2}$
$0.0125$	$7.41\cdot {10}^{-4}$	$8.26\cdot {10}^{-4}$	$3.62\cdot {10}^{-4}$	$7.90\cdot {10}^{-3}$
$0.00625$	$2.23\cdot {10}^{-4}$	$2.45\cdot {10}^{-4}$	$1.08\cdot {10}^{-4}$	$9.55\cdot {10}^{-4}$

Table 9. Table 9. Maximum residual in the semi-implicit scheme on the fine mesh

$\Delta t$	BDF1	BDF2	CNm	CNt
$0.05$	$3.59\cdot {10}^{-2}$	$3.59\cdot {10}^{-2}$	$3.43\cdot {10}^{-2}$	$3.81\cdot {10}^{-2}$
$0.025$	$1.03\cdot {10}^{-2}$	$1.07\cdot {10}^{-2}$	$8.19\cdot {10}^{-3}$	$1.02\cdot {10}^{-2}$
$0.0125$	$5.28\cdot {10}^{-3}$	$5.87\cdot {10}^{-3}$	$1.54\cdot {10}^{-3}$	$2.27\cdot {10}^{-3}$
$0.00625$	$1.46\cdot {10}^{-3}$	$1.44\cdot {10}^{-3}$	$4.33\cdot {10}^{-4}$	$5.08\cdot {10}^{-4}$

Equations

$${\mathbf{\sigma }}_f=-p_f\mathbf{I}+{\nu }_f{\nabla }_{\mathrm{sym}}{\mathbf{u}}_f,$$

$${\dot{\mathbf{u}}}_f=\frac{\partial {\mathbf{u}}_f}{\partial t}+{\mathbf{u}}_f\cdot \nabla {\mathbf{u}}_f.$$

$$\mathbf{u}_{s}(\mathbf{x},t)=\frac{\partial\mathbf{X}(\mathbf{s},t)}{\partial t}\Big{|}_{\mathbf{x}=\mathbf{X}(\mathbf{s},t)}$$

$${\mathbf{\sigma }}_s^v=-p_s\mathbb{I}+{\nu }_s{\nabla }_{\mathrm{sym}}{\mathbf{u}}_s,$$

$$\mathbb{P}(\mathbb{F}(\mathbf{s},t))=|\mathbb{F}(\mathbf{s},t)|{\mathbf{\sigma }}_s^e(\mathbf{x},t){\mathbb{F}}^{-\top }(\mathbf{s},t)\text{for }\mathbf{x}=\mathbf{X}(\mathbf{s},t).$$

$$(\mathbb{P}(\mathbb{F}(\mathbf{s},t){)}_{\alpha i}=\frac{\partial W}{\partial {\mathbb{F}}_{\alpha i}}(\mathbb{F}(\mathbf{s},t))={\left(\frac{\partial W}{\partial \mathbb{F}}(\mathbb{F}(\mathbf{s},t))\right)}_{\alpha i},$$

$$E\left(\mathbf{X}(t)\right)={\int }_{\mathcal{B}}W(\mathbb{F}(s,t))\mathrm{d}\mathbf{s}.$$

$${\rho }_f{\dot{\mathbf{u}}}_f=\mathrm{div}{\mathbf{\sigma }}_f$$

$$\mathrm{div}{\mathbf{u}}_f=0$$

$$\displaystyle\rho_{s}\frac{\partial^{2}\mathbf{X}}{\partial t^{2}}={\operatorname*{div}}_{s}(|\mathbb{F}|\bm{\sigma}_{s}^{v}\mathbb{F}^{-\top}+\mathbb{P}(\mathbb{F}))\$$

$$\mathrm{div}{\mathbf{u}}_s=0$$

$${\mathbf{u}}_f=\frac{\partial \mathbf{X}}{\partial t}$$

$${\mathbf{\sigma }}_f{\mathbf{n}}_f=-({\mathbf{\sigma }}_s^f+|\mathbb{F}{|}^{-1}\mathbb{P}{\mathbb{F}}^{\top }){\mathbf{n}}_s$$

$${\mathbf{u}}_f(0)={\mathbf{u}}_{f0}$$

$${\mathbf{u}}_s(0)={\mathbf{u}}_{s0}$$

$$\mathbf{X}(0)={\mathbf{X}}_0$$

$${\mathbf{u}}_f(t)=0$$

$$\mathbf{u}=\left\{\begin{array}[]{ll}\mathbf{u}_{f}&\text{ in }\Omega^{f}_{t}\\ \mathbf{u}_{s}&\text{ in }\Omega^{s}_{t}\end{array}\right.\quad p=\left\{\begin{array}[]{ll}p_{f}&\text{ in }\Omega^{f}_{t}\\ p_{s}&\text{ in }\Omega^{s}_{t}.\end{array}\right.$$

$${\rho }_f\dot{\mathbf{u}}=\mathrm{div}(\nu {\nabla }_{\mathrm{sym}}\mathbf{u})-\nabla p$$

$$\mathrm{div}\mathbf{u}=0$$

$$\nu=\left\{\begin{array}[]{ll}\nu_{f}&\text{ in }\Omega^{f}_{t}\\ \nu_{s}&\text{ in }\Omega^{s}_{t}.\end{array}\right.$$

$${\rho }_f\frac{{\partial }^2\mathbf{X}}{\partial t^2}={\mathrm{div}}_s(|\mathbb{F}|\mathbf{\sigma })s^v{\mathbb{F}}^{-\top }\text{in }\mathcal{B}$$

$$({\rho }_s-{\rho }_f)\frac{{\partial }^2\mathbf{X}}{\partial t^2}={\mathrm{div}}_s\mathbb{P}(\mathbb{F})\text{in }\mathcal{B}$$

$${\rho }_f\dot{\mathbf{u}}=\mathrm{div}(\nu {\nabla }_{\mathrm{sym}}\mathbf{u})-\nabla p$$

$$\mathrm{div}\mathbf{u}=0$$

$$({\rho }_s-{\rho }_f)\frac{{\partial }^2\mathbf{X}}{\partial t^2}={\mathrm{div}}_s\mathbb{P}(\mathbb{F})$$

$$\mathbf{u}(\mathbf{X}(\mathbf{s},t),t)=\frac{\partial \mathbf{X}(\mathbf{s},t)}{\partial t}$$

$${\mathbf{\sigma }}_f{\mathbf{n}}_f=-({\mathbf{\sigma }}_s^v+|\mathbb{F}{|}^{-1}\mathbb{P}{\mathbb{F}}^{\top }){\mathbf{n}}_s$$

$$\mathbf{u}(0)={\mathbf{u}}_0$$

$$\mathbf{X}(0)={\mathbf{X}}_0$$

$$\mathbf{u}(t)=0$$

$${\int }_{\Omega }{\rho }_f\dot{\mathbf{u}}\mathbf{v}\mathrm{d}\mathbf{x}$$

$$-{\int }_{{\Gamma }_t}p\mathbf{v}\cdot {\mathbf{n}}_f\mathbf{d}\gamma +{\int }_{{\Gamma }_t}{\nu }_s{\nabla }_{\mathrm{sym}}\mathbf{u}{\mathbf{n}}_s\mathbf{vd}\gamma -{\int }_{{\Gamma }_t}p\mathbf{v}\cdot {\mathbf{n}}_s\mathbf{d}\gamma$$

$$=-{\int }_{\Omega }\nu {\nabla }_{\mathrm{sym}}\mathbf{u}:{\nabla }_{\mathrm{sym}}\mathbf{v}\mathrm{d}\mathbf{x}+{\int }_{\Omega }p\mathrm{div}\mathbf{v}\mathrm{d}\mathbf{x}+{\int }_{{\Gamma }_t}|\mathbb{F}{|}^{-1}\mathbb{P}(\mathbb{F}){\mathbb{F}}^{\top }{\mathbf{n}}_s\mathbf{vd}\gamma .$$

$$({\rho }_s-{\rho }_f){\int }_{\mathcal{B}}\frac{{\partial }^2\mathbf{X}}{\partial t^2}\mathbf{z}\mathrm{d}\mathbf{s}=-{\int }_{\mathcal{B}}\mathbb{P}(\mathbb{F}):{\nabla }_{\mathrm{s}}\mathbf{z}\mathrm{d}\mathbf{s}+{\int }_{\partial \mathcal{B}}\mathbb{P}(\mathbb{F})\mathbf{Nzd}\gamma ,$$

$$\mathbf{c}(\mu ,\mathbf{z})=0\forall \mu \in \mathbf{\Lambda }\text{ implies }\mathbf{z}=0.$$

$$\mathbf{c}(\lambda ,\mathbf{z})={\int }_{\partial \mathcal{B}}\mathbb{P}(\mathbb{F})\mathbf{Nzd}\gamma \forall \mathbf{z}\in H^1(\mathcal{B}{)}^d.$$

$$\mathbf{c}\left(\mu ,\mathbf{u}(\mathbf{X}(\mathbf{s},t),t)-\frac{\partial (\mathbf{X}(\mathbf{s},t))}{\partial t}\right)=0$$

$${\rho }_f{\left(\frac{\partial }{\partial t}\mathbf{u}(t),\mathbf{v}\right)}_{\Omega }+b\left(\mathbf{u}(t),\mathbf{u}(t),\mathbf{v}\right)+a\left(\mathbf{u}(t),\mathbf{v}\right)$$

Full text

∗ Corresponding author: Daniele Boffi

Daniele Boffi∗

Dipartimento di Matematica “F. Casorati”, University of Pavia

Pavia, Italy

Lucia Gastaldi

DICATAM, University of Brescia

Brescia, Italy

Sebastian Wolf

Technische Universität München (TUM)

München, Germany

1. Introduction

We discuss a scheme involving a fictitious domain approach with a distributed Lagrange multiplier for the modeling of fluid-structure interaction problems [3, 5], which has originated as a natural evolution of the Finite Element Immersed Boundary Method introduced and studied in [6, 8, 2, 4]. This investigation started from the framework of Peskin’s research [19] who introduced the (finite difference) Immersed Boundary Method for the modeling of fluid-structure interaction problems.

The project described in this paper has been carried on during the Master thesis of the third author who spent a semester in Pavia within an exchange program between TUM and Pavia. The aim of the project was to investigate and analyze higher order time schemes for our fluid-structure interaction numerical approach which, so far, had been presented only in combination with low order Euler schemes. On the other hand, in the case of thick (codimension zero) solids, the regularity of the solution allows for a convergence in space higher than first order, so that it may pay off to make use of higher order schemes.

The main result of this paper consists in the implementation and in the stability analysis for the second order Backward Differentiation Formula BDF2 and for the Crank–Nicolson scheme. We present two possible variants of the Crank–Nicolson scheme, which differ in the treatment of the nonlinear terms.

The structure of our paper is as follows: after introducing the model and its finite element discretization in Section 2, we describe different time stepping schemes in Section 3: Backward Euler BDF1, BDF2, Crank–Nicolson (version based on midpoint rule CNm or based on trapezoidal rule CNt). Finally, in Section 4, we present several numerical experiments confirming the convergence and the stability proved in the previous section.

2. Setting of the FSI problem

We consider a fluid-structure interaction problem consisting of a visco-elastic solid immersed in a fluid. The solid is initially distorted from its equilibrium configuration, so that it tends to return to its equilibrium position. In doing so, the region occupied by the fluid changes its shape, thus inducing a flow which in turn produces a force on the solid, which deforms accordingly. We assume that both the fluid and the solid are incompressible. An extension to our model to compressible solids has been studied in [7] but is not going to be considered in this paper.

Let $\Omega\subset\mathbb{R}^{d}$, with $d=2,3$, be a connected, open, and bounded domain with Lipschitz continuous boundary $\partial\Omega$. For simplicity, we assume that $\Omega$ is a polyhedron. The domain $\Omega$ is split into two non intersecting open domains $\Omega^{f}_{t}$ and $\Omega^{s}_{t}$ which represent the regions occupied at time $t$ by fluid and solid, respectively. Hence we have $\overline{\Omega}=\overline{\Omega^{f}_{t}}\cup\overline{\Omega^{s}_{t}}$. We denote by $\Gamma_{t}$ the interface between $\Omega^{f}_{t}$ and $\Omega^{s}_{t}$ and assume that it has empty intersection with the exterior boundary $\partial\Omega$. Let $\mathcal{B}$ be the reference domain of $\Omega^{s}_{t}$, and let $\mathbf{X}:\mathcal{B}\to\Omega^{s}_{t}$ represent the corresponding deformation mapping. Hence a point $x\in\Omega^{s}_{t}$ is the image at time $t$ of a point $\mathbf{s}\in\mathcal{B}$, that is $\mathbf{x}=\mathbf{X}(\mathbf{s},t)$. For simplicity, we assume that $\mathcal{B}=\Omega^{s}_{0}$ is the initial position of the solid. We denote by $\mathbb{F}=\operatorname*{\bm{\nabla}}_{s}\mathbf{X}$ the deformation gradient and by $J=\det(\mathbb{F})$ its Jacobian.

We are going to use the following notation: $\mathbf{u}_{f}$, $p_{f}$, $\bm{\sigma}_{f}$, and $\rho_{f}$ denote, respectively, velocity, pressure, stress tensor and density in the fluid. We consider a Newtonian fluid characterized by the usual Navier–Stokes stress tensor

[TABLE]

where $\nu_{f}$ is the fluid viscosity and $\operatorname*{\bm{\nabla}_{\rm sym}}\mathbf{u}=(1/2)\left(\operatorname*{\bm{\nabla}}\mathbf{u}_{f}+(\operatorname*{\bm{\nabla}}\mathbf{u}_{f})^{\top}\right)$. In the fluid, we use an Eulerian description so that the material derivative is given by

[TABLE]

In the solid, $\mathbf{u}_{s}$, $p_{s}$, and $\rho_{s}$ stand, respectively, for velocity, pressure, and density. In the solid the Lagrangian framework is preferred, and the spatial description of the material velocity reads

[TABLE]

so that $\dot{\mathbf{u}}_{s}=\partial^{2}\mathbf{X}/\partial t^{2}$. Moreover, we assume that the solid material is viscous-hyperelastic, so that the Cauchy stress tensor is given by the sum of a viscous part

[TABLE]

where $\nu_{s}$ is the viscosity, and an elastic part $\bm{\sigma}_{s}^{e}$, which can be expressed in terms of the Piola–Kirchhoff stress tensor $\mathbb{P}$:

[TABLE]

The Piola–Kirchhoff stress tensor is related to the positive energy density $W(\mathbb{F})$, which characterizes hyperelastic materials, as follows:

[TABLE]

where $i=1,\ldots,m$ and $\alpha=1,\ldots,d$. The elastic potential energy of the body is given by:

[TABLE]

Assuming that both the fluid and the solid material are incompressible, we have the following mathematical model for the fluid-structure system:

[TABLE]

The last two equations in (7) represent the transmission condition at the interface $\Gamma_{t}$. The model is completed with initial and boundary conditions:

[TABLE]

Following [5], we apply a fictitious domain approach by extending the first equation in (7) to the whole domain $\Omega$ and by using the following new unknowns:

[TABLE]

The extended velocity and pressure satisfy the following equation all over the domain:

[TABLE]

where we set

[TABLE]

The first equation in (10) enforces that $\rho_{f}\dot{\mathbf{u}}_{s}=\operatorname*{div}\bm{\sigma}_{s}^{v}$ in $\mathcal{B}_{t}$. By taking into account (2) and changing variable, this is equivalent to

[TABLE]

which subtracted from the third equation in (7) gives

[TABLE]

and the model problem can be rewritten as follows: find $u$, $p$, $\mathbf{X}$ such that

[TABLE]

We observe that the second equation in (12) enforces the divergence free constraint both for fluid and solid. In order to arrive to our weak formulation we multiply the first equation by a test function $v\in H^{1}_{0}(\Omega)^{d}$ and integrate by parts the right hand side taking into account (9) and the fact that the normal derivatives of $\mathbf{u}$ might jump across $\Gamma_{t}$, hence we have

[TABLE]

On the other hand, multiplying by a test function $\mathbf{z}\in H^{1}(\mathcal{B})^{d}$ the third equation in (12) and integrating by parts, we obtain

[TABLE]

where $\mathbf{N}$ is the outward normal unit vector to $\partial\mathcal{B}$. By change of variables, the integral on the boundary is equal to the last integral on the right hand side of (13) if we choose $\mathbf{v}(\mathbf{X}(\mathbf{s},t))=\mathbf{z}(s)$ for $\mathbf{s}\in\partial\mathcal{B}$.

Let $\bm{\Lambda}$ be a functional space to be defined later on and $\mathbf{c}:\bm{\Lambda}\times H^{1}(\mathcal{B})^{d}\to\mathbb{R}$ a continuous bilinear form such that

[TABLE]

Possible definitions for $\bm{\Lambda}$ and $\mathbf{c}$ are:

• •

$\mathbf{c}$ is the duality pairing between $H^{1}(\mathcal{B})^{d}$ its dual, that is $\bm{\Lambda}=(H^{1}(\mathcal{B})^{d})^{\prime}$ and $\mathbf{c}(\mu,\mathbf{z})=\langle\mu,\mathbf{z}\rangle_{\mathcal{B}}$ for $\mu\in\bm{\Lambda}$, $\mathbf{z}\in H^{1}(\mathcal{B})^{d}$;

• •

$\mathbf{c}$ is the scalar product in $H^{1}(\mathcal{B})^{d}$, that is $\bm{\Lambda}=H^{1}(\mathcal{B})^{d}$ and $\mathbf{c}(\mu,\mathbf{z})=(\mu,\mathbf{z})_{\mathcal{B}}+(\operatorname{\bm{\nabla}_{s}}\mu,\operatorname{\bm{\nabla}_{s}}\mathbf{z})_{\mathcal{B}}$ for $\mu,\mathbf{z}\in H^{1}(\mathcal{B})^{d}$.

Then we introduce $\lambda\in\bm{\Lambda}$ satisfying the following relation:

[TABLE]

Using the bilinear form $\mathbf{c}$ we can formulate the kinematic constraint in the fourth equation in (12) as

[TABLE]

and our problem can be rewritten in the following weak form (see, also, [3, 5]).

Problem 1**.**

For given $\mathbf{u}_{0}\in H^{1}_{0}(\Omega)^{d}$ and $\mathbf{X}_{0}\in W^{1,\infty}(\mathcal{B})$, find $\mathbf{u}(t)\in H^{1}_{0}(\Omega)^{d}$, $p(t)\in L^{2}_{0}(\Omega)$, $\mathbf{X}(t)\in H^{1}(\mathcal{B})^{d}$, and $\lambda(t)\in\bm{\Lambda}$ such that for almost all $t\in(0,T)$:

[TABLE]

In the above problem we have used the following notation: $\delta\rho=\rho_{s}-\rho_{f}$, the scalar product in $L^{2}(D)$ is denoted by $(\cdot,\cdot)_{D}$, and

[TABLE]

Remark*.*

We remark that the unknown $\lambda$ in Problem 1 plays the role of a Lagrange multiplier associated with the condition which enforces the kinematic constraint, that is the equality of the velocity $\mathbf{u}$ with the solid velocity in the region occupied by the structures, see also (2).

Remark*.*

We can introduce an alternative to the third equation in (16) by splitting it into a system of two equations of first order in time, and by introducing a new unknown $\dot{\mathbf{X}}(t)=\frac{\partial\mathbf{X}(t)}{\partial t}$:

[TABLE]

This formulation is more suited when a second order time marching scheme is used, since we do not need to introduce a second order approximation of the second time derivative.

By choosing properly the test functions in (16), one can obtain the following energy estimate [5] :

Proposition 2.1*.*

Let us assume that $\delta\rho\geq 0$, that the potential energy density $W$ is a $C^{1}$ convex function over the set of second order tensors, and that for almost every $t\in[0,T]$ the solution of Problem 1 is such that $\mathbf{X}(t)\in(W^{1,\infty}(\mathcal{B}))^{d}$ with $\frac{\partial\mathbf{X}}{\partial t}(t)\in L^{2}(\mathcal{B})^{d}$, then the following equality holds true

[TABLE]

In the above proposition $\|\cdot\|_{D}$ stands for the norm in $L^{2}(D)$. We observe that the above proposition still holds true when the alternative formulation in (18) is used.

2.1. Finite element discretization

Let $\mathcal{T}_{h}$ and $\mathcal{T}_{h}^{\mathcal{B}}$ be regular meshes in $\Omega$ and $\mathcal{B}$, respectively, which are independent one from each other. The corresponding mesh sizes will be denoted by $h_{f}$ and $h_{s}$, respectively. We consider two finite element spaces $\mathbf{V}_{h}\subset H^{1}_{0}(\Omega)^{d}$ and $Q_{h}\subset L^{2}_{0}(\Omega)$ such that the pair $(\mathbf{V}_{h},Q_{h})$ satisfies the usual inf-sup condition for the Stokes equations. We assume that $\mathcal{T}_{h}^{\mathcal{B}}$ contains only simplices and introduce the space of continuous piecewise affine functions on $\mathcal{T}_{h}^{\mathcal{B}}$

[TABLE]

In order to discretize $\bm{\Lambda}$, we set $\bm{\Lambda}_{h}=\mathbf{S}_{h}$. With this definition we have that when $\bm{\Lambda}=(H^{1}(\mathcal{B})^{d})^{\prime}$ and $\mathbf{c}$ is the duality pairing, we can compute easily $\mathbf{c}$ using the scalar product in $L^{2}(\mathcal{B})$.

Then the discrete counterpart of Problem 1 reads as follows.

Problem 2**.**

For given $\mathbf{u}_{0h}\in\mathbf{V}_{h}$ and $\mathbf{X}_{0h}\in\mathbf{S}_{h}$, find $\mathbf{u}_{h}(t)\in\mathbf{V}_{h}$, $p_{h}(t)\in Q_{h}$, $\mathbf{X}_{h}(t)\in\mathbf{S}_{h}$, $\lambda_{h}(t)\in\bm{\Lambda}_{h}$ such that for almost all $t\in(0,T)$:

[TABLE]

This semi discrete problem inherits the same energy estimate as the continuous one, which can be proved with the same technique.

The system (18) can be discretized as follows: find $\dot{\mathbf{X}}_{h}(t)\in\mathbf{S}_{h}$ and $\mathbf{X}_{h}(t)\in\mathbf{S}_{h}$ such that

[TABLE]

We observe that the first equation can be solved exactly. In the following we shall use these two equations instead of the third equation in (21) when we apply higher order time marching schemes.

3. Higher-order time-stepping method

The system (21) is a system of ODEs with some algebraic constraints in $\mathbb{R}^{n}$, where $n=\dim(\mathbf{V}_{h})+\dim(Q_{h})+\dim(\mathbf{S}_{h})+\dim(\bm{\Lambda}_{h})$. The numerical solution of Problem 2 can be computed by using some ODE/DAE solver. In order to avoid to use excessively small time steps, in [3, 5] the Backward Euler formula was applied for the time-integration, and the unconditional stability of the scheme has been proved. The aim of this paper is to introduce methods which achieve second order convergence in time and are unconditionally stable. In particular, we shall consider two one-step methods - based on the midpoint and trapezoidal rules - and the BDF2 method and analyze their stability. We observe that the resulting fully discrete scheme is nonlinear, and we shall discuss how the solution can be obtained.

3.1. Backward Euler

Before entering into the details of higher order methods, we recall the Backward Euler method analyzed in [3]. We subdivide the time interval $[0,T]$ into $N$ equal parts with size $\Delta t=T/N$ and subdivision points $t_{n}=n\Delta t$. Moreover, for a certain function $y(t)$, we set $y^{n}=y(t_{n})$ and use the following finite difference in order to approximate the time derivatives:

[TABLE]

Notice that both approximations are of first order.

The fully discrete version of Problem 1 using the Backward Euler scheme is the following one.

Problem 3**.**

Given $\mathbf{u}_{0h}\in\mathbf{V}_{h}$ and $\mathbf{X}_{0h}\in\mathbf{S}_{h}$, for all $n=1,\dots,N$ find $\mathbf{u}_{h}^{n}\in\mathbf{V}_{h}$, $p_{h}^{n}\in Q_{h}$, $\mathbf{X}_{h}^{n}\in\mathbf{S}_{h}$, and $\lambda_{h}^{n}\in\bm{\Lambda}_{h}$ fulfilling:

[TABLE]

We recall the stability estimate proven in [3].

Proposition 3.1*.*

Let the material behavior be governed by an energy density $W$ which is $\mathcal{C}^{1}$ and convex. Let $\mathbf{u}_{h}^{n}\in\mathbf{V}_{h}$ and $\mathbf{X}_{h}^{n}\in\mathbf{S}_{h}$, $n=1,\dots N$ be solutions of (24). Then the following estimate holds true:

[TABLE]

This system is highly nonlinear and hence not immediate to solve. It can be solved by a fixed-point method or by replacing some quantities at time $t_{n+1}$ with their known value at time $t_{n}$. For example, in [3] the latter approach has been preferred; the value of $\mathbf{v}(\mathbf{X}_{h}^{n+1})$ in the first equation has been replaced by $\mathbf{v}(\mathbf{X}^{n}_{h})$ and, similarly, $\mathbf{u}^{n+1}_{h}(\mathbf{X}^{n+1}_{h})$ by $\mathbf{u}^{n+1}_{h}(\mathbf{X}^{n}_{h})$. Moreover, as it is usual in the discretization of the Navier–Stokes equations, the first argument in the trilinear form $b$ has been evaluated at time $t_{n}$. Notice that the same stability estimate as for the fully implicit scheme has been proved. We observe also that this modification introduces an approximation of first order, which should not affect the accuracy of the Backward Euler scheme, since it is of first order, too.

At the first step $n=0$, the approximation of the second time derivative requires the knowledge of $\mathbf{X}_{h}^{-1}$. Following [5], we can compute this value using the kinematic constraint and the initial data

[TABLE]

3.2. BDF2 method

Backwards differentiation formulae (BDF) are popular multistep methods to solve stiff ODE problems. They can be derived by approximating the time derivative at time $t_{n+1}$ with finite differences with order greater than or equal to 1. In particular, using a first order finite difference we have the Backward Euler method. In this paper we focus on the second order formula BDF2, hence, for a certain function $y(t)$, we approximate the time derivative with

[TABLE]

When applied to ODEs, it is well-known that BDF2 is convergent of order 2 and A-stable, (see for more detail [10, Chapter 7]). BDF methods have been successfully applied to the Navier–Stokes equations [16], to nonlinear structural mechanics [11], to electro-magnetic problems [18], and to Stokes–Darcy flow [9].

The application of the BDF2 method to Problem 2 gives the following formulation.

Problem 4**.**

Given $\mathbf{u}_{0h}\in\mathbf{V}_{h}$, $\mathbf{X}_{0h}\in\mathbf{S}_{h}$, for all $n=1,\dots,N$ find $\mathbf{u}_{h}^{n}\in\mathbf{V}_{h}$, $p_{h}^{n}\in Q_{h}$, $\mathbf{X}_{h}^{n}\in\mathbf{S}_{h}$, $\dot{\mathbf{X}}_{h}^{n}\in\mathbf{S}_{h}$ and $\lambda_{h}^{n}\in\bm{\Lambda}_{h}$ fulfilling:

[TABLE]

To start the computations, we additionally need the quantities $\mathbf{u}_{h}^{1}$ and $\mathbf{X}_{h}^{1}$ which are usually obtained by a one-step method. To achieve second order consistency theoretically, we need to apply a method of second order to the first step. This can be done, for instance, by using the Crank–Nicolson scheme, which is analyzed in this paper as well. In our numerical experiments we also tried a start-up with the backward Euler method, which also produced second order convergence.

In the following proposition, we prove an energy estimate for the scheme described by Equations (26a)-(26f).

Proposition 3.2*.*

Assume that $\delta\rho\geq 0$ and that the solids behavior is linear: $\mathbb{P}(\mathbb{F})=\kappa\mathbb{F}$. Let $\mathbf{u}_{h}^{n}$, $\mathbf{X}_{h}^{n}$, $n=1,\dots N$, be solutions of Problem 4, then the following estimate holds true:

[TABLE]

Proof.

We mimic the proof of [3, Prop. 3] and use some useful tricks presented in [16]. We test (26a) with $\mathbf{u}_{h}^{n+1}$:

[TABLE]

We remind that the trilinear form $b$ vanishes when the last two arguments coincide. The third term reduces to $a\left(\mathbf{u}_{h}^{n+1},\mathbf{u}_{h}^{n+1}\right)=\nu\left\|\operatorname*{\bm{\nabla}_{\rm sym}}\mathbf{u}_{h}^{n+1}\right\|_{\Omega}^{2}$. Then by testing (26b) with $p_{h}^{n+1}$, we see that the divergence term vanishes. Applying the following formula to the first term

[TABLE]

with $a=\mathbf{u}_{h}^{n+1}$, $b=\mathbf{u}_{h}^{n}$, $c=\mathbf{u}_{h}^{n-1}$, we obtain

[TABLE]

We are now left with the treatment of the term involving $\mathbf{c}$. We take $\mu_{h}=\lambda_{h}^{n+1}$ in (26e) and $\mathbf{z}_{h}=\frac{3\mathbf{X}_{h}^{n+1}-4\mathbf{X}_{h}^{n}+\mathbf{X}_{h}^{n-1}}{2\Delta t}$ in (26d), and use (26c) to obtain:

[TABLE]

Using the formula (28) we have

[TABLE]

Similarly we bound the term involving the Piola–Kirchhoff stress-tensor $\mathbb{P}$ using again (28):

[TABLE]

Putting all these relations together yields the result and concludes our proof. ∎

Summing up $n=2,\dots,m-1\leq N$ in (27) and dropping out positive terms, yields the following unconditional stability bound:

[TABLE]

We note that the scheme presented in (26a)-(26f) is a fully implicit system which involves the solution of a nonlinear equation, stemming from the convection part in the Navier-Stokes equation as well as from the kinematic coupling term. The problem can be presented in the following matrix form:

[TABLE]

with

[TABLE]

Here $\phi_{i}$, $\psi_{k}$, $\chi_{i}$, and $\zeta_{l}$ denote the basis functions in $\mathbf{V}_{h}$, $Q_{h}$, $\mathbf{S}_{h}$, and $\bm{\Lambda}_{h}$, respectively.

Thanks to the theory developed in [5], we know that the linearization of the system above is associated to a steady saddle point problem which admits a unique solution. Moreover, the finite element discretization is stable, thus giving optimal convergence rates depending on the regularity of the solution.

We can either solve this system by a solver for nonlinear systems of equations like a fixed point iteration or Newton like methods, or make this semi-implicit by replacing the implicit terms with explicit ones. In the numerical experiments reported in the last section, we used a fixed point iteration or the well known extrapolation formula $x^{n+1}=2x^{n}-x^{n-1}+\mathcal{O}(\Delta t^{2})$, in order to replace $A(\mathbf{u}_{h}^{n+1})$ by $A(2\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1})$ and $L_{f}(\mathbf{X}_{h}^{n+1})$ by $L_{f}(2\mathbf{X}_{h}^{n}-\mathbf{X}_{h}^{n-1})$.

3.3. Crank–Nicolson scheme

The Crank–Nicolson scheme is another second order method widely used for the discretization of evolutionary equations thanks to the fact that it is A-stable and one-step. In the literature, we can find two different schemes referred to as Crank–Nicolson scheme, which can be obtained by applying either the midpoint (CNm) or the trapezoidal (CNt) rule to integrate the ODE from $t_{n}$ to $t_{n+1}$. Notice that the two methods coincide when the problem is linear. As far as the Navier–Stokes equations are concerned, in [15] the CNm approach is used and the stability properties are analyzed, while the CNt formula has been introduced, for example, in [17, Chapter 7]. We are going to apply both formulations in order to compare their numerical performances.

Let us start with the CNm version of the method.

Problem 5**.**

Given $\mathbf{u}_{0h}\in\mathbf{V}_{h}$ and $\mathbf{X}_{0h}\in\mathbf{S}_{h}$, for all $n=1,\dots,N$ find $\mathbf{u}_{h}^{n}\in\mathbf{V}_{h}$, $p_{h}^{n}\in Q_{h}$, $\mathbf{X}_{h}^{n}\in\mathbf{S}_{h}$, $\dot{\mathbf{X}}_{h}^{n}\in\mathbf{S}_{h}$, and $\lambda_{h}^{n}\in\bm{\Lambda}_{h}$ such that

[TABLE]

Equation (30e) is obtained by applying the midpoint rule to the kinematic constraint. Hence we approximate the value $\mathbf{u}\left(\mathbf{X}\left(t+\frac{\Delta t}{2}\right),t+\frac{\Delta t}{2}\right)$ and we average both $\mathbf{u}$ and $\mathbf{X}$ providing an approximation which is second order accurate. The following proposition states a stability estimate similar to the one for the BDF2 method.

Proposition 3.3*.*

Assume $\delta\rho\geq 0$ and that the material behavior is governed by an energy density $W$ which is $\mathcal{C}^{1}$ and convex. Let $\mathbf{u}_{h}^{n}\in\mathbf{V}_{h}$ and $\mathbf{X}_{h}^{n}\in\mathbf{S}_{h}$, for $n=1,\dots N$, be the solutions of Problem 5. Then the following estimate holds true:

[TABLE]

Proof.

We mimic again the proof of [3, Prop. 3]. We take $\mathbf{v}_{h}=\left(\mathbf{u}_{h}^{n+1}+\mathbf{u}_{h}^{n}\right)/2$ in (30a) to get:

[TABLE]

The trilinear form $b$ vanishes since the second two arguments coincide. Equation (30b) implies that the terms including the divergence vanish. If the initial condition $\mathbf{u}^{0}_{h}$ is discretely divergence free, the theorem even hold for $n=0$. Hence the above equality reduces to

[TABLE]

We set $\mu_{h}=\lambda^{n+1}$ in (30e), then, using the equations in the solid (30c) with $\mathbf{z}_{h}=(\mathbf{X}_{h}^{n+1}-\mathbf{X}_{h}^{n})/2$ and (30d) with $\mathbf{w}_{h}=(\dot{\mathbf{X}}_{h}^{n+1}-\dot{\mathbf{X}}_{h}^{n})/2$, we arrive at

[TABLE]

Now we want to see how the last term relates to the energy (6). Let us define the function $\mathcal{W}:[0,1]\rightarrow\mathbb{R}$ by

[TABLE]

The convexity of $\mathcal{W}$ is inherited from $W$, so we see $\mathcal{W}^{\prime}(1)\geq\mathcal{W}(1)-\mathcal{W}(0)$. Using the chain rule, we obtain $\mathcal{W}^{\prime}(1)=\mathbb{P}(\mathbb{F}^{n+1}_{h}):(\mathbb{F}^{n+1}_{h}-\mathbb{F}^{n}_{h})$. This is exactly the integrand in the term we want estimate. We obtain

[TABLE]

Putting all equations together yields the result. ∎

Again we can sum the equations with respect to $n$ from [math] to $m-1\leq N$ and get the unconditional stability bound:

[TABLE]

If we use the CNt method, we are led to consider the problem:

Problem 6**.**

Given $\mathbf{u}_{0h}\in\mathbf{V}_{h}$ and $\mathbf{X}_{0}\in\mathbf{S}_{h}$, for all $n=1,\dots,N$ find $\mathbf{u}_{h}^{n}\in\mathbf{V}_{h}$, $p_{h}^{n}\in Q_{h}$, $\mathbf{X}_{h}^{n}\in S_{h}$, $\dot{\mathbf{X}}_{h}^{n}\in\mathbf{S}_{h}$ and $\lambda_{h}^{n}\in\bm{\Lambda}_{h}$, such that:

[TABLE]

We omit the stability analysis of this problem which is not straightforward not even for the Navier–Stokes equations. Nevertheless this scheme has been used in our numerical tests and instabilities seem not to occur.

Both Problems (5) and (6) can be presented in matrix form with the a structure similar to that of BDF2, with some terms properly modified. The solution of the resulting nonlinear system can be obtained with a Newton like solver, Picard iterations, or by linearization. In our numerical experiments we adopted the latter two approaches.

4. Numerical Experiments

In this section we present some numerical results, with the aim of verifying the accuracy of the higher-order schemes presented in the previous sections. We consider fluid and solid with the same density $\rho_{f}=\rho_{s}=\rho$ and same viscosity $\nu_{f}=\nu_{s}=\nu$. In our numerical experiments we adopt either Picard iterations or the semi-implicit schemes obtained by linearization of the nonlinear terms with second order extrapolations. More precisely, we substitute the value of a certain quantity $y(t_{n+1})$ with $2y(t_{n})-y(t_{n-1})$.

The first test case is the deformed annulus considered in [3] where convergence tests illustrate the first order rate for the implicit Euler method. The second test case is the floating disk for which we want to analyze the volume conservation properties of our methods.

In all our experiments, we use triangular meshes both in $\Omega$ and $\mathcal{B}$, enhanced Bercovier–Pironneau elements (i.e. $P1$iso$P2$/$(P1+P0)$) for the $(\mathbf{u},p)$ discretization (see [1]) and $P1$ elements for the discretization of the body’s deformation $\mathbf{X}$, body’s velocity $\dot{\mathbf{X}}$ and the distributed Lagrange multipliers $\lambda$. We evaluate the bilinear form $\mathbf{c}$ by means of the $L^{2}$ scalar product.

At each time step or at each fixed point iteration we have to assemble the matrix in (29) and to solve the resulting algebraic system. First of all, we observe that the matrix is neither symmetric nor positive definite. The sparsity pattern of the matrix is shown in Fig. 2.

Another aspect that influences the computational time is the assembly of the coupling terms. In order to construct the contribution of the matrix $L_{f}(\mathbf{X}_{h}^{n})$, we have to integrate on an element $T_{s}$ in the Lagrangian mesh $\mathcal{T}_{h}^{\mathcal{B}}$ quantities living on the Eulerian mesh such as $\mathbf{v}(\mathbf{X}(\cdot,t),t)$. Numerically, this is done using some quadrature formula on $T_{s}$. Consequently, we have to find all intersecting elements $T_{f}\in\mathcal{T}_{h}$ of the fluid mesh such that $\mathbf{X}(T_{s},t)\cap T_{f}\neq\emptyset$. Essentially this results in checking every deformed solid element against each fluid element. Bounding boxes are a good tool to detect when two cells surely do not intersect. This reduces the computational cost for each cell to cell comparison. Anyway, the computational cost for finding the intersecting cells grows linearly with the number of cells:

[TABLE]

Although, the theoretical analysis allows to choose meshes for the fluid and the solid independently one from each other, it has been observed, for example in [14], that the mesh parameter $h_{s}$ of the solid should be approximately half the size of the mesh parameter $h_{f}$ of the fluid. This can be explained as a trade-off between accuracy and computational effort.

4.1. The deformed annulus

In this test, an annulus is placed at the center of a square filled with some fluid. Initially the annulus is deformed from its initial configuration and the fluid is at rest. The internal forces steer the annulus back into its undeformed configuration and set the surrounding fluid into motion. The material of the annulus is hyperelastic and described by the identity $\mathbb{P}(\mathbb{F})=\kappa\mathbb{F}$. Thanks to the symmetry of the geometry, we run the simulation in the upper right quarter of the domain. The fluid’s domain is the unit square $\Omega=(0,1)^{2}$. On the upper and right boundary we impose no-slip boundary conditions $\mathbf{u}(\mathbf{x},t)=0$ for the fluid velocity. The reference domain for the immersed structure is a section of the annulus: $\mathcal{B}=\{\mathbf{x}\in\mathbb{R}^{2}:x_{1},x_{2}\geq 0,0.3\leq|\mathbf{x}|\leq 0.5\}$. On the remaining part of the boundary, the fluid and the solid are allowed to move along the tangential direction, while the normal component is set to [math]. The initial conditions are

[TABLE]

The meshes for the fluid and the structure are schematically reported in Fig. 3.

We use two meshes whose number of degrees of freedom can be found in Table 1.

In Fig. 1 we report the position of the annulus and the streamlines of the velocity corresponding to the following choice of parameters $\kappa=10$, $\nu=0.1$, $\rho=1$, $T=1$. The BDF2 method is used with $\Delta t=0.05$. The snapshots are taken at $t=0$, $t=0.1$, $t=0.5$, and $t=1$.

In what follows, we report some tests to investigate the higher-order convergence of the BDF2 method and the two version of the Crank–Nicolson scheme CNm and CNt and compare them with the backward Euler method BDF1. The convergence tests are done using the two meshes plotted in Fig. 3 and parameters $\rho=1$, $\nu=1$, $\kappa=10$, and $T=0.2$. Solutions have been calculated for $\Delta t=T/4,\ T/8,\ T/16,\ T/32$. Since in this case no analytic solution is available, for each mesh we computed a reference solution using the BDF2 method with time-step $\Delta t=0.001$ and we report the relative errors with respect to it. When we use Picard iterations to obtain the solution of fully implicit schemes, we calculate the residual in the $L^{2}$-norm and stop the iteration when it is less than the tolerance $\epsilon=10^{-6}$.

We first investigate the convergence of the fully implicit scheme. The relative errors in the $L^{2}$ norm are displayed in Tables 2 and 3.

The backward Euler method shows a clean first order rate of convergence, while the BDF2 method gives a clean second order rate of convergence only on the coarse mesh and in the fine mesh for the fluid velocity. The results for the structure displacement are not very clear on the fine mesh. The second order rate of convergence is achieved only for larger value of the time step and deteriorates as the time step decreases. On the other hand, we see that in this case the error is about $10^{-5}$, which is likely close to the accuracy that our fine mesh can provide. The CNt method achieves second order convergence for fluid and structure on both meshes, while the CNm scheme provides a second order convergence only for the fluid velocity. In general, the CNt scheme seems to perform better from the point of view of the rate of convergence and of the size of the error.

In order to investigate the behavior of the nonlinear solver, we track the maximum number of Picard iterations needed to reach the given tolerance. The results can be found in Tables 4 and 5. We observe that smaller time steps give smaller number of iterates. This can be motivated by the fact that we use the solution at the previous time step as the initial value for the iteration, therefore for smaller time step it is closer to the value at the current time.

In the case of fluid structure interaction, not only the nonlinear convection term has to be assembled at each iteration, but also the one coupling fluid and structure. Therefore, the semi-implicit scheme is attractive as it requires the solution of only one big system at each time steps. Results for the convergence study of the semi-implicit versions of the methods can be found in Tables 6 and 7.

We see a similar behavior as the one given by Picard iterations. The backward Euler method gives a good first order convergence. The BDF2 method shows a higher convergence order, which slows down for smaller time-steps in the fluid domain, in particular it is almost first order for the structure deformation. The behavior of the CNm is not very clear, in fact only the error in the fluid velocity on the coarse mesh achieves second order convergence. The CNt method provides better results which seem to respect the second order of convergence both for fluid velocity and structure deformation independently of the mesh.

Tables 8 and 9 report the residual obtained with the semi-implicit scheme which is higher of two orders of magnitude than the tolerance prescribed in the fixed point iterations. However, this fact does not influence too much the size of the errors.

4.2. The floating disk

In [12] it is pointed out that the Immersed Boundary Method may show poor volume conservation properties. Roy, Heltai, and Costanzo in [20] have compared the volume preservation properties of their numerical method and propose it as a benchmark problem for fluid-structure interaction codes. In the following we investigate the volume preserving properties of the higher-order time-stepping schemes applied to a circular disk placed in a lid-driven cavity. The movement of the fluid induces a motion to the disk, which is being deformed and transported.

The computational domain is a square with unit length: $\Omega=(0,1)^{2}$. The disk has a diameter of $0.2$ and its center is initially placed at $(0.6,0.5)$. We describe Dirichlet boundary conditions on the fluid velocity: on the left, right, and bottom part of the boundary no-slip boundary conditions are enforced; on the upper part the fluid is set to $\mathbf{u}=(1,0)^{\top}$. The fluid and the solid have the same density $\rho_{f}=\rho_{s}=1$. The viscosity is $\nu=0.01$ and the material behaves as $\mathbb{P}(\mathbb{F})=\kappa\mathbb{F}$ with $\kappa=0.1$. The final time of the simulation is $T=4$. For the space discretization the same code as in the previous test is used on finer meshes with $18,818$ DOFs for the fluid velocity and $7,009$ DOFs for the pressure. The dimension of the space $\mathbf{S}_{h}$ is equal to $4,402$. The maximal diameter is $h_{f}=0.029$ for fluid elements, and $h_{s}=0.012$ for structure cells.

The simulation has been run for all the four time marching schemes with a time-step of $\Delta t=0.01$. At each time-step the volume of the immersed solid $\mathcal{B}_{t}$ is calculated and compared to the volume of the non-deformed solid. The percentage of volume change is plotted over time in Fig. 5. Throughout the literature on Immersed Boundary Method this problem has been addressed. Griffith and Luo [13] use a combination of finite differences for the fluid part and finite elements for the structure part. They report a volume conservation, which is more or less equal to results reported here. Wang and Zhan [21] use finite element discretizations for fluid and solid, while the coupling is enforced via some interpolation function, which mimics the Dirac $\delta$ in the continuous problem. They report a much larger volume change than what we have found. Even their volume preserving scheme performs worse. In our numerical experiments, we note that the one step methods give better volume preservation, whereas BDF2 method tends to reduce more the volume of the immersed solid.

The BDF1, the CNm, and the CNt schemes produce the same volume preservation pattern. In the beginning the volume decreases, from $t=1.5$ to $t=3.0$ the volume almost stays constant. In the last part of the time interval, the volume decreases again faster. At the end the volume is decreased by $0.5\%$. The BDF2 method produces the same pattern, but the volume change is larger, at the final time the volume is decreased by $0.7\%$. It is interesting to note that the BDF2 method, which is more accurate than the Backward Euler, produces a bigger volume change.

To investigate further the influence of the discretization parameters $\Delta t$, $h_{f}$, and $h_{s}$ on the volume preservation, the equations have been solved again, once on the same fine mesh with the doubled time-step $\Delta t=0.02$ and once on a coarser mesh with the same time-step $\Delta t=0.01$. The coarse mesh consists of $8,450$ velocity DOFs, $3,137$ pressure DOFs and $1,986$ DOFs for the structure deformation and the Lagrangian multiplier. The volume preservation is plotted over time in Fig. LABEL:VolumePreservationCoarse for the coarse mesh (left) and the coarser time-step (right). One can observe that the qualitative behavior is similar to that reported in Fig. 5 for the fine mesh and double time step $\Delta t=0.02$. The behavior on the coarser mesh, reported in Fig. LABEL:VolumePreservationCoarse (left), presents a larger kink at the end of the interval with respect to that in Fig. 5 and the absolute volume change is approximately doubled, at the end the volume is decreased by $-1.0\%$.

We conclude that in absolute numbers the volume change is the same for both simulations with less accuracy. The spatial discretization has a bigger influence on the qualitative behavior than the time discretization. This results is perfectly compatible with the fact that the area loss is strictly related to the approximation of the divergence free condition which, in turn, depends on the spatial discretization.

With a maximal area change of $-0.7\%$ on the fine mesh with the fine time-step our method performs very well in comparison to other methods. A direct comparison is done with the method developed by Heltai and Costanzo in [14]. Their approach is quite similar to ours. They use a finite element approximation in space based on rectangular meshes and Backward Euler method for the time discretization. Their solver has been run with $33,282$ velocity DOFs, $12,288$ pressure DOFs and $10,370$ DOFs for the displacement field. This resulted in a maximal volume change of $-2.6\%$. A plot over time is shown in Figure 7. It is interesting to note that the volume change is monotone for this method, while our method shows some oscillations. This could be due to the different type of meshes and the fact that intersection of a mapped structural element with the elements in the fluid mesh are simpler to detect.

5. Conclusion

In this paper we discussed a finite element discretization of fluid-structure interaction problems based in the use of a distributed Lagrange multiplier. We introduced three higher-order time-stepping methods: BDF2, CNm, and CNt. We have proved unconditional stability estimates and we have performed a series of numerical tests. In the numerical experiments the BDF2 and CNt methods showed-higher order convergence in the fully implicit version. Using a semi-implicit approach, the CNt method showed second order convergence.

Acknowledgments

The first and the second authors are members of the INdAM Research group GNCS and their research is partially supported by IMATI/CNR and by PRIN/MIUR