# Fully Eulerian models for the numerical simulation of capsules with an elastic bulk nucleus

Florian Desmons, Thomas Milcent, Anne-Virginie Salsac, Mirco Ciallella

► To cite this version:

Florian Desmons, Thomas Milcent, Anne-Virginie Salsac, Mirco Ciallella. Fully Eulerian models for the numerical simulation of capsules with an elastic bulk nucleus. 2023. hal-04272881

HAL Id: hal-04272881

https://hal.science/hal-04272881v1

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Fully Eulerian models for the numerical simulation of
capsules with an elastic bulk nucleus

Florian Desmons<sup>a</sup>, Thomas Milcent<sup>a</sup>, Anne-Virginie Salsac<sup>b</sup>, Mirco Ciallella<sup>a,\*</sup>

<sup>a</sup>École Nationale Supérieure d'Arts et Métiers, Institut de Mécanique et d'Ingénierie, 33400 Talence, France

<sup>b</sup>Biomechanics and Bioengineering Laboratory, Université de Technologie de Compiègne, 60203 Compiègne, France

**Abstract**

In this paper, we present a computational framework based on fully Eulerian models for fluid-
structure interaction for the numerical simulation of biological capsules. The flexibility of such
models, given by the Eulerian treatment of the interface and deformations, allows us to easily deal
with the large deformations experienced by the capsule. The modeling of the membrane is based
on the full membrane elasticity model introduced in (Milcent, T., Maitre, E. (2016)) that is capable
of capturing both area and shear variations thanks to the so-called backward characteristics. In
the validation section several test cases are presented with the goal of comparing our results to
others present in the literature. In this part, the comparisons are done with different well-known
configurations (capsule in shear flow and square-section channel), and by deepening the effect of
the elastic constitutive law and capillary number on the membrane dynamics. Finally, to show the
potential of this framework we introduce a new test case that describes the relaxation of a capsule
in an opening channel. In order to increase the challenges of this test we study the influence of an
internal nucleus, modeled as a hyperelastic solid, on the membrane evolution. Several numerical
simulations are presented to deeply study its influence by modifying the characteristic parameters
of the nucleus (size and elastic parameter).

***Keywords:*** Biological capsules, Fluid-structure interaction, Eulerian elasticity, full membrane elasticity, Incompressible Navier-Stokes

**1. Introduction**

Biological capsules have revolutionized the fields of medicine and biotechnology, offering re-
markable advancements and opportunities for research, development, and practical applications.
For example, in medicine, the use of biological capsules have greatly improved drug delivery sys-
tems [9]. By encapsulating drugs within protective coatings, scientists can enhance their stability,
solubility, and targeted release. Another field of application in biotechnology concerns the en-
capsulation for the efficient confinement of enzymes within protective matrices, enhancing their
reusability, and enabling their application in diverse industrial processes [27].

In the context of studying the interaction between such capsules and a surrounding fluid, nu-
merical simulation of fluid-structure interaction (FSI) is a very promising tool to understand their
behavior and deformation. FSI simulations involve the coupled analysis of fluid flow and struc-
tural deformation, allowing researchers to investigate the complex interactions between structures
and their surrounding fluid environments. When applied to biological capsules, FSI simulations

<footnote>*Corresponding author (mirco.ciallella@ensam.eu)</footnote>
provide valuable insights into their mechanical behavior, fluid dynamics, and overall performance. Nowadays, several methods are used to tackle the coupled problem coming from the interaction between fluid and structure. In the seminal work [35], Peskin introduced an innovative approach where Lagrangian markers are used to track the membrane deformation. The elastic force appears as a source term in the fluid equations and is spread on the fluid grid with discretized Dirac mass. This work was established as a reference for future developments of the so-called immersed boundary method [32, 26, 36, 19, 42]. Other popular approaches to solve such problems are based on arbitrary Lagrangian–Eulerian (ALE) methods [29, 41, 11, 14, 15], which consider body-fitted grids to follow the displacement of the interface. Although ALE method has been widely used in many FSI problems, they become extremely cumbersome to apply to problems with large deformations, due to the challenging re-meshing process.

More recently, Eulerian models for fluid-structure interaction have been derived for computational purposes as they enable complex multi-dimensional applications, with large deformation, to be simulated quite easily on fixed Cartesian meshes [4, 5, 6]. In this framework, both the fluid and the solid are described by an Eulerian approach. The interface and deformations of the media are captured by Eulerian fields, which are advected by the fluid velocity in a level-set fashion [34]. The fluid-structure problem is recast as a complex flow: the Navier-Stokes equations with an elastic source term is coupled with a transport equation on the Eulerian interface and deformation. In the last decade, this approach has been applied to both incompressible and compressible materials [30, 3, 18, 7] proving its high flexibility, especially when dealing with complex geometries.

When dealing with microcapsules, Stokes flows are sometimes considered due to the low Reynolds numbers, which corresponds to slow and highly viscous flows. Furthermore, the Stokes equations are simpler than the Navier-Stokes equations and can be solved more efficiently using the boundary integral method [39, 38, 43]. The boundary integral method directly solves the fluid flow problem on the surface of the capsule, eliminating the need for a volumetric discretization of the entire fluid domain. This reduction in dimensionality significantly reduces the computational cost. Difficulties remain for the simulation of large deformations, interactions with complex geometries and applications where the viscosity of the inner fluid is different from that of the outer fluid [16].

The novelty of this paper is to present a framework based on Eulerian models for FSI for the realistic simulation of capsules undergoing large deformations. In particular, the objective is to study the influence of an internal nucleus, modeled as an elastic solid [8], within the capsule for long-time simulations. The introduction of an internal nucleus is extremely promising for numerous applications since many biological cells, e.g. white blood cells, are often modeled as a thin membrane with a nucleus, also modeled as a thin elastic membrane [1]. Once the proposed model and method are validated against the literature of microcapsules, a new configuration of a relaxation phenomenon [17] is simulated. Due to the difficulties in simulating both experimentally and numerically this phenomenon, up to our knowledge this is the first simulation of this kind, where the nuclues is modeled as an elastic bulk, and it was possible thanks to the flexibility of the fully Eulerian framework.

The paper is organized as follows. First, in Section 2 we describe the physical problem, represented by a capsule immersed within an incompressible fluid, and its main parameters. Next, in Section 3, we briefly recall the formulation of the Eulerian model used to simulated this complex fluid-structure problem. In Section 4, we present the general numerical scheme used to discretize the FSI model along with a novel method to mitigate the degradation of the deformation vector in the case of large deformation and long-time simulations. The mathematical model, along with the proposed numerical method, is then validated with respect to other approaches present in the literature, in Section 5. In Section 6, the relaxation phenomenon is deeply analyzed considering the
<figure><img src="image_2.png" /><figcaption>Figure 1: Characteristic variables of the problem.</figcaption></figure>

impact of an internal solid, and its stiffness, on the capsule elastic behavior. Finally, in Appendix A, we present a test case that is a modification of the sheared elastic sphere relaxation proposed in [31].

## 2. Problem statement

We consider an elastic capsule (see Figure 1) modeled as a thin elastic membrane, with an internal nucleus of radius $a_b$ and stiffness $\chi$, immersed within an incompressible Newtonian fluid with viscosity $\mu$ and density $\rho$. The same viscosity and density are considered inside and outside the membrane. The radius $a_m$ of the capsule membrane is taken as a characteristic length of the problem. The capsule immersed in the fluid undergoes large deformations and its elastic properties are characterized by a non linear constitutive law depending on the shear modulus $G_s$ and the dilatation modulus $K_s$. The fluid and the elastic capsule are strongly coupled and thus this phenomenon is a fluid-structure interaction problem. Therefore, the phenomena arising are governed by two dimensionless numbers: the Reynolds number

