# Robin-Robin loose coupling for incompressible fluid-structure interaction: non-linear setting and nearly-optimal error analysis

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Oscar Ruz

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

## To cite this version:

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Oscar Ruz. Robin-Robin loose coupling for incompressible fluid-structure interaction: non-linear setting and nearly-optimal error analysis. Journal of Scientific Computing, 2025, 104 (34), 10.1007/s10915-025-02939-5. hal-04258861v3

HAL Id: hal-04258861

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

Submitted on 1 Feb 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>
Robin-Robin loose coupling for
incompressible fluid-structure interaction:
non-linear setting and nearly-optimal error
analysis

Erik Burman¹, Rebecca Durst², Miguel A. Fernández³*, Johnny Guzmán² and Oscar Ruz³

¹Department of Mathematics, University College London,
London, United Kingdom.

²Division of Applied Mathematics, Brown University, Providence,
USA.

³Sorbonne Université, Inria & CNRS, UMR 7598 Laboratoire
Jacques-Louis Lions, Paris, France.

*Corresponding author(s). E-mail(s): miguel.fernandez@inria.fr;

Contributing authors: e.burman@ucl.ac.uk;

rebecca_durst@alumni.brown.edu; johnny_guzman@brown.edu;

oscar.ruz@inria.fr;

Abstract

The numerical simulation of incompressible fluid-structure interaction
systems with loosely coupled schemes is a delicate problem. Indeed,
the splitting method must both be stable for the full nonlinear sys-
tem and have sufficient accuracy to be of use in practice. In the
case of the coupling of an incompressible fluid with thick-walled
solids the error analyses reported in the literature are limited to
O(τ½) accuracy in time, where τ denotes the time step. The ob-
jective of the present work is to show two important extensions of
the analysis of the Robin-Robin loosely coupled scheme recently re-
ported in [*Numer. Math.*, **151**(4):807-840, 2022]. First, we give a
formulation of the scheme in a general non-linear setting and prove
its unconditional energy stability. Then we show that nearly-optimal
2 Robin-Robin loose coupling for FSI

$O(\tau\sqrt{1 + \log \tau^{-1}})$ accuracy can be achieved in the linear case. These
theoretical findings are illustrated in a series of numerical examples.

**Keywords:** Fluid-structure interaction, loosely coupled scheme, moving domains, Robin conditions, optimal accuracy

# 1 Introduction

The numerical approximation of the interaction of an elastic solid with an in-
compressible viscous fluid over a moving interface (FSI) remains a challenging
problem in scientific computing. For a global overview of the framework and
different approaches with special focus on cardiovascular flow which is a rel-
evant application for the discussion herein, we refer to [1, 2]. Typically the
problem is nonlinear both through the separate subproblems and the mov-
ing interface. The coupling is very stiff due to the different time scales of the
subproblems. In particular, the incompressibility and the related compatibil-
ity conditions present difficulties. For computational efficiency it is appealing
to use a loosely coupled scheme to advance the coupled problem in time. In
such an approach one solve is performed in each subdomain and the coupling
data is tranferred across the interface in a staggered manner. Although this is
feasible in elasto-acoustic, where the fluid is compressible (see, e.g., [3]), it is
much more delicate in the incompressible regime due to the added-mass effect
(see [4, 5]) unless special care is taken in the design of the coupling conditions.

The first stable loosely coupled scheme for FSI in the presence of incom-
pressibility was proposed in [6], using Nitsche's method for FSI (see [7]) in
combination with a pressure stabilization that relaxes the incompressibility on
the interface. For this type of method, however, the stiffness of the Nitsche
penalty term introduces the parabolic condition $\tau \lesssim h^2$ for overall first-order
accuracy (but not for stability), where $\tau$ stands for the time-step length and $h$
for the spatial grid parameter. It was then noted in [8] that the Nitsche cou-
pling method could be transformed into Robin type coupling scheme without
loss of stability. Robin conditions had been shown to perform well in the par-
titioned iterative solution of strong coupling in [9, 10] (see also [11]). A number
of studies have since then been reported in the literature on Robin-like loose
coupling for incompressible FSI involving thick-walled solids (see, e.g., [12–
18]), but without rigorous error analysis nor convergence rates. In [19], it is
shown that the generalized Robin-Neumann scheme of [13] introduces a split-
ting error scaling as $O(\tau/h^{3/2})$, which requires a $\tau \lesssim h^{3/2}$ condition for accuracy
(not for stability).

In [8], the Robin coefficient inherited the $O(h^{-1})$ scaling from Nitsche's
method, but it was noted numerically that this scaling appeared not to be
necessary for the stability of the algorithm. In a recent development [20], it
was shown for a semi-discretization in time that the coupling scheme indeed
is stable independent of the added mass effect and an error estimate of order $O(\tau^{\frac{1}{2}})$, where $\tau$ denotes the time step, was shown. Another result in this spirit was reported in [21], for a fully discrete Robin-Robin method and under a CFL-like condition on the discretization parameters. A stability estimate for the time semi-discretized problem with moving fluid domain is also provided, but the treatment of the convective term involves an artificial modification of the fluid velocity on the interface. An alternative fully discrete formulation is introduced in [22], where the stresses are represented in terms of Lagrange multipliers or variational residuals, and the exponential growth of the stability factor was moderated to sub-linear growth without conditions on the discretization parameters. Finally, both for the cases of parabolic-parabolic and parabolic-hyperbolic couplings, it was shown in [23] that in most situations one can expect the nearly-optimal order $O(\tau\sqrt{1+\log\tau^{-1}})$. Depending on the physical parameters however it was noted that this may not be attained, except for very small time steps unless the parameter in the Robin condition is chosen carefully.

The objective of the present paper is to continue our work on the Robin method introduced in [8] by revisiting the results of [22, 23] and extending them to more realistic settings. In particular we show that the Robin-Robin loose coupling is indeed expected to result in a first order time discretization scheme also in physically relevant fluid-structure interaction problems. The scientific contributions of this work are threefold. First, in a fully non-linear setting (coupling viscous incompressible flow in moving domains with the dynamics of an hyper-elastic solid), we provide a rigorous fully discrete formulation with suitable stabilization terms which guarantee unconditional energy stability. Second, in the linear case, we present the extension of the nearly-optimal error analysis from [23] to the case of the Stokes-elasticity system. To the best of our knowledge, this is the first time that such accuracy is proved for a loosely coupled scheme in incompressible fluid-structure interaction involving thick-walled solids. The latter results were first reported in [24, Chapter 5]. Finally, the accuracy properties of the proposed method are illustrated in a series of well-known numerical examples. We in particular investigate the dependence of the accuracy on the physical parameters and the choice of the free parameter in the Robin coupling.

The rest of the paper is organized as follows. In Section 2, we introduce the considered physically realistic setting. We show how to extend the fully discrete formulation of [22] to this nonlinear case and prove the fundamental stability estimate in Section 3. Section 4 is devoted to the derivation of the nearly-optimal error estimate in a simplified linear setting. Finally, the paper is ended with some numerical examples which are reported in Section 5.
4 Robin-Robin loose coupling for FSI

## 2 Problem setting

We consider a fluid-structure interaction system in which the fluid is modeled by the incompressible Navier-Stokes equations, in arbitrary Lagrangian-Eulerian (ALE) formalism, and the solid by the non-linear elastodynamics equations. The reference configuration of system is given by the domain $\Omega \stackrel{\text{def}}{=} \Omega^f \cup \Omega^s \subset \mathbb{R}^d$, with $d = 2, 3$. For all time $t > 0$, the current solid configuration is denoted by $\Omega^s(t) \subset \mathbb{R}^d$, whereas $\Omega^f(t) \subset \mathbb{R}^d$ stands for the moving fluid control volume. We denote by $\Sigma(t) \stackrel{\text{def}}{=} \partial\Omega^f(t) \cap \partial\Omega^s(t)$ the current configuration of the fluid-structure interface, whose reference configurations corresponds to $\Sigma \stackrel{\text{def}}{=} \partial\Omega^f \cap \partial\Omega^s$. The remaining parts of the fluid and solid boundaries are supposed to be time independent and partitioned as $\partial\Omega^f = \Gamma^f \cup \Sigma$ and $\partial\Omega^s = \Gamma^s \cup \Sigma$.

The moving fluid domain $\Omega^f(t)$ can be parametrized in terms of the so-called ALE map, $\mathcal{A}: \Omega^f \times \mathbb{R}^+ \to \mathbb{R}^d$, in such a way that $\Omega^f(t) = \mathcal{A}(\Omega^f, t)$. The ALE deformation map $\mathcal{A}$ is often written in terms of fluid domain displacement $d^f: \Omega^f \times \mathbb{R}^+ \to \mathbb{R}^d$ via the relation $\mathcal{A} = I_{\Omega^f} + d^f$, where $I_{\Omega^f}$ stands for the identity operator in $\Omega^f$. We can hence introduce the notation $\mathcal{A}_t \stackrel{\text{def}}{=} \mathcal{A}(\cdot, t)$, the fluid domain velocity $w \stackrel{\text{def}}{=} \partial_t \mathcal{A} = \partial_t d^f$, the fluid gradient of deformation $\mathbf{F} \stackrel{\text{def}}{=} \nabla \mathcal{A}$ (spatial gradient) and its Jacobian $J \stackrel{\text{def}}{=} \det \mathbf{F}$. For a given field $f$ defined in the current configuration $\Omega^f(t)$, we shall use the notation $\hat{f}$ to defined its corresponding ALE description in $\Omega^f$, by composing with the ALE map, namely,

$$ \hat{f}(x, t) \stackrel{\text{def}}{=} f(\mathcal{A}(x, t), t) \quad \forall x \in \Omega^f. $$

Conversely, a given field defined in the reference configuration $\Omega^f$ can also be transported into the current configuration by composition with $\mathcal{A}_t^{-1}$. For the sake of simplicity, this composition is omitted. For instance, $w$ in (1)$_1$ below has to be read as $w \circ \mathcal{A}_t^{-1}$.

## Coupled problem in strong form

The considered fluid-structure interaction model reads as follows: find the fluid domain displacement $d^f: \Omega^f \times \mathbb{R}^+ \to \mathbb{R}^d$, the fluid velocity $\hat{u}: \Omega^f \times \mathbb{R}^+ \to \mathbb{R}^d$, the fluid pressure $\hat{p}: \Omega^f \times \mathbb{R}^+ \to \mathbb{R}$, the solid displacement $d: \Omega^s \times \mathbb{R}^+ \to \mathbb{R}^d$, and the structure velocity $\dot{d}: \Omega^s \times \mathbb{R}^+ \to \mathbb{R}^d$ such that, for all $t > 0$,

