# Convergence analysis of an unfitted mesh semi-implicit coupling scheme for incompressible fluid-structure interaction

Erik Burman, Miguel Angel Fernández, Fannie Gerosa

► To cite this version:

Erik Burman, Miguel Angel Fernández, Fannie Gerosa. Convergence analysis of an unfitted mesh semi-implicit coupling scheme for incompressible fluid-structure interaction. Vietnam Journal of Mathematics, 2023, 51, pp.37-69. 10.1007/s10013-022-00589-w. hal-03065803v3

HAL Id: hal-03065803

https://inria.hal.science/hal-03065803v3

Submitted on 27 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Convergence analysis of an unfitted mesh semi-implicit coupling scheme for incompressible fluid-structure interaction

Erik Burman,\* Miguel A. Fernández,† Fannie M. Gerosa‡

October 27, 2022

## Abstract

The complete numerical analysis of time splitting schemes which avoid strong coupling has rarely been addressed in the literature of unfitted mesh methods for incompressible fluid-structure interaction. In this paper, an error analysis of the semi-implicit scheme recently reported in [*Int. J. Numer. Methods Eng.* 2021; 122:5384–5408] is performed for a linear fluid-structure interaction system. The analysis shows that, under a hyperbolic-CFL condition, the leading term in the energy error scales as $O(h^{r-1/2})$, where $r = 1, 2$ stands for the extrapolation order of the solid velocity in the viscous fluid substep. The theoretical findings are illustrated via some numerical experiments which show, in particular, that the considered method avoids the spatial non-uniformity issues of standard loosely coupled schemes and that it delivers practically the same accuracy as the strongly coupled scheme.

## 1 Introduction

This paper is devoted to the error analysis of a numerical method for a linear fluid-structure coupled system involving the transient Stokes equations (in a fixed domain) and a thin-walled solid elastodynamics model. This system is often used as model problem for the analysis of time-splitting schemes for incompressible fluid-structure interaction (see, e.g., [4, 26, 19, 28, 13, 9]). Indeed, it retains the fundamental numerical difficulties that have to be faced also in general incompressible fluid-structure systems. A large amount of added-mass in the system is known to severely compromise stability and accuracy in standard explicit coupling schemes (i.e., those which invoke the fluid and solid solvers only once per time-step, see, e.g., [44, 22, 32, 53]). The simplest approach to overcome these issues is to resort to a strongly coupled scheme (i.e., one in which the interface coupling is implicitly discretized in time), but at the expense of a higher computational complexity.

The development and the analysis of time splitting schemes which avoid strong coupling without compromising stability and accuracy has been a very active field of research during the last fifteen years. The vast majority of the studies have been devoted to the case of spatial approximations based on fluid meshes which are fitted to the interface (see, e.g., [27, 49, 5, 4, 26, 12, 13, 36, 46, 6, 34, 14]). For many applications, such a mesh compatibility can however be cumbersome to maintain in practice (see, e.g., [48, 33, 51, 8, 19, 42, 43]).

The earliest explicit coupling schemes with unfitted meshes have been reported in [8, 43], using the immersed boundary method, and in [19, 42], using unfitted Nitsche approximations with overlapping meshes. Nevertheless, these methods suffer from major stability/accuracy issues which either require severe time-step restrictions (see [8, 19]) or are limited by the amount of added-mass in the system (see [42, 43]). A new class of semi-implicit schemes with unfitted meshes has been recently reported in [30, 1] for the case of the coupling with thin-walled solids. These methods robustly avoid strong coupling but at the expense of introducing additional unknowns in the fluid sub-problem (intermediate solid velocity). Fully explicit variants of these approaches have been derived in [30] and in [10]. Nevertheless, the formulation of the former in the case of immersed solids remains open and the accuracy of the latter relies on a grad-div penalty stabilization (for enhanced mass conservation) which spoils the conditioning of the fluid subsystem.

<footnote>\*Department of Mathematics, University College London, London, UK-WC1E 6BT, United Kingdom</footnote>

<footnote>†Sorbonne Université & CNRS, UMR 7598 LJLL, 75005 Paris France – Inria Paris, 75012 Paris, France</footnote>

<footnote>‡Sorbonne Université & CNRS, UMR 7598 LJLL, 75005 Paris France – Inria Paris, 75012 Paris, France</footnote>
Additional difficulties, which are not considered in this paper, are the cases of curved and dynamic interfaces. Both subjects have not been particularly analyzed in the literature. We refer to [45] for stability and error analysis accounting also geometric errors. When dynamic interfaces are involved, a further complication derives from the variation of the approximation spaces in time. For the study of parabolic problems with moving interfaces we refer to [54, 37].

This paper is devoted to the numerical analysis of the unfitted mesh semi-implicit coupling scheme recently introduced in [29]. The method combines a Nitsche based unfitted mesh spatial approximation with a fractional-step time-marching in the fluid. The viscous part of the coupling is treated in an explicit fashion (which avoids strong coupling), while the remaining fluid pressure and solid contributions are treated implicitly (which guarantees added-mass free stability). Robust a priori error estimates are derived for two extrapolated variants ($r = 1, 2$ stands for the extrapolation of the solid velocity). The analysis highlights the fundamental role played by the time discretization of the Nitsche's penalty term in the stability and accuracy of the splitting. In particular, an $O(\tau^r/h^{\frac{1}{2}})$ splitting error is obtained instead of the standard $O(\tau/h)$ for the stabilized explicit coupling scheme of [19]. The superior accuracy of the method is also supported by numerical experiments in an academic benchmark.

The rest of the paper is organized as follows. Section 2 presents the continuous setting. Its numerical approximation is discussed in Section 3. The numerical analysis of the semi-implicit scheme is reported in Section 4. Numerical evidence on the convergence properties of the methods is given in Section 5. Finally, Section 6 summarizes the main conclusions of this work.

## 2 Problem Setting

We consider a linear fluid-structure interaction problem in which the fluid is described by the Stokes equations in a fixed polyhedral bounded domain $\Omega^f \subset \mathbb{R}^d$, with $d = 2, 3$ and the structure by a linear thin-walled solid model with mid-surface given by $\Sigma$, also assumed polyhedral. Let the boundary of $\Omega^f$ be partitioned as $\partial\Omega^f = \Sigma \cup \Gamma$ and denote the outward unit normal to $\partial\Omega^f$ by $\mathbf{n}$. In this framework, the considered coupled problem reads as follows: find the fluid velocity and pressure $\mathbf{u}: \Omega^f \times \mathbb{R}^+ \to \mathbb{R}^d$, $p: \Omega^f \times \mathbb{R}^+ \to \mathbb{R}$, the solid displacement and velocity $\mathbf{d}: \Sigma \times \mathbb{R}^+ \to \mathbb{R}^d$, $\dot{\mathbf{d}}: \Sigma \times \mathbb{R}^+ \to \mathbb{R}^d$ such that for all $t \in \mathbb{R}^+$ we have

