# An unfitted mesh semi-implicit coupling scheme for fluid-structure interaction with immersed solids

Miguel Angel Fernández, Fannie Gerosa

► To cite this version:

Miguel Angel Fernández, Fannie Gerosa. An unfitted mesh semi-implicit coupling scheme for fluid-structure interaction with immersed solids. International Journal for Numerical Methods in Engineering, 2021, 122, pp.5384-5408. 10.1002/nme.6449. hal-02288723

HAL Id: hal-02288723

https://inria.hal.science/hal-02288723

Submitted on 15 Sep 2019

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE TYPE

# An unfitted mesh semi-implicit coupling scheme for fluid-structure interaction with immersed solids

Miguel A. Fernández<fnref n="1" />,<fnref n="2" /> | Fannie M. Gerosa<fnref n="1" />,<fnref n="2" />

Correspondence

*Email: miguel.fernandez@inria.fr

Summary

Unfitted mesh finite element approximations of immersed incompressible fluid-
structure interaction problems which efficiently avoid strong coupling without com-
promising stability and accuracy are rare in the literature. Moreover, most of the
existing approaches introduce additional unknowns or are limited by penalty terms
which yield ill conditioning issues. In this paper, we introduce a new unfitted
mesh semi-implicit coupling scheme which avoids these issues. To this purpose, we
provide a consistent generalization of the projection based semi-implicit coupling
paradigm of [Int. J. Num. Meth. Engrg.,69(4):794-821, 2007] to the unfitted mesh
Nitsche-XFEM framework.

KEYWORDS:

Fluid-structure interaction, semi-implicit coupling scheme, projection method, unfitted meshes, Nitsche-
XFEM.

1 | INTRODUCTION

Numerical methods for the approximation of mathematical models describing the mechanical interaction of an incompressible viscous fluid with an immersed elastic structure are an essential ingredient in the computer simulation of many living and engineering systems (see, e.g.,<fnref n="1" />,<fnref n="2" />). These coupled problems often feature large interface displacements, with potential contact between solids, so that the favored numerical approaches are mainly based on unfitted mesh approximations (the fluid mesh is not fitted to the fluid-solid interface). Among these methods, the most popular are the immersed boundary method (see, e.g.,³,⁴,⁵) and the fictitious domain method (see, e.g.,⁶,⁷,⁸,⁹,¹⁰,¹¹), which treat the solid in Lagrangian form. We can also mention the methods based on fully Eulerian descriptions of the coupled system (see, e.g.,¹²,¹³).

In general, unfitted mesh methods have the reputation of being inaccurate in space. This is due to the approximation of the interface conditions in an unfitted framework and to the fact that the fluid spatial discretization does not generally allows for discontinuities across the interface (which often yields severe interfacial mass loss). Mesh adaptation can improve the situation (see, e.g.,¹⁴), but it does not cure the problem. The extended-FEM (XFEM) method, which combines a local enrichment with a cut-FEM approach (see, e.g.,¹⁵,¹⁶,¹⁷), fixes these issues but at the expense of introducing Lagrange multipliers (additional unknowns) and deteriorating the robustness (ill-conditioning). The Nitsche-XFEM method (see¹⁸,¹⁹,²⁰) circumvents these difficulties through a Nitsche's treatment of the interface coupling (with overlapping meshes) and the addition of suitable stabilization in the vicinity of the interface. The superior accuracy of Nitsche-XFEM with respect to the traditional immersed boundary or fictitious domain methods (see²¹ for a recent comparative study) comes, however, at the price of a much more involved computer implementation and a superior computational complexity. The latter is particularly due to the fact that accurate time-splitting schemes for Nitsche-XFEM are mainly of strongly coupled nature.

<footnote n="1">¹Inria, 75012 Paris, France</footnote>

<footnote n="2">²Sorbonne Université & CNRS, UMR 7598
LJLL, 75005 Paris, France</footnote>
Time splitting is generally difficult to marry with unfitted meshes without compromising stability and/or accuracy. This is a direct consequence of the weak treatment of the kinematic interface coupling. To the best of our knowledge, the sole available approaches are the splitting methods introduced in$^{5,22,23,24}$ for the immersed boundary or fictitious domain methods, and in$^{18,19,25,26}$ for unfitted Nitsche based methods. The loosely coupled schemes reported in$^{5,18,19,25,23}$ are known to enforce severe time-step restrictions for stability/accuracy or to be sensitive to the amount of added-mass effect. In the case of the coupling with thin-walled solids, these issues are circumvented by the semi-implicit and loosely coupled schemes reported in$^{22,19,26}$ and in$^{24}$, respectively. Nevertheless, the former introduces additional unknowns in the fluid sub-problem (intermediate solid velocity) and the accuracy of the latter relies on a grad-div penalty term (for enhanced mass conservation) which spoils the conditioning of the fluid problem.

In this paper, we introduce and analyze a new semi-implicit coupling scheme for unfitted mesh approximations of fluid-structure interaction problems with immersed solids which overcomes the above mentioned drawbacks. To this purpose, we propose to generalize the projection based semi-implicit splitting paradigm reported in$^{27}$ with fitted meshes, to the unfitted Nitsche based framework of$^{18,19}$. The basic idea consists in the explicit treatment of the geometrical non-linearities, convective and viscous fluid contributions (which avoids strong coupling), whereas the remain fluid pressure and solid contributions are coupled in a fully implicit fashion (which guarantees added-mass free stability). In contrast to alternative immersed boundary and fictitious domain methods involving fractional-step time-marching in the fluid (see, e.g.,$^{4,28}$), the Nitsche-XFEM approximation guarantees the spatial consistency of the Laplacian operator in the projection step. For a model problem with static interface, we prove a stability result which states that the conditionally stability of the coupling scheme in the energy norm. Numerical evidence in a series of well-known two-dimensional examples, involving large interface displacements and solid contact, highlights the stability and accuracy properties of the proposed method.

The rest of the paper is organized as follows. Section 2 presents the derivation of the proposed coupling scheme in a linear setting with static interfaces. The energy stability of this method is addressed in Section 2.3. In Section 3, the coupling scheme is formulated within a fully non-linear setting involving dynamic interfaces. The numerical experiments are reported in Section 4. Finally, a summary of the results of the present work are discussed in Section 5. Through this paper and without loss of generality, the solid is assumed to be thin-walled, which corresponds to the most difficult case. The methods and theoretical results presented in this paper remain valid in the case of the coupling with a thick-walled solid, by simply limiting the fluid problem to a single side of the interface.

## 2 | LINEAR MODEL PROBLEM: STATIC INTERFACES

We first consider a linear fluid-structure interaction problem in which the fluid is described by the Stokes equations in a fixed domain and the structure by a linear immersed thin-walled solid model. We denote by $\Sigma \subset \mathbb{R}^d$, with $d = 2, 3$, the reference configuration of the solid mid-surface. The structure is supposed to be immersed within a fixed domain $\Omega \subset \mathbb{R}^d$, with boundary $\Gamma = \partial\Omega$ (see Figure 1). In this section, we assume that the solid undergoes infinitesimal displacements so that the fluid flows within the fixed domain $\Omega^f = \Omega \setminus \Sigma \subset \mathbb{R}^d$. The immersed interface $\Sigma$ is supposed to divide $\Omega^f$ into two subdomais $\Omega_1^f = \Omega_1 \cup \Omega_2^f$, with respective unit normals $n_1 = n_\Sigma$ and $n_2 = -n_\Sigma$. Here, $n_\Sigma$ the normal unit vector given by the orientation of the surface $\Sigma$. For a given field $f$ defined in $\Omega^f$ (possibly discontinuous across the interface), we can then define its sided-restrictions, denoted

<figure><img src="image_2.png" /><figcaption>**Figure 1** Geometric configuration of the fluid domain and the immersed solid.</figcaption></figure>
by $f_1$ and $f_2$, as

$$f_1(\mathbf{x}) \stackrel{\text{def}}{=} \lim_{\xi \to 0^-} f(\mathbf{x} + \xi \mathbf{n}_1), \quad f_2(\mathbf{x}) \stackrel{\text{def}}{=} \lim_{\xi \to 0^-} f(\mathbf{x} + \xi \mathbf{n}_2),$$

for all $\mathbf{x} \in \Sigma$, and the following jump and average operators across $\Sigma$:

$$[f] \stackrel{\text{def}}{=} f_1 - f_2 \quad [fn] \stackrel{\text{def}}{=} f_1 n_1 + f_2 n_2, \quad [f] \stackrel{\text{def}}{=} \frac{1}{2}(f_1 + f_2).$$

In this framework, the considered coupled problem reads as follows: find the fluid velocity and pressure $\mathbf{u} : \Omega \times \mathbb{R}^+ \to \mathbb{R}^d$, $p : \Omega \times \mathbb{R}^+ \to \mathbb{R}$, the solid displacement and velocity $\mathbf{d} : \Sigma \times \mathbb{R}^+ \to \mathbb{R}^d$, $\dot{\mathbf{d}} : \Sigma \times \mathbb{R}^+ \to \mathbb{R}^d$ such that for all $t \in \mathbb{R}^+$ we have

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