$$ \left\{ \begin{array}{ll} \rho^f \partial_t u|_{\mathcal{A}} + \rho^f(u-w) \cdot \nabla u - \operatorname{div} \sigma(u,p) = 0 & \text{in } \Omega^f(t), \\ \operatorname{div} u = 0 & \text{in } \Omega^f(t), \\ u = 0 & \text{on } \Gamma^f, \end{array} \right. \quad (1) $$

$$ \left\{ \begin{array}{ll} \rho^s \partial_t \dot{d} - \operatorname{div} (\mathbf{F}^s \boldsymbol{\Sigma}(d)) = 0 & \text{in } \Omega^s, \\ \dot{d} = \partial_t d & \text{in } \Omega^s, \\ d = 0 & \text{on } \Gamma^s, \end{array} \right. \quad (2) $$
$$
\left\{
\begin{array}{l}
\displaystyle \mathbf{d}^f = \mathcal{L}(\mathbf{d}|_{\Sigma}), \quad \mathbf{w} = \partial_t \mathbf{d}^f, \quad \mathcal{A} = \mathbf{I}_{\Omega^f} + \mathbf{d}^f, \quad \Omega^f(t) = \mathcal{A}(\Omega^f, t), \\
\displaystyle \hat{\mathbf{u}} = \dot{\mathbf{d}} \quad \text{on } \Sigma, \qquad (3) \\
\displaystyle \mathbf{F}^s \Sigma(\mathbf{d}) \mathbf{n}^s = -J \hat{\sigma}(\mathbf{u}, p) \mathbf{F}^{-T} \hat{\mathbf{n}} \quad \text{on } \Sigma.
\end{array}
\right.
$$

The constants $\rho^f, \rho^s > 0$ stand, respectively, for the fluid and solid densities, $\partial_t|_{\mathcal{A}} \stackrel{\text{def}}{=} \partial_t + \mathbf{w} \cdot \nabla$ for the ALE time derivative and $\sigma(\mathbf{u}, p) \stackrel{\text{def}}{=} 2\mu\epsilon(\mathbf{u}) - p\mathbf{I}$ for the fluid Cauchy stress tensor, with $\mu > 0$ the fluid dynamic viscosity (supposed constant) and $\epsilon(\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + (\nabla\mathbf{u})^T)$ the strain rate tensor. The symbol $\mathbf{F}^s \stackrel{\text{def}}{=} \mathbf{I} + \nabla\mathbf{d}$ corresponds to the gradient of deformation of the solid and $\Sigma(\mathbf{d})$ denotes its second Piola-Kirchoff stress tensor which, for an hyper-elastic material, is given in terms of the internal energy functional $W : \mathbb{R}_{\text{sym}}^{d \times d} \rightarrow \mathbb{R}^+$ as

$$
\Sigma = \frac{\partial W}{\partial E}(E),
$$

where the symbol $E$ denotes the Green-Lagrange strain tensor given by the relation $E \stackrel{\text{def}}{=} \frac{1}{2}((F^s)^T F^s - \mathbf{I})$. The relation (3)$_1$ enforces the geometrical compatibility between the fluid and solid domains. Here, the symbol $\mathcal{L}$ represents a smooth lifting operator (possibly non-linear) from $\Sigma$ over $\Omega^f$, which vanishes on $\Gamma^f$. The remaining interface conditions (3)$_{2,3}$ enforce, respectively, the so-called interface kinematic and dynamic coupling.

Finally, the coupled system (1)-(3) is complemented with the following initial conditions: $\hat{\mathbf{u}}(0) = \hat{\mathbf{u}}_0$ in $\Omega^f$ and $\mathbf{d}(0) = \mathbf{d}_0$, $\dot{\mathbf{d}}(0) = \dot{\mathbf{d}}_0$ in $\Omega^s$.

Energy identity and weak form

By taking the scalar product of (1)$_1$ with $\mathbf{u}$ and of (2)$_1$ with $\dot{\mathbf{d}}$, integrating by parts, using the boundary conditions and by noting that the interface integrals cancel due to the coupling conditions (3), we get the following classical energy identity (see, e.g., [1]):

$$
\begin{align}
& \frac{\rho^f}{2} \|\mathbf{u}\|_{0,\Omega^f(t)}^2 + \frac{\rho^s}{2} \|\dot{\mathbf{d}}\|_{0,\Omega^s}^2 + \int_{\Omega^s} W(\mathbf{E}(\mathbf{d})) + 2\mu \int_0^t \|\epsilon(\mathbf{u}(s))\|_{0,\Omega^f(s)}^2 \mathrm{d}s \nonumber \\
& = \frac{\rho^f}{2} \|\mathbf{u}_0\|_{0,\Omega^f(0)}^2 + \frac{\rho^s}{2} \|\dot{\mathbf{d}}_0\|_{0,\Omega^s}^2 + \int_{\Omega^s} W(\mathbf{E}(\mathbf{d}_0)). \tag{4}
\end{align}
$$

Here, the first two terms represent the kinetic energy of the fluid and solid, respectively. The third term takes into account the hyperelastic energy of the solid, and the fourth term corresponds to the viscous dissipation of the fluid.

We consider the following functional spaces:

$$
\begin{align*}
\mathbf{V}^f &\subset \{ \mathbf{v} \in \mathbf{H}^1(\Omega^f) : \mathbf{v}|_{\Gamma^f} = \mathbf{0} \}, & M^f &\stackrel{\text{def}}{=} L^2(\Omega^f), \\
\mathbf{V}^s &\subset \{ \mathbf{v} \in \mathbf{H}^1(\Omega^s) : \mathbf{v}|_{\Gamma^s} = \mathbf{0} \}.
\end{align*}
$$
6 Robin-Robin loose coupling for FSI

Let $(\hat{u}, \hat{p}, d, \dot{d})$ be solution of the coupled system (1)-(3), one can show that the following monolithic variational formulation holds (see, e.g., [1]), for $t > 0$,

$$ \rho^f \frac{d}{dt} \int_{\Omega^f(t)} \mathbf{u} \cdot \mathbf{v} - \rho^f \int_{\Omega^f(t)} (\mathrm{div} \, \mathbf{w}) \mathbf{u} \cdot \mathbf{v} + \rho^f \int_{\Omega^f(t)} (\mathbf{u} - \mathbf{w}) \cdot \nabla \mathbf{u} \cdot \mathbf{v} + 2\mu \int_{\Omega^f(t)} \epsilon(\mathbf{u}) : \epsilon(\mathbf{v}) \\ - \int_{\Omega^f(t)} p \mathrm{div} \, \mathbf{v} + \int_{\Omega^f(t)} q \mathrm{div} \, \mathbf{u} + \rho^s \int_{\Omega^s} \partial_t \dot{\mathbf{d}} \cdot \boldsymbol{\xi} + \int_{\Omega^s} \boldsymbol{\Sigma} : \partial_d(\mathbf{E})\boldsymbol{\xi} = 0 $$

for all $(\hat{v}, \hat{q}, \boldsymbol{\xi}) \in \mathbf{V}^f \times \mathbf{M}^f \times \mathbf{V}^s$ with $\hat{v}|_{\Sigma} = \boldsymbol{\xi}|_{\Sigma}$, and where we have used the notation

$$ \partial_d(\mathbf{E})\boldsymbol{\xi} \stackrel{\mathrm{def}}{=} \frac{1}{2} [ (\mathbf{F}^s)^{\mathrm{T}} \nabla \boldsymbol{\xi} + (\nabla \boldsymbol{\xi})^{\mathrm{T}} \mathbf{F}^s ], \quad (5) $$

for the differential of $\mathbf{E}$ applied to $\boldsymbol{\xi}$.

## 3 Robin-Robin explicit coupling scheme

The purpose of this section is to formulate the Robin-Robin loosely coupled scheme of [20, 22] in the framework of the non-linear coupled problem (1)-(3).

### Time semi-discretization

In what follows, the scalar $\tau > 0$ denotes the time-step length and $\{t_n\}_{n \in \mathbb{N}}$ represents the temporal grid. We shall also make use of the following notations:

$$ \partial_{\tau} x^{n} \stackrel{\mathrm{def}}{=} \frac{1}{\tau}(x^{n} - x^{n-1}), \quad x^{n-\frac{1}{2}} \stackrel{\mathrm{def}}{=} \frac{1}{2}(x^{n} + x^{n-1}), $$

for the first-order backward difference and the mid-point value, respectively. Moreover, for a given time-dependent quantity $x(t)$, we will use the notation $x^n \stackrel{\mathrm{def}}{=} x(t_n)$.

The bulk terms of the fluid sub-system (1) will be discretized in time with a backward-Euler semi-implicit scheme, whereas for the solid sub-system (2) we consider a mid-point scheme with the so-called Gonzalez correction (see [25]). Finally, for the time discretization of the interface conditions (3), we combine an explicit treatment of the geometrical compatibility with a Robin-Robin type splitting of the kinematic and dynamic interface conditions, originally proposed in [20] for the linear case. By introducing the following notation for the interfacial fluid stresses on the reference configuration,

$$ \boldsymbol{\lambda}^n \stackrel{\mathrm{def}}{=} J^n \hat{\boldsymbol{\sigma}}(\mathbf{u}^n, p^n) (\mathbf{F}^n)^{-\mathrm{T}} \hat{\mathbf{n}} \quad \text{on } \Sigma, \quad (6) $$

the resulting time-stepping procedure can formally be written as follows:
• Solve solid:

$$
\left\{
\begin{array}{ll}
\rho^s \partial_\tau \dot{\mathbf{d}}^n - \mathrm{div} \left( \mathbf{F}^{s,n-\frac{1}{2}} \bar{\mathbf{\Sigma}}^{n-\frac{1}{2}} \right) = \mathbf{0} & \text{in } \Omega^s, \\
\dot{\mathbf{d}}^{n-\frac{1}{2}} = \partial_\tau \mathbf{d}^n & \text{in } \Omega^s, \\
\mathbf{d}^n = \mathbf{0} & \text{on } \Gamma^s, \\
\mathbf{F}^{s,n-\frac{1}{2}} \bar{\mathbf{\Sigma}}^{n-\frac{1}{2}} \mathbf{n}^s + \alpha (\dot{\mathbf{d}}^{n-\frac{1}{2}} - \hat{\mathbf{u}}^{n-1}) = -\lambda^{n-1} & \text{on } \Sigma.
\end{array}
\right.
\tag{7}
$$

• Fluid domain update:

$$
\mathbf{d}^{f,n} = \mathcal{L}(\mathbf{d}^{n-1}|_{\Sigma}), \quad \mathbf{w}^n = \partial_{\tau} \mathbf{d}^{f,n}, \quad \mathcal{A}^n = \mathbf{I}_{\Omega^{f}} + \mathbf{d}^{f,n}, \quad \Omega^{f,n} = \mathcal{A}^n(\Omega^{f}). \tag{8}
$$

• Solve fluid:

$$
\left\{
\begin{array}{ll}
\rho^f \partial_\tau \mathbf{u}^n|_{\mathcal{A}} + \rho^f (\mathbf{u}^{n-1} - \mathbf{w}^n) \cdot \nabla \mathbf{u}^n - \mathrm{div} \, \boldsymbol{\sigma}(\mathbf{u}^n, p^n) = \mathbf{0} & \text{in } \Omega^{f,n}, \\
\mathrm{div} \, \mathbf{u}^n = 0 & \text{in } \Omega^{f,n}, \\
\mathbf{u}^n = \mathbf{0} & \text{on } \Gamma^f, \\
\boldsymbol{\lambda}^n + \alpha (\hat{\mathbf{u}}^n - \dot{\mathbf{d}}^{n-\frac{1}{2}}) = \boldsymbol{\lambda}^{n-1} & \text{on } \Sigma.
\end{array}
\right.
\tag{9}
$$

Here, the symbol $\bar{\Sigma}^{n-\frac{1}{2}}$ denotes the mid-point correction of the solid stress given by (see [25]):

$$
\bar{\Sigma}^{n-\frac{1}{2}} \stackrel{\text{def}}{=} \Sigma^{n-\frac{1}{2}} + \frac{\partial_{\tau} W(E^n) - \Sigma^{n-\frac{1}{2}} : \partial_{\tau} E^n}{|\partial_{\tau} E^n|^2} \partial_{\tau} E^n, \quad (10)
$$

which guarantees energy stability in the solid. The mid-point stress $\bar{\Sigma}^{n-\frac{1}{2}}$ can
be taken as the standard half-sum or (with a little abuse of notation) as

$$
\Sigma^{n-\frac{1}{2}} \stackrel{\text{def}}{=} \frac{\partial W}{\partial \mathbf{E}}(\mathbf{E}^{n-\frac{1}{2}}), \quad \mathbf{E}^{n-\frac{1}{2}} \stackrel{\text{def}}{=} \mathbf{E}(\mathbf{d}^{n-\frac{1}{2}}),
$$

which is often preferred for implementation purposes (as in this paper). Finally,
the scalar $\alpha > 0$ stands for the user-defined Robin parameter of the fluid-solid
splitting. Note that the above interfacial Robin-Robin splitting, viz.,

$$
\left\{
\begin{array}{ll}
\boldsymbol{F}^{s,n-\frac{1}{2}} \bar{\boldsymbol{\Sigma}}^{n-\frac{1}{2}} \boldsymbol{n}^s + \alpha (\dot{\boldsymbol{d}}^{n-\frac{1}{2}} - \hat{\boldsymbol{u}}^{n-1}) = -\boldsymbol{\lambda}^{n-1} & \text{on } \Sigma, \\
\boldsymbol{\lambda}^n + \alpha (\hat{\boldsymbol{u}}^n - \dot{\boldsymbol{d}}^{n-\frac{1}{2}}) = \boldsymbol{\lambda}^{n-1} & \text{on } \Sigma,
\end{array}
\right.
\tag{11}
$$

is written in the reference configuration of the interface.
8 Robin-Robin loose coupling for FSI

Space discretization, fully discrete scheme

For the spatial discretization of (7)–(9), we adopt a finite element approximation. Let $\mathcal{T}_h^f$ and $\mathcal{T}_h^s$ be, respectively, simplicial triangulations of $\Omega^f$ and $\Omega^s$. We assume that $\mathcal{T}_h^f$ and $\mathcal{T}_h^s$ match on the interface $\Sigma$. We then define the following finite element spaces:

$$ 
\begin{align*}
\mathbf{V}_h^s &\stackrel{\text{def}}{=} \{\mathbf{v} \in \mathbf{V}^s : \mathbf{v}|_K \in \mathbb{P}^1(K), \forall K \in \mathcal{T}_h^s\}, \\
\mathbf{V}_h^f &\stackrel{\text{def}}{=} \{\mathbf{v} \in \mathbf{V}^f : \mathbf{v}|_K \in \mathbb{P}^1(K), \forall K \in \mathcal{T}_h^f\}, \\
M_h^f &\stackrel{\text{def}}{=} \{v \in M^f : v|_K \in \mathbb{P}^1(K), \forall K \in \mathcal{T}_h^f\}.
\end{align*}
$$

Since the fluid and solid meshes are assumed to match at the interface, we can define the discrete traces space

$$ \mathbf{V}_h^g \stackrel{\text{def}}{=} \text{Tr}_{\Sigma} \mathbf{V}_h^f = \text{Tr}_{\Sigma} \mathbf{V}_h^s. \qquad (12) $$

The fluid bi-linear form is given by

$$ 
\begin{align}
a_{\Omega^{f,n}}(\mathbf{u}_h^{n-1}, \mathbf{w}_h^n; (\mathbf{u}_h^n, p_h^n), (\mathbf{v}_h, q_h)) &\stackrel{\text{def}}{=} \rho^f \int_{\Omega^{f,n}} (\mathbf{u}_h^{n-1} - \mathbf{w}_h^n) \cdot \nabla \mathbf{u}_h^n \cdot \mathbf{v}_h \nonumber \\
&\quad + \frac{\rho^f}{2} \int_{\Omega^{f,n}} \text{div}(\mathbf{u}_h^{n-1} - 2\mathbf{w}_h^n) \mathbf{u}_h^n \cdot \mathbf{v}_h \nonumber \\
&\quad - \frac{\rho^f}{2} \int_{\Sigma^n} (\mathbf{u}_h^{n-1} - \dot{\mathbf{d}}_h^{n-\frac{1}{2}}) \cdot n \mathbf{u}_h^n \cdot \mathbf{v}_h \nonumber \\
&\quad + 2\mu \int_{\Omega^{f,n}} \epsilon_h(\mathbf{u}_h^n) : \epsilon_h(\mathbf{v}_h) \nonumber \\
&\quad - \int_{\Omega^{f,n}} p_h^n \text{div} \mathbf{v}_h + \int_{\Omega^{f,n}} q_h \text{div} \mathbf{u}_h^n \nonumber \\
&\quad + s_h^f (\mathbf{u}_h^{n-1} - \mathbf{w}_h^n; (\mathbf{u}_h^n, p_h^n), (\mathbf{v}_h, q_h)). \tag{13}
\end{align}
$$

The second term in (13) involves the so-called Temam’s trick, which copes with the fact that the discrete fluid velocities are not divergence free. The third term is a weakly consistent stabilization term which allows to control the energy contributions of the convective term (first term of (13)) induced by the fact that

$$ \hat{\mathbf{u}}_h^{n-1}|_{\Sigma} \neq \mathbf{w}_h^n|_{\Sigma} = \dot{\mathbf{d}}_h^{n-\frac{1}{2}}|_{\Sigma}. $$

The last term in (13) stands for a given velocity and pressure stabilization operator (see, e.g., [26, 27]), whose purpose is to circumvent the issues related to the lack of robustness for high local Reynolds numbers and the lack of inf-sup stability for the pair $\mathbf{V}_h^f/M^f$. In the stability analysis of Section 3.1, we only use the positiveness of $s_h^f$, viz.,

$$ s_h^f(\mathbf{z}_h; (\mathbf{v}_h, q_h), (\mathbf{v}_h, q_h)) \geq 0 \qquad (14) $$
for all $(z_h, v_h, q_h) \in \mathbf{V}_h^f \times \mathbf{V}_h^f \times M_h^f$. Following [28, Chapter 5], we introduce a second weakly consistent stabilization term, namely,

$$s_{\Omega^{f,n}}(\mathbf{u}_h^n, \mathbf{v}_h) \stackrel{\text{def}}{=} -\frac{\rho^f}{2\tau} \left( \int_{\Omega^{f,n}} \mathbf{u}_h^n \cdot \mathbf{v}_h - \int_{\Omega^{f,n-1}} \mathbf{u}_h^n \cdot \mathbf{v}_h \right) + \frac{\rho^f}{2} \int_{\Omega^{f,n}} (\text{div} \, \mathbf{w}_h^n) \mathbf{u}_h^n \cdot \mathbf{v}_h, \quad (15)$$

whose purpose is to cope with the issues related to the so-called discrete geometric conservation law (see, e.g., [29]). Note that the weak consistency of (15) is a consequence of the Reynolds transport formula:

$$\frac{d}{dt} \int_{\Omega^f(t)} q = \int_{\Omega^f(t)} q \, \text{div} \, \mathbf{w}$$

for any smooth enough field $\hat{q} : \Omega^f \to \mathbb{R}^+$. Finally, we let $\mathcal{L}_h : \mathbf{V}_h^g \to \mathbf{V}_h^f$ denote a given discrete counterpart of the continuous lifting operator $\mathcal{L}$.

With all the above ingredients, the proposed fully discrete counterpart of (7)–(9) is detailed in Algorithm 1. A few remarks are now in order. The relation (18) in Algorithm 1 is essentially a rewriting of (11)$_2$, in which the intermediate variable $\lambda_h^n$ represents the discrete counterpart of the interfacial fluid stress on the reference configuration $\lambda^n$, given by (6). Indeed, by combining (17) and (18), the update of the fluid-stress in the fourth step of Algorithm 1 can be formulated, in a more standard form as a fluid variational residual (see, e.g., [3]):

$$\int_{\Sigma} \lambda_h^n \cdot \mu_h = \frac{\rho^f}{\tau} \left( \int_{\Omega^{f,n}} \mathbf{u}_h^n \cdot \mathcal{L}_h^f \mu_h - \int_{\Omega^{f,n-1}} \mathbf{u}_h^{n-1} \cdot \mathcal{L}_h^f \mu_h \right) + s_{\Omega^{f,n}}(\mathbf{u}_h^n, \mathcal{L}_h^f \mu_h) + a_{\Omega^{f,n}}(\mathbf{u}_h^{n-1}, \mathbf{w}_h^n; (\mathbf{u}_h^n, p_h^n), (\mathcal{L}_h^f \mu_h, 0)) \quad (19)$$

for all $\mu_h \in \mathbf{V}_h^g$ and where $\mathcal{L}_h^f : \mathbf{V}_h^g \to \mathbf{V}_h^f$ represents a linear fluid-sided lifting operator, such that, the nodal values of $\mathcal{L}_h^f \mu_h$ vanish out of $\Sigma$ and $(\mathcal{L}_h^f \mu_h)|_{\Sigma} = \mu_h$, for all $\mu_h \in \mathbf{V}_h^g$. We can hence see that, as indicated above, $\lambda_h^n$ is nothing but the discrete counterpart of the interfacial fluid stress on the reference configuration. It is worth noting that the interfacial stabilization term in (13) introduces a weakly consistent perturbation of the fluid stresses approximation provided by (19). Indeed, this relation formally writes

$$\int_{\Sigma} \lambda_h^n \cdot \mu_h \equiv \int_{\Sigma^n} \sigma(\mathbf{u}_h^n, p_h^n) \mathbf{n} \cdot \mu_h - \frac{\rho^f}{2} \int_{\Sigma^n} (\mathbf{u}_h^{n-1} - \dot{\mathbf{d}}_h^{n-\frac{1}{2}} \circ (\mathcal{A}_h^n)^{-1}) \cdot \mathbf{n} (\mathbf{u}_h^n \cdot \mu_h)$$

for all $\mu \in \mathbf{V}_h^g$, or point-wise

$$\lambda^n \equiv J^n \hat{\sigma}(\mathbf{u}^n, p^n) (\mathbf{F}^n)^{-\mathrm{T}} \hat{\mathbf{n}} - \frac{\rho^f J^n}{2} (\hat{\mathbf{u}}^{n-1} - \dot{\mathbf{d}}^{n-\frac{1}{2}}) \cdot ((\mathbf{F}^n)^{-\mathrm{T}} \hat{\mathbf{n}}) \hat{\mathbf{u}}^n \quad \text{on } \Sigma.$$
**Algorithm 1** Explicit Robin-Robin scheme (fully discrete)

Assume $u_h^0, d_h^{f,0}, d_h^0, \dot{d}_h^0$ and $\lambda_h^0$ to be given and, for $n \ge 1$, perform:

* Solve solid: find $(d_h^n, \dot{d}_h^n) \in V_h^s \times V_h^s$, with $\dot{d}_h^{n-\frac{1}{2}} = \partial_\tau d_h^n$ such that

$$ \rho^s \int_{\Omega^s} \partial_\tau \dot{d}_h^n \cdot \xi_h + \int_{\Omega^s} \bar{\Sigma}_h^{n-\frac{1}{2}} : \partial_d(E_h^{n-\frac{1}{2}}) \xi_h + \alpha \int_{\Sigma} (\dot{d}_h^{n-\frac{1}{2}} - \hat{u}_h^{n-1}) \cdot \xi_h \\ = - \int_{\Sigma} \lambda_h^{n-1} \cdot \xi_h \quad (16) $$

for all $\xi_h \in V_h^s$.

* Update fluid domain:

$$ d_h^{f,n} = \mathcal{L}_h(d_h^{n-1}|_{\Sigma}), \quad w_h^n = \partial_\tau d_h^{f,n}, \quad \mathcal{A}_h^n = I_{\Omega^f} + d_h^{f,n}, \quad \Omega^{f,n} = \mathcal{A}_h^n(\Omega^f). $$

* Solve fluid: find $(\hat{u}^n, \hat{p}^n) \in V_h^f \times Q_h$ such that

$$ \frac{\rho^f}{\tau} \left( \int_{\Omega^{f,n}} \mathbf{u}_h^n \cdot \mathbf{v}_h - \int_{\Omega^{f,n-1}} \mathbf{u}_h^{n-1} \cdot \mathbf{v}_h \right) + s_{\Omega^{f,n}} (\mathbf{u}_h^n, \mathbf{v}_h) \\ + a_{\Omega^{f,n}} (\mathbf{u}_h^{n-1}, \mathbf{w}_h^n; (\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{v}_h, q_h)) + \alpha \int_{\Sigma} (\hat{u}_h^n - \dot{\mathbf{d}}_h^{n-\frac{1}{2}}) \cdot \hat{\mathbf{v}}_h \\ = \int_{\Sigma} \lambda_h^{n-1} \cdot \hat{\mathbf{v}}_h \quad (17) $$

for all $(\hat{\mathbf{v}}_h, \hat{q}_h) \in V_h^f \times Q_h$.

* Fluid-stress update: set $\lambda_h^n \in V_h^g$ as

$$ \lambda_h^n = \lambda_h^{n-1} + \alpha(\dot{\mathbf{d}}_h^{n-\frac{1}{2}} - \hat{\mathbf{u}}_h^n) \quad \text{on } \Sigma. \quad (18) $$

It should also be noted that these last two expressions are merely formal and, therefore, can not be used in the energy stability analysis of Algorithm 1, presented in Section 3.1 below.

Algorithm 1 shares some similarities with the time-semidiscrete Robin-Robin scheme proposed in [21], but with several fundamental differences. First, the method is here formulated in a fully discrete setting, with a rigorous definition of the interfacial fluid-stresses. Second, a suitable interfacial stabilization is added to control the energy contributions of the convective term. This avoids the artificial modifications of the fluid velocity on the interface advocated in [21, Section 3]. Third, the stability of the method reported in the next paragraph also holds in 3D, irrespectively of the discrete geometric conservation law, via the stabilization term (15). Finally, the stability analysis (Section 3.1)
and the numerical examples (of Section 5) cover the case of the coupling with a general non-linear solid model.

## 3.1 Energy stability analysis

For $n \ge 1$, let the discrete energy and dissipation be respectively defined at time $t_n$ as

$$E_h^n \stackrel{\text{def}}{=} \frac{\rho^f}{2} \|\mathbf{u}_h^n\|_{0,\Omega^{f,n}}^2 + \frac{\rho^s}{2} \|\dot{\mathbf{d}}_h^n\|_{0,\Omega^s}^2 + \int_{\Omega^s} W(\mathbf{E}_h^n) + \frac{\tau\alpha}{2} \|\hat{\mathbf{u}}_h^n\|_{0,\Sigma}^2 + \frac{\tau}{2\alpha} \|\boldsymbol{\lambda}_h^n\|_{0,\Sigma}^2,$$

$$G_h^n \stackrel{\text{def}}{=} 2\mu \|\epsilon(\mathbf{u}_h^n)\|_{0,\Omega^{f,n}}^2 + \frac{\rho^f}{2\tau} \|\mathbf{u}_h^n - \mathbf{u}_h^{n-1}\|_{0,\Omega^{f,n-1}}^2 + \frac{\alpha}{2} \|\dot{\mathbf{d}}_h^{m-\frac{1}{2}} - \hat{\mathbf{u}}_h^{m-1}\|_{0,\Sigma}^2 \\ + s_h^f (\mathbf{u}_h^{n-1} - \mathbf{w}_h^n; (\mathbf{u}_h^n, p_h^n), (\mathbf{u}_h^n, p_h^n)).$$

Note that the positiveness of $s_h^f$ (i.e., relation (14)) is used in the definition of $G_h^n$.

**Theorem 1** Let $\{(\hat{\mathbf{u}}_h^n, \hat{p}_h^n, \mathbf{d}_h^n, \dot{\mathbf{d}}_h^n, \boldsymbol{\lambda}_h^n)\}_{n \ge 1}$ be given by Algorithm 1. The following energy identity holds for $n \ge 1$

$$E_h^n + \sum_{m=1}^{n} \tau G_h^m = E_h^0. \qquad (20)$$

As a result, Algorithm 1 is unconditionally stable in the energy norm.

*Proof* We first proceed by testing (16) with $\xi_h = \dot{\mathbf{d}}_h^{n-\frac{1}{2}}$, which yields

$$\frac{\rho^s}{2} \partial_{\tau} \|\dot{\mathbf{d}}_h^n\|_{0,\Omega^s}^2 + \int_{\Omega^s} \bar{\Sigma}_h^{n-\frac{1}{2}} : \partial_d(\mathbf{E}_h^{n-\frac{1}{2}}) \dot{\mathbf{d}}_h^{n-\frac{1}{2}} + \alpha \int_{\Sigma} (\dot{\mathbf{d}}_h^{n-\frac{1}{2}} - \hat{\mathbf{u}}_h^{n-1}) \cdot \dot{\mathbf{d}}_h^{n-\frac{1}{2}} \\ = - \int_{\Sigma} \boldsymbol{\lambda}_h^{n-1} \cdot \dot{\mathbf{d}}_h^{n-\frac{1}{2}}. \quad (21)$$

From (5), we have

$$\begin{aligned} \partial_d(\mathbf{E}_h^{n-\frac{1}{2}}) \dot{\mathbf{d}}_h^{n-\frac{1}{2}} &= \frac{1}{2} \left( (\mathbf{F}_h^{s,n-\frac{1}{2}})^{\mathrm{T}} \partial_{\tau} \mathbf{F}_h^{s,n} + (\partial_{\tau} \mathbf{F}_h^{s,n})^{\mathrm{T}} \mathbf{F}_h^{s,n-\frac{1}{2}} \right) \\ &= \frac{1}{2\tau} \left( (\mathbf{F}_h^{s,n})^{\mathrm{T}} \mathbf{F}_h^{s,n} - (\mathbf{F}_h^{s,n-1})^{\mathrm{T}} \mathbf{F}_h^{s,n-1} \right) \\ &= \partial_{\tau} \mathbf{E}_h^n. \end{aligned}$$

So that, by inserting the previous relation into (21) and using the definition of the stress correction (10), we have

$$\int_{\Omega^s} \bar{\Sigma}_h^{n-\frac{1}{2}} : \partial_{\tau} \mathbf{E}_h^n = \int_{\Omega^s} \partial_{\tau} W(\mathbf{E}_h^n). \qquad (22)$$

By inserting (22) into (21), we ultimately arrive at:

$$\frac{\rho^s}{2} \partial_{\tau} \|\dot{\mathbf{d}}_h^n\|_{0,\Omega^s}^2 + \partial_{\tau} \int_{\Omega^s} W(\mathbf{E}_h^n) + \alpha \int_{\Sigma} (\dot{\mathbf{d}}_h^{n-\frac{1}{2}} - \hat{\mathbf{u}}_h^{n-1}) \cdot \dot{\mathbf{d}}_h^{n-\frac{1}{2}} = - \int_{\Sigma} \boldsymbol{\lambda}_h^{n-1} \cdot \dot{\mathbf{d}}_h^{n-\frac{1}{2}}. \qquad (23)$$
12 Robin-Robin loose coupling for FSI

For the fluid, testing (17) with $(\hat{v}_h, \hat{q}_h) = (\hat{u}_h^n, \hat{p}_h^n)$ yields

$$ \frac{\rho^f}{\tau} \left( \int_{\Omega^{f,n}} \mathbf{u}_h^n \cdot \mathbf{u}_h^n - \int_{\Omega^{f,n-1}} \mathbf{u}_h^{n-1} \cdot \mathbf{u}_h^n \right) + a_{\Omega^{f,n}} (\mathbf{u}_h^{n-1}, \mathbf{w}_h^n; (\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{u}_h^n, \mathbf{p}_h^n)) \\ + s_{\Omega^{f,n}} (\mathbf{u}_h^n, \mathbf{u}_h^n) + \alpha \int_{\Sigma} (\hat{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{n-\frac{1}{2}}) \cdot \hat{\mathbf{u}}_h^n = \int_{\Sigma} \boldsymbol{\lambda}_h^{n-1} \cdot \hat{\mathbf{u}}_h^n. \quad (24) $$

The bulk terms are treated in a standard fashion, with the sole differences in the control of the ALE time derivative and of the convective terms, with the considered stabilizations. For the first term in (24), we have

$$ \frac{\rho^f}{\tau} \left( \int_{\Omega^{f,n}} \mathbf{u}_h^n \cdot \mathbf{u}_h^n - \int_{\Omega^{f,n-1}} \mathbf{u}_h^{n-1} \cdot \mathbf{u}_h^n \right) = \frac{\rho^f}{2} \partial_\tau \|\mathbf{u}_h^n\|_{0,\Omega^{f,n}}^2 \\ + \frac{\rho^f}{2\tau} \|\mathbf{u}_h^n - \mathbf{u}_h^{n-1}\|_{0,\Omega^{f,n-1}}^2 + \frac{\rho^f}{2\tau} \left( \|\mathbf{u}_h^n\|_{0,\Omega^{f,n}}^2 - \|\mathbf{u}_h^n\|_{0,\Omega^{f,n-1}}^2 \right). $$

On the other hand, using integration by parts and the fact that (by construction) $\mathbf{w}_h^n|_{\Sigma} = \dot{\mathbf{d}}_h^{n-\frac{1}{2}}|_{\Sigma}$, we get

$$ \begin{aligned} a_{\Omega^{f,n}} (\mathbf{u}_h^{n-1}, \mathbf{w}_h^n; (\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{u}_h^n, \mathbf{p}_h^n)) = & 2\mu \|\epsilon(\mathbf{u}_h^n)\|_{0,\Omega^{f,n}}^2 \\ & + \frac{\rho^f}{2} \int_{\Omega^{f,n}} (\mathbf{u}_h^{n-1} - \mathbf{w}_h^n) \cdot \nabla |\mathbf{u}_h^n|^2 + \frac{\rho^f}{2} \int_{\Omega^{f,n}} \operatorname{div}(\mathbf{u}_h^{n-1} - 2\mathbf{w}_h^n) |\mathbf{u}_h^n|^2 \\ & - \frac{\rho^f}{2} \int_{\Sigma^n} (\mathbf{u}_h^{n-1} - \dot{\mathbf{d}}_h^{n-\frac{1}{2}}) \cdot n |\mathbf{u}_h^n|^2 + s_h^f (\mathbf{u}_h^{n-1} - \mathbf{w}_h^n; (\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{u}_h^n, \mathbf{p}_h^n)) \\ = & 2\mu \|\epsilon(\mathbf{u}_h^n)\|_{0,\Omega^{f,n}}^2 - \frac{\rho^f}{2} \int_{\Omega^{f,n}} (\operatorname{div} \mathbf{w}_h^n) |\mathbf{u}_h^n|^2 + s_h^f (\mathbf{u}_h^{n-1} - \mathbf{w}_h^n; (\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{u}_h^n, \mathbf{p}_h^n)). \end{aligned} $$

By inserting the above two identities into (24), and by using (15) and (18), it follows that

$$ \begin{aligned} & \frac{\rho^f}{2} \partial_\tau \|\mathbf{u}_h^n\|_{0,\Omega^{f,n}}^2 + 2\mu \|\epsilon(\mathbf{u}_h^n)\|_{0,\Omega^{f,n}}^2 + s_h^f (\mathbf{u}_h^{n-1} - \mathbf{w}_h^n; (\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{u}_h^n, \mathbf{p}_h^n)) \\ & + \frac{\rho^f}{2\tau} \|\mathbf{u}_h^n - \mathbf{u}_h^{n-1}\|_{0,\Omega^{f,n-1}}^2 = \int_{\Sigma} \boldsymbol{\lambda}_h^{n-1} \cdot \hat{\mathbf{u}}_h^n - \alpha \int_{\Sigma} (\hat{\mathbf{u}}_h^n - \dot{\mathbf{d}}_h^{n-\frac{1}{2}}) \cdot \hat{\mathbf{u}}_h^n \\ & = \int_{\Sigma} \boldsymbol{\lambda}_h^n \cdot \hat{\mathbf{u}}_h^n. \end{aligned} \quad (25) $$

Finally, by adding (23) and (25), we get the following total energy balance

$$ \begin{aligned} & \frac{\rho^f}{2} \partial_\tau \|\mathbf{u}_h^n\|_{0,\Omega^{f,n}}^2 + \frac{\rho^s}{2} \partial_\tau \|\dot{\mathbf{d}}_h^n\|_{0,\Omega^s}^2 + \partial_\tau \int_{\Omega^s} W(E_h^n) + 2\mu \|\epsilon(\mathbf{u}_h^n)\|_{0,\Omega^{f,n}}^2 \\ & + \frac{\rho^f}{2\tau} \|\mathbf{u}_h^n - \mathbf{u}_h^{n-1}\|_{0,\Omega^{f,n-1}}^2 + s_h^f (\mathbf{u}_h^{n-1} - \mathbf{w}_h^n; (\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{u}_h^n, \mathbf{p}_h^n)) = R_h^n, \end{aligned} \quad (26) $$

where

$$ R_h^n \stackrel{\text{def}}{=} \alpha \int_{\Sigma} (\hat{\mathbf{u}}_h^{n-1} - \dot{\mathbf{d}}_h^{n-\frac{1}{2}}) \cdot \dot{\mathbf{d}}_h^{n-\frac{1}{2}} - \int_{\Sigma} \boldsymbol{\lambda}_h^{n-1} \cdot \dot{\mathbf{d}}_h^{n-\frac{1}{2}} + \int_{\Sigma} \boldsymbol{\lambda}_h^n \cdot \hat{\mathbf{u}}_h^n. $$

In order to control $R_h^n$, we can proceed as in the case of a fixed fluid domain (see [22]). By using twice the relation (18), we get

$$ R_h^n = \alpha \int_{\Sigma} (\hat{\mathbf{u}}_h^{n-1} - \hat{\mathbf{u}}_h^n) \cdot \dot{\mathbf{d}}_h^{n-\frac{1}{2}} + \frac{1}{\alpha} \int_{\Sigma} \boldsymbol{\lambda}_h^n \cdot (\boldsymbol{\lambda}_h^{n-1} - \boldsymbol{\lambda}_h^n). \quad (27) $$
On the other hand, by using the relation

$$ (a-b)c = \frac{1}{2} \left[ (b-c)^2 - (a-c)^2 + a^2 - b^2 \right], \qquad (28) $$

we have

$$ R_h^n = \frac{\alpha}{2} \left( \|\hat{u}_h^{n-1}\|_{0,\Sigma}^2 - \|\hat{u}_h^n\|_{0,\Sigma}^2 - \|\dot{\hat{d}}_h^{n-\frac{1}{2}} - \hat{u}_h^{n-1}\|_{0,\Sigma}^2 + \|\dot{\hat{d}}_h^{n-\frac{1}{2}} - \hat{u}_h^n\|_{0,\Sigma}^2 \right) \\ + \frac{1}{2\alpha} \left( \|\lambda_h^{n-1}\|_{0,\Sigma}^2 - \|\lambda_h^n\|_{0,\Sigma}^2 - \|\lambda_h^n - \lambda_h^{n-1}\|_{0,\Sigma}^2 \right), $$

so that, by using (18) once more, we get

$$ R_h^n = -\frac{\tau\alpha}{2} \partial_{\tau} \|\hat{u}_h^n\|_{0,\Sigma}^2 - \frac{\tau}{2\alpha} \partial_{\tau} \|\lambda_h^n\|_{0,\Sigma}^2 - \frac{\alpha}{2} \|\dot{\hat{d}}_h^{n-\frac{1}{2}} - \hat{u}_h^{n-1}\|_{0,\Sigma}^2. $$

Therefore, by inserting this identity into (26), we finally have

$$ \frac{\rho^f}{2} \partial_{\tau} \|\mathbf{u}_h^n\|_{0,\Omega^{f,n}}^2 + \frac{\rho^s}{2} \partial_{\tau} \|\dot{\mathbf{d}}_h^n\|_{0,\Omega^s}^2 + \partial_{\tau} \int_{\Omega^s} W(\mathcal{E}_h^n) + 2\mu \|\epsilon(\mathbf{u}_h^n)\|_{0,\Omega^{f,n}}^2 \\ + \frac{\rho^f}{2\tau} \|\mathbf{u}_h^n - \mathbf{u}_h^{n-1}\|_{0,\Omega^{f,n-1}}^2 + s_h^f(\mathbf{u}_h^{n-1} - \mathbf{w}_h^n; (\mathbf{u}_h^n, p_h^n), (\mathbf{u}_h^n, p_h^n)) \\ + \frac{\alpha\tau}{2} \partial_{\tau} \|\hat{u}_h^n\|_{0,\Sigma}^2 + \frac{\tau}{2\alpha} \partial_{\tau} \|\lambda_h^n\|_{0,\Sigma}^2 + \frac{\alpha}{2} \|\dot{\hat{d}}_h^{n-\frac{1}{2}} - \hat{u}_h^{n-1}\|_{0,\Sigma}^2 = 0. $$

The energy estimate (20) then follows by multiplying this expression by $\tau$ and then summing over $m = 1, \dots, n$, which completes the proof. $\square$

## 4 Error estimate in a linearized setting

In this section, we provide an error analysis of Algorithm 1 in a simplified setting. To this purpose, we assume that the solid undergoes infinitesimal displacements and that convective effects are negligible in the fluid. Under these assumptions, the non-linear coupled problem (1)-(3) can be replaced by the following linear Stokes-elasticity system in the reference configuration $\Omega^f \cup \Omega^s$:

$$ \left\{ \begin{array}{ll} \rho^f \partial_t \mathbf{u} - \mathrm{div} \, \boldsymbol{\sigma}(\mathbf{u}, p) = \mathbf{0} & \text{in } \Omega^f \times (0, T), \\ \mathrm{div} \, \mathbf{u} = 0 & \text{in } \Omega^f \times (0, T), \\ \mathbf{u} = \mathbf{0} & \text{on } \Gamma^f \times (0, T), \end{array} \right. \qquad (29) $$

$$ \left\{ \begin{array}{ll} \rho^s \partial_t \dot{\mathbf{d}} - \mathrm{div} \, \boldsymbol{\sigma}(\mathbf{d}) = \mathbf{0} & \text{in } \Omega^s \times (0, T), \\ \dot{\mathbf{d}} = \partial_t \mathbf{d} & \text{in } \Omega^s \times (0, T), \\ \mathbf{d} = \mathbf{0} & \text{on } \Gamma^s \times (0, T), \end{array} \right. \qquad (30) $$

$$ \left\{ \begin{array}{ll} \mathbf{u} = \dot{\mathbf{d}} & \text{on } \Sigma \times (0, T), \\ \boldsymbol{\sigma}(\mathbf{d})\mathbf{n} = -\boldsymbol{\sigma}(\mathbf{u}, p)\mathbf{n} & \text{on } \Sigma \times (0, T), \end{array} \right. \qquad (31) $$

complemented with the initial conditions: $\mathbf{d}(0, \cdot) = \mathbf{d}_0$, $\dot{\mathbf{d}}(0, \cdot) = \dot{\mathbf{d}}_0$ in $\Omega^s$, and $\mathbf{u}(0, \cdot) = \mathbf{u}_0$ in $\Omega^f$ and where $T > 0$ denotes the final time. The symbol $\boldsymbol{\sigma}(\mathbf{d})$
14 Robin-Robin loose coupling for FSI

stands for the linearized solid stress tensor, which is assumed to be given in terms of the solid displacement **d** by the following relation:

$$ \sigma(\mathbf{d}) \stackrel{\text{def}}{=} 2L_1\epsilon(\mathbf{d}) + L_2(\operatorname{div} \mathbf{d})\mathbf{I}, $$

where $L_1, L_2 > 0$ are the so called Lamé coefficients of the solid, given by the relations $L_1 = E/(1 + \nu)$ and $L_2 = \nu E/((1 + \nu)(1 - 2\nu))$ in terms of the Young's modulus $E$ and the Poisson's ratio $\nu$.

Assuming that $\lambda \stackrel{\text{def}}{=} \sigma(\mathbf{u}, p)\mathbf{n} \in \mathbf{L}^2(\Sigma)$, the solution of (29)–(31) satisfies the weak formulation

$$ \left\{ \begin{array}{l} \rho^s \int_{\Omega^s} \partial_t \dot{\mathbf{d}} \cdot \boldsymbol{\xi} + a^s(\mathbf{d}, \boldsymbol{\xi}) + \int_{\Sigma} \boldsymbol{\lambda} \cdot \boldsymbol{\xi} = 0, \\ \\ \rho^f \int_{\Omega^f} \partial_t \mathbf{u} \cdot \mathbf{v} + a^f((\mathbf{u}, p), (\mathbf{v}, q)) - \int_{\Sigma} \boldsymbol{\lambda} \cdot \mathbf{v} = 0, \\ \\ \int_{\Sigma} (\mathbf{u} - \dot{\mathbf{d}}) \cdot \boldsymbol{\mu} = 0 \end{array} \right. \qquad (32) $$

for all $(\boldsymbol{\xi}, \mathbf{v}, q) \in \mathbf{V}^s \times \mathbf{V}^f \times \mathbf{M}^f$. Here, the bilinear forms of the fluid, denoted as $a^f$, and of the solid, denoted as $a^s$, are respectively given by

$$ \begin{aligned} a^f((\mathbf{u}, p), (\mathbf{v}, q)) &\stackrel{\text{def}}{=} 2\mu \int_{\Omega^f} \boldsymbol{\epsilon}(\mathbf{u}) : \boldsymbol{\epsilon}(\mathbf{v}) - \int_{\Omega^f} p \operatorname{div} \mathbf{v} + \int_{\Omega^f} q \operatorname{div} \mathbf{u}, \\ a^s(\mathbf{d}, \boldsymbol{\xi}) &\stackrel{\text{def}}{=} 2L_1 \int_{\Omega^s} \boldsymbol{\epsilon}(\mathbf{d}) : \boldsymbol{\epsilon}(\boldsymbol{\xi}) + L_2 \int_{\Omega^s} (\operatorname{div} \mathbf{d})(\operatorname{div} \boldsymbol{\xi}). \end{aligned} \qquad (33) $$

We also introduce the elastic energy norm of the solid given by the relation

$$ \|\boldsymbol{\xi}\|_{\mathrm{s}} \stackrel{\mathrm{def}}{=} \sqrt{a^{\mathrm{s}}(\boldsymbol{\xi}, \boldsymbol{\xi})} \qquad (34) $$

for all $\boldsymbol{\xi} \in \mathbf{V}^s$.

For the numerical approximation of (32) we consider the Robin-Robin loosely coupled scheme introduced in [22] and detailed in Algorithm 2. As in Section 3, we employ a backward Euler method for the fluid and a mid-point scheme for the solid. In addition, for the spatial discretization, we consider a piecewise affine continuous finite elements, where we assume that the fluid and solid meshes are fitted. In other words, Algorithm 2 is essentially the linear counterpart of Algorithm 1. Here, we accommodate for the lack of inf-sup stability by using a pressure stabilized bilinear form for the fluid (see [30]),

$$ 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)) + \frac{\gamma_p h^2}{\mu} \int_{\Omega^f} \nabla p_h \cdot \nabla q_h, \quad (35) $$

with $\gamma_p > 0$ a user-defined parameter.
**Algorithm 2** Robin-based loosely coupled scheme (linear case).

Assume $\mathbf{u}_h^0$, $\mathbf{d}_h^0$, $\dot{\mathbf{d}}_h^0$ and $\boldsymbol{\lambda}_h^0$ to be given and, for $n \ge 1$, perform:

1. Solid subproblem: Find $\dot{\mathbf{d}}_h^n, \mathbf{d}_h^n \in \mathbf{V}_h^s$ such that $\dot{\mathbf{d}}_h^{n-\frac{1}{2}} = \partial_\tau \mathbf{d}_h^n$ and

$$ \rho^s \int_{\Omega^s} \partial_\tau \dot{\mathbf{d}}_h^n \cdot \boldsymbol{\xi}_h + a^s(\mathbf{d}_h^{n-\frac{1}{2}}, \boldsymbol{\xi}_h) + \alpha \int_{\Sigma} (\dot{\mathbf{d}}_h^{n-\frac{1}{2}} - \mathbf{u}_h^{n-1}) \cdot \boldsymbol{\xi}_h = - \int_{\Sigma} \boldsymbol{\lambda}_h^{n-1} \cdot \boldsymbol{\xi}_h \quad (36) $$

for all $\boldsymbol{\xi}_h \in \mathbf{V}_h^s$.

2. Fluid subproblem: Find $(\mathbf{u}_h^n, \mathbf{p}_h^n) \in \mathbf{V}_h^f \times \mathbf{M}_h^f$ such that

$$ \rho^f \int_{\Omega^f} \partial_\tau \mathbf{u}_h^n \cdot \mathbf{v}_h + a_h^f((\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{v}_h, q_h)) + \alpha \int_{\Sigma} (\mathbf{u}_h^n - \dot{\mathbf{d}}_h^{n-\frac{1}{2}}) \cdot \mathbf{v}_h = \int_{\Sigma} \boldsymbol{\lambda}_h^{n-1} \cdot \mathbf{v}_h \quad (37) $$

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

3. Fluid-stress update: set $\boldsymbol{\lambda}_h^n \in \mathbf{V}_h^g$ as

$$ \boldsymbol{\lambda}_h^n = \boldsymbol{\lambda}_h^{n-1} + \alpha(\dot{\mathbf{d}}_h^{n-\frac{1}{2}} - \mathbf{u}_h^n) \quad \text{on } \Sigma. \quad (38) $$

In [22], we proved that Algorithm 2 was optimal in space and sub-optimal in time, specifically $O(h + \tau^{\frac{1}{2}})$. However, following new numerical experiments and the results of [23], it has become clear that the result in [22] is not a strict bound. Following the techniques developed in [23], we are able to prove that the Robin-Robin scheme for the linear case is, in fact, nearly-optimal in time. This result is presented below (see also [24, Chapter 5]), along with the main features of the analysis.

## 4.1 Nearly-optimal error estimate

Let $N \in \mathbb{N}$ be such that $N\tau \le T$ and $\{(\mathbf{u}_h^n, \mathbf{p}_h^n, \boldsymbol{\lambda}_h^n, \mathbf{d}_h^n, \dot{\mathbf{d}}_h^n)\}_{1 \le n \le N} \subset \mathbf{V}_h^f \times \mathbf{M}_h^f \times \mathbf{V}_h^g \times \mathbf{V}_h^s \times \mathbf{V}_h^s$ be given by Algorithm 2. We introduce the error terms

$$ \left\{ \begin{aligned} e_{\mathbf{d}}^n &\stackrel{\text{def}}{=} \dot{\mathbf{d}}_h^n - \mathbf{d}^n = e_{\mathbf{d},h}^n + e_{\mathbf{d},\Pi}^n, & e_{\dot{\mathbf{d}}}^n &\stackrel{\text{def}}{=} \dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}^n = e_{\dot{\mathbf{d}},h}^n + e_{\dot{\mathbf{d}},\Pi}^n, \\ e_{\mathbf{u}}^n &\stackrel{\text{def}}{=} \mathbf{u}_h^n - \mathbf{u}^n = e_{\mathbf{u},h}^n + e_{\mathbf{u},\Pi}^n, & e_{\mathbf{p}}^n &\stackrel{\text{def}}{=} \mathbf{p}_h^n - \mathbf{p}^n = e_{\mathbf{p},h}^n + e_{\mathbf{p},\Pi}^n, \\ e_{\boldsymbol{\lambda}}^n &\stackrel{\text{def}}{=} \boldsymbol{\lambda}_h^n - \boldsymbol{\lambda}^n = e_{\boldsymbol{\lambda},h}^n + e_{\boldsymbol{\lambda},\Pi}^n, \end{aligned} \right. \quad (39) $$

$$ \left\{ \begin{aligned} e_{\mathbf{d},h}^n &\stackrel{\text{def}}{=} \dot{\mathbf{d}}_h^n - R_h^s \mathbf{d}^n, & e_{\dot{\mathbf{d}},h}^n &\stackrel{\text{def}}{=} \dot{\mathbf{d}}_h^n - R_h^s \dot{\mathbf{d}}^n, \\ e_{\mathbf{u},h}^n &\stackrel{\text{def}}{=} \mathbf{u}_h^n - R_h^f \mathbf{u}^n, & e_{\mathbf{p},h}^n &\stackrel{\text{def}}{=} \mathbf{p}_h^n - S_h \mathbf{p}^n, & e_{\boldsymbol{\lambda},h}^n &\stackrel{\text{def}}{=} \boldsymbol{\lambda}_h^n - \pi_h \boldsymbol{\lambda}^n, \end{aligned} \right. \quad (40) $$

$$ \left\{ \begin{aligned} e_{\mathbf{d},\Pi}^n &\stackrel{\text{def}}{=} R_h^s \mathbf{d}^n - \mathbf{d}^n, & e_{\dot{\mathbf{d}},\Pi}^n &\stackrel{\text{def}}{=} R_h^s \dot{\mathbf{d}}^n - \dot{\mathbf{d}}^n, \\ e_{\mathbf{u},\Pi}^n &\stackrel{\text{def}}{=} R_h^f \mathbf{u}^n - \mathbf{u}^n, & e_{\mathbf{p},\Pi}^n &\stackrel{\text{def}}{=} S_h \mathbf{p}^n - \mathbf{p}^n, & e_{\boldsymbol{\lambda},\Pi}^n &\stackrel{\text{def}}{=} \pi_h \boldsymbol{\lambda}^n - \boldsymbol{\lambda}^n, \end{aligned} \right. \quad (41) $$
16 Robin-Robin loose coupling for FSI

with $\lambda^n = \sigma(\mathbf{u}^n, \mathbf{p}^n)\mathbf{n}|_{\Sigma}$ and where $R_h^f, R_h^s$ respectively denote the Scott-Zhang interpolants (see, e.g., [31]) projecting onto the finite element spaces $\mathbf{V}_h^s$ and $\mathbf{V}_h^f$, and $S_h$ is the Scott-Zhang interpolant into $M_h^f$ modified by a global constant so the average on $\Omega^f$ is zero. For projecting into $\mathbf{V}_h^g$, we use the notation $\pi_h$ for the $L^2$ orthogonal projection from $\mathbf{V}^g$ to $\mathbf{V}_h^g$. We employ the Scott-Zhang interpolant for this analysis because it allows significant flexibility when selecting the degrees of freedom. Indeed, we may choose the degrees of freedom so that we will have $R_h^s\mathbf{v} = R_h^f\mathbf{w}$ on $\Sigma$ if $\mathbf{v} \in [H^1(\Omega^s)]^2$ and $\mathbf{w} \in [H^1(\Omega^f)]^2$ and $\mathbf{v} = \mathbf{w}$ on $\Sigma$. Furthermore, for simplicity, we assume that Algorithm 2 is initialized as follows

$$ 
\begin{aligned}
\dot{\mathbf{d}}_h^0 &= R_h^s \dot{\mathbf{d}}^0, \quad \dot{\mathbf{d}}_h^0 = R_h^s \dot{\mathbf{d}}^0, \quad \mathbf{u}_h^0 = R_h^f \mathbf{u}^0, \\
\boldsymbol{\lambda}_h^0 &= \pi_h \boldsymbol{\lambda}^0 = -\pi_h (\boldsymbol{\sigma}(\mathbf{d}_0) \mathbf{n}^s|_{\Sigma}),
\end{aligned} 
\quad (42) $$

which immediately yields

$$ e_{\mathbf{d},h}^0 = \mathbf{0}, \quad e_{\mathbf{d},h}^0 = \mathbf{0}, \quad e_{\mathbf{u},h}^0 = \mathbf{0}, \quad e_{\boldsymbol{\lambda},h}^0 = \mathbf{0}. \quad (43) $$

We may then define the following quantities to construct our error estimate:

$$ 
\begin{aligned}
S_h^n &\stackrel{\text{def}}{=} \frac{\rho^f}{2} \|e_{\mathbf{u},h}^n\|_{0,\Omega^f}^2 + \frac{\rho^s}{2} \|e_{\mathbf{d},h}^n\|_{0,\Omega^s}^2 + \frac{1}{2} \|e_{\mathbf{d},h}^n\|_{\mathrm{s}}^2, \\
\mathcal{E}_h^n &\stackrel{\text{def}}{=} \frac{\tau}{2} \alpha \|e_{\mathbf{u},h}^n\|_{0,\Sigma}^2 + \frac{\tau}{2\alpha} \|e_{\boldsymbol{\lambda},h}^n\|_{0,\Sigma}^2, \\
W_h^n &\stackrel{\text{def}}{=} \frac{\rho^f}{2} \|e_{\mathbf{u},h}^n - e_{\mathbf{u},h}^{n-1}\|_{0,\Omega^f}^2 + 2\mu\tau \|\boldsymbol{\varepsilon}(e_{\mathbf{u},h}^n)\|_{0,\Omega^f}^2 + \tau \frac{\gamma_p h^2}{\mu} \|\nabla e_{p,h}^n\|_{0,\Omega^s}^2, \\
Z_h^n &\stackrel{\text{def}}{=} \frac{\tau\alpha}{2} \|e_{\mathbf{d},h}^{n-\frac{1}{2}} - e_{\mathbf{u},h}^{n-1}\|_{0,\Sigma}^2
\end{aligned} 
\quad (44) $$

for $1 \le n \le N$. Finally, the quasi-optimal error result requires two key assumptions, which we state here. The motivation behind these assumptions will be explained further below. First, we assume that the normal $\mathbf{n}^s$ can be extended from $\Sigma$ into $\Omega^s$ in such a way that its gradient is bounded.

**Assumption 1** There exists $\tilde{\mathbf{n}}^s \in [\mathbf{W}^{1,\infty}(\Omega^s)]^2$ such that $\tilde{\mathbf{n}}^s|_{\Sigma} = \mathbf{n}^s$.

This assumption is necessary because the proof of the error estimate requires a way to evaluate the outward facing normal $\mathbf{n}^s$ (a quantity only defined on boundaries) at points in the interior. To do this, we must find an extension $\tilde{\mathbf{n}}^s$ of the normal vector from $\Sigma$ to $\Omega^s$ that is smooth enough to play well with the analysis. For many domains, this extension may be trivial. However, this assumption is made to acknowledge the possibility of more complicated domains. For example, if the domain is an annulus, then $\tilde{\mathbf{n}}^s$ will need to be a function that gradually changes direction from $\mathbf{n}^s$ on the outer boundary to $\mathbf{n}^s$ on the interior boundary.
The second assumption is the existence of a function $\tilde{\phi}: \Omega^s \to \mathbb{R}$ that serves as a pseudo-lifting operator from the interface $\Sigma$ into the interior of the solid domain $\Omega^s$.

**Assumption 2** Assume that $\tau < \frac{1}{2}$. There exists a function $\tilde{\phi}: \Omega^s \to \mathbb{R}$ satisfying:

(i) $0 \le \tilde{\phi} \le 1$;

(ii) $\tilde{\phi} \in \mathbf{V}^s$;

(iii) $|\{x \in \Sigma: \tilde{\phi}(x) \neq 1\}| \le C\tau$;

(iv) $\|\nabla \tilde{\phi}\|_{0,\Omega^s}^2 \le C(1 + \log \frac{1}{\tau})$;

where each $C > 0$ represents a general constant independent of $\tau$ and of the physical parameters.

The function $\tilde{\phi}$ should be viewed as a smooth cutoff function that is 1 on the interface $\Sigma$ and 0 on $\Gamma^s$. Here, smooth refers to the necessary control of the gradient (condition (iv)) that prevents a rapid transition from 1 to 0 near the interface. We refer to this as a pseudo-lifting operator because it allows to take functions on the interface and lift them into the bulk function space $\mathbf{V}^s$. In this regard, condition (iii) may be interpreted as allowing the lifted functions to go to 0 on $\Gamma^s$ while still matching the unlifted function on most of $\Sigma$.

**Remark 1** It should be noted that, while the above two assumptions may seem restrictive, an example has been constructed on a rectangular domain in [24]. This example may serve as a building block to construct examples on more complicated domains. Furthermore, it should be noted that the function $\tilde{\phi}$ is only a tool used for the error analysis below and does not need to be known in order to implement Algorithm 2.

With Assumptions 1 and 2 established, the quasi-optimal error results proceed as follows.

**Theorem 2** Under Assumptions 1 and 2, and assuming that (43) holds and that the solution of (29)-(31) is smooth enough so that $Y_h$ defined below is bounded, we have the following bound on the discrete error

$$ \max_{1 \le m \le n} \{S_h^m + \mathcal{E}_h^m\} + \sum_{m=1}^n (Z_h^m + W_h^m) \le \tau^2 \frac{T(L_1 + L_2)^3}{\alpha^2} \left(1 + \log \frac{1}{\tau}\right) \|\dot{\mathbf{d}}\|_{L^2(0,T; W^{1,\infty}(\Omega^s))} + Y_{\tau,h} \quad (45) $$

for $1 \le n \le N$ and where $Y_{\tau,h} \stackrel{\text{def}}{=} \tau^2 Y_1 + h^2 Y_2$ is given by

$$ Y_1 \stackrel{\text{def}}{=} \left[ \frac{T(L_1 + L_2)^3}{\alpha^2} \|\dot{\mathbf{d}}\|_{L^2(0,T; H^2(\Omega^s))}^2 + \frac{T(L_1 + L_2)^2}{\alpha} \|\nabla \dot{\mathbf{d}}\|_{L^2(0,T; H^1(\Sigma))}^2 \right] \quad (46) $$
18 Robin-Robin loose coupling for FSI

$$
\begin{align*}
&+ \tau^2 T (L_1 + L_2) \|\partial_t \boldsymbol{d}\|_{L^2(0,T; H^1(\Omega^s))}^2 + \left(\alpha + \frac{1}{\alpha} + \frac{1+\alpha^2}{\mu}\right) \|\partial_t \boldsymbol{u}\|_{L^2(0,T; L^2(\Sigma))}^2 \\
&+ \left(1 + \frac{1}{\alpha}\right) (L_1 + L_2)^2 \|\boldsymbol{\nabla} \dot{\boldsymbol{d}}\|_{L^2(0,T; L^2(\Sigma))}^2 + \tau^2 \rho^s T \|\partial_t^3 \boldsymbol{d}\|_{L^2(0,T; L^2(\Omega^s))}^2 \\
&+ \tau^2 \rho^f T \|\partial_t^3 \boldsymbol{u}\|_{L^2(0,T; L^2(\Omega^f))}^2 + (L_1 + L_2)^2 \|\boldsymbol{\nabla} \dot{\boldsymbol{d}}\|_{L^2(0,T; L^\infty(\Sigma))}^2 \\
&+ \rho^s T \|\partial_t^2 \dot{\boldsymbol{d}}\|_{L^2(0,T; L^2(\Omega^s))}^2 + \alpha \|\partial_t \boldsymbol{u}\|_{L^2(0,T; L^2(\Sigma))}^2 + \frac{\rho^s T}{\alpha} \|\partial_t^2 \boldsymbol{d}\|_{L^2(0,T; H^1(\Omega^s))}^2 \\
&+ \frac{\rho^s}{T} \|\partial_t \dot{\boldsymbol{d}}\|_{L^2(0,T; L^2(\Omega^s))}^2 + \frac{\rho^s (L_1 + L_2)^2}{T \alpha^2} \|\dot{\boldsymbol{d}}\|_{L^2(0,T; H^1(\Omega^s))}^2 \\
&\rho^s \|\partial_t \dot{\boldsymbol{d}}\|_{L^\infty(0,T; L^2(\Omega^s))}^2 + \frac{\rho^s (L_1 + L_2)^2}{\alpha^2} \|\dot{\boldsymbol{d}}\|_{L^\infty(0,T; H^1(\Omega^s))}^2
\end{align*}
$$

$$
Y_2 \stackrel{\text{def}}{=} \left[ (L_1 + L_2) \|\boldsymbol{d}\|_{L^\infty(0,T; H^2(\Omega^s))}^2 + h^2 \frac{T(L_1 + L_2)^2}{\alpha} \|\boldsymbol{\nabla} \dot{\boldsymbol{d}}\|_{L^2(0,T; H^1(\Sigma))}^2 \right. \\
+ \mu T \|\boldsymbol{u}\|_{L^\infty(0,T; H^2(\Omega^f))}^2 + \frac{T}{\mu} \|p\|_{L^\infty(0,T; H^1(\Omega^f))}^2 \\
+ T(L_1 + L_2) \|\dot{\boldsymbol{d}}\|_{L^2(0,T; H^2(\Omega^s))}^2 + (L_1 + L_2) \|\boldsymbol{d}\|_{L^\infty(0,T; H^2(\Omega^s))}^2 \\
+ \mu T \left( \|\boldsymbol{u}\|_{L^\infty(0,T; H^2(\Omega^f))}^2 + \|\boldsymbol{u}\|_{L^\infty(0,T; H^3(\Omega^f))}^2 \right) \\
+ \left( \alpha + \frac{1}{\alpha} + \frac{1+\alpha^2}{\mu} \right) T \|\boldsymbol{u}\|_{L^\infty(0,T; H^3(\Omega^f))}^2 + \rho^s T \|\partial_t \dot{\boldsymbol{d}}\|_{L^2(0,T; H^1(\Omega^s))}^2 \\
+ \rho^f T \|\partial_t \boldsymbol{u}\|_{L^2(0,T; H^1(\Omega^f))}^2 + \frac{T}{\mu} \|p\|_{L^\infty(0,T; H^1(\Omega^f))}^2 + \alpha T \|\boldsymbol{u}\|_{L^\infty(0,T; H^3(\Omega^f))}^2 \bigg].
$$

We conclude this section with a few fundamental observations. From (45),
it follows that the time splitting error in Algorithm 2 scales as

$$
\tau \left( \alpha^{-1} \sqrt{(L_1 + L_2)^3 (1 + \log \tau^{-1})} + \rho^s (L_1 + L_2)^2 + \sqrt{1 + \alpha^{-1}} (L_1 + L_2) \right. \\
\left. \qquad + \sqrt{\alpha + \alpha^{-1}} + \frac{\alpha}{\sqrt{\mu}} \right), \quad (46)
$$

which yields the nearly-optimal first-order accuracy $O(\tau\sqrt{1+\log\tau^{-1}})$ in time.
Owing to the dependence of (46) on $\alpha$ and on $\alpha^{-1}$, we can anticipate that
the accuracy Algorithm 2 is spoiled for (relatively) large or small values of
$\alpha$. Finally, it is also worth noting that the leading term in the splitting er-
ror (46) degrades for large values of the elastic parameters of the solid. This
can be mitigated by balancing the contributions of the error bound provided
by Theorem 2, which suggests the choice $\alpha = \gamma\sqrt{\rho^s E}$, with $\gamma > 0$ a user-
defined dimensionless parameter. Numerical evidence of these observations will
be provided in the numerical examples of Section 5.

4.2 Proof to Theorem 2

The argument consists of three parts. First, we begin with a preliminary bound.
Lemma 1 For $1 \le n \le N$, it holds

$$ \max_{1 \le m \le n} \{S_h^m + \varepsilon_h^m\} + \sum_{m=1}^n (W_h^m + Z_h^m) \lesssim D_h^n + \sum_{m=1}^n \frac{\tau}{\alpha} \int_{\Sigma} e_{\lambda,h}^{m-1} \cdot g_{3,h}^m, \quad (47) $$

where

$$ g_{3,h}^n \stackrel{\text{def}}{=} \frac{\alpha}{2} g_{1,h}^n - \pi_h g_{2,h}^n, \quad g_{1,h}^n \stackrel{\text{def}}{=} R_h^f (u^n - u^{n-1})|_{\Sigma}, \quad g_{2,h}^n \stackrel{\text{def}}{=} \lambda^n - \lambda^{n-1} \quad (48) $$

and

$$ D_h^n \stackrel{\text{def}}{=} \|e_{d,\Pi}^n\|_S^2 + \tau \sum_{m=1}^n \left[ \begin{aligned} & \mu \|\epsilon(e_{u,\Pi}^m)\|_{0,\Omega^f}^2 + \frac{1}{\mu} \|e_{p,\Pi}^m\|_{0,\Omega^f}^2 \\ & + T \|R_h^s (\dot{d}^{n-\frac{1}{2}} - \partial_\tau d^n)\|_S^2 + T \|\partial_\tau e_{d,\Pi}^m\|_S^2 + \frac{1}{T} \|e_{d,\Pi}^{m-\frac{1}{2}}\|_S^2 \\ & + \frac{\mu}{h^2} \left( \|e_{u,\Pi}^m\|_{0,\Omega^f}^2 + \|e_{u,\Pi}^m\|_{0,\Sigma}^2 \right) \\ & + \left(\alpha + \frac{\alpha^2}{\mu}\right) \|g_{1,h}^m\|_{0,\Sigma}^2 + \left(\frac{1}{\alpha} + \frac{1}{\mu}\right) \|g_{2,h}^m\|_{0,\Sigma}^2 \\ & + \rho^s T \|\partial_t \dot{d}^{m-\frac{1}{2}} - \partial_\tau R_h^s \dot{d}^m\|_{0,\Omega^s}^2 + \rho^f T \|\partial_t u^m - \partial_\tau R_h^f u^m\|_{0,\Omega^f}^2 \\ & + \frac{h^2}{\mu} \|\nabla S_h p^m\|_{0,\Omega^f}^2 + \left(\frac{1}{\alpha} + \frac{1}{\mu}\right) \|g_{3,h}^m\|_{0,\Sigma}^2 \end{aligned} \right]. $$

*Proof* The first step of the proof consists in showing that the discrete errors (40) satisfy Algorithm 2 with appropriate perturbations of the right-hand sides. To this purpose, we first note that since $\dot{d}_h^{n-\frac{1}{2}} = \partial_\tau d_h^n$, from (40) we get

$$ e_{\bullet, d,h}^{n-\frac{1}{2}} = \partial_\tau e_{d,h}^n - \underbrace{R_h^s (\dot{d}^{n-\frac{1}{2}} - \partial_\tau d^n)}_{g_{4,h}^n}. \quad (49) $$

On the other hand, from (32) it follows that

$$ \rho^s \int_{\Omega^s} \partial_\tau \dot{d}^n \cdot \xi_h + a^s (d^{n-\frac{1}{2}}, \xi_h) + \alpha \int_{\Sigma} (\dot{d}^{n-\frac{1}{2}} - u^{n-1}) \cdot \xi_h = - \int_{\Sigma} \lambda^{n-1} \cdot \xi_h \\ + \rho^s \int_{\Omega^s} (\partial_\tau \dot{d}^n - \partial_t \dot{d}^{n-\frac{1}{2}}) \cdot \xi_h + \alpha \int_{\Sigma} (\dot{d}^{n-\frac{1}{2}} - u^{n-1}) \cdot \xi_h - \int_{\Sigma} (\lambda^{n-\frac{1}{2}} - \lambda^{n-1}) \cdot \xi_h $$

for all $\xi_h \in V_h^s$. Therefore, by subtracting (36) from this relation and using the error expressions (39), we have

$$ \rho^s \int_{\Omega^s} \partial_\tau e_{\bullet, d,h}^n \cdot \xi_h + a^s (e_{d,h}^{n-\frac{1}{2}}, \xi_h) + \alpha \int_{\Sigma} (e_{\bullet, d,h}^{n-\frac{1}{2}} - e_{u,h}^{n-1}) \cdot \xi_h + \int_{\Sigma} e_{\lambda,h}^{n-1} \cdot \xi_h \\ = \rho^s \int_{\Omega^s} (\partial_t \dot{d}^{n-\frac{1}{2}} - \partial_\tau R_h^s \dot{d}^n) \cdot \xi_h - a^s (e_{d,\Pi}^{n-\frac{1}{2}}, \xi_h) - \alpha \int_{\Sigma} (R_h^s \dot{d}^{n-\frac{1}{2}} - R_h^f u^{n-1}) \cdot \xi_h \\ + \underbrace{\frac{1}{2} \int_{\Sigma} (\lambda^n - \lambda^{n-1}) \cdot \xi_h}_{T_4^s(\xi_h)} \quad (50) $$
20 Robin-Robin loose coupling for FSI

for all $\xi_h \in V_h^s$. Similarly, from (32) and (35), for the fluid we have

$$ \rho^f \int_{\Omega^f} \partial_\tau \boldsymbol{u}^n \cdot \boldsymbol{v}_h + a_h^f ((\boldsymbol{u}^n, p^n), (\boldsymbol{v}_h, q_h)) - \int_\Sigma \boldsymbol{\lambda}^n \cdot \boldsymbol{v}_h \\ = \rho^f \int_{\Omega^f} (\partial_\tau \boldsymbol{u}^n - \partial_\tau \boldsymbol{v}^n) \cdot \boldsymbol{v}_h + \frac{\gamma_p h^2}{\mu} \int_{\Omega^f} \nabla p^n \cdot \nabla q_h, $$

for all $(\boldsymbol{v}_h, q_h) \in V_h^f \times M_h^f$, so that by subtracting (37) from this relation and using the error expressions (39), we get

$$ \rho^f \int_{\Omega^f} \partial_\tau \boldsymbol{e}_{u,h}^n \cdot \boldsymbol{v}_h + a_h^f ((\boldsymbol{e}_{u,h}^n, e_{p,h}^n), (\boldsymbol{v}_h, q_h)) - \int_\Sigma \boldsymbol{e}_{\lambda,h}^n \cdot \boldsymbol{v}_h \\ = \underbrace{\rho^f \int_{\Omega^f} (\partial_\tau \boldsymbol{u}^n - \partial_\tau \boldsymbol{R}_h^f \boldsymbol{u}^n) \cdot \boldsymbol{v}_h}_{T_1^f(\boldsymbol{v}_h, q_h)} - \underbrace{\frac{\gamma_p h^2}{\mu} \int_{\Omega^f} \nabla S_h p^n \cdot \nabla q_h}_{T_3^f(\boldsymbol{v}_h, q_h)} - \underbrace{-a^f((\boldsymbol{e}_{u,\Pi}^n, e_{p,\Pi}^n), (\boldsymbol{v}_h, q_h))}_{T_2^f(\boldsymbol{v}_h, q_h)} \quad (51) $$

for all $(\boldsymbol{v}_h, q_h) \in V_h^f \times M_h^f$. Finally, using the fact that $\boldsymbol{u}^n|_{\Sigma} = \dot{\boldsymbol{d}}^n|_{\Sigma}$, the discrete error counterpart of the Robin relation (38) becomes

$$ \boldsymbol{e}_{\lambda,h}^n - \boldsymbol{e}_{\lambda,h}^{n-1} = \alpha(\boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}} - \boldsymbol{e}_{u,h}^n) + \boldsymbol{g}_{3,h}^n \quad \text{on } \Sigma. \quad (52) $$

We now proceed as in the proof of Theorem 1. By testing (50) with $\xi_h = \tau \boldsymbol{e}_{\bar{d},h}^{n-\frac{1}{2}}$ and using (49), we get

$$ \frac{\rho^s}{2} \|\boldsymbol{e}_{\cdot, \bar{d}, h}^n\|_{0, \Omega^s}^2 + \frac{1}{2} \|\boldsymbol{e}_{\bar{d}, h}^n\|_{\bar{d}, h}^2 = \frac{\rho^s}{2} \|\boldsymbol{e}_{\cdot, \bar{d}, h}^{n-1}\|_{0, \Omega^s}^2 + \frac{1}{2} \|\boldsymbol{e}_{\bar{d}, h}^{n-1}\|_{\bar{d}, h}^2 - \tau a^s (\boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}}, \boldsymbol{g}_{4,h}^n) \\ + \tau \sum_{i=1}^4 T_i^s (\boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}}) - \alpha \tau \int_{\Sigma} (\boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}} - \boldsymbol{e}_{u,h}^{n-1}) \cdot \boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}} - \tau \int_{\Sigma} \boldsymbol{e}_{\lambda,h}^{n-1} \cdot \boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}}. \quad (53) $$