$$Re = \frac{\rho U a_m}{\mu}, \qquad (1)$$

and the capillary number

$$Ca = \frac{\mu U}{G_s} \qquad (2)$$

where $U$ is the fluid reference velocity.

The additional dimensionless parameters $a_b/a_m$ and $\chi/G_s$ are considered to study the effect of the internal nucleus on the capsule profile and deformation. In this work, we propose the original idea of using a fully Eulerian model for the simulation of capsules. The elastic deformations of the capsule are described by Eulerian fields that are advected by the fluid velocity and the associated elastic force appears as an Eulerian source term in the fluid equations.

## 3. Fully Eulerian fluid-structure model

### 3.1. Eulerian description of the fluid

Herein, the incompressible fluid modeled through the incompressible Navier-Stokes equations:

$$ \left\{ \begin{array}{l} \rho(\partial_t u + (u \cdot \nabla)u) - \text{div}(2\mu D(u)) + \nabla p = F \\ \text{div}(u) = 0 \end{array} \right. \qquad (3) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Forward and backward characteristics.</figcaption></figure>

where $D(u) = \frac{1}{2}([\nabla u] + [\nabla u]^T)$ is the strain rate tensor, $p$ the pressure and $F$ is a source term. When no structure is considered, the source term $F$ is set to zero retrieving the original Navier-Stokes equations. In later sections, more details are given on the formulation of the elastic force term for both thin membranes and elastic solids.

### 3.2. Eulerian description of solid deformations

Let $\Omega_0 \subset \mathbb{R}^3$ be the reference configuration of a continuous medium and $\Omega_t \subset \mathbb{R}^3$ the deformed configuration at time $t$. In order to describe the evolution of this medium in the Lagrangian frame we define the forward characteristics $X(t, \xi)$ as the image at time $t$ in the deformed configuration $\Omega_t = X(t, \Omega_0)$ of a material point $\xi$ belonging to the initial configuration, i.e., $X : \mathbb{R}^3 \times \mathbb{R}^3 \to \mathbb{R}^3$ (see Fig. 2). The corresponding Eulerian velocity field is defined as $u : \mathbb{R}^3 \times \mathbb{R}^+ \to \mathbb{R}^3$ where

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

To describe the continuous medium in the Eulerian frame, we introduce the backward characteristics $Y(x, t)$ that for a time $t$ and a point $x$ in the deformed configuration, gives the corresponding initial point $\xi$ in the initial configuration, i.e., $Y : \mathbb{R}^3 \times \mathbb{R}^+ \to \mathbb{R}^3$ (see Fig. 2). Since $Y(X(t, \xi), t) = \xi$, differentiating with respect to time and space in turn we have:

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

and

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

The relation (5) is the Eulerian equivalent of the characteristic equation (4). In addition, equation (6) allows to compute the gradient of the deformation in the Eulerian frame via $Y$. The next sections are devoted to the description of bulk and membrane elastic deformations with the backward characteristics $Y$.

### 3.3. Hyperelastic bulk models

For hyperelastic bulk material, the internal energy $E_b$ is a function of the deformation tensor $\nabla_{\xi} X$. We focus in this paper on materials that are Galilean invariant and isotropic. With these
assumptions it can be proven (see [21]) that $\mathcal{E}_b$ is expressed as a function of the invariants of the left Cauchy-Green tensor $B = [\nabla_\xi X][\nabla_\xi X]^T$. This tensor is written in the Eulerian description with (6)

$$B(x,t) = [\nabla_x Y]^{-1} [\nabla_x Y]^{-T}. \qquad (7)$$

In this work, we consider the neo-Hookean constitutive law ($\mathcal{E}_b$ depends only on $\mathrm{Tr}(B)$). The associated Cauchy stress tensor is given by $\sigma_b = 2\chi B$ where $\chi$ is the elastic modulus coefficient. In this work, we aim at developing a numerical tool able to simulate elastic capsules with a bulk kernel inside. The neo-Hookean model allows to describe, with one parameter, elastic media subjected to large but moderate deformations. To model various aspects of the resulting non linear relation between stress and deformation, other elastic constitutive laws can be considered: Mooney-Rivlin, Saint Venant Kirchhoff, Ogden. In the final model the bulk solid will be immersed in a incompressible fluid and the fluid-structure interface will be captured by a level set function $\phi_b : \mathbb{R}^3 \times \mathbb{R}^+ \to \mathbb{R}$ which is advected by the Eulerian velocity field $u$:

$$\partial_t \phi_b + u \cdot \nabla \phi_b = 0. \qquad (8)$$

The force associated to the neo-Hookean solid is then given by

$$F_b = \mathrm{div} \left( H \left( \frac{\phi_b}{\varepsilon} \right) \sigma_b \right), \quad \sigma_b = 2\chi B, \qquad (9)$$

where $H$ is a smoothed Heaviside function and $2\varepsilon$ is the width of the interface. Therefore $\phi_b > \varepsilon$ corresponds to the fluid domain and $\phi_b < -\varepsilon$ to the solid domain.

### 3.4. Hyperelastic membrane models

The notations and results summarized in this section are detailed in [31]. We consider a surface $\Gamma_t = \{x \in \mathbb{R}^3 / \phi_m(x,t) = 0\}$ captured by a level set function $\phi_m : \mathbb{R}^3 \times \mathbb{R}^+ \to \mathbb{R}$ and advected by the Eulerian velocity field $u$:

$$\partial_t \phi_m + u \cdot \nabla \phi_m = 0. \qquad (10)$$

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

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

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

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

where $B$ is given by (7). This tensor measures surface deformations by projecting the 3D 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$. We denote by $(\lambda_1)^2$ and $(\lambda_2)^2$ the two other eigenvalues. The other invariants are used to define the following quantities

$$Z_1 = \sqrt{\mathrm{Tr}(\mathrm{Cof}(\mathcal{A}))} = |\lambda_1 \lambda_2|, \qquad (13)$$

$$Z_2 = \frac{\mathrm{Tr}(\mathcal{A})}{2\sqrt{\mathrm{Tr}(\mathrm{Cof}(\mathcal{A}))}} = \frac{1}{2} \left( \left| \frac{\lambda_1}{\lambda_2} \right| + \left| \frac{\lambda_2}{\lambda_1} \right| \right). \qquad (14)$$
The eigenvalues $\lambda_1$ and $\lambda_2$ correspond to the local deformation of the surface, therefore it is intuitive that $Z_1$ measures the local area variation whereas $Z_2$ measures the local shear variation. These results are demonstrated in details in [31]. Note also that $Z_1 \ge 0$ whereas $Z_2 \ge 1$. Now we introduce a constitutive law E that depends on the invariants $Z_1$ and $Z_2$ and the associated energy which is localized in a neighborhood of the membrane,

$$ \mathcal{E}_m = \int_Q E(Z_1, Z_2) \frac{1}{\varepsilon} \zeta \left( \frac{\phi_m}{\varepsilon} \right) dx. \qquad (15) $$

Here $Q$ is a box containing the membrane, $\varepsilon$ the width of the interface and $\zeta$ is a cut-off function used to spread the interface near $\{\phi_m = 0\}$. Elastic energy is conservative so a variation of the membrane shape induces an elastic force (principle of virtual power) and is given by

$$ F_m = \operatorname{div} (E_{,1}(Z_1, Z_2) Z_1 \mathcal{C}_1 + E_{,2}(Z_1, Z_2) Z_2 \mathcal{C}_2) \frac{1}{\varepsilon} \zeta \left( \frac{\phi_m}{\varepsilon} \right), \qquad (16) $$

where

$$ \mathcal{C}_1 = \mathbb{I} - n \otimes n, \qquad \mathcal{C}_2 = \frac{2\mathcal{A}}{\operatorname{Tr}(\mathcal{A})} - (\mathbb{I} - n \otimes n), \qquad (17) $$

and $E_{,i}$ represents the derivative of E with respect to $Z_i$. Three constitutive laws to model membranes are used in this article:

* Evan-Skalak model:

$$ E_{,1}(r_1, r_2) = K_s(r_1 - 1), \qquad E_{,2}(r_1, r_2) = G_s, \qquad (18) $$

* Skalak model:

$$ E_{,1}(r_1, r_2) = \frac{G_s}{2}(-r_1^3 + 4r_1r_2^2 - r_1 - 2r_2) + \frac{K_s}{2}(r_1^3 - r_1), \qquad E_{,2}(r_1, r_2) = G_s r_1(2r_1r_2 - 1), \qquad (19) $$

* Membrane Neo-Hookean model:

$$ E_{,1}(r_1, r_2) = G_s \left( r_2 - \frac{1}{r_1^3} \right), \qquad E_{,2}(r_1, r_2) = G_s r_1. \qquad (20) $$

### 3.5. Overall model

The elastic membrane and bulk media are immersed in a incompressible fluid modeled by the Navier-Stokes equations. The overall fully Eulerian model is given by

$$ \left\{ \begin{array}{l} \rho(\partial_t u + (u \cdot \nabla)u) - \operatorname{div}(2\mu D(u)) + \nabla p = F_b(\phi_b, Y_b) + F_m(\phi_m, Y_m) \\ \partial_t Y_\ell + (u \cdot \nabla)Y_\ell = 0 \qquad \ell = b, m \\ \partial_t \phi_\ell + u \cdot \nabla \phi_\ell = 0 \qquad \ell = b, m \end{array} \right. \qquad (21) $$

The membrane force $F_m$ is given by (16) and the bulk force $F_b$ is given by (9). Note that a level set function $\phi_\ell$ and a backward characteristic field $Y_\ell$ is needed for each media (bulk (b) and membrane (m)). These equations are completed with appropriate initial and boundary conditions that will be detailed in the numerical results sections.
<figure><img src="image_4.png" /><figcaption>Figure 3: Staggered grid with position of unknowns</figcaption></figure>

# 4. Numerical Scheme

## 4.1. General scheme

The equations (21) are discretized with finite volume schemes on a staggered grid (see Figure 3 for a 2D configuration). In the following, the notation $\phi$ (respectively $Y$) denote either $\phi_b$ or $\phi_m$ (respectively $Y_b$ or $Y_m$) because the same discretization is used for each variable. 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 semi-discretization in time is given by the projection method (see [20] for an overview on projection methods) for the Navier Stokes equations (Steps 1-2-3) and an Euler explicit scheme for the advection equations (Step 4):

$$ \text{Step 1: } \rho^n \left( \frac{u^\star - u^n}{\Delta t} + \text{div}(u^n \otimes u^\star) \right) - \text{div}(2\mu^n D(u^\star)) + \nabla p^n = F(\phi^n, Y^n) $$

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

$$ \text{Step 3: } u^{n+1} = u^\star - \frac{\Delta t}{\rho^n} \nabla \psi^{n+1}, \qquad p^{n+1} = p^n + \psi^{n+1} $$

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

In Step 1 a prediction of the velocity $u^\star$ is computed with an Euler implicit scheme for the viscous term. The convection term is treated explicitly in time with NSSP3 Runge-Kutta [24] and in space with WENO5 [23], following an high-order momentum procedure [10], and the source term explicitly with a standard second-order discretization. In step 2, the Poisson equation for the pressure increment $\psi^{n+1}$ is solved with appropriate boundary conditions [20]. The resulting linear systems are solved with the GMRES algorithm of the HYPRE library [13, 12] with Jacobi preconditioning for Step 1 and with multigrid preconditioning for Step 2. In Step 3 the velocity is corrected to enforce the incompressibility condition and the pressure is updated. In Step 4 the transport equations are discretized with an explicit NSSP3 Runge-Kutta scheme in time and a WENO5 scheme in space for $\phi, Y$.

For the cut-off function, we considered the following expression $\zeta(r) = \frac{1}{2}(1 + \cos(\pi r))$ on $[-1, 1]$ and $\zeta(r) = 0$ elsewhere. The Heaviside function is defined as $H(r) = \int_{-\infty}^{r} \zeta(x)dx$. In
<figure><img src="image_5.png" /><figcaption>Figure 4: Deformation of the isolines of Y at early stage rotation in the linear shear illustration: (a) without Aslam extrapolations and inner diffusion; (b) with Aslam extrapolations and inner diffusion.</figcaption></figure>

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

## 4.2. Extrapolations, reinitilization and inner diffusion

It was already discussed in section 3.2 of [8] that the numerical scheme associated to the fully Eulerian fluid-structure model with bulk (corresponding here to the terms indexed by $b$ in (21)) can be unstable in severe situations. Indeed the backward characteristics $Y$ and the level set $\phi$ are computed on the whole domain and these fields can become distorted and irregular in the vicinity of the fluid-structure interface. These perturbations can lead to numerical errors on the location and magnitude of the elastic force and therefore create instabilities. These limitations were overcome using a redistancing technique for the level set function [33, 40] together with a linear Aslam extrapolation [2] on the backward characteristics of the external fluid.

In numerical simulations of fluid-structure interaction with bulk materials, the backward characteristics can become distorted only in the outside fluid region. In the inside region, the characteristics are smooth because the elastic forces tend to bring them back to their initial configuration. However these characteristics can become irregular in the interior fluid in the case of a thin elastic membrane surrounded from either side by a fluid. To illustrate this property we consider a representative test case in capsule applications where an elastic membrane is immersed within a fluid where a linear shear configuration is considered. More information about the numerical simulation of this test case are given in Section 5.1. The membrane will elongate in one direction and then will rotate along itself in a stationary velocity field: this behavior is called the tank-treading motion. The isolines of $Y$ are plotted in Figure 4a at the early stage of rotation. It is clear that for larger physical times, i.e when the membrane will perform several full rotation, the characteristics will become even more distorted. This eventually leads also to a bad representation of the membrane.

The first idea to address this problem was to use an Aslam extrapolation of $Y$ inside the fluid with the information coming from the membrane. Since the membrane is closed, the characteristics will cross together and then give rise to a non smooth field. Instead, we propose in this paper to perform a method, we called *inner diffusion*: the backward characteristics are smoothed by using a diffusion equation in the internal region. In Figure 4b we present the same illustration where
Aslam extrapolations, and the inner diffusion were used. We clearly see that the characteristics are smooth everywhere and this allows us to run more severe configurations for longer times.

For the reinitilization technique we solve the equation for a fictitious time $\tau$,

$$ \partial_{\tau} \varphi + \operatorname{sgn}(\varphi_0)(|\nabla \varphi| - 1) = 0 $$

with a WENO5 scheme. For the linear Aslam extrapolation,

$$ \partial_{\tau} Y_n + H(\varphi)n \cdot \nabla Y_n = 0, \quad \partial_{\tau} Y + H(\varphi)(n \cdot \nabla Y - Y_n) = 0 $$

the normal derivative of $Y$ is extrapolated firstly in a constant, and then in a linear, manner in the region $\varphi > 0$ ($H$ is a Heaviside function) with the initial condition $Y_n(t=0) = n \cdot \nabla Y$. Both equations are discretized by an Euler explicit scheme in time and a WENO5 scheme in space.

For the inner diffusion approach, we solve the equation $\partial_{\tau} Y = \Delta Y$, in the region $\varphi < 0$, by discretizing it with an Euler explicit scheme in time and a centered second-order scheme in space. The fictitious time $\tau$ has been set for every test case depending on the configurations: more severe tests may need more diffusion.

# 5. Validation

In this section, two test cases have been setup with different physical parameters (size, elasticity) of the membrane in order to validate the proposed method with respect to other approaches present in the literature. One additional test case, relevant for the grid convergence of the proposed approach has been proposed in Appendix A starting from a similar test introduced in [31]. In particular, this academic test is shown to give an idea about the required mesh refinement needed to obtain a good enough resolution of the deformed capsule.

## 5.1. Simple shear flow

In this section, we present the numerical simulations for the simple shear flow test case. Thanks to the results available in the literature [25, 28, 44], this test case allows us to validate the proposed method with complex configurations, typical for capsule applications, arising when a membrane is immersed within a linear shear velocity field. The complexity comes from the fact that the steady state consists in a constant shear occurring on both the inside and the outside of the membrane. The numerical challenge of such test cases is related to evolution of the deformation vector $Y$ close to the interface. For this reason, it is mandatory for the long-time simulation of this test case to use both the Aslam extrapolation and the inner diffusion in order to capture properly the tank-treading motion (mentioned in Section 4.2).

The simulation is setup by considering a thin membrane of radius $a_m$, modeled with several constitutive laws, immersed within a viscous flow, characterized by the physical parameters in Table 1, in the domain $[-4, 4] \times [-2, 2] \times [-2, 2]$. The dimensionless parameters for this configuration read

$$ Re = \frac{\rho a_m^2 \dot{\gamma}}{\mu} = 0.0625, \quad Ca = \frac{\mu \dot{\gamma} a_m}{G_s}, \qquad (22) $$

where $Re$ is the Reynolds number, and $Ca$ is the capillary number, which is modified according to the configuration chosen. The initial conditions are given by Equation (A.1) to initialize the level-set of a sphere, and

$$ u_0(x, y, z) = \begin{pmatrix} \dot{\gamma}y \\ 0 \\ 0 \end{pmatrix}, \quad Y_0(x, y, z) = \begin{pmatrix} x \\ y \\ z \end{pmatrix}, \qquad (23) $$
<table><caption>Table 1: Simple shear flow: physical parameters.</caption><thead><tr><th>ρ</th><th>μ</th><th>a<sub>m</sub></th><th>shear rate γ̇</th></tr></thead><tbody><tr><td>1</td><td>4.0</td><td>0.5</td><td>1.0</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Simple shear flow: steady state configuration of the deformed capsule, with $L_1$ and $L_2$ represent the principal directions.</figcaption></figure>

to impose the linear shear velocity field and no pre-deformation on the initial configuration. We impose Neumann boundary conditions along the x-axis, moving wall along the y-axis, and slip wall along the z-axis.

The computational domain is discretized using a uniform grid of 512 × 256 × 256 = 33 554 432 cells. The time step chosen for these simulations is set to $\Delta t = 2 \times 10^{-3}$s. Further details about the configurations simulated in this section are given in Table 2, where different constitutive law types and multiple values of the capillary number are tested. The constitutive laws employed for these simulations are the Neo-Hookean (20) and the Skalak (19) laws. We remind here that, in general, the stretching modulus can be set as $K_s = 3G_s$. In these contributions, the shape of the membrane is assumed to be an ellipsoid once the steady state is reached and, by taking the slice $z = 0$ of this geometry, it is possible to introduce a dimensionless parameter depending on the characteristics of the ellipse,

$$D_{12} = \frac{L_1 - L_2}{L_1 + L_2}. \qquad (24)$$

This parameter, computed by taking the major and minor axes of the ellipse (see Figure 5), has been mainly introduced for validation purposes. Here, the proposed approach is compared to the literature [25, 28, 44]. In Figure 6 we show two three-dimensional visualizations of the tank-treading motion for different capillary numbers by plotting the velocity vectors of the fluid on the membrane surface. Here we show the flow velocity field over the membrane surface to point out that even if

<table><caption>Table 2: Simple shear flow: stretching and shear moduli, $K_s$ and $G_s$, when varying the capillary number for the Neo-Hookean (NK) and Skalak (SK) constitutive laws.</caption><thead><tr><th>test case</th><th>Ca</th><th>G<sub>s</sub></th><th>K<sub>s</sub></th><th>law type</th></tr></thead><tbody><tr><td>TC1</td><td>0.15</td><td>13.33</td><td>-</td><td>NH</td></tr><tr><td>TC2</td><td>0.4</td><td>5</td><td>-</td><td>NH</td></tr><tr><td>TC3</td><td>0.6</td><td>3.33</td><td>-</td><td>NH</td></tr><tr><td>TC4</td><td>0.9</td><td>2.22</td><td>-</td><td>NH</td></tr><tr><td>TC5</td><td>0.15</td><td>13.33</td><td>40</td><td>SK</td></tr><tr><td>TC6</td><td>0.4</td><td>5</td><td>15</td><td>SK</td></tr><tr><td>TC7</td><td>0.9</td><td>2.22</td><td>6.66</td><td>SK</td></tr><tr><td>TC8</td><td>2.0</td><td>1</td><td>3</td><td>SK</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6: Simple shear flow: three-dimensional visualization of the tank-treading motion at steady state.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Simple shear flow: comparisons with respect to the dimensionless parameter $D_{12}$ with Neo-Hookean (left) and Skalak (right) law types.</figcaption></figure>

the geometry of the capsule is no longer changing, the structure dynamics is still evolving. It is then clear from the visualization that the steady state configuration consists in a fixed but rotating membrane. Figure 7 presents the results of $D_{12}$ when varying the capillary number. The results presented in this section are in good agreement with those coming from the literature showing a reasonable increasing trend of $D_{12}$, when a larger capillary number is considered. In agreement with the results in Figure 7, it is possible to observe also in Figure 6 that using smaller capillary numbers leads to a stiffer membrane that tends to elongate less.

## 5.2. Capsule in a square-section channel

The test cases presented in this section have the goal of simulating capsules flowing in a square-section channel, even when the capsule size at rest is larger than the channel dimension. In these situations, an interaction between the membrane and the boundary layer occurs and the capsule deforms into a parachute shape.

The simulation is setup by considering a thin membrane, modeled with several constitutive laws, immersed within a viscous flow in a channel of dimensions $[-1, 1] \times [-0.5, 0.5] \times [-0.5, 0.5]$. The
<table><caption>Table 3: Capsule in a square-section channel: stretching and shear moduli, $K_s$ and $G_s$, when varying the ratio $a_m/ℓ$ for the Neo-Hookean (NK) and Skalak (SK) constitutive laws.</caption><thead><tr><th>test case</th><th>a<sub>m</sub>/ℓ</th><th>t<sub>0</sub></th><th>μ</th><th>G<sub>s</sub></th><th>K<sub>s</sub></th><th>law type</th></tr></thead><tbody><tr><td>TC9</td><td>0.85</td><td>0.0</td><td>8.5</td><td>85</td><td>-</td><td>NH</td></tr><tr><td>TC10</td><td>1.00</td><td>0.1</td><td>10</td><td>100</td><td>-</td><td>NH</td></tr><tr><td>TC11</td><td>0.85</td><td>0.0</td><td>8.5</td><td>85</td><td>255</td><td>SK</td></tr><tr><td>TC12</td><td>1.00</td><td>0.1</td><td>10</td><td>100</td><td>300</td><td>SK</td></tr><tr><td>TC13</td><td>1.10</td><td>0.2</td><td>11</td><td>110</td><td>330</td><td>SK</td></tr></tbody></table>

initial conditions for both the level-set field $\phi$ and deformation vector $Y$ are computed so that the capsule has a pre-deformed shape. This choice is related to the fact that, for some test cases, the capsule is assumed to be larger than the channel cross section. In this case, the initial conditions read

$$ \phi_0(x, y, z) = \sqrt{(xe^{-2t_0})^2 + (ye^{t_0})^2 + (ze^{t_0})^2} - a_m, \quad (25) $$

where $a_m$ is the radius of the relaxed sphere, and

$$ u_0(x, y, z) = \begin{pmatrix} V_p - V \\ 0 \\ 0 \end{pmatrix}, \quad Y_0(x, y, z) = \begin{pmatrix} xe^{-2t_0} \\ ye^{t_0} \\ ze^{t_0} \end{pmatrix}, \quad (26) $$

where $t_0$ is the pre-deformation time, $V$ is the mean capsule velocity and $V_p$ is a two-dimensional Poiseuille on a square section [37] such that,

$$ V_p = \frac{\pi V}{2\beta} \sum_{i=0}^{\infty} \left( \frac{1}{(2i+1)^3} - \frac{\cosh\left(\frac{(2i+1)\pi z}{\ell}\right)}{(2i+1)^3 \cosh\left(\frac{(2i+1)\pi}{2}\right)} \right) \sin\left( n\pi \left( \frac{y}{\ell} + \frac{1}{2} \right) \right), $$

where

$$ \beta = \frac{\pi^4}{96} - \sum_{i=0}^{\infty} \frac{\tan\left(\frac{(2i+1)\pi}{2}\right)}{\frac{(2i+1)^5\pi}{2}}. $$

The equations on $\phi_0$ and $Y_0$ represent those of a deformed sphere into an ellipsoid elongated, by a virtual divergence-free velocity field, along the $x$-axis and compressed along the $y$ and $z$ axes. We impose inlet and Neumann conditions along the $x$-axis, and moving wall on the other boundaries. Multiple configurations are studied by varying the constitutive laws used to model the elastic properties of the membrane and the ratio $a/ℓ$, where $2ℓ = 1$ represents the channel's size. Further details about the numerical parameters chosen in this section are given in Table 3. All simulations are set up by considering the following dimensionless parameters,

$$ Re = \frac{\rho a_m V}{\mu} = 0.05, \quad Ca = \frac{\mu V}{G_s} = 0.1, \quad (27) $$

where $\rho$ is set to one. The goal of this benchmark is analyzing the deformation of the parachute shape with respect to the influence of the constitutive law and the ratio $a_m/ℓ$. The computational domain is discretized using a uniform grid of 128 × 64 × 64 = 524 288 cells and the time step is set to $\Delta t = 10^{-3}$s. In Figure 8 we present the numerical results obtained when modeling the elastic membrane with Neo-Hookean constitutive law and varying the size of the
<figure><img src="image_9.png" /><figcaption>Figure 8: Capsule in a square-section channel: the membrane profiles obtained by the proposed method with Neo-Hookean constitutive law (dashed blue line) are compared to those presented in [22] (black dotted line).</figcaption></figure>

immersed capsule. Figure 9 shows similar numerical results obtained by using the Skalak law, also with a capsule larger than the channel (TC14) to assess our method for highly compressed capsules. Due to the small size of the channel, the capsule is highly influenced by the higher velocity profile occurring in the middle of the channel that brings the structure to take an elongated shape that resembles that of a phantom. Increasing the size of the capsule leads to an even more elongated parachute, whose shape even features corners close to the boundary layer (see Figure 8b). Little differences in the computation of the membrane profile with respect to the literature are probably related to the different numerical approaches used for the simulations. Moreover, in our approach the forces are spread near the interface, which is not the case in [22], hence regularizing the corners behind the parachute. In general, for all configurations a very good agreement with the literature is observed.

## 6. Numerical simulations of a capsule with an internal nucleus with complex geometries

In this section, we present the numerical simulations of a highly compressed capsule travelling within a square-section channel until it opens to a much wider space bringing the capsule to feature a relaxation phenomenon. A novelty with respect to previous works is the introduction of a solid bulk structure within the thin membrane. In particular, we perform thorough numerical studies on the influence of the size and stiffness of the bulk on the capsule relaxation.

The domain is divided in two part: the first one being a square-section channel $[-1, 3] \times [-0.5, 0.5] \times [-0.5, 0.5]$ where the compressed capsule is initialized and the second one being the wider part $[3, 7] \times [-1, 1] \times [-0.5, 0.5]$ where the relaxation happens (see Figure 10). The simulation is setup by considering a thin membrane, of radius $a_m$, with an internal bulk of radius $a_b$, immersed within a viscous flow. The internal square-section channel is treated as an immersed boundary by imposing wall boundary conditions through a classical penalization approach.
<figure><img src="image_10.png" /><figcaption>Figure 9: Capsule in a square-section channel: the membrane profiles obtained by the proposed method with Skalak constitutive law (red bold line) are compared to those presented in [22] (black dotted line).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Capsule relaxation phenomenon: initial configuration of a membrane with an internal solid.</figcaption></figure>
<table><caption>Table 4: Capsule relaxation phenomenon: simulation setup for the membrane.</caption><thead><tr><th>ρ</th><th>μ</th><th>V</th><th>a<sub>m</sub></th><th>G<sub>s</sub></th><th>K<sub>s</sub></th><th>law type</th></tr></thead><tbody><tr><td>1</td><td>11</td><td>1</td><td>0.55</td><td>110</td><td>330</td><td>SK</td></tr></tbody></table>

The dimensionless parameters used for this simulation are computed considering,

$$Re = \frac{\rho a_m V}{\mu} = 0.05, \quad Ca = \frac{\mu V}{G_s} = 0.1, \qquad (28)$$

and are shown in Table 4, where $2\ell$ represents the channel dimension with $\ell = 0.5$. The initial condition on the velocity is given by the two-dimensional Poiseuille on a square-section given in Section 5.2. The level-set $\phi_m$ and deformation vector $Y_m$ are initialized as

$$\phi_{m,0}(x, y, z) = \sqrt{(xe^{-2t_0})^2 + (ye^{t_0})^2 + (ze^{t_0})^2} - a_m, \quad Y_{m,0}(x, y, z) = \begin{pmatrix} xe^{-2t_0} \\ ye^{t_0} \\ ze^{t_0} \end{pmatrix}, \qquad (29)$$

where $t_0 = 0.2$. The level-set $\phi_b$ and deformation vector $Y_b$ for the elastic bulk are initialized as a sphere with no pre-deformation,

$$\phi_{b,0}(x, y, z) = \sqrt{x^2 + y^2 + z^2} - a_b, \quad Y_{b,0}(x, y, z) = \begin{pmatrix} x \\ y \\ z \end{pmatrix}. \qquad (30)$$

As mentioned above, several simulations have been performed by varying the size $a_b$ and the elastic modulus $\chi$ of the internal bulk. Further details about the setup of both $a_b$ and $\chi$ are given in Table 5.

The computational domain is discretized using a uniform grid of 512 × 128 × 64 = 4194304 cells, and the time step chosen is $\Delta t = 5 \times 10^{-4}$s.

<table><caption>Table 5: Capsule relaxation phenomenon: simulation setup for different configurations of the internal bulk.</caption><thead><tr><th>test case</th><th>χ</th><th>a<sub>b</sub>/a<sub>m</sub></th></tr></thead><tbody><tr><td>TC14</td><td>-</td><td>No bulk</td></tr><tr><td>TC15</td><td>2200</td><td>0.6</td></tr><tr><td>TC16</td><td>2200</td><td>0.5</td></tr><tr><td>TC17</td><td>1100</td><td>0.5</td></tr><tr><td>TC18</td><td>220</td><td>0.5</td></tr><tr><td>TC19</td><td>110</td><td>0.5</td></tr><tr><td>TC20</td><td>2200</td><td>0.4</td></tr><tr><td>TC21</td><td>2200</td><td>0.2</td></tr></tbody></table>

In Figure 11, we plotted the evolution of the membrane profile, with and without an internal bulk. At the beginning, as presented in Section 5.2, when flowing within the square-section channel, the compressed capsule takes a phantom shape. Then, once the capsule is at the end of the channel, it starts relaxing and modifying its shape until it becomes again similar to a sphere. To give more insights on the simulations, in Figure 12 we define three characteristic lengths to study the evolution of the parachute during the relaxation phenomenon. Figure 13 shows the evolution in
<figure><img src="image_12.png" /><figcaption>Figure 11: Capsule relaxation phenomenon: evolution of the capsule profile at different times.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Capsule relaxation phenomenon: $L_1, L_2, L_3$ are the characteristic length of deformed capsule.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Capsule relaxation phenomenon: time evolution of $L_1/\ell$ and $L_3/\ell$ when varying the size of the internal bulk, while keeping a constant elastic modulus $\chi = 2200$. The numerical results of TC14 are depicted using a dashed line, TC15 in blue, and TC16 in red.</figcaption></figure>

non-dimensional time $tV/\ell$ of $L_1/\ell$ and $L_3/\ell$ when varying the bulk's size, while keeping a constant elastic modulus $\chi = 2200$ (TC14, TC15, and TC16 in Table 5). As expected, in all cases $L_1/\ell$ suddenly decreases and features a small jump when the capsule is exiting the channel. After that, the relaxation process begins and $L_1/\ell$ increases again. Figure 13 shows that the presence of an elastic bulk slightly influences the evolution of $L_1/\ell$. Similar trends are observed for $L_3/\ell$, however it should be noticed that in this case the presence of an elastic bulk introduces remarkable changes in the parachute shape, which appears much more elongated for larger bulks. In particular, it is observed that when increasing the bulk's size, the minimum of $L_3/\ell$ increases too. Although the dynamics of the relaxation process is highly influenced by the internal bulk, it is interesting to observe that the minimum value of $L_3/\ell$ always occurs around $tV/\ell = 6.2$. In Figure 14 we plot the evolution of $L_2/\ell$ where the strong correlation with $L_1/\ell$ is observed: when $L_1/\ell$ decreases, $L_2/\ell$ increases and vice versa.

Figure 15 shows the time evolution of the $L_3/\ell$ characteristic length when varying the elastic modulus of the internal bulk, while keeping a constant radius $a_b/a_m = 0.5$ (TC14, TC16, TC17, TC18, and TC19). Similar trends to Figure 13b are followed when performing analysis on the bulk stiffness. In particular, it should be noticed that the stiffer the internal bulk is the higher the minimum of $L_3/\ell$ gets. As already observed in Figure 13b, even though the capsule dynamics is highly affected by the presence of the internal solid, the minimum value of $L_3/\ell$ always occurs at the same time.

More details about the influence of an internal bulk on the capsule's shape is given in Figure 16, where the capsule dynamics is depicted at different times for different configurations (TC15, TC16, TC20, and TC21). For this case, an internal bulk with $\chi = 2200$ is considered for the simulations. This brings the solid not to deform that much during the relaxation process. However, the bulk's size remarkably influences the shape of the parachute that appears differently. In particular, the smaller the bulk is, the more the capsule resembles the phantom shape shown in Figure 9a. Instead, when a bigger bulk is considered the shape of the parachute takes a more deformed profile.

Finally, in Figure 17 a fixed ratio $a_b/a_m = 0.5$ is considered, and the stiffness of the bulk varies from 110 to 2200. In this case, the bulk now shows larger deformations because a lower elastic modulus is chosen. This brings the parachute to slightly modify its shape assuming a more, or
<figure><img src="image_15.png" /><figcaption>Figure 14: Temporal evolution of the length $L_2/\ell$ for three test case: no elastic volume (TC15), $a_b/a = 0.6$ (TC16) and $a_b/a = 0.5$ (TC17). The elastic volume coefficient is $\chi = 2200$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Capsule relaxation phenomenon: time evolution of $L_3$ when varying the elastic modulus of the internal bulk, while keeping a constant radius ratio $a_b/a = 0.5$. The numerical results of TC14 are depicted using a dashed line, TC16 in blue, TC17 in red, TC18 in orange, and TC19 in green.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: Capsule relaxation phenomenon: time evolution of the membrane and bulk profiles during the relaxation process when varying the size of the internal bulk, while keeping a constant elastic modulus $\chi = 2200$. The numerical results of TC15 are depicted in blue, TC16 in red, TC20 in orange, and TC21 in green. The walls of the square-section channel are represented in grey.</figcaption></figure>

less, elongated appearance. This also seems to be in line with the time evolution of $L_3/\ell$ presented in Figure 15.

## 7. Conclusion and discussion

In this paper, we have presented a numerical approach for the simulation of elastic capsules. The numerical model is based on a fully Eulerian formulation for both the fluid and the deformable structures, membrane and solid volume, that are here taken into account by the Navier-Stokes equation by means of a source term describing the elastic forces, which act on the fluid. The fluid-structure system of equations considered herein describes both the forces related to the area variation and shear by following the full membrane model introduced in [31]. As mentioned above, the evaluation of the backward characteristics used to compute the elastic deformations can degrade when complex simulations of capsules are addressed, especially when membrane are considered. In this case, distortions of backward characteristics may occur within the membrane, causing a less reliable prediction of the capsule dynamics. The novelty of our method is based on coupling the Aslam extrapolation to improve the isolines of the backward characteristics outside the capsule (as done in [8]) with the *inner diffusion* approach to smooth the internal isolines. Unfortunately, the good performances provided by the Aslam extrapolation in the outer zone does not solve the same problem in the internal part.

The new algorithm is validated on several complex 3D configurations where the capsule dynamics is simulated and compared to other numerical experiments present in the literature. In particular, the challenges of these test cases are such that the proposed algorithm is essential to simulate the correct evolution of the capsule. First, we focus on the validation of the proposed model and method for classical test cases related to capsule dynamics. We studied before the deformation of a capsule when immersed into a linear shear velocity field, which gives rise to the tank-treading motion. Then, we performed several simulations of a capsule in a square-section channel, which allows to study different parachute shapes that the capsule takes when highly compressed. The numerical
<figure><img src="image_18.png" /><figcaption>Figure 17: Capsule relaxation phenomenon: time evolution of the membrane and bulk profiles at different times for fixed radius ratio $a_b/a_m = 0.5$ and variable elastic modulus of the bulk. The numerical results of TC16 are depicted in blue, TC17 in red, TC18 in orange, and TC19 in green. The walls of the square-section channel are represented in grey.</figcaption></figure>

results obtained for different capillary numbers and constitutive laws have been compared with other numerical techniques present in the literature, obtaining a good agreement in all cases. In Appendix A, we also present a grid refinement analysis for an academic test case to understand the needed level of refinement to achieve satisfying results. In particular, a modification of the circular shear benchmark introduced in [31] is proposed by introducing a pressurization coefficient that allows, with the right constitutive law, the correct recovery of the spherical shape after its relaxation. Finally, an original test case is presented where an elastic nucleus is considered within the thin membrane. We provide interesting insights about the influence of a nucleus on the deformation of a capsule during its relaxation after a compressed phase.

## Appendix A. Relaxation of a sheared elastic sphere

This academic test case is a modification of the sheared elastic sphere relaxation proposed in [31]. Herein, we propose to use again this test case to validate the new method described in this paper by performing a mesh refinement analysis. As it was shown in [31], the relaxation of the pre-deformed sphere ends with a final configuration that does no longer resemble an actual sphere. This is mainly related to the fact that the sphere was initialized with no area variation ($Z_1 = 1$), which eventually brings the membrane to modify its aspect.

The simulation is setup by considering a thin membrane, modeled with the Evans-Skalak constitutive law (18) immersed within a viscous flow in the domain $[-1, 1]^3$. It should be noticed that also the constitutive law has changed with respect to the test case presented in [31]. The physical parameters used for this test case are given in Table A.6. We impose zero velocity as initial condition and Neumann boundary conditions on all boundaries of the domain. The initial conditions for the level-set field $\phi$ is given by a sphere of radius $a_m$,

$$ \phi_0(x, y, z) = \sqrt{x^2 + y^2 + z^2} - a_m. \qquad (A.1) $$
<table><caption>Table A.6: Relaxation of a sheared elastic sphere: physical parameters.</caption><thead><tr><th>ρ</th><th>μ</th><th>G<sub>s</sub></th><th>K<sub>s</sub></th><th>a<sub>m</sub></th></tr></thead><tbody><tr><td>1</td><td>0.01</td><td>0.1</td><td>1</td><td>0.5</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>Figure A.18: Relaxation of a sheared elastic sphere: initial condition of $Z_2$.</figcaption></figure>

The membrane is pre-deformed with a 3D circular shear until a fictitious time $t_0 = \pi$. This corresponds to the following initial conditions on $Y$,

$$Y_0(x, y, z) = \frac{1}{1+\alpha} \begin{pmatrix} x \cos(t_0 z) + y \sin(t_0 z) \\ y \cos(t_0 z) - x \sin(t_0 z) \\ z \end{pmatrix} \qquad (A.2)$$

where $\alpha$ is the pressurization coefficient. Following the Equation (14), the initial area and shear variations reads,

$$Z_{1,0} = (1 + \alpha)^2, \quad Z_{2,0} = 1 + \frac{t_0^2(x^2 + y^2)^2}{2(x^2 + y^2 + z^2)}. \qquad (A.3)$$

This introduces a shear variation that is equal to zero at the poles and varies along the z-axis (see Figure A.18). It should be noticed that, $\alpha$ does not influence the shear variation which stays the same, but introduce an area variation on the initial conditions. The coefficient $\alpha$ has been introduced to obtain a sphere shape after the relaxation process. Moreover, when taking $\alpha = 0$ we obtain the same configuration analyzed in [31].

The computational domain is discretized using a Cartesian mesh with $N = \{64, 128, 256\}$ cells in each direction. The time step chosen for the simulations varies depending on the mesh refinement: $\Delta t = 4 \times 10^{-3}$ for $N = 64$, $\Delta t = 2 \times 10^{-3}$ for $N = 128$ and $\Delta t = 10^{-3}$s for $N = 256$. In this simulation the pressurization coefficient $\alpha$ is set to 0.05.

The numerical results at different time steps are presented in Figure A.19, showing the deformation isocontours until the new equilibrium is reached. At the beginning of the simulation, only the force related to the shear coefficient is acting on the membrane. However, the membrane inertia introduces a shape modification that brings into play also the force term related to the area variation. As expected, the shear deformation gives rise to a notable deformation and inertia, along both the $x$ and $z$ axes, such that the membrane begins to swing to one side and then to the other. From $t = 3s$, the membrane starts stabilizing and slightly oscillating around the steady state
<figure><img src="image_20.png" /><figcaption>Figure A.19: Relaxation of a sheared elastic sphere: surface profile with the deformation isocontours $Y_1$ (red) and $Y_3$ (black) at different time steps.</figcaption></figure>

where $Z_2$ is uniformly equal to one, and $Z_1 = Z_{1,0}$. Because of the incompressibility condition, the internal fluid is going to be characterized by the same volume, that translates into an area variation at the equilibrium equal to the initial one. Having introduced the pressure coefficient $\alpha$, the equilibrium configuration has an area variation greater than one that implies an elastic force at the interface. The pressure jump resulting from the elastic force is the reason why we are able to return to an actual sphere configuration at the equilibrium (which was not achieved in [31] where $\alpha$ was virtually set to zero). In order to provide more insights about these oscillations, in Figure A.20, we plotted the evolution of the horizontal and vertical radii. These curves also give an idea about the mesh refinement needed to achieve good enough results. The dimensionless pressure along the $x$ and the $z$ axes is also studied, in Figures A.21 and A.22 to provide more quantitative results on the mesh convergence. It is possible to notice that different results, with respect to those in [31], are found here because of aforementioned differences in the initial pressurization of the membrane and elastic constitutive law.

## References

[1] Alizad Banaei, A., Loiseau, J.C., Lashgari, I., Brandt, L., 2017. Numerical simulations of elastic capsules with nucleus in shear flow. European Journal of Computational Mechanics 26, 131–153.

[2] Aslam, T.D., 2003. A partial differential equation approach to multidimensional extrapolation. Journal of Computational Physics 193, 349–355.
<figure><img src="image_21.png" /><figcaption>Figure A.20: Relaxation of a sheared elastic sphere: evolution of the x-radius and z-radius for $N = \{64, 128, 256\}$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure A.21: Relaxation of a sheared elastic sphere: dimensionless pressure $pa/G_s$ along the x-axis at time $t = 0.1s$ (left), and $t = 1.2$ (right) for $N = \{64, 128, 256\}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure A.22: Relaxation of a sheared elastic sphere: dimensionless pressure $pa/G_s$ along the z-axis at time $t = 0.1s$ (left), and $t = 1.2$ (right) for $N = \{64, 128, 256\}$.</figcaption></figure>
[3] Bergmann, M., Fondanèche, A., Iollo, A., 2022. An eulerian finite-volume approach of fluid-structure interaction problems on quadtree meshes. Journal of Computational Physics 471, 111647.

[4] Cottet, G.H., Maitre, E., 2006. A level set method for fluid-structure interactions with immersed surfaces. Mathematical models and methods in applied sciences 16, 415-438.

[5] Cottet, G.H., Maitre, E., Milcent, T., 2008. Eulerian formulation and level set models for incompressible fluid-structure interaction. ESAIM: Mathematical Modelling and Numerical Analysis 42, 471-492.

[6] Cottet, G.H., Maitre, E., Milcent, T., 2021. Méthodes Level Set pour l'interaction fluide-structure. volume 86. Springer.

[7] De Brauer, A., Iollo, A., Milcent, T., 2017. A cartesian scheme for compressible multimaterial hyperelastic models with plasticity. Communications in Computational Physics 22, 1362-1384.

[8] Deborde, J., Milcent, T., Lubin, P., Glockner, S., 2020. Numerical simulations of the interaction of solitary waves and elastic structures with a fully eulerian method. Water Waves 2, 433-466.

[9] Del Mercato, L.L., Ferraro, M.M., Baldassarre, F., Mancarella, S., Greco, V., Rinaldi, R., Leporatti, S., 2014. Biological applications of lbl multilayer capsules: From drug delivery to sensing. Advances in colloid and interface science 207, 139-154.

[10] Desmons, F., Coquerelle, M., 2021. A generalized high-order momentum preserving (homp) method in the one-fluid model for incompressible two phase flows with high density ratio. Journal of Computational Physics 437, 110322.

[11] Donea, J., Giuliani, S., Halleux, J.P., 1982. An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions. Computer methods in applied mechanics and engineering 33, 689-723.

[12] Falgout, R., Jones, J., Yang, U., 2006. 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, 267-294.

[13] Falgout, R., Yang, U., 2002. 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, 632-641.

[14] Fernández, M.A., Gerbeau, J.F., Grandmont, C., 2007. A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid. International Journal for Numerical Methods in Engineering 69, 794-821.

[15] Fernández, M.A., Moubachir, M., 2005. A newton method using exact jacobians for solving fluid-structure coupling. Computers & Structures 83, 127-142.
[16] Foessel, E., Walter, J., Salsac, A.V., Barthès-Biesel, D., 2011. Influence of internal viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow. Journal of Fluid Mechanics 672, 477-486.

[17] Gires, P.Y., Barthès-Biesel, D., Leclerc, E., Salsac, A.V., 2016. Transient behavior and relaxation of microcapsules with a cross-linked human serum albumin membrane. Journal of the mechanical behavior of biomedical materials 58, 2-10.

[18] Gorsse, Y., Iollo, A., Milcent, T., Telib, H., 2014. A simple cartesian scheme for compressible multimaterials. Journal of Computational Physics 272, 772-798.

[19] Griffith, B.E., Patankar, N.A., 2020. Immersed methods for fluid-structure interaction. Annual review of fluid mechanics 52, 421-448.

[20] Guermond, J.L., Minev, P., Shen, J., 2006. An overview of projection methods for incompressible flows. Computer Methods in Applied Mechanics and Engineering 195, 6011-6045.

[21] Holzapfel, G., 2000. Nonlinear Solid Mechanics. A continuum approach for engineering. J. Wiley and Sons.

[22] Hu, X.Q., Sévénié, B., Salsac, A.V., Leclerc, E., Barthès-Biesel, D., 2013. Characterizing the membrane properties of capsules flowing in a square-section microfluidic channel: Effects of the membrane constitutive law. Phys. Rev. E 87, 063008. URL: https://link.aps.org/doi/10.1103/PhysRevE.87.063008, doi:10.1103/PhysRevE.87.063008.

[23] Jiang, G., Shu, C., 1996. Efficient implementation of weighted eno schemes. Journal of Computational Physics.

[24] Ketcheson, D.I., 2008. Highly efficient strong stability-preserving runge-kutta methods with low-storage implementations. SIAM Journal on Scientific Computing 30, 2113-2136.

[25] Lac, E., Barthès-Biesel, D., Pelekasis, N.A., Tsamopoulos, J., 2004. Spherical capsules in three-dimensional unbounded stokes flows: effect of the membrane constitutive law and onset of buckling. Journal of Fluid Mechanics 516, 303-334. doi:10.1017/S002211200400062X.

[26] Lee, L., LeVeque, R.J., 2003. An immersed interface method for incompressible navier-stokes equations. SIAM Journal on Scientific Computing 25, 832-856.

[27] Li, M., Qiao, S., Zheng, Y., Andaloussi, Y.H., Li, X., Zhang, Z., Li, A., Cheng, P., Ma, S., Chen, Y., 2020. Fabricating covalent organic framework capsules with commodious microenvironment for enzymes. Journal of the American Chemical Society 142, 6675-6681.

[28] Li, X., Sarkar, K., 2008. Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane. Journal of Computational Physics 227, 4998-5018. URL: https://www.sciencedirect.com/science/article/pii/S0021999108000569, doi:https://doi.org/10.1016/j.jcp.2008.01.034.

[29] Liu, H., Kawachi, K., 1998. A numerical study of insect flight. Journal of computational physics 146, 124-156.

[30] Maitre, E., Milcent, T., Cottet, G.H., Raoult, A., Usson, Y., 2009. Applications of level set methods in computational biophysics. Mathematical and Computer Modelling 49, 2161-2169.
[31] Milcent, T., Maitre, E., 2016. Eulerian model of immersed elastic surfaces with full membrane elasticity. Communications in Mathematical Sciences 14, 857-881. doi:https://dx.doi.org/10.4310/CMS.2016.v14.n3.a11.

[32] Mittal, R., Iaccarino, G., 2005. Immersed boundary methods. Annu. Rev. Fluid Mech. 37, 239-261.

[33] M.Sussman, Smereka, P., Osher, S., 1994. A level set approach for computing solutions to incompressible two-phase flow. Journal of Computational Physics 114, 146-159.

[34] Osher, S., Fedkiw, R., Piechor, K., 2004. Level set methods and dynamic implicit surfaces. Appl. Mech. Rev. 57, B15-B15.

[35] Peskin, C.S., 1972. Flow patterns around heart valves: a numerical method. Journal of computational physics 10, 252-271.

[36] Peskin, C.S., 2002. The immersed boundary method. Acta numerica 11, 479-517.

[37] Pozrikidis, C., 1997. Introduction to Theoretical and Computational Fluid Dynamics. EngineeringPro collection, OUP USA. URL: https://books.google.fr/books?id=Wys9n5yB18EC.

[38] Pozrikidis, C., 2001. Interfacial dynamics for stokes flow. Journal of Computational Physics 169, 250-301.

[39] Ramanujan, S., Pozrikidis, C., 1998. Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of fluid viscosities. Journal of fluid mechanics 361, 117-143.

[40] Russo, G., Smereka, P., 2000. A remark on computing distance functions. Journal of Computational Physics 163, 51-67. URL: https://www.sciencedirect.com/science/article/pii/S0021999100965537, doi:https://doi.org/10.1006/jcph.2000.6553.

[41] Sahin, M., Mohseni, K., 2009. An arbitrary lagrangian-eulerian formulation for the numerical simulation of flow patterns generated by the hydromedusa aequorea victoria. Journal of Computational Physics 228, 4588-4605.

[42] Sotiropoulos, F., Yang, X., 2014. Immersed boundary methods for simulating fluid-structure interaction. Progress in Aerospace Sciences 65, 1-21.

[43] Walter, J., Salsac, A.V., Barthès-Biesel, D., Le Tallec, P., 2010a. Coupling of finite element and boundary integral methods for a capsule in a stokes flow. International journal for numerical methods in engineering 83, 829-850.

[44] Walter, J., Salsac, A.V., Barthès-Biesel, D., Le Tallec, P., 2010b. Coupling of finite element and boundary integral methods for a capsule in a stokes flow. International Journal for Numerical Methods in Engineering 83, 829-850. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2859, doi:https://doi.org/10.1002/nme.2859.