# Semi-implicit Eulerian method for the fluid structure interaction of elastic membranes

Mirco Ciallella, Thomas Milcent

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

## To cite this version:

Mirco Ciallella, Thomas Milcent. Semi-implicit Eulerian method for the fluid structure interaction of elastic membranes. 2023. [hal-04272884](https://hal.science/hal-04272884v1)

HAL Id: [hal-04272884](https://hal.science/hal-04272884v1)

Preprint submitted on 6 Nov 2023

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# SEMI-IMPLICIT EULERIAN METHOD FOR THE FLUID STRUCTURE INTERACTION OF ELASTIC MEMBRANES *

MIRCO CIALLELLA† AND THOMAS MILCENT‡

**Abstract.** In this paper we propose a novel and general approach to design semi-implicit methods for the simulation of fluid-structure interaction problems in a fully Eulerian framework. In order to properly present the new method, we focus on the two-dimensional version of the general model developed to describe full membrane elasticity. The approach consists in treating the elastic source term by writing an evolution equation on the structure stress tensor, even if it is nonlinear. Then, it is possible to show that its semi-implicit discretization allows us to add to the linear system of the Navier-Stokes equations some consistent dissipation terms that depend on the local deformation and stiffness of the membrane. Due to the linearly implicit discretization, the approach does not need iterative solvers and can be easily applied to any Eulerian framework for fluid-structure interaction. Its stability properties are studied by performing a Von Neumann analysis on a simplified one-dimensional model and proving that, thanks to the additional dissipation, the discretized coupled system is unconditionally stable. Several numerical experiments are shown for two-dimensional problems by comparing the new method to the original explicit scheme and studying the effect of structure stiffness and mesh refinement on the membrane dynamics. The newly designed scheme is able to relax the time step restrictions that affect the explicit method and reduce crucially the computational costs, especially when very stiff membranes are under consideration.

**Key words.** Eulerian elasticity, Fluid-structure interaction, Semi-implicit, Immersed boundary, Incompressible Navier-Stokes

**MSC codes.** 65M08, 65M12, 65M85, 76D05, 74B20

**1. Introduction.** Fluid-structure interaction problems arise in a wide range of applications and their numerical simulation is extremely challenging. In biomechanics, the correct modeling and computation of membranes immersed in an incompressible fluid is a key point to study biological capsules and cells, like red blood cells (RBCs). Capsules are in general modeled as a liquid drop protected by a thin elastic membrane and their use is widely spread in food, cosmetics and pharmaceutical industries. Simulations of such fluid-structure problems are expensive and require the correct coupling between the flow behavior and the structure behavior, which often comes with the imposition of some kind of interface conditions. One of the most classical approaches is the so-called Arbitrary Lagrangian-Eulerian (ALE) method [31, 30, 42, 13, 16]. This approach consists in solving the fluid-structure problem on a moving mesh that follows the displacement of the interface, and enforce the coupling conditions right on the grid points that describe the interface. Although this philosophy may bring accurate results for some problems, its limitations are clear when dealing with large deformations or topology changes.

In [39], Peskin introduced a new approach that simplifies the fluid-structure coupling for applications with structures modeled as lower dimensional manifolds of the problem. For immersed boundaries [40, 29, 35], the coupling interface conditions reduce to a forcing term in the Navier-Stokes equations that modifies the flow behavior by taking into account the deformation of the immersed membrane. In this case, the

<footnote>*Submitted to the editors 27/10/2023.</footnote>

<footnote>**Funding:** This work was funded by the ANR CAPSEULERIAN-FSI.</footnote>

<footnote>†École Nationale Supérieure d'Arts et Métiers, Institut de Mécanique et d'Ingénierie, 33400 Talence, France (mirco.ciallella@ensam.eu).</footnote>

<footnote>‡École Nationale Supérieure d'Arts et Métiers, Institut de Mécanique et d'Ingénierie, 33400 Talence, France (thomas.milcent@ensam.eu).</footnote>
structure is tracked using Lagrangian markers and the source term is spread on the fluid with a discretized Dirac delta. When large deformations are under consideration the need to insert or delete markers for stretching and shrinking of the membrane remains a delicate aspect, which requires the introduction of additional parameters to tune.

Both ALE and immersed boundary (IB) methods belong to the class of front-tracking methods that need an explicit definition and displacement of the interface. Fully Eulerian models for fluid-structure interaction have been recently introduced [7] for computational purposes as they allow to simulate quite easily multi-dimensional problems with large deformations thanks to an implicit treatment of interfaces and deformations through level-set fields [38, 36]. In the last decade, Eulerian elasticity has drawn a lot of attention and its implementation turned out to be very successful for both incompressible [9, 32, 41, 34, 12, 4] and compressible [21, 10, 11] applications. These models rely on the so-called backward characteristics which are defined as the inverse map of the deformation, and correlates the location of a material point to its initial position. As mentioned above, this vector is also defined with Eulerian variables that are advected by the Eulerian fluid velocity, and its spatial derivatives are used to compute the elastic stress tensor. Rather than following the backward characteristics, it was proposed in [44] an Eulerian approach based on following the evolution of the symmetrical left Cauchy-Green tensor through an inhomogeneous advection equation with extra terms depending on the velocity gradient. In this case, six components in three dimensions (three in two dimensions) must be advected, rather than three in three dimensions (two in two dimensions) for the backward characteristics.

Whether one deals with ALE, IB or fully Eulerian models for incompressible fluid-structure interaction problems, stability is a well-known issue of monolithic frameworks in which both the fluid and the solid are solved at the same time. For Eulerian models this was proven to be related to the velocity of elastic waves that (for stiff materials) are order of magnitude higher than that of the fluid, which leads to dramatic time step restrictions [10]. In order to overcome such limitations related to monolithic frameworks, several works have been dedicated to the derivation of implicit schemes [47, 33, 43, 37]. However, due to the strong nonlinearity of the problem, these schemes require computationally expensive and cumbersome iterative approaches. With the goal of moving towards methods that are easier to introduce in different computational frameworks, semi-implicit or approximate implicit schemes are a very interesting and widespread research topic since their goal is to provide a less expensive and efficient alternative to standard implicit methods [47, 43, 17, 18]. Similar interest in such schemes was also shown in the context of multi-phase flows for the computation of surface tension [25, 45] to avoid the time step restriction proportional to $\Delta x^{3/2}$.

In the context of fully Eulerian models for fluid-structure interaction, the attempts to relax the time step restrictions have been proposed in [26, 8] with methods that seem to be tailored to the application under consideration. In particular, the first reference [26] addresses the simulation of solid bulks immersed in an incompressible fluid by introducing a linear expansion for the nonlinear elastic stress with a Jacobian tensor, to improve the performances of the scheme initially proposed in [44]. Whereas in [8] a semi-implicit approach for thin membranes subjected to area variation is developed by solving a diffusion equation to predict the interface position, in the context of the model proposed in [7]. In this work, we aim at proposing a general way of deriving semi-implicit schemes for Eulerian models for fluid-structure interaction that can be directly applicable to general (nonlinear) elastic models without any need
of linearization. This approach consists in discretizing the evolution equation of the stress tensor in a semi-implicit way, and writing the value of the stress at the new time level as a function of that at the previous time plus some additional terms. In order to show the generality and potential of this approach, we applied it to the two-dimensional version of the full membrane model introduced in [34], which is capable of dealing with both area and shear variations. While the level-set approach proposed in [8] only works for the elastic force that models the area variation introduced in [7], our idea can be ideally applied to more general and nonlinear elastic models, because it depends on the deformation vector. Indeed, the model treated in this paper is the starting point to develop semi-implicit schemes for full membrane elasticity that includes both area and shear variations [34]. The complete three-dimensional analysis will be tackled in a future work.

The paper is organized as follows. In Section 2, we recall the basic definitions to derive the fully Eulerian model for fluid-structure interaction focusing on backward characteristics, level-set and the membrane elastic model treated herein. In Section 3 we recall and summarize the main stability results obtained for a simplified one-dimensional model, and derive our semi-implicit scheme in this framework. Here it is shown that the new scheme consists in adding a consistent dissipation to the system which allows us to prove its unconditional stability. In Section 4, we begin by describing the classical explicit coupling of the elastic source term and the general algorithm to solve the full model. Here, we also derive an evolution equation on the stress tensor that will be then used to design the semi-implicit approach. After that, we focus on the derivation of the semi-implicit coupling which is performed by introducing the discretized evolution equation of the stress in the momentum equation. Here we show that also the approach for the full model consists in adding a diffusion that depends on stiffness and area variation. In Section 5, we test the convergence properties of the discretized semi-implicit operators to show that they are second-order accurate in space. In Section 6, we provide numerical illustrations on the ability of the newly developed semi-implicit scheme to introduce an important relaxation on the time step restriction to simulate a membrane immersed in a simple shear flow. In particular, we focus on the influence of membrane stiffness and mesh refinement on the computational gain and quality of the solution. In Section 7, we conclude with a summary of the work and give some perspectives on aspects that need further investigation. Finally, in Appendix A we give a proof on the parameter chosen to measure the area variation, in Appendix B we give the proof of the evolution equation of the stress tensor, and in Appendix C we provide examples on how to discretize the semi-implicit operators on Cartesian staggered meshes.

## 2. Fully Eulerian Model.

**2.1. Forward and backward characteristics.** Let $\Omega_0 \subset \mathbb{R}^2$ be the reference configuration of a continuous medium and assume that this medium is deformed by a smooth map $X: \mathbb{R}^+ \times \mathbb{R}^2 \to \mathbb{R}^2$ (the forward characteristics) to $\Omega_t = X(t, \Omega_0)$. A velocity field $u: \mathbb{R}^2 \times \mathbb{R}^+ \to \mathbb{R}^2$ is naturally associated with $X$:

$$ (2.1) \qquad \partial_t X(t, \xi) = u(X(t, \xi), t), \quad X(0, \xi) = \xi, \quad \xi \in \Omega_0. $$

We introduce the backward characteristics $Y: \mathbb{R}^2 \times \mathbb{R}^+ \to \mathbb{R}^2$ by the formula $Y(X(t, \xi), t) = \xi$. The physical interpretation of $Y(x, t)$ is the position at time 0 of a material particle lying in $x$ at time $t$ and moving at speed $u$. The derivative of this
<figure><img src="image_3.png" /><figcaption>Fig. 1: Forward and backward characteristics.</figcaption></figure>

relation with respect to $t$ and $\xi$ in turn gives with (2.1)

$$ (2.2) \qquad \partial_t Y + (u \cdot \nabla_x) Y = 0, \quad Y(x, 0) = x, \quad x \in \Omega_t, $$

and

$$ (2.3) \qquad [\nabla_{\xi} X(t, \xi)] = [\nabla_x Y(x, t)]^{-1}, \quad \text{for } x = X(t, \xi). $$

The relation (2.2) is the Eulerian equivalent of the characteristic equation (2.1). In addition, equation (2.3) allows to compute the gradient of the deformation in the Eulerian frame via $Y$. The next sections are devoted to the description of the mathematical model used to treat the membrane elastic deformations with the backward characteristics $Y$. To simplify the notation, we drop the subscript of the operator $\nabla_x$.

**2.2. Hyperelastic membrane models.** The notations and results summarized in this section are detailed in [34] for the general three-dimensional model. We consider a surface $\Gamma_t = \{x \in \mathbb{R}^2 / \phi(x, t) = 0\}$ captured by a level set function $\phi: \mathbb{R}^2 \times \mathbb{R}^+ \to \mathbb{R}$ and advected by the Eulerian velocity field $u$:

$$ (2.4) \qquad \partial_t \phi + u \cdot \nabla \phi = 0. $$

The normal $n(x, t)$ for $x \in \Gamma_t$ is then expressed in terms of the normalized gradient of the level set:

$$ (2.5) \qquad n(x, t) = \frac{\nabla \phi(x, t)}{|\nabla \phi(x, t)|}. $$

To measure the deformations on the surface $\Gamma_t$ we introduce the tensor

$$ (2.6) \qquad \mathcal{A} = B - \frac{(Bn) \otimes (Bn)}{(Bn) \cdot n}. $$

where $B = [\nabla Y]^{-1}[\nabla Y]^{-T}$ is the left Cauchy-Green tensor. Tensor $\mathcal{A}$ measures the surface deformations by projecting the deformations (measured by $B$) on the surface
$\Gamma_t$ (represented locally by $n$). The vector $n$ is an eigenvector of $\mathcal{A}$ associated to the eigenvalue 0 so $\det(\mathcal{A}) = 0$. In two dimensions, the other invariant (the trace) is used to define

$$ (2.7) \qquad Z = \sqrt{\mathrm{Tr}(\mathcal{A})}. $$

Following the same approach in [34], it is possible to prove that $Z$ measures the local area variation for two-dimensional problems (see Appendix A for the proof).

**Remark 2.1.** Using the backward characteristics $Y$, rather than only the level-set $\phi$ as in [7], allows us to have the deformations that can also be used to measure the shear variation for tree-dimensional membranes. Indeed, the transport equation for the level-set function only records information on the surface area variations, and ignores any tangential component of $u$. This is due to the fact that $\nabla\phi$ is normal to the interface.

In order to compute the elastic force, we then introduce the regularized membrane energy,

$$ (2.8) \qquad \mathcal{E} = \int_Q E(Z) \frac{1}{\varepsilon} \zeta \left( \frac{\phi}{\varepsilon} \right) dx. $$

Here $Q$ is a box containing the membrane, $E(Z)$ is the constitutive law, $\varepsilon$ the width of the interface and $\zeta$ is a cut-off function used to spread the interface near $\{\phi = 0\}$. The associated force, computed with the principle of virtual power, is given by

$$ (2.9) \qquad F = \mathrm{div}(\sigma) \delta_{\varepsilon}(\phi), \quad \text{where} \quad \sigma = E'(Z)Z(\mathbb{I} - n \otimes n) $$

where

$$ \delta_{\varepsilon}(\phi) = \frac{1}{\varepsilon} \zeta \left( \frac{\phi}{\varepsilon} \right) $$

is introduced to simplify the notation.

The Evan-Skalak constitutive law is used in this article:

$$ (2.10) \qquad E'(Z) = K(Z - 1), $$

where $K$ is the elastic modulus of the membrane.

**2.3. Fluid-structure interaction model.** The elastic membrane is immersed in a incompressible fluid modeled by the Navier-Stokes equations. The overall fully Eulerian model is given by

$$ (2.11) \qquad \left\{ \begin{array}{l} \rho(\partial_t u + (u \cdot \nabla)u) + \nabla p - \mathrm{div}(2\mu D(u)) = F(\phi, Y), \\ \mathrm{div}(u) = 0, \\ \partial_t \phi + u \cdot \nabla \phi = 0, \\ \partial_t Y + (u \cdot \nabla)Y = 0, \end{array} \right. $$

where $u$ is the fluid velocity, $D(u) = ([\nabla u] + [\nabla u]^T)/2$ the strain rate tensor, $\rho$ the density, $p$ the pressure and $\mu$ the viscosity. Different values of density and viscosity, outside (1) and inside (2) the membrane, are taken into account by defining

$$ \rho = \mathcal{H} \left( \frac{\phi}{\varepsilon} \right) \rho_1 + \left( 1 - \mathcal{H} \left( \frac{\phi}{\varepsilon} \right) \right) \rho_2, \quad \mu = \mathcal{H} \left( \frac{\phi}{\varepsilon} \right) \mu_1 + \left( 1 - \mathcal{H} \left( \frac{\phi}{\varepsilon} \right) \right) \mu_2, $$
where $\mathcal{H}$ represents a smooth Heaviside function. The membrane force $F$ is given by (2.9). These equations are completed with appropriate initial and boundary conditions that will be detailed in the numerical validation section.

The fluid-structure problem is governed by two parameters, the Reynolds number $Re$ and the capillary number $Ca$, that can be recovered by rewriting the momentum equation in a non-dimensional form. Therefore, we introduce the reference variables $\varphi^*$ and the non-dimensional variables $\bar{\varphi}$, which give the following operators and physical variables:

$$ \rho = \rho^* \bar{\rho}, \quad u = u^* \bar{u}, \quad p = p^* \bar{p}, \quad \mu = \mu^* \bar{\mu}, \quad F = \frac{K^*}{(L^*)^2} \bar{F}, \quad \nabla = \frac{\bar{\nabla}}{L^*}. $$

By replacing all variables with their counterpart, and taking $p^* = \mu^* u^*/L^*$, the non-dimensional momentum equation reads,

$$ (2.12) \qquad \mathrm{Re} \left[ \bar{\rho} \left( \partial_t \bar{u} + (\bar{u} \cdot \bar{\nabla}) \bar{u} \right) \right] + \bar{\nabla} \bar{p} - \mathrm{div} \left( \bar{\mu} \left[ (\bar{\nabla} \bar{u}) + [\bar{\nabla} \bar{u}]^T \right) \right) = \frac{1}{\mathrm{Ca}} \bar{F}, $$

where $\mathrm{Re} = \rho^* u^* L^*/\mu^*$, and $\mathrm{Ca} = \mu^* u^*/K^*$.

### 2.4. Reinitialization and extrapolation.
It is a well-known problem that the level-set field $\phi$ stops being a signed distance function after its time evolution with the associated advection equation. In this work, for its reinitialization, we follow the idea proposed in [36] and solve for a fictitious time $\tau$

$$ (2.13) \qquad \partial_{\tau} \phi + \operatorname{sgn}(\phi_0)(|\nabla \phi| - 1) = 0. $$

Moreover, as discussed in [12, 4] also the backward characteristics $Y$ can get extremely distorted with time. Being the elastic force computed by means of $Y$ and $\phi$, if not treated, these distortions are likely to cause spurious disturbances close to the fluid-structure interface.

To summarize previous works on this topic, we perform the linear extrapolation proposed in [3] for the outer characteristics, i.e. $\phi > 0$,

$$ (2.14) \qquad \partial_{\tau} Y_n + \mathcal{H} \left( \frac{\phi}{\varepsilon} \right) (n \cdot Y_n) = 0, \quad \text{with } Y_n = n \cdot \nabla Y, $$

$$ (2.15) \qquad \partial_{\tau} Y + \mathcal{H} \left( \frac{\phi}{\varepsilon} \right) (n \cdot Y - Y_n) = 0. $$

### 3. Stability analysis for a simplified one-dimensional model.
In this section, we recall recent stability results in the context of explicit and implicit coupling for a simplified model, and we introduce the new semi-implicit approach in this framework. We consider here the same linearization and assumptions taken for similar models [5, 20]. To summarize we set $\rho = 1$ and we decompose the solution as the small perturbation of a stationary solution $(u_0, Y_0) = (0, x)$.

By avoiding the treatment of nonlinear terms we end up with the following linearized one-dimensional model:

$$ (3.1) \qquad \left\{ \begin{aligned} \frac{\partial u}{\partial t} - \mu \frac{\partial^2 u}{\partial x^2} &= -\frac{K}{\varepsilon} \frac{\partial^2 Y}{\partial x^2}, \\ \frac{\partial Y}{\partial t} + u &= 0. \end{aligned} \right. $$
Here we recall the stability conditions for such models when using a classical explicit coupling similarly to what was introduced in [5]. For simplicity we consider a classical centered discretization of second order derivatives, with grid size $\Delta x$, with an implicit treatment for the viscous term:

$$ (3.2) \qquad \begin{cases} \frac{u_j^{n+1} - u_j^n}{\Delta t} - \mu \frac{u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}}{(\Delta x)^2} = -\frac{K Y_{j+1}^n - 2Y_j^n + Y_{j-1}^n}{\Delta x}, \\ \frac{Y_j^{n+1} - Y_j^n}{\Delta t} + u_j^{n+1} = 0. \end{cases} $$

For this system, the stability condition was proven in [5] to be

$$ (3.3) \qquad \Delta t < \frac{\mu\varepsilon + \max(\mu\varepsilon, \sqrt{K\varepsilon\Delta x})}{K}. $$

Instead, when treating implicitly the elastic term, the discretized system reads:

$$ (3.4) \qquad \begin{cases} \frac{u_j^{n+1} - u_j^n}{\Delta t} - \mu \frac{u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}}{(\Delta x)^2} = -\frac{K Y_{j+1}^{n+1} - 2Y_j^{n+1} + Y_{j-1}^{n+1}}{\Delta x}, \\ \frac{Y_j^{n+1} - Y_j^n}{\Delta t} + u_j^{n+1} = 0. \end{cases} $$

Once again in [5], this system was proven to be unconditionally stable. However, for real applications, this scheme is very cumbersome to design and implement due to the full fluid-structure Navier-Stokes system with nonlinear source term and advection equations.

For this reason, following the same motivation of [8], we introduce and study the stability analysis of the new semi-implicit scheme in this simplified framework. The new approach developed for the full nonlinear model is described in Section 4.2.

We start by introducing an auxiliary variable $\tilde{Y}$ for the implicit elastic term,

$$ (3.5) \qquad \begin{cases} \frac{u_j^{n+1} - u_j^n}{\Delta t} - \mu \frac{u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}}{(\Delta x)^2} = -\frac{K \tilde{Y}_{j+1}^{n+1} - 2\tilde{Y}_j^{n+1} + \tilde{Y}_{j-1}^{n+1}}{(\Delta x)^2}, \\ \frac{Y_j^{n+1} - Y_j^n}{\Delta t} + u_j^{n+1} = 0. \end{cases} $$

For the simplified model, writing an evolution equation on the stress tensor, discretized in a semi-implicit manner, entails that $\tilde{Y}_j^{n+1} = Y_j^n - \Delta t u_j^{n+1}$. This gives us an explicit part and an implicit one that depends on the second order derivatives of $u$,

$$ \frac{\tilde{Y}_{j+1}^{n+1} - 2\tilde{Y}_j^{n+1} + \tilde{Y}_{j-1}^{n+1}}{(\Delta x)^2} = \frac{Y_{j+1}^n - 2Y_j^n + Y_{j-1}^n}{(\Delta x)^2} - \Delta t \frac{u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}}{(\Delta x)^2}. $$