On the other hand, by testing (51) with $(\boldsymbol{v}_h, q_h) = \tau(\boldsymbol{e}_{u,h}^n, \boldsymbol{e}_{p,h}^n)$ we get

$$ \frac{1}{2} \|\boldsymbol{e}_{u,h}^n\|_{0,\Omega^f}^2 + \frac{1}{2} \|\boldsymbol{e}_{u,h}^n - \boldsymbol{e}_{u,h}^{n-1}\|_{0,\Omega^f}^2 + 2\mu\tau \|\boldsymbol{\epsilon}(\boldsymbol{e}_{u,h}^n)\|_{0,\Omega^f}^2 + \tau s_h (e_{p,h}^n, e_{p,h}^n) \\ = \frac{1}{2} \|\boldsymbol{e}_{u,h}^{n-1}\|_{0,\Omega^f}^2 + \tau \sum_{i=1}^3 T_i^f (\boldsymbol{e}_{u,h}^n, e_{p,h}^n) + \tau \int_{\Sigma} \boldsymbol{e}_{\lambda,h}^n \cdot \boldsymbol{e}_{u,h}^n. $$

Hence, by adding this relation to (53) and using (44), we get

$$ S_h^n + W_h^n = S_h^{n-1} - \tau a^s (\boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}} \cdot \boldsymbol{g}_{4,h}^n) + \tau \sum_{i=1}^4 T_i^s (\boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}}) + \tau \sum_{i=1}^3 T_i^f (\boldsymbol{e}_{u,h}^n, e_{p,h}^n) \\ + \tau \tilde{R}_h^n, \quad (54) $$