$$ \left\{ \begin{array}{ll} \rho^f \partial_t \mathbf{u} - \mathrm{div} \sigma(\mathbf{u}, p) = \mathbf{0} & \text{in } \Omega^f \times \mathbb{R}^+, \\ \mathrm{div} \mathbf{u} = 0 & \text{in } \Omega^f \times \mathbb{R}^+, \\ \mathbf{u} = \mathbf{0} & \text{on } \Gamma \times \mathbb{R}^+ \end{array} \right. \quad (1) $$

$$ \left\{ \begin{array}{ll} \rho^s \varepsilon \partial_t \dot{\mathbf{d}} + L \mathbf{d} = \mathbf{T} & \text{in } \Sigma \times \mathbb{R}^+, \\ \dot{\mathbf{d}} = \partial_t \mathbf{d} & \text{in } \Sigma \times \mathbb{R}^+, \\ \mathbf{d} = \mathbf{0} & \text{on } \partial \Sigma \times \mathbb{R}^+ \end{array} \right. \quad (2) $$

$$ \left\{ \begin{array}{ll} \mathbf{u} = \dot{\mathbf{d}} & \text{on } \Sigma \times \mathbb{R}^+, \\ \mathbf{T} = -\sigma(\mathbf{u}, p)\mathbf{n} & \text{on } \Sigma \times \mathbb{R}^+, \end{array} \right. \quad (3) $$

with the initial conditions $\mathbf{u}(0) = \mathbf{u}_0$, $\mathbf{d}(0) = \mathbf{d}_0$ and $\dot{\mathbf{d}}(0) = \dot{\mathbf{d}}_0$. Here, the symbols $\rho^f$ and $\rho^s$ stand, respectively, for the fluid and solid densities. The thickness of the solid is denoted by $\varepsilon$ and the fluid Cauchy stress tensor is given by

$$ \sigma(\mathbf{u}, p) \stackrel{\text{def}}{=} 2\mu\epsilon(\mathbf{u}) - p\mathbf{I}, \quad \epsilon(\mathbf{u}) \stackrel{\text{def}}{=} \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T), $$

where $\mu$ denotes the fluid dynamic viscosity and $\mathbf{I}$ is the identity matrix. The relations in (3) enforce, respectively, the kinematic and dynamic interface coupling conditions. The abstract differential operator $\mathbf{L}$ in (2) describes the elastic behaviour of the solid.

In the following, we will make use of the usual Sobolev's spaces $H^m(\Omega)(m \ge 0)$, with norm $\| \cdot \|_{m,\Omega}$ and seminorm $| \cdot |_{m,\Omega}$, along with the closed spaces $H_\Gamma^1(\Omega)$, of functions in $H^1(\Omega)$ with zero trace on $\Gamma$, and $L_0^2(\Omega)$, of functions in $L^2(\Omega)$ with zero mean in $\Omega$. The scalar product in $L^2(\Omega)$ is denoted by $(\cdot, \cdot)_\Omega$.

For the weak formulation of the fluid problem we consider $V = [H_\Gamma^1(\Omega^f)]^d$ and $Q = L^2(\Omega^f)$ as the fluid velocity and pressure functional spaces, respectively. The standard Stokes bilinear forms are given by

$$ a(\mathbf{u}, \mathbf{v}) \stackrel{\text{def}}{=} 2\mu(\epsilon(\mathbf{u}), \epsilon(\mathbf{v}))_{\Omega^f}, \quad b(q, \mathbf{v}) \stackrel{\text{def}}{=} -(q, \mathrm{div}\mathbf{v})_{\Omega^f}. $$
For the solid weak problem we suppose that $W \subset [H_{\partial\Sigma}^1(\Sigma)]^d$ is the space of admissible displacements and we assume that $L: D \subset [L^2(\Sigma)]^d \to [L^2(\Sigma)]^d$ is a self-adjoint second-order differential operator symmetric, coercive and continuous on $W$. Associated to the operator $L$, we define the elastic bilinear form

$$a^s(d, w) \stackrel{\text{def}}{=} (Ld, w)_{\Sigma} \qquad (4)$$

for all $d \in D$ and $w \in W$. We define the following elastic energy norm on

$$\|w\|_s^2 \stackrel{\text{def}}{=} a^s(w, w)$$

and we assume that the following continuity estimate holds

$$\|w\|_s^2 \le \beta^s \|w\|_{1,\Sigma}^2 \qquad (5)$$

for all $w \in W$, with $\beta^s > 0$.

# 3 Numerical methods

In this section, we discuss three numerical methods for the approximations of the coupled problem (1)-(3). These methods involve an unfitted mesh spatial discretization and different levels of fluid-solid time splitting.

## 3.1 Unfitted mesh spatial approximation

In a standard conforming finite element approximation, typically based on fitted meshes (see, e.g., [23, 26]), the kinematic coupling condition (3)1 is strongly enforced. This is no longer feasible in the unfitted mesh setting. We consider the robust and optimal unfitted mesh method with overlapping meshes proposed in [19] (see also the related works [47, 2]). Therein, the interfacial fluid-solid coupling is treated in a fully weak fashion via a Nitsche's type mortaring.

Let be $\{\mathcal{T}_h^s\}_{0<h<1}$ a family of triangulations of $\Sigma$, such that $\Sigma = \bigcup_{K \in \mathcal{T}_h^s}$. We then consider the standard space of continuous piecewise affine functions associated to $\mathcal{T}_h^s$, namely,

$$X_h^s \stackrel{\text{def}}{=} \{v_h \in C^0(\bar{\Sigma}) \mid v_{h|K} \in \mathbb{P}_1(K), \forall K \in \mathcal{T}_h^s\}.$$

For the approximation of the solid discrete space for the displacement and velocity we consider the following space $W_h = [X_h^s]^d \cap W$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Unfitted meshes.</figcaption></figure>

We denote with $\{\mathcal{T}_h^f\}_{0<h<1}$ a family of triangulations that cover the fluid domain $\Omega^f$ such that:

1. Every $\mathcal{T}_h^f$ is fitted to $\Gamma$ but, in general, not to $\Sigma$;

2. For every simplex $K \in \mathcal{T}_h^f$, we have $K \cap \Omega^f \neq \emptyset$.

In what follows, $\Omega_h^f$ stands for the domain covered by $\mathcal{T}_h^f$ (i.e., the fluid computational domain). We denote by $G_h$ the set of elements of $\mathcal{T}_h^f$ intersected by $\Sigma$, by $F_h$ the set of the internal edges or faces of $\mathcal{T}_h^f$, and by $F_{G_h}$ the set of edges or faces of the elements of $G_h$ that do not lie on $\partial\Omega_h^f$, namely,

$$\begin{align*}
\Omega_h^f &\stackrel{\text{def}}{=} \text{int}(\cup_{K \in \mathcal{T}_h^f} K), \quad G_h &\stackrel{\text{def}}{=} \{K \in \mathcal{T}_h^f \mid K \cap \Sigma \neq \emptyset\}, \\
F_h &\stackrel{\text{def}}{=} \{F \in \partial K \mid K \in \mathcal{T}_h^f, F \cap \partial \Omega_h^f \neq F\}, \quad F_{G_h} &\stackrel{\text{def}}{=} \{F \in \partial K \mid K \in G_h, F \cap \partial \Omega_h^f \neq F\}.
\end{align*}$$
We consider the following space of continuous piecewise affine functions defined over $\mathcal{T}_h^f$

$$X_h^f \stackrel{\text{def}}{=} \{v_h \in C^0(\bar{\Omega}_h) \mid v_{h|K} \in \mathbb{P}_1(K) \quad \forall K \in \mathcal{T}_h^f\}. \quad (6)$$

For the approximation of the fluid velocity and pressure, we consider the following spaces

$$V_h^f \stackrel{\text{def}}{=} \{v_h \in [X_h^f]^d \mid v_h|_{\Gamma} = \mathbf{0}\}, \quad Q_h \stackrel{\text{def}}{=} X_h^f.$$

Since the discrete pair $V_h/Q_h$ lacks inf-sup stability, we consider the classical Brezzi-Pitkäranta symmetric pressure stabilization (see, e.g., [11]):

$$s_h(p_h, q_h) \stackrel{\text{def}}{=} \gamma_p \int_{\Omega_h^f} \frac{h^2}{\mu} \nabla p_h \cdot \nabla q_h, \quad \gamma_p > 0 \quad (7)$$

Note that the pressure stabilization is defined over the computational fluid domain $\Omega_h^f$. In order to guarantee robustness of the method with respect to the way the interface is cutting the fluid mesh, we consider the following ghost-penalty operator (see [15]):

$$g_h(\mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} \gamma_g \mu h \sum_{F \in \mathcal{F}_{g_h}} (\llbracket \nabla \mathbf{u}_h \rrbracket_F, \llbracket \nabla \mathbf{v}_h \rrbracket_F)_F. \quad (8)$$

We can hence introduce the following total stabilization operator $S_h$ and associated semi-norm:

$$S_h((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) \stackrel{\text{def}}{=} s_h(p_h, q_h) + g_h(\mathbf{u}_h, \mathbf{v}_h), \\ |(\mathbf{u}_h, p_h)|_S \stackrel{\text{def}}{=} S_h((\mathbf{u}_h, p_h), (\mathbf{u}_h, p_h))^{\frac{1}{2}}, \quad (9)$$

so that the fluid discrete bi-linear form is given by

$$a_h^f((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) \stackrel{\text{def}}{=} a^f((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) + S_h((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)),$$

with

$$a^f((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) \stackrel{\text{def}}{=} a(\mathbf{u}_h, \mathbf{v}_h) - b(p_h, \mathbf{v}_h) + b(q_h, \mathbf{u}_h).$$

Finally, the considered space semi-discrete unfitted mesh approximation of (1)-(3) reads as follows: for $t > 0$, find $(\mathbf{u}_h(t), p_h(t), \dot{\mathbf{d}}_h(t), \mathbf{d}_h(t)) \in V_h \times Q_h \times W_h \times W_h$, such that $\dot{\mathbf{d}}_h = \partial_t \mathbf{d}_h$ and

$$\left\{ \begin{array}{l} \rho^f(\partial_t \mathbf{u}_h, \mathbf{v}_h)_\Omega + a_h^f((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) + \rho^s \varepsilon(\partial_t \dot{\mathbf{d}}_h, \mathbf{w}_h)_\Sigma + a^s(\mathbf{d}_h, \mathbf{w}_h) \\ - (\sigma(\mathbf{u}_h, p_h)\mathbf{n}, \mathbf{v}_h - \mathbf{w}_h)_\Sigma - (\mathbf{u}_h - \dot{\mathbf{d}}_h, \sigma(\mathbf{v}_h, -q_h)\mathbf{n})_\Sigma + \frac{\gamma\mu}{h}(\mathbf{u}_h - \dot{\mathbf{d}}_h, \mathbf{v}_h - \mathbf{w}_h)_\Sigma = 0 \end{array} \right. \quad (10)$$

for all $(\mathbf{v}_h, q_h, \mathbf{w}_h) \in V_h \times Q_h \times W_h$. Here, $\gamma > 0$ denotes the Nitsche's penalty parameter.

**Remark 3.1.** It should be noted that in (10) the Nitsche mortaring is taken only from the fluid side (as in [39]). The fundamental reason for this is that there is no dissipative mechanism in the solid model (2)$_{1,2}$ that can control the perturbation induced by the solid stresses on the interface.

**Remark 3.2.** In the numerical experiments of Section 5, the second assumption on $\mathcal{T}_h^f$, that is, all the elements of $\Omega_h^f$ intersect the physical domain $\Omega^f$ is relaxed. This is achieved by extending the ghost-penalty operator (8) to $\mathcal{F}_h$ (all internal edges or faces of $\mathcal{T}_h^f$), i.e.,

$$g_h(\mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} \gamma_g \mu h \sum_{F \in \mathcal{F}_h} (\llbracket \nabla \mathbf{u}_h \rrbracket_F, \llbracket \nabla \mathbf{v}_h \rrbracket_F)_F. \quad (11)$$

This ensures the invertibility of the fluid stiffness matrix. It should be noted that the results of the numerical analysis reported in Section 4 below also hold in this case.

## 3.2 Time splitting schemes

This section is devoted to the time-discretization of (10). We first discuss the strongly coupled and the stabilized explicit coupling schemes reported in [19]. Particular attention is paid to the well-known accuracy issues of the latter. We then discuss the semi-implicit projection based coupling scheme reported in [29], whose purpose was precisely to circumvent such difficulties without resorting to strong coupling.

In the following, the parameter $\tau > 0$ denotes the time-step length, $\partial_\tau x^n$ stands for the first-order backward difference formula and $x^{*,r}$ represents the r-th order explicit extrapolations to $x^n$, viz.,

$$\partial_\tau x^n \stackrel{\text{def}}{=} \frac{1}{\tau} (x^n - x^{n-1}), \quad x^{*,r} \stackrel{\text{def}}{=} \begin{cases} x^{n-1} & \text{if } r=1, \\ 2x^{n-1} - x^{n-2} & \text{if } r=2. \end{cases} \quad (12)$$
### 3.2.1 Strongly coupled scheme

Traditionally, the natural way to achieve numerical stability has been to consider a strongly coupled scheme, that is, a fully implicit time-discretization of (10). An example of such an approach is reported in Algorithm 1. The method is also known to deliver an optimal $O(\tau) + O(h)$ accuracy in the energy norm (see [30]). The price to pay for this robustness is the resolution (at each time-step) of the hybrid coupled system (13), which can be computationally demanding in practice, particularly, due to its hybrid nature. Indeed, this monolithic system often yields ill-conditioned matrices (see, e.g., [2, 52]) which require dedicated solvers.

**Algorithm 1** Strongly coupled scheme (from [19]).

For $n \ge 1$:

find $(\mathbf{u}_h^n, p_h^n, \dot{\mathbf{d}}_h^n, \dot{d}_h^n) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{W}_h \times \mathbf{W}_h$ with $\dot{\mathbf{d}}_h^n = \partial_\tau \mathbf{d}_h^n$ and such that

$$ \rho^f(\partial_\tau \mathbf{u}_h^n, \mathbf{v}_h)_{\Omega^f} + a_h^f((\mathbf{u}_h^n, p_h^n), (\mathbf{v}_h, q_h)) + \rho^s \varepsilon(\partial_\tau \dot{\mathbf{d}}_h^n, \mathbf{w}_h)_{\Sigma} + a^s(\dot{\mathbf{d}}_h^n, \mathbf{w}_h) \\ - (\sigma(\mathbf{u}_h^n, p_h^n)\mathbf{n}, \mathbf{v}_h - \mathbf{w}_h)_{\Sigma} - (\mathbf{u}_h^n - \dot{\mathbf{d}}_h^n, \sigma(\mathbf{v}_h, -q_h))\mathbf{n})_{\Sigma} + \frac{\gamma\mu}{h}(\mathbf{u}_h^n - \dot{\mathbf{d}}_h^n, \mathbf{v}_h - \mathbf{w}_h)_{\Sigma} = 0 \quad (13) $$

for all $(\mathbf{v}_h, q_h, \mathbf{w}_h) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{W}_h$.

### 3.2.2 Stabilized explicit coupling scheme

The stabilized explicit coupling scheme reported in Algorithm 2 enables a fully sequential decoupled time-marching of (10). Energy stability is achieved under a mild CFL-like condition (see [19]), thanks to the specific explicit treatment of the Nitsche penalty interface term and to the addition of an interface pressure stabilization in time (weakly consistent interfacial compressibility). The stability of the method is independent of the added-mass effect and of the considered local fluid and solid time-marching schemes. These features come however at a price: the sub-optimality of the error with respect to Algorithm 1. The interface time splitting (the explicit treatment of the penalty term in the right-hand side of (14)) introduces a truncation error whose leading contribution scales as $O(\tau/h)$, which prevents convergence when $\tau = O(h)$ (see Remark 4.3 and [18]). The pressure stabilization is introduced to control the pressure fluctuations at the interface, produced by the splitting in time of pressure and it has a contribution of the order $O(\tau^{1/2}h^{1/2})$. Correction iterations are thus needed to enhance accuracy, under restrictive constraints on the discretization parameters.

**Algorithm 2** Stabilized explicit coupling scheme (from [19]).

For $n \ge 1$:

* Solid sub-step: find $(\dot{\mathbf{d}}_h^n, \mathbf{d}_h^n) \in \mathbf{W}_h \times \mathbf{W}_h$ with $\dot{\mathbf{d}}_h^n = \partial_\tau \mathbf{d}_h^n$ and such that

$$ \rho^s \varepsilon(\partial_\tau \dot{\mathbf{d}}_h^n, \mathbf{w}_h)_{\Sigma} + a^s(\dot{\mathbf{d}}_h^n, \mathbf{w}_h) + \frac{\gamma\mu}{h}(\dot{\mathbf{d}}_h^n, \mathbf{w}_h)_{\Sigma} = \frac{\gamma\mu}{h}(\mathbf{u}_h^{n-1}, \mathbf{w}_h)_{\Sigma} - (\sigma(\mathbf{u}_h^{n-1}, p_h^{n-1})\mathbf{n}, \mathbf{w}_h)_{\Sigma} \quad (14) $$

for all $\mathbf{w}_h \in \mathbf{W}_h$.

* Fluid sub-step: find $(\mathbf{u}_h^n, p_h^n) \in \mathbf{V}_h \times \mathbf{Q}_h$, such that

$$ \rho^f(\partial_\tau \mathbf{u}_h^n, \mathbf{v}_h)_{\Omega^f} + a_h^f((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) + \frac{\gamma_0 h}{\gamma\mu}(p_h^n - p_h^{n-1}, q_h)_{\Sigma} \\ - (\mathbf{u}_h^n - \dot{\mathbf{d}}_h^n, q_h \mathbf{n})_{\Sigma} + \frac{\gamma\mu}{h}(\mathbf{u}_h^n, \mathbf{v}_h)_{\Sigma} = \frac{\gamma\mu}{h}(\dot{\mathbf{d}}_h^n, \mathbf{v}_h)_{\Sigma} + (\sigma(\mathbf{u}_h^{n-1}, p_h^{n-1})\mathbf{n}, \mathbf{v}_h)_{\Sigma} \quad (15) $$

for all $(\mathbf{v}_h, q_h) \in \mathbf{V}_h \times \mathbf{Q}_h$.

Roughly speaking, the lack of spatial uniformity in the splitting error of Algorithm 2 can be explained as follows. After spatial refinement, i.e., whenever $h \to 0$, the solid sub-problem (14) forces the solid velocity $\dot{\mathbf{d}}_h^n$ to be close to $\mathbf{u}_h^{n-1}|_{\Sigma}$, whereas in the fluid sub-problem (15) the fluid velocity $\mathbf{u}_h^n|_{\Sigma}$
approximates $\dot{\mathbf{u}}_h^n$. In summary, the spatial discretization forces $\|\mathbf{u}_h^n - \mathbf{u}_h^{n-1}\|_{0,\Sigma}$ to be small as $h \to 0$, by amplifying the time-splitting error. This is an essential ingredient of the scheme that guarantees numerical stability but it degrades accuracy.

### 3.2.3 Projection based semi-implicit coupling scheme

Algorithm 3 reports the projection based semi-implicit scheme of [29]. The fundamental idea of this method, borrowed from [27] in the case of fitted mesh approximations, consists in combining a fractional-step time-marching in the fluid (1) (see, e.g., [35]) with a semi-implicit treatment of the interface coupling (3). In Algorithm 3, the fluid is discretized in time with an incremental pressure-correction method and a backward-Euler method is considered for the solid. We introduce the fluid discrete viscous bi-linear form

$$a_h^f(\mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} 2\mu(\epsilon(\mathbf{u}_h), \epsilon(\mathbf{v}_h))_{\Omega^f}.$$

Note that the fluid viscous-step (16) is explicitly coupled with the solid, hence avoiding strong coupling (i.e., reducing computational complexity), whereas the coupled pressure-displacement system (17) guarantees added-mass free stability through the implicit treatment of the fluid incompressibility and solid inertial effects. For $r = 2$, the Algorithm 3 can be initialized with one step of the scheme with $r = 1$.

**Algorithm 3** Projection-based semi-implicit scheme (from [29]).

For $n \ge r$:

1. Fluid viscous step: find $\tilde{\mathbf{u}}_h^n \in V_h$ such that

$$\frac{\rho^f}{\tau} (\tilde{\mathbf{u}}_h^n - \mathbf{u}_h^{n-1}, \tilde{\mathbf{v}}_h)_{\Omega^f} + a_h^f (\tilde{\mathbf{u}}_h^n, \tilde{\mathbf{v}}_h) + g_h(\tilde{\mathbf{u}}_h^n, \mathbf{v}_h) + (\nabla p_h^{n-1}, \tilde{\mathbf{v}}_h)_{\Omega^f} \\ - (\sigma(\tilde{\mathbf{u}}_h^n, 0)\mathbf{n}, \tilde{\mathbf{v}}_h)_\Sigma + \frac{\gamma\mu}{h} (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{*,r}, \tilde{\mathbf{v}}_h)_\Sigma - (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{*,r}, \sigma(\tilde{\mathbf{v}}_h, 0)\mathbf{n})_\Sigma = 0 \quad (16)$$

for all $\tilde{\mathbf{v}}_h \in V_h$.

2. Pressure-displacement step: find $(p_h^n, d_h^n) \in Q_h \times W_h$ with $\dot{d}_h^n = \partial_\tau d_h^n$, such that

$$\left\{ \begin{aligned} & \frac{\tau}{\rho^f} (\nabla(p_h^n - p_h^{n-1}), \nabla q_h)_{\Omega_h^f} + s_h(p_h^n, q_h) = (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^n, q_h \mathbf{n})_\Sigma - (\text{div} \tilde{\mathbf{u}}_h^n, q_h)_{\Omega_h^f}, \\ & \rho^s \varepsilon (\partial_\tau \dot{\mathbf{d}}_h^n, \mathbf{w}_h)_\Sigma + a^s (\mathbf{d}_h^n, \mathbf{w}_h) = \frac{\gamma \mu}{h} (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{*,r}, \mathbf{w}_h)_\Sigma - (\sigma(\tilde{\mathbf{u}}_h^n, p_h^n)\mathbf{n}, \mathbf{w}_h)_\Sigma \end{aligned} \right. \quad (17)$$

for all $(q_h, \mathbf{w}_h) \in Q_h \times W_h$.

3. Intermediate velocity step: find $\mathbf{u}_h^n \in V_h$ such that

$$\frac{\rho^f}{\tau} (\mathbf{u}_h^n, \mathbf{v}_h)_{\Omega^f} = \frac{\rho^f}{\tau} (\tilde{\mathbf{u}}_h^n, \mathbf{v}_h)_{\Omega^f} - (\nabla(p_h^n - p_h^{n-1}), \mathbf{v}_h)_{\Omega^f} \quad (18)$$

for all $\mathbf{v}_h \in V_h$.

It is worth noting that the discrete interface stresses in the $(17)_2$ involve the same penalty term as in (16). In other words, the viscous stresses in $(17)_2$ correspond to the variationally consistent residual of (16). This constitutes a fundamental difference with respect to Algorithm 2 (and also with respect to [3] with fitted meshes).

The next section provides an error estimate for Algorithm 3 which shows superior accuracy with respect to Algorithm 2, namely: $O(\tau^r/h^{1/2})$, with $r = 1, 2$, instead of $O(\tau/h)$. Furthermore, the numerical evidence reported in Section 5 suggests Algorithm 3 delivers practically the same accuracy as Algorithm 1, which is uniform with respect to $h$. The price to pay for this superior accuracy is threefold:

* An additional CFL-like condition for stability (see [29, Theorem 1] and Theorem 4.1 below), with respect to Algorithm 1;

* The solution of the coupled pressure-displacement system (17), not required in Algorithm 2;

* A limited flexibility in the choice of the time-stepping for the fluid and solid sub-systems.
An additional appealing characteristic of Algorithm 3, is the reduced computational cost with respect to Algorithm 1. We refer to [29] for a run time study in the context of Nitsche-XFEM for immersed thin-walled structures.

**Remark 3.3.** In practice, we can avoid solving the third step, by inserting (18) into the viscous step (16), obtaining

$$ \rho^f(\partial_\tau \tilde{\mathbf{u}}_h^n, \tilde{\mathbf{v}}_h)_{\Omega^f} + a_h^f(\tilde{\mathbf{u}}_h^n, \tilde{\mathbf{v}}_h) + g_h(\tilde{\mathbf{u}}_h^n, \mathbf{v}_h) + (\nabla(2p^{n-1} - p^{n-2}), \tilde{\mathbf{v}}_h)_{\Omega^f} - (\sigma(\tilde{\mathbf{u}}_h^n, 0)\mathbf{n}, \tilde{\mathbf{v}}_h)_{\Sigma} + \frac{\gamma\mu}{h}(\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_{h}^{*,r}, \tilde{\mathbf{v}}_h)_{\Sigma} - (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_{h}^{*,r}, \sigma(\tilde{\mathbf{v}}_{h,i}, 0)\mathbf{n}_i)_{\Sigma} = 0 \\ \text{for all } \tilde{\mathbf{v}}_h \in \mathbf{V}_h. $$

With the purpose of further investigating the space uniformity of the splitting error, we introduce a variant of the projection-based semi-implicit scheme with second-order backward difference discretization. The resulting method is shown in Algorithm 4 by choosing $l = 2$ and with the following definition of the $l$-order backward difference formulas:

$$ \partial_{\tau}^{l} x^{n} \stackrel{\text{def}}{=} \frac{(l+1)x^{n} - 2lx^{n-1} + (l-1)x^{n-2}}{2\tau}, \quad D_{\tau}^{l} \tilde{x}^{n} \stackrel{\text{def}}{=} \frac{(l+1)\tilde{x}^{n} - 2lx^{n-1} + (l-1)x^{n-2}}{2\tau}. $$

Observe that, for $l = 1$ we retrieve Algorithm 3 from [29].

In the following, Algorithm 4 is introduced only for comparison purposes conducted in Section 5. We do not convey any stability analysis or error estimates with second-order in time discretization.

**Algorithm 4** Projection-based semi-implicit scheme with $l$-th order backward difference.

For $n \ge r$:

1. Fluid viscous step: find $\tilde{\mathbf{u}}_h^n \in \mathbf{V}_h$ such that

$$ \rho^f (D_{\tau}^l \tilde{\mathbf{u}}_h^n, \tilde{\mathbf{v}}_h)_{\Omega^f} + a_h^f (\tilde{\mathbf{u}}_h^n, \tilde{\mathbf{v}}_h) + g_h (\tilde{\mathbf{u}}_h^n, \mathbf{v}_h) + (\nabla p_h^{n-1}, \tilde{\mathbf{v}}_h)_{\Omega^f} - (\sigma(\tilde{\mathbf{u}}_h^n, 0)\mathbf{n}, \tilde{\mathbf{v}}_h)_{\Sigma} + \frac{\gamma\mu}{h} (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_{h}^{*,r}, \tilde{\mathbf{v}}_h)_{\Sigma} - (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_{h}^{*,r}, \sigma(\tilde{\mathbf{v}}_h, 0)\mathbf{n})_{\Sigma} = 0 $$

for all $\tilde{\mathbf{v}}_h \in \mathbf{V}_h$.

2. Pressure-displacement step: find $(p_h^n, d_h^n) \in Q_h \times W_h$ with $\dot{d}_h^n = \partial_\tau d_h^n$, such that

$$ \left\{ \begin{array}{l} \frac{\tau}{\rho^f}(\nabla(p_h^n - p_h^{n-1}), \nabla q_h)_{\Omega_h^f} + s_h(p_h^n, q_h) = \frac{l+1}{2}(\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^n, q_h \mathbf{n})_{\Sigma} - \frac{l+1}{2}(\text{div} \tilde{\mathbf{u}}_h^n, q_h)_{\Omega^f}, \\ \rho^s \varepsilon(\partial_\tau^l \dot{\mathbf{d}}_h^n, \mathbf{w}_h)_{\Sigma} + a^s(d_h^n, \mathbf{w}_h) = \frac{\gamma\mu}{h}(\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{*,r}, \mathbf{w}_h)_{\Sigma} - (\sigma(\tilde{\mathbf{u}}_h^n, p_h^n)\mathbf{n}, \mathbf{w}_h)_{\Sigma} \end{array} \right. $$

for all $(q_h, \mathbf{w}_h) \in Q_h \times W_h$.

3. Intermediate velocity step: find $\mathbf{u}_h^n \in \mathbf{V}_h$ such that

$$ \frac{l+1}{2\tau} \rho^f (\mathbf{u}_h^n, \mathbf{v}_h)_{\Omega^f} = \frac{l+1}{2\tau} \rho^f (\tilde{\mathbf{u}}_h^n, \mathbf{v}_h)_{\Omega^f} - (\nabla(p_h^n - p_h^{n-1}), \mathbf{v}_h)_{\Omega^f} $$

for all $\mathbf{v}_h \in \mathbf{V}_h$.

# 4 Numerical analysis

This section is devoted to the numerical analysis of Algorithm 3. We first recall the main ingredients for the energy stability of the method reported in [29] and extend the proof to cover the case of a second-order extrapolation ($r = 2$). An a priori error estimate is derived in Section 4.2.
## 4.1 Energy stability

We assume (see, e.g., [38, 20]) that the following trace inequality holds

$$ \|v\|_{0,\Sigma \cap K}^2 \le C_T (h^{-1} \|v\|_{0,K}^2 + h \|\nabla v\|_{0,K}^2) \quad (19) $$

for all $v \in H^1(K)$, $K \in \mathcal{T}_h^f$ and $C_T$ depending only on $\Sigma$. By combining (19) with a discrete inverse inequality, it follows

$$ \|\boldsymbol{\varepsilon}(\mathbf{v}_h) \boldsymbol{n}\|_{0,\Sigma}^2 \le \sum_{K \in \mathcal{G}_h} \|\boldsymbol{\varepsilon}(\mathbf{v}_h)\|_{0,\Sigma \cap K}^2 \le C_T \sum_{K \in \mathcal{G}_h} \left(h^{-1} \|\boldsymbol{\varepsilon}(\mathbf{v}_h)\|_{0,K}^2 + h \|\nabla \boldsymbol{\varepsilon}(\mathbf{v}_h)\|_{0,K}^2\right) \le \frac{C_{\text{TI}}}{h} \sum_{K \in \mathcal{G}_h} \|\boldsymbol{\varepsilon}(\mathbf{v}_h)\|_{0,K}^2 \quad (20) $$

for all $v_h \in V_h$. Hence,

$$ h \|\boldsymbol{\varepsilon}(\mathbf{v}_h) \boldsymbol{n}\|_{0,\Sigma}^2 \le C_{\text{TI}} \|\boldsymbol{\varepsilon}(\mathbf{v}_h)\|_{0,\Omega_h^f}^2 \quad (21) $$

for all $v_h \in V_h$. This estimates are fundamental for the energy stability of the method.

Let define the discrete total energy $E_h^n$ by the following expression:

$$ E_h^n \stackrel{\text{def}}{=} \frac{\rho^f}{2} \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + \frac{\rho^s \varepsilon}{2} \|\dot{\mathbf{d}}_h^n\|_{\Sigma}^2 + \frac{1}{2} \|\mathbf{d}_h^n\|_{\mathbf{s}}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^n\|_{0,\Omega_h^f}^2 \quad (22) $$

The following result states the conditional energy based stability of the approximation provided by Algorithm 3.

**Theorem 4.1.** Let $\{(\mathbf{u}_h^n, p_h^n, \dot{\mathbf{d}}_h^n, \mathbf{d}_h^n)\}_{n \ge 1}$ be given by Algorithm 3 with $r = 1, 2$. Under the following conditions

$$ \gamma \ge \frac{2(3 + 4\alpha)}{\alpha} \frac{C_{\text{TI}}}{c_g}, \quad (23) $$

$$ \gamma \mu \tau \le \frac{2^{2-r}}{1 + 4\alpha} \rho^s \varepsilon h, \quad (24) $$

with $\alpha > 0$, the discrete energy estimate presented below holds:

$$ E_h^n \lesssim E_h^0, \quad (25) $$

for all $n \ge 1$. As a result, Algorithm 3 is conditionally stable in the energy norm (22).

*Proof.* The proof for $r = 1$ was reported in [29]. We recall here some of the steps and provide some details that include also the case $r = 2$, which follows similarly. We proceed by testing (16)-(18) with $(\mathbf{v}_h, \tilde{\mathbf{v}}_h, q_h, w_h) = (\mathbf{u}_h^n, \tilde{\mathbf{u}}_h^n, p_h^n, \dot{\mathbf{d}}_h^n)$. By proceeding like in [29] (the sole difference lies in the choice of $\dot{\mathbf{d}}_h^{*,r}$), we get the following energy estimate

$$ \frac{\rho^f}{2} \partial_\tau \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + 2c_g \mu \|\epsilon(\tilde{\mathbf{u}}_h^n)\|_{0,\Omega_h^f}^2 + \frac{\tau}{2\rho^f} \left[ \|\nabla p_h^n\|_{0,\Omega_h^f}^2 - \|\nabla p_h^{n-1}\|_{0,\Omega_h^f}^2 \right] + \|(\tilde{\mathbf{u}}_h^n, p_h^n)\|_{\mathbf{s}}^2 \\ \underbrace{-2\mu (\epsilon(\tilde{\mathbf{u}}_h^n)\mathbf{n}, \tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^n)_{\Sigma}}_{T_1} - 2\mu (\epsilon(\tilde{\mathbf{u}}_h^n)\mathbf{n}, \tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{*,r})_{\Sigma} + \underbrace{\frac{\gamma\mu}{h} (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{*,r}, \tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^n)_{\Sigma}}_{T_2} \le 0 \quad (26) $$

for $n \ge r$. Term $T_1$ can be bounded by adding and subtracting $\dot{\mathbf{d}}_h^n$, using the Cauchy–Schwarz, Young's and trace inequalities (20), as follows:

$$ \begin{aligned} T_1 = & -2\mu (\epsilon(\tilde{\mathbf{u}}_h^n)\mathbf{n}, \dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{*,r})_{\Sigma} - 4\mu (\epsilon(\tilde{\mathbf{u}}_h^n)\mathbf{n}, \tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^n)_{\Sigma} \\ & \ge -\left(\frac{1}{\alpha_1} + \frac{2}{\alpha_2}\right) \frac{\mu C_{\text{TI}}}{\gamma} \|\epsilon(\tilde{\mathbf{u}}_h^n)\|_{0,\Omega_h^f}^2 - r\alpha_1 \frac{\gamma\mu}{h} \|\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}\|_{0,\Sigma}^2 - 2(r-1)\alpha_1 \frac{\gamma\mu}{h} \|\dot{\mathbf{d}}_h^{n-1} - \dot{\mathbf{d}}_h^{n-2}\|_{0,\Sigma}^2 \\ & \quad - 2\alpha_2 \frac{\gamma\mu}{h} \|\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2, \end{aligned} \quad (27) $$
with $\alpha_1, \alpha_2 > 0$. Similarly, for the second term, we have

$$ T_2 = \frac{\gamma\mu}{h} \|\tilde{u}_h^n - \dot{\tilde{u}}_h^n\|_{0,\Sigma}^2 + \frac{\gamma\mu}{h} (\dot{\tilde{u}}_h^n - \dot{\tilde{u}}_h^{n,r}, \tilde{u}_h^n - \dot{\tilde{u}}_h^n)_{\Sigma} \\ \geq \left(1 - \frac{1}{2\alpha_3}\right) \frac{\gamma\mu}{h} \|\tilde{u}_h^n - \dot{\tilde{u}}_h^n\|_{0,\Sigma}^2 - r \frac{\alpha_3\gamma\mu}{2h} \|\dot{\tilde{u}}_h^n - \dot{\tilde{u}}_h^{n-1}\|_{0,\Sigma}^2 - (r-1) \frac{\alpha_3\gamma\mu}{h} \|\dot{\tilde{u}}_h^{n-1} - \dot{\tilde{u}}_h^{n-2}\|_{0,\Sigma}^2 \quad (28) $$

with $\alpha_3 > 0$. By inserting (27) and (28) into (26) we get

$$ \frac{\rho^f}{2} \partial_\tau \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + \frac{\tau}{2\rho^f} \left[ \|\nabla p_h^n\|_{0,\Omega_h^f}^2 - \|\nabla p_h^{n-1}\|_{0,\Omega_h^f}^2 \right] + \frac{\rho^s \varepsilon}{2} \partial_\tau \|\dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 + \frac{1}{2\tau} (\|\mathbf{d}_h^n\|_s^2 - \|\mathbf{d}_h^{n-1}\|_s^2) \\ + 2\mu \left[ c_g - \frac{C_{\text{TI}}}{\gamma} \left( \frac{1}{2\alpha_1} + \frac{1}{\alpha_2} \right) \right] \|\epsilon(\tilde{u}_h^n)\|_{0,\Omega_h^f}^2 + \left[ \frac{\rho^s \varepsilon}{2\tau} - r \left( \alpha_1 + \frac{\alpha_3}{2} \right) \frac{\gamma\mu}{h} \right] \|\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}\|_{0,\Sigma}^2 \\ + \left(1 - \frac{1}{2\alpha_3} - 2\alpha_2\right) \frac{\gamma\mu}{h} \|\tilde{u}_h^n - \dot{\tilde{u}}_h^n\|_{0,\Sigma}^2 - (r-1)(2\alpha_1 + \alpha_3) \frac{\gamma\mu}{h} \|\dot{\tilde{u}}_h^{n-1} - \dot{\tilde{u}}_h^{n-2}\|_{0,\Sigma}^2 \le 0 \quad (29) $$

for $n \ge r$. We now set

$$ \alpha_1 = \frac{\alpha}{2}, \quad \alpha_2 = \frac{\alpha}{2(1+2\alpha)}, \quad \alpha_3 = \frac{1}{2} + \alpha, \quad \alpha > 0, $$

so that (29) yields

$$ \frac{\rho^f}{2} \partial_\tau \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + \frac{\tau}{2\rho^f} \left[ \|\nabla p_h^n\|_{0,\Omega_h^f}^2 - \|\nabla p_h^{n-1}\|_{0,\Omega_h^f}^2 \right] + \frac{\rho^s \varepsilon}{2} \partial_\tau \|\dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 + \frac{1}{2\tau} (\|\mathbf{d}_h^n\|_s^2 - \|\mathbf{d}_h^{n-1}\|_s^2) \\ + \left[ \frac{\rho^s \varepsilon}{2\tau} - r \left( \frac{1}{4} + \alpha \right) \frac{\gamma\mu}{h} \right] \|\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}\|_{0,\Sigma}^2 - (r-1) \left( \frac{1}{2} + 2\alpha \right) \frac{\gamma\mu}{h} \|\dot{\mathbf{d}}_h^{n-1} - \dot{\mathbf{d}}_h^{n-2}\|_{0,\Sigma}^2 \\ + 2\mu \left[ c_g - \frac{C_{\text{TI}}}{\gamma} \frac{3+4\alpha}{\alpha} \right] \|\epsilon(\tilde{u}_h^n)\|_{0,\Omega_h^f}^2 + \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \|\tilde{u}_h^n - \dot{\tilde{u}}_h^n\|_{0,\Sigma}^2 \le 0 $$

for $n \ge r$. We now replace the upper index $n$ by $m$ and sum over $m = 2...n$ and multiply by $\tau$. This yields

$$ \frac{\rho^f}{2} \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^n\|_{0,\Omega_h^f}^2 + \frac{\rho^s \varepsilon}{2} \|\dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 + \frac{1}{2\tau} \|\mathbf{d}_h^n\|_s^2 \\ + \left[ \frac{\rho^s \varepsilon}{2} - \tau r \left( \frac{1}{4} + \alpha \right) \frac{\gamma\mu}{h} \right] \sum_{m=1}^n \|\dot{\mathbf{d}}_h^m - \dot{\mathbf{d}}_h^{m-1}\|_{0,\Sigma}^2 - \tau(r-1) \left( \frac{1}{2} + 2\alpha \right) \frac{\gamma\mu}{h} \sum_{m=2}^n \|\dot{\mathbf{d}}_h^{m-1} - \dot{\mathbf{d}}_h^{m-2}\|_{0,\Sigma}^2 \\ + 2\tau\mu \left[ c_g - \frac{C_{\text{TI}}}{\gamma} \frac{3+4\alpha}{\alpha} \right] \sum_{m=2}^n \|\epsilon(\tilde{u}_h^m)\|_{0,\Omega_h^f}^2 + \tau \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \sum_{m=2}^n \|\tilde{u}_h^m - \dot{\tilde{u}}_h^m\|_{0,\Sigma}^2 \\ \leq \frac{\rho^f}{2} \|\mathbf{u}_h^{r-1}\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^{r-1}\|_{0,\Omega_h^f}^2 + \frac{\rho^s \varepsilon}{2} \|\dot{\mathbf{d}}_h^{r-1}\|_{0,\Sigma}^2 + \frac{1}{2} \|\mathbf{d}_h^{r-1}\|_s^2 $$

for $n \ge r$. By rearranging the terms in the first sums, we get

$$ \frac{\rho^f}{2} \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^n\|_{0,\Omega_h^f}^2 + \frac{\rho^s \varepsilon}{2} \|\dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 + \frac{1}{2\tau} \|\mathbf{d}_h^n\|_s^2 \\ + \left[ \frac{\rho^s \varepsilon}{2} - \tau 4^{r-1} \left( \frac{1}{4} + \alpha \right) \frac{\gamma\mu}{h} \right] \sum_{m=2}^n \|\dot{\mathbf{d}}_h^m - \dot{\mathbf{d}}_h^{m-1}\|_{0,\Sigma}^2 - \tau(r-1) \left( \frac{1}{2} + 2\alpha \right) \frac{\gamma\mu}{h} \|\dot{\mathbf{d}}_h^1 - \dot{\mathbf{d}}_h^0\|_{0,\Sigma}^2 \\ + 2\tau\mu \left[ c_g - \frac{C_{\text{TI}}}{\gamma} \frac{3+4\alpha}{\alpha} \right] \sum_{m=2}^n \|\epsilon(\tilde{u}_h^m)\|_{0,\Omega_h^f}^2 + \tau \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \sum_{m=2}^n \|\tilde{u}_h^m - \dot{\tilde{u}}_h^m\|_{0,\Sigma}^2 \\ \leq \frac{\rho^f}{2} \|\mathbf{u}_h^{r-1}\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^{r-1}\|_{0,\Omega_h^f}^2 + \frac{\rho^s \varepsilon}{2} \|\dot{\mathbf{d}}_h^{r-1}\|_{0,\Sigma}^2 + \frac{r}{2} \|\mathbf{d}_h^{r-1}\|_s^2 \quad (30) $$

for $n \ge r$.

In the case $r = 1$, the previous bound yields the energy estimate provided in [29]. For $r = 2$, we need to control the contributions coming from the initialization step, namely,

$$ \frac{\rho^f}{2} \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^n\|_{0,\Omega_h^f}^2 + \frac{\rho^s \varepsilon}{2} \|\dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 + \frac{1}{2} \|\mathbf{d}_h^n\|_s^2, $$
which are obtained from one step of Algorithm 3 with $r = 1$. We hence consider (30) with ($n = 1, r = 1$), which yields

$$ \begin{aligned} & \frac{\rho^f}{2} \|\mathbf{u}_h^1\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^1\|_{0,\Omega_h^f}^2 + \frac{\rho^s\epsilon}{2} \|\dot{\mathbf{d}}_h^1\|_{0,\Sigma}^2 + \frac{1}{2} \|\mathbf{d}_h^1\|_s^2 \\ & \qquad + \left[ \frac{\rho^s\epsilon}{2} - \tau \left(\frac{1}{4} + \alpha\right) \frac{\gamma\mu}{h} \right] \|\dot{\mathbf{d}}_h^1 - \dot{\mathbf{d}}_h^0\|_{0,\Sigma}^2 + 2\tau\mu \left[ c_g - \frac{C_{TI}}{\gamma} \frac{3+4\alpha}{\alpha} \right] \|\epsilon(\tilde{\mathbf{u}}_h^1)\|_{0,\Omega_h^f}^2 \\ & \qquad + \tau \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \|\tilde{\mathbf{u}}_h^1 - \dot{\mathbf{d}}_h^1\|_{0,\Sigma}^2 \le \frac{\rho^f}{2} \|\mathbf{u}_h^0\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^0\|_{0,\Omega_h^f}^2 + \frac{\rho^s\epsilon}{2} \|\dot{\mathbf{d}}_h^0\|_{0,\Sigma}^2 + \frac{1}{2} \|\mathbf{d}_h^0\|_s^2. \end{aligned} $$

Hence, by adding this expression to (30), we finally get

$$ \begin{aligned} & \frac{\rho^f}{2} \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^n\|_{0,\Omega_h^f}^2 + \frac{\rho^s\epsilon}{2} \|\dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 + \frac{1}{2} \|\mathbf{d}_h^n\|_s^2 \\ & \qquad + \left[ \frac{\rho^s\epsilon}{2} - \tau (1+4\alpha) \frac{\gamma\mu}{h} \right] \sum_{m=1}^n \|\dot{\mathbf{d}}_h^m - \dot{\mathbf{d}}_h^{m-1}\|_{0,\Sigma}^2 \\ & \qquad + 2\tau\mu \left[ c_g - \frac{C_{TI}}{\gamma} \frac{3+4\alpha}{\alpha} \right] \sum_{m=1}^n \|\epsilon(\tilde{\mathbf{u}}_h^m)\|_{0,\Omega_h^f}^2 + \tau \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \sum_{m=1}^n \|\tilde{\mathbf{u}}_h^m - \dot{\mathbf{d}}_h^m\|_{0,\Sigma}^2 \\ & \qquad \le \frac{\rho^f}{2} \|\mathbf{u}_h^0\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla p_h^0\|_{0,\Omega_h^f}^2 + \frac{\rho^s\epsilon}{2} \|\dot{\mathbf{d}}_h^0\|_{0,\Sigma}^2 + \frac{1}{2} \|\mathbf{d}_h^0\|_s^2 \end{aligned} \quad (31) $$

for $n \ge 1$. The energy estimate (22) hence follows from (31) under the conditions (23)-(24), which completes the proof. $\square$

**Remark 4.1.** A similar stability analysis can be derived in the case of a thick-walled solid. The solid quantities appearing on the interface $\Sigma$, such as $\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}$, are controlled on the whole solid domain using element-wise trace inequalities. This yields to a parabolic CFL-type stability condition, namely,

$$ \gamma\mu\tau \lesssim \rho^s h^2, $$

which is more restrictive than in the case of a thin-walled solid. An analogous stability result is reported in [27] for the non-incremental version of Algorithm 3 within the framework of fitted mesh. Stability is guaranteed under the CFL-like condition $\rho^f h^2 + 2\mu\tau \lesssim \rho^s \epsilon h$ for a thin-walled solid and $\rho^f h^2 + 2\mu\tau \lesssim \rho^s h^2$ in the case of thick-walled solid.

## 4.2 A priori error estimate

In the following we use the notation $v^n \stackrel{\text{def}}{=} v(t_n)$ (without subscript $h$) to denote the temporal value $v(t_n)$ of a time dependent function $v$. For conciseness, an abuse of notation will be committed by denoting $(\partial_t v)^n$ with $\partial_t v^n$. Furthermore, the symbol $\lesssim$ indicates inequalities up to a multiplicative constant (independent of the discretization parameter $h$ and of the physical parameters). We consider the following mesh-dependent seminorms for functions defined on the interface $\Sigma$:

$$ \|f\|_{\frac{1}{2},h,\Sigma}^2 = \sum_{K \in \mathcal{G}_h} h^{-1} \|f\|_{0,\Sigma_K}^2, \quad \|f\|_{-\frac{1}{2},h,\Sigma}^2 = \sum_{K \in \mathcal{G}_h} h \|f\|_{0,\Sigma_K}^2 $$

where $\Sigma_K$ denotes the part of the interface intersecting the simplex $K$, i.e., $\Sigma_K \stackrel{\text{def}}{=} K \cap \Sigma$.

For the sake of simplicity, in the error analysis we assume that the interface $\Sigma$ is flat. Furthermore, the elements of the solid mesh are supposed to be grouped in disjoint macropatches $P_i$, with meas($P_i$) = $O(h^d)$. Each $(d-1)$-dimensional macro patch $P_i$ is assumed to contain at least one interior node and the union of the $P_i$ is assumed to cover $\Sigma$, viz., $\cup_i F_i = \Sigma$.

The discrete interpolation operators are those introduced in [19] (see also [30]). For the solid displacement, we consider the elastic Ritz-projection operator $\pi_h^s: W \to W_h$ defined by the relation

$$ a^s(w - \pi_h^s w, w_h) = 0 \quad (32) $$

for all $w_h \in W_h$, and for which there holds

$$ \|w - \pi_h^s w\|_{0,\Sigma} + h \|\nabla(w - \pi_h^s w)\|_{0,\Sigma} \lesssim h^2 |w|_{2,\Sigma} \quad (33) $$
for all $w \in [H^2(\Sigma)]^d \cap W$. For the solid velocity, we consider the operator $I_h : W \to W_h$ which is defined as a correction of the operator $\pi_h^s$ by the relation

$$I_h w \stackrel{\text{def}}{=} \pi_h^s w + \sum_i \alpha_i \varphi_i,$$

with $\alpha_i \in \mathbb{R}$ to be fixed with a constraint detailed below. The $\varphi_i$ are functions with support in the macropatches $P_i$, such that

$$0 \le \varphi_i \le 1, \quad \|\varphi_i\|_{0,P_i} \lesssim h^{\frac{d-1}{2}}$$

and take the value 1, component-wise, in the interior nodes of the associated patch $P_i$. The scalars $\alpha_i$ are chosen so that the following condition holds:

$$\int_{P_i} (w - I_h w) \cdot n = 0. \qquad (34)$$

This orthogonality condition is used in the error analysis to control the interface terms coupling the fluid pressure and the solid velocity. We refer to [7] for the detailed construction of such an operator. It can be shown (see [19, Lemma 3.3]) that

$$\|w - I_h w\|_{0,\Sigma} + h \|\nabla (w - I_h w)\|_{0,\Sigma} \lesssim h^2 |w|_{2,\Sigma} \qquad (35)$$

for all $w \in [H^2(\Sigma)]^d \cap W$.

Since the physical solution and the discrete one, are defined on different domains, namely $\Omega^f$ and $\Omega_h^f$, with $\Omega^f \subset \Omega_h^f$, we assume that there exist two linear continuous lifting operators $E_2 : H^2(\Omega^f) \to H^2(\mathbb{R}^d)$ and $E_1 : H^1(\Omega^f) \to H^1(\mathbb{R}^d)$, satisfying the bounds $\|E_1 v\|_{H^1(\mathbb{R}^d)} \lesssim \|v\|_{H^1(\Omega^f)}$ for all $v \in H^1(\Omega^f)$ and $\|E_2 v\|_{H^2(\mathbb{R}^d)} \lesssim \|v\|_{H^2(\Omega^f)}$ for all $v \in H^2(\Omega^f)$, (see, e.g., [25, 50]). To interpolate the resulting extended fluid solution we consider the Scott-Zhang operator $i_{sz}$, see [24] for extra details. Then it holds (see [19, Lemma 3.3]):

$$\begin{aligned}
& \|v - i_{sz} E_2 v\|_{0,\Omega^f} + h \|\nabla (v - i_{sz} E_2 v)\|_{0,\Omega^f} \lesssim h^2 |v|_{2,\Omega^f}, \\
& \|q - i_{sz} E_1 q\|_{0,\Omega^f} + h \|\nabla (q - i_{sz} E_1 q)\|_{0,\Omega^f} \lesssim h |q|_{1,\Omega^f}, \\
& \|\sigma (v - i_{sz} E_2 v, q - i_{sz} E_1 q) n\|_{-\frac{1}{2},h,\Sigma} \lesssim h (\|v\|_{2,\Omega^f} + \|q\|_{1,\Omega^f}), \\
& \|v - i_{sz} E_2 v\|_{\frac{1}{2},h,\Sigma} \lesssim h \|v\|_{2,\Omega^f}, \\
& \|w - I_h w\|_{\frac{1}{2},h,\Sigma} \lesssim h^{\frac{3}{2}} \|w\|_{2,\Sigma}
\end{aligned} \qquad (36)$$

for all $v \in [H^2(\Omega^f)]^d$, $q \in H^1(\Omega^f)$ and $w \in [H^2(\Sigma)]^d \cap W$. Moreover, using an inverse inequality, (36) and the stability of the extension operator we have the following stability result for the gradient projection

$$\|\nabla i_{sz} E_1 q\|_{0,\Omega_h^f}^2 \le h^{-1} \|i_{sz} E_1 q - E_1 q\|_{0,\Omega_h^f}^2 + \|\nabla E_1 q\|_{0,\Omega_h^f}^2 \lesssim \|q\|_{1,\Omega^f}^2. \qquad (37)$$

For the pressure and ghost-penalty stabilization operators (7)-(8), the following consistency properties hold (see, e.g., [19, 21]):

$$s_h (i_{sz} E_1 q, i_{sz} E_1 q)^{\frac{1}{2}} \lesssim \mu^{-\frac{1}{2}} h |q|_{1,\Omega^f} \qquad (38)$$

and

$$g_h (i_{sz} E_2 v, i_{sz} E_2 v)^{\frac{1}{2}} \lesssim h \mu^{\frac{1}{2}} |v|_{2,\Omega^f}. \qquad (39)$$

In the following we will make use of the discrete Grönwall lemma (see, e.g., [41]), which we recall here without a proof.

**section 4.2.** Let $\tau, B$ and $a_m, b_m, c_m, \eta_m$ (for integers $m \ge 1$) be non-negative numbers such that

$$a_n + \tau \sum_{m=1}^{n} b_m \le \tau \sum_{m=1}^{n} \eta_m a_m + \tau \sum_{m=1}^{n} c_m + B$$

for $n \ge 1$. Suppose that $\tau \eta_m < 1$ for all $m \ge 1$. Then there holds

$$a_n + \tau \sum_{m=1}^{n} b_m \le \exp \left( \tau \sum_{m=1}^{n} \frac{\eta_m}{1 - \tau \eta_m} \right) \left( \tau \sum_{m=1}^{n} c_m + B \right)$$

for $n \ge 1$.
For a priori error estimate, we assume that the exact solution of problem (1)-(3) has the following regularity, for a given final time $T \ge \tau$:

$$
\begin{align}
& \mathbf{u} \in [H^1(0, T; H^2(\Omega))]^d, \quad \mathbf{u}|_{\Sigma} \in [H^1(0, T; H^2(\Sigma))]^d, \nonumber \\
& \partial_{tt}\mathbf{u} \in [L^2(0, T; L^2(\Omega))]^d, \partial_{tt}\mathbf{u}|_{\Sigma} \in [L^2(0, T; L^2(\Sigma))]^d, \nonumber \\
& p \in H^1(0, T; H^1(\Omega)), \quad \mathbf{u}|_{\Sigma} \in [H^2(0, T; H^2(\Sigma))]^d, \nonumber \\
& \partial_t^{(r)}\mathbf{u}|_{\Sigma} \in [L^2(0, T; L^2(\Sigma))]^d, \quad \mathbf{d} \in [H^1(0, T; H^2(\Sigma))]^d. \tag{40}
\end{align}
$$

We define the energy norm of the error at time $t_n$ as

$$
\begin{align}
Z_h^n \stackrel{\text{def}}{=} (\rho^f)^{\frac{1}{2}} \|\mathbf{u}^n - \mathbf{u}_h^n\|_{0,\Omega^f} &+ \frac{\tau}{(\rho^f)^{\frac{1}{2}}} \|\nabla(p^n - p_h^n)\|_{0,\Omega_h^f} \nonumber \\
&+ (\rho^s \varepsilon)^{\frac{1}{2}} \|\dot{\mathbf{d}}^n - \dot{\mathbf{d}}_h^n\|_{0,\Sigma} + \|\mathbf{d}^n - \mathbf{d}_h^n\|_{\mathrm{s}} \nonumber \\
&+ \left( \sum_{m=1}^n \tau \frac{\alpha}{1+2\alpha} \gamma \mu \|\tilde{\mathbf{u}}_h^m - \dot{\mathbf{d}}_h^m\|_{\frac{2}{h},\Sigma}^2 \right)^{\frac{1}{2}} \nonumber \\
&+ \left( \sum_{m=1}^n \tau |(\tilde{\mathbf{u}}_h^m, p_h^m)|_{S}^2 \right)^{\frac{1}{2}}.
\end{align}
$$

We can then state the following a priori error estimate.

**Theorem 4.3.** Let $(\mathbf{u}, p, \mathbf{d}, \dot{\mathbf{d}})$ be the solution of (1)-(2) and, $\{(\mathbf{u}_h^n, \tilde{\mathbf{u}}_h^n, p_h^n, \mathbf{d}_h^n, \dot{\mathbf{d}}_h^n)\}_{n \ge r}$ be given by Algorithm 3, with initial data

$$ (\tilde{\mathbf{u}}_h^{r-1}, p_h^{r-1}, \mathbf{d}_h^{r-1}, \dot{\mathbf{d}}_h^{r-1}) = (i_{sz}E_2\mathbf{u}^{r-1}, i_{sz}E_1p^{r-1}, \pi_h^s\mathbf{d}^{r-1}, I_h\dot{\mathbf{d}}^{r-1}) $$

for $r=1,2$. Suppose that the exact solution has the regularity (40) and that the stability conditions (23)-(24) hold. Then, for $n \ge r$ and $n\tau < T$, we have the following discrete error estimate:

$$ Z_h^n \lesssim c_1 h + c_2 \tau + c_3 \frac{\tau^r}{h^{\frac{1}{2}}}, \quad (41) $$

where $\{c_i\}_{i=1}^3$ denote positive constants independent of $h$ and $\tau$, but which depend on the physical parameters and on the regularity of the exact solution.

*Proof.* The proof combines some of the arguments reported in [16, 19]. Note however that analysis of [19] focuses on the spatial semi-discrete problem (10) and the work of [16] is limited to a pure fluid problem. Multiplying (1)-(2) by $(\mathbf{v}_h, q_h) \in \mathbf{V}_h \times Q_h$ and $\mathbf{w}_h \in \mathbf{W}_h$, integrating by parts over $\Omega^f$ and using (1)3 and (3)2 we obtain

1. Fluid sub-problem:

$$
\left\{
\begin{array}{l}
\rho^f (\partial_t \mathbf{u}, \mathbf{v}_h)_{\Omega^f} + 2\mu (\epsilon(\mathbf{u}), \epsilon(\mathbf{v}_h))_{\Omega^f} + (\nabla p, \mathbf{v}_h)_{\Omega^f} - (\sigma(\mathbf{u}, 0) \mathbf{n}, \mathbf{v}_h)_{\Sigma} = 0, \\
(q_h, \text{div } \mathbf{u})_{\Omega^f} = 0.
\end{array}
\right.
\qquad (42)
$$

2. Solid sub-problem:

$$ \rho^s \varepsilon (\partial_t \dot{\mathbf{d}}, \mathbf{w}_h)_\Sigma + a^s (\mathbf{d}, \mathbf{w}_h) + (\sigma (\mathbf{u}, p) \mathbf{n}, \mathbf{w}_h)_\Sigma = 0. $$

Note that only the viscous term has been integrated by parts in the fluid.

On the other hand, owing to the kinematic coupling condition (3)1, we also have

1. Fluid sub-problem:

$$
\left\{
\begin{array}{l}
\rho^f (\partial_t \mathbf{u}, \mathbf{v}_h)_{\Omega^f} + 2\mu (\epsilon(\mathbf{u}), \epsilon(\mathbf{v}_h))_{\Omega^f} + (\nabla p, \mathbf{v}_h)_{\Omega^f} - (\sigma(\mathbf{u}, 0) \mathbf{n}, \mathbf{v}_h)_{\Sigma} \\
\qquad + \frac{\gamma\mu}{h} (\mathbf{u} - \dot{\mathbf{d}}, \mathbf{v}_h)_{\Sigma} - (\mathbf{u} - \dot{\mathbf{d}}, \sigma(\mathbf{v}_h, -q_h) \mathbf{n})_{\Sigma} = 0, \hfill (43) \\
(q_h, \text{div } \mathbf{u})_{\Omega^f} = 0
\end{array}
\right.
$$

for all $(\mathbf{v}_h, q_h) \in \mathbf{V}_h \times Q_h$.

2. Solid sub-problem:

$$ \rho^s \varepsilon (\partial_t \dot{\mathbf{d}}, \mathbf{w}_h)_\Sigma + a^s (\mathbf{d}, \mathbf{w}_h) + \frac{\gamma\mu}{h} (\dot{\mathbf{d}} - \mathbf{u}, \mathbf{w}_h)_\Sigma + (\sigma (\mathbf{u}, p) \mathbf{n}, \mathbf{w}_h)_\Sigma = 0 \quad (44) $$

for all $\mathbf{w}_h \in \mathbf{W}_h$.
Thereafter, using the lifting operators (component-wise) we introduce the following decomposition of the errors for the fluid:

$$
\begin{align}
E_2 \boldsymbol{u}^n - \boldsymbol{u}_h^n &= \underbrace{E_2 \boldsymbol{u}^n - i_{sz} E_2 \boldsymbol{u}^n}_{\stackrel{\text{def}}{=} \boldsymbol{\theta}_{\pi}^n} + \underbrace{i_{sz} E_2 \boldsymbol{u}^n - \boldsymbol{u}_h^n}_{\stackrel{\text{def}}{=} \boldsymbol{\theta}_{h}^n} \quad \text{in } \Omega_h^f, \nonumber \\
E_2 \boldsymbol{u}^n - \tilde{\boldsymbol{u}}_h^n &= \underbrace{E_2 \boldsymbol{u}^n - i_{sz} E_2 \boldsymbol{u}^n}_{\stackrel{\text{def}}{=} \boldsymbol{\theta}_{\pi}^n} + \underbrace{i_{sz} E_2 \boldsymbol{u}^n - \tilde{\boldsymbol{u}}_h^n}_{\stackrel{\text{def}}{=} \tilde{\boldsymbol{\theta}}_{h}^n} \quad \text{in } \Omega_h^f, \tag{45} \\
E_1 p^n - p_h^n &= \underbrace{E_1 p^n - i_{sz} E_1 p^n}_{\stackrel{\text{def}}{=} y_{\pi}^n} + \underbrace{i_{sz} E_1 p^n - p_h^n}_{\stackrel{\text{def}}{=} y_{h}^n} \quad \text{in } \Omega_h^f \nonumber
\end{align}
$$

and for the solid:

$$
\begin{equation}
\begin{aligned}
\boldsymbol{d}^n - \boldsymbol{d}_h^n &= \underbrace{\boldsymbol{d}^n - \boldsymbol{\pi}_h^s \boldsymbol{d}^n}_{\stackrel{\text{def}}{=} \boldsymbol{\xi}_{\pi}^n} + \underbrace{\boldsymbol{\pi}_h^s \boldsymbol{d}^n - \boldsymbol{d}_h^n}_{\stackrel{\text{def}}{=} \boldsymbol{\xi}_{h}^n} && \text{in } \Sigma, \\
\dot{\boldsymbol{d}}^n - \dot{\boldsymbol{d}}_h^n &= \underbrace{\dot{\boldsymbol{d}}^n - \boldsymbol{I}_h \dot{\boldsymbol{d}}^n}_{\stackrel{\text{def}}{=} \dot{\boldsymbol{\xi}}_{\pi}^n} + \underbrace{\boldsymbol{I}_h \dot{\boldsymbol{d}}^n - \dot{\boldsymbol{d}}_h^n}_{\stackrel{\text{def}}{=} \dot{\boldsymbol{\xi}}_{h}^n} && \text{in } \Sigma.
\end{aligned}
\tag{46}
\end{equation}
$$

By adding and subtracting $\partial_{\tau} \pi_h^s d^n$, we can rewrite $\dot{\xi}_h^n$ as

$$
\begin{equation}
\begin{aligned}
\dot{\xi}_h^n = \mathcal{I}_h \dot{d}^n - \dot{d}_h^n &= \underbrace{\mathcal{I}_h \dot{d}^n - \partial_{\tau} \pi_h^s \dot{d}^n}_{\stackrel{\text{def}}{=} \mathcal{Z}_h^n} + \underbrace{\partial_{\tau} \pi_h^s \dot{d}^n - \partial_{\tau} \dot{d}_h^n}_{\stackrel{\text{def}}{=} \partial_{\tau} \xi_h^n} = \mathcal{Z}_h^n + \partial_{\tau} \xi_h^n. \quad (47)
\end{aligned}
\end{equation}
$$

We also introduce the following notations:

$$
\dot{\chi}_{\pi}^{*,r} \stackrel{\text{def}}{=} \dot{d}^n - I_h \dot{d}^{*,r} = \begin{cases} \dot{d}^n - I_h \dot{d}^{n-1} & \text{if } r=1, \\ \dot{d}^n - 2I_h \dot{d}^{n-1} + I_h \dot{d}^{n-2} & \text{if } r=2, \end{cases} \quad (48)
$$

$$
\psi_h^{n-1} \stackrel{\text{def}}{=} y_h^{n-1} + i_{sz} E_1 p^n - i_{sz} E_1 p^{n-1}. \quad (49)
$$

In particular, owing to (48), we have

$$
\mathbf{u}^n - \tilde{\mathbf{u}}_h^n - (\dot{\mathbf{d}}^n - \dot{\mathbf{d}}_h^{*,r}) = \theta_{\pi}^n + \tilde{\theta}_h^n - (\chi_{\pi}^{*,r} + \mathcal{I}_h \dot{\mathbf{d}}^{*,r} - \dot{\mathbf{d}}_h^{n,r}) = \theta_{\pi}^n + \tilde{\theta}_h^n - (\chi_{\pi}^{*,r} + \dot{\xi}_h^{n,r}), \quad (50)
$$

Similar, from (49), one straightforwardly gets the following useful relations:

$$
\begin{align}
p^n - p_h^{n-1} &= \psi_h^{n-1} + y_\pi^n, \tag{51} \\
p_h^n - p_h^{n-1} &= \psi_h^{n-1} - y_h^n. \nonumber
\end{align}
$$

The essential part of the proof focuses on deriving a priori estimate for the discrete errors

$$
\{ (\boldsymbol{\theta}_h^n, \tilde{\boldsymbol{\theta}}_h^n, y_h^n, \boldsymbol{\xi}_h^n, \dot{\boldsymbol{\xi}}_h^n) \}_{n \ge r},
$$

in terms of the following energy norm

$$
\begin{align}
\mathcal{E}_h^n \stackrel{\text{def}}{=} (\rho^f)^{\frac{1}{2}} \|\boldsymbol{\theta}_h^n\|_{0,\Omega_h^f} + \left(\frac{\tau}{\rho^f}\right)^{\frac{1}{2}} \|\nabla y_h^n\|_{0,\Omega_h^f} &+ (\rho^s \varepsilon)^{\frac{1}{2}} \|\dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma} + \|\boldsymbol{\xi}_h^n\|_{\mathrm{s}} \nonumber \\
&+ \left(\sum_{m=r}^{n} \tau \mu \|\nabla \tilde{\boldsymbol{\theta}}_h^m\|_{0,\Omega_h^f}^2\right)^{\frac{1}{2}} + \left(\sum_{m=r}^{n} \tau \frac{\alpha}{1+2\alpha} \gamma \mu \|\tilde{\boldsymbol{\theta}}_h^m - \dot{\boldsymbol{\xi}}_h^m\|_{\frac{2}{1},h,\Sigma}^2\right)^{\frac{1}{2}} + \left(\sum_{m=r}^{n} \tau c_g |(\tilde{\boldsymbol{\theta}}_h^m, y_h^m)|_{\Sigma}^2\right)^{\frac{1}{2}}. \nonumber
\end{align}
$$

To this purpose, we first focus on the fluid subsystem. By subtracting (16) from the momentum equation of (43) at $t = t_n$, with $n \ge r$, we get

$$
\begin{align}
& \rho^f (\partial_t \mathbf{u}^n, \tilde{\mathbf{v}}_h)_{\Omega^f} - \frac{\rho^f}{\tau} (\tilde{\mathbf{u}}_h^n - \mathbf{u}_h^{n-1}, \tilde{\mathbf{v}}_h)_{\Omega^f} + 2\mu (\epsilon(\mathbf{u}^n - \tilde{\mathbf{u}}_h^n), \epsilon(\tilde{\mathbf{v}}_h))_{\Omega^f} + (\nabla(p^n - p_h^{n-1}), \tilde{\mathbf{v}}_h)_{\Omega^f} \nonumber \\
& \qquad + \frac{\gamma\mu}{h} ((\mathbf{u}^n - \tilde{\mathbf{u}}_h^n) - (\dot{\mathbf{d}}^n - \dot{\mathbf{d}}_h^{*,r}), \tilde{\mathbf{v}}_h)_{\Sigma} - (\mathbf{u}^n - \dot{\mathbf{d}}^n, \sigma(\tilde{\mathbf{v}}_h, -q_h)n)_{\Sigma} + (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{*,r}, \sigma(\tilde{\mathbf{v}}_h, 0)n)_{\Sigma}
\end{align}
$$
$$ - (\sigma(\mathbf{u}^n - \tilde{\mathbf{u}}_h^n, 0)\mathbf{n}, \tilde{\mathbf{v}}_h)_{\Sigma} - g_h(\tilde{\mathbf{u}}_h^n, \mathbf{v}_h^n) = 0. \quad (52) $$

Owing to the error decompositions (45)-(46) and using (50)-(51), the identity (52) can be rewritten as

$$ \begin{aligned} & \frac{\rho^f}{\tau} (\tilde{\boldsymbol{\theta}}_h^n - \boldsymbol{\theta}_h^{n-1}, \tilde{\mathbf{v}}_h)_{\Omega^f} + 2\mu(\epsilon(\tilde{\boldsymbol{\theta}}_h^n), \epsilon(\tilde{\mathbf{v}}_h))_{\Omega^f} + (\nabla \psi_h^{n-1}, \tilde{\mathbf{v}}_h)_{\Omega^f} \\ & \qquad + \frac{\gamma\mu}{h} (\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^{*,r}, \tilde{\mathbf{v}}_h)_{\Sigma} - (\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^{*,r}, \sigma(\tilde{\mathbf{v}}_h, 0)\mathbf{n})_{\Sigma} - (\sigma(\tilde{\boldsymbol{\theta}}_h^n, 0)\mathbf{n}, \tilde{\mathbf{v}}_h)_{\Sigma} + g_h(\tilde{\boldsymbol{\theta}}_h^n, \mathbf{v}_h^n) \\ & = \rho^f (-\partial_t \mathbf{u}^n + \partial_\tau \mathbf{u}^n - \partial_\tau \boldsymbol{\theta}_\pi^n, \tilde{\mathbf{v}}_h)_{\Omega^f} - 2\mu(\epsilon(\boldsymbol{\theta}_\pi^n), \epsilon(\tilde{\mathbf{v}}_h))_{\Omega^f} - (\nabla y_\pi^n, \tilde{\mathbf{v}}_h)_{\Omega^f} \\ & \qquad - \frac{\gamma\mu}{h} (\boldsymbol{\theta}_\pi^n - \dot{\boldsymbol{\chi}}_\pi^{*,r}, \tilde{\mathbf{v}}_h)_{\Sigma} + (\boldsymbol{\theta}_\pi^n - \dot{\boldsymbol{\chi}}_\pi^{*,r}, \sigma(\tilde{\mathbf{v}}_h, 0)\mathbf{n})_{\Sigma} + (\sigma(\boldsymbol{\theta}_\pi^n, 0)\mathbf{n}, \tilde{\mathbf{v}}_h)_{\Sigma} + g_h(i_{sz}E_2\mathbf{u}^n, \mathbf{v}_h^n) \end{aligned} \quad (53) $$

for $n \ge r$.

For the pressure, subtracting the pressure-projection step of (17) from the mass conservation equation (43) at $t = t_n$, with $n \ge r$, we get the following relation

$$ (\operatorname{div}(\mathbf{u}^n - \tilde{\mathbf{u}}_h^n), q_h)_{\Omega^f} - \frac{\tau}{\rho^f} (\nabla(p_h^n - p_h^{n-1}), \nabla q_h)_{\Omega_h^f} - s_h(p_h^n, q_h) + (\tilde{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^n, q_h \mathbf{n})_{\Sigma} = 0. \quad (54) $$

Again, using the definition of error decomposition (45)-(46), the coupling kinematic condition (3)$_1$ and (51), from (54) we obtain

$$ \begin{aligned} & \frac{\tau}{\rho^f} (\nabla(y_h^n - \psi_h^{n-1}), \nabla q_h)_{\Omega_h^f} + (\operatorname{div}\tilde{\boldsymbol{\theta}}_h^n, q_h)_{\Omega^f} + s_h(y_h^n, q_h) - (\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n, q_h \mathbf{n})_{\Sigma} \\ & = -(\operatorname{div}\boldsymbol{\theta}_\pi^n, q_h)_{\Omega^f} + s_h(i_{sz}E_1 p_h^n, q_h) + (\boldsymbol{\theta}_\pi^n - \dot{\boldsymbol{\chi}}_\pi^n, q_h \mathbf{n})_{\Sigma} \end{aligned} \quad (55) $$

for $n \ge r$.

Finally, adding and subtracting $i_{sz}E_2\mathbf{u}^n, i_{sz}E_1 p^n, i_{sz}E_1 p^{n-1}$ in (18) and using (51), we obtain the following relation for the end-of-step velocity error

$$ \frac{\rho^f}{\tau} (\boldsymbol{\theta}_h^n - \tilde{\boldsymbol{\theta}}_h^n, \mathbf{v}_h)_{\Omega^f} + (\nabla(y_h^n - \psi_h^{n-1}), \mathbf{v}_h)_{\Omega^f} = 0 \quad (56) $$

for $n \ge r$.

Subtracting the solid problem of (17) from (44) $t = t_n$, with $n \ge r$, and using the relation (50), we obtain

$$ \begin{aligned} & \rho^s \varepsilon (\partial_t \dot{\mathbf{d}}^n - \partial_\tau^* \dot{\mathbf{d}}_h^n, \mathbf{w}_h)_{\Sigma} + a^s (\mathbf{d}^n - \mathbf{d}_h^n, \mathbf{w}_h) - \frac{\gamma\mu}{h} (\mathbf{u}^n - \tilde{\mathbf{u}}_h^n - (\dot{\mathbf{d}}^n - \dot{\mathbf{d}}_h^{*,r}), \mathbf{w}_h)_{\Sigma} \\ & \qquad + (\boldsymbol{\sigma}(\mathbf{u}^n - \tilde{\mathbf{u}}_h^n, p^n - p_h^n)\mathbf{n}, \mathbf{w}_h)_{\Sigma} = 0. \end{aligned} $$

Thus, using (45)-(46) and (50), we finally get the equation for the solid discrete errors:

$$ \begin{aligned} & \rho^s \varepsilon (\partial_\tau \dot{\boldsymbol{\xi}}_h^n, \mathbf{w}_h)_{\Sigma} + a^s (\boldsymbol{\xi}_h^n, \mathbf{w}_h) - \frac{\gamma\mu}{h} (\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^{*,r}, \mathbf{w}_h)_{\Sigma} + (\boldsymbol{\sigma}(\tilde{\boldsymbol{\theta}}_h^n, y_h^n)\mathbf{n}, \mathbf{w}_h)_{\Sigma} \\ & = \rho^s \varepsilon (\partial_\tau \dot{\mathbf{d}}^n - \partial_t \dot{\mathbf{d}}^n - \partial_\tau \dot{\boldsymbol{\xi}}_\pi^n, \mathbf{w}_h)_{\Sigma} - \underbrace{a^s (\boldsymbol{\xi}_\pi^n, \mathbf{w}_h)}_{=0} + \frac{\gamma\mu}{h} (\boldsymbol{\theta}_\pi^n - \dot{\boldsymbol{\chi}}_\pi^{*,r}, \mathbf{w}_h)_{\Sigma} - (\boldsymbol{\sigma}(\boldsymbol{\theta}_\pi^n, y_\pi^n)\mathbf{n}, \mathbf{w}_h)_{\Sigma} \end{aligned} \quad (57) $$

for $n \ge r$. Note that term $a^s(\boldsymbol{\xi}_\pi^n, \mathbf{w}_h)$ vanishes due to the definition of the solid velocity projection operator (32).

Owing to (53), (55), (56) and (57), we have that the discrete errors $(\boldsymbol{\theta}_h^n, \tilde{\boldsymbol{\theta}}_h^n, y_h^n, \boldsymbol{\xi}_h^n, \dot{\boldsymbol{\xi}}_h^n)$ satisfy a time-stepping scheme similar to Algorithm 3, but with a modified right-hand side and pressure increment (i.e., we have $y_h^n - \psi_h^{n-1}$ instead of $p_h^n - p_h^{n-1}$). Therefore, we can leverage the stability arguments of Theorem 4.1 to derive an a priori error estimate. We proceed by testing (53), (55), (56) and (57) with

$$ (\tilde{\mathbf{v}}_h, \mathbf{v}_h, q_h, \mathbf{w}_h) = \tau (\tilde{\boldsymbol{\theta}}_h^n, \boldsymbol{\theta}_h^n, y_h^n, \dot{\boldsymbol{\xi}}_h^n). $$

By adding the resulting expressions, using the steps of Theorem 4.1 under condition (23), we obtain the following energy inequality for the discrete errors:

$$ \frac{\rho^f}{2} (\|\boldsymbol{\theta}_h^n\|_{0,\Omega^f}^2 - \|\boldsymbol{\theta}_h^{n-1}\|_{0,\Omega^f}^2) + \frac{\tau^2}{2\rho^f} (\|\nabla y_h^n\|_{0,\Omega_h^f}^2 - \|\nabla \psi_h^{n-1}\|_{0,\Omega_h^f}^2) + \frac{\rho^s \varepsilon}{2} (\|\dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2 - \|\dot{\boldsymbol{\xi}}_h^{n-1}\|_{0,\Sigma}^2) $$
$$
+ \frac{1}{2} (\|\boldsymbol{\xi}_h^n\|_S^2 - \|\boldsymbol{\xi}_h^{n-1}\|_S^2) + \left[ \frac{\rho^s \varepsilon}{2} - \tau \left(\frac{r}{4} + r\alpha\right) \frac{\gamma\mu}{h} \right] \|\dot{\boldsymbol{\xi}}_h^n - \dot{\boldsymbol{\xi}}_h^{n-1}\|_{0,\Sigma}^2 - (r-1)\tau \left(\frac{1}{2} + 2\alpha\right) \frac{\gamma\mu}{h} \|\dot{\boldsymbol{\xi}}_h^{n-1} - \dot{\boldsymbol{\xi}}_h^{n-2}\|_{0,\Sigma}^2 \\
+ \tau \|(\tilde{\boldsymbol{\theta}}_h^n, y_h^n)\|_S^2 + 2\mu\tau \left(c_g - \frac{C_{TI}}{\gamma} \frac{3+4\alpha}{\alpha}\right) \|\epsilon(\tilde{\boldsymbol{\theta}}_h^n)\|_{0,\Omega_h^f}^2 + \tau \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \|\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2 \leq \mathcal{G}_h^{\star,r} \quad (58)
$$

for $n \ge r$, and with the right-hand side $\mathcal{G}_h^{\star,r}$ defined by

$$
\begin{equation}
\begin{aligned}
\mathcal{G}_h^{\star,r} &\stackrel{\text{def}}{=} \tau \rho^f \left( -\partial_t \mathbf{u}^n + \partial_\tau \mathbf{u}^n - \partial_\tau \boldsymbol{\theta}_\pi^n, \tilde{\boldsymbol{\theta}}_h^n \right)_{\Omega^f} \\
&\quad + \tau \left( y_\pi^n, \operatorname{div} \tilde{\boldsymbol{\theta}}_h^n \right)_{\Omega^f} - 2\tau \mu \left( \epsilon(\boldsymbol{\theta}_\pi^n), \epsilon(\tilde{\boldsymbol{\theta}}_h^n) \right)_{\Omega^f} \\
&\stackrel{\text{def}}{=} T_1 \qquad \stackrel{\text{def}}{=} T_2 \qquad \stackrel{\text{def}}{=} T_3 \\
&\quad + \tau \left( \boldsymbol{\theta}_\pi^n - \dot{\boldsymbol{\chi}}_{\pi}^{\star,r}, \boldsymbol{\sigma}(\tilde{\boldsymbol{\theta}}_h^n, 0) \mathbf{n} \right)_{\Sigma} - \tau \frac{\gamma\mu}{h} \left( \boldsymbol{\theta}_\pi^n - \dot{\boldsymbol{\chi}}_{\pi}^{\star,r}, \tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n \right)_{\Sigma} \\
&\stackrel{\text{def}}{=} T_4 \qquad \stackrel{\text{def}}{=} T_5 \qquad \stackrel{\text{def}}{=} T_6 \\
&\quad - \tau \left( \operatorname{div} \boldsymbol{\theta}_\pi^n, y_h^n \right)_{\Omega^f} + \tau \left( \boldsymbol{\theta}_\pi^n - \dot{\boldsymbol{\xi}}_\pi^n, y_h^n \mathbf{n} \right)_{\Sigma} \\
&\stackrel{\text{def}}{=} T_7 \qquad \stackrel{\text{def}}{=} T_8 \\
&\quad + \tau \rho^s \varepsilon \left( \partial_\tau \dot{\mathbf{d}}^n - \partial_\tau \dot{\mathbf{d}}(t_n) - \partial_\tau \dot{\boldsymbol{\xi}}_\pi^n, \dot{\boldsymbol{\xi}}_h^n \right)_{\Sigma} - \tau a^s (\boldsymbol{\xi}_h^n, z_h^n) \\
&\stackrel{\text{def}}{=} T_9 \qquad \stackrel{\text{def}}{=} T_{10}
\end{aligned}
\tag{59}
\end{equation}
$$

Considering condition (23), equation (58) can be written as:

$$
\begin{aligned}
& \frac{\rho^f}{2} (\|\boldsymbol{\theta}_h^n\|_{0,\Omega_h^f}^2 - \|\boldsymbol{\theta}_h^{n-1}\|_{0,\Omega_h^f}^2) + \frac{\tau^2}{2\rho^f} (\|\nabla y_h^n\|_{0,\Omega_h^f}^2 - \|\nabla \psi_h^{n-1}\|_{0,\Omega_h^f}^2) + \frac{\rho^s \varepsilon}{2} (\|\dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2 - \|\dot{\boldsymbol{\xi}}_h^{n-1}\|_{0,\Sigma}^2) \\
& + \frac{1}{2} (\|\boldsymbol{\xi}_h^n\|_S^2 - \|\boldsymbol{\xi}_h^{n-1}\|_S^2) + \left[ \frac{\rho^s \varepsilon}{2} - \tau \left(\frac{r}{4} + r\alpha\right) \frac{\gamma\mu}{h} \right] \|\dot{\boldsymbol{\xi}}_h^n - \dot{\boldsymbol{\xi}}_h^{n-1}\|_{0,\Sigma}^2 \\
& - \tau(r-1) \left(\frac{1}{2} + 2\alpha\right) \frac{\gamma\mu}{h} \|\dot{\boldsymbol{\xi}}_h^{n-1} - \dot{\boldsymbol{\xi}}_h^{n-2}\|_{0,\Sigma}^2 \\
& + \tau \|(\tilde{\boldsymbol{\theta}}_h^n, y_h^n)\|_S^2 + \mu c_g \tau \|\epsilon(\tilde{\boldsymbol{\theta}}_h^n)\|_{0,\Omega_h^f}^2 + \tau \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \|\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2 \leq \mathcal{G}_h^{\star,r}.
\end{aligned}
$$

The lack of telescoping sum on the pressure terms $\|\nabla y_h^n\|_{0,\Omega_h^f}^2 - \|\nabla \psi_h^{n-1}\|_{0,\Omega_h^f}^2$ is not an issue (see, e.g., [16]). Indeed, using (37) we have

$$
\begin{aligned}
\|\nabla \psi_h^{n-1}\|_{0,\Omega_h^f}^2 &= \|\nabla y_h^{n-1} + \tau \nabla \partial_\tau i_{sz} E_1 p^n\|_{0,\Omega_h^f}^2 \\
&\le \left(1 + \frac{\tau}{T}\right) \|\nabla y_h^{n-1}\|_{0,\Omega_h^f}^2 + \left(1 + \frac{T}{\tau}\right) \tau^2 \|\nabla \partial_\tau i_{sz} E_1 p^n\|_{0,\Omega_h^f}^2 \\
&\lesssim \left(1 + \frac{\tau}{T}\right) \|\nabla y_h^{n-1}\|_{0,\Omega_h^f}^2 + (\tau + T) \|\partial_t p\|_{L^2(t_{n-1}, t_n; H^1(\Omega^f))}^2,
\end{aligned}
$$

so that by inserting this expression into (58), we have

$$
\begin{aligned}
& \frac{\rho^f}{2} (\|\boldsymbol{\theta}_h^n\|_{0,\Omega_h^f}^2 - \|\boldsymbol{\theta}_h^{n-1}\|_{0,\Omega_h^f}^2) + \frac{\tau^2}{2\rho^f} (\|\nabla y_h^n\|_{0,\Omega_h^f}^2 - \|\nabla y_h^{n-1}\|_{0,\Omega_h^f}^2) + \frac{\rho^s \varepsilon}{2} (\|\dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2 - \|\dot{\boldsymbol{\xi}}_h^{n-1}\|_{0,\Sigma}^2) \\
& + \frac{1}{2} (\|\boldsymbol{\xi}_h^n\|_S^2 - \|\dot{\boldsymbol{\xi}}_h^{n-1}\|_S^2) + \left[ \frac{\rho^s \varepsilon}{2} - \tau \left(\frac{r}{4} + r\alpha\right) \frac{\gamma\mu}{h} \right] \|\dot{\boldsymbol{\xi}}_h^n - \dot{\boldsymbol{\xi}}_h^{n-1}\|_{0,\Sigma}^2 - \tau(r-1) \left(\frac{1}{2} + 2\alpha\right) \frac{\gamma\mu}{h} \|\dot{\boldsymbol{\xi}}_h^{n-1} - \dot{\boldsymbol{\xi}}_h^{n-2}\|_{0,\Sigma}^2 \\
& + \tau \|(\tilde{\boldsymbol{\theta}}_h^n, y_h^n)\|_S^2 + \mu c_g \tau \|\epsilon(\tilde{\boldsymbol{\theta}}_h^n)\|_{0,\Omega_h^f}^2 + \tau \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \|\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2 \\
&\lesssim \mathcal{G}_h^{\star,r} + \frac{\tau^3}{2\rho^f T} \|\nabla y_h^{n-1}\|_{0,\Omega_h^f}^2 \\
&\quad + (\tau+T) \frac{\tau^2}{2\rho^f} \|\partial_t p\|_{L^2(t_{n-1}, t_n; H^1(\Omega^f))}^2
\end{aligned}
$$

for $n \ge r$. We now replace the upper index $n$ by $m$ and sum over $m = r \dots n$, this yields

$$
\frac{\rho^f}{2} \|\boldsymbol{\theta}_h^n\|_{0,\Omega_h^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla y_h^n\|_{0,\Omega_h^f}^2 + \frac{\rho^s \varepsilon}{2} \|\dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2 + \frac{1}{2} \|\boldsymbol{\xi}_h^n\|_S^2 + \mu c_g \sum_{m=r}^{n} \tau \|\epsilon(\tilde{\boldsymbol{\theta}}_h^m)\|_{0,\Omega_h^f}^2 + \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \sum_{m=r}^{n} \tau \|\tilde{\boldsymbol{\theta}}_h^m - \dot{\boldsymbol{\xi}}_h^m\|_{0,\Sigma}^2
$$
$$
+ \sum_{m=r}^{n} \tau \| (\tilde{\theta}_h^m, y_h^m) \|_S^2 + \left[ \frac{\rho^s \varepsilon}{2} - \tau 4^r \left( \frac{1}{4} + \alpha \right) \frac{\gamma \mu}{h} \right] \sum_{m=r}^{n} \| \dot{\xi}_h^m - \dot{\xi}_h^{m-1} \|_{0,\Sigma}^2 - \tau(r-1) \left( \frac{1}{2} + 2\alpha \right) \frac{\gamma \mu}{h} \| \dot{\xi}_h^{r-1} - \dot{\xi}_h^{r-2} \|_{0,\Sigma}^2 \\
\lesssim \frac{\tau^2}{2\rho^f T} \sum_{m=r-1}^{n-1} \tau \| \nabla y_h^m \|_{0,\Omega_h^f}^2 + \frac{\rho^f}{2} \| \theta_h^{r-1} \|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \| \nabla y_h^{r-1} \|_{0,\Omega_h^f}^2 + \frac{\rho^s \varepsilon}{2} \| \dot{\xi}_h^{r-1} \|_{0,\Sigma}^2 + \frac{1}{2} \| \xi_h^{r-1} \|_S^2 + \sum_{m=r}^{n} G_h^{m,r} \\
+ (\tau + T) \frac{\tau^2}{2\rho^f} \| \partial_t p \|_{L^2(t_{r-1},T;H^1(\Omega^f))}^2 \quad (60)
$$

for $n \ge r$.

We proceed by estimating $G_h^{*,r}$, by treating each term in (59) separately. The first term can be bound in a standard fashion using a Taylor expansion, (36), the Cauchy–Schwarz and the Poincaré’s inequalities with constant $C_P$. This yields

$$
\begin{align}
T_1 &\le \rho^f \tau (\|\partial_t \mathbf{u}^n - \partial_\tau \mathbf{u}^n\|_{0,\Omega^f} + \|\partial_\tau \boldsymbol{\theta}_\pi^n\|_{0,\Omega^f}) \|\tilde{\boldsymbol{\theta}}_h^n\|_{0,\Omega^f} \nonumber \\
&\le \rho^f \tau (\tau^{\frac{1}{2}} \|\partial_{tt} \mathbf{u}^n\|_{L^2(t_{n-1},t_n;L^2(\Omega^f))} + \tau^{-\frac{1}{2}} \|\partial_t \boldsymbol{\theta}_\pi\|_{L^2(t_{n-1},t_n;L^2(\Omega^f))}) \|\tilde{\boldsymbol{\theta}}_h^n\|_{0,\Omega^f} \nonumber \\
&\le \frac{(\rho^f C_P)^2}{2\varepsilon_1 \mu} (\tau^2 \|\partial_{tt} \mathbf{u}\|_{L^2(t_{n-1},t_n;L^2(\Omega^f))} + \|\partial_t \boldsymbol{\theta}_\pi\|_{L^2(t_{n-1},t_n;L^2(\Omega^f))}^2) + \varepsilon_1 \tau \mu \|\nabla \tilde{\boldsymbol{\theta}}_h^n\|_{0,\Omega_h^f}^2 \tag{61} \\
&\lesssim \frac{(\rho^f C_P)^2}{2\varepsilon_1 \mu} \tau^2 \|\partial_{tt} \mathbf{u}\|_{L^2(t_{n-1},t_n;L^2(\Omega^f))}^2 + \frac{(\rho^f C_P)^2}{2\varepsilon_1 \mu} h^2 \|\partial_t \mathbf{u}\|_{L^2(t_{n-1},t_n;H^2(\Omega^f))}^2 \nonumber \\
&\quad + \varepsilon_1 \tau \mu \|\nabla \tilde{\boldsymbol{\theta}}_h^n\|_{0,\Omega_h^f}^2 \nonumber
\end{align}
$$

with $\varepsilon_1 > 0$. Observe that the last term can be absorbed in the left-hand side of (60) with $\varepsilon_1$ small enough. For term $T_9$ we proceed in a similar fashion. Using (36), we get

$$
\begin{align}
T_9 &\le \rho^s \varepsilon \tau (\| (\partial_t - \partial_\tau) \dot{\mathbf{d}}^n \|_{0,\Sigma} + \| \partial_\tau \dot{\boldsymbol{\xi}}_\pi^n \|_{0,\Sigma} ) \| \dot{\boldsymbol{\xi}}_h^n \|_{0,\Sigma} \nonumber \\
&\le \rho^s \varepsilon \tau (\tau^{\frac{1}{2}} \| \partial_{tt} \mathbf{u} \|_{L^2(t_{n-1},t_n;L^2(\Sigma))} + \tau^{-\frac{1}{2}} \| \partial_t \dot{\boldsymbol{\xi}}_\pi^n \|_{L^2(t_{n-1},t_n;L^2(\Sigma))} ) \| \dot{\boldsymbol{\xi}}_h^n \|_{0,\Sigma} \tag{62} \\
&\le \frac{\rho^s \varepsilon T}{2\varepsilon_9} (\tau^2 \| \partial_{tt} \mathbf{u} \|_{L^2(t_{n-1},t_n;L^2(\Sigma))}^2 + h^2 \| \partial_t \mathbf{u} \|_{L^2(t_{n-1},t_n;H^2(\Sigma))}^2) + \tau \frac{\varepsilon_9 \rho^s \varepsilon}{2T} \| \dot{\boldsymbol{\xi}}_h^n \|_{0,\Sigma}^2, \nonumber
\end{align}
$$

with $\varepsilon_9 > 0$ and where the last term can be controlled in (60) using a Grönwall argument (Lemma 4.2).

For term $T_2$, using (36), we have

$$
T_2 = \tau (y_{\pi}^{n}, \operatorname{div} \tilde{\theta}_{h}^{n})_{\Omega^{f}} \le \tau h^{2} \frac{1}{2\varepsilon_{2}\mu} \|p^{n}\|_{1,\Omega^{f}}^{2} + \tau \frac{\varepsilon_{2}}{2}\mu \|\nabla \tilde{\theta}_{h}^{n}\|_{0,\Omega_{h}^{f}}^{2}. \quad (63)
$$

The second term can be absorbed in the left-hand side of (60) for $\varepsilon_2 > 0$ small enough.

$$
T_3 \lesssim \tau \frac{\mu}{\varepsilon_3} h^2 \|\mathbf{u}^n\|_{2,\Omega}^2 + \tau \varepsilon_3 \mu \|\nabla \tilde{\boldsymbol{\theta}}_h^n\|_{0,\Omega_h^f}^2. \quad (64)
$$

The last term can be, once again, absorbed in the left-hand side of (60), for $\varepsilon_3 > 0$ sufficiently small.

Terms $T_4$ and $T_5$ involve the Nitsche splitting error, namely $\|\boldsymbol{\theta}_\pi^n - \dot{\boldsymbol{\chi}}_\pi^{*,r}\|_{\frac{1}{2},h,\Sigma}$. Using (35), (36) and a Taylor expansion we have

$$
\begin{align}
\| \boldsymbol{\theta}_{\pi}^{n} - \dot{\boldsymbol{\chi}}_{\pi}^{*,r} \|_{\frac{1}{2},h,\Sigma} &\le \| \boldsymbol{\theta}_{\pi}^{n} - \dot{\boldsymbol{\xi}}_{\pi}^{n} \|_{\frac{1}{2},h,\Sigma} + \| \boldsymbol{I}_h \dot{\mathbf{d}}^{n} - \boldsymbol{I}_h \dot{\mathbf{d}}^{*,r} \|_{\frac{1}{2},h,\Sigma} \nonumber \\
&\lesssim h \| \mathbf{u}^n \|_{2,\Omega^f} + h^{\frac{3}{2}} \left( \| \dot{\mathbf{d}}^n \|_{2,\Sigma} + \sum_{i=1}^{r} \| \dot{\mathbf{d}}^{n-r} \|_{2,\Sigma} \right) + \frac{\tau^{\frac{2r-1}{2}}}{h^{\frac{1}{2}}} \| \partial_t^r \dot{\mathbf{d}} \|_{L^2(t_{n-r},t_n;L^2(\Sigma))}. \tag{65}
\end{align}
$$

To estimate $T_4$, we follow the same idea of [19]. Using the robust trace inequality (21) combined with (65), we get

$$
\begin{align}
T_4 &\le \tau \frac{\mu}{\varepsilon_4} \|\boldsymbol{\theta}_{\pi}^n - \dot{\boldsymbol{\chi}}_{\pi}^{*,r}\|_{\frac{1}{2},h,\Sigma}^2 + \tau \mu \varepsilon_4 \|\boldsymbol{\epsilon}(\tilde{\boldsymbol{\theta}}_h^n)\boldsymbol{n}\|_{-\frac{1}{2},h,\Sigma}^2 \nonumber \\
&\lesssim \tau \frac{\mu}{\varepsilon_4} h^2 \left( \|\mathbf{u}^n\|_{2,\Omega}^2 + \sum_{i=1}^r \|\dot{\mathbf{d}}^{n-r}\|_{2,\Sigma} \right) + \frac{\mu}{\varepsilon_4} \frac{\tau^{2r}}{h} \|\partial_t^r \dot{\mathbf{d}}\|_{L^2(t_{n-r},t_n;L^2(\Sigma))}^2 + \tau \varepsilon_4 \mu \|\nabla \tilde{\boldsymbol{\theta}}_h^n\|_{0,\Omega_h^f}^2. \tag{66}
\end{align}
$$
Once more, the last term can be absorbed in the left-hand side of (60), for $\epsilon_4 > 0$ sufficiently small. Similarly, for $T_5$ we have

$$ T_5 \le \frac{\tau}{2\epsilon_5} \gamma\mu \|\boldsymbol{\theta}_{\pi}^n - \dot{\boldsymbol{\chi}}_{\pi}^{*,r}\|_{\frac{2}{3},h,\Sigma}^2 + \frac{\tau\epsilon_5}{2} \gamma\mu \|\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n\|_{\frac{2}{3},h,\Sigma}^2 \\ \le \tau \frac{\gamma\mu}{\epsilon_5^2} h^2 \left( \|\boldsymbol{u}^n\|_{2,\Omega}^2 + \sum_{i=1}^r \|\dot{\boldsymbol{d}}^{n-r}\|_{2,\Sigma} \right) + \frac{\gamma\mu}{\epsilon_6} \frac{\tau^{2r}}{h} \|\partial_t^r \dot{\boldsymbol{d}}\|_{L^2(t_{n-r},t_n;L^2(\Sigma))}^2 + \frac{\tau\epsilon_5}{2} \gamma\mu \|\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n\|_{\frac{2}{3},h,\Sigma}^2. \quad (67) $$

Note that the last term can be included in the left-hand side of (60) for $\epsilon_5 > 0$ small enough.

Term $T_6$ can be handled using (36) as follows:

$$ T_6 \le \tau \frac{1}{2\epsilon_6} \frac{1}{\gamma\mu} \|\boldsymbol{\sigma}(\boldsymbol{\theta}_{\pi}^n, y_{\pi}^n) \boldsymbol{n}\|_{-\frac{1}{2},h,\Sigma}^2 + \tau \frac{\epsilon_6}{2} \gamma\mu \|\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n\|_{\frac{2}{3},h,\Sigma}^2 \\ \le \tau h^2 \frac{1}{2\epsilon_6\gamma\mu} (\|\boldsymbol{u}^n\|_{2,\Omega^f}^2 + \|p^n\|_{1,\Omega^f}^2) + \tau \frac{\epsilon_6}{2} \gamma\mu \|\tilde{\boldsymbol{\theta}}_h^n - \dot{\boldsymbol{\xi}}_h^n\|_{\frac{2}{3},h,\Sigma}^2. \quad (68) $$

Again, the last term can be absorbed in the left-hand side of (60), for $\epsilon_6 > 0$ small enough.

Integrating by parts in $T_7$, we have

$$ T_7 = \underbrace{\tau(\boldsymbol{\theta}_{\pi}^n, \nabla y_h^n)_{\Omega^f}}_{\text{def}} + \underbrace{\tau(\dot{\boldsymbol{\xi}}_{\pi}^n, \tau y_h^n \boldsymbol{n})_{\Sigma}}_{\text{def}}. \quad (69) $$

Term $T_{7,1}$ can be easily handled by using (36) as follows:

$$ T_{7,1} \le \tau \frac{1}{2\epsilon_7} \|\boldsymbol{\theta}_{\pi}^n\|_{0,\Omega^f}^2 + \tau \frac{\epsilon_7}{2} \|\nabla y_h^n\|_{0,\Omega^f}^2 \le \tau h^2 \frac{\mu}{2\epsilon_7} \|\boldsymbol{u}^n\|_{0,\Omega^f}^2 + \tau \frac{\epsilon_7}{2} |(0, y_h^n)|_S^2. \quad (70) $$

The second term can be absorbed in the left-hand side of (60), for $\epsilon_7 > 0$ sufficiently small. For the second term of (69), we proceed as in [19] (see also [30]). We denote by $\tilde{y}_i^n \in \mathbb{R}$ the average of $y_h^n$ over the interface patch $P_i$. Combining the trace inequality (19) with the orthogonality property (34) of the interpolation operator $I_h$ and the standard estimate

$$ \|y_h^n - \tilde{y}_i^n\|_{0,P_i} \lesssim h \|\nabla y_h^n\|_{0,P_i}, $$

term $T_{7,2}$ can be estimated as follows:

$$ T_{7,2} = -\tau \sum_i (y_h^n - \tilde{y}_i^n, \dot{\xi}_\pi^n \cdot \boldsymbol{n})_{P_i} \lesssim \tau \sum_i h \|\nabla y_h^n\|_{0,P_i} \|\dot{\xi}_\pi^n\|_{0,P_i} \\ \lesssim \tau h^3 \frac{\mu}{2\epsilon_7} \|\dot{\boldsymbol{d}}^n\|_{2,\Sigma}^2 + \tau h^2 \frac{\epsilon_7}{2\mu} \|\nabla y_h^n\|_{0,\Omega_h}^2 \lesssim \tau h^3 \frac{\mu}{2\epsilon_7} \|\dot{\boldsymbol{d}}^n\|_{2,\Sigma}^2 + \tau \frac{\epsilon_7}{2} |(0, y_h^n)|_S^2. \quad (71) $$

It should be noted here we have assumed that the solid mesh step has an asymptotic regime similar to the fluid mesh step, namely, $h^s = O(h^f)$. As for $T_{7,1}$, the last term in (71) can be, once again, absorbed in the left-hand side of (60), for $\epsilon_7 > 0$ sufficiently small.

For term $T_8$, using the weak consistency of the stabilization operators (38) and (39), we have

$$ T_8 \le \tau \frac{1}{2\epsilon_8} |(i_{sz} E_2 \boldsymbol{u}^n, i_{sz} E_1 p^n)|_S^2 + \tau \frac{\epsilon_8}{2} |(\tilde{\boldsymbol{\theta}}_h^n, y_h^n)|_S^2 \\ \le \tau h^2 \frac{1}{\epsilon_8 \mu} \left( \mu \|\boldsymbol{u}^n\|_{2,\Omega^f}^2 + \mu^{-1} \|p^n\|_{1,\Omega^f}^2 \right) + \tau \frac{\epsilon_8}{2} |(\tilde{\boldsymbol{\theta}}_h^n, y_h^n)|_S^2. \quad (72) $$

Again, the last term can be absorbed in the left-hand side of (60), for $\epsilon_8 > 0$ small enough.

Term $T_{10}$ can be bounded using the continuity estimate for the elastic bilinear form (5), (32), (36) and a triangular inequality. This yields

$$ \begin{align} T_{10} &= -\tau a^s (\xi_h^n, I_h \dot{\boldsymbol{d}}^n - \partial_\tau \pi_h^s \boldsymbol{d}^n) \nonumber \\ &= -\tau a^s (\xi_h^n, I_h \dot{\boldsymbol{d}}_s^n - \partial_\tau \boldsymbol{d}^n) \nonumber \\ &\le \tau \|\xi_h^n\|_S \|I_h \dot{\boldsymbol{d}}^n - \partial_\tau \boldsymbol{d}^n\|_S \tag{73} \\ &\le \tau \|\xi_h^n\|_S (\|I_h \dot{\boldsymbol{d}}^n - \dot{\boldsymbol{d}}^n\|_S + \|\dot{\boldsymbol{d}}^n - \partial_\tau \boldsymbol{d}^n\|_S) \nonumber \\ &\le \tau \frac{\epsilon_{11}}{2T} \|\xi_h^n\|_S^2 + \tau \frac{h^2 \beta^s T}{2\epsilon_{10}} \|\boldsymbol{u}^n\|_{2,\Sigma}^2 + \tau^2 \frac{\beta^s T}{2\epsilon_{10}} \|\partial_t \boldsymbol{u}\|_{L^2(t_{n-1},t_n;H^1(\Sigma))}^2. \nonumber \end{align} $$
Note that the first term can controlled via Lemma 4.2 in (60).

In summary, the term $\sum_{m=r}^{n} G_h^{m,r}$ in the right-hand side of (60) can be estimated by collecting the estimates (61)-(73) and by inserting them into (59), for $n \ge r$. The desired estimate (41) hence follows from (60) together with the stability condition (24) and Lemma 4.2 with

$$a_m = \frac{\rho^f}{2} \|\boldsymbol{\theta}_h^m\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla y_h^n\|_{0,\Omega_h^f}^2 + \frac{\rho^s\varepsilon}{2} \|\dot{\boldsymbol{\xi}}_h^m\|_{0,\Sigma}^2 + \frac{1}{2} \|\boldsymbol{\xi}_h^n\|_s^2, \quad \eta_m = \frac{1}{T} \quad (74)$$

and by noting that, owing to the initial data, we have

$$\boldsymbol{\theta}_h^{r-1} = \mathbf{0}, \quad \tilde{\boldsymbol{\theta}}_h^{r-1} = \mathbf{0}, \quad y_h^{r-1} = \mathbf{0}, \quad \boldsymbol{\xi}_h^{r-1} = \mathbf{0}, \quad \dot{\boldsymbol{\xi}}_h^{r-1} = \mathbf{0} \quad (75)$$

for $r = 1, 2$. This completes proof. $\square$

**Corollary 4.1.** Assume that Algorithm 3 with $r = 2$ is initialized with one step of the method with $r = 1$. Then, under the assumptions of Theorem 4.3, for $n \ge 1$ and $n\tau < T$, the following discrete error estimate holds for the scheme with $r = 2$:

$$\mathcal{Z}_h^n \lesssim c_1 h + c_2 \tau + c_3 \frac{\tau^2}{h^{\frac{3}{2}}} + c_4 \frac{\tau^{\frac{3}{2}}}{h^{\frac{1}{2}}}, \quad (76)$$

where $\{c_i\}_{i=1}^4$ denote positive constants independent of $h$ and $\tau$, but which depend on the physical parameters and on the regularity of the exact solution.

*Proof.* For $r = 2$, we have to bound the contributions from the initialization step in the right-hand side of (60), namely,

$$\frac{\rho^f}{2} \|\boldsymbol{\theta}_h^1\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla y_h^1\|_{0,\Omega_h^f}^2 + \frac{\rho^s\varepsilon}{2} \|\dot{\boldsymbol{\xi}}_h^1\|_{0,\Sigma}^2 + \frac{1}{2} \|\boldsymbol{\xi}_h^1\|_s^2 + \mu c_g \tau \|\epsilon(\tilde{\boldsymbol{\theta}}_h^1)\|_{0,\Omega_h^f}^2 + \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \tau \|\tilde{\boldsymbol{\theta}}_h^1 - \dot{\boldsymbol{\xi}}_h^1\|_{0,\Sigma}^2 + \tau \|(\tilde{\boldsymbol{\theta}}_h^1, y_h^1)\|_S^2 + \left[\frac{\rho^s\varepsilon}{2} - \tau \left(\frac{1}{4} + \alpha\right) \frac{\gamma\mu}{h}\right] \|\dot{\boldsymbol{\xi}}_h^1 - \dot{\boldsymbol{\xi}}_h^0\|_{0,\Sigma}^2 \lesssim \mathcal{G}_h^{1,1} + \frac{\tau^3}{2\rho^f} \|\partial_t p\|_{L^2(0,\tau; H^1(\Omega^f))}^2. \quad (77)$$

To this purpose, we use the fact that the initialization of Algorithm 3 with $r = 2$ is provided by the first step of the scheme with $r = 1$. We can hence use the estimate provided by (60) with $r = 1$ and $n = 1$ to control (77). More precisely, using (75), this yields

$$\frac{\rho^f}{2} \|\boldsymbol{\theta}_h^1\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla y_h^1\|_{0,\Omega_h^f}^2 + \frac{\rho^s\varepsilon}{2} \|\dot{\boldsymbol{\xi}}_h^1\|_{0,\Sigma}^2 + \frac{1}{2} \|\boldsymbol{\xi}_h^1\|_s^2 + \mu c_g \tau \|\epsilon(\tilde{\boldsymbol{\theta}}_h^1)\|_{0,\Omega_h^f}^2 + \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \tau \|\tilde{\boldsymbol{\theta}}_h^1 - \dot{\boldsymbol{\xi}}_h^1\|_{0,\Sigma}^2 + \tau \|(\tilde{\boldsymbol{\theta}}_h^1, y_h^1)\|_S^2 + \left[\frac{\rho^s\varepsilon}{2} - \tau \left(\frac{1}{4} + \alpha\right) \frac{\gamma\mu}{h}\right] \|\dot{\boldsymbol{\xi}}_h^1 - \dot{\boldsymbol{\xi}}_h^0\|_{0,\Sigma}^2 \lesssim \mathcal{G}_h^{1,1} + \frac{\tau^3}{2\rho^f} \|\partial_t p\|_{L^2(0,\tau; H^1(\Omega^f))}^2. \quad (78)$$

Hence, by inserting this estimate in (60), we get

$$\frac{\rho^f}{2} \|\boldsymbol{\theta}_h^n\|_{0,\Omega^f}^2 + \frac{\tau^2}{2\rho^f} \|\nabla y_h^n\|_{0,\Omega_h^f}^2 + \frac{\rho^s\varepsilon}{2} \|\dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2 + \frac{1}{2} \|\boldsymbol{\xi}_h^n\|_s^2 + \mu c_g \sum_{m=1}^n \tau \|\epsilon(\tilde{\boldsymbol{\theta}}_h^m)\|_{0,\Omega_h^f}^2 + \frac{\alpha}{1+2\alpha} \frac{\gamma\mu}{h} \tau \|\tilde{\boldsymbol{\theta}}_h^m - \dot{\boldsymbol{\xi}}_h^m\|_{0,\Sigma}^2 + \sum_{m=1}^n \tau \|(\tilde{\boldsymbol{\theta}}_h^m, y_h^m)\|_S^2 + \left[\frac{\rho^s\varepsilon}{2} - \tau(1+4\alpha)\frac{\gamma\mu}{h}\right] \sum_{m=1}^n \|\dot{\boldsymbol{\xi}}_h^m - \dot{\boldsymbol{\xi}}_h^{m-1}\|_{0,\Sigma}^2 \lesssim \frac{\tau^2}{2\rho^f T} \sum_{m=1}^{n-1} \tau \|\nabla y_h^m\|_{0,\Omega_h^f}^2 + \mathcal{G}_h^{1,1} + \sum_{m=2}^n \mathcal{G}_h^{m,2} + (\tau+T) \frac{\tau^2}{2\rho^f} \|\partial_t p\|_{L^2(0,\tau; H^1(\Omega^f))}^2 \quad (79)$$

for $n \ge 1$. Owing to the initialization procedure, the bounds provided in (66)-(67) for terms $T_4$ and $T_5$ of $\mathcal{G}_h^{1,1}$ yield a $O(\tau^3/h)$ splitting error, by noting that

$$\gamma\mu \frac{\tau^2}{h} \|\partial_t \dot{\boldsymbol{d}}\|_{L^2(0,\tau; L^2(\Sigma))}^2 \le \gamma\mu \frac{\tau^3}{h} \|\partial_t \dot{\boldsymbol{d}}\|_{L^\infty(0,\tau; L^2(\Sigma))}^2. \quad (79)$$

The estimate (41) for $r = 2$ hence follows from (78) together with the stability condition (24) and Lemma 4.2 with (74). This completes proof. $\square$

We conclude this section with a series of remarks.

**Remark 4.2.** For $r = 2$, the last term in (41) comes from the bound of the first step of Algorithm 3 with $r = 1$, that is, the estimate given by (79). This bound is quasi-optimal in time because the Taylor expansions are evaluated in $L^2(0, T)$ instead of $L^1(0, T)$. Alternatively, one could avoid this term by initializing Algorithm 3 with the first-step of Algorithm 1.
**Remark 4.3.** Note that Algorithm 2 introduces the following perturbations terms in the discrete error equation

$$ \frac{\gamma\mu}{h}(\mathbf{u}_h^n - \mathbf{u}_h^{n-1}, \mathbf{w}_h)_{0,\Sigma} + (\sigma(\mathbf{u}_h^n, p_h^n)\mathbf{n} - \sigma(\mathbf{u}_h^{n-1}, p_h^{n-1})\mathbf{n}, \mathbf{v}_h - \mathbf{w}_h)_{0,\Sigma}, $$
with $\mathbf{w}_h = \tau\dot{\boldsymbol{\xi}}_h^n$. The first term leads to the following bound:

$$ \tau \frac{\gamma\mu}{h} (\mathbf{u}^n - \mathbf{u}^{n-1}, \dot{\boldsymbol{\xi}}_h^n)_{0,\Sigma} \leq \frac{(\gamma\mu)^2}{2} \frac{\tau^2}{h^2} \|\partial_t \mathbf{u}\|_{L^2(t_{n-1}, t_n; L^2(\Sigma))}^2 + \frac{\tau}{2} \|\dot{\boldsymbol{\xi}}_h^n\|_{0,\Sigma}^2. $$

The second term can be controlled via Lemma 4.2 while the first yields the above mentioned $O(\tau/h)$ sub-optimal splitting error.

**Remark 4.4.** As shown in Theorem 4.3, the discrete error estimates of Algorithm 3 contains terms of order $O(\tau^r/h^{1/2})$, which are not visible numerically (see Section 5). To fully understand the impact of selecting the same penalty term in the viscous step as in the solid sub-step in Algorithm 3, we consider the coupling of a parabolic equation with an hyperbolic one. The considered coupled problem reads as follow: find $u: \Omega^f \times \mathbb{R}^+ \to \mathbb{R}$, $d: \Sigma \times \mathbb{R}^+ \to \mathbb{R}$, $\dot{d}: \Sigma \times \mathbb{R}^+ \to \mathbb{R}$, such that for all $t \in \mathbb{R}^+$ we have

$$ \left\{ \begin{array}{ll} \partial_t u - \Delta u = 0 & \text{in } \Omega^f \times \mathbb{R}^+, \\ u = 0 & \text{on } \Gamma \times \mathbb{R}^+, \end{array} \right. \quad \left\{ \begin{array}{ll} \partial_t \dot{d} - \Delta d = T & \text{in } \Sigma \times \mathbb{R}^+, \\ \dot{d} = \partial_t d & \text{in } \Sigma \times \mathbb{R}^+, \\ d = 0 & \text{on } \partial\Sigma \times \mathbb{R}^+, \end{array} \right. \quad \left\{ \begin{array}{ll} u = \dot{d} & \text{on } \Sigma \times \mathbb{R}^+, \\ T = -\frac{\partial u}{\partial n} & \text{on } \Sigma \times \mathbb{R}^+, \end{array} \right. $$

with the respective initial conditions. We propose to discretize the problem via a loosely coupled scheme, inspired by the semi-implicit scheme of Algorithm 3. The fully discrete approximation results in the following (explicit) scheme:

For $n \ge 1$:

* Parabolic step: find $u_h^n \in V_h$, such that

$$ (\partial_{\tau} u_h^n, v_h)_{\Omega^f} + (\nabla u_h^n, \nabla v_h)_{\Omega^f} + \frac{\gamma}{h} (u_h^n - \dot{d}_h^{n-1}, v_h)_{\Sigma} - \left(\frac{\partial u_h^n}{\partial n}, v_h\right)_{\Sigma} = 0 $$

for all $v_h \in V_h$.

* Hyperbolic step: find $(d_h^n, \dot{d}_h^n) \in W_h \times W_h$ with $\dot{d}_h^n = \partial_{\tau} d_h^n$ and such that

$$ (\partial_{\tau} \dot{d}_h^n, w_h)_{\Sigma} + (\nabla d_h^n, \nabla w_h)_{\Sigma} - \frac{\gamma}{h} (u_h^n - \dot{d}_h^{n-1}, w_h)_{\Sigma} + \left(\frac{\partial u_h^n}{\partial n}, w_h\right)_{\Sigma} = 0 $$

for all $w_h \in W_h$.

When considering loosely coupled schemes with Nitsche's coupling, the sub-optimal terms come typically from the fact that we introduce a time-splitting error inside the Nitsche's penalty term, which is scaled with an $h^{-1}$. A possible way to overcome this issue, is to remove the time-splitting error from the Nitsche's penalty term, by introducing an error in time within the definition of the projection errors. Thus, considering the following decomposition of the errors for the parabolic-hyperbolic explicit scheme:

$$ \begin{align} E_2 u^n - u_h^n &= \underbrace{E_2 u^n - i_{sz} E_2 u^n}_{\stackrel{\text{def}}{=} \theta_{\pi}^n} + \underbrace{i_{sz} E_2 u^n - u_h^n}_{\stackrel{\text{def}}{=} \theta_h^n} \quad \text{in } \Omega_h^f, \notag \\ d^n - d_h^n &= \underbrace{d^n - \pi_h d^n}_{\stackrel{\text{def}}{=} \xi_{\pi}^n} + \underbrace{\pi_h d^n - d_h^n}_{\stackrel{\text{def}}{=} \xi_h^n} \quad \text{in } \Sigma, \tag{80} \\ \dot{d}^n - \dot{d}_h^n &= \underbrace{\dot{d}^n - \pi_h \dot{d}^n}_{\stackrel{\text{def}}{=} \dot{\xi}_{\pi}^n} + \underbrace{\pi_h \dot{d}^n - \pi_h \dot{d}^{n+1}}_{\stackrel{\text{def}}{=} \dot{\xi}_{\tau}^n} + \underbrace{\pi_h \dot{d}^{n+1} - \dot{d}_h^n}_{\stackrel{\text{def}}{=} \dot{\xi}_h^n} \quad \text{in } \Sigma, \notag \end{align} $$

it can be proven that the scheme delivers optimal space and time accuracy. More in detail, using similar arguments of the proof of Theorem 4.3, we will get the following terms inside the Nitsche's penalty part:

$$ u^n - u_h^n - (\dot{d}^n - \dot{d}_h^{n-1}) = \theta_{\pi}^n + \theta_h^n - \underbrace{(\dot{d}^n - \pi_h \dot{d}^n)}_{\stackrel{\text{def}}{=} \dot{\xi}_{\pi}^n} + \underbrace{\pi_h \dot{d}^n - \dot{d}_h^{n-1}}_{\stackrel{\text{def}}{=} \dot{\xi}_h^{n-1}}, $$

which does not contain error in time, in fact the arising terms involving $\theta_h^n - \dot{\xi}_h^{n-1}$ are controlled via the stability result and terms involving $\theta_{\pi}^n - \dot{\xi}_{\pi}^n$ have optimal convergence order. The only terms which contain $\dot{\xi}_{\tau}^n$ are the corresponding $T_9$ and $T_{10}$ terms of (59) and their optimality can be proved.

A similar strategy fails when considered for the semi-implicit scheme of Algorithm 3. In particular we will retrieve terms of order $O(\tau/h^{1/2})$ when controlling the pressure term $T_{7,2}$ of (59).
# 5 Numerical experiments

In this section, we illustrate via numerical experiments the convergence properties of Algorithm 3 with $r = 1, 2$ (semi-implicit scheme) in an academic numerical example. The obtained results are compared with those of Algorithms 1 (strongly coupled scheme) and Algorithm 2 (stabilized explicit coupling scheme).

<figure><img src="image_3.png" /><figcaption>Figure 2: Geometric configuration.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Fluid and solid meshes for $h = 0.1$.</figcaption></figure>

The considered test case is the well-known fluid-structure interaction benchmark describing the propagation of a pressure wave within a straight two-dimensional elastic tube (see, e.g., [31, 17, 19]). In the following, all the units are given in the CGS system.

The fluid domain is defined as $\Omega^f = [0, L] \times [0, R]$, with $L = 6$ and $R = 0.5$, while the fluid computational domain is given by $\Omega_h^f = [0, 6] \times [0, 0.8]$. The solid domain is $\Sigma = [0, L] \times \{R\}$, as shown in Figure 2. In the sequel, the solid is described by a one-dimensional string model, viz.,

$$d = \begin{pmatrix} 0 \\ \eta \end{pmatrix}, \quad Ld = \begin{pmatrix} 0 \\ -\lambda_1 \partial_{xx}\eta + \lambda_0 \eta \end{pmatrix}, \quad \lambda_1 \stackrel{\text{def}}{=} \frac{E\epsilon}{2(1+v)}, \quad \lambda_0 \stackrel{\text{def}}{=} \frac{E\epsilon}{R^2(1-\nu^2)},$$

hence in (4), taking $w = (0, w)^T$, we have

$$a^s(d, w) \stackrel{\text{def}}{=} \lambda_1(\partial_x\eta, \partial_x w)_\Sigma + \lambda_0(\eta, w)_\Sigma.$$

The fluid physical parameters are given by $\rho^f = 1.0, \mu = 0.035$. For the solid we have $\rho^s = 1.1$ and $\epsilon = 0.1$, with Young's modulus $E = 0.75 \times 10^6$ and Poisson's ratio $\nu = 0.5$. Regarding the boundary conditions, we consider both fluid and structure to be initially at rest and we impose a sinusoidal normal traction of maximal amplitude $2 \times 10^4$ for $5 \times 10^{-3}$ time instants at the inlet $\Gamma_I$. A symmetry condition is applied on the lower wall $\Gamma_W$ and zero traction is enforced at $\Gamma_O$. All the computations have been performed with FreeFem++ [40]. An example of the unfitted meshes is shown in Figure 3 with fluid space discretization parameter $h = 0.1$. The Nitsche parameter is set to $\gamma = 10^3$ and for the pressure and ghost-penalty stabilization terms (7)-(8) we consider $\gamma_p = 10^{-3}$ and $\gamma_g = 1$, respectively.

As expected, all the considered methods deliver a numerical solution with a stable pressure-wave propagation. For illustration purposes, Figure 4 provides the snapshots of the fluid pressure and solid deformation at time $t = 5 \times 10^{-3}, 10^{-2}$ and $1.5 \times 10^{-2}$, obtained with $\tau = 10^{-4}$ and $h = 0.05$ using respectively Algorithms 1-3. The solid displacement has been amplified by a factor 20. A very good agreement between Algorithm 1 and Algorithm 3 ($r = 1, 2$) is clearly visible, while a difference on the solid displacement is noticeable with Algorithm 2.

In order to quantify the accuracy properties of each coupling scheme we have evaluated the convergence histories by uniformly refining in space and in time

$$(h, \tau) \in \{0.1/2^i, 2 \times 10^{-4}/2^i\}_{i=0}^4. \qquad (81)$$

Figure 5 shows the corresponding solid displacement at $t = 1.5 \times 10^{-2}$ for $i = 0, \dots, 3$ and the different coupling schemes. As in Figure 4, a very good fit is observed between Algorithm 1 and Algorithm 3 ($r = 1, 2$), while a degradation of accuracy is visible for Algorithm 2 under space-time refinement. The
<figure><img src="image_5.png" /><figcaption>Figure 4: Snapshots of the fluid pressure and deformation (magnified) at different time instants.</figcaption></figure>

depicted reference solution has been generated using the strongly coupled fitted method with a high space-time grid resolution ($h = 3.125 \times 10^{-3}$ and $\tau = 10^{-6}$).

Aiming to investigate the space uniformity of the splitting error, we compare Algorithms 1–3 with the variant of the projection-based semi-implicit scheme with second-order backward difference discretization, shown in Algorithm 4. Figure 6 reports the convergence history of the solid displacement at time $t = 1.5 \times 10^{-2}$, in the relative elastic energy-norm. Note that, from the choice of space and time discretization parameters in (81), we have $\tau = O(h)$. The results show that Algorithm 3 with $r = 1, 2$ and Algorithm 1 retrieve the overall optimal first-order accuracy $O(h)$ of Algorithm 1. As expected, Algorithm 2 shows a non-convergence behaviour. This points out the expected $O(\tau/h)$ sub-optimality of the splitting error (see Remark 4.3). Algorithm 4 exhibits better accuracy than Algorithm 3. Full second-order accuracy is not achieved due to the first-order extrapolation of the pressure. Finally, it is worth noting that no effect from the $O(\tau/h^{\frac{1}{2}})$ and $O(\tau^2/h^{\frac{1}{2}})$, anticipated by Theorem 4.3, is visible on the convergence history of Algorithm 3.

# 6 Conclusion

In this paper, an a priori error analysis of the unfitted mesh based semi-implicit coupling scheme introduced in [29] has been performed in the context of linear-fluid structure interaction system involving a thin-walled solid. The considered method combines a Nitsche based unfitted mesh spatial approximation with a fractional-step time-splitting in the fluid. Strong coupling is avoided by treating explicitly the
<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison of the solid displacements at $t = 1.5 \times 10^{-2}$ for different levels of $(\tau, h)$-refinement, given by (81) with $i = 0, \dots, 3$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Convergence history of the solid displacements at $t = 1.5 \times 10^{-2}$ in the elastic-energy norm.</figcaption></figure>
fluid viscous part, while the added-mass free stability (see [29, Theorem 1] and Theorem 4.1 herein) is achieved by treating implicitly the fluid incompressibility and the solid inertia.

The stability analysis of [29] has been extended to the case of second order extrapolation of the solid velocity ($r = 2$) under a similar hyberbolic CFL condition. The proposed error analysis predicts superior accuracy of the semi-implicit scheme of [29] with respect to the stabilized explicit scheme of [19], namely $O(\tau^r/h^{\frac{1}{2}})$, $r = 1, 2$, instead of $O(\tau/h)$. Extensive numerical evidence indicated that the semi-implicit algorithm and the strongly coupled (from [19]) deliver practically the same accuracy. In particular, no effect from the non-uniformity in space is visible, for both the first and second-order backward difference variants.

# Acknowledgment

This work has been partially supported by Inria through the IMFIBIO associated team.

# References

[1] ANNESE, M., FERNÁNDEZ, M. A., AND GASTALDI, L. Splitting schemes for a lagrange multiplier formulation of fsi with immersed thin-walled structure: stability and convergence analysis. (2020), arXiv:2007.04699.

[2] ANTONIETTI, P. AND VERANI, M. AND VERGARA, C. AND ZONCA, S. Numerical solution of fluid-structure interaction problems by means of a high order Discontinuous Galerkin method on polygonal grids. *Finite Elements in Analysis and Design* **159** (2019), 1–14.

[3] ASTORINO, M., CHOULY, F., AND FERNÁNDEZ, M. A. Robin based semi-implicit coupling in fluid-structure interaction: Stability analysis and numerics. *SIAM J. Sci. Comput.* **31**, 6 (2009), 4041–4065.

[4] ASTORINO, M., AND GRANDMONT, C. Convergence analysis of a projection semi-implicit coupling scheme for fluid-structure interaction problems. *Numer. Math.* **116** (2010), 721–767.

[5] BADIA, S., QUAINI, A., AND QUARTERONI, A. Splitting methods based on algebraic factorization for fluid-structure interaction. *SIAM J. Sci. Comput.* **30**, 4 (2008), 1778–1805.

[6] BANKS, J. W. AND HENSHAW, W. D. AND SCHWENDEMAN, D. W. An analysis of a new stable partitioned algorithm for FSI problems. Part I: Incompressible flow and elastic solids. *Journal of Computational Physics* **269** (2014), 108–137.

[7] BECKER, R., BURMAN, E., AND HANSBO, P. A nitsche extended finite element method for incompressible elasticity with discontinuous modulus of elasticity. *Computer Methods in Applied Mechanics and Engineering* **198** (09 2009).

[8] BOFFI, D., CAVALLINI, N., AND GASTALDI, L. Finite element approach to immersed boundary method with different fluid and solid densities. *Math. Models Methods Appl. Sci.* **21**, 12 (2011), 2523–2550.

[9] BOFFI, D., AND GASTALDI, L. A fictitious domain approach with lagrange multiplier for fluid-structure interactions. *Numer. Math.* **135**, 3 (2017), 711–732.

[10] BOILEVIN-KAYL, L., FERNÁNDEZ, M. A., AND GERBEAU, J.-F. A loosely coupled scheme for fictitious domain approximations of fluid-structure interaction problems with immersed thin-walled structures. *SIAM J. Sci. Comput.* **41**, 2 (2019), B351–B374.

[11] BREZZI, F., AND PITKÄRANTA, J. On the stabilization of finite element approximations of the Stokes equations. In *Efficient solutions of elliptic systems* (Kiel, 1984), vol. 10 of *Notes Numer. Fluid Mech*. Vieweg, 1984, pp. 11–19.

[12] BUKAC, M., CANIC, C., GLOWINSKI, R., TAMBACA, T., AND QUAINI, A. Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement. *J. Comp. Phys.* **235**, 0 (2013), 515–541.

[13] BUKAC, M., AND MUHA, B. Stability and convergence analysis of the extensions of the kinematically coupled scheme for the fluid-structure interaction. *SIAM J. Numer. Anal.* **54**, 5 (2016), 3032–3061.

[14] BUKAČ, MARTINA AND ČANIĆ, SUNČICA. A partitioned numerical scheme for fluid-structure interaction with slip. *Mathematical Modelling of Natural Phenomena* **16** (2021), 8.
[15] BURMAN, E. Ghost penalty. *C. R. Math. Acad. Sci. Paris 348*, 21-22 (2010), 1217-1220.

[16] BURMAN, E., ERN, A., AND FERNÁNDEZ, M. A. Fractional-step methods and finite elements with symmetric stabilization for the transient Oseen problem. *ESAIM Math. Model. Numer. Anal. 51*, 2 (2017), 487-507.

[17] BURMAN, E., AND FERNÁNDEZ, M. Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility. *Comput. Methods Appl. Mech. Engrg. 198*, 5-8 (2009), 766-784.

[18] BURMAN, E., AND FERNÁNDEZ, M. Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin-Robin coupling. *Int. J. Num. Meth. Engrg. 97*, 10 (2014), 739-758.

[19] BURMAN, E., AND FERNÁNDEZ, M. An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes. *Comput. Methods Appl. Mech. Engrg. 279* (2014), 497-514.

[20] BURMAN, E., AND HANSBO, P. Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method. *Appl. Numer. Math. 62*, 4 (2012), 328-341.

[21] BURMAN, E., AND HANSBO, P. Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes' problem. *ESAIM Math. Model. Numer. Anal. 48*, 3 (2014), 859-874.

[22] CAUSIN, P., GERBEAU, J.-F., AND NOBILE, F. Added-mass effect in the design of partitioned algorithms for fluid-structure problems. *Comput. Methods Appl. Mech. Engrg. 194*, 42-44 (2005), 4506-4527.

[23] DU, Q., GUNZBURGER, M., HOU, L., AND LEE, J. Semidiscrete finite element approximations of a linear fluid-structure interaction problem. *SIAM J. Numer. Anal. 42*, 1 (2004), 1-29.

[24] ERN, A., AND GUERMOND, J.-L. *Theory and Practice of Finite Elements*, vol. 159 of *Applied Mathematical Sciences*. Springer, New York, 2004.

[25] EVANS, L. C. Partial differential equations, ams. *Graduate Studies in Mathematics 19* (2002).

[26] FERNÁNDEZ, M. Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis. *Numer. Math. 123*, 1 (2013), 21-65.

[27] FERNÁNDEZ, M., GERBEAU, J., AND GRANDMONT, C. A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid. *Int. J. Num. Meth. Engrg. 69*, 4 (2007), 794-821.

[28] FERNÁNDEZ, M., AND MULLAERT, J. Convergence and error analysis for a class of splitting schemes in incompressible fluid-structure interaction. *IMA J. Numer. Anal. 36*, 4 (2016), 1748-1782.

[29] FERNÁNDEZ, M. A., AND GEROSA, F. An unfitted mesh semi-implicit coupling scheme for fluid-structure interaction with immersed solids. *International Journal for Numerical Methods in Engineering* (2020). DOI: 10.1002/nme.6449.

[30] FERNÁNDEZ, M. A., AND LANDAJUELA, M. Splitting schemes and unfitted-mesh methods for the coupling of an incompressible fluid with a thin-walled structure. *IMA J. Numer. Anal. 40*, 2 (2020), 1407-1453.

[31] FORMAGGIA, L., GERBEAU, J., NOBILE, F., AND QUARTERONI, A. On the coupling of 3d and 1d navier-stokes equations for flow problems in compliant vessels. *Computer Methods in Applied Mechanics and Engineering 191* (12 2001), 561-582.

[32] FÖRSTER, C., WALL, W., AND RAMM, E. Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows. *Comput. Methods Appl. Mech. Engrg. 196*, 7 (2007), 1278-1293.

[33] GERSTENBERGER, A., AND WALL, W. An extended finite element method/Lagrange multiplier based approach for fluid-structure interaction. *Comput. Methods Appl. Mech. Engrg. 197*, 19-20 (2008), 1699-1714.

[34] GIGANTE, G. AND VERGARA, C. On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction. *Computers & Mathematics with Applications 96* (2021), 109-119.

[35] GUERMOND, J. L., MINEV, P., AND SHEN, J. An overview of projection methods for incompressible flows. *Comput. Methods Appl. Mech. Engrg. 195*, 44-47 (2006), 6011-6045.

[36] GUIDOBONI, G. AND GLOWINSKI, R. AND CAVALLINI, N. AND CANIC, S. Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow. *Journal of Computational Physics 228*, 18 (2009), 6916-6937.
[37] GUO, R. Solving parabolic moving interface problems with dynamical immersed spaces on unfitted meshes: Fully discrete analysis. *SIAM Journal on Numerical Analysis* 59, 2 (2021), 797–828.

[38] HANSBO, A., AND HANSBO, P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. *Comput. Methods Appl. Mech. Engrg.* 193, 33-35 (2004), 3523–3540.

[39] HANSBO, P., HERMANSSON, J., AND SVEDBERG, T. Nitsche's method combined with space-time finite elements for ALE fluid-structure interaction problems. *Comput. Methods Appl. Mech. Engrg.* 193, 39-41 (2004), 4195–4206.

[40] HECHT, F. New development in freefem++. *J. Numer. Math.* 20, 3-4 (2012), 251–265.

[41] HEYWOOD, J. G., AND RANNACHER, R. Finite-element approximation of the nonstationary navier–stokes problem. part iv: Error analysis for second-order time discretization. *SIAM Journal on Numerical Analysis* 27, 2 (1990), 353–384.

[42] KADAPA, C., DETTMER, W., AND PERIĆ, D. A stabilised immersed framework on hierarchical b-spline grids for fluid-flexible structure interaction with solid-solid contact. *Comput. Methods Appl. Mech. Engrg.* 335 (2018), 472–489.

[43] KIM, W., AND LEE, I. CHOI, H. A weak-coupling immersed boundary method for fluid-structure interaction with low density ratio of solid to fluid. *J. Comput. Phys.* 359 (2018), 296–311.

[44] LE TALLEC, P., AND MOURO, J. Fluid structure interaction with large structural displacements. *Comput. Meth. Appl. Mech. Engrg.* 190 (2001), 3039–3067.

[45] LEHRENFELD, C. AND OLSHANSKII, M. An Eulerian finite element method for PDEs in time-dependent domains. *ESAIM: Mathematical Modelling and Numerical Analysis* 53, 2 (2019), 585–614.

[46] LUKÁČOVÁ-MEDVID'OVÁ, M. AND RUSNÁKOVÁ, G. AND HUNDERTMARK-ZAUŠKOVÁ, A. Kinematic splitting algorithm for fluid-structure interaction in hemodynamics. *Computer Methods in Applied Mechanics and Engineering* 265 (2013), 83–106.

[47] MASSING, ANDRÉ AND LARSON, MATS G. AND LOGG, ANDERS. Efficient implementation of finite element methods on nonmatching and overlapping meshes in three dimensions. *SIAM J. Sci. Comput.* 35, 1 (2013), C23–C47.

[48] PESKIN, C. The immersed boundary method. *Acta Numer.* 11 (2002), 479–517.

[49] QUAINI, A., AND QUARTERONI, A. A semi-implicit approach for fluid-structure interaction based on an algebraic fractional step method. *Math. Models Methods Appl. Sci.* 17, 6 (2007), 957–983.

[50] SALSA, S. *Partial Differential Equations in Action: From Modelling to Theory*. Springer, Milan, 2009.

[51] SAWADA, T., AND TEZUKA, A. LLM and X-FEM based interface modeling of fluid-thin structure interactions on a non-interface-fitted mesh. *Comput. Mech.* 48, 3 (2011), 319–332.

[52] SCHOTT, B. AND AGER, C. AND WALL, W. A. Monolithic cut finite element-based approaches for fluid-structure interaction. *International Journal for Numerical Methods in Engineering* 119, 8 (2019), 757–796.

[53] VAN BRUMMELEN, E. Added mass effects of compressible and incompressible flows in fluid-structure interaction. *J. Appl. Mech.* 76, 2 (2009), 021206–7.

[54] ZUNINO, P. Analysis of backward Euler/extended finite element discretization of parabolic problems with moving interfaces. *Computer Methods in Applied Mechanics and Engineering* 258 (2013), 152–165.