When introducing the previous relation in System (3.5), the overall system reads

$$ (3.6) \qquad \begin{cases} \frac{u_j^{n+1} - u_j^n}{\Delta t} - \left(\mu + \Delta t \frac{K}{\varepsilon}\right) \frac{u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}}{(\Delta x)^2} = -\frac{K Y_{j+1}^n - 2Y_j^n + Y_{j-1}^n}{\Delta x}, \\ \frac{Y_j^{n+1} - Y_j^n}{\Delta t} + u_j^{n+1} = 0. \end{cases} $$
It should be noticed that this approach reduces to solving a PDE similar to the one arising from the explicit discretization, but with an additional consistent viscosity term $\Delta t \frac{K}{\epsilon}$ that depends on the time step and the elastic parameter of the structure.

**PROPOSITION 3.1.** *The semi-implicit scheme (3.6) is unconditionally stable.*

*Proof.* With the assumption of periodic boundary conditions we perform a Von-Neumann stability analysis. Therefore, we decompose the solution $(u, Y)$ as Fourier expansions of the grid values,

$$u_j^n = \hat{u}^n(k)e^{2i\pi kj\Delta x}, \quad Y_j^n = \hat{Y}^n(k)e^{2i\pi kj\Delta x} \xrightarrow{\theta:=2\pi k\Delta x} u_j^n = \hat{u}^n(\theta)e^{ij\theta}, \quad Y_j^n = \hat{Y}^n(\theta)e^{ij\theta}.$$