with

$$ \tilde{R}_h^n \stackrel{\text{def}}{=} \alpha \int_{\Sigma} (\boldsymbol{e}_{u,h}^{n-1} - \boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}}) \cdot \boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}} - \int_{\Sigma} \boldsymbol{e}_{\lambda,h}^{n-1} \cdot \boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}} + \int_{\Sigma} \boldsymbol{e}_{\lambda,h}^n \cdot \boldsymbol{e}_{u,h}^n. $$

We can now use (52) twice, to get

$$ \tilde{R}_h^n = \alpha \int_{\Sigma} (\boldsymbol{e}_{u,h}^{n-1} - \boldsymbol{e}_{u,h}^n) \cdot \boldsymbol{e}_{\cdot, \bar{d}, h}^{n-\frac{1}{2}} + \frac{1}{\alpha} \int_{\Sigma} \boldsymbol{e}_{\lambda,h}^n \cdot (\boldsymbol{e}_{\lambda,h}^{n-1} - \boldsymbol{e}_{\lambda,h}^n) $$
$$ + \frac{1}{\alpha} \int_{\Sigma} g_{3,h}^{n} \cdot e_{\lambda,h}^{n} + \int_{\Sigma} g_{3,h}^{n} \cdot e_{\cdot,d,h}^{n-\frac{1}{2}}. $$

Therefore, using the relation (28), we have