$$\left\{ \begin{array}{ll} \rho^s \epsilon \partial_t \dot{\mathbf{d}} + L \mathbf{d} = \mathbf{T} & \text{in } \Sigma, \\ \dot{\mathbf{d}} = \partial_t \mathbf{d} & \text{in } \Sigma, \end{array} \right. \qquad (2)$$

$$\mathbf{d} = \mathbf{0} \quad \text{on} \quad \partial\Sigma,$$

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

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

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

where $\mu$ denotes the fluid dynamic viscosity. The operator $L$ describes the elastic behavior of the solid. The relations in (3) enforce, respectively, the kinematic and dynamic interface coupling conditions. Note that the former enforces two conditions since it has to be seen as $\mathbf{u}_1 = \mathbf{u}_2 = \dot{\mathbf{d}}$ on $\Sigma$.

## 2.1 | Time-discretization: semi-implicit coupled scheme

In what follows, we will use the following notation for the first-order backward difference: $\partial_{\tau} x^n \stackrel{\text{def}}{=} (x^n - x^{n-1})/\tau$, where $\tau > 0$ denotes the time-step length. For the time discretization of the coupled problem (1)-(3) we consider the projection based semi-implicit splitting scheme proposed in$^{27,29}$ for the case of fitted-mesh spatial approximations (see also$^{30,31,32,33}$). The scheme avoids strong coupling without compromising stability and accuracy. The fundamental idea consists in combining a fractional-step time-marching in the fluid with a semi-implicit treatment of the interface coupling (3). The resulting time semi-discrete method reads as follows (see$^{27,29}$) for $n \ge 1$:

1. Explicit fluid viscous step: Find $\tilde{\mathbf{u}}^n : \Omega^f \to \mathbb{R}^d$ such that

$$\left\{ \begin{array}{ll} \frac{\rho^f}{\tau}(\tilde{\mathbf{u}}^n - \mathbf{u}^{n-1}) - \text{div}\sigma(\tilde{\mathbf{u}}^n, p^{n-1}) = \mathbf{0} & \text{in } \Omega^f, \\ \tilde{\mathbf{u}}^n = \mathbf{0} & \text{on } \Gamma, \\ \tilde{\mathbf{u}}^n = \dot{\mathbf{d}}^{n-1} & \text{on } \Sigma. \end{array} \right. \qquad (4)$$

2. Implicit pressure-displacement step: Find $\mathbf{u}^n : \Omega^f \to \mathbb{R}^d$, $p^n : \Omega^f \to \mathbb{R}$, $\mathbf{d}^n : \Sigma \to \mathbb{R}^d$ and $\dot{\mathbf{d}}^n : \Sigma \to \mathbb{R}^d$ such that

$$\left\{ \begin{array}{ll} \frac{\rho^f}{\tau}(\mathbf{u}^n - \tilde{\mathbf{u}}^n) + \nabla (p^n - p^{n-1}) = 0 & \text{in } \Omega^f, \\ \text{div}\mathbf{u}^n = 0 & \text{in } \Omega^f, \\ \mathbf{u}^n \cdot \mathbf{n} = 0 & \text{on } \Gamma, \end{array} \right. \qquad (5)$$

$$\left\{ \begin{array}{ll} \rho^s \epsilon \partial_{\tau} \dot{\mathbf{d}}^n + L \mathbf{d}^n = \mathbf{T}^n & \text{on } \Sigma, \\ \dot{\mathbf{d}}^n = \partial_{\tau} \mathbf{d}^n & \text{on } \Sigma, \\ \mathbf{d}^n = \mathbf{0} & \text{on } \partial\Sigma, \end{array} \right. \qquad (6)$$
$$ \left\{ \begin{array}{ll} \mathbf{u}_i^n \cdot \mathbf{n}_i = \dot{\mathbf{d}}^n \cdot \mathbf{n}_i & \text{on } \Sigma, \quad i = 1, 2, \\ T^n = -[\sigma(\tilde{\mathbf{u}}^n, p^n)\mathbf{n}] & \text{on } \Sigma. \end{array} \right. \qquad (7) $$

The viscous-step (4) is loosely coupled with the solid, which avoids strong coupling, whereas the step (5)-(7) guarantees added-mass free stability by the implicit treatment of the fluid pressure and solid inertia. From a computational point of view, the scheme can be reformulated exclusively in terms of $\tilde{\mathbf{u}}^n$, $p^n$, $\mathbf{d}^n$ and $\dot{\mathbf{d}}^n$ as shown in Algorithm 1, where the end-of-step velocity $\mathbf{u}^n$ has been eliminated by inserting (5)$_1$ into (4)$_1$ and by applying (5)$_2$ to (5)$_1$ in each sub-domain $\Omega_i$. Here, we have used the notation $p^{n, \star} \stackrel{\text{def}}{=} 2p^{n-1} - p^{n-2}$ for the second-order temporal extrapolation of the pressure.

**Algorithm 1** Time semi-discrete projection based semi-implicit scheme (from$^{27,29}$).

For $n \ge 2$:

1. Fluid viscous step: Find $\tilde{\mathbf{u}}^n : \Omega^f \to \mathbb{R}^d$ such that

$$ \left\{ \begin{array}{ll} \rho^f \partial_{\tau} \tilde{\mathbf{u}}^n - \mathbf{div} \sigma(\tilde{\mathbf{u}}^n, p^{n, \star}) = \mathbf{0} & \text{in } \Omega^f, \\ \tilde{\mathbf{u}}^n = \mathbf{0} & \text{on } \Gamma, \\ \tilde{\mathbf{u}}^n = \dot{\mathbf{d}}^{n-1} & \text{on } \Sigma. \end{array} \right. \qquad (8) $$

2. Pressure-displacement step: Find $p^n : \Omega^f \to \mathbb{R}$, $\mathbf{d}^n : \Sigma \to \mathbb{R}^d$ and $\dot{\mathbf{d}}^n : \Sigma \to \mathbb{R}^d$ such that

$$ \left\{ \begin{array}{ll} -\frac{\tau}{\rho^f} \Delta (p^n - p^{n-1}) = -\mathbf{div} \tilde{\mathbf{u}}^n & \text{in } \Omega^f, \\ \frac{\tau}{\rho^f} \nabla (p^n - p^{n-1}) \cdot \mathbf{n} = 0 & \text{on } \Gamma, \end{array} \right. \qquad (9) $$

$$ \left\{ \begin{array}{ll} \rho^s \varepsilon \partial_{\tau} \dot{\mathbf{d}}^n + L \mathbf{d}^n = \mathbf{T}^n & \text{on } \Sigma, \\ \dot{\mathbf{d}}^n = \partial_{\tau} \mathbf{d}^n & \text{on } \Sigma, \\ \mathbf{d}^n = \mathbf{0} & \text{on } \partial \Sigma, \end{array} \right. \qquad (10) $$

$$ \left\{ \begin{array}{ll} \frac{\tau}{\rho^f} \nabla (p_i^n - p_i^{n-1}) \cdot \mathbf{n}_i = (\tilde{\mathbf{u}}^n - \dot{\mathbf{d}}^n) \cdot \mathbf{n}_i & \text{on } \Sigma, \quad i = 1, 2, \\ \mathbf{T}^n = -[\sigma(\tilde{\mathbf{u}}^n, p^n)\mathbf{n}] & \text{on } \Sigma. \end{array} \right. \qquad (11) $$

An energy based stability result for the non-incremental version of (4)-(7) (i.e., with $p^{n-1} = 0$) with a fitted mesh based finite element approximation in space, has been reported in$^{27}$. Therein, stability is guaranteed under the CFL-like condition

$$ \rho^f h^2 + 2\mu\tau \lesssim \rho^s \varepsilon h, \qquad (12) $$

where $h > 0$ stands for the spatial grid parameter. Numerical evidence in$^{27}$ indicates that (12) is not necessary for stability. It is also worth noting that unconditional stability was achieved in$^{32}$ via a specific Nitsche's treatment of the viscous coupling. Unfortunately, the splitting error of the resulting scheme is known to be non uniform with respect to $h$, namely, to scales as $\mathcal{O}(\tau/h)$, so that suitable correction iterations are needed to enhance accuracy under restrictive constraints on the discretization parameters (as in$^{19}$).

**Remark 1.** Note that from the relation (11)$_1$ we get the continuity of flux on the pressure increment across $\Sigma$, namely,

$$ [\nabla (p^n - p^{n-1}) \cdot \mathbf{n}] = 0 \quad \text{on } \Sigma. $$

However, both the pressure $p^n$ and the pressure increment $p^n - p^{n-1}$ are generally discontinuous across $\Sigma$, so that the pressure-Poisson equation (9)$_1$ is not valid across $\Sigma$, only in $\Omega^f$. Nevertheless, most of the immersed boundary and fictitious domain methods involving fractional-step time-marching in the fluid assume that this relation is valid in $\Omega$ (see, e.g.,$^{4,28}$).
## 2.2 | Unfitted mesh approximation: fully discrete scheme

In the following, the closed spaces $H_Γ^1(\omega)$, of functions in $H^1(\omega)$ with zero trace on $Γ$, and $L_0^2(\omega)$, of functions in $L^2(\omega)$ with zero mean in $\omega$, will be considered. The scalar product in $L^2(\omega)$ is denoted by $(\cdot, \cdot)_\omega$. In this section we introduce a consistent unfitted mesh spatial approximation of the time semi-discrete scheme given by Algorithm 1. The fluid fields $(\tilde{u}^n, p^n)$ will hence be approximated in triangulations of $\Omega$ which are independent of $\Sigma$. To this purpose, it is important to note that:

* The velocity gradient $\nabla \tilde{u}^n$ and the pressure $p^n$ are discontinuous across $\Sigma$;

* The bulk relations $(8)_1$ and $(9)_1$ are not valid across $\Sigma$, only in $\Omega^f$ (see Remark 1).

In the case of fitted mesh approximations, these discontinuous features of the solution can be introduced in the discrete approximation in simple fashion (e.g., by considering cracked meshes with duplicated nodes on the interface). However, if the fluid mesh does not have a geometrical representation of the interface $\Sigma$, guaranteeing consistency of the approximations requires a specific treatment. In this paper, we build on the unifted Nitsche's based method for incompressible fluid-structure interaction with overlapping meshes reported in$^{18,19}$. The fundamental reasons for this choice are: (i) it is Lagrange multipliers free and robust; (ii) it is mathematically sound (i.e., optimal error estimates are obtained for spatial semi-discrete approximations of linear model problems such as (1)-(3) under reasonable regularity assumptions$^{18}$); and (iii) it naturally provides a consistent of the pressure-Poisson problem in step 2 of Algorithm 1.

In order to simplify the presentation, we assume that both $\Omega$ and $\Sigma$ are polyhedral. Let be $\{\mathcal{T}_h^s\}_{0<h<1}$ a family of triangulations of $\Sigma$. We then consider the standard space of continuous piecewise affine functions:

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

The weak form of the abstract solid elastic operator $L$ in (2) is assumed to be given by a positive and symmetric bi-linear form $a^s: W \times W \to \mathbb{R}$, where $W = [H^1(\Sigma)]_0^d$ denotes the space of admissible displacements. The discrete space for the solid displacement and velocity approximations is hence defined as $W_h = [X_h^s]^d \cap W$. For the fluid, we introduce two family of meshes $\{\mathcal{T}_{h,i}\}_{0<h<1}$, $i = 1, 2$, where each $\mathcal{T}_{h,i}$ covers the $i$-th fluid region $\Omega_i^f$ defined by $\Sigma$. Each mesh $\mathcal{T}_{h,i}$ is fitted to the exterior boundary $\Gamma$ but not to $\mathcal{T}_h^s$. Furthermore for every element $K \in \mathcal{T}_{h,1} \cap \mathcal{T}_{h,2}$ we assume that $K \cap \Sigma_h^n \neq \emptyset$. We denote by $\Omega_{h,i}$ the

<figure><img src="image_3.png" /><figcaption>Figure 2 One dimensional illustration of the construction of the discrete spaces $X_{h,i}$.</figcaption></figure>

domain covered by $\mathcal{T}_{h,i}$, viz.,

$$\Omega_{h,i} \stackrel{\text{def}}{=} \text{int} \left( \cup_{K \in \mathcal{T}_{h,i}} K \right).$$

We shall also make use of the following notation for the broken $L^2$-product in the whole computational domain

$$(\cdot, \cdot)_h \stackrel{\text{def}}{=} \sum_{i=1}^{2} \sum_{K \in \mathcal{T}_{h,i}} (\cdot, \cdot)_K.$$

For $i = 1, 2$, we can hence introduce the following spaces of continuous piecewise affine functions:

$$X_{h,i} \stackrel{\text{def}}{=} \left\{ v_h \in C^0(\overline{\Omega_{h,i}}) / v_{h|K} \in \mathbb{P}_1(K), \forall K \in \mathcal{T}_{h,i} \right\},$$
For the approximation of the fluid velocity and pressure we will consider the following discrete product spaces

$$V_h \stackrel{\text{def}}{=} ([X_{h,1}]^2 \times [X_{h,2}]^2) \cap [H_\Gamma^1(\Omega)]^d, \quad Q_h \stackrel{\text{def}}{=} (X_{h,1} \times X_{h,2}) \cap L^2(\Omega)_0, \quad (13)$$

which guarantee that interfacial (strong and weak) discontinuities are included in the discrete approximation of both the fluid velocity and pressure. Indeed, the functions of (13) are continuous in the physical fluid domain $\Omega^f$ but discontinuous across the interface $\Sigma$ (see Figure 2). Since the discrete pair $V_h/Q_h$ is not inf-sup stable, we consider a symmetric stabilization operator, such as, the one given by Continuous Interior Penalty method (see$^{34}$) over the whole computational domain:

$$s_h(p_h, q_h) = \frac{\gamma_p h^3}{\mu} \sum_{i=1}^{2} \sum_{F \in \mathcal{F}_{h,i}} (\llbracket \nabla p_h \rrbracket_F, \llbracket \nabla q_h \rrbracket_F)_F,$$

where $\mathcal{F}_{h,i}$ denotes the set of interior edges or faces of $\mathcal{T}_{h,i}$. Finally, we introduce the fluid discrete viscous bi-linear form

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

where the ghost-penalty operator is given by (see$^{35}$)

$$g_h(\mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} \gamma_g \mu h \sum_{i=1}^{2} \sum_{F \in \mathcal{F}_{i,h}^{\Sigma}} (\llbracket \nabla \mathbf{u}_{i,h} \rrbracket_F, \llbracket \nabla \mathbf{v}_{i,h} \rrbracket_F)_F \quad (14)$$

and where $\mathcal{F}_{i,h}^{\Sigma}$ denotes the set of interior edges or faces of the elements intersected by $\Sigma$. This operator guarantees robustness irrespectively to the way the interface is cutting the fluid mesh, by extending the coercivity of the viscous bi-linear form to the whole computational domain.

We can now introduce the unfitted mesh approximation of Algorithm 1 detailed in Algorithm 2. Its main ingredients are the following:

* Unfitted Nitsche's mortaring for the spatial discretization of the kinematic/dynamic viscous coupling in (8), (10) and (11), which is Lagrange multipliers free (i.e., without additional unknowns) and guarantees accuracy and robustness;

* For robustness, the Laplacian operator in the projection-step (9) is integrated over the whole computational domain, whereas for consistency the remaining fluid bulk terms in (8) and in (9) are integrated in the whole physical domain.

Note that the price to pay for consistency in the last point is a specific track of the interface intersections and the integration over cut elements (see, e.g.,$^{36,19}$ and the references therein). As regards the first point, it should be noted that in (16) the discrete interface stresses are exactly the variationally consistent viscous stress of (8). This constitutes a fundamental difference with respect to the Robin based semi-implicit and explicit coupling schemes respectively reported in$^{32,19}$ with fitted meshes. The main reason is to avoid the accuracy loss observed with this methods (see Section 2.1). The next section is devoted the energy based stability analysis of Algorithm 2.

*Remark 2.* In the case in which the interface has a boundary inside the fluid domain (the tip), we consider the construction of the fluid and solid spaces proposed in$^{19}$, which basically consists in introduce a virtual boundary $\tilde{\Sigma}$ which closes the fluid domain (see Figure 3).

<figure><img src="image_4.png" /><figcaption>**Figure 3** Case in which the $\Sigma$ has a boundary inside the fluid domain.</figcaption></figure>

Then, we enforce the kinematic/dynamic continuity of the fluid on $\tilde{\Sigma}$ in a discontinuous Galerkin fashion (see, e.g.,$^{37}$). More precisely, the following terms are added

$$- (\llbracket \sigma(\tilde{\mathbf{u}}_h^n, 0) \rrbracket_F, \llbracket \tilde{\mathbf{v}}_h \rrbracket_F)_{\tilde{\Sigma}} - (\llbracket \sigma(\tilde{\mathbf{v}}_h, 0) \rrbracket_F, \llbracket \tilde{\mathbf{u}}_h^n \rrbracket_F)_{\tilde{\Sigma}} + \frac{\gamma \mu}{h} (\llbracket \tilde{\mathbf{u}}_h^n \rrbracket_F, \llbracket \tilde{\mathbf{v}}_h \rrbracket_F)_{\tilde{\Sigma}}, \quad (\llbracket p_h^{n,\star} \rrbracket_F, \llbracket \tilde{\mathbf{v}}_h \rrbracket_F)_{\tilde{\Sigma}}, \quad (\llbracket p_h^{n,\star} \rrbracket_F, \llbracket \tilde{\mathbf{v}}_h \rrbracket_F)_{\tilde{\Sigma}}, \quad (17)$$
**Algorithm 2** Projection based semi-implicit scheme with unfitted meshes and static interface.

For $n \ge 1$:

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

$$ \rho^f (\partial_\tau \tilde{u}_h^n, \tilde{v}_h)_{\Omega^f} + a_h^f (\tilde{u}_h^n, \tilde{v}_h) - \sum_{i=1}^2 (\sigma(\tilde{u}_h^n, 0)n_i, \tilde{v}_{h,i})_{\Sigma} + \frac{\gamma\mu}{h} \sum_{i=1}^2 (\tilde{u}_{h,i}^n - \dot{d}_h^{n-1}, \tilde{v}_{h,i})_{\Sigma} \\ - \sum_{i=1}^2 (\tilde{u}_{h,i}^n - \dot{d}_h^{n-1}, \sigma(\tilde{v}_{h,i}, 0)n_i)_{\Sigma} = -(\nabla p_h^{n,\star}, \tilde{v}_h)_{\Omega^f} \quad (15) $$

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

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

$$ \left\{ \begin{array}{l} \frac{\tau}{\rho^f} (\nabla(p_h^n - p_h^{n-1}), \nabla q_h)_h + s_h(p_h^n, q_h) = \sum_{i=1}^2 (\tilde{u}_{h,i}^n - \dot{d}_{h,i}^n, q_{h,i} n_i)_{\Sigma} - (\text{div} \tilde{u}_h^n, q_h)_{\Omega^f}, \\ \\ \rho^s \epsilon (\partial_\tau \dot{d}_h^n, w_h)_\Sigma + a^s(d_h^n, w_h)_\Sigma = \frac{\gamma \mu}{h} \sum_{i=1}^2 ((\tilde{u}_{h,i}^n - \dot{d}_{h,i}^{n-1}), w_h)_\Sigma - \sum_{i=1}^2 (\sigma(\tilde{u}_{h,i}^n, p_h^n) n_i, w_h)_\Sigma \end{array} \right. \quad (16) $$

for all $(q_h, w_h) \in Q_h \times W_h$.

into the left- and right-hand side of step (15), respectively, where as in (16)$_1$ we add

$$ -\frac{\tau}{\rho^f} (\llbracket \nabla p_h^n \cdot \mathbf{n} \rrbracket, \llbracket q_h \rrbracket)_{\tilde{\Sigma}} - \frac{\tau}{\rho^f} (\llbracket p_h^n \rrbracket, \llbracket \nabla q_h \cdot \mathbf{n} \rrbracket)_{\tilde{\Sigma}} + \frac{\tau}{\rho^f} \frac{\gamma}{h} (\llbracket p_h^n \rrbracket, \llbracket q_h \rrbracket)_{\tilde{\Sigma}}, \quad (\llbracket \tilde{u}_h^n \rrbracket \cdot \mathbf{n}, \llbracket q_h \rrbracket)_{\tilde{\Sigma}} \quad (18) $$

to the left- and right-hand side, respectively.

## 2.3 | Energy based stability analysis

For the purpose of the analysis below, we recall the following estimate from$^{35}$:

$$ c_g (2\mu \|\epsilon(v_h)\|_{0,\Omega_h}^2 + g_h(v_h, v_h)) \le 2\mu \|\epsilon(v_h)\|_{0,\Omega^f}^2 + g_h(v_h, v_h) \quad (19) $$

for all $v_h \in V_h$ with $c_g > 0$ and the notation

$$ \| \cdot \|_{0,\Omega_h}^2 \stackrel{\text{def}}{=} (\cdot, \cdot)_h. $$

We shall also make use of the following (robust) discrete trace inequality

$$ h \|\epsilon(v_h)\mathbf{n}\|_{0,\Sigma}^2 \le C_T \|\epsilon(v_h)\|_{0,\Omega_h}^2 \quad (20) $$

for all $v_h \in V_h$. Let define the discrete total energy $E_h^n$ by the following expression:

$$ E_h^n \stackrel{\text{def}}{=} \frac{\rho_f}{2} \|\mathbf{u}_h^n\|_{0,\Omega^f}^2 + \frac{\rho_s \epsilon}{2} \|\dot{\mathbf{d}}_h^n\|_{\Sigma}^2 + \frac{1}{2} a^s(d_h^n, d_h^n) + \frac{\tau^2}{2 \rho_f} \|\nabla p_h^n\|_{0,\Omega_h}^2. $$

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

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

$$ c_g \gamma \ge 10C_{\text{TI}}, \quad \frac{3\gamma}{2} \mu \tau \le \rho_s \epsilon h, \quad (21) $$

the discrete energy estimate presented below holds:

$$ E_h^n \lesssim E_h^0 \quad (22) $$

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

*Proof.* We first introduce the $L^2$-projection operator $\pi_h : [L^2(\Omega)]^d \to V_h$ given by

$$ (\pi_h s, v_h)_{\Omega^f} = (s, v_h)_{\Omega^f} \quad (23) $$
for all $v_h \in V_h^n$. Note that, depending on how the solid mesh $T_h^s$ intersects the fluid overlapping meshes $T_{h,i}$, $\pi_h s$ may not be uniquely defined in the whole computation domain. However, a simple argument show that $\pi_h s$ is uniquely defined in the physical domain $\Omega^f$ (it suffices to remove the indetermination by blocking appropriate nodes outside the physical domain). This feature will be enough for the purpose of the present proof. In a similar fashion, we introduce the intermediate velocity $u_h^n \in V_h^n$ given by

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

for all $v_h \in V_h^n$, so that $u_h^n$ is uniquely defined in $\Omega^f$. In particular, owing to (23) and (24), we have

$$ u_h^n = \tilde{u}_h^n - \frac{\tau}{\rho_f} \pi_h \nabla (p_h^n - p_h^{n-1}) \quad \text{in } \Omega^f. \quad (25) $$

Now, since the fluid bulk terms of the viscous step (15) are integrated (only) in the physical domain and using (24), it can alternatively written as

$$ \frac{\rho_f}{\tau} (\tilde{u}_h^n, \tilde{v}_h)_{\Omega} + a_h^{f,n} (\tilde{u}_h^n, \tilde{v}_h)_{\Sigma} + \sum_{i=1}^{2} (\sigma(\tilde{u}_h^n, 0)n_i, \tilde{v}_{h,i})_{\Sigma} + \frac{\gamma \mu}{h} \sum_{i=1}^{2} (\tilde{u}_{h,i}^n - \dot{\tilde{a}}_h^{n-1}, \tilde{v}_{h,i})_{\Sigma} \\ - \sum_{i=1}^{2} (\tilde{u}_{h,i}^n - \dot{\tilde{a}}_h^{n-1}, \sigma(\tilde{v}_{h,i}, 0)n_i)_{\Sigma} = \frac{\rho_f}{\tau} (u_h^{n-1}, \tilde{v}_h)_{\Omega} - (\nabla p_h^{n-1}, \tilde{v}_h)_{\Omega} \quad (26) $$

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

We then proceed by testing the relation (24) with $v_h^n = u_h^n$, with yields

$$ \frac{\rho_f}{2\tau} [ \|u_h^n\|_{0,\Omega^f}^2 - \|\tilde{u}_h^n\|_{0,\Omega^f}^2 + \|u_h^n - \tilde{u}_h^n\|_{0,\Omega^f}^2 ] + (\nabla(p_h^n - p_h^{n-1}), u_h^n)_{\Omega^f} = 0. \quad (27) $$

By inserting (25) into the last equality and by rearranging the terms we get

$$ \frac{\rho_f}{2\tau} [ \|u_h^n\|_{0,\Omega^f}^2 - \|\tilde{u}_h^n\|_{0,\Omega^f}^2 + (\nabla(p_h^n - p_h^{n-1}), \tilde{u}_h^n)_{\Omega^f} - \frac{\tau}{2\rho_f} \|\pi_h \nabla(p_h^n - p_h^{n-1})\|_{0,\Omega^f}^2 = 0. \quad (28) $$

On the other hand, by testing (26) with $\tilde{v}_h^n = \tilde{u}_h^n$ and using (19) we have

$$ \frac{\rho_f}{2\tau} [ \|\tilde{u}_h^n\|_{0,\Omega^f}^2 - \|u_h^{n-1}\|_{0,\Omega^f}^2 + \|\tilde{u}_h^n - u_h^{n-1}\|_{0,\Omega^f}^2 ] + 2c_g\mu \|\epsilon(\tilde{u}_h^n)\|_{0,\Omega_h}^2 + (\nabla p_h^{n-1}, \tilde{u}_h^n)_{\Omega^f} \\ + \frac{\gamma\mu}{h} \sum_{i=1}^{2} (\tilde{u}_{h,i}^n - \dot{\tilde{a}}_h^{n-1}, \tilde{u}_{h,i}^n)_{\Sigma} - 2\mu \sum_{i=1}^{2} (\epsilon(\tilde{u}_{h,i}^n)n_i, \tilde{u}_{h,i}^n)_{\Sigma} - 2\mu \sum_{i=1}^{2} (\epsilon(\tilde{u}_{h,i}^n)n_i, \tilde{u}_{h,i}^n - \dot{\tilde{a}}_h^{n-1})_{\Sigma} \le 0 \quad (29) $$

and by testing the fluid projection-step (16)$_1$ with $q_h = p_h^n$ and by integrating by parts the divergence term we get

$$ \frac{\tau}{2 \rho_f} [ \| \nabla p_h^n \|_{0,\Omega_h}^2 - \| \nabla p_h^{n-1} \|_{0,\Omega_h}^2 + \| \nabla (p_h^n - p_h^{n-1}) \|_{0,\Omega_h}^2 ] - (\nabla p_h^n, \tilde{u}_h^n)_{\Omega^f} + \sum_{i=1}^{2} (\dot{\tilde{a}}_h^n, n_i p_{h,i}^n)_{\Sigma} + s_h(p_h^n, p_h^n) = 0. \quad (30) $$

Finally, by adding the relations (28)-(30) we get the following energy estimate for the fluid

$$ \frac{\rho_f}{2} \partial_\tau \|u_h^n\|_{0,\Omega^f}^2 + 2c_g \mu \|\epsilon(\tilde{u}_h^n)\|_{0,\Omega_h}^2 + \frac{\tau}{2\rho_f} [ \|\nabla p_h^n\|_{0,\Omega_h}^2 - \|\nabla p_h^{n-1}\|_{0,\Omega_h}^2 ] \\ + \frac{\tau}{2\rho_f} \|(\mathcal{I} - \pi_h)\nabla(p_h^n - p_h^{n-1})\|_{\Omega^f}^2 + \frac{\tau}{2\rho_f} \|\nabla(p_h^n - p_h^{n-1})\|_{\Omega_h \setminus \Omega^f}^2 + \sum_{i=1}^2 (\dot{\tilde{a}}_h^n, n_i p_{h,i}^n)_{\Sigma} \\ + \frac{\gamma\mu}{h} \sum_{i=1}^2 (\tilde{u}_{h,i}^n - \dot{\tilde{a}}_h^{n-1}, \tilde{u}_{h,i}^n)_{\Sigma} - 2\mu \sum_{i=1}^2 (\epsilon(\tilde{u}_{h,i}^n)n_i, \tilde{u}_{h,i}^n)_{\Sigma} - 2\mu \sum_{i=1}^2 (\epsilon(\tilde{u}_{h,i}^n)n_i, \tilde{u}_{h,i}^n - \dot{\tilde{a}}_h^{n-1})_{\Sigma} \le 0. \quad (31) $$

We now proceed by testing the solid equation (16)$_2$ with $w_h = \dot{a}_h^n$, which yields

$$ \frac{\rho_s \epsilon}{2 \tau} [ \| \dot{a}_h^n \|_{0,\Sigma}^2 - \| \dot{a}_h^{n-1} \|_{0,\Sigma}^2 + \| \dot{a}_h^n - \dot{a}_h^{n-1} \|_{0,\Sigma}^2 ] + \frac{1}{2\tau} [ a_s(\dot{d}_h^n, \dot{d}_h^n) - a_s(\dot{d}_h^{n-1}, \dot{d}_h^{n-1}) + a_s(\dot{d}_h^n - \dot{d}_h^{n-1}, \dot{d}_h^n - \dot{d}_h^{n-1}) ] \\ + 2\mu \sum_{i=1}^2 (\epsilon(\tilde{u}_{h,i}^n)n_i, \dot{a}_h^n)_{\Sigma} - \sum_{i=1}^2 (p_{h,i}^n n_i, \dot{a}_h^n)_{\Sigma} - \frac{\gamma\mu}{h} \sum_{i=1}^2 (\tilde{u}_{h,i}^n - \dot{\tilde{a}}_h^{n-1}, \dot{a}_h^n)_{\Sigma} = 0. $$
By adding this relation to (31) we get the following total energy estimate

$$
\begin{equation}
\begin{aligned}
& \frac{\rho_f}{2} \partial_\tau \|\mathbf{u}_h^n\|_{0,\Omega^t}^2 + 2c_g \mu \|\epsilon(\tilde{\mathbf{u}}_h^n)\|_{0,\Omega_h}^2 + \frac{\tau}{2\rho_f} \left[ \|\nabla p_h^n\|_{0,\Omega_h}^2 - \|\nabla p_h^{n-1}\|_{0,\Omega_h}^2 \right] \\
& \qquad + \frac{\rho_s \epsilon}{2\tau} \left[ \|\dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 - \|\dot{\mathbf{d}}_h^{n-1}\|_{0,\Sigma}^2 + \|\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}\|_{0,\Sigma}^2 \right] \\
& \qquad + \frac{1}{2\tau} \left[ a_s(\mathbf{d}_h^n, \mathbf{d}_h^n) - a_s(\mathbf{d}_h^{n-1}, \mathbf{d}_h^{n-1}) \right] \\
& - 2\mu \sum_{i=1}^2 \left( \epsilon(\tilde{\mathbf{u}}_{h,i}^n) \mathbf{n}_i, \tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^n \right)_\Sigma - 2\mu \sum_{i=1}^2 \left( \epsilon(\tilde{\mathbf{u}}_{h,i}^n) \mathbf{n}_i, \tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^{n-1} \right)_\Sigma \\
& \qquad + \underbrace{\frac{\gamma \mu}{h} \sum_{i=1}^2 (\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^{n-1}, \tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^n)_\Sigma}_{T_1} \le 0. \quad (32)
\end{aligned}
\end{equation}
$$

Terms $T_1$ can be bounded from every side of the interface by adding and subtracting $\dot{\mathbf{d}}_h^n$, using the Cauchy–Schwarz, Young's and trace inequalities (20), as follows:

$$
\begin{align}
T_1 &= -2\mu (\epsilon(\tilde{\mathbf{u}}_{h,i}^n) \mathbf{n}_i, \dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1})_\Sigma - 4\mu (\epsilon(\tilde{\mathbf{u}}_{h,i}^n) \mathbf{n}_i, \tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^n)_\Sigma \nonumber \\
    &\geq -\frac{20 \mu C_{TI}}{\gamma} \|\epsilon(\tilde{\mathbf{u}}_h^n)\|_{0,\Omega_h}^2 - \frac{\gamma \mu}{4h} \|\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 - \frac{\gamma \mu}{4h} \|\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}\|_{0,\Sigma}^2. \tag{33}
\end{align}
$$

Similarly, for the second term we have

$$
\begin{align}
T_2 &= \frac{\gamma \mu}{h} \|\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 + \frac{\gamma \mu}{h} (\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}, \tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^n)_\Sigma \nonumber \\
    &\geq \frac{\gamma \mu}{2h} \|\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 - \frac{\gamma \mu}{2h} \|\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}\|_{0,\Sigma}^2. \tag{34}
\end{align}
$$

By applying (33)-(34) we get

$$
\begin{equation}
\begin{aligned}
& \frac{\rho_f}{2} \partial_\tau \|\mathbf{u}_h^n\|_{0,\Omega^t}^2 + 2\mu \left(c_g - \frac{10 C_{TI}}{\gamma}\right) \|\epsilon(\tilde{\mathbf{u}}_h^n)\|_{0,\Omega_h}^2 + \frac{\tau}{2\rho_f} \left[ \|\nabla p_h^n\|_{0,\Omega_h}^2 - \|\nabla p_h^{n-1}\|_{0,\Omega_h}^2 \right] \\
& \qquad + \frac{\rho_s \epsilon}{2} \partial_\tau \|\dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 + \left(\frac{\rho_s \epsilon}{2\tau} - \frac{3\gamma \mu}{4h}\right) \|\dot{\mathbf{d}}_h^n - \dot{\mathbf{d}}_h^{n-1}\|_{0,\Sigma}^2 + \frac{1}{2\tau} [a_s(\mathbf{d}_h^n, \mathbf{d}_h^n) - a_s(\mathbf{d}_h^{n-1}, \mathbf{d}_h^{n-1})] \\
& \qquad + \frac{\gamma \mu}{4h} \|\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_h^n\|_{0,\Sigma}^2 \le 0. \quad (35)
\end{aligned}
\end{equation}
$$

Finally, the energy estimate (22) follows from (35) under de assumption (21), which completes the proof. □

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

It is worth noting that the stability condition (21) provided by Theorem 1 is similar to the stability condition (12) obtained in²⁷ for Algorithm 1 with fitted meshes.

## 3 | NON-LINEAR MODEL: MOVING INTERFACE

In this section, we propose an extension of the semi-implicit coupling scheme given by Algorithm 2 to the case of non-linear fluid–structure interaction problems involving an incompressible viscous fluid and a moving immersed thin-walled structure. The fluid is described by the Navier-Stokes equations (in Eulerian form) and the structure by a possibly non-linear (beam or shell) solid model (in Lagrangian form).

### 3.1 | Problem setting

Let $\Omega \subset \mathbb{R}^d$ be the reference configuration of the fluid domain with boundary $\Gamma \stackrel{\text{def}}{=} \partial\Omega$, and $\Sigma \subset \mathbb{R}^2$ be the reference solid mid-surface. In contrast to Section 2, the structure is now supposed to move within the fluid domain. The current position of the interface $\Sigma(t)$ is described in terms of a deformation map $\phi : \Sigma \times \mathbb{R}^+ \to \mathbb{R}^d$ as $\Sigma(t) = \phi(\Sigma, t)$, with $\phi \stackrel{\text{def}}{=} I_\Sigma + \mathbf{d}$ and where $\mathbf{d}$ denotes the solid displacement. To simplify the notation we will refer to $\phi_t \stackrel{\text{def}}{=} \phi(\cdot, t)$, so that we also have $\Sigma(t) = \phi_t(\Sigma)$. Note that the fluid control volume is now time-dependent, namely $\Omega^f(t) \stackrel{\text{def}}{=} \Omega \setminus \Sigma(t) \subset \mathbb{R}^d$ with boundary $\partial\Omega(t) = \Sigma(t) \cup \Gamma$. The notations introduced in Section 2 for the surface normal vector $\mathbf{n}_\Sigma$, jumps and average operators remain valid with the sole difference that they refer to the current interface position $\Sigma(t)$. The considered coupled problem reads therefore as follow: find the fluid velocity and pressure $\mathbf{u} : \Omega \times \mathbb{R}^+ \to \mathbb{R}^d$, $\mathbf{p} : \Omega \times \mathbb{R}^+ \to \mathbb{R}$, the solid displacement and velocity $\mathbf{d} : \Sigma \times \mathbb{R}^+ \to \mathbb{R}^d$, $\dot{\mathbf{d}} : \Sigma \times \mathbb{R}^+ \to \mathbb{R}^d$ such that for all $t \in \mathbb{R}^+$ we have
$$ egin{cases} \rho^f (\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u}) - \mathbf{div} \sigma(\mathbf{u}, p) = \mathbf{0} & \text{in } \Omega^f(t), \\ \mathbf{div} \mathbf{u} = \mathbf{0} & \text{in } \Omega^f(t), \\ \mathbf{u} = \mathbf{0} & \text{on } \Gamma, \end{cases} \quad (36) $$

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

$$ \left\{ \begin{array}{ll} \boldsymbol{\phi} = \mathbf{I}_{\Sigma} + \mathbf{d}, & \Sigma(t) = \boldsymbol{\phi}_t(\Sigma), \quad \Omega^f(t) = \Omega \setminus \Sigma(t), \\ \mathbf{u} = \dot{\mathbf{d}} \circ \boldsymbol{\phi}_t^{-1} & \text{on } \Sigma(t), \\ \displaystyle\int_{\Sigma} T \cdot \mathbf{w} = - \int_{\Sigma(t)} [\boldsymbol{\sigma}(\mathbf{u}, p)\mathbf{n}] \cdot \mathbf{w} \circ \boldsymbol{\phi}_t^{-1}. & \end{array} \right. \quad (38) $$

The relations in (38) respectively enforce the geometrical compatibility, the kinematic and the dynamic coupling at the interface between the fluid and the solid. In the next section, we propose a numerical method for the coupled system (36)-(38) based on Algorithm 2.

## 3.2 | Numerical methods

With the purpose of avoiding geometrical non-linearities in the fluid, we will discretize the geometric compatibility condition (38)$_1$, namely $\Omega^f(t) = \Omega \setminus \Sigma(t)$, in an explicit fashion. For a given displacement approximation $\mathbf{d}_h^n \in W_h$, we define by $\boldsymbol{\phi}_h^n$ its associated deformation map as $\boldsymbol{\phi}_h^n \stackrel{\text{def}}{=} \mathbf{I}_\Sigma + \mathbf{d}_h^n$. This map characterizes the interface position, at time level $n$, as $\Sigma^n \stackrel{\text{def}}{=} \boldsymbol{\phi}_h^n(\Sigma)$. We hence propose to explicitly update the physical fluid domain as

$$ \Omega^{f,n} \stackrel{\text{def}}{=} \Omega \setminus \Sigma^{n-1}. \quad (39) $$

For the fluid discrete spaces, $V_h^n$ and $Q_h^n$, we proceed as in Section 2.2, with the difference that they are now given in term of the time-dependent overlapping meshes $\mathcal{T}_{h,i}^n$ covering each side of $\Omega^{f,n}$. We recall that this functional spaces are made of functions that are continuous in each side of $\Omega^{f,n}$ but discontinuous across $\Sigma^{n-1}$. Finally, we introduce the broken $L^2$-product in the moving computational domain as

$$ (\cdot, \cdot)_{n,h} \stackrel{\text{def}}{=} \sum_{i=1}^{2} \sum_{K \in \mathcal{T}_{h,i}^n} (\cdot, \cdot)_K. $$

The approximation space for the solid $W_h$ is the same as in Section 2.2.

For the spatial approximation of the fluid, we introduce the following discrete tri-linear form associated to the convective term

$$ c^n(\mathbf{z}_h, \mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} \rho^f(\mathbf{z}_h \cdot \nabla \mathbf{u}_h, \mathbf{v}_h)_{\Omega^{f,n}} + \frac{\rho^f}{2} \left( (\mathbf{div}\mathbf{z}_h)\mathbf{u}_h, \mathbf{v}_h \right)_{\Omega^{f,n}} - \rho^f(\left[\mathbf{z}_h\right] \cdot \mathbf{n})[\mathbf{u}_h], \left[\mathbf{v}_h\right]_{\Sigma^{n-1}} - \frac{\rho^f}{2} (\left[\mathbf{z}_h \cdot \mathbf{n}\right], \left[\mathbf{u}_h \cdot \mathbf{v}_h\right]_{\Sigma^{n-1}}), \quad (40) $$

where the three last terms are added in order to guarantee that $c^n(\mathbf{v}_h, \mathbf{z}_h, \mathbf{z}_h) = 0$ for all $\mathbf{z}_h \in V_h^n$ (see $^{19}$). Numerical instabilities due to the lack of inf-sup compatibility of the discrete spaces and to large local Reynolds number, will be handled by the continuous interior penalty stabilization method (CIP) of$^{38,39}$. The associated symmetric velocity and pressure stabilization operators are given by:

$$ s_{v,h}^n(\mathbf{z}_h; \mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} \gamma_v h^2 \sum_{i=1}^2 \sum_{F \in \mathcal{F}_{h,i}^n} \xi(\text{Re}_F(\mathbf{z}_h)) \|\mathbf{z}_h \cdot \mathbf{n}\|_{L^\infty(F)} (\llbracket \nabla \mathbf{u}_h \rrbracket_F, \llbracket \nabla \mathbf{v}_h \rrbracket_F)_F, $$

$$ s_{p,h}^n(\mathbf{z}_h; p_h, q_h) \stackrel{\text{def}}{=} \gamma_p h^2 \sum_{i=1}^2 \sum_{F \in \mathcal{F}_{h,i}^n} \frac{\xi(\text{Re}_F(\mathbf{z}_h))}{\|\mathbf{z}_h\|_{L^\infty(F)}} (\llbracket \nabla p_h \rrbracket_F, \llbracket \nabla q_h \rrbracket_F)_F, $$

where $\mathcal{F}_{h,i}^n$ denotes the set of interior edges or faces of $\mathcal{T}_{h,i}^n$, $\text{Re}_F(\mathbf{z}_h) \stackrel{\text{def}}{=} \rho^f \|\mathbf{z}_h\|_{L^\infty(F)} h \mu^{-1}$ denotes the local Reynolds number, $\xi(x) \stackrel{\text{def}}{=} \min\{1, x\}$ is a cut-off function and $\gamma_p, \gamma_v > 0$ are user-defined parameters. In order to guarantee robustness with respect
to the way the interface $\Sigma^n$ is cutting the fluid domain $\Omega$, we introduce the time-dependent ghost-penalty operator, given by

$$g_h^n(\mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} \gamma_g \mu h \sum_{i=1}^{2} \sum_{F \in \mathcal{F}_{i,h}^{\Sigma^{n-1}}} (\llbracket \nabla \mathbf{u}_{i,h} \rrbracket_F, \llbracket \nabla \mathbf{v}_{i,h} \rrbracket_F)_F,$$

where $\mathcal{F}_{i,h}^{\Sigma^{n-1}}$ denotes the set of interior edges or faces of the elements intersected by $\Sigma^{n-1}$. Finally, we collect all the above terms in a single contribution

$$a_h^{f,n}(\mathbf{z}_h; \mathbf{u}_h, \mathbf{v}_h) \stackrel{\text{def}}{=} c^n(\mathbf{z}_h, \mathbf{u}_h, \mathbf{v}_h) + 2\mu(\epsilon(\mathbf{u}_h), \epsilon(\mathbf{v}_h))_{\Omega^{f,n}} + s_{\tilde{\nu}_h}^n(\mathbf{z}_h; \mathbf{u}_h, \mathbf{v}_h) + g_h^n(\mathbf{u}_h, \mathbf{v}_h). \quad (41)$$

With all the above ingredients, we propose to approximate (36)-(38) by the semi-implicit coupling scheme reported in Algorithm 3. The basic idea consists in combining the interface kinematic/dynamic coupling of Algorithm 2 with the explicit treatment of the geometrical compatibility (39).

**Algorithm 3** Projection based semi-implicit scheme with unfitted meshes and moving interfaces.

For $n \ge 1$:

1. Interface update:

$$\Sigma^{n-1} = \phi_h^{n-1}(\Sigma), \quad \Omega^{f,n} = \Omega \setminus \Sigma^{n-1}.$$

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

$$\rho^f (\partial_{\tau} \tilde{\mathbf{u}}_h^n, \tilde{\mathbf{v}}_h)_{\Omega^{f,n}} + a_h^{f,n} (\tilde{\mathbf{u}}_h^{n-1}; \tilde{\mathbf{u}}_h^n, \tilde{\mathbf{v}}_h) - \sum_{i=1}^{2} (\sigma(\tilde{\mathbf{u}}_{h,i}^n, 0)n_i, \tilde{\mathbf{v}}_{h,i})_{\Sigma^{n-1}} + \frac{\gamma \mu}{h} \sum_{i=1}^{2} (\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_{h}^{n-1}, \tilde{\mathbf{v}}_{h,i})_{\Sigma^{n-1}} \\ - \sum_{i=1}^{2} (\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_{h}^{n-1}, \sigma(\tilde{\mathbf{v}}_{h,i}, 0)n_i)_{\Sigma^{n-1}} = -(\nabla p_h^{n,\star}, \tilde{\mathbf{v}}_h)_{\Omega^{f,n}} \quad (42)$$

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

3. Pressure-displacement step: Find $(p_h^n, d_h^n) \in \mathcal{Q}_h^n \times \mathbf{W}_h$ with $\dot{d}_h^n = \partial_{\tau} d_h^n$, such that

$$\left\{ \begin{array}{l} \frac{\tau}{\rho^f} (\nabla(p_h^n - p_h^{n-1}), \nabla q_h)_{n,h} + s_{p,h}^n (\tilde{\mathbf{u}}_h^n; p_h^n, q_h) = \sum_{i=1}^{2} (\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_{h,i}^n, q_{h,i} n_i)_{\Sigma^{n-1}} - (\text{div} \tilde{\mathbf{u}}_h^n, q_h)_{\Omega^{f,n}}, \\ \\ \rho^s \epsilon (\partial_{\tau} \dot{\mathbf{d}}_h^n, \mathbf{w}_h)_\Sigma + a^s (d_h^n, \mathbf{w}_h) = \frac{\gamma \mu}{h} \sum_{i=1}^{2} ((\tilde{\mathbf{u}}_{h,i}^n - \dot{\mathbf{d}}_{h,i}^{n-1}), \mathbf{w}_h)_{\Sigma^{n-1}} - \sum_{i=1}^{2} (\sigma(\tilde{\mathbf{u}}_{h,i}^n, p_{h,i}^n)n_i, \mathbf{w}_h)_{\Sigma^{n-1}} \end{array} \right. \quad (43)$$

for all $(q_h, \mathbf{w}_h) \in \mathcal{Q}_h^n \times \mathbf{W}_h$. Then, set $\phi_h^n = I_\Sigma + d_h^n$.

Note that steps (42) and (43) of Algorithm 3 involves integrals of functions associated with different time levels, namely,

$$(\tilde{\mathbf{u}}_h^{n-1}, \tilde{\mathbf{v}}_h)_{\Omega^n}, \quad (\nabla p_h^{n-1}, \tilde{\mathbf{v}}_h)_{\Omega^n}, \quad (\nabla p_h^{n-2}, \tilde{\mathbf{v}}_h)_{\Omega^n}, \quad (\nabla p_h^{n-1}, \nabla q_h)_{n,h}.$$

with $\tilde{\mathbf{u}}_h^n \in \mathbf{V}_h^n$, $p_h^{n-1} \in \mathcal{Q}_h^{n-1}$ and $p_h^{n-2} \in \mathcal{Q}_h^{n-2}$. This requires the integration of products of functions that might be discontinuous at different locations in the same element. In order avoid the simultaneous intersection of different interface locations with the same fluid element, we consider the approach introduced in ¹⁹ (see also ⁴⁰), which basically consists in locally shifting the discontinuity at time $t^*$ to the structure location at time $t^n$, where $t^*$ refers to $t^{n-1}$ and $t^{n-2}$ respectively. In the case where the discontinuities are located in different elements, the quadrature is performed in a standard fashion since we keep track of the (previous) intersections at different times and we can treat each discontinuity separately.

**Remark 3.** In the case of partially intersected fluid domain with dynamic interface, we proceed as in Remark 2. The terms in (17) and (18) are now evaluated on $\tilde{\Sigma}^{n-1}$ and we add the following convective Discontinuous Galerkin contributions (see, e.g., ³⁷) to the tri-linear form (40):

$$-\rho^f (\llbracket \mathbf{z}_h \rrbracket \cdot \mathbf{n} [\llbracket \mathbf{u}_h \rrbracket], \llbracket \mathbf{v}_h \rrbracket)_{\tilde{\Sigma}^{n-1}} - \frac{\rho^f}{2} (\llbracket \mathbf{z}_h \cdot \mathbf{n} \rrbracket, \llbracket \mathbf{u}_h \cdot \mathbf{v}_h \rrbracket)_{\tilde{\Sigma}^{n-1}}.$$
Remark 4. Note that whenever $\dot{\mathbf{d}}_{h,i}^n, \dot{\mathbf{d}}_{h,i}^{n-1}$ and $\mathbf{w}_h$ are integrated on $\Sigma^{n-1}$, one has to understand $\dot{\mathbf{d}}_{h,i}^n \circ (\phi_h^{n-1})^{-1}, \dot{\mathbf{d}}_{h,i}^{n-1} \circ (\phi_h^{n-1})^{-1}$ and $\mathbf{w}_h \circ (\phi_h^{n-1})^{-1}$ respectively. This abuse of notation is simply made to ease the presentation.

# 4 | NUMERICAL EXPERIMENTS

In this section we illustrate the stability and accuracy of the proposed semi-implicit scheme (Algorithm 3) in different 2D test cases motivated by the simulation of heart valves and of micro-encapsulation. To this purpose we compare the results obtained with Algorithm 3 and those obtained with the strongly coupled and loosely coupled (stabilized explicit coupling) schemes proposed in¹⁹. The implicit step in Algorithm 3 is solved in a partitioned fashion by a Dirichlet-Neumann based Newton-GMRES iterative algorithm. In all the tests, the solid is described by a non-linear Reissner-Mindlin curved beam model with a MITC spatial discretization. All the units are given in the CGS units system

## 4.1 | Idealized valve without contact

The first example is the heart-valve-inspired benchmark proposed in⁴¹,⁴²,⁴³,⁴⁴,²¹. The considered geometry is shown in Figure 4(a). The fluid domain is a rectangle $\Omega = [0, 8] \times [0, 0.805]$, while the immersed solid reference configuration $\Sigma$ is the straight segment AB, with $A = (2, 0)$ and $B = (2, 0.7)$.

<figure><img src="image_6.png" /><figcaption>Figure 4 (a) Geometric configuration of the idealized valve without contact, (b) Zoom of the fluid and solid meshes.</figcaption></figure>

The physical parameters used for the fluid in this test are $\rho_f = 100$, $\mu = 10$. While for the solid we have $\rho_s = 100$, $\epsilon_s = 0.0212$, the Young's modulus $E = 5.6 \cdot 10^7$ and Poisson's ratio $\nu = 0.4$. Concerning the boundary conditions, no-slip boundary condition is applied on $\Gamma_w$, a symmetry condition is imposed on $\Gamma_{sym}$, zero traction on $\Gamma_{out}$ and finally on $\Gamma_{in}$ the following half parabolic profile is applied:

$$u_{\max}(t) = 5(0.805)^2(\sin(2\pi t) + 1.1), \quad t \in \mathbb{R}^+.$$

The solid rotation and displacement are set to zero at the bottom endpoint A and rest initial conditions are considered for both fluid and solid.

The solid mesh is made of 64 edges while the fluid unfitted mesh is made of 18662 triangles (see Figure 4(b)). The Nitsche parameter is set to $\gamma = 10$ and the Ghost penalty parameter to $\gamma_g = 1$. The CIP stabilization parameters are $\gamma_v = \gamma_p = 10^{-2}$. Three different levels of time-step refinement, $\tau \in \{(10^{-3}/2^i)\}_{i=0}^2$, are considered in order to compare results from Algorithms 3 and the loosely coupled and strongly coupled schemes. The final time is $T = 3$, which corresponds to 3 full oscillations cycles for the structure.

For illustration purposes, snapshots of the fluid velocity magnitude and the position of the interface, computed with $\tau = 10^{-3}$, are shown in Figure 5 at time $t = 0.5$ and $t = 0.8$ respectively. A very good agreement is obtained for the three methods already with the larger time step and all algorithms reproduce very well the vortex induced after the leaflet. The two times selected correspond to a situation of opening of the valve at $t = 0.5$ and closing at $t = 0.8$. At opening state, there is an increasing velocity magnitude on top of the channel, while the velocity is decreasing at closing state. Figure 6 presents the pressure elevation computed with the coarsest time step, $\tau = 10^{-3}$, obtained with the semi-implicit coupling scheme (Algorithm 3), the loosely coupled case.
<figure><img src="image_7.png" /><figcaption>Figure 5 Velocity magnitude snapshots at $\tau = 10^{-3}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6 Pressure snapshots at $\tau = 10^{-3}$.</figcaption></figure>

coupled and the strongly coupled schemes at the same time instants as before. The discontinuity of the pressure is well captured with all methods. A very good agreement can be observed between Algorithm 3 and the strongly coupled scheme, while some differences are clear visible in the loosely coupled scheme.

Figures 7 and 8 report the displacement history of the upper structure endpoint $B$ as function of time, in terms of $x$-displacement and $y$-displacement respectively. Algorithm 3 delivers practically the same results as the strongly coupled scheme (the two curves are indistinguishable already with the larger time step), whereas some differences are clearly visible with the loosely coupled scheme. This mismatch is reduced with the time-step refinement.
<figure><img src="image_9.png" /><figcaption>Figure 7 Time evolution of the x-displacement for the structure endpoint B.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8 Time evolution of the y-displacement for the structure endpoint B.</figcaption></figure>

## 4.2 | Idealized valve with contact

The second numerical example corresponds to the idealized valve test with contact introduced in$^{45}$. It is an extension of the previous one in which the structure is sufficiently long to get in contact with $\Gamma_{sym}$. The geometry is shown in Figure 9(a). The fluid domain $\Omega$ is the same as in the previous example, while as reference configuration for the solid $\Sigma$ we consider a curve segment of extreme $A = (4, 0)$ and $B = (5.112, 0.483)$, parametrized by the analytical function

$$y(x) = \frac{1}{2} \sqrt{1 - \frac{(x - 11/2)^2}{(3/2)^2}}, \quad x \in [4, 5.112].$$

The physical parameters used for the fluid in this test are $\rho_f = 100$, $\mu = 10$. While for the solid we have $\rho_s = 100$, $\epsilon_s = 0.0212$, the Young's modulus $E = 5.6 \cdot 10^7$ and Poisson's ratio $\nu = 0.4$.

<figure><img src="image_11.png" /><figcaption>Figure 9 (a) Geometric configuration of the idealized valve with contact, (b) Zoom of the leaflet mesh and fluid mesh.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10 Velocity magnitude snapshots with $\tau = 10^{-3}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11 Pressure snapshots with $\tau = 10^{-3}$.</figcaption></figure>

Regarding the boundary condition, a symmetry condition is enforced on $\Gamma_{sym}$, a no-slip condition on $\Gamma_{w}$, zero traction on the outflow boundary $\Gamma_{out}$ and a traction condition is imposed on $\Gamma_{in}$ in terms of the following time-dependent pressure:

$$p_{\text{in}}(t) = \begin{cases} -200 \text{ atanh}(100t) & \text{if } 0 < t < 0.7, \\ 200 & \text{if } t \ge 0.7. \end{cases}$$

The final time is $T = 1$ and it corresponds to one full valve oscillation cycle. The fluid and the solid are initially at rest and the beam is pinched at the bottom tip A.
<figure><img src="image_14.png" /><figcaption>Figure 12 Time evolution of the x-displacement for the structure endpoint B.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13 Time evolution of the y-displacement for the structure endpoint B.</figcaption></figure>

In order to avoid penetration on $\Gamma_{\text{top}}$, we enforce the following contact condition:

$$ d \cdot n_{\Gamma_{\text{sym}}} - g \le 0 \quad \text{on} \quad \Sigma, \qquad (44) $$

where $n_{\Gamma_{\text{sym}}}$ denotes the exterior unit normal to $\Gamma_{\text{sym}}$ (see Figure 9 (a)) and $g: \Sigma \to \mathbb{R}^+$ refer to the gap function between $\Sigma$ and $\Gamma_{\text{sym}}$. The inequality constraint (44) is approximated via a penalty method (see, e.g.,$^{46}$), by inserting the following term in the solid discrete problem:

$$ \frac{\gamma_c E \epsilon_s}{h^2} ([d_h^n \cdot n_{\Gamma_{\text{sym}}} - g + \epsilon_h]_+, w_h)_\Sigma, \qquad (45) $$

where $[x]_+ \stackrel{\text{def}}{=} \max\{0, x\}$, $\gamma_c > 0$ is a (dimensionless) user-defined parameter and $\epsilon_h > 0$ is a contact tolerance. The contact parameters are given by $\epsilon_h = 0.01$ and $\gamma_c = 5 \cdot 10^{-3}$.

The fluid mesh has 16384 triangles and the solid 50 edges. The zoom on the both meshes is presented in Figure 9(b). The Nitsche parameter $\gamma = 100$ and the Ghost penalty coefficient $\gamma_g = 1$. The CIP stabilization parameters are $\gamma_v = \gamma_p = 10^{-2}$. As in the previous example, we consider three levels of time refinement $\tau \in \{(10^{-3}/2^i)\}_{i=0}^2$ for the comparisons.

For illustration purposes, we report in Figure 10 the velocity magnitude, with the corresponding contour lines, at two different instants. In Figures 10 (a), (b) and (c) are reported the solutions obtained at time $t = 0.6$, when the valve is supposed to get into contact with the upper wall and the fluid velocity decreases globally as consequence of the closing of the valve. The same comparison is performed at time $t = 1$ in Figures 10 (d), (e) and (f) in a situation where the valve is open and far from contact. In this case the flux is reestablishes and the velocity increases in the channel. Again, a very good agreement is observed between Algorithm 3 and the strongly coupled scheme. On the contrary, the loosely coupled scheme delivers an approximation that is far from the previous ones. Similar observations can be made from Figure 11, where we compare the pressure at the same instants the pressure. We can see the high pressure jump when the valve is getting in contact with the wall (Figures 11(a) and (b)), while at $t = 1$ the discontinuity between the two sides of the interface is weaker (see Figures 11(d) and (e)). Algorithm 3 reproduces
very well the pressure jump obtained with the strongly coupled scheme, though a small difference is visible at $t = 0.6$. Figure 11 (c) and (f) show the results obtained with the loosely coupled algorithm, which is unable to deliver reasonable approximations.

Finally, Figures 12 and 13 present the time history of the horizontal and vertical displacement, respectively, at the upper solid point B for the different levels of time refinement. The contact condition with the wall can be seen in Figures 13, whereas Figures 12 shows that the structure is sliding and bouncing over the top wall. These results clearly show that Algorithm 3 is able to capture the dynamics of the interface before and after contact with the upper wall. Only slightly differences are observed, in particular close to the contact instant, but which decrease with the time refinement. On the contrary, the loosely coupled is not able to reproduce the dynamics obtained with the strongly coupled scheme, even with finest time refinement. This illustrates the limitations of the loosely coupled scheme.

<table><caption>**Table 1** L²-error displacement between Strongly coupled and Semi-implicit algorithms at time $t = 1$.</caption><thead><tr><th>Time-step</th><th>L<sup>2</sup>-difference</th></tr></thead><tbody><tr><td>τ = 10<sup>-3</sup></td><td>5.98 · 10<sup>-4</sup></td></tr><tr><td>τ = 5 · 10<sup>-4</sup></td><td>3.87 · 10<sup>-4</sup></td></tr><tr><td>τ = 2.5 · 10<sup>-4</sup></td><td>2.99 · 10<sup>-4</sup></td></tr></tbody></table>

A more quantitative comparison is given in Table 1 which shows L²-difference of the displacement obtained with the strongly coupled scheme and Algorithm 3 at time $t = 0.9$. We can observe that the difference reduces with the time refinement.

## 4.3 | Vesicle in lid-driven cavity flow

The last example is the well-known lid-driven cavity test with an immersed capsule (see, e.g., 47,48,49,50). The fluid geometry is shown in Figure 14 (a) and consists in three rigid wall $Γ_w$ and a lid $Γ_{top}$ moving with tangential velocity. The domain $Ω$ is a square given by $[0, 0]×[0, 1]$. The reference solid configuration is instead a circle of center $(0.6, 0.5)$ and radius $0.2$. The physical parameters used for the fluid in this test are $ρ_f = 100$, $μ = 10$. For the solid we have $ρ_s = 100$, $ε_s = 0.0212$, the Young's modulus $E = 5.6 · 10^3$ and Poisson's ratio $ν = 0.4$. Both fluid and solid are initially at rest. For the boundary conditions, we impose zero velocity on $Γ_w$ and tangential velocity of magnitude $\bar{u} = 1$ on $Γ_{top}$.

<figure><img src="image_16.png" /><figcaption>**Figure 14** (a) Geometric configuration and external boundary condition, (b) Fluid and solid meshes at time $t=0$.</figcaption></figure>

Since the fluid is entirely enclosed by Dirichlet-type boundary condition, standard Dirichlet-Neumann partitioned procedures for the solution of the implicit step (43) are known to diverge. This is due to the ill-posedness of the fluid system which enforces a volumetric constraint on the interfacial solid velocity. In order to avoid this issues, we consider the approach proposed in⁵¹ which consists in enforcing the volumetric constraint into the structure equation using a scalar Lagrange multiplier, in order to avoid the incompressibility incompatibilities with the subsequent fluid problem. More precisely, in step (43) of Algorithm 3 the
<figure><img src="image_17.png" /><figcaption>Figure 15 Velocity magnitude snapshots with $\tau = 5 \cdot 10^{-3}$.</figcaption></figure>

solid problem is replaced by the following mixed system

$$ \left\{ \begin{array}{l} \rho^s \epsilon (\partial_\tau \dot{d}_h^n, w_h)_\Sigma + a^s (d_h^n, w_h) + \lambda^n \int_{\Sigma^{n-1}} \hat{w}_h \cdot n = \frac{\gamma \mu}{h} \sum_{i=1}^2 (\tilde{u}_{h,i}^n - \dot{d}_{h,i}^{n-1}), \hat{w}_h)_{\Sigma^{n-1}} - \sum_{i=1}^2 (\sigma(\tilde{u}_{h,i}^n, p_{h,i}^n)n_i, \hat{w}_h)_{\Sigma^{n-1}}, \\ \int_{\Sigma^{n-1}} \dot{d}_h^n \cdot n = 0 \end{array} \right. $$

for all $w_h \in W_h$. The scalar Lagrange multiplier $\lambda^n \in \mathbb{R}$ represents the unknown constant pressure in the fluid.

The fluid mesh has 3200 elements while the solid mesh is composed by 80 edges. Both meshes are presented in Figure 14 (b). For the remaining parameters, we take the Nitsche parameter $\gamma = 1$, the Ghost penalty coefficient $\gamma_g = 1$ and the CIP stabilization parameters $\gamma_v = \gamma_p = 10^{-2}$. The chosen time steps are $\tau \in \{(5 \cdot 10^{-3}/2^i)\}_{i=0}^2$ and the final time is $T = 10$.

Always for illustration purposes, Figure 15 shows the snapshots of the fluid magnitude for both schemes and at different time, i.e., two different positions of the structure. At $t = 2$, Figures 15(a) and (b), the vesicle is starting the upper region of the cavity, while at $t = 5$, Figures 15(d) and (e), is moving away this region. Even with the coarsest time step, Algorithm 3 is able to predict the same location as the strongly coupled scheme. Figures 15(c) and (f) show the results of the loosely coupled scheme, which is clearly unable to reproduce the previous dynamics. Similar observations can be made from Figure 16 which shows the elevation of the pressure for the same time instants. No notable differences can be seen between Algorithm 3 and the strongly coupled scheme, whereas a major mismatch is obtained with the loosely coupled algorithm, particularly at time $t = 5$.

Finally, in Figure 17, we present the trajectory of the vesicle rightmost node obtained with the three numerical methods for the different levels of time refinement. Algorithm 3 is able to capture the dynamics of the strongly coupled scheme at all the discretization levels, whereas the loosely coupled scheme requires a sufficiently small time-step to deliver a minimally reasonable approximation.
<figure><img src="image_18.png" /><figcaption>Figure 16 Pressure snapshots with $\tau = 5 \cdot 10^{-3}$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17 Trajectory of the extreme right node of the vesicle from $t = 0$ to $t = 10$.</figcaption></figure>

## 5 | CONCLUSION

In this paper, we have introduced a new semi-implicit coupling scheme for the numerical approximations of incompressible fluid-structure interaction problems involving immersed solids. The proposed method generalizes the projection based semi-implicit coupling paradigm of$^{27}$ to the Nitsche-XFEM framework, with the following main ingredients:

* The traditional accuracy issues of previous splitting schemes based on a Nitsche's interface treatment (see$^{32,19}$) are circumvented through a variationally consistent transfer of the fluid viscous stresses to the solid problem;

* Consistent spatial approximation of the pressure-Poisson problem through the Nitsche-XFEM unfitted framework.
Moreover, Theorem 1 has shown that the method preserves the stability properties of the original splitting in the fitted mesh framework (see$^{27}$). The comprehensive numerical study reported in Section 4 confirmed these findings and showed a very good performance, in terms of stability and accuracy, with respect to the previous strongly coupled and loosely coupled schemes reported in$^{19}$. As a result, the present semi-implicit coupling scheme can be considered as a robust approach to avoid strong coupling in unfitted meshes without compromising stability and accuracy.

## References

1. Pozrikidis C.. *Computational hydrodynamics of capsules and biological cells*. Chapman & Hall/CRC Mathematical and Computational Biology CRC Press; 2010.

2. Hou G., Wang J., Layton A.. Numerical methods for fluid-structure interaction—a review. *Commun. Comput. Phys.* 2012;12(2):337–377.

3. Peskin C.S.. The immersed boundary method. *Acta Numer.* 2002;11:479–517.

4. Newren E.P., Fogelson A.L., Guy R.D., Kirby R.M.. Unconditionally stable discretizations of the immersed boundary equations. *J. Comput. Phys.* 2007;222(2):702–719.

5. Boffi D., Cavallini N., Gastaldi L.. Finite element approach to immersed boundary method with different fluid and solid densities. *Math. Models Methods Appl. Sci.* 2011;21(12):2523–2550.

6. Glowinski R., Pan T.W., Hesla T.I., Joseph D.D.. A distributed Lagrange multiplier/fictitious domain method for particulate flows. *Int. J. of Multiphase Flow*. 1999;25:755-794.

7. Baaijens F.. A fictitious domain/mortar element method for fluid-structure interaction. *Int. Jour. Num. Meth. Fluids*. 2001;35:743-761.

8. De Hart J., Peters G.W.M., Schreurs P.J.G., Baaijens F.P.T.. A three-dimensional computational analysis of fluid-structure interaction in the aortic valve. *Journal of Biomechanics*. 2003;36:103–112.

9. Astorino M., Gerbeau J.-F., Pantz O., Traoré K.-F.. Fluid-structure interaction and multi-body contact: application to aortic valves. *Comput. Methods Appl. Mech. Engrg.* 2009;198(45-46):3603–3612.

10. Boffi Daniele, Cavallini Nicola, Gastaldi Lucia. The finite element immersed boundary method with distributed Lagrange multiplier. *SIAM Journal on Numerical Analysis*. 2015;53(6):2584–2604.

11. Boffi Daniele, Gastaldi Lucia. A fictitious domain approach with Lagrange multiplier for fluid-structure interactions. *Numer. Math.* 2017;135(3):711–732.

12. Cottet G.-H., Maitre E., Milcent T.. Eulerian formulation and level set models for incompressible fluid-structure interaction. *M2AN Math. Model. Numer. Anal.* 2008;42(3):471–492.

13. Richter T.. A fully Eulerian formulation for fluid-structure-interaction problems. *J. Comput. Phys.* 2013;233:227–240.

14. Hachem E., Feghali S., Codina R., Coupez T.. Immersed stress method for fluid-structure interaction using anisotropic mesh adaptation. *Internat. J. Numer. Methods Engrg.* 2013;94(9):805–825.

15. Zilian A., Legay A.. The enriched space-time finite element method (EST) for simultaneous solution of fluid-structure interaction. *Internat. J. Numer. Methods Engrg.* 2008;75(3):305–334.

16. Gerstenberger Axel, Wall W.A.. An extended finite element method/Lagrange multiplier based approach for fluid-structure interaction. *Comput. Methods Appl. Mech. Engrg.* 2008;197(19-20):1699–1714.

17. Sawada T., Tezuka A.. LLM and X-FEM based interface modeling of fluid-thin structure interactions on a non-interface-fitted mesh. *Comput. Mech.* 2011;48(3):319–332.
18. Burman E., Fernández M.A. An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes. *Comput. Methods Appl. Mech. Engrg.* 2014;279:497-514.

19. Alauzet F., Fabrèges B., Fernández M.A, Landajuela M., Nitsche-XFEM for the coupling of an incompressible fluid with immersed thin-walled structures. *Computer Methods in Applied Mechanics and Engineering*. 2015;301.

20. Zonca S., Vergara C., Formaggia L. An unfitted formulation for the interaction of an incompressible fluid with a thick structure via an XFEM/DG approach. *SIAM J. Sci. Comput.* 2018;40(1):B59-B84.

21. Boilevin-Kayl L., Fernández M.A, Gerbeau J.F. Numerical methods for immersed FSI with thin-walled structures. *Computers & Fluids*. 2018;

22. Annese M. Time integration schemes for fluid-structure interaction problems: non-fitted FEMs for immersed thin structures. PhD thesisUniversità degli studi di Brescia2017.

23. Kim W., Lee H. A weak-coupling immersed boundary method for fluidâĂŞstructure interaction with low density ratio of solid to fluid. *J. Comput. Phys.* 2018;359:296-311.

24. Boilevin-Kayl L., Fernández M. A., Gerbeau J.-F. A loosely coupled scheme for fictitious domain approximations of fluid-structure interaction problems with immersed thin-walled structures. *SIAM J. Sci. Comput.* 2019;41(2):B351-B374.

25. Kadapa C., Dettmer W.G., Perić D. A stabilised immersed framework on hierarchical b-spline grids for fluid-flexible structure interaction with solid-solid contact. *Comput. Methods Appl. Mech. Engrg.* 2018;335:472-489.

26. Fernández M.A., Landajuela M. Splitting schemes for incompressible fluid/thin-walled structure interaction with unfitted meshes. *C. R. Math. Acad. Sci. Paris, 2015*. 2015; DOI: 10.1016/j.crma.2015.04.003.

27. Fernández M.A., Gerbeau J.F., Grandmont C. A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid. *Int. J. Num. Meth. Engrg.* 2007;69(4):794-821.

28. Roshchenko A., Minev P. D., Finlay W. H. A time splitting fictitious domain algorithm for fluidâĂŞstructure interaction problems (A fictitious domain algorithm for FSI). *Journal of Fluids and Structures*. 2015;58:109-126.

29. Landajuela M., Vidrascu M., Chapelle D., Fernández M. A. Coupling schemes for the FSI forward prediction challenge: comparative study and validation. *Int. J. Numer. Methods Biomed. Eng.* 2017;33(4):e02813, 23.

30. Quaini A., Quarteroni A. A semi-implicit approach for fluid-structure interaction based on an algebraic fractional step method. *Math. Models Methods Appl. Sci.* 2007;17(6):957-983.

31. Badia S., Quaini A., Quarteroni A. Splitting methods based on algebraic factorization for fluid-structure interaction. *SIAM J. Sci. Comput.* 2008;30(4):1778-1805.

32. Astorino M., Chouly F., Fernández M. A. Robin Based Semi-Implicit Coupling in Fluid-Structure Interaction: Stability Analysis and Numerics. *SIAM J. Sci. Comput.* 2009;31(6):4041-4065.

33. Astorino M., Grandmont C. Convergence analysis of a projection semi-implicit coupling scheme for fluid-structure interaction problems. *Numer. Math.* 2010;116:721-767.

34. Burman E., Hansbo P. Edge stabilization for the generalized Stokes problem: A continuous interior penalty method. *Computer Methods in Applied Mechanics and Engineering*. 2006;195:2393-2410.

35. Burman E. Ghost penalty. *C. R. Math. Acad. Sci. Paris*. 2010;348(21-22):1217-1220.

36. Massing André, Larson Mats G., Logg Anders. Efficient implementation of finite element methods on nonmatching and overlapping meshes in three dimensions. *SIAM J. Sci. Comput.* 2013;35(1):C23-C47.

37. Di Pietro D.A., Ern A. *Mathematical aspects of discontinuous Galerkin methods* Mathematics & Applications, vol. 69: Springer, Heidelberg; 2012.
38. Burman E., Fernández M.A., Hansbo P.. Continuous interior penalty finite element method for Oseen's equations. *SIAM J. Numer. Anal.* 2006;44(3):1248-1274.

39. Burman E., Fernández M.A.. Continuous interior penalty finite element method for the time-dependent Navier-Stokes equations: space discretization and convergence. *Numer. Math.* 2007;107(1):39–77.

40. Fries Thomas-Peter, Zilian Andreas. On time integration in the XFEM. *International Journal for Numerical Methods in Engineering*. 2009;79:69 - 93.

41. Gil A., Bonet J., Hassan O.. An enhanced Immersed Structural Potential Method for fluid-structure interaction. *Journal of Computational Physics*. 2013;250:178–205.

42. Hesch C., Gil A., Arranz Carreño A., Bonet J.. On continuum immersed strategies for Fluid-Structure Interaction. *Computer Methods in Applied Mechanics and Engineering*. 2012;247–248:51–64.

43. Wick T.. Flapping and Contact FSI Computations with the Fluid-Solid Interface-Tracking/Interface-Capturing Technique and Mesh Adaptivity. *Computational Mechanics*. 2013;53.

44. Kamensky D., Hsu M., Schillinger D., et al. An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves. *Computer methods in applied mechanics and engineering*. 2015;284:1005-1053.

45. Boilevin-Kayl L., Fernández M.A, Gerbeau J.F.. A Loosely Coupled Scheme for Fictitious Domain Approximations of Fluid-Structure Interaction Problems with Immersed Thin-Walled Structures. *SIAM Journal on Scientific Computing*. 2019;41.

46. Scholz R. Numerical solution of the obstacle problem by the penalty method. *Computing*. 1984;32:297-306.

47. Roy S., Heltai L., Costanzo F.. Benchmarking the Immersed Finite Element Method for Fluid-Structure Interaction Problems. *Computers and Mathematics with Applications*. 2013;69.

48. Wang X., Zhang L.. Interpolation functions in the immersed boundary and finite element methods. *Computational Mechanics*. 2010;45:321-334.

49. E Griffith Boyce. On the Volume Conservation of the Immersed Boundary Method. *Communications in Computational Physics*. 2012;12.

50. Griffith Boyce, Luo Xiaoyu, H Charney Leon. Hybrid finite difference/finite element version of the immersed boundary method. *International journal for numerical methods in biomedical engineering*. 2016;33.

51. Küttler Ulrich, Förster Christiane, Wall Wolfgang A.. A Solution for the Incompressibility Dilemma in Partitioned Fluid-Structure Interaction with Pure Dirichlet Fluid Domains. *Computational Mechanics*. 2006;38(4):417–429.

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