We then obtain from (3.6)

$$\left\{ \begin{array}{l} \hat{u}^{n+1}(\theta) \left[ 1 + 4 \frac{\Delta t}{(\Delta x)^2} \left( \mu + \Delta t \frac{K}{\epsilon} \right) \sin^2 \left( \frac{\theta}{2} \right) \right] = \hat{u}^n(\theta) + \hat{Y}^n(\theta) \left[ 4 \frac{\Delta t}{(\Delta x)^2} \frac{K}{\epsilon} \sin^2 \left( \frac{\theta}{2} \right) \right], \\ \\ \Delta t \hat{u}^{n+1}(\theta) + \hat{Y}^{n+1}(\theta) = \hat{Y}^n(\theta). \end{array} \right.$$

By setting

$$\alpha_{\theta} = 1 + 4 \frac{\Delta t}{(\Delta x)^2} \left( \mu + \Delta t \frac{K}{\epsilon} \right) \sin^2 \left( \frac{\theta}{2} \right), \quad \beta_{\theta} = 4 \frac{\Delta t}{(\Delta x)^2} \frac{K}{\epsilon} \sin^2 \left( \frac{\theta}{2} \right),$$

we can recast the system in matrix form

$$\underbrace{\begin{pmatrix} \alpha_{\theta} & 0 \\ \Delta t & 1 \end{pmatrix}}_{A_{\theta}} \begin{pmatrix} u^{n+1}(\theta) \\ \hat{Y}^{n+1}(\theta) \end{pmatrix} = \underbrace{\begin{pmatrix} 1 & \beta_{\theta} \\ 0 & 1 \end{pmatrix}}_{B_{\theta}} \begin{pmatrix} u^n(\theta) \\ \hat{Y}^n(\theta) \end{pmatrix}.$$

The eigenvalues $\lambda_1$ and $\lambda_2$ are the solution of the characteristic polynomial $\det(A_{\theta}^{-1}B_{\theta} - \lambda I) = 0$, which is

$$\det \begin{pmatrix} \frac{1}{\alpha_{\theta}} - \lambda & \frac{\beta_{\theta}}{\alpha_{\theta}} \\ -\frac{\Delta t}{\alpha_{\theta}} & 1 - \Delta t \frac{\beta_{\theta}}{\alpha_{\theta}} - \lambda \end{pmatrix} = \alpha_{\theta} \lambda^2 - \lambda(1 + \alpha_{\theta} - \Delta t \beta_{\theta}) + 1 = 0,$$

whose invariants are

$$\lambda_1 \lambda_2 = \frac{1}{\alpha_{\theta}} < 1, \quad \lambda_1 + \lambda_2 = \frac{1 + \alpha_{\theta} - \Delta t \beta_{\theta}}{\alpha_{\theta}} > 0.$$

If $(\lambda_1 + \lambda_2)^2 < 4\lambda_1\lambda_2$, the eigenvalues are complex conjugates, and their modulus is always less than 1. If $(\lambda_1 + \lambda_2)^2 > 4\lambda_1\lambda_2$, the eigenvalues are real and always positive. In particular,

$$\lambda_1 = \frac{1 + \alpha_{\theta} - \Delta t \beta_{\theta} + \sqrt{(1 + \alpha_{\theta} - \Delta t \beta_{\theta})^2 - 4\alpha_{\theta}}}{2\alpha_{\theta}} \\ \le \frac{1 + \alpha_{\theta} + \sqrt{(1 + \alpha_{\theta})^2 - 4\alpha_{\theta}}}{2\alpha_{\theta}} = 1,$$

which proves that the spectral radius of $A_{\theta}^{-1}B_{\theta}$ is always less than one. □

**Remark 3.2** (Stability of the fully nonlinear semi-implicit schemes). Observe that the above stability condition has been proven for a linear simplified one-dimensional model and it cannot be generally applied to the fully nonlinear fluid-structure problem. In particular, for realistic applications, a relaxation of the time step restriction is experienced thanks to the additional dissipation added depending on the local deformation of the media.
<figure><img src="image_4.png" /><figcaption>Fig. 2: Staggered grid with position of unknowns. ■: scalar field (p, ρ, Y, φ). ▶: vector field – component along x. ▲: vector field – component along y.</figcaption></figure>

**4. Numerical schemes.** In this section we present how to develop both explicit and semi-implicit couplings for the full nonlinear model.

**4.1. Explicit method.** In general, The equations (2.11) are discretized with finite volume schemes on a staggered grid (see Figure 4.1 for a two-dimensional configuration). It should be noticed that the force term $F(\phi, Y)$ is here written as the divergence of a stress tensor multiplied by a cut-off function.

Let $\Delta t$ be the time step and $u^n, p^n, \phi^n, Y^n, \rho^n, \mu^n$ the time discretization of the variables at $t_n = n\Delta t$. The discretization is given by the projection method (see [22] for an overview on projection methods) for the Navier Stokes equations and an explicit scheme for the advection equations for the interface and backward characteristics:

$$ (4.1) \quad \text{Step 1:} \quad \rho^n \left( \frac{u^* - u^n}{\Delta t} + \text{div}(u^n \otimes u^*) \right) + \nabla p^n - \text{div}(2\mu^n D(u^*)) = \delta_{\varepsilon}(\phi^n) \text{div}(\sigma^n), $$

$$ (4.2) \quad \text{Step 2:} \quad \text{div} \left( \frac{\Delta t}{\rho^n} \nabla \psi^{n+1} \right) = \text{div}(u^*), $$

$$ (4.3) \quad \text{Step 3:} \quad u^{n+1} = u^* - \frac{\Delta t}{\rho^n} \nabla \psi^{n+1}, \quad p^{n+1} = p^n + \psi^{n+1}, $$

$$ (4.4) \quad \text{Step 4:} \quad \frac{\phi^{n+1} - \phi^n}{\Delta t} + u^{n+1} \cdot \nabla \phi^n = 0, \quad \frac{Y^{n+1} - Y^n}{\Delta t} + (u^{n+1} \cdot \nabla)Y^n = 0. $$

In Step 1 (4.1) a prediction of the velocity $u^*$ is computed with an Euler implicit scheme for the viscous term, and an explicit treatment of the elastic source term $\sigma^n$. The convection term is treated following a semi-implicit approach, with a standard second-order discretization. In Step 2 (4.2), the Poisson equation for the pressure increment $\psi^{n+1}$ is solved with appropriate boundary conditions [22]. The resulting
linear systems are solved with the GMRES algorithm of the HYPRE library [15, 14] with preconditioning. In Step 3 (4.3) the velocity is corrected to enforce the incompressibility condition and the pressure is updated. In Step 4 (4.4) the transport equations are discretized with an explicit Runge-Kutta3 scheme in time and a WENO5 [27] scheme in space for $\phi$ and $Y$.

For the cut-off function, we considered the following expression:

$$ (4.5) \qquad \zeta(r) = \begin{cases} \frac{1}{2}(1 + \cos(\pi r)), & \text{on } [-1, 1] \\ 0, & \text{elsewhere.} \end{cases} $$

In our simulations $\varepsilon$ is fixed at $2\Delta x$ which is the standard value used in the literature to spread the interface.

**4.2. Semi-implicit method.** In the context of fluid-structure interaction problems, a classical explicit treatment of the elastic source term introduces very strict time step restrictions, especially when the structure is stiff. This restriction leads to long and cumbersome simulations, particularly when working with incompressible flows for which long times are needed to simulate realistic configurations, and the speed of elastic waves is much faster than the fluid speed.

For such problems, a good approximation of the time step is given by the surface tension condition introduced in [6], which limits the time step to $\sim \sqrt{\rho\Delta x^3/K}$. Although capturing properly all elastic waves becomes crucial for compressible materials in which transients are under study [10], in incompressible simulations one is more interested in the long-time dynamics of deformable structures. For instance when dealing with very stiff structures, the admissible time step will be extremely low even if the structure will not deform that much.

For this reason here we present the new semi-implicit method for thin elastic membranes, following the approach introduced in Section 3 for which it was possible to prove unconditional stability in a simplified framework. The general idea consists in writing the evolution equation of the stress tensor, and then discretizing it in a semi-implicit way with respect to the Eulerian velocity $u$. Therefore, we introduce here an evolution equation for the nonlinear membrane stress tensor $\sigma$ defined in (2.9).

**PROPOSITION 4.1.** *Under the smoothness assumption made on $u$, the membrane stress tensor $\sigma$ verifies*

$$ (4.6) \qquad \partial_t \sigma + u \cdot \nabla \sigma = (f'(Z)Z[\nabla u] : C)C + f(Z) ([\nabla u]^T(n \otimes n) + (n \otimes n)[\nabla u] - 2([\nabla u]n \cdot n)(n \otimes n)). $$

where we defined $f(Z) = E'(Z)Z$, and $C = I - n \otimes n$.

*Proof.* See Appendix B. □

In order to develop the new method for the full model, we start by replacing the explicit term $\sigma^n$ by $\sigma^*$, within the prediction step (4.1),

$$ (4.7) \qquad \rho^n \left( \frac{u^* - u^n}{\Delta t} + \text{div}(u^n \otimes u^*) \right) + \nabla p^n - \text{div}(2\mu^n D(u^*)) = \delta_{\varepsilon}(\phi^n) \text{div}(\sigma^*). $$

Once one has obtained the evolution equation for $\sigma$ (4.6), we discretize it in a semi-implicit way for all the terms that depend on $u$ as,

$$ (4.8) \qquad \frac{\sigma^* - \sigma^n}{\Delta t} + (u^* \cdot \nabla)\sigma^n = T(Z^n, \nabla u^*, n) + f(Z^n) ([\nabla u^*]^T(n \otimes n) + (n \otimes n)[\nabla u^*]) $$
where

$$ \mathcal{T}(Z^n, \nabla u^*, n) = (f'(Z^n)Z^n[\nabla u^*] : \mathcal{C}) \mathcal{C} - 2f(Z^n)([\nabla u^*]n \cdot n)(n \otimes n). $$

When replacing Equation (4.8) into (4.7) the fluid-structure momentum equation becomes

$$ \begin{aligned} & (4.9) \rho^n \left( \frac{u^* - u^n}{\Delta t} + \operatorname{div}(u^n \otimes u^*) \right) + \nabla p^n + \delta_{\varepsilon}(\phi^n) \Delta t \operatorname{div} \left[ (u^* \cdot \nabla) \sigma^n - \mathcal{T}(Z^n, \nabla u^*, n) \right] + \\ & \operatorname{div} \left[ \left( \mu^n \mathbb{I} + \delta_{\varepsilon}(\phi^n) \Delta t f(Z^n)(n \otimes n) \right) [\nabla u^*] + [\nabla u^*]^T \left( \mu^n \mathbb{I} + \delta_{\varepsilon}(\phi^n) \Delta t f(Z^n)(n \otimes n) \right) \right] \\ = & \delta_{\varepsilon}(\phi^n) \operatorname{div}(\sigma^n). \end{aligned} $$

It should be noticed that the consistent terms coming from the evolution equation of $\sigma$ add some viscosity to the classical one. Moreover, the considered elastic model also gives additional viscosity terms appearing in the $\mathcal{T}$ term of (4.9). All these additional discretized terms depend on the time step $\Delta t$, the geometric features of the membrane described by $n$, and its area variation measured by $Z$. These new viscosity terms differ significantly from $2\mu D(u)$ because $\mu$ is scalar, while the new terms, like $f(Z^n)(n \otimes n)$, are matrices. For this reason they are going to be referred to as *tensorial viscosity terms*.

**Remark 4.2** (Computational costs of the semi-implicit method). It should be noticed that the resolution of the momentum equation for the explicit scheme already deals with a linear system to treat implicitly and semi-implicitly the strain rate tensor and the advection term, respectively. Therefore, the additional computational cost needed to perform the semi-implicit computations is given by adding the tensorial viscosity terms to the linear system. As shown in Appendix C a wider stencil is needed with respect to the classical viscous terms. These differences do not impact the costs that much also because only those elements around the zero level-set require the tensorial viscosity, meaning when $\delta_{\varepsilon}(\phi^n) \neq 0$.

**5. Manufactured solution.** For validation purposes in this section we analyze the convergence properties of the semi-implicit terms, discussed in Section 4.2, discretized with standard second-order accurate formulas. In Appendix C we give further details on how to discretize these operators.

The following equation is solved

$$ 10u + \operatorname{div}((u \cdot \nabla)M - M[\nabla u] - [\nabla u]^T M - ([\nabla u] : M)M) = S, $$

which reduces to solving a linear system $\tilde{M}u_h = S$ in a square domain $[0, 1] \times [0, 1]$, where $S$ is computed by following the method of manufactured solutions and considering the exact solution

$$ (5.1) \quad u = \begin{bmatrix} \sin(\pi x) \sin(\pi y) \\ \sin(\pi x) \sin(\pi y) \end{bmatrix} \quad \text{with} \quad M = \begin{bmatrix} 1 + \sin(\pi x y) & 1 + \sin(\pi x y) \\ 1 + \sin(\pi x y) & 1 + \sin(\pi x y) \end{bmatrix}. $$

**Remark 5.1** (Manufactured solution). It should be noticed that an additional term $10u$ is added to the diagonal because solving only the semi-implicit operators gives rise to an ill-conditioned linear system which is never considered in practical
<figure><img src="image_5.png" /><figcaption>Fig. 3: Validation: convergence analysis with the discretized semi-implicit operators. $u_h$ and $u_{ex}$ represent the numerical and exact solution, respectively.</figcaption></figure>

applications. As a matter of fact, when solving the Navier-Stokes equation, we always have some additional terms to the diagonal given by the advection and diffusion parts of the system. This allows us to study the convergence of the discretized operators without having conditioning issues.

In Figure 5.1 the discretization error is computed with a $L_2$ norm on a set of nested meshes: 20 × 20, 40 × 40, 80 × 80, 160 × 160, 320 × 320. The convergence analysis shows that the discretized operators needed to design the semi-implicit scheme are indeed second-order accurate.

**6. Numerical experiments.** In this section, we study the effectiveness of the proposed scheme by comparing it to the original explicit one for a thin membrane immersed in a simple shear flow. The simulation is setup by considering a membrane of radius $a = 0.5$ immersed within a computational domain $[-4, 4] \times [-2, 2]$ discretized with uniform Cartesian meshes. From (2.12) the dimensionless parameters are modified by taking the reference speed equal to $a\dot{\gamma}$, and read

$$ (6.1) \qquad \mathrm{Re} = \frac{\rho a^2 \dot{\gamma}}{\mu} = 0.1, \qquad \mathrm{Ca} = \frac{\mu a \dot{\gamma}}{K}, $$

where the capillary number Ca, is modified according to the configuration chosen. For all cases, $\rho = \rho_1 = \rho_2$ and the angular speed $\dot{\gamma}$ are both set equal to 1. From the non-dimensional momentum Equation (2.12) it is possible to notice that when Ca tends to zero, the elastic source term becomes stiff. This also relates to what happens when dealing with low Mach flows or materials with the compressible Euler equations [1, 2, 46]. For these experiments, four capillary numbers are considered: 0.02, 0.01, 0.008 and 0.001. It should be noticed that when dealing with Ca = 0.001, the explicit scheme no longer provide any physical result (when taking $\mu_2/\mu_1 = 1$) due to the stiff nature of the system, and therefore no comparison can be done with the
<figure><img src="image_6.png" /><figcaption>Fig. 4: Membrane in simple shear flow: initial and final configurations for soft membranes.</figcaption></figure>

new semi-implicit method. Therefore, we take the viscosity ratio equal to 10 for the stiffest membrane to perform an adequate comparison, and 1 for the others. This also shows the ability of the method to deal easily with a viscosity ratio that may differ from one, as it may be essential for practical applications [19]. The initial conditions are given by

$$ (6.2) \qquad u_0(x, y) = \begin{bmatrix} \dot{\gamma}y \\ 0 \end{bmatrix}, \qquad \phi_0(x, y) = \sqrt{x^2 + y^2} - a, \qquad Y_0(x, y) = \begin{bmatrix} x \\ y \end{bmatrix}, $$

to impose the linear shear velocity field, a circular membrane, and no pre-deformation on the initial configuration. We set Neumann boundary conditions on the left and on the right of the domain, and moving wall at the top and at the bottom.

When a thin membrane is immersed into a shear flow, it starts elongating in the direction of the velocity vectors until it reaches a steady state position in which its shape no longer changes but its dynamics is still evolving. Indeed, in this configuration, the velocity vectors are perfectly tangent to the membrane and there is a constant stress applied on it, which makes it rotating around itself: this phenomenon is called *tank-treading* [28]. An overview of the initial and final position of the membrane for this experiment is shown in Figure 6.1. To compare properly the final shapes of the membrane, the reference time to stop the simulation is set for all cases equal to 1.5, which was observed to be enough reach the configuration described above.

As mentioned before, we are interested in studying the dynamics of the membrane when the new semi-implicit scheme is used, and specifically how the additional viscosity influences the computational cost and the quality of the solution when compared to the original explicit scheme. In particular, we focus on the simulation of the tank-treading phenomenon by changing the stiffness coefficient of the structure, ranging from a soft membrane to a very stiff one, and performing a grid refinement study to analyze its influence.

In Table 6.1 we present the $\Delta t$ needed to perform several experiments with varying capillary number, from the softest to the stiffest, on a set of three nested uniform meshes 128 × 64, 256 × 128, and 512 × 256. It can be noticed that in all cases the semi-implicit method provides higher maximum time steps, proving also numerically speaking that the additional dissipation introduces consistent viscosity where it is
<table><caption>Table 1: Membrane in simple shear flow: comparisons between explicit (EX) and semi-implicit (SI) schemes in terms of maximum $\Delta t$ to reach a final time $t_f = 1.5$.</caption><thead><tr><th rowspan="2">N<sub>x</sub> × N<sub>y</sub></th><th colspan="2">Ca = 0.001</th><th colspan="2">Ca = 0.008</th><th colspan="2">Ca = 0.01</th><th colspan="2">Ca = 0.02</th></tr><tr><th>EX</th><th>SI</th><th>EX</th><th>SI</th><th>EX</th><th>SI</th><th>EX</th><th>SI</th></tr></thead><tbody><tr><td>128 × 64</td><td>6.0 10<sup>-3</sup></td><td>6.0 10<sup>-2</sup></td><td>2.5 10<sup>-2</sup></td><td>1.0 10<sup>-1</sup></td><td>3.0 10<sup>-2</sup></td><td>1.0 10<sup>-1</sup></td><td>5.0 10<sup>-2</sup></td><td>1.5 10<sup>-1</sup></td></tr><tr><td>256 × 128</td><td>1.0 10<sup>-3</sup></td><td>2.5 10<sup>-2</sup></td><td>1.0 10<sup>-2</sup></td><td>8.5 10<sup>-2</sup></td><td>1.5 10<sup>-2</sup></td><td>9.5 10<sup>-2</sup></td><td>2.5 10<sup>-2</sup></td><td>1.0 10<sup>-1</sup></td></tr><tr><td>512 × 256</td><td>2.0 10<sup>-4</sup></td><td>5.0 10<sup>-3</sup></td><td>6.0 10<sup>-3</sup></td><td>6.0 10<sup>-2</sup></td><td>5.0 10<sup>-3</sup></td><td>5.0 10<sup>-2</sup></td><td>1.0 10<sup>-2</sup></td><td>7.5 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>Table 2: Membrane in simple shear flow: time step ratio $\Delta t_{SI}/\Delta t_{EX}$ with data taken from Table 6.1.</caption><thead><tr><th rowspan="2">N<sub>x</sub> × N<sub>y</sub></th><th colspan="2">Ca = 0.001</th><th colspan="2">Ca = 0.008</th><th colspan="2">Ca = 0.01</th><th colspan="2">Ca = 0.02</th></tr><tr><th colspan="8">Δt<sub>SI</sub>/Δt<sub>EX</sub></th></tr></thead><tbody><tr><td>128 × 64</td><td>10</td><td></td><td>4</td><td></td><td>3.33</td><td></td><td>3</td><td></td></tr><tr><td>256 × 128</td><td>25</td><td></td><td>8.5</td><td></td><td>6.33</td><td></td><td>4</td><td></td></tr><tr><td>512 × 256</td><td>25</td><td></td><td>10</td><td></td><td>10</td><td></td><td>7.5</td><td></td></tr></tbody></table>

needed allowing us to use larger time steps. Moreover, the new method provides larger time steps when stiffer structures are under study, which is extremely promising especially for real applications. We believe that this numerical behavior can be explained by considering that when dealing with stiffer structure the elastic waves seem to become less relevant to capture the general dynamics of the problem. This entails larger time steps for stiff structures that do not deform that much, for which the CFL condition given by the elastic waves is very strict. Instead, since the deformations are more important for softer membranes, smaller time step ratios seem to be crucial to have good predictions.

The mesh refinement analysis also provides interesting results regarding the effect of the new method when refining the mesh. In particular, it is observed that a larger gain in terms of speedup is experienced on finer grids with respect to coarser grids, which is very promising also in view of the future development of the method for three-dimensional models. Table 6.2 shows the ratio between the semi-implicit time step $\Delta t_{SI}$ and the explicit time step $\Delta t_{EX}$, showing its observed increasing trend when using finer grids and stiffer membranes.

The increased dissipation experienced for softer structures can also be observed when studying the membrane final shape once the final time is reached. In the following, we discuss the effect of the numerical scheme on the physics of the problem, which is extremely important to understand if the maximum admissible time steps given in Table 6.1 are enough to provide a good approximated solution. For this reason, we perform a set of numerical simulations and compare the final shape obtained through the original explicit scheme with those simulated with the semi-implicit scheme. Moreover, to investigate further the impact of the time step on the diffusion we compare the explicit result with a set of three simulations, for each problem, for which the time step has been adapted ranging from the explicit upper bound to the semi-implicit one. For all cases, the semi-implicit simulations are run with larger time steps with respect to the explicit one. For instance, in Figure 6.2 where we show the numerical results
<figure><img src="image_7.png" /><figcaption>Fig. 5: Membrane in simple shear flow: final configuration on the finest mesh for Ca = 0.02 with both explicit (EX) and semi-implicit (SI) methods. For the solution computed with the explicit method, the maximum Δt shown in Table 6.1 is taken. For the semi-implicit method, a set of three simulations have been performed with Δt ranging from the explicit bound to semi-implicit one shown in Table 6.1.</figcaption></figure>

for the softest membrane Ca = 0.02 on the finest grid, we compare the zero level-set of the membrane obtained with the explicit method with $\Delta t = 1.0 \ 10^{-2}$ to the same curve given by the semi-implicit scheme run with $\Delta t = 2.5 \ 10^{-2}, 5.0 \ 10^{-2}, 7.5 \ 10^{-2}$. For all experiments we perform the same study and show the results on the finest grid. As expected since the proposed scheme is consistent we obtain that all shapes follow the trace of the explicit one. For the very soft membrane shown in Figure 6.2 we experience a greater difference in the local error with respect to the explicit solution which increases for larger time steps. It is observed that, when deformations are important, damping the elastic waves with the tensorial viscosity terms comes with less accurate solutions. Moreover, we conclude that for soft structures CFL restrictions are not that important, and therefore the use of more sophisticated methods is not needed.

Indeed, it is well-known that stricter restrictions occur for stiffer structures, and it is for those applications that an implicit or semi-implicit approach is crucial. In Figure 6.3, two intermediate test cases are considered with lower capillary numbers to study the effect of the new scheme in these configurations. For simplicity, the experiments performed with Ca = 0.01 and Ca = 0.008 will be discussed together due to similarities in the final shapes, but also for the maximum $\Delta t$ and gains shown in Tables 6.1 and 6.2. For both cases, the semi-implicit scheme introduces much less dissipation giving a final shape of the membrane that is much closer to the reference one with a higher speedup $\Delta t_{SI}/\Delta t_{EX}$, especially when compared to the results
<figure><img src="image_8.png" /><figcaption>Fig. 6: Membrane in simple shear flow: final configuration on the fine mesh for Ca = 0.01, 0.008 with both explicit (EX) and semi-implicit (SI) methods. For the solution computed with the explicit method, the maximum Δt shown in Table 6.1 is taken. For the semi-implicit method, a set of three simulations have been performed with Δt ranging from the explicit bound to semi-implicit one shown in Table 6.1.</figcaption></figure>

observed in Figure 6.2. In particular, in this case a speedup of 10 is observed on the finest mesh. As expected from previous results and assumptions, the best results are obtained for the stiffest membrane simulated with Ca = 0.001. Figure 6.4 shows that all semi-implicit results with higher time step are almost superimposed on the reference shape computed with Δt = 2.0 10⁻⁴. In particular, for this experiment, we observe a speedup of 10 on the coarse mesh and 25 on the fine one that translates into a dramatic reduction of the computational time which also comes with a very precise prediction of the membrane shape. Although elastic waves are extremely fast for this experiment, the additional viscosity seem to be able to dampen them and relax the time step restriction, still providing a good dynamics of the membrane with almost no additional costs.

We conclude by pointing out that the proposed approach behaves the best for the stiffest structure on the finest mesh. This is critical for both real applications that often concern stiff materials and computational efficiency, which one would like to improve especially when refining the mesh.

**7. Summary and future perspectives.** We presented a novel approach for the construction of semi-implicit schemes for Eulerian models for fluid-structure interaction. This approach is very interesting for the simulation of incompressible flows in which the speed of elastic waves is much higher than that of the fluid causing the need of using extremely low time steps when compared to the simulation time. The focus of this work was on hyperelastic membranes and on the two-dimensional version of the model for full membrane elasticity introduced in [34]. The tensorial viscosity terms provide additional but consistent dissipation that activates only on the few cells where the membrane interface is spread. This was also demonstrated for a simplified
<figure><img src="image_9.png" /><figcaption>Fig. 7: Membrane in simple shear flow: final configuration on the fine mesh for Ca = 0.001 with both explicit (EX) and semi-implicit (SI) methods. For the solution computed with the explicit method, the maximum Δt shown in Table 6.1 is taken. For the semi-implicit method, a set of three simulations have been performed with Δt ranging from the explicit bound to semi-implicit one shown in Table 6.1.</figcaption></figure>

one-dimensional model in Section 3, where this additional term allowed us to prove
that the scheme is unconditionally stable. As shown in Section 4, in the full nonlinear
fluid-structure system these consistent terms appear in the momentum equation (4.9)
and depend on the time step, the stiffness parameter, and the local deformation of the
membrane. In Section 6 a set of numerical experiments is performed for a membrane
immersed in a simple shear flow. It was shown that the dependences mentioned before
influence the performances of the new scheme, which is observed to be more efficient
and accurate for stiffer structures and finer grids. In particular, the numerical method
proposed herein performed best on the experiment with the lowest capillary number
showing a remarkable speedup of 10 on the coarse mesh and 25 on the fine one while
preserving a very good agreement with the reference solution.

The perspectives of this work are several and range from the development of the
method for other hyperelastic models to its application to advance simulations of
biological systems.

Regarding the investigations on the method, in this work we considered the two-
dimensional version of the full membrane model introduced in [34], and future studies
will also focus on the development of the method for the full three-dimensional model.
Following the same reasoning, new evolution equations will be written to take into
account the modeling of shear variation. Therefore, this will introduce viscosity terms
that depend on both the area and shear variation. Other models will also be considered
to study the effect of this approach on surface tension (special case of the membrane
model with area variation) and, anisotropic elastic bodies [9].

About the possible applications, the proposed method is very promising for the simulation of biological capsules and red blood cells (RBCs), modelled as thin membranes characterized by shear effects, and their interaction with complex geometries [24, 23].

## Appendix A. Z is the local area variation.

In this appendix, we prove that Z introduced in (2.7) measures the local area variation as it was firstly introduced in [7]. In particular, in [7] it was proven that for incompressible media the gradient of the level-set $|∇φ|$ transports information about the local area variation. Afterwards, in [34] a new theory was introduced with the goal of having a single formulation (for both compressible and incompressible media) to deal with both area and shear variations. For this reason, two parameters have been introduced for three-dimensional membranes as functions only of the backward characteristics Y.

Here, we prove that to capture the area variation for two-dimensional membranes a different parameter is needed with respect to [34], and we also show how to relate it to the theory introduced in [7].

**PROPOSITION A.1.** Given that $J = \det(∇_ξX) = \det(∇Y)^{-1}$ is the volume ratio, different than one for general compressible media,

$$ (A.1) \qquad Z = J \frac{|∇φ|}{|∇φ_0(Y)|}. $$

*Proof.* The local area variation defined in (2.7) is measured by

$$ (A.2) \qquad Z = \sqrt{\mathrm{Tr}(\mathcal{A})}. $$

When using definition (2.6), we can write

$$ (A.3) \qquad \mathrm{Tr}(\mathcal{A}) = \mathrm{Tr}(B) - \frac{\mathrm{Tr}((Bn) \otimes (Bn))}{(Bn) \cdot n} = \frac{((\mathrm{Tr}(B)B - B^2)n) \cdot n}{(Bn) \cdot n} = \frac{\det(B)}{(Bn) \cdot n}, $$

where, for $B \in \mathcal{M}_2(\mathbb{R})$, the Cayley-Hamilton theorem $B^2 - \mathrm{Tr}(B)B + \det(B)I = 0$ was used in the last step.

With Equation (2.5), and $∇φ = [∇Y]^T∇φ_0(Y)$, we can write

$$ (A.4) \qquad (Bn) \cdot n = \left( [∇Y]^{-1} [∇Y]^{-T} \frac{∇φ}{|∇φ|} \right) \cdot \frac{∇φ}{|∇φ|} = \left| [∇Y]^{-T} \frac{∇φ}{|∇φ|} \right|^2 = \frac{|∇φ_0(Y)|^2}{|∇φ|^2}. $$

With the relation $\det(B) = J^2$ and (A.4), Equation (A.2) reduces to (A.1). □

**Remark A.2 (Volume ratio).** It should be noticed that Equation (A.2) measures the area variation for both compressible and incompressible flows. In particular, for incompressible flows, J is exactly equal to one and $Z = |∇φ|/|∇φ_0(Y)|$ as it was introduced in [7].

## Appendix B. Proof of Proposition 4.1.

*Proof.* By performing the derivative in time of $σ$ we obtain

$$ (B.1) \qquad \partial_t \sigma = \partial_t f(Z) \mathcal{C} + f(Z) \partial_t \mathcal{C}. $$
We first focus on the equation related to $f(Z)$. Considering the time derivative of $f(Z)$ gives us

$$ (B.2) \qquad \partial_t f(Z) + u \cdot \nabla f(Z) = f'(Z)Z[\nabla u] : \mathcal{C}, $$

where the last step is obtained with the equation on $Z$ taken from [34], recalled here

$$ (B.3) \qquad \partial_t Z + u \cdot \nabla Z = Z[\nabla u] : \mathcal{C}. $$

To compute the equation on $\mathcal{C}$, an equation on $n$ is needed. The evolution of $n$ is computed by deriving Equation (2.5) as follows

$$ \partial_t n = \frac{1}{|\nabla \phi|} \left[ \partial_t (\nabla \phi) - \left( \partial_t (\nabla \phi) \cdot \frac{\nabla \phi}{|\nabla \phi|} \right) \frac{\nabla \phi}{|\nabla \phi|} \right]. $$

This relation together with the evolution equation on the gradient of $\phi$ (2.6)

$$ \partial_t (\nabla \phi) + u \cdot \nabla (\nabla \phi) = -[\nabla u]^T \nabla \phi, $$

gives us

$$ (B.4) \qquad \partial_t n + u \cdot \nabla n = -[\nabla u]^T n + (([\nabla u]n) \cdot n)n. $$

With the identities $(Ab) \otimes c = A(b \otimes c)$ and $b \otimes (Ac) = (b \otimes c)A^T$, it is possible to write the equation on $\mathcal{C}$, which reads

$$ (B.5) \qquad \partial_t \mathcal{C} + u \cdot \nabla \mathcal{C} = [\nabla u]^T (n \otimes n) + (n \otimes n)[\nabla u] - 2(([\nabla u]n) \cdot n)(n \otimes n). $$

Using (B.5) and (B.2) in (B.1) we get the evolution equation (4.6) stated in the proposition. □

## Appendix C. Discretization of semi-implicit operators.

In this appendix, we briefly describe how we discretize the semi-implicit operators for the staggered mesh framework used by the incompressible Navier-Stokes solver. In particular, for a two-dimensional problem, up to 16 terms need to be discretized: $\partial_t(m \partial_j u_k)$, with $i, j, k = x, y$, for both $u_x$ and $u_y$ located on vertical and horizontal faces, respectively (see Figure 4.1). Indeed, since the position of the $x-$ and $y-$components of a vector field are located in different points the discretization of the same operator on one direction is going to be different to that on the other direction. For compactness, here we do not present the discretization of all terms but, we mainly focus on few of them to give an idea on how to discretize these operators. There are of course several ways to discretize such terms but as long as they are consistent no remarkable difference is observed. For simplicity a uniform discretization is considered, meaning that each cell has a $\Delta x \times \Delta y$ area. To simplify the notation, only in this appendix we take the following abuse of notation: $u = u_x$ and $v = u_y$.

For both components, we discretize two trivial operators and a third trickier one. For the $x$-direction (on vertical faces) we consider the discretization of $\partial_x(m \partial_x u)$, $\partial_y(m \partial_y u)$, and $\partial_y(m \partial_y v)$:

$$ \partial_x(m \partial_x u)|_{i-1/2,j} = \frac{1}{\Delta x} \left( m_{i,j} \frac{u_{i+1/2,j} - u_{i-1/2,j}}{\Delta x} - m_{i-1,j} \frac{u_{i-1/2,j} - u_{i-3/2,j}}{\Delta x} \right), $$

$$ \partial_y(m \partial_y u)|_{i-1/2,j} = \frac{1}{\Delta y} \left( m^*_{i-1/2,j+1/2} \frac{u_{i-1/2,j+1} - u_{i-1/2,j}}{\Delta y} - m^*_{i-1/2,j-1/2} \frac{u_{i-1/2,j} - u_{i-1/2,j-1}}{\Delta y} \right), $$

$$ \partial_y(m \partial_y v)|_{i-1/2,j} = \frac{1}{\Delta y} \left( m^*_{i-1/2,j+1/2} \frac{v^*_{i-1/2,j+1} - v^*_{i-1/2,j}}{\Delta y} - m^*_{i-1/2,j-1/2} \frac{v^*_{i-1/2,j} - v^*_{i-1/2,j-1}}{\Delta y} \right), \blacksquare $$
where $(\cdot)^*$ identifies an interpolated variable. In this case, we can approximate the multiplicative factor $m$ as

$$m_{i-1/2, j+1/2}^* = \frac{1}{4}(m_{i,j} + m_{i,j+1} + m_{i-1,j+1} + m_{i-1,j}),$$

$$m_{i-1/2, j-1/2}^* = \frac{1}{4}(m_{i,j-1} + m_{i,j} + m_{i-1,j} + m_{i-1,j-1}),$$

and the velocity component $v$ as

$$v_{i-1/2, j+1}^* = \frac{1}{4}(v_{i,j+1/2} + v_{i,j+3/2} + v_{i-1,j+3/2} + v_{i-1,j+1/2}),$$

$$v_{i-1/2, j}^* = \frac{1}{4}(v_{i,j-1/2} + v_{i,j+1/2} + v_{i-1,j+1/2} + v_{i-1,j-1/2}),$$

$$v_{i-1/2, j-1}^* = \frac{1}{4}(v_{i,j-3/2} + v_{i,j-1/2} + v_{i-1,j-1/2} + v_{i-1,j-3/2}).$$

For the $y$-direction (on horizontal faces) we consider the discretization of $\partial_y(m \partial_y v)$, $\partial_x(m \partial_x v)$, and $\partial_x(m \partial_x u)$:

$$\partial_y(m \partial_y v)|_{i,j-1/2} = \frac{1}{\Delta y} \left( m_{i,j} \frac{v_{i,j+1/2} - v_{i,j-1/2}}{\Delta y} - m_{i,j-1} \frac{v_{i,j-1/2} - v_{i,j+3/2}}{\Delta y} \right),$$

$$\partial_x(m \partial_x v)|_{i,j-1/2} = \frac{1}{\Delta x} \left( m_{i+1/2,j-1/2}^* \frac{v_{i+1,j-1/2} - v_{i,j-1/2}}{\Delta x} - m_{i-1/2,j-1/2}^* \frac{v_{i,j-1/2} - v_{i-1,j-1/2}}{\Delta x} \right),$$

$$\partial_x(m \partial_x u)|_{i,j-1/2} = \frac{1}{\Delta x} \left( m_{i+1/2,j-1/2}^* \frac{u_{i+1,j-1/2}^* - u_{i,j-1/2}^*}{\Delta x} - m_{i-1/2,j-1/2}^* \frac{u_{i,j-1/2}^* - u_{i-1,j-1/2}^*}{\Delta x} \right).$$

where we can approximate the remaining $m$ term as

$$m_{i+1/2, j-1/2}^* = \frac{1}{4}(m_{i+1,j-1} + m_{i+1,j} + m_{i,j} + m_{i,j-1}),$$

and the velocity component $u$ as

$$u_{i+1,j-1/2}^* = \frac{1}{4}(u_{i+3/2,j-1} + u_{i+3/2,j} + u_{i+1/2,j} + u_{i+1/2,j-1}),$$

$$u_{i,j-1/2}^* = \frac{1}{4}(u_{i+1/2,j-1} + u_{i+1/2,j} + u_{i-1/2,j} + u_{i-1/2,j-1}),$$

$$u_{i-1,j-1/2}^* = \frac{1}{4}(u_{i-1/2,j-1} + u_{i-1/2,j} + u_{i-3/2,j} + u_{i-3/2,j-1}).$$

**Acknowledgments.** MC would like to acknowledge the support of Antoine Lemoine for the development of this project within the open-source massively parallel software notus CFD (https://notus-cfd.org/).

REFERENCES

[1] E. ABBATE, A. IOLLO, AND G. PUPPO, An all-speed relaxation scheme for gases and compressible materials, Journal of Computational Physics, 351 (2017), pp. 1–24.

[2] E. ABBATE, A. IOLLO, AND G. PUPPO, An asymptotic-preserving all-speed scheme for fluid dynamics and nonlinear elasticity, SIAM Journal on Scientific Computing, 41 (2019), pp. A2850–A2879.

[3] T.-D. ASLAM, A partial differential equation approach to multidimensional extrapolation, Journal of Computational Physics, 193 (2003), pp. 349–355.
[4] M. BERGMANN, A. FONDANÈCHE, AND A. IOLLO, *An eulerian finite-volume approach of fluid-structure interaction problems on quadtree meshes*, Journal of Computational Physics, 471 (2022), p. 111647.

[5] C. BOST, G.-H. COTTET, AND E. MAITRE, *Linear stability analysis of a level set model of immersed elastic membrane*, (2009).

[6] J. U. BRACKBILL, D. B. KOTHE, AND C. ZEMACH, *A continuum method for modeling surface tension*, Journal of computational physics, 100 (1992), pp. 335-354.

[7] G.-H. COTTET AND E. MAITRE, *A level set method for fluid-structure interactions with immersed surfaces*, Mathematical models and methods in applied sciences, 16 (2006), pp. 415-438.

[8] G.-H. COTTET AND E. MAITRE, *A semi-implicit level set method for multiphase flows and fluid-structure interaction problems*, Journal of Computational Physics, 314 (2016), pp. 80-92.

[9] G.-H. COTTET, E. MAITRE, AND T. MILCENT, *Eulerian formulation and level set models for incompressible fluid-structure interaction*, ESAIM: Mathematical Modelling and Numerical Analysis, 42 (2008), pp. 471-492.

[10] A. DE BRAUER, A. IOLLO, AND T. MILCENT, *A cartesian scheme for compressible multimaterial models in 3d*, Journal of Computational Physics, 313 (2016), pp. 121-143.

[11] A. DE BRAUER, A. IOLLO, AND T. MILCENT, *A cartesian scheme for compressible multimaterial hyperelastic models with plasticity*, Communications in Computational Physics, 22 (2017), pp. 1362-1384.

[12] J. DEBORDE, T. MILCENT, P. LUBIN, AND S. GLOCKNER, *Numerical simulations of the interaction of solitary waves and elastic structures with a fully eulerian method*, Water Waves, 2 (2020), pp. 433-466.

[13] J. DONEA, S. GIULIANI, AND J.-P. HALLEUX, *An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions*, Computer methods in applied mechanics and engineering, 33 (1982), pp. 689-723.

[14] R. FALGOUT, J. JONES, AND U. YANG, *The design and implementation of hypre, a library of parallel high performance preconditioners*, chapter in Numerical Solution of Partial Differential Equations on Parallel Computers, A.M. Bruaset and A. Tveito, eds., Springer-Verlag, UCRL-JRNL-205459, 51 (2006), pp. 267-294.

[15] R. FALGOUT AND U. YANG, *Hypre: a library of high performance preconditioners*, Computational Science - ICCS 2002 Part III, P.M.A. Sloot, C.J.K. Tan. J.J. Dongarra, and A.G. Hoekstra, eds., vol. 2331 of Lecture Notes in Computer Science, Springer-Verlag, UCRL-JC-146175, (2002), pp. 632-641.

[16] T. FANION, M. FERNÁNDEZ, AND P. LE TALLEC, *Deriving adequate formulations for fluid-structure interaction problems: from ale to transpiration*, Revue européenne des éléments finis, 9 (2000), pp. 681-708.

[17] M. A. FERNÁNDEZ, J.-F. GERBEAU, AND C. GRANDMONT, *A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid*, International Journal for Numerical Methods in Engineering, 69 (2007), pp. 794-821.

[18] M. A. FERNÁNDEZ, J. MULLAERT, AND M. VIDRASCU, *Explicit robin-neumann schemes for the coupling of incompressible fluids with thin-walled structures*, Computer Methods in Applied Mechanics and Engineering, 267 (2013), pp. 566-593.

[19] E. FOESSEL, J. WALTER, A.-V. SALSAC, AND D. BARTHÈS-BIESEL, *Influence of internal viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow*, Journal of Fluid Mechanics, 672 (2011), pp. 477-486.

[20] A. FONDANÈCHE, *Interaction fluide-structure dans un dispositif vasculaire actif*, PhD thesis, Université de Bordeaux, 2021.

[21] Y. GORSSE, A. IOLLO, T. MILCENT, AND H. TELIB, *A simple cartesian scheme for compressible multimaterials*, Journal of Computational Physics, 272 (2014), pp. 772-798.

[22] J.-L. GUERMOND, P. MINEV, AND J. SHEN, *An overview of projection methods for incompressible flows*, Computer Methods in Applied Mechanics and Engineering, 195 (2006), p. 6011-6045.

[23] E. HÄNER, D. VESPERINI, A.-V. SALSAC, A. LE GOFF, AND A. JUEL, *Sorting of capsules according to their stiffness: from principle to application*, Soft Matter, 17 (2021), pp. 3722-3732.

[24] X.-Q. HU, B. SÉVÉNIÉ, A.-V. SALSAC, E. LECLERC, AND D. BARTHÈS-BIESEL, *Characterizing the membrane properties of capsules flowing in a square-section microfluidic channel: Effects of the membrane constitutive law*, Phys. Rev. E, 87 (2013), p. 063008, https://doi.org/10.1103/PhysRevE.87.063008, https://link.aps.org/doi/10.1103/PhysRevE.87.063008.

[25] S. HYSING, *A new implicit surface tension implementation for interfacial flows*, International Journal for Numerical Methods in Fluids, 51 (2006), pp. 659-672.
[26] S. II, K. SUGIYAMA, S. TAKEUCHI, S. TAKAGI, AND Y. MATSUMOTO, *An implicit full eulerian method for the fluid-structure interaction problem*, International Journal for Numerical Methods in Fluids, 65 (2011), pp. 150–165.

[27] G. JIANG AND C. SHU, *Efficient implementation of weighted eno schemes*, Journal of Computational Physics, (1996).

[28] S. R. KELLER AND R. SKALAK, *Motion of a tank-treading ellipsoidal particle in a shear flow*, Journal of Fluid Mechanics, 120 (1982), pp. 27–47.

[29] L. LEE AND R. J. LEVEQUE, *An immersed interface method for incompressible navier–stokes equations*, SIAM Journal on Scientific Computing, 25 (2003), pp. 832–856.

[30] A. LEGAY, J. CHESSA, AND T. BELYTSCHKO, *An eulerian–lagrangian method for fluid–structure interaction based on level sets*, Computer Methods in Applied Mechanics and Engineering, 195 (2006), pp. 2070–2087.

[31] H. LIU AND K. KAWACHI, *A numerical study of insect flight*, Journal of computational physics, 146 (1998), pp. 124–156.

[32] E. MAITRE, T. MILCENT, G.-H. COTTET, A. RAOULT, AND Y. USSON, *Applications of level set methods in computational biophysics*, Mathematical and Computer Modelling, 49 (2009), pp. 2161–2169.

[33] A. A. MAYO AND C. S. PESKIN, *An implicit numerical method for fluid dynamics problems with immersed elastic boundaries*, Contemporary Mathematics, 141 (1992), pp. 261–261.

[34] T. MILCENT AND E. MAITRE, *Eulerian model of immersed elastic surfaces with full membrane elasticity*, Communications in Mathematical Sciences, 14 (2016), pp. 857–881, https://doi.org/https://dx.doi.org/10.4310/CMS.2016.v14.n3.a11.

[35] R. MITTAL AND G. IACCARINO, *Immersed boundary methods*, Annu. Rev. Fluid Mech., 37 (2005), pp. 239–261.

[36] M. SUSSMAN, P. SMEREKA, AND S. OSHER, *A level set approach for computing solutions to incompressible two-phase flow*, Journal of Computational Physics, 114 (1994), pp. 146–159.

[37] E. P. NEWREN, *Enhancing the immersed boundary method: stability, volume conservation, and implicit solvers*, vol. 68, Citeseer, 2007.

[38] S. OSHER AND J. A. SETHIAN, *Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations*, Journal of Computational Physics, 79 (1988), pp. 12–49, https://doi.org/10.1016/0021-9991(88)90002-2, https://www.sciencedirect.com/science/article/pii/0021999188900022.

[39] C. S. PESKIN, *Flow patterns around heart valves: a numerical method*, Journal of computational physics, 10 (1972), pp. 252–271.

[40] C. S. PESKIN, *The immersed boundary method*, Acta numerica, 11 (2002), pp. 479–517.

[41] T. RICHTER, *A fully eulerian formulation for fluid-structure-interaction problems*, Journal of Computational Physics, 233 (2013), pp. 227–240.

[42] M. SAHIN AND K. MOHSENI, *An arbitrary lagrangian–eulerian formulation for the numerical simulation of flow patterns generated by the hydromedusa aequorea victoria*, Journal of Computational Physics, 228 (2009), pp. 4588–4605.

[43] J. M. STOCKIE AND B. R. WETTON, *Analysis of stiffness in the immersed boundary method and implications for time-stepping schemes*, Journal of Computational Physics, 154 (1999), pp. 41–64.

[44] K. SUGIYAMA, S. II, S. TAKEUCHI, S. TAKAGI, AND Y. MATSUMOTO, *A full eulerian finite difference approach for solving fluid-structure coupling problems*, Journal of Computational Physics, 230 (2011), pp. 596–627.

[45] M. SUSSMAN AND M. OHTA, *A stable and efficient method for treating surface tension in incompressible two-phase flow*, SIAM Journal on Scientific Computing, 31 (2009), pp. 2447–2471.

[46] A. THOMANN, A. IOLLO, AND G. PUPPO, *Implicit relaxed all mach number schemes for gases and compressible materials*, SIAM Journal on Scientific Computing, 45 (2023), pp. A2632–A2656.

[47] C. TU AND C. S. PESKIN, *Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods*, SIAM Journal on Scientific and Statistical Computing, 13 (1992), pp. 1361–1376.