$$ \tilde{R}_h^n = \frac{\alpha}{2} \left( \|e_{u,h}^{n-1}\|_{0,\Sigma}^2 - \|e_{u,h}^n\|_{0,\Sigma}^2 - \|e_{\cdot,d,h}^{n-\frac{1}{2}} - e_{u,h}^{n-1}\|_{0,\Sigma}^2 + \|e_{\cdot,d,h}^{n-\frac{1}{2}} - e_{u,h}^n\|_{0,\Sigma}^2 \right) \\ + \frac{1}{2\alpha} \left( \|e_{\lambda,h}^{n-1}\|_{0,\Sigma}^2 - \|e_{\lambda,h}^n\|_{0,\Sigma}^2 - \|e_{\lambda,h}^n - e_{\lambda,h}^{n-1}\|_{0,\Sigma}^2 \right) + \frac{1}{\alpha} \int_{\Sigma} g_{3,h}^n \cdot e_{\lambda,h}^n + \int_{\Sigma} g_{3,h}^n \cdot e_{\cdot,d,h}^{n-\frac{1}{2}}. $$

We conclude by taking the square of (52), which yields

$$ \frac{\alpha}{2} \|e_{\cdot,d,h}^{n-\frac{1}{2}} - e_{u,h}^n\|_{0,\Sigma}^2 = \frac{1}{2\alpha} \|e_{\lambda,h}^n - e_{\lambda,h}^{n-1}\|_{0,\Sigma}^2 + \frac{1}{2\alpha} \|g_{3,h}^n\|_{0,\Sigma}^2 - \frac{1}{\alpha} \int_{\Sigma} (e_{\lambda,h}^n - e_{\lambda,h}^{n-1}) \cdot g_{3,h}^n, $$

so that, by inserting this in the previous expression of $\tilde{R}_h^n$ we finally get the relation

$$ \tilde{R}_h^n = \frac{\alpha}{2} \left( \|e_{u,h}^{n-1}\|_{0,\Sigma}^2 - \|e_{u,h}^n\|_{0,\Sigma}^2 - \|e_{\cdot,d,h}^{n-\frac{1}{2}} - e_{u,h}^{n-1}\|_{0,\Sigma}^2 \right) \\ + \frac{1}{2\alpha} \left( \|e_{\lambda,h}^{n-1}\|_{0,\Sigma}^2 - \|e_{\lambda,h}^n\|_{0,\Sigma}^2 \right) + \frac{1}{\alpha} \int_{\Sigma} g_{3,h}^n \cdot e_{\lambda,h}^{n-1} + \int_{\Sigma} g_{3,h}^n \cdot e_{\cdot,d,h}^{n-\frac{1}{2}} + \frac{1}{2\alpha} \|g_{3,h}^n\|_{0,\Sigma}^2. $$

Therefore by inserting this identity into (54) and using (44), we get

$$ S_h^n + W_h^n + \mathcal{E}_h^n + Z_h^n = S_h^{n-1} + \mathcal{E}_h^{n-1} - \tau a^s (e_{\cdot,d,h}^{n-\frac{1}{2}}, g_{4,h}^m) + \tau \sum_{i=1}^4 T_i^s (e_{\cdot,d,h}^{n-\frac{1}{2}}) \\ + \tau \sum_{i=1}^3 T_i^f (e_{u,h}^m, e_{p,h}^m) + \frac{\tau}{\alpha} \int_{\Sigma} g_{3,h}^n \cdot e_{\lambda,h}^{n-1} + \tau \int_{\Sigma} g_{3,h}^n \cdot e_{\cdot,d,h}^{n-\frac{1}{2}} + \frac{\tau}{2\alpha} \|g_{3,h}^n\|_{0,\Sigma}^2. $$

Finally, by summing over $m = 1, \dots, n$ and by noting that $S_h^0 = \mathcal{E}_h^0 = 0$ (owing to (43)), we have

$$ \begin{aligned} & \max_{1 \le m \le n} S_h^m + \max_{1 \le m \le n} \mathcal{E}_h^m + \sum_{m=1}^n (W_h^m + Z_h^m) \\ & \le - \underbrace{\sum_{m=1}^n \tau a^s (e_{\cdot,d,h}^{m-\frac{1}{2}}, g_{4,h}^m)}_{T_1} + \underbrace{\sum_{m=1}^n \tau T_1^s (e_{\cdot,d,h}^{m-\frac{1}{2}})}_{T_2} + \underbrace{\sum_{m=1}^n \tau T_2^s (e_{\cdot,d,h}^{m-\frac{1}{2}})}_{T_3} \\ & \quad + \underbrace{\sum_{m=1}^n \tau T_3^s (e_{\cdot,d,h}^{m-\frac{1}{2}})}_{T_4} + \underbrace{\sum_{m=1}^n \tau T_4^s (e_{\cdot,d,h}^{m-\frac{1}{2}})}_{T_5} + \underbrace{\sum_{m=1}^n \tau T_1^f (e_{u,h}^m, e_{p,h}^m)}_{T_6} \\ & \quad + \underbrace{\tau T_2^f (e_{u,h}^m, e_{p,h}^m)}_{T_7} + \underbrace{\sum_{m=1}^n \tau T_3^f (e_{u,h}^m, e_{p,h}^m)}_{T_8} + \underbrace{\sum_{m=1}^n \tau \int_{\Sigma} g_{3,h}^m \cdot e_{\cdot,d,h}^{m-\frac{1}{2}}}_{T_9} \\ & \quad + \sum_{m=1}^n \frac{\tau}{2\alpha} \|g_{3,h}^m\|_{0,\Sigma}^2 + \sum_{m=1}^n \frac{\tau}{\alpha} \int_{\Sigma} g_{3,h}^m \cdot e_{\lambda,h}^{m-1}. \end{aligned} \quad (55) $$

We proceed by estimating separately each term $T_i$, $i = 1, \dots, 9$ in (55). In what follows, $\epsilon_i$, with $i = 1, \dots$, denote generic positive constants independent of the physical and discretization parameters. For the first term, using (43), we have
22 Robin-Robin loose coupling for FSI

$$
\begin{aligned}
T_1 &\le \sum_{m=1}^{n} \tau \|\mathbf{e}_{d,h}^{m-\frac{1}{2}}\|_{\mathrm{s}} \|\mathbf{g}_{4,h}^{m}\|_{\mathrm{s}} \le \frac{1}{2} \sum_{m=1}^{n} \tau \left( \|\mathbf{e}_{d,h}^{m}\|_{\mathrm{s}} + \|\mathbf{e}_{d,h}^{m-1}\|_{\mathrm{s}} \right) \|\mathbf{g}_{4,h}^{m}\|_{\mathrm{s}} \\
&\le \frac{\epsilon_1 \tau}{2T} \sum_{m=1}^{n} \|\mathbf{e}_{d,h}^{m}\|_{\mathrm{s}}^2 + \frac{\tau 2T}{\epsilon_1} \sum_{m=1}^{n} \|\mathbf{g}_{4,h}^{m}\|_{\mathrm{s}}^2 \le \epsilon_1 \max_{1 \le m \le n} S_h^m + \frac{\tau 2T}{\epsilon_1} \sum_{m=1}^{n} \|\mathbf{g}_{4,h}^{m}\|_{\mathrm{s}}^2.
\end{aligned}
$$

For the second, we have

$$
\begin{aligned}
T_2 &\le \sum_{m=1}^{n} \tau \rho^s \|\partial_t \dot{\mathbf{d}}^{m-\frac{1}{2}} - \partial_\tau R_h^s \dot{\mathbf{d}}^m\|_{0,\Omega^s} \|\mathbf{e}_{\dot{\mathbf{d}},h}^{m-\frac{1}{2}}\|_{0,\Omega^s} \\
&\le \frac{\tau \rho^s \epsilon_1}{2T} \sum_{m=1}^{n} \|\mathbf{e}_{\dot{\mathbf{d}},h}^{m-\frac{1}{2}}\|_{0,\Omega^s}^2 + \frac{\tau \rho^s T}{2\epsilon_1} \sum_{m=1}^{n} \|\partial_t \dot{\mathbf{d}}^{m-\frac{1}{2}} - \partial_\tau R_h^s \dot{\mathbf{d}}^m\|_{0,\Omega^s}^2 \\
&\le \epsilon_2 \max_{1 \le m \le n} S_h^m + \frac{\tau \rho^s T}{2\epsilon_2} \sum_{m=1}^{n} \|\partial_t \dot{\mathbf{d}}^{m-\frac{1}{2}} - \partial_\tau R_h^s \dot{\mathbf{d}}^m\|_{0,\Omega^s}^2.
\end{aligned}
$$

For the third term, we use the relation (49) to get

$$
T_3 = - \underbrace{\sum_{m=1}^{n} a^s (\mathbf{e}_{d,\Pi}^{m-\frac{1}{2}}, \mathbf{e}_{d,h}^m - \mathbf{e}_{d,h}^{m-1})}_{T_{3,1}} + \underbrace{\sum_{m=1}^{n} \tau a^s (\mathbf{e}_{d,\Pi}^{m-\frac{1}{2}}, \mathbf{g}_{4,h}^m)}_{T_{3,2}}.
$$

Using summation by parts and (43), one gets

$$
\begin{aligned}
T_{3,1} &= \frac{1}{2} \sum_{m=1}^{n} a^s (\mathbf{e}_{d,\Pi}^m, \mathbf{e}_{d,h}^m - \mathbf{e}_{d,h}^{m-1}) + \frac{1}{2} \sum_{m=1}^{n} a^s (\mathbf{e}_{d,\Pi}^{m-1}, \mathbf{e}_{d,h}^m - \mathbf{e}_{d,h}^{m-1}) \\
&= a^s (\mathbf{e}_{d,\Pi}^n, \mathbf{e}_{d,h}^n) - a^s (\mathbf{e}_{d,\Pi}^0, \mathbf{e}_{d,h}^0) - \sum_{m=1}^{n} \tau a^s (\partial_\tau \mathbf{e}_{d,\Pi}^m, \mathbf{e}_{d,h}^{m-\frac{1}{2}}) \\
&\le \|\mathbf{e}_{d,\Pi}^n\|_{\mathrm{s}} \|\mathbf{e}_{d,h}^n\|_{\mathrm{s}} + \frac{1}{2} \sum_{m=1}^{n} \tau \|\partial_\tau \mathbf{e}_{d,\Pi}^m\|_{\mathrm{s}} \left( \|\mathbf{e}_{d,h}^m\|_{\mathrm{s}} + \|\mathbf{e}_{d,h}^{m-1}\|_{\mathrm{s}} \right) \\
&\le \epsilon_3 \max_{1 \le m \le n} S_h^m + \frac{1}{2\epsilon_3} \|\mathbf{e}_{d,\Pi}^n\|_{\mathrm{s}}^2 + \frac{\epsilon_4 \tau}{2T} \sum_{m=1}^{n} \|\mathbf{e}_{d,h}^m\|_{\mathrm{s}}^2 + \frac{\tau 2T}{\epsilon_4} \sum_{m=1}^{n} \|\partial_\tau \mathbf{e}_{d,\Pi}^m\|_{\mathrm{s}}^2 \\
&\le (\epsilon_3 + \epsilon_4) \max_{1 \le m \le n} S_h^m + \frac{1}{2\epsilon_3} \|\mathbf{e}_{d,\Pi}^n\|_{\mathrm{s}}^2 + \frac{\tau 2T}{\epsilon_4} \sum_{m=1}^{n} \|\partial_\tau \mathbf{e}_{d,\Pi}^m\|_{\mathrm{s}}^2.
\end{aligned}
$$

On the other hand, the term $T_{3,2}$ can simply be handled as follows

$$
T_{3,2} \le \sum_{m=1}^{n} \tau \| \mathbf{e}_{d,\Pi}^{m-\frac{1}{2}} \|_{\mathrm{s}} \| \mathbf{h}_{3,h}^{m} \|_{\mathrm{s}} \le \frac{\tau}{2T} \sum_{m=1}^{n} \| \mathbf{e}_{d,\Pi}^{m-\frac{1}{2}} \|_{\mathrm{s}}^2 + \frac{\tau T}{2} \sum_{m=1}^{n} \| \mathbf{h}_{3,h}^{m} \|_{\mathrm{s}}^2.
$$

For the fourth term in (55), using (43) and the trace and Korn's inequalities, we have

$$
\begin{aligned}
T_4 ={}& \sum_{m=1}^{n} \tau \alpha \int_{\Sigma} (R_h^f \mathbf{u}^{m-1} - R_h^s \dot{\mathbf{d}}^{m-\frac{1}{2}}) \cdot (\mathbf{e}_{\dot{\mathbf{d}},h}^{m-\frac{1}{2}} - \mathbf{e}_{\mathbf{u},h}^{m-1}) \\
& + \sum_{m=1}^{n} \tau \alpha \int_{\Sigma} (R_h^f \mathbf{u}^{m-1} - R_h^s \dot{\mathbf{d}}^{m-\frac{1}{2}}) \cdot \mathbf{e}_{\mathbf{u},h}^{m-1}
\end{aligned}
$$
$$
\begin{aligned}
&\le \frac{\epsilon_5 \tau \alpha}{2} \sum_{m=1}^{n} \| e_{\dot{\mathbf{d}},h}^{m-\frac{1}{2}} - e_{\mathbf{u},h}^{m-1} \|_{0,\Sigma}^2 + 2\epsilon_6 \tau \mu \sum_{m=1}^{n} \| \epsilon(e_{\mathbf{u},h}^m) \|_{0,\Omega^f}^2 \\
&\qquad + \frac{\tau}{2} \left( \frac{\alpha}{\epsilon_5} + \frac{\alpha^2 C_{\text{TK}}^2}{4\epsilon_6 \mu} \right) \sum_{m=1}^{n} \| R_h^f \mathbf{u}^{m-1} - R_h^s \dot{\mathbf{d}}^{m-\frac{1}{2}} \|_{0,\Sigma}^2 \\
&\le \epsilon_5 \sum_{m=1}^{n} Z_h^n + \epsilon_6 \sum_{m=1}^{n} W_h^n + \frac{\tau}{2} \left( \frac{\alpha}{\epsilon_5} + \frac{\alpha^2 C_{\text{TK}}^2}{4\epsilon_6 \mu} \right) \sum_{m=1}^{n} \| R_h^f \mathbf{u}^{m-1} - R_h^s \dot{\mathbf{d}}^{m-\frac{1}{2}} \|_{0,\Sigma}^2.
\end{aligned}
$$

Similarly, for $T_5$, we have

$$
\begin{aligned}
T_5 &= \frac{\tau}{2} \sum_{m=1}^{n} \int_{\Sigma} (\lambda^m - \lambda^{m-1}) \cdot (e_{\dot{\mathbf{d}},h}^{m-\frac{1}{2}} - e_{\mathbf{u},h}^{m-1}) + \frac{\tau}{2} \sum_{m=1}^{n} \int_{\Sigma} (\lambda^m - \lambda^{m-1}) \cdot e_{\mathbf{u},h}^{m-1} \\
&\le \frac{\epsilon_7 \tau \alpha}{2} \sum_{m=1}^{n} \| e_{\dot{\mathbf{d}},h}^{m-\frac{1}{2}} - e_{\mathbf{u},h}^{m-1} \|_{0,\Sigma}^2 + 2\epsilon_8 \tau \mu \sum_{m=1}^{n} \| \epsilon(e_{\mathbf{u},h}^m) \|_{0,\Omega^f}^2 \\
&\qquad + \frac{\tau}{8} \left( \frac{1}{\alpha \epsilon_7} + \frac{C_{\text{TK}}^2}{4 \epsilon_8 \mu} \right) \sum_{m=1}^{n} \| \lambda^m - \lambda^{m-1} \|_{0,\Sigma}^2 \\
&\le \epsilon_7 \sum_{m=1}^{n} Z_h^n + \epsilon_8 \sum_{m=1}^{n} W_h^n + \frac{\tau}{8} \left( \frac{1}{\alpha \epsilon_7} + \frac{C_{\text{TK}}^2}{4 \epsilon_8 \mu} \right) \sum_{m=1}^{n} \| \lambda^m - \lambda^{m-1} \|_{0,\Sigma}^2.
\end{aligned}
$$

For $T_6$, we have

$$
\begin{aligned}
T_6 &\le \frac{\tau \rho^f \epsilon_9}{2T} \sum_{m=1}^{n} \| e_{\mathbf{u},h}^m \|_{0,\Omega^f}^2 + \frac{\tau \rho^f T}{2\epsilon_9} \sum_{m=1}^{n} \| \partial_t \mathbf{u}^m - \partial_\tau R_h^f \mathbf{u}^m \|_{0,\Omega^f}^2 \\
&\le \epsilon_9 \max_{1 \le m \le n} S_h^m + \frac{\tau \rho^f T}{2\epsilon_9} \sum_{m=1}^{n} \| \partial_t \mathbf{u}^m - \partial_\tau R_h^f \mathbf{u}^m \|_{0,\Omega^f}^2.
\end{aligned}
$$

Similarly, for $T_7$, we have

$$
\begin{aligned}
T_7 &\le \epsilon_{10} \tau \frac{\gamma_p h^2}{\mu} \sum_{m=1}^{n} \| \nabla e_{p,h}^m \|_{0,\Omega^f}^2 + \frac{\tau}{4\epsilon_{10}} \frac{\gamma_p h^2}{\mu} \sum_{m=1}^{n} \| \nabla S_h p^m \|_{0,\Omega^f}^2 \\
&\le \epsilon_{10} \sum_{m=1}^{n} W_h^n + \frac{\tau}{4\epsilon_{10}} \frac{\gamma_p h^2}{\mu} \sum_{m=1}^{n} \| \nabla S_h p^m \|_{0,\Omega^f}^2.
\end{aligned}
$$

For $T_8$, using Korn's inequality, we have

$$
\begin{aligned}
T_8 ={}& -\tau 2\mu \sum_{m=1}^{n} \int_{\Omega^f} \epsilon(e_{\mathbf{u},\Pi}^m) : \epsilon(e_{\mathbf{u},h}^m) + \tau \sum_{m=1}^{n} \int_{\Omega^f} e_{p,\Pi}^m \operatorname{div} e_{\mathbf{u},h}^m \\
& -\tau \sum_{m=1}^{n} \underbrace{\int_{\Omega^f} \operatorname{div} e_{\mathbf{u},\Pi}^m e_{p,h}^m}_{T_{8,1}} \le 2(\epsilon_{11}+\epsilon_{12})\tau\mu \sum_{m=1}^{n} \|\epsilon(e_{\mathbf{u},h}^m)\|_{0,\Omega^f}^2 + \frac{\tau\mu}{2\epsilon_{11}} \sum_{m=1}^{n} \|\epsilon(e_{\mathbf{u},\Pi}^m)\|_{0,\Omega^f}^2 \\
& + \frac{\tau C_K^2}{8\epsilon_{12}\mu} \sum_{m=1}^{n} \|e_{p,\Pi}^m\|_{0,\Omega^f}^2 + T_{8,1} \le (\epsilon_{11}+\epsilon_{12}) \sum_{m=1}^{n} W_h^n + \frac{\tau\mu}{2\epsilon_{11}} \sum_{m=1}^{n} \|\epsilon(e_{\mathbf{u},\Pi}^m)\|_{0,\Omega^f}^2 \\
& + \frac{\tau C_K^2}{8\epsilon_{12}\mu} \sum_{m=1}^{n} \|e_{p,\Pi}^m\|_{0,\Omega^f}^2 + T_{8,1}.
\end{aligned}
$$
24 Robin-Robin loose coupling for FSI

For the term $T_{8,1}$ we proceed by integration by parts and using the trace and Poincare's inequalities on the zero-mean function $e_{p,h}^m$, which yields

$$
\begin{align*}
T_{8,1} ={}& \tau \sum_{m=1}^{n} \int_{\Omega^f} \boldsymbol{e}_{\boldsymbol{u},\Pi}^{m} \cdot \boldsymbol{\nabla} e_{p,h}^{m} - \tau \sum_{m=1}^{n} \int_{\Sigma} e_{p,h}^{m} \boldsymbol{e}_{\boldsymbol{u},\Pi}^{m} \cdot \boldsymbol{n} \\
\le{}& \tau \epsilon_{13} \frac{\gamma_p h^2}{\mu} \sum_{m=1}^{n} \|\boldsymbol{\nabla} e_{p,h}^{m}\|_{0,\Omega^f}^2 + \frac{\tau}{4\epsilon_{13}} \frac{\mu}{\gamma_p h^2} \sum_{m=1}^{n} \|\boldsymbol{e}_{\boldsymbol{u},\Pi}^{m}\|_{0,\Omega^f}^2 \\
& + \tau C_P \sum_{m=1}^{n} \|\boldsymbol{\nabla} e_{p,h}^{m}\|_{0,\Omega^f} \|\boldsymbol{e}_{\boldsymbol{u},\Pi}^{m}\|_{0,\Sigma} \\
\le{}& (\epsilon_{13} + \epsilon_{14}) \tau \frac{\gamma_p h^2}{\mu} \sum_{m=1}^{n} \|\boldsymbol{\nabla} e_{p,h}^{m}\|_{0,\Omega^f}^2 + \frac{\tau}{4\epsilon_{13}} \frac{\mu}{\gamma_p h^2} \sum_{m=1}^{n} \|\boldsymbol{e}_{\boldsymbol{u},\Pi}^{m}\|_{0,\Omega^f}^2 \\
& + \frac{\tau C_P^2}{4\epsilon_{14}} \frac{\mu}{\gamma_p h^2} \sum_{m=1}^{n} \|\boldsymbol{e}_{\boldsymbol{u},\Pi}^{m}\|_{0,\Sigma}^2 \\
\le{}& (\epsilon_{13} + \epsilon_{14}) \sum_{m=1}^{n} W_h^n + \frac{\tau}{4\epsilon_{13}} \frac{\mu}{\gamma_p h^2} \sum_{m=1}^{n} \|\boldsymbol{e}_{\boldsymbol{u},\Pi}^{m}\|_{0,\Omega^f}^2 \\
& + \frac{\tau C_P^2}{4\epsilon_{14}} \frac{\mu}{\gamma_p h^2} \sum_{m=1}^{n} \|\boldsymbol{e}_{\boldsymbol{u},\Pi}^{m}\|_{0,\Sigma}^2.
\end{align*}
$$

Finally, for the term $T_9$, we proceed as for term $T_4$, using (43) and the trace and Korn's inequalities, which yields

$$
\begin{align*}
T_9 ={}& \tau \sum_{m=1}^{n} \int_{\Sigma} \boldsymbol{g}_{3,h}^{m} \cdot (\boldsymbol{e}_{\dot{\boldsymbol{d}}_h}^{m-\frac{1}{2}} - \boldsymbol{e}_{\boldsymbol{u},h}^{m-1}) + \tau \sum_{m=1}^{n} \int_{\Sigma} \boldsymbol{g}_{3,h}^{m} \cdot \boldsymbol{e}_{\boldsymbol{u},h}^{m-1} \\
\le{}& \frac{\tau \alpha \epsilon_{15}}{2} \sum_{m=1}^{n} \|\boldsymbol{e}_{\dot{\boldsymbol{d}}_h}^{m-\frac{1}{2}} - \boldsymbol{e}_{\boldsymbol{u},h}^{m-1}\|_{0,\Sigma}^2 + 2\tau \mu \epsilon_{16} \sum_{m=1}^{n} \|\boldsymbol{\epsilon}(\boldsymbol{e}_{\boldsymbol{u},h}^m)\|_{0,\Omega^f}^2 \\
& + \frac{\tau}{2} \left( \frac{1}{\alpha \epsilon_{15}} + \frac{C_{\mathrm{TK}}^2}{2\mu \epsilon_{16}} \right) \sum_{m=1}^{n} \|\boldsymbol{g}_{3,h}^m\|_{0,\Sigma}^2 \\
\le{}& \epsilon_{15} \sum_{m=1}^{n} Z_h^m + \epsilon_{16} \sum_{m=1}^{n} W_h^m + \frac{\tau}{2} \left( \frac{1}{\alpha \epsilon_{15}} + \frac{C_{\mathrm{TK}}^2}{2\mu \epsilon_{16}} \right) \sum_{m=1}^{n} \|\boldsymbol{g}_{3,h}^m\|_{0,\Sigma}^2.
\end{align*}
$$

The estimate (47) then follows by inserting into (55) all the above bounds on $T_i$, $i = 1, \dots, 9$, and by taking $\epsilon_i$ small enough for $i = 1, \dots, 15$. This concludes the proof. □

The key takeaway of Lemma 1 is that all of the terms in $D_h^n$ are either the norms of splitting error terms such as $g_{1,h}^n$ or interpolation errors such as $S_h p^n - p^n$. In both cases, the terms may be optimally bounded in terms of $h$ and $\tau$ in the final step of the proof. However, the remaining unbounded term in (47), namely,

$$
\sum_{m=1}^{n} \frac{\tau}{\alpha} \int_{\Sigma} \boldsymbol{e}_{\lambda,h}^{m-1} \cdot \boldsymbol{g}_{3,h}^{m}.
$$
requires special care. Indeed, if we attempt to bound this term by balancing $e_{\lambda,h}^{m-1}$ with terms appearing on the left-hand side of (47), we ultimately sacrifice a factor of $\tau^{\frac{1}{2}}$, resulting in suboptimal convergence, as done in [22]. The second part of the proof focuses on deriving an optimal bound for this term. Essentially, the difficulty resides on the interface $\Sigma$, so that to avoid a loss in accuracy we must somehow avoid the interface as much as possible. To do this, the key idea consists in leveraging the error equations of the solid (50) to pull the term $e_{\lambda,h}^{n-1}$ from the interface $\Sigma$ into the solid domain $\Omega^s$.

However, this requires the construction of an extension of $g_{3,h}^n$ into $\Omega^s$ which is in $V_h^s$. Critically, this means that the extension needs to vanish on $\Gamma^s$. Owing to (31)$_1$ and (48), for the first contribution of $g_{3,h}^n$ we have

$$g_{1,h}^n = R_h^s (\dot{d}^n - \dot{d}^{n-1}) \quad \text{on } \Sigma,$$

so that we can simply take

$$\tilde{g}_{1,h}^n \stackrel{\text{def}}{=} R_h^s (\dot{d}^n - \dot{d}^{n-1}), \qquad (56)$$

as extension of $g_{1,h}^n$ onto $\Omega^s$. For the remaining term $\pi_h g_2^n$ in $g_{3,h}^n$ the situation is more delicate. Thus, we will utilize the extended normal $\tilde{n}^s$ and the function $\tilde{\phi}$ introduced in Assumptions 1 and 2 to continuously lift $\pi_h g_2^n$ into $V_h^s$. To this purpose, we define the following three terms:

$$\tilde{\ell}_s \stackrel{\text{def}}{=} \sigma(d)\tilde{n}^s, \quad \mathcal{L}_s \stackrel{\text{def}}{=} \tilde{\phi}\tilde{\ell}_s, \\ \tilde{g}_2^n \stackrel{\text{def}}{=} -(\mathcal{L}_s^n - \mathcal{L}_s^{n-1}) = -\tilde{\phi}\sigma(d^n - d^{n-1})\tilde{n}^s.$$

Note that, from (31)$_2$ we have $\sigma(d^n)n^s = -\lambda^n$ on $\Sigma$, so that $\tilde{g}_2^n = \tilde{\phi}g_2^n$ represents an extension of $g_2^n$ into $\Omega^s$. We can then define

$$\tilde{g}_{2,h}^n \stackrel{\text{def}}{=} R_h^s \tilde{g}_2^n = -R_h^s (\mathcal{L}_s^n - \mathcal{L}_s^{n-1}) = -(\mathcal{L}_{s,h}^n - \mathcal{L}_{s,h}^{n-1}). \qquad (57)$$

Additionally, using Assumptions 1 and 2, we can show that (see [24, Lemma 5.2.4]):

$$\|\tilde{g}_{2,h}^n - g_2^n\|_{0,\Sigma}^2 \lesssim \tau \|g_2^n\|_{L^\infty(\Sigma)}^2 + h^2 \|g_2^n\|_{1,\Sigma}^2. \qquad (58)$$

The result of this extension step yields the bound stated in the following lemma.

**Lemma 2** Under Assumptions 1 and 2, for any $\epsilon > 0$ we have

$$\frac{\tau}{\alpha} \sum_{m=1}^{n} \int_{\Sigma} e_{\lambda,h}^{m-1} \cdot g_{3,h}^{m} \lesssim \epsilon \max_{1 \le m \le n} \{S_h^m + \mathcal{E}_h^m\} + \epsilon \sum_{m=1}^{n} Z_h^m + \frac{1}{\epsilon} \Psi_h^n, \qquad (59)$$
26 Robin-Robin loose coupling for FSI

with

$$
\begin{align*}
\Psi_h^n &\stackrel{\text{def}}{=} \frac{T h^2}{\alpha} \sum_{m=1}^{n} \|\mathbf{g}_2^m\|_{1,\Sigma}^2 + \frac{\rho^s}{\alpha^2} \|\tilde{\mathbf{g}}_{5,h}^n\|_{0,\Omega^s}^2 \\
&\quad + \tau \sum_{m=1}^{n} \left( \|\mathbf{g}_2^m\|_{L^\infty(\Sigma)}^2 + \frac{\rho^s T \tau^2}{\alpha^2} \|\partial_\tau^2 \mathcal{L}_{s,h}^n\|_{0,\Omega^s}^2 + \rho^s T \tau^2 \|\partial_\tau^2 R_h^s \dot{\mathbf{d}}^m\|_{0,\Omega^s}^2 \right. \\
&\quad \left. + \frac{T}{\alpha^2} \|\tilde{\mathbf{g}}_{5,h}^m\|_s^2 + \frac{1}{\alpha} \|\tilde{\mathbf{g}}_{5,h}^m\|_{0,\Sigma}^2 + T \rho^s \|\partial_t \dot{\mathbf{d}}^{m-\frac{1}{2}} - \partial_\tau R_h^s \dot{\mathbf{d}}^m\|_{0,\Omega^s}^2 + \frac{\rho^s}{T \alpha^2} \|\tilde{\mathbf{g}}_{5,h}^m\|_{0,\Omega^s}^2 \right. \\
&\quad \left. + \frac{1}{T} \|\mathbf{e}_{\mathbf{d},\Pi}^{m-\frac{1}{2}}\|_s^2 + \alpha \|R_h^s \dot{\mathbf{d}}^{m-\frac{1}{2}} - R_h^f \mathbf{u}^{m-1}\|_{0,\Sigma}^2 + \frac{1}{\alpha} \|\mathbf{g}_2^m\|_{0,\Sigma}^2 \right),
\end{align*}
$$

and $\mathbf{g}_{5,h}^m \stackrel{\text{def}}{=} \frac{\alpha}{2} \tilde{\mathbf{g}}_{1,h}^m - \tilde{\mathbf{g}}_{2,h}^m$, where $\tilde{\mathbf{g}}_{1,h}^m$ and $\tilde{\mathbf{g}}_{2,h}^m$ are given by (56) and (57), respectively.

Proof We first note that, owing to the definitions (48) and (56), we have the decomposition

$$
\mathbf{g}_{3,h}^{m} = \tilde{\mathbf{g}}_{2,h}^{m} - \pi_{h} \mathbf{g}_{2}^{m} + \mathbf{g}_{5,h}^{m} \quad \text{on } \Sigma.
$$

Therefore,

$$
\begin{equation}
\begin{split}
\frac{\tau}{\alpha} \sum_{m=1}^{n} \int_{\Sigma} \mathbf{e}_{\lambda,h}^{m-1} \cdot \mathbf{g}_{3,h}^{m} &= \underbrace{\frac{\tau}{\alpha} \sum_{m=1}^{n} \int_{\Sigma} \mathbf{e}_{\lambda,h}^{m-1} \cdot (\tilde{\mathbf{g}}_{2,h}^{m} - \mathbf{g}_{2}^{m})}_{T_1} + \underbrace{\frac{\tau}{\alpha} \sum_{m=1}^{n} \int_{\Sigma} \mathbf{e}_{\lambda,h}^{m-1} \cdot \mathbf{g}_{5,h}^{m}}_{T_2} \\
&= \underbrace{\frac{\tau}{\alpha} \sum_{m=1}^{n} \int_{\Sigma} \mathbf{e}_{\lambda,h}^{m-1} \cdot \mathbf{g}_{5,h}^{m}}_{T_2}.
\end{split}
\tag{60}
\end{equation}
$$

Using (58) and (43), the first term in the right-hand side can be bounded as

$$
\begin{align*}
T_1 &\le \frac{T}{2\epsilon\alpha} \sum_{m=1}^{n} \|\tilde{\mathbf{g}}_{2,h}^m - \mathbf{g}_2^m\|_{0,\Sigma}^2 + \frac{\epsilon_1 \tau^2}{2T\alpha} \sum_{m=1}^{n} \|\mathbf{e}_{\lambda,h}^m\|_{0,\Sigma}^2 \\
&\le \frac{CT}{2\epsilon\alpha} \sum_{m=1}^{n} \left( \tau \|\mathbf{g}_2^m\|_{L^\infty(\Sigma)}^2 + h^2 \|\mathbf{g}_2^m\|_{1,\Sigma}^2 \right) + \epsilon \max_{1 \le m \le n} \mathcal{E}_h^m
\end{align*}
$$

for $\epsilon > 0$. In order to bound the second term in (60), we test the solid error equation
(50) with $\boldsymbol{\xi}_h = \frac{\tau}{\alpha} \tilde{\boldsymbol{g}}_{5,h}^n$ and sum the resulting expression over $m = 1 \dots n$, which yields

$$
\begin{align}
T_2 ={}& \underbrace{\frac{\rho^s}{\alpha} \sum_{m=1}^n \int_{\Omega^s} (\mathbf{e}_{\mathbf{d},h}^m - \mathbf{e}_{\mathbf{d},h}^{m-1}) \cdot \tilde{\mathbf{g}}_{5,h}^m}_{T_{2,1}} + \underbrace{\frac{\tau}{\alpha} \sum_{m=1}^n a^s (\mathbf{e}_{\mathbf{d},h}^{m-\frac{1}{2}}, \tilde{\mathbf{g}}_{5,h}^m)}_{T_{2,2}} \nonumber \\
& + \underbrace{\tau \sum_{m=1}^n \int_{\Sigma} (\mathbf{e}_{\mathbf{d},h}^{m-\frac{1}{2}} - \mathbf{e}_{\mathbf{u},h}^{m-1}) \cdot \tilde{\mathbf{g}}_{5,h}^m}_{T_{2,3}} - \underbrace{\frac{\tau}{\alpha} \sum_{m=1}^n \sum_{i=1}^4 T_i^s (\tilde{\mathbf{g}}_{5,h}^m)}_{T_{2,4}}. \tag{61}
\end{align}
$$

We now proceed by estimating each term separately. For the first term, using summation by parts and (43), we get

$$
T_{2,1} = \frac{\rho^s}{\alpha} \int_{\Omega^s} \mathbf{e}_{\mathbf{d},h}^m \cdot \tilde{\mathbf{g}}_{5,h}^n - \frac{\rho^s}{\alpha} \sum_{m=1}^n \int_{\Omega^s} (\tilde{\mathbf{g}}_{5,h}^m - \tilde{\mathbf{g}}_{5,h}^{m-1}) \cdot \mathbf{e}_{\mathbf{d},h}^{m-1} \le \frac{\epsilon \rho^s}{2} \|\mathbf{e}_{\mathbf{d},h}^n\|_{0,\Omega^s}^2
$$
$$ 
\begin{aligned}
&+ \frac{\rho^s}{2\epsilon\alpha^2} \|\tilde{g}_{5,h}^n\|_{0,\Omega^s}^2 + \frac{\epsilon\rho^s\tau}{2T} \sum_{m=1}^n \|e_{\cdot, \tilde{d}, h}^{m-1}\|_{0,\Omega^s}^2 + \frac{\rho^s T}{2\tau\epsilon\alpha^2} \sum_{m=1}^n \|\tilde{g}_{5,h}^m - \tilde{g}_{5,h}^{m-1}\|_{0,\Omega^s}^2 \\
&\le 2\epsilon \max_{1 \le m \le n} S_h^m + \frac{\rho^s}{2\epsilon\alpha^2} \|\tilde{g}_{5,h}^n\|_{0,\Omega^s}^2 + \frac{\rho^s T \tau^3}{\epsilon\alpha^2} \sum_{m=1}^n \|\partial_\tau^2 \mathcal{L}_{s,h}^n\|_{0,\Omega^s}^2 \\
&\quad + \frac{\rho^s T \tau^3}{\epsilon} \sum_{m=1}^n \|\partial_\tau^2 R_h^s \dot{\tilde{d}}^m\|_{0,\Omega^s}^2.
\end{aligned}
$$

For the second term in (61), we have

$$ T_{2,2} \le \frac{\tau \epsilon}{2T} \sum_{m=1}^{n} \|e_{\cdot, \tilde{d}, h}^{n}\|_{s}^{2} + \frac{\tau T}{2\epsilon\alpha^2} \sum_{m=1}^{n} \|\tilde{g}_{5,h}^{m}\|_{s}^{2} \le \epsilon \max_{1 \le m \le n} S_{h}^{m} + \frac{\tau T}{2\epsilon\alpha^2} \sum_{m=1}^{n} \|\tilde{g}_{5,h}^{m}\|_{s}^{2}. $$

Similarly, for the third term, we have

$$ 
\begin{aligned}
T_{2,3} \le \frac{\tau \alpha \epsilon}{2} \sum_{m=1}^{n} \|e_{\cdot, \tilde{d}, h}^{m-\frac{1}{2}} - e_{\cdot, u, h}^{m-1}\|_{0,\Sigma}^2 &+ \frac{\tau}{2\epsilon\alpha} \sum_{m=1}^{n} \|\tilde{g}_{5,h}^m\|_{0,\Sigma}^2 \\
&\le \epsilon \sum_{m=1}^{n} Z_h^m + \frac{\tau}{2\epsilon\alpha} \sum_{m=1}^{n} \|\tilde{g}_{5,h}^m\|_{0,\Sigma}^2.
\end{aligned}
$$

Finally, for the last term we have

$$ 
\begin{aligned}
T_{2,4} \le \frac{\tau \rho^s T}{2} \sum_{m=1}^{n} \|\partial_t \dot{\tilde{d}}^{m-\frac{1}{2}} - \partial_\tau R_h^s \dot{\tilde{d}}^m\|_{0,\Omega^s}^2 &+ \frac{\tau \rho^s}{2T\alpha^2} \sum_{m=1}^{n} \|\tilde{g}_{5,h}^m\|_{0,\Omega^s}^2 \\
&+ \frac{\tau}{2T} \sum_{m=1}^{n} \|e_{\cdot, \tilde{d}, \Pi}^{m-\frac{1}{2}}\|_{s}^{2} + \frac{\tau T}{2\alpha^2} \sum_{m=1}^{n} \|\tilde{g}_{5,h}^{m}\|_{s}^{2} + \frac{\tau \alpha}{2} \sum_{m=1}^{n} \|R_h^s \dot{\tilde{d}}^{m-\frac{1}{2}} - R_h^f u^{m-1}\|_{0,\Sigma}^2 \\
&+ \frac{\tau}{\alpha} \sum_{m=1}^{n} \|\tilde{g}_{5,h}^{m}\|_{0,\Sigma}^2 + \frac{\tau}{4\alpha} \sum_{m=1}^{n} \|g_2^m\|_{0,\Sigma}^2.
\end{aligned}
$$

The estimate (59) then follows by inserting all the above bounds into (60), which completes the proof. □

Finally, in the third step of the proof, it remains to bound the terms $D_h^n$ and $\Psi_h^n$ involved in Lemmas 1 and 2. This is the purpose of the next result, whose proof relies on the properties of the operators $R_h^f$, $R_h^s$, $S_h$ and $\pi_h$, Taylor expansions in time and Assumptions 1 and 2.

Lemma 3 Under Assumptions 1 and 2, there follows

$$ D_h^n + \Psi_h^n \lesssim \tau^2 \frac{T(L_1 + L_2)^3}{\alpha^2} \left(1 + \log \frac{1}{\tau}\right) \|\dot{\tilde{d}}\|_{L^2(0,T; W^{1,\infty}(\Omega^s))}^2 + Y_{\tau,h} \quad (62) $$

for $1 \le n \le N$.

*Proof* See Appendix B. □

The error estimate (45) then follows by combining the results of Lemmas 1–3.
# 5 Numerical examples

In this section, we illustrate the accuracy properties of Algorithms 1 in a series of numerical examples.

## 5.1 Convergence study with a manufactured solution

The objective of this first example is to provide numerical evidence of the nearly-optimal convergence result of Theorem 2 in an academic test case using a manufactured solution. To this purpose, we consider the simplified coupled problem (29)-(31) with $\Omega^f = [0, 1] \times [0, \frac{3}{4}]$, $\Omega^s = [0, 1] \times [\frac{3}{4}, 1]$. All the physical parameters $\rho^s$, $\rho^f$, $\mu$, $L_1$, $L_2$ are set to 1.

We then take the exact solution of (29)-(31) given by

$$ \mathbf{u}(x, y, t) = 10^{-3} e^t \begin{pmatrix} 2x(1-x)y(1-y) \\ x(1-x)(1-y) \end{pmatrix}, $$

$$ p(x, y, t) = 10^{-3} e^t L_2 (2(1 - 2x)y(1 - y) + x(1 - x)(1 - 2y)), $$

$$ \mathbf{d}(x, y, t) = 10^{-3} e^t \begin{pmatrix} 2x(1 - x)y(1 - y) \\ x(1 - x)(1 - y) \end{pmatrix}, $$

with the corresponding non-homogeneous right-hand side (in both $\Omega^f$ and in $\Omega^s$):

$$ 10^{-3} e^t \left[ \begin{pmatrix} 2x(1-x)y(1-y) + 12y(1-y) \\ x(1-x)y(1-y) + 6x(1-x) \end{pmatrix} + \begin{pmatrix} 4x(1-x) - 2(1-2x)(1-2y) \\ 2y(1-y) - 4(1-2x)(1-2y) \end{pmatrix} \right]. $$

Note that, the divergence free condition is not preserved, so that the following right-hand side has to be included in the continuity equation (29)$_2$

$$ 10^{-3} e^t (2(1 - 2x)y(1 - y) + x(1 - x)(1 - 2y)). $$

In order to illustrate the accuracy of Algorithm 2, we evaluate the $L^2$-velocity error $\|\dot{\mathbf{d}}_h^N - \dot{\mathbf{d}}(T)\|_{0,\Omega^s}$ at the final time $T = 0.5$ for different values of the discretization parameters $\tau = h \in \{0.5^i\}_{i=2}^7$. The numerical computations were performed using the FEniCS automated finite element solver [32].

In Figure 1, we report the error convergence history with respect to $\tau$ (in log-log scale) for different values of the Robin parameter $\alpha$. We can observe that Algorithm 2 delivers first-order accuracy for almost all values of $\alpha$. This is especially true for $\alpha = 10, 15, 20$, whereas for smaller or larger values of $\alpha$, smaller time-steps are required to observe the first-order accuracy provided by the theory (Theorem 2).

## 5.2 Convergence study in a simplified framework

In this second example we provide numerical evidence of the nearly-optimal convergence result of Theorem 2 in a more realistic configuration, namely, the
<figure><img src="image_3.png" /><figcaption>Figure 1 Convergence history of the solid velocity error.</figcaption></figure>

pressure wave propagation benchmark of [13, Section 6.1.1], We also investigate the effect on the accuracy of the solid parameters $L_1, L_2$ involved in the error estimate (45). In (29)-(31), we take $\Omega^f = [0, L] \times [0, R]$, $\Omega^s = [0, L] \times [R, R+\epsilon]$, with $L = 6$, $R = 0.5$ and $\epsilon = 0.1$. All the units are in the CGS system. At the fluid inlet boundary $x = 0$ a sinusoidal pressure of maximal amplitude $2 \cdot 10^4$ is enforced during $5 \cdot 10^{-3}$ second. Zero traction is prescribed at the outlet $x = 6$ and a symmetry condition is enforced on the lower wall. Transverse membrane effects are included in the solid equation (30)$_1$ through the additional zeroth-order term $c_0d$. The solid is clamped at its extremities and zero traction is enforced on its top boundary.

<figure><img src="image_4.png" /><figcaption>Figure 2 Convergence histories for $E = 1.58 \cdot 10^6$.</figcaption></figure>

The fluid and solid physical parameters are $\rho^f = 1$, $\mu = 0.035$, $\rho^s = 1.1$, $L_1 = 1.15 \cdot 10^6$, $L_2 = 1.7 \cdot 10^6$ and $c_0 = 4 \cdot 10^6$. The pressure stabilization parameter in (35) is set to $\gamma_p = 10^{-3}$. The computations have been carried out with the FreeFem++ software (see [33]).

In order to illustrate the accuracy of Algorithm 2, we evaluate the relative displacement error $\|d_h^N - d(T)\|_e / \|d(T)\|_e$ at the final time $T = 0.015$ for different values of the discretization parameters under $\tau = O(h)$, namely,

$$ \tau \in \{5 \cdot 10^{-4} / 2^i\}_{i=0}^4, \quad h \in \{0.1 / 2^i\}_{i=0}^4. $$
30 Robin-Robin loose coupling for FSI

<figure><img src="image_5.png" /><figcaption>**Figure 3** Convergence histories for $E = 3.16 \cdot 10^5$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4** Convergence histories for $E = 1.58 \cdot 10^5$.</figcaption></figure>

The relative errors are calculated with a reference solution obtained with an implicit coupling scheme and a fine space-time grid ($h = 0.003125$ and $\tau = 10^{-6}$). The resulting convergence histories are reported in Figures 2–4 for three decreasing values of the Young modulus $E$. For comparison purposes, the results obtained with the implicit coupling scheme are also provided.

Several remarks are in order. The results indicate that there exists an optimal value of the Robin parameter $\alpha$ in terms of accuracy, and that accuracy tends to degrade for large or small values of $\alpha$. Figures 2–4 also show that overall first-order accuracy $O(h)$ is obtained for some values of $\alpha$. This is more noticeable in Figures 3 and 4. The results show that accuracy tends to degrade when increasing the elastic parameters $L_1$ and $L_2$. All these observations are in agreement with the error estimate provided by Theorem 2. The results of Figures 2–4 indicate also that the optimal value of the Robin parameter $\alpha$ is proportional to $\sqrt{E}$. Such scaling is consistent with the choice $\alpha = O(\sqrt{\rho^s E})$ advocated in Section 4.1, with the purpose of equilibrating the contributions in the error bound provided by Theorem 2.

## 5.3 Lid-driven cavity with flexible bottom

As second numerical example, we consider the classical shear-driven cavity problem with a flexible bottom (see, e.g., [5]). The fluid domain is the unit
square $\Omega^f = [0, 1] \times [0, 1]$ with the fluid-solid interface on the lower boundary
$\Sigma = [0, 1] \times \{0\}$. All the units are given in the SI system.

<figure><img src="image_7.png" /><figcaption>**Figure 5** Geometric description.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6** Snapshots of the fluid velocity in the deformed configuration at three different time instants obtained with Algorithm 1 and $\alpha = 1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7** Time histories of the mid-point displacement magnitude obtained with the strongly coupled scheme and with Algorithm 1 for different values of $\alpha$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8** Time histories of the interface mid-point displacement magnitude obtained the strongly coupled scheme and with Algorithm 1 for $\alpha = 1$ and different correction iterations.</figcaption></figure>

The system is described by the non-linear coupled problem (1)-(3) in which the solid model (2) is replaced by a non-linear Timoshenko beam (see, e.g., [34]). The overall system is initially at rest and an oscillating shear velocity profile is imposed on the upper boundary $u(t) = (1 - \cos(0.4\pi t), 0)^T$ on $[0, 1] \times \{1\}$. A no-slip condition is enforced on the portions $\{0, 1\} \times [0, 0.9]$ of the lateral cavity walls, and zero traction is prescribed on the remaining parts $\{0, 1\} \times [0.9, 1]$. The fluid physical parameters are given by $\rho^f = 1.0$ and $\mu = 0.01$, and for the solid $\rho^s = 250$, thickness $\epsilon = 0.002$, $E = 250$ and $\nu = 0$.

The spatial fluid approximation in (17) is made of $Q_1/Q_1$ finite elements with a SUPG/PSPG stabilized formulation. Linear MITC beam elements are used for the approximation of the Timoshenko beam in (16). In this numerical example, the computational meshes for the fluid and the solid meshes are, respectively, made of 400 quadrilaterals and of 20 segments and a time-step length of $\tau = 0.1$ is considered. The numerical computations have been performed with the FELiScE library<fnref n="1" />.

For illustration purposes, we have reported in Figure 6 some snapshots of the fluid velocity magnitude in the current configuration at different time-instants, obtained with Algorithm 1 and $\alpha = 1$. We now focus on the accuracy of this loosely coupled scheme. To this purpose, we have first reported in Figure 7 the time history of the interfacial mid-point displacement magnitude obtained with Algorithm 1, for different values of $\alpha$, and with a strongly coupled scheme. The impact of $\alpha$ on the accuracy is clearly visible. As in the previous numerical example (Section 5.2), the accuracy of Algorithm 1 degrades when $\alpha$ is both smaller and larger than the optimal value $\alpha = 1$. It is also worth noting the very good accuracy of the results obtained with Algorithm 1 and $\alpha = 1$. Indeed, the proposed loosely coupled scheme delivers practically the same accuracy as the strongly coupled scheme. This observation is in line with the results of the previous section (Figures 2–4) and the error estimate of Theorem 2, which indicate that small values of the solid elastic parameters yield a low splitting error. Furthermore, the choice $\alpha = 1$ is consistent with the scaling $\alpha = O(\epsilon\sqrt{\rho^s E})$ advocated in Section 4.1 (the factor

<footnote n="1">¹https://gitlab.inria.fr/felisce/felisce</footnote>
$\epsilon$ commes from the thin-walled nature of the solid model). Figure 8 compares the results obtained with Algorithm 1 for $\alpha = 1$ and different correction iterations $\mathcal{K}$. We recall that each correction iteration amounts to perform once more Steps 1 and 3 of Algorithm 1 with updated Robin conditions (see Appendix A). Since the accuracy is already very good without any correction, it only slight improves as we increase the number of corrections.

<table><caption>**Table 1** L²-difference between the time histories of the interfacial mid-point displacement magnitude obtained with the strongly coupled scheme and with Algorithm 1, for different values of $\alpha$ and $\mathcal{K}$.</caption><thead><tr><th>α</th><th>0.1</th><th>1</th><th>10</th></tr></thead><tbody><tr><th>0</th><td>0.0718</td><td>0.0161</td><td>0.0656</td></tr><tr><th>1</th><td>0.0419</td><td>0.0087</td><td>0.0175</td></tr><tr><th>2</th><td>0.0296</td><td>0.0030</td><td>0.0080</td></tr><tr><th>3</th><td>0.0231</td><td>0.0020</td><td>0.0038</td></tr></tbody></table>

Finally, in order to quantitatively assess the impact of the free parameters $\alpha, \mathcal{K}$ on the accuracy of the approximations, Table 1 provides the L²-difference of the time histories of the interfacial mid-point displacement magnitude obtained with the strongly coupled scheme and with Algorithm 1, for different values of $\alpha$ and $\mathcal{K}$. We can clearly observe that each correction iteration enhances the accuracy of the approximation. It is also worth noting that the number of corrections necessary to achieve a certain level of accuracy is sensitive to $\alpha$.

## 5.4 Pressure wave in an elastic tube

We consider the 3D non-linear counterpart of the numerical example considered in Section 5.2, viz., the propagation of a pressure wave within an elastic tube (see, e.g., [2, Chapter 12]). The fluid domain is a straight tube of radius 0.5 and length 5. The vessel wall has a thickness of 0.1. All the units are given in the CGS system. The fluid-structure system is modeled by the non-linear coupled problem (1)-(3), with a Saint Venant-Kirchhoff law for the solid. The overall system is initially at rest and an over pressure of $1.3332 \cdot 10^4$ is imposed on the inlet boundary during the time interval $[0, 0.005]$. The solid wall is clamped at its extremities. The physical parameters for the fluid are $\rho^f = 1$ and $\mu = 0.035$, and for the solid $\rho^s = 1.2$, $E = 3 \cdot 10^6$ and $\nu = 0.3$. In (16) and (17) continuous $\mathbb{P}_1$ finite elements are used (a SUPG/PSPG stabilized formulation is considered in the fluid). The numerical computations have been performed with the FELiScE library.

For the sake of illustration, Figure 5.4 shows some snapshots of the fluid pressure and of the deformed solid configuration at different time-instants, obtained with Algorithm 1 and $\alpha = 500$.

Figure 10 reports the time-histories of the interfacial mid-point displacement magnitude obtained with Algorithm 1 for several values of $\alpha$. The
<figure><img src="image_11.png" /><figcaption>**Figure 9** Snapshots of the fluid pressure and of the deformed solid configuration at three different time instants obtained with Algorithm 1 and $\alpha = 500$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10** Evolution of the interface mid-point displacement magnitude obtained with the strongly coupled scheme and Algorithm 1 for different values of $\alpha$.</figcaption></figure>

numerical solution provided by the strongly coupled scheme is also displayed, for comparison purposes. As observed in the previous numerical examples, the accuracy of Algorithm 1 is sensitive to the choice of $\alpha$. The best result with the proposed loosely coupled scheme is obtained for around $\alpha = 500$ and $\alpha = 1000$. These values are consistent with the optimal scaling $\alpha = O(\sqrt{\rho^s E})$ suggested in Section 4.1. It is also worth noting that, in contrast to the numerical example of Section 5.3, a significant deviation is obtained with respect to the approximation provided by the strongly coupled scheme. This observation is in agreement with the error estimate of Section 4.1 and with the numerical results of Section 5.2, which indicate that larger values of the solid elastic parameters amplify the splitting error.

Figure 11 reports the results obtained with Algorithm 1 for $\alpha = 500$ by performing a few correction iterations $K$ (see Appendix A). We can clearly observe that the accuracy of the approximation is enhanced as we increase the number of correction iterations. It is also worth noting that a single correction significantly improves the results.

As in Section 5.3, we can quantitatively assess the impact of $\alpha$, $K$ on the accuracy of the approximations by evaluating the $L^2$-difference of the time histories of the interfacial mid-point displacement magnitude obtained with the strongly coupled scheme and with Algorithm 1. The results are reported in Table 2. Each correction iteration enhances the accuracy of the approximation,
<figure><img src="image_13.png" /><figcaption>**Figure 11** Time histories of the interface mid-point displacement magnitude obtained with the strongly coupled scheme and Algorithm 1 for $\alpha = 500$ and different correction iterations.</figcaption></figure>

<table><caption>**Table 2** L²-difference between the time histories of the interfacial mid-point displacement magnitude obtained with the strongly coupled scheme and with Algorithm 1, for different values of $\alpha$ and $K$.</caption><thead><tr><th>α</th><th>125</th><th>250</th><th>500</th><th>1000</th><th>2000</th></tr></thead><tbody><tr><th>0</th><td>0.3962</td><td>0.3284</td><td>0.2700</td><td>0.3802</td><td>6.7457</td></tr><tr><th>1</th><td>0.2609</td><td>0.1990</td><td>0.0926</td><td>0.1342</td><td>4.0078</td></tr><tr><th>2</th><td>0.2361</td><td>0.1381</td><td>0.0585</td><td>0.0536</td><td>2.5329</td></tr><tr><th>3</th><td>0.1977</td><td>0.1014</td><td>0.0407</td><td>0.0224</td><td>1.6369</td></tr></tbody></table>

and this is particularly striking for the optimal values $\alpha = 500$ and $\alpha = 1000$. We can also observe that taking larger or smaller values of $\alpha$ yields larger deviations, so that further correction iterations are necessary to achieve a certain level of accuracy. Finally, it is also worth noting that the errors reported in Table 2 are larger than those of Table 1, which is also an indication that the splitting error induced by Algorithm 1 increases with increasing solid elastic parameters.

# 6 Conclusion

We have extended the stability results of the loosely coupled Robin-Robin coupling from [22] to the fully nonlinear case. A nearly-optimal error estimate was also reported for the linear case. The numerical section validated the theoretical results, but also pointed to the importance of making a judicious choice of the Robin parameter $\alpha$. The best results are obtained with the scaling $\alpha = O(\sqrt{\rho^8 E})$ suggested by the error analysis. Furthermore, large values of the elastic parameters amplify the splitting error, so that a few correction iterations might be necessary to enhance accuracy in the case of the coupling with stiffer solids.
# Declarations

## Conflict of interest

The authors declare that they have no conflict of interest.

## Data availability

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

# Acknowledgments

EB was partially supported by EPSRC grant EP/W007460/1. MF was partially supported by the French National Research Agency (ANR) through the SIMR project (ANR-19-CE45-0020). EB and MF were partially supported by the IMFIBIO associated team. OR was financially supported by the Vicerrectoria de Investigación, Desarrollo e Innovación (VRIDEI) of the Universidad de Santiago de Chile, and the National Agency for Research and Development (ANID) Doctorado Becas Chile/2020 - 72210101 of the Government of Chile.

# A Algorithm 1 with correction iterations

We introduce the following notations

$$ \partial_{\tau} x^{n,k} \stackrel{\text{def}}{=} \frac{1}{\tau}(x^{n,k} - x^{n-1}), \quad x^{n-\frac{1}{2},k} \stackrel{\text{def}}{=} \frac{1}{2}(x^{n,k} + x^{n-1}), $$

where $x^{n,k}$ denotes the approximation of $x^n$ at the correction iteration $k-1$, for $1 \le k \le K$. Assume that $u_h^0$, $d_h^0$ and $d_h^{f,0}$, $\lambda_h^0$ are given. Algorithm 1 with $K-1$ correction iterations, reads as follows for $n \ge 1$:

1. Update fluid domain:

$$ d_h^{f,n} = \mathcal{L}_h(d_h^{n-1}|_{\Sigma}), \quad w_h^n = \partial_{\tau} d_h^{f,n}, \quad \mathcal{A}_h^n = I_{\Omega^f} + d_h^{f,n}, \quad \Omega^{f,n} = \mathcal{A}_h^n(\Omega^f). $$

2. Set $(d_h^{n,0}, \dot{d}_h^{n,0}, \lambda_h^{n,0}) = (d_h^{n-1}, \dot{d}_h^{n-1}, \lambda_h^{n-1})$ and for $1 \le k \le K$ perform:

* Solve solid: find $(d_h^{n,k}, \dot{d}_h^{n,k}) \in V_h^s \times V_h^s$, with $\dot{d}_h^{n-\frac{1}{2},k} = \frac{1}{\tau}(d_h^{n,k} - d_h^{n-1})$ such that

$$ \rho^s \int_{\Omega^s} \partial_{\tau} \dot{d}_h^{n,k} \cdot \xi_h + \int_{\Omega^s} \bar{\Sigma}_h^{n-\frac{1}{2},k} : \partial_d(E_h^{n-\frac{1}{2},k}) \xi_h + \alpha \int_{\Sigma} (\dot{d}_h^{n-\frac{1}{2},k} - \hat{u}_h^{n,k-1}) \cdot \xi_h \\ = - \int_{\Sigma} \lambda_h^{n,k-1} \cdot \xi_h $$

for all $\xi_h \in V_h^s$.
• Solve fluid: find $(\hat{u}^{n,k}, \hat{p}^{n,k}) \in V_h^f \times Q_h$ such that

$$ \frac{\rho}{\tau} \left( \int_{\Omega^{f,n}} \mathbf{u}_h^{n,k} \cdot \mathbf{v}_h - \int_{\Omega^{f,n-1}} \mathbf{u}_h^{n-1} \cdot \mathbf{v}_h \right) + s_{\Omega^{f,n}} (\mathbf{u}_h^{n,k}, \mathbf{v}_h) \\ + a_{\Omega^{f,n}} (\mathbf{u}_h^{n-1}, \mathbf{w}_h^n; (\mathbf{u}_h^{n,k}, \mathbf{p}_h^{n,k}), (\mathbf{v}_h, q_h)) + \alpha \int_{\Sigma} (\hat{\mathbf{u}}_h^{n,k} - \dot{\mathbf{d}}_h^{n-\frac{1}{2},k}) \cdot \hat{\mathbf{v}}_h \\ = \int_{\Sigma} \lambda_h^{n,k-1} \cdot \hat{\mathbf{v}}_h $$

for all $(\hat{\mathbf{v}}_h, \hat{q}_h) \in V_h^f \times Q_h$.

• Fluid-stress update: set $\lambda_h^{n,k} \in V_h^g$ as

$$ \lambda_h^{n,k} = \lambda_h^{n,k-1} + \alpha (\dot{\mathbf{d}}_h^{n-\frac{1}{2},k} - \hat{\mathbf{u}}_h^{n,k}) \quad \text{on } \Sigma. $$

3. Update approximation:

$$ (\mathbf{d}_h^n, \dot{\mathbf{d}}_h^n, \hat{\mathbf{u}}^n, \hat{\mathbf{p}}^n, \lambda_h^n) = (\mathbf{d}_h^{n,\mathcal{K}}, \dot{\mathbf{d}}_h^{n,\mathcal{K}}, \hat{\mathbf{u}}^{n,\mathcal{K}}, \hat{\mathbf{p}}^{n,\mathcal{K}}, \lambda_h^{n,\mathcal{K}}). $$

Note that for $\mathcal{K} = 1$ (i.e., with 0 correction iterations) this numerical procedure is nothing but Algorithm 1. On the other hand, if we iterate in $k$ until convergence we retrieve the strongly coupled scheme.
B Proof of Lemma 3

By simply gathering the contributions of $D_h^n$ and $\Psi_h^n$, we have

$$
\begin{align}
D_h^n + \Psi_h^n \lesssim & \|e_{\mathbf{d},\Pi}^n\|_s^2 + \frac{Th^2}{\alpha} \sum_{m=1}^n \|g_2^m\|_{1,\Sigma}^2 \nonumber \\
& + \tau \sum_{m=1}^n \left[ \mu \|\epsilon(e_{\mathbf{u},\Pi}^m)\|_{0,\Omega^f}^2 + \frac{1}{\mu} \|e_{p,\Pi}^m\|_{0,\Omega^f}^2 + T \|R_h^s(\dot{\mathbf{d}}^{n-\frac{1}{2}} - \partial_\tau \mathbf{d}^m)\|_s^2 \right. \nonumber \\
& \left. + T \|\partial_\tau e_{\mathbf{d},\Pi}^m\|_s^2 + \frac{1}{T} \|e_{\mathbf{d},\Pi}^{m-\frac{1}{2}}\|_s^2 + \frac{\mu}{h^2} \left( \|e_{\mathbf{u},\Pi}^m\|_{0,\Omega^f}^2 + \|e_{\mathbf{u},\Pi}^m\|_{0,\Sigma}^2 \right) \right. \nonumber \\
& \left. + \left( \alpha + \frac{1}{\alpha} + \frac{1+\alpha^2}{\mu} \right) \|g_{1,h}^m\|_{0,\Sigma}^2 + \left( \frac{1}{\alpha} + \frac{1}{\mu} \right) \|g_2^m\|_{0,\Sigma}^2 \right. \nonumber \\
& + \rho^s T \|\partial_t \dot{\mathbf{d}}^{m-\frac{1}{2}} - \partial_\tau R_h^s \dot{\mathbf{d}}^m\|_{0,\Omega^s}^2 + \rho^f T \|\partial_t \mathbf{u}^m - \partial_\tau R_h^f \mathbf{u}^m\|_{0,\Omega^f}^2 \nonumber \\
& + \frac{h^2}{\mu} \|\nabla S_h p^m\|_{0,\Omega^f}^2 + \|g_2^m\|_{L^\infty(\Sigma)}^2 \nonumber \\
& + \rho^s T \tau^2 \left\| \partial_\tau^2 R_h^s \dot{\mathbf{d}}^m \right\|_{0,\Omega^s}^2 + \frac{\rho^s T \tau^2}{\alpha^2} \left\| \partial_\tau^2 \mathcal{L}_{s,h}^n \right\|_{0,\Omega^s}^2 \nonumber \\
& + \frac{T}{\alpha^2} \left\| \tilde{g}_{5,h}^m \right\|_s^2 + \frac{1}{\alpha} \left\| \tilde{g}_{5,h}^m \right\|_{0,\Sigma}^2 + \frac{\rho^s}{T\alpha^2} \left\| \tilde{g}_{5,h}^m \right\|_{0,\Omega^s}^2 \right] + \frac{\rho^s}{\alpha^2} \left\| \tilde{g}_{5,h}^n \right\|_{0,\Omega^s}^2. \tag{63}
\end{align}
$$

We now proceed by estimating each term separately.

Using (33)–(34) and the approximation properties of $R_h^s$, we have

$$
\|e_{\mathbf{d},\Pi}^n\|_s^2 \lesssim (L_1 + L_2)h^2 \|\mathbf{d}^n\|_{2,\Omega^s}^2 \le (L_1 + L_2)h^2 \|\mathbf{d}\|_{L^\infty(0,T; H^2(\Omega^s))}^2
$$

For the second term, we have

$$
\begin{align*}
\frac{Th^2}{\alpha} \sum_{m=1}^{n} \|g_2^m\|_{1,\Sigma}^2 &\le \frac{Th^2\tau}{\alpha} \|\partial_t \lambda\|_{L^2(0,T; H^1(\Sigma))}^2 \\
&\lesssim \frac{(h^4 + \tau^2)T(L_1 + L_2)^2}{\alpha} \|\nabla \dot{\mathbf{d}}\|_{L^2(0,T; H^1(\Sigma))}^2.
\end{align*}
$$

Similarly, for the next two terms, we have

$$
\mu\tau \sum_{m=1}^{n} \|\epsilon(e_{\mathbf{u},\Pi}^{m})\|_{0,\Omega^f}^2 \lesssim \mu h^2 \tau \sum_{m=1}^{n} \|\mathbf{u}^n\|_{2,\Omega^f}^2 \le \mu h^2 T \|\mathbf{u}\|_{L^\infty(0,T; H^2(\Omega^f))}^2
$$

and

$$
\frac{\tau}{\mu} \sum_{m=1}^{n} \|e_{p,\Pi}^{m}\|_{0,\Omega^{f}}^{2} \lesssim \frac{h^{2}\tau}{\mu} \sum_{m=1}^{n} \|p^{n}\|_{1,\Omega^{f}}^{2} \le \frac{h^{2}T}{\mu} \|p\|_{L^{\infty}(0,T; H^{1}(\Omega^{f}))}^{2}.
$$
On the other hand, using the $H^1$-stability of $R_h^s$ and a Taylor expansion, we have

$$ 
\begin{aligned}
\tau T \sum_{m=1}^{n} \left( \| R_h^s (\dot{\mathbf{d}}^{m-\frac{1}{2}} - \partial_{\tau} \mathbf{d}^m) \|_{\mathrm{s}}^2 + \| \partial_{\tau} \mathbf{e}_{\mathbf{d},\Pi}^m \|_{\mathrm{s}}^2 \right) & \\
& \lesssim \tau T \sum_{m=1}^{n} \left( \| \partial_t \mathbf{d}^{m-\frac{1}{2}} - \partial_{\tau} \mathbf{d}^m \|_{\mathrm{s}}^2 + \| \partial_{\tau} \mathbf{e}_{\mathbf{d},\Pi}^m \|_{\mathrm{s}}^2 \right) \\
& \lesssim T(L_1 + L_2) \left( \tau^4 \| \partial_t^3 \mathbf{d} \|_{L^2(0,T; H^1(\Omega^s))}^2 + h^2 \| \dot{\mathbf{d}} \|_{L^2(0,T; H^2(\Omega^s))}^2 \right).
\end{aligned}
$$

Similarly,

$$ \frac{\tau}{T} \sum_{m=1}^{n} \| e_{\mathbf{d},\Pi}^{m-\frac{1}{2}} \|_{\mathrm{s}}^2 \lesssim \frac{\tau}{T} \sum_{m=1}^{n} \| e_{\mathbf{d},\Pi}^{m} \|_{\mathrm{s}}^2 + \frac{\tau}{T} \| e_{\mathbf{d},\Pi}^{0} \|_{\mathrm{s}}^2 \lesssim (L_1+L_2)h^2 \| \mathbf{d} \|_{L^{\infty}(0,T; H^2(\Omega^s))}^2. $$

Using the approximation properties of $R_h^f$, we have

$$ \frac{\mu \tau}{h^2} \sum_{m=1}^{n} \left( \| e_{\mathbf{u},\Pi}^{m} \|_{0,\Omega^f}^2 + \| e_{\mathbf{u},\Pi}^{m} \|_{0,\Sigma} \right) \lesssim \mu T h^2 \left( \| \mathbf{u} \|_{L^{\infty}(0,T; H^2(\Omega^f))}^2 + \| \mathbf{u} \|_{L^{\infty}(0,T; H^3(\Omega^f))}^2 \right). $$

On the other hand,

$$ 
\begin{aligned}
\left( \alpha + \frac{1}{\alpha} + \frac{1+\alpha^2}{\mu} \right) \tau \sum_{m=1}^{n} \| g_{1,h}^{m} \|_{0,\Sigma}^2 & \\
& \lesssim \left( \alpha + \frac{1}{\alpha} + \frac{1+\alpha^2}{\mu} \right) \tau \sum_{m=1}^{n} \left( \| e_{\mathbf{u},\Pi}^{n} \|_{0,\Sigma}^2 + \| e_{\mathbf{u},\Pi}^{n-1} \|_{0,\Sigma}^2 + \| \mathbf{u}^n - \mathbf{u}^{n-1} \|_{0,\Sigma}^2 \right) \\
& \lesssim \left( \alpha + \frac{1}{\alpha} + \frac{1+\alpha^2}{\mu} \right) \left( h^2 T \| \mathbf{u} \|_{L^{\infty}(0,T; H^3(\Omega^f))}^2 + \tau^2 \| \partial_t \mathbf{u} \|_{L^2(0,T; L^2(\Sigma))}^2 \right).
\end{aligned}
$$

Similarly,

$$ 
\begin{aligned}
\left( \frac{1}{\alpha} + \frac{1}{\mu} \right) \tau \sum_{m=1}^{n} \| g_2^m \|_{0,\Sigma}^2 & \lesssim \left( \frac{1}{\alpha} + \frac{1}{\mu} \right) \tau^2 \| \partial_t \boldsymbol{\lambda} \|_{L^2(0,T; L^2(\Sigma))}^2 \\
& \lesssim \left( \frac{1}{\alpha} + \frac{1}{\mu} \right) (L_1 + L_2)^2 \tau^2 \| \nabla \dot{\mathbf{d}} \|_{L^2(0,T; L^2(\Sigma))}^2.
\end{aligned}
$$

As above, we also have

$$ \rho^s T \tau \sum_{m=1}^{n} \| \partial_t \dot{\mathbf{d}}^{m-\frac{1}{2}} - \partial_{\tau} R_h^s \dot{\mathbf{d}}^m \|_{0,\Omega^s}^2 $$
40 Robin-Robin loose coupling for FSI

$$
\begin{aligned}
&\lesssim \rho^s T \tau \sum_{m=1}^{n} \left( \| \partial_t \dot{\mathbf{d}}^{m-\frac{1}{2}} - \partial_\tau \dot{\mathbf{d}}^m \|_{0,\Omega^s}^2 + \| \partial_\tau e_{\dot{\mathbf{d}},\Pi}^m \|_{0,\Omega^s}^2 \right) \\
&\lesssim \rho^s T \left( \tau^4 \| \partial_t^3 \dot{\mathbf{d}} \|_{L^2(0,T; L^2(\Omega^s))}^2 + h^2 \| \partial_t \dot{\mathbf{d}} \|_{L^2(0,T; H^1(\Omega^s))}^2 \right)
\end{aligned}
$$

and

$$
\begin{aligned}
& \rho^f T \tau \sum_{m=1}^{n} \| \partial_t \mathbf{u}^m - \partial_\tau R_h^f \mathbf{u}^m \|_{0,\Omega^f}^2 \\
& \lesssim \rho^f T \left( \tau^4 \| \partial_t^3 \mathbf{u} \|_{L^2(0,T; L^2(\Omega^f))}^2 + h^2 \| \partial_t \mathbf{u} \|_{L^2(0,T; H^1(\Omega^f))}^2 \right).
\end{aligned}
$$

Using the $H^1$-stability of $S_h$, we have

$$
\frac{h^2 \tau^2}{\mu} \sum_{m=1}^{n} \| \nabla S_h p^m \|_{0, \Omega^f}^2 \lesssim \frac{h^2 T}{\mu} \| p \|_{L^{\infty}(0, T; H^1(\Omega^f))}^2.
$$

Similarly,

$$
\tau \sum_{m=1}^{n} \| \boldsymbol{g}_2^m \|_{L^\infty(\Sigma)}^2 \lesssim \tau^2 \| \partial_t \boldsymbol{\lambda} \|_{L^2(0,T; L^\infty(\Sigma))}^2 \lesssim \tau^2 (L_1 + L_2)^2 \| \nabla \dot{\boldsymbol{d}} \|_{L^2(0,T; L^\infty(\Sigma))}^2.
$$

It now remains to control the terms with the extensions onto $\Omega^s$. For the first term, we use the stability of $R_h^s$:

$$
\rho^s T \tau^3 \sum_{m=1}^{n} \| \partial_{\tau}^2 R_h^s \dot{\mathbf{d}}^m \|_{0, \Omega^s}^2 \lesssim \rho^s T \tau^2 \| \partial_t^2 \dot{\mathbf{d}} \|_{L^2(0, T; L^2(\Omega^s))}^2.
$$

On the other hand, using Assumptions 1 and 2, we have

$$
\begin{aligned}
& \frac{\rho^s T \tau^3}{\alpha^2} \sum_{m=1}^{n} \| \partial_{\tau}^2 \mathcal{L}_{s,h}^n \|_{0,\Omega^s}^2 \lesssim \frac{\rho^s T \tau^3}{\alpha^2} \sum_{m=1}^{n} \| \partial_{\tau}^2 \mathcal{L}_{s}^n \|_{0,\Omega^s}^2 \\
& \lesssim \frac{\rho^s T \tau^2}{\alpha} \| \tilde{\phi} \tilde{\sigma} (\partial_t^2 \mathbf{d}) \tilde{n}^s \|_{L^2(0,T; L^2(\Omega^s))}^2 \lesssim \frac{\rho^s T \tau^2}{\alpha} \| \partial_t^2 \mathbf{d} \|_{L^2(0,T; H^1(\Omega^s))}^2.
\end{aligned}
$$

For the next term,

$$
\frac{T \tau}{\alpha^2} \sum_{m=1}^{n} \| \tilde{\mathbf{g}}_{5,h}^m \|_{\mathrm{s}}^2 \lesssim \underbrace{T \tau \sum_{m=1}^{n} \| \tilde{\mathbf{g}}_{1,h}^m \|_{\mathrm{s}}^2}_{T_1} + \underbrace{\frac{T \tau}{\alpha^2} \sum_{m=1}^{n} \| \tilde{\mathbf{g}}_{2,h}^m \|_{\mathrm{s}}^2}_{T_2}.
$$

For the first term, using the $H^1$-stability of $R_h^s$, we have
$$
\begin{aligned}
T_1 &= T\tau \sum_{m=1}^{n} \|R_h^s(\dot{\mathbf{d}}^n - \dot{\mathbf{d}}^{n-1})\|_s^2 \lesssim T\tau(L_1 + L_2) \sum_{m=1}^{n} \|\dot{\mathbf{d}}^n - \dot{\mathbf{d}}^{n-1}\|_{1,\Omega^s}^2 \\
&\lesssim T\tau^2(L_1 + L_2) \|\partial_t \dot{\mathbf{d}}\|_{L^2(0,T; H^1(\Omega^s))}^2.
\end{aligned}
$$

On the other hand, for the second, using Assumptions 1 and 2, we have

$$
\begin{aligned}
T_2 &\lesssim \frac{T\tau}{\alpha^2}(L_1 + L_2) \sum_{m=1}^{n} \|\nabla \tilde{g}_2^n\|_{0,\Omega^s}^2 \\
&\lesssim \frac{T\tau}{\alpha^2}(L_1 + L_2) \sum_{m=1}^{n} \left( \|\nabla \tilde{\phi}\|_{0,\Omega^s}^2 \|\tilde{\ell}_s^m - \tilde{\ell}_s^{m-1}\|_{L^\infty(\Omega^s)}^2 \right. \\
&\qquad \left. + \|\tilde{\phi}\|_{L^\infty(\Omega^s)}^2 \|\nabla(\tilde{\ell}_s^m - \tilde{\ell}_s^{m-1})\|_{0,\Omega^s}^2 \right) \\
&\lesssim \frac{T\tau}{\alpha^2}(L_1 + L_2) \sum_{m=1}^{n} \left( \|\nabla \tilde{\phi}\|_{0,\Omega^s}^2 \|\boldsymbol{\sigma}(\mathbf{d}^m - \mathbf{d}^{m-1})\|_{L^\infty(\Omega^s)}^2 \right. \\
&\qquad \left. + \|\nabla \tilde{n}^s\|_{L^\infty(\Omega^s)}^2 \|\boldsymbol{\sigma}(\mathbf{d}^m - \mathbf{d}^{m-1})\|_{0,\Omega^s}^2 + \|\nabla \boldsymbol{\sigma}(\mathbf{d}^m - \mathbf{d}^{m-1})\|_{0,\Omega^s}^2 \right) \\
&\lesssim \frac{T\tau}{\alpha^2}(L_1 + L_2)^3 \sum_{m=1}^{n} \left( \|\nabla \tilde{\phi}\|_{0,\Omega^s}^2 \|\mathbf{d}^m - \mathbf{d}^{m-1}\|_{W^{1,\infty}(\Omega^s)}^2 + \|\mathbf{d}^m - \mathbf{d}^{m-1}\|_{2,\Omega^s}^2 \right) \\
&\lesssim \frac{T\tau^2}{\alpha^2}(L_1 + L_2)^3 \left( \|\nabla \tilde{\phi}\|_{0,\Omega^s}^2 \|\dot{\mathbf{d}}\|_{L^2(0,T; W^{1,\infty}(\Omega^s))}^2 + \|\dot{\mathbf{d}}\|_{L^2(0,T; H^2(\Omega^s))}^2 \right) \\
&\lesssim \frac{T\tau^2}{\alpha^2}(L_1 + L_2)^3 \left[ \left(1 + \log \frac{1}{\tau}\right) \|\dot{\mathbf{d}}\|_{L^2(0,T; W^{1,\infty}(\Omega^s))}^2 + \|\dot{\mathbf{d}}\|_{L^2(0,T; H^2(\Omega^s))}^2 \right].
\end{aligned}
$$

Similarly, we have

$$
\begin{aligned}
\frac{\rho^s \tau}{T \alpha^2} \sum_{m=1}^{n} \|\tilde{g}_{5,h}^m\|_{0,\Omega^s}^2 &\lesssim \frac{\rho^s \tau}{T} \sum_{m=1}^{n} \|\tilde{g}_{1,h}^m\|_{0,\Omega^s}^2 + \frac{\rho^s \tau}{T \alpha^2} \sum_{m=1}^{n} \|\tilde{g}_{2,h}^m\|_{0,\Omega^s}^2 \\
&\lesssim \frac{\rho^s \tau}{T} \sum_{m=1}^{n} \|\dot{\mathbf{d}}^m - \dot{\mathbf{d}}^{m-1}\|_{0,\Omega^s}^2 + \frac{\rho^s \tau}{T \alpha^2} \sum_{m=1}^{n} \|\boldsymbol{\sigma}(\mathbf{d}^m - \mathbf{d}^{m-1})\|_{0,\Omega^s}^2 \\
&\lesssim \frac{\rho^s \tau^2}{T} \|\partial_t \dot{\mathbf{d}}\|_{L^2(0,T; L^2(\Omega^s))}^2 + \frac{\rho^s (L_1 + L_2)^2 \tau^2}{T \alpha^2} \|\dot{\mathbf{d}}\|_{L^2(0,T; H^1(\Omega^s))}^2,
\end{aligned}
$$

and

$$
\begin{aligned}
\frac{\rho^s}{\alpha^2} \|\tilde{g}_{5,h}^n\|_{0,\Omega^s}^2 &\lesssim \rho^s \|\tilde{g}_{1,h}^m\|_{0,\Omega^s}^2 + \frac{\rho^s}{\alpha^2} \|\tilde{g}_{2,h}^m\|_{0,\Omega^s}^2 \\
&\lesssim \rho^s \tau^2 \|\partial_t \dot{\mathbf{d}}\|_{L^\infty(0,T; L^2(\Omega^s))}^2 + \frac{\rho^s (L_1 + L_2)^2 \tau^2}{\alpha^2} \|\dot{\mathbf{d}}\|_{L^\infty(0,T; H^1(\Omega^s))}^2.
\end{aligned}
$$
42 Robin-Robin loose coupling for FSI

Finally,

$$
\begin{aligned}
\frac{\tau}{\alpha} \sum_{m=1}^{n} \| \tilde{g}_{5,h}^{m} \|_{0,\Sigma}^{2} &\lesssim \tau \alpha \sum_{m=1}^{n} \| R_{h}^{f}(\boldsymbol{u}^{n} - \boldsymbol{u}^{n-1}) \|_{0,\Sigma}^{2} + \frac{\tau}{\alpha} \sum_{m=1}^{n} \| g_{2}^{m} \|_{0,\Sigma}^{2} \\
&\lesssim \alpha \left( h^{2} T \| \boldsymbol{u} \|_{L^{\infty}(0,T; H^{3}(\Omega^{f}))}^{2} + \tau^{2} \| \partial_{t} \boldsymbol{u} \|_{L^{2}(0,T; L^{2}(\Sigma))}^{2} \right) \\
&\qquad + \frac{\tau^{2}(L_{1} + L_{2})^{2}}{\alpha} \| \boldsymbol{\nabla} \dot{\boldsymbol{d}} \|_{L^{2}(0,T; L^{2}(\Sigma))}^{2}.
\end{aligned}
$$

The bound (62) then follows by inserting all the above estimates into (63), which completes the proof.

References

[1] Fernández, M.A., Gerbeau, J.-F.: Algorithms for fluid-structure interaction problems. In: Cardiovascular Mathematics. MS&A. Model. Simul. Appl., vol. 1, pp. 307–346. Springer, ??? (2009)

[2] Formaggia, L., Quarteroni, A., Veneziani, A. (eds.): Cardiovascular Mathematics. Modeling and Simulation of the Circulatory System. Modeling, Simulation and Applications, vol. 1. Springer, ??? (2009)

[3] Farhat, C., Lesoinne, M., LeTallec, P.: Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: momentum and energy conservation, optimal discretization and application to aeroelasticity. Comput. Methods Appl. Mech. Engrg. 157(1-2), 95–114 (1998). https://doi.org/10.1016/S0045-7825(97)00216-8

[4] Causin, P., Gerbeau, J.F., Nobile, F.: Added-mass effect in the design of partitioned algorithms for fluid-structure problems. Comput. Methods Appl. Mech. Engrg. 194(42-44), 4506–4527 (2005). https://doi.org/10.1016/j.cma.2004.12.005

[5] Förster, C., Wall, W.A., Ramm, E.: Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows. Comput. Methods Appl. Mech. Engrg. 196(7), 1278–1293 (2007)

[6] Burman, E., Fernández, M.A.: Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility. Comput. Methods Appl. Mech. Engrg. 198(5-8), 766–784 (2009). https://doi.org/10.1016/j.cma.2008.10.012

[7] Hansbo, P., Hermansson, J., Svedberg, T.: Nitsche’s method combined
with space-time finite elements for ALE fluid-structure interaction prob-
lems. Comput. Methods Appl. Mech. Engrg. **193**(39-41), 4195–4206
(2004). https://doi.org/10.1016/j.cma.2003.09.029

[8] Burman, E., Fernández, M.A.: Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin-Robin coupling. Internat. J. Numer. Methods Engrg. **97**(10), 739–758 (2014). https://doi.org/10.1002/nme.4607

[9] Badia, S., Nobile, F., Vergara, C.: Fluid-structure partitioned procedures based on Robin transmission conditions. J. Comput. Phys. **227**(14), 7027–7051 (2008)

[10] Badia, S., Nobile, F., Vergara, C.: Robin-Robin preconditioned Krylov methods for fluid-structure interaction problems. Comput. Methods Appl. Mech. Engrg. **198**(33-36), 2768–2784 (2009). https://doi.org/10.1016/j.cma.2009.04.004

[11] Bukač, M., Fu, G., Seboldt, A., Trenchea, C.: Time-adaptive partitioned method for fluid-structure interaction problems with thick structures. J. Comput. Phys. **473**, 111708–19 (2023)

[12] Banks, J.W., Henshaw, W.D., Schwendeman, D.W.: An analysis of a new stable partitioned algorithm for FSI problems. Part I: Incompressible flow and elastic solids. J. Comput. Phys. **269**, 108–137 (2014). https://doi.org/10.1016/j.jcp.2014.03.006

[13] Fernández, M.A., Mullaert, J., Vidrascu, M.: Generalized Robin-Neumann explicit coupling schemes for incompressible fluid-structure interaction: stability analysis and numerics. Internat. J. Numer. Methods Engrg. **101**(3), 199–229 (2015)

[14] Bukač, M., Čanić, S., Glowinski, R., Muha, B., Quaini, A.: A modular, operator-splitting scheme for fluid-structure interaction problems with thick structures. Internat. J. Numer. Methods Fluids **74**(8), 577–604 (2014). https://doi.org/10.1002/fld.3863

[15] Serino, D.A., Banks, J.W., Henshaw, W.D., Schwendeman, D.W.: A stable added-mass partitioned (amp) algorithm for elastic solids and incompressible flow: Model problem analysis. SIAM Journal on Scientific Computing **41**(4), 2464–2484 (2019)

[16] Gigante, G., Vergara, C.: On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction. Comput. Math. Appl. **96**, 109–119 (2021). https://doi.org/10.1016/j.camwa.2021.05.012
44 Robin-Robin loose coupling for FSI

[17] Gigante, G., Vergara, C.: On the choice of interface parameters in Robin–Robin loosely coupled schemes for fluid–structure interaction. Fluids **6**(6) (2021)

[18] Bucelli, M., Gabriel, M.G., Quarteroni, A., Gigante, G., Vergara, C.: A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction. J. Comput. Phys. **490**, 112326–24 (2023)

[19] Fernández, M.A., Mullaert, J.: Convergence and error analysis for a class of splitting schemes in incompressible fluid-structure interaction. IMA J. Numer. Anal. **36**(4), 1748–1782 (2016). https://doi.org/10.1093/imanum/drv055

[20] Burman, E., Durst, R., Guzmán, J.: Stability and error analysis of a splitting method using Robin-Robin coupling applied to a fluid-structure interaction problem. Numer. Methods Partial Differential Equations **38**(5), 1396–1406 (2022). https://doi.org/10.1002/num.22840

[21] Seboldt, A., Bukač, M.: A non-iterative domain decomposition method for the interaction between a fluid and a thick structure. Numer. Methods Partial Differential Equations **37**(4), 2803–2832 (2021)

[22] Burman, E., Durst, R., Fernández, M.A., Guzmán, J.: Fully discrete loosely coupled Robin-Robin scheme for incompressible fluid-structure interaction: stability and error analysis. Numer. Math. **151**(4), 807–840 (2022). https://doi.org/10.1007/s00211-022-01295-y

[23] Burman, E., Durst, R., Fernández, M.A., Guzmán, J.: Loosely coupled, non-iterative time-splitting scheme based on Robin-Robin coupling: unified analysis for parabolic/parabolic and parabolic/hyperbolic problems. J. Numer. Math. **31**(1), 59–77 (2023). https://doi.org/10.1515/jnma-2021-0119

[24] Durst, R.: Recent advances in splitting methods based on robin-robin coupling conditions. PhD thesis, Brown University (2022). https://repository.library.brown.edu/studio/item/bdr:v34deskr/

[25] Gonzalez, O.: Exact energy and momentum conserving algorithms for general models in nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering **190**(13), 1763–1783 (2000)

[26] Tezduyar, T.E.: Stabilized finite element formulations for incompressible flow computations. Advances in Applied Mechanics, vol. 28, pp. 1–44. Elsevier (1991). https://doi.org/10.1016/S0065-2156(08)70153-4

[27] Burman, E., Fernández, M.A., Hansbo, P.: Continuous interior penalty finite element method for Oseen's equations. SIAM J. Numer. Anal. **44**(3),
1248–1274 (2006)

[28] Smaldone, S.: Numerical analysis and simulations of coupled problems for the cardiovascular system. PhD thesis, Université Pierre et Marie Curie (2014). https://tel.archives-ouvertes.fr/tel-01287506

[29] Farhat, C., Geuzaine, P., Grandmont, C.: The discrete geometric conservation law and the nonlinear stability of ale schemes for the solution of flow problems on moving grids. Journal of Computational Physics 174(2), 669–694 (2001)

[30] Brezzi, F., Pitkäranta, J.: On the stabilization of finite element approximations of the Stokes equations. In: Efficient Solutions of Elliptic Systems (Kiel, 1984). Notes Numer. Fluid Mech., vol. 10, pp. 11–19. Vieweg, ??? (1984)

[31] Scott, L.R., Zhang, S.: Finite element interpolation of nonsmooth functions satisfying boundary conditions. Mathematics of Computation 54(190), 483–493 (1990)

[32] Logg, A., Mardal, K.-A., Wells, G.N., et al.: Automated Solution of Differential Equations by the Finite Element Method. Springer, ??? (2012). https://doi.org/10.1007/978-3-642-23099-8

[33] Hecht, F.: New development in FreeFem++. J. Numer. Math. 20(3-4), 251–265 (2012)

[34] Bathe, K.J.: Finite Element Procedures. Prentice Hall, ??? (1996)