# Numerical Study of a 3D Eulerian Monolithic Formulation for Incompressible Fluid-Structures Systems

Chen-Yu Chiang, Olivier Pironneau, Tony W H Sheu, Marc Thiriet

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

## To cite this version:

Chen-Yu Chiang, Olivier Pironneau, Tony W H Sheu, Marc Thiriet. Numerical Study of a 3D Eulerian Monolithic Formulation for Incompressible Fluid-Structures Systems. Fluids, 2017, Mechanics of Fluid-Particles Systems and Fluid-Solid Interactions, 2 (2), pp.34-53. 10.3390/fluids2020034. hal-01558747

HAL Id: hal-01558747

https://hal.sorbonne-universite.fr/hal-01558747

Submitted on 10 Jul 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Numerical Study of a 3D Eulerian Monolithic Formulation for Incompressible Fluid-Structures Systems

Chen-Yu Chiang<sup>1,2,†</sup>, Olivier Pironneau<sup>2,*†</sup>, Tony W. H. Sheu<sup>3</sup> and Marc Thiriet<sup>2,4,†</sup>

Received: 27 March 2017; Accepted: 8 June 2017; Published: date

**Abstract:** An algorithm is derived for a hyperelastic incompressible solid coupled with a Newtonian fluid. It is based on a Eulerian formulation of the full system in which the main variables are the velocities. After a fully implicit discretization in time it is possible to eliminate the displacements and solve a variational equation for the velocities and pressures only. The stability of the method depends heavily on the use of characteristic-Galerkin discretization of the total derivatives. For comparison with previous works, the method is tested on a three-dimensional (3D) clamped beam in a pipe filled with a fluid. Convergence is studied numerically on an axisymmetric case.

**Keywords:** fluid-structure-interaction; Eulerian formulation; monolithic; finite element; characteristic-Galerkin; large displacement; energy compatible algorithm.

## 1. Introduction

Computer simulation of fluid-structure systems in interaction (FSI) is important for aerospace [1], hemodynamics [2,3], automobile industries [4], etc.

A popular approach to numerically match the velocities and stresses at the fluid-solid interfaces is the arbitrary Lagrangian Eulerian formulation (ALE) whereby at every time step the fluid equations are mapped back into the solid domain [5]. Difficulties with the mesh used in the fluid part arise in the case of large displacements of the solids [6,7].

The immersed boundary method (IBM) proposed by Charles Peskin [2] and analyzed in [8] is efficient for shells in fluid but is more difficult to implement for thick structures [9]. Particles in suspension in a fluid have been simulated by IBM-like methods but mostly for hard particles (see for instance [10,11]).

A third way is the fully Eulerian approach, advocated by Liu et al. [12] and Dunne-Rannacher [13,14] (see also [15,16]). This approach is well suited to problems with large displacements.

In [17,18] a similar monolithic numerical method was proposed for the fluid-structure equations in two dimensions. It was also shown to be energy-stable. For more details on the motivations behind such a radically different approach to FSI, see [18]. The basic idea is to write the equations in terms of velocity, pressure, and displacement, as usual but in an Eulerian frame, and then eliminate

<footnote n="1">1 Department of Engineering Science and Ocean Engineering, National Taiwan University; Taipei 10617, Taiwan; f01525007@ntu.edu.tw</footnote>

<footnote n="2">2 Laboratoire Jacques-Louis-Lions (LJLL), Sorbonne University, UPMC (Paris VI), 75005 Paris, France; marc.thiriet@inria.fr</footnote>

<footnote n="3">3 Institute of Applied Mathematical Sciences and Center of Advanced Study in Theoretical Sciences (CASTS), National Taiwan University, Taipei 10617, Taiwan; twhsheu@ntu.edu.tw</footnote>

<footnote n="4">4 INRIA-Paris-Rocquencourt, EPC-REO project, Domaine de Voluceau, BP105, 78153 Le Chesnay Cedex, France</footnote>

<footnote>* Correspondence: olivier.pironneau@upmc.fr; Tel.: +33-1-4427-4411</footnote>

<footnote>† Current address: LJLL-UPMC, boite 187, 4 Place Jussieu, F75005 Paris, France.</footnote>
the displacement in the time-discretized equations; the result is a monolithic formulation which is well-posed at each time step. Furthermore, to obtain energy stability it was suggested in [18] to update the structure with its own velocity and remesh the fluid domain at every time step.

In this paper we keep the same approach but in 3D. The Cayley–Hamilton theorem was used in 2D and it is different in 3D. Hence, the equations for a structure in Eulerian variables need to be derived again. The analysis is performed for hyper-elastic Mooney–Rivlin incompressible material. The result is not very different except that one coefficient of the constitutive relation is now a nonlinear function of the deformation tensor.

The article has four sections.

* The fully Eulerian formulation is derived in Section 1.

* In Section 2 the equations are discretized implicitly in time as in [18,19]; then the displacements are eliminated and a non-linear system for the velocities and pressures remains. Spatial discretization is obtained by using stable finite element spaces like linear pressures and quadratic velocities on a tetrahedral mesh.

* In Section 3 it is shown that the energy decays at each time step, which is an indication that the method is robust.

* Finally, in Section 4, the method is implemented in 2D for axisymmetric systems and in 3D for general systems. Energy decrease, mass conservation and convergence are analyzed on an axisymmetric case: an elastic torus in a cylindrical canister filled with a fluid at a Reynolds number of a few hundred. Then, the method is carefully evaluated on the test case proposed in [20] and comparisons with previous publications are made.

## 2. Derivation of the Formulation from the General Laws of Continuum Mechanics

### 2.1. Notations

The computational domain is denoted $\Omega^t$; it is time-dependent and is the union of the fluid part $\Omega_f^t$ and solid part $\Omega_s^t$. Technically $\Omega^t$ must be an open set and so we denote by $\overline{\Omega}^t$ its closure. The fluid and solid domains must not overlap: $\Omega_f^t \cap \Omega_s^t = \emptyset$. Initially, $\Omega_f^0$ and $\Omega_s^0$ are prescribed. Let:

* $\Sigma^t = \overline{\Omega_f^t} \cap \overline{\Omega_s^t}$ be the fluid-structure interface, and $\partial\Omega^t$ be the boundary of $\Omega^t$,

* $\Gamma$ be the part of $\partial\Omega^t$ where the structure is clamped or the fluid does not slip. It is assumed to be independent of $t$.

The following standard notations are used:

* $\mathbf{X} : \Omega^0 \times (0, T), \mapsto \Omega^t : \mathbf{X}(x^0, t)$ is the Lagrangian position at $t$ of $x^0$,

* $\mathbf{d} = \mathbf{X}(x^0, t) - x^0$ is the displacement,

* $\mathbf{u}(X(x^0, t), t) = \partial_t X(x^0, t)$ is the Eulerian velocity of the deformation at $t$ and $x = X(x^0, t)$,

* $\mathbf{F}_{ji} = \partial_{x^i} \mathbf{X}_j$ is the transposed gradient of the deformation,

* $J = \det \mathbf{F}$ is the Jacobian of the deformation.

Let $\rho_f$, $\sigma_f$ be the density and stress tensor in the fluid, as for the solid with $\rho_s$, $\sigma_s$, respectively. It is convenient to define a unique density and stress tensor by using the set function indicators $1_{\Omega_f^t}$ and $1_{\Omega_s^t}$:

* $\rho(x, t) = 1_{\Omega_f^t}(x, t)\rho_f(x, t) + 1_{\Omega_s^t}(x, t)\rho_s(x, t)$, the density,

* $\sigma(x, t) = 1_{\Omega_f^t}(x, t)\sigma_f(x, t) + 1_{\Omega_s^t}(x, t)\sigma_s(x, t)$, the stress tensor.

For readability, vectors, tensors and matrices are noted in bold, except $x$ and $x^0$. Unless specified otherwise, all spatial derivatives are with respect to $x \in \Omega^t$ and not with respect to $x^0 \in \Omega^0$. If $\phi$ is a function of $x = \mathbf{X}(x^0, t)$,

$$ \nabla_{x^0} \phi = [\partial_{x_i^0} \phi] = [\partial_{x_i^0} \mathbf{X}_j \partial_{x_j} \phi] = \mathbf{F}^T \nabla \phi \quad (1) $$
When $\mathbf{X}$ is invertible, $\mathbf{d}$ and $\mathbf{F}$ can be seen as functions of $(x, t)$ instead of $(x^0, t)$. They are related by:

$$ \mathbf{F}^T = \nabla_{x^0} \mathbf{X} = \nabla_{x^0} (\mathbf{d} + x^0) = \nabla_{x^0} \mathbf{d} + \mathbf{I} = \mathbf{F}^T \nabla \mathbf{d} + \mathbf{I} \Rightarrow \mathbf{F} = (\mathbf{I} - \nabla \mathbf{d})^{-T} \quad (2) $$

Time derivatives are related by:

$$ \mathbb{D}_t \phi(x, t) := \frac{d}{dt} \phi \left( \mathbf{X} \left( x^0, t \right), t \right) \bigg|_{\mathbf{X}(x^0, t)=x} = \partial_t \phi(x, t) + \mathbf{u} \cdot \nabla \phi(x, t) \quad (3) $$

It is convenient to introduce a notation,

$$ \mathbf{D} \mathbf{u} = \nabla \mathbf{u} + \nabla^T \mathbf{u}. $$

## 2.2. Conservation Laws

Conservation of momentum and conservation of mass take the same form for the fluid and the solid. With $\mathbf{f}(x, t)$, the volumic force in the system is:

$$ \rho \mathbb{D}_t \mathbf{u} = \mathbf{f} + \nabla \cdot \boldsymbol{\sigma}, \quad \frac{d}{dt} (J\rho) = 0, \text{ i.e., } J\rho = \rho^0 $$

with continuity of $\mathbf{u}$ and of $\boldsymbol{\sigma} \cdot \mathbf{n}$ at the fluid-structure interface $\boldsymbol{\sigma}^t$ in the absence of external surface force, like surface tension. There are also unwritten constraints pertaining to the realizability of the map $\mathbf{X}$. Finally, incompressibility implies $J = 1$ and so $\rho = \rho^0$ along the Lagrangian trajectories. In particular, if $\rho^0$ is piecewise constant and equal to $\rho_f$, $\rho_s$ in the fluid and the solid, at initial time, then it remains so at later times.

## 2.3. Constitutive Equations

Three-dimensional hyperelastic incompressible Mooney-Rivilin/neo Hookean materials and Newtonian incompressible fluids are considered.

* For a Newtonian incompressible fluid, $\boldsymbol{\sigma}_f = -p_f \mathbf{I} + \mu \mathbf{D} \mathbf{u}$

* For an hyperelastic incompressible material, $\boldsymbol{\sigma}_s = -p_s \mathbf{I} + \partial_F \Psi \mathbf{F}^T$

where $\Psi$ is an Helmholtz potential, like,

$$ \Psi(\mathbf{F}) = c_1 \mathrm{tr}_{\mathbf{F}^T \mathbf{F}} + c_2 \left( \mathrm{tr}_{(\mathbf{F}^T \mathbf{F})^2} - \mathrm{tr}_{\mathbf{F}^T \mathbf{F}}^2 \right). \quad (4) $$

## 2.4. The Mooney-Rivilin 3D Stress Tensor

Note that $\partial_F \mathrm{tr}_{\mathbf{F}^T \mathbf{F}} = 2\mathbf{F}$ and $\partial_F \mathrm{tr}_{(\mathbf{F}^T \mathbf{F})^2} = 4\mathbf{F} \mathbf{F}^T \mathbf{F}$. Hence,

$$ \partial_F \Psi \mathbf{F}^T = (2c_1 - 4c_2 \mathrm{tr}_B) \mathbf{B} + 4c_2 \mathbf{B}^2, \quad \text{where } \mathbf{B} = \mathbf{F} \mathbf{F}^T. \quad (5) $$

In order to compute the tensors at Eulerian points $(x, t)$ rather than Lagrangian points $(x^0, t)$, it is necessary to compute $x \to \mathbf{F}(x, t)$ or $x \to \mathbf{d}(x, t)$ in terms of $\Psi$. By (2)

$$ \mathbf{B}^{-1} = \mathbf{F}^{-T} \mathbf{F}^{-1} = (\mathbf{I} - \nabla \mathbf{d})(\mathbf{I} - \nabla \mathbf{d})^T = \mathbf{I} - \mathbf{D} \mathbf{d} + \nabla \mathbf{d} \nabla^T \mathbf{d}. \quad (6) $$

By Cayley-Hamilton theorem in 3D, for an invertible symmetric matrix,

$$ \mathbf{B}^3 - \mathrm{tr}_{\mathbf{B}} \mathbf{B}^2 + \gamma \mathbf{B} - \det_{\mathbf{B}} \mathbf{I} = 0 \quad \text{with} \quad \gamma = \frac{1}{2} (\mathrm{tr}_{\mathbf{B}}^2 - \mathrm{tr}_{\mathbf{B}^2}), \quad (7) $$
and by incompressibility det$_{B}$ = 1. Multiplying (7) by B$^{-2}$, and the result by B$^{-1}$ leads to polynomial expressions in terms of B$^{-1}$:

$$
\mathbf{B} = \mathrm{tr}_{\mathbf{B}} \mathbf{I} - \gamma \mathbf{B}^{-1} + \mathbf{B}^{-2}, \quad \mathbf{B}^2 = (\mathrm{tr}_{\mathbf{B}}^2 - \gamma) \mathbf{I} + (1 - \mathrm{tr}_{\mathbf{B}} \gamma) \mathbf{B}^{-1} + \mathrm{tr}_{\mathbf{B}} \mathbf{B}^{-2}. \quad (8)
$$

With the help of (8), (5) becomes:

$$
\begin{align*}
\partial_F \Psi \mathbf{F}^T &= (2c_1 - 4c_2 \mathrm{tr}_B) \mathbf{B} + 4c_2 \mathbf{B}^2 \\
&= (2c_1 - 4c_2 \mathrm{tr}_B) (\mathrm{tr}_B \mathbf{I} - \gamma \mathbf{B}^{-1} + \det_B \mathbf{B}^{-2}) \\
&\quad + 4c_2 [(\mathrm{tr}_B^2 - \gamma) \mathbf{I} + (1 - \mathrm{tr}_B \gamma) \mathbf{B}^{-1} + \mathrm{tr}_B \mathbf{B}^{-2}] \\
&= 2c_1 \mathbf{B}^{-2} + (4c_2 - 2c_1 \gamma) \mathbf{B}^{-1} + (2c_1 \mathrm{tr}_B - 4c_2 \gamma) \mathbf{I}.
\end{align*}
$$

Consequently, as a function of **d**, the 3D Mooney–Rivlin stress tensor is obtained from (6),

$$
\partial_F \Psi \mathbf{F}^T = 2c_1 (\mathbf{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d})^2 + 2c_3 (\mathbf{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) + \alpha \mathbf{I}, \quad (9)
$$

$$
c_3 := \frac{c_1}{2} (\mathrm{tr}_{\mathrm{B}}^2 - \mathrm{tr}_{\mathrm{B}}^2 - 4) - 4c_2 \quad (10)
$$

with a complex expression for α in (10) which is unimportant for the time being because it is absorbed
by the pressure.

This is compared with the 2D model of [19] obtained with the same Ψ:

$$
\partial_F \Psi \mathbf{F}^T = \begin{cases} 2c_1 (\mathbf{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) + \alpha' \mathbf{I}, & \text{2D model.} \\ 2c_1 (\mathbf{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d})^2 + 2c_3(x,t) (\mathbf{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) + \alpha \mathbf{I}, & \text{3D model.} \end{cases}
$$

with $c_3$ given by (10).

2.5. From 3D to 2D

In the 3D model, $c_3$ varies with $\mathrm{tr}_B$ and $\mathrm{tr}_{B^2}$ instead of being constant, like $c_1$ in 2D model. The way to calculate $c_3$ is to compute directly $\mathrm{tr}_B$ and $\mathrm{tr}_{B^2}$ by (6) which involves a 3 by 3 inverse matrix at each point of the solid domain.

The 3D model degenerates into a 2D system when the geometry and the data are invariant with
respect to one coordinate, like a translation invariance with respect to z for instance, in a Cartesian
coordinate system, or a rotation invariance in θ, in a cylindrical coordinate system.

Assume invariance with respect to the third coordinate. Then,

$$
\mathbf{F} = \begin{pmatrix} \mathbf{I} - \nabla_2 \mathbf{d} & \mathbf{0} \\ \mathbf{0} & \mathbf{1} \end{pmatrix}^{-T} = \begin{pmatrix} \tilde{\mathbf{F}} & \mathbf{0} \\ \mathbf{0} & \mathbf{1} \end{pmatrix} \quad \text{with} \quad \tilde{\mathbf{F}} = (\mathbf{I} - \nabla_2 \mathbf{d})^{-T}, \quad \nabla_2 \mathbf{d} = \begin{pmatrix} \partial_1 d_1 & \partial_1 d_2 \\ \partial_2 d_1 & \partial_2 d_2 \end{pmatrix}.
$$

Consequently tr$_{(F^T F)_j}$ = tr$_{(\tilde{F}^T \tilde{F})_j}$ + 1, j = 1, 2, (tr$_{F^T F}$)$_2$ = (tr$_{\tilde{F}^T \tilde{F}}$)$_2$ + 2tr$_{\tilde{F}^T \tilde{F}}$ + 1,

$$
\Psi(\mathbf{F}) = c_1 \mathrm{tr}_{\mathbf{F}^T \mathbf{F}} + c_2 (\mathrm{tr}_{(\mathbf{F}^T \mathbf{F})^2} - \mathrm{tr}_{\mathbf{F}^T \mathbf{F}}^2) = (c_1 - 2c_2) \mathrm{tr}_{\tilde{\mathbf{F}}^T \tilde{\mathbf{F}}} + c_2 (\mathrm{tr}_{(\tilde{\mathbf{F}}^T \tilde{\mathbf{F}})^2} - \mathrm{tr}_{\tilde{\mathbf{F}}^T \tilde{\mathbf{F}}}^2) + \text{constant.} \quad (11)
$$

Hence a comparison between the 2D model and the 3D model on a 2D configuration requires
replacing the $c_1$ of the 2D Helmholtz potential by $c_1 - 2c_2$ from the 3D Helmholtz potential.

2.6. Variational Formulation

For simplicity, we shall consider only the case of homogeneous boundary conditions on Γ ⊂ ∂Ω$^t$, i.e., clamped or no-slip, and homogeneous Neumann conditions on ∂Ω$^t$ \ Γ.
So, given $\Omega_f^0, \Omega_s^0, \mathbf{d}$, and $\mathbf{u}$ at $t=0$, we must find $(\mathbf{u}, p, \mathbf{d}, \Omega_f^t, \Omega_s^t)$ with $\mathbf{u}|_{\Gamma} = 0$ and,

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \int_{\Omega_f^t} (\rho \mathbb{D}_t \mathbf{u} \cdot \hat{\mathbf{u}} - p \nabla \cdot \hat{\mathbf{u}} - \hat{p} \nabla \cdot \mathbf{u}) + \int_{\Omega_f^t} \frac{\mu}{2} \mathbf{D} \mathbf{u} : \mathbf{D} \hat{\mathbf{u}} \\
\displaystyle + \int_{\Omega_s^t} \left[ c_1 (\mathbf{D} \mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d})^2 + c_3 (\mathbf{D} \mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) \right] : \mathbf{D} \hat{\mathbf{u}} = \int_{\Omega_f^t} \mathbf{f} \cdot \hat{\mathbf{u}}, \\
\mathbb{D}_t \mathbf{d} = \mathbf{u},
\end{array}
\right. 
\quad (12) 
$$

for all $(\hat{\mathbf{u}}, \hat{p})$ with $\hat{\mathbf{u}}|_{\Gamma} = 0$ and where $\Omega_s^t$ and $\Omega_f^t$ are defined by:

$$ \frac{d\chi}{d\tau} = \mathbf{u}(\chi(\tau), \tau), \quad \chi(t) \in \Omega_r^t \quad \Rightarrow \quad \chi(\tau) \in \Omega_r^\tau \quad \forall \tau \in (0, T), \quad r=s, f. $$

### 3. Numerical Schemes

#### 3.1. Characteristic-Galerkin Derivatives

The characteristic-Galerkin method is applied to (12) to discretize the total derivatives.

Let $\mathbf{u}: \Omega^t \times (0, T) \to \mathbb{R}^3$ be given, Lipschitz in space and continuous in time. Let $\chi_{\mathbf{u},x}^t(\tau)$ be the solution at time $\tau < t$ of:

$$ \dot{\chi}(\tau) = \mathbf{u}(\chi(\tau), \tau) \quad \text{with} \quad \chi(t) = x. $$

So $\chi$ is the particle path in a velocity field $\mathbf{u}$ and $\chi_{\mathbf{u},x}^t(\tau)$ is the position of a particle at time $\tau$ which will be at position $x$ at (future) time $t$. The method of characteristics relies on the concept of total derivative: for any differentiable function $\mathbf{w}: \Omega^t \times (0, T) \to \mathbb{R}^m, m = 1, 2, 3, \dots$,

$$ \mathbb{D}_t \mathbf{w}(x, t) := \frac{d}{d\tau} \mathbf{w}(\chi(\tau), \tau) \Big|_{\tau=t} = \partial_t \mathbf{w} + \mathbf{u} \cdot \nabla \mathbf{w}. \quad (13) $$

Given a time step $\delta t$, the position of a particle at time $n\delta t$, which will be at position $x$ at time $(n+1)\delta t$, is approximated by $\mathbb{Y}^{n+1}(x)$:

$$ \chi_{\mathbf{u}^{n+1},x}^{(n+1)\delta t}(n\delta t) \approx \mathbb{Y}^{n+1}(x) := x - \mathbf{u}^{n+1}(x) \delta t. $$

Then a first-order consistent scheme for the total derivative of $\mathbf{w}$ is:

$$ \partial_t \mathbf{w} + \mathbf{u} \cdot \nabla \mathbf{w} = \mathbb{D}_t \mathbf{w}(x, t) \approx \frac{1}{\delta t} (\mathbf{w}^{n+1}(x) - \mathbf{w}^n(\mathbb{Y}^{n+1}(x))). \quad (14) $$

Combined with a Galerkin approximation in space the scheme has been analyzed in depth in [21] in the variational context of the finite element methods.

#### 3.2. A Monolithic Time-Discrete Variational Formulation

$$ \text{Define } L_0^2(\Omega) = \{q \in L^2(\Omega) : \int_{\Omega} q = 0\}. $$

**Proposition 1.** The following variational problem is a first order in time-consistent approximation of (12): Find, $\Omega^{n+1}, \mathbf{u}^{n+1} \in \mathbf{H}_0^1(\Omega^{n+1})$, $p^{n+1} \in L_0^2(\Omega^{n+1})$ such that for all $\hat{\mathbf{u}} \in \mathbf{H}_0^1(\Omega^{n+1})$, $\hat{p} \in L_0^2(\Omega^{n+1})$,
$$
\left\{
\begin{aligned}
& \int_{\Omega^{n+1}} \left( \rho^{n+1} \frac{\mathbf{u}^{n+1} - \mathbf{u}^n \circ \mathbb{Y}^{n+1}}{\delta t} \cdot \hat{\mathbf{u}} - p^{n+1} \nabla \cdot \hat{\mathbf{u}} - \hat{\rho} \nabla \cdot \mathbf{u}^{n+1} \right) + \int_{\Omega_f^{n+1}} \frac{\mu}{2} \mathbf{D} \mathbf{u}^{n+1} : \mathbf{D} \hat{\mathbf{u}} \\
& \qquad + \delta t \int_{\Omega_s^{n+1}} \left[ 2c_1 \left( \mathbf{D} \mathbf{u}^{n+1} - \nabla \mathbf{u}^{n+1} \nabla^T \tilde{\mathbf{d}}^n - \nabla \tilde{\mathbf{d}}^n \nabla^T \mathbf{u}^{n+1} \right) \left( \mathbf{D} \tilde{\mathbf{d}}^n - \nabla \tilde{\mathbf{d}}^n \nabla^T \tilde{\mathbf{d}}^n \right) \right. \\
& \qquad \left. + c_3^n \left( \mathbf{D} \mathbf{u}^{n+1} - \nabla \mathbf{u}^{n+1} \nabla^T \tilde{\mathbf{d}}^n - \nabla \tilde{\mathbf{d}}^n \nabla^T \mathbf{u}^{n+1} \right) \right] : \mathbf{D} \hat{\mathbf{u}} \\
& \qquad + \int_{\Omega_s^{n+1}} \left[ c_1 \left( \mathbf{D} \tilde{\mathbf{d}}^n - \nabla \tilde{\mathbf{d}}^n \nabla^T \tilde{\mathbf{d}}^n \right)^2 + c_3^n \left( \mathbf{D} \tilde{\mathbf{d}}^n - \nabla \tilde{\mathbf{d}}^n \nabla^T \tilde{\mathbf{d}}^n \right) \right] : \mathbf{D} \hat{\mathbf{u}} = \int_{\Omega_t} \mathbf{f} \cdot \hat{\mathbf{u}}, \\
& \mathbf{d}^{n+1} = \mathbf{d}^n \circ \mathbb{Y}^{n+1} + \delta t \mathbf{u}^{n+1}.
\end{aligned}
\right.
\quad (15)
$$

where $\tilde{\mathbf{d}}^n$ stands for $\mathbf{d}^n(\mathbb{Y}^{n+1})$.

## Some comments

1. Formula (14) has been used to approximate $D_t \mathbf{u}$ and $D_t \mathbf{d}$ in (12).

2. One may wonder why the scheme is applied to $\mathbf{u}$ and not to $\rho \mathbf{u}$. Note that $\rho = \rho_f \mathbf{1}_{\Omega_f^t} + \rho_s \mathbf{1}_{\Omega_s^t}$ is convected by the velocity $\mathbf{u}$. Hence $\rho^{n+1}(x) = \rho^n \circ \mathbb{Y}^{n+1}(x)$. This shows that discretizing the total derivative of $\mathbf{u}$ or the total derivative of $\rho \mathbf{u}$ gives the same scheme:

$$ \frac{1}{\delta t} \left( \mathbf{w}^{n+1}(x) - \mathbf{w}(\mathbb{Y}^{n+1}(x)) \right) = (\partial_t \mathbf{w} + \mathbf{u} \cdot \nabla \mathbf{w})|_{x,t^{n+1}} + O(\delta t) \text{ with } \mathbf{w} = \mathbf{u} \text{ or } \mathbf{w} = \rho \mathbf{u}. \quad (16) $$

**Proof.** With the characteristic Galerkin method, a consistent time discretization of (12) would be such that at each time step one must:

Find $\mathbf{u}^{n+1} \in \mathbf{H}_0^1(\Omega^{n+1})$, $p \in L_0^2(\Omega^{n+1})$, $\Omega^{n+1} = \Omega_f^{n+1} \cup \Omega_s^{n+1}$, such that $\forall \hat{\mathbf{u}} \in \mathbf{H}_0^1(\Omega^{n+1})$, $\forall \hat{p} \in L_0^2(\Omega^{n+1})$; the three relations below hold:

$$ \int_{\Omega^{n+1}} \left( \rho^{n+1} \frac{\mathbf{u}^{n+1} - \mathbf{u}^n \circ \mathbb{Y}^{n+1}}{\delta t} \cdot \hat{\mathbf{u}} - p^{n+1} \nabla \cdot \hat{\mathbf{u}} - \hat{\rho} \nabla \cdot \mathbf{u} \right) + \int_{\Omega_f^{n+1}} \frac{\mu}{2} \mathbf{D} \mathbf{u}^{n+1} : \mathbf{D} \hat{\mathbf{u}} \\ + \int_{\Omega_s^{n+1}} \left[ c_1 \left( \mathbf{D} \mathbf{d}^{n+1} - \nabla \mathbf{d}^{n+1} \nabla^T \mathbf{d}^{n+1} \right)^2 + c_3^{n+1} \left( \mathbf{D} \mathbf{d}^{n+1} - \nabla \mathbf{d}^{n+1} \nabla^T \mathbf{d}^{n+1} \right) \right] : \mathbf{D} \hat{\mathbf{u}} = \int_{\Omega_t} \mathbf{f} \cdot \hat{\mathbf{u}}, \quad (17) $$

$$ \mathbf{d}^{n+1} = \tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1}, \quad \text{where } \tilde{\mathbf{d}}^n = \mathbf{d}^n \circ \mathbb{Y}^{n+1} \quad (18) $$

$$ \Omega^{n+1} = (\mathbb{Y}^{n+1})^{-1}(\Omega^n) = \{x : \mathbb{Y}^{n+1}(x) := x - \delta t \mathbf{u}^{n+1}(x) \in \Omega^n\}. \quad (19) $$

A fully implicit monolithic formulation with variables $(\mathbf{u}^{n+1}, p^{n+1})$ can be derived by substituting $\mathbf{d}^{n+1}$ in (17) with (18):

$$ \int_{\Omega^{n+1}} \left( \rho^{n+1} \frac{\mathbf{u}^{n+1} - \mathbf{u}^n \circ \mathbb{Y}^{n+1}}{\delta t} \cdot \hat{\mathbf{u}} - p^{n+1} \nabla \cdot \hat{\mathbf{u}} - \hat{\rho} \nabla \cdot \mathbf{u}^{n+1} \right) + \int_{\Omega_f^{n+1}} \frac{\mu}{2} \mathbf{D} \mathbf{u}^{n+1} : \mathbf{D} \hat{\mathbf{u}} \\ + \int_{\Omega_s^{n+1}} \left[ c_1 \left( \mathbf{D} \left( \tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1} \right) - \nabla \left( \tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1} \right) \nabla^T \left( \tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1} \right) \right)^2 \right. \\ \left. + c_3^{n+1} \left( \mathbf{D} \left( \tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1} \right) - \nabla \left( \tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1} \right) \nabla^T \left( \tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1} \right) \right) \right] : \mathbf{D} \hat{\mathbf{u}} = \int_{\Omega_t} \mathbf{f} \cdot \hat{\mathbf{u}}. \quad (20) $$

System (15) is found by expanding the nonlinear terms in (20), and keeping only the zero order terms and the terms of order one in $\delta t$. $\square$
### 3.3. Spacial Discretization with Finite Elements

Let $T_h^0$ be a triangulation of the initial domain. Spatial discretization can be done with the most popular finite element for fluids: the Lagrangian triangular elements of degree 2 for the space $V_h$ of velocities and displacements and Lagrangian triangular elements of degree 1 for the pressure space $Q_h$, provided that the pressure is different in the structure and the fluid because the pressure is discontinuous at the interface $\sigma$; therefore $Q_h$ is the space of piecewise linear functions on the triangulations, continuous in $\Omega_h^{n+1}$, $r = s$, $f$. Appropriate couples like $\{V_h, Q_h\}$ are chosen to satisfy the inf-sup condition to avoid checker-board oscillations (see for example [8]). A small penalization with parameter $\epsilon$ must be added to impose uniqueness of the pressure when $L_h^2$ is replaced by $Q_h \approx L^2$.

This leads us to find $\mathbf{u}_h^{n+1} \in V_{0h}$, $p_h^{n+1} \in Q_h$, $\Omega^{n+1}$, such that for all $\hat{\mathbf{u}}_h$, $\hat{p}_h \in V_{0h} \times Q_h$ with,

$$ \tilde{\mathbf{d}}_h^n := \mathbf{d}_h^n \circ \mathbb{Y}^{n+1}, \quad \text{where} \quad \mathbb{Y}^{n+1}(x) = x - \mathbf{u}_h^{n+1}(x) \delta t, $$

the following holds:

$$ \left\{ \begin{aligned} & \int_{\Omega^{n+1}} \left( \rho^{n+1} \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n \circ \mathbb{Y}^{n+1}}{\delta t} \cdot \hat{\mathbf{u}}_h - p_h^{n+1} \nabla \cdot \hat{\mathbf{u}}_h - \hat{p}_h \nabla \cdot \mathbf{u}_h^{n+1} + \epsilon \nabla p_h^{n+1} \nabla \hat{p}_h \right) \\ & + \int_{\Omega_f^{n+1}} \frac{\mu}{2} \mathbf{D} \mathbf{u}_h^{n+1} : \mathbf{D} \hat{\mathbf{u}}_h \\ & + \delta t \int_{\Omega_s^{n+1}} \left[ 2c_1 (\mathbf{D} \mathbf{u}_h^{n+1} - \nabla \mathbf{u}_h^{n+1} \nabla^T \tilde{\mathbf{d}}_h^n - \nabla \tilde{\mathbf{d}}_h^n \nabla^T \mathbf{u}_h^{n+1}) (\mathbf{D} \tilde{\mathbf{d}}_h^n - \nabla \tilde{\mathbf{d}}_h^n \nabla^T \tilde{\mathbf{d}}_h^n) \right. \\ & \qquad \left. + c_3^{n+1} (\mathbf{D} \mathbf{u}_h^{n+1} - \nabla \mathbf{u}_h^{n+1} \nabla^T \tilde{\mathbf{d}}_h^n - \nabla \tilde{\mathbf{d}}_h^n \nabla^T \mathbf{u}_h^{n+1}) \right] : \mathbf{D} \hat{\mathbf{u}}_h \\ & + \int_{\Omega_s^{n+1}} \left[ c_1 (\mathbf{D} \tilde{\mathbf{d}}_h^n - \nabla \tilde{\mathbf{d}}_h^n \nabla^T \tilde{\mathbf{d}}_h^n)^2 + c_3^{n+1} (\mathbf{D} \tilde{\mathbf{d}}_h^n - \nabla \tilde{\mathbf{d}}_h^n \nabla^T \tilde{\mathbf{d}}_h^n) \right] : \mathbf{D} \hat{\mathbf{u}}_h = \int_{\Omega^{n+1}} \mathbf{f} \cdot \hat{\mathbf{u}}_h \\ & \Omega^{n+1} = (\mathbb{Y}^{n+1})^{-1} (\Omega^n) = \{x: \mathbb{Y}^{n+1}(x) \in \Omega^n\}. \\ & \mathbf{d}_h^{n+1} = \tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1}. \end{aligned} \right. \quad (21) $$

### 3.4. Solution Algorithm

Equation (21) must be solved iteratively because $\Omega^{n+1}$ is updated by (19). The most natural method is to freeze some coefficients, $c_3$ included, so as to obtain a well-posed linear problem and iterate.

To explain the iteration process let us rename, in (19), $\rho^{n+1}$, $\Omega^{n+1}$, $\Omega_f^{n+1}$, $\Omega_s^{n+1}$, $\mathbb{Y}^{n+1}$, $c_3^{n+1}$, into $\Omega$, $\Omega^f$, $\Omega^s$, $\mathbb{Y}$, $c_3$.

1. Set $\rho = \rho_n$, $c_3 = c_3^n$, $\Omega^r = \Omega_r^n$, $r = f, s$, $\mathbf{u} = \mathbf{u}_h^n$, $\mathbb{Y}(x) = x - \mathbf{u}\delta t$,

2. Solve system (21). In this study a direct solver is used for the linear system.

3. Set $\mathbf{u} = \mathbf{u}_h^{n+1}$, $\mathbb{Y}(x) = x - \mathbf{u}\delta t$, $\Omega^r = \mathbb{Y}^{-1}(\Omega_r^n)$, $r = s, f$; update $c_3$ and $\rho$.

4. If not converged return to step 2.

**Remark 1.** $\Omega_s^{n+1}$ we move the vertices $q^j$ in the solid part by its velocity $\mathbf{u}(q^j)$. As $\tilde{\mathbf{d}}_h^n$ is moved by $\mathbb{Y}$, it is naturally convected on the new mesh. Similarly, $\rho$ being constant on $\Omega^s$, it is well defined on the new mesh. Finally, $c_3$ being defined in terms of $\mathbf{u}$, it is also well defined on the new mesh.

In the fluid region we generate a new mesh with a Delauney-Voronoi automatic mesh generator from the new position of the interface $\sigma$. There are many other possibilities like solving a Laplace equation in the fluid part with Dirichlet conditions equal to the velocity of the fluid-solid interface. It works equally well but what we advocate is, in our experience, more robust, even if it is at the cost of a possibly greater interpolation error at every step because $\tilde{\mathbf{u}}_h^n$ needs to be projected on the new mesh. A $P^2$ polynomial interpolation is used.
**Remark 2.** An earlier study in [17] showed that 2 iterations are necessary but 3 or more make invisible improvements; this result is confirmed in 3D by the numerical tests below.

By inspection of (23) below, notice that this iterative algorithm computes the stationary point of $\frac{d}{dt} \int_{\Omega_s^0} \Psi (\mathbf{I} + \nabla_{x_0} \mathbf{d}^T)$. As $\Psi$ is convex, it should be a convergent process, but it requires a more rigorous proof.

# 4. Stability of the Scheme

## 4.1. Conservation of Energy

**Proposition 2.** In the continuous case, it holds that:

$$ \frac{d}{dt} \int_{\Omega^t} \frac{\rho}{2} |\mathbf{u}|^2 + \frac{\mu^f}{2} \int_{\Omega_f^t} |\mathbf{D}\mathbf{u}|^2 + \frac{d}{dt} \int_{\Omega_s^0} \Psi (\mathbf{I} + \nabla_{x_0} \mathbf{d}^T) = \int_{\Omega^t} \mathbf{f} \cdot \mathbf{u}. \quad (22) $$

**Proof.** Choosing $\hat{\mathbf{u}} = \mathbf{u}$, $\hat{\rho} = -p$ in Equation (12), then one term becomes (9):

$$ \int_{\Omega_s^t} \left[ c_1 (\mathbf{D}\mathbf{d} - \nabla\mathbf{d}\nabla^T\mathbf{d})^2 + c_3 (\mathbf{D}\mathbf{d} - \nabla\mathbf{d}\nabla^T\mathbf{d}) \right] : \nabla\mathbf{u} = \int_{\Omega_s^t} (\partial_{\mathbf{F}}\Psi\mathbf{F}^T - \alpha\mathbf{I}) : \nabla\mathbf{u}. \quad (23) $$

Now, $\frac{d}{dt}\Psi(\mathbf{F}) = \partial_{\mathbf{F}}\Psi(\mathbf{F}) : \partial_t \mathbf{F}$ and $\mathbf{F}^T\nabla\mathbf{u} = \nabla_{x_0}\mathbf{u}(x_0) = \partial_t\nabla_{x_0}\mathbf{d}(x_0) = \partial_t\mathbf{F}^T(x_0)$, so conservation of the Helmholtz potential $\Psi$ can be derived as follows:

$$ \begin{aligned} \int_{\Omega_s^t} (\partial_{\mathbf{F}}\Psi\mathbf{F}^T - \alpha\mathbf{I}) : \nabla\mathbf{u} &= \int_{\Omega_s^0} \partial_{\mathbf{F}}\Psi : \nabla_{x_0}\mathbf{u} = \int_{\Omega_s^0} \partial_{\mathbf{F}}\Psi : \partial_t\nabla_{x_0}\mathbf{d} = \int_{\Omega_s^0} \partial_{\mathbf{F}}\Psi : \partial_t\mathbf{F}^T \\ &= \int_{\Omega_s^0} \frac{d}{dt}\Psi(\mathbf{F}) = \frac{d}{dt}\int_{\Omega_s^0} \Psi(\mathbf{F}) = \frac{d}{dt}\int_{\Omega_s^0} \Psi(\mathbf{I} + \nabla_{x_0}\mathbf{d}^T). \end{aligned} $$

The other terms in (12) are standard, in particular, with enough regularity on $\Omega^t$,

$$ \int_{\Omega^t} (\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{u} = \frac{d}{dt} \int_{\Omega^t} \frac{|\mathbf{u}|^2}{2}. $$

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

**Remark 3.** When $\Psi$ is convex, some regularity can be gained from this energy conservation which can lead to existence of solution up to time $T$ which, loosely speaking, in the most optimistic perspective, is the first time of contact of two boundaries which were not in contact initially (see [22–24] etc.).

## 4.2. Stability of the Scheme Discretized in Time

**Lemma 1.**

$$ \mathbf{F}^{n+1} = \mathbf{F}^n + \delta t \nabla_{x^0} \mathbf{u}^{n+1}. $$

**Proof.** By construction, $\mathbb{X}^{n+1} = (\mathbb{Y}^1 \circ \mathbb{Y}^2 \dots \mathbb{Y}^{n+1})^{-1} = (\mathbb{Y}^{n+1})^{-1} \circ \mathbb{X}^n$; therefore $\mathbf{d}_h^{n+1} = \tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1}$ implies that:

$$ \nabla_{x^0} \mathbf{d}^{n+1} = \nabla_{x^0} \mathbf{d}^n + \delta t \nabla_{x^0} \mathbf{u}^{n+1}. $$

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

**Theorem 1.** When (17) is used and if $\mathbf{f} = 0$ and $\rho_r$ is constant in each domain $\Omega_r^n$, $r = s, f$,

$$ \int_{\Omega^{n+1}} \frac{\rho^{n+1}}{2} |\mathbf{u}^{n+1}|^2 + \delta t \int_{\Omega_f^{n+1}} \frac{\mu}{2} |\mathbf{D}\mathbf{u}^{n+1}|^2 + \int_{\Omega_s^0} \Psi(\mathbf{F}^{n+1}) \le \int_{\Omega^n} \frac{\rho^n}{2} |\mathbf{u}^n|^2 + \int_{\Omega_s^0} \Psi(\mathbf{F}^n) \quad (24) $$
**Proof.** The proof is similar to that given in [19] for the 2D case.

With $\hat{u} = u^{n+1}$, $\hat{p} = -p^{n+1}$, and $f = 0$, (17) becomes:

$$ \int_{\Omega^{n+1}} \left( \rho^{n+1} \frac{u^{n+1} - u^n \circ \mathbb{Y}^{n+1}}{\delta t} \cdot u^{n+1} \right) + \frac{\mu}{2} \int_{\Omega_f^{n+1}} |\mathbf{D}u^{n+1}|^2 + \int_{\Omega_s^{n+1}} \left[ c_1 (\mathbf{D}\mathbf{d}^{n+1} - \nabla\mathbf{d}^{n+1}\nabla^T\mathbf{d}^{n+1})^2 + c_3^{n+1} (\mathbf{D}\mathbf{d}^{n+1} - \nabla\mathbf{d}^{n+1}\nabla^T\mathbf{d}^{n+1}) \right] : \mathbf{D}u^{n+1} = 0 \quad (25) $$

We note that $\rho_r^n \circ \mathbb{Y}^{n+1}(x) = \rho_r^{n+1}(x)$, $x \in \Omega_r^{n+1}$; so:

$$ \int_{\Omega_r^{n+1}} (\rho_r^n \mathbf{u}^n) \circ \mathbb{Y}^{n+1} \cdot \mathbf{u}^{n+1} = \int_{\Omega_r^{n+1}} (\sqrt{\rho_r^n} \mathbf{u}^n) \circ \mathbb{Y}^{n+1} \cdot \sqrt{\rho_r^{n+1}} \mathbf{u}^{n+1} $$

Schwartz inequality applied to the right side followed by Young's inequality leads to:

$$ \int_{\Omega_r^{n+1}} (\rho_r^n \mathbf{u}^n) \circ \mathbb{Y}^{n+1} \cdot \mathbf{u}^{n+1} \le \left( \int_{\Omega_r^{n+1}} (\sqrt{\rho_r^n} \mathbf{u}^n)^2 \circ \mathbb{Y}^{n+1} \right)^{\frac{1}{2}} \left( \int_{\Omega_r^{n+1}} \rho_r^{n+1} (\mathbf{u}^{n+1})^2 \right)^{\frac{1}{2}} \\ = \left( \int_{\Omega_r^{n+1}} (\sqrt{\rho_r^n} \mathbf{u}^n)^2 \right)^{\frac{1}{2}} \left( \int_{\Omega_r^{n+1}} \rho_r^{n+1} (\mathbf{u}^{n+1})^2 \right)^{\frac{1}{2}} \le \frac{1}{2} \int_{\Omega_r^{n+1}} \rho_r^{n+1} |\mathbf{u}^{n+1}|^2 + \frac{1}{2} \int_{\Omega_r^{n+1}} \rho_r^n |\mathbf{u}^n|^2 $$

Because $(\mathbb{Y}^{n+1})^{-1}(\Omega^{n+1}) = \Omega^n$, and $\det_{\nabla\mathbb{Y}^{n+1}} = 1$, for any $f$,

$$ \int_{\Omega^{n+1}} f (\mathbb{Y}^{n+1}(x)) dx = \int_{(\mathbb{Y}^{n+1})^{-1}(\Omega^{n+1})} f(y) \det_{(\nabla\mathbb{Y}^{n+1})^{-1}} dy = \int_{\Omega^n} f. \quad (26) $$

The third integral in (25) can be transformed with Equation (9):

$$ \int_{\Omega_s^{n+1}} \left[ c_1 (\mathbf{D}\mathbf{d}^{n+1} - \nabla\mathbf{d}^{n+1}\nabla^T\mathbf{d}^{n+1})^2 + c_3^{n+1} (\mathbf{D}\mathbf{d}^{n+1} - \nabla\mathbf{d}^{n+1}\nabla^T\mathbf{d}^{n+1}) \right] : \mathbf{D}u^{n+1} \quad (27) \\ = \int_{\Omega_s^{n+1}} \partial_{\mathbf{F}} \Psi \mathbf{F}^T |^{n+1} : \nabla \mathbf{u}^{n+1} = \frac{1}{\delta t} \int_{\Omega_s^0} \partial_{\mathbf{F}} \Psi (\mathbf{F}^{n+1}) : (\mathbf{F}^{n+1} - \mathbf{F}^n) $$

The last equality is due to Lemma 1. Now, by the convexity of $\Psi$:

$$ \int_{\Omega_s^0} \partial_{\mathbf{F}} \Psi (\mathbf{F}^{n+1}) : (\mathbf{F}^{n+1} - \mathbf{F}^n) \ge \int_{\Omega_s^0} \Psi (\mathbf{F}^{n+1}) - \int_{\Omega_s^0} \Psi (\mathbf{F}^n) $$

Finally,

$$ \int_{\Omega^{n+1}} \frac{\rho^{n+1}}{2} |\mathbf{u}^{n+1}|^2 + \delta t \int_{\Omega_f^{n+1}} \frac{\mu}{2} |\mathbf{D}\mathbf{u}^{n+1}|^2 + \int_{\Omega_s^0} \Psi (\mathbf{F}^{n+1}) \le \int_{\Omega^n} \frac{\rho^n}{2} |\mathbf{u}^n|^2 + \int_{\Omega_s^0} \Psi (\mathbf{F}^n) \quad (28) $$

□

### 4.3. Energy Inequality for the Fully Discrete Scheme

It was argued in [18] that the same proof holds for scheme discretized with $P^1$ triangular elements for velocity and pressure because then $\mathbb{Y}$ transforms a tetrahedron into a tetrahedron and so (27) holds also in the discrete case. Compatible linear elements are either the $P^1 - P^1$ stabilized element or the $P^1$ element for the pressure with $P^1$ also for the velocity but on a mesh where each tetrahedron is divided in sub-tetrahedra from an inner additional vertex. In the first case the above stability analysis needs to be modified to include the stabilization term. Numerical tests in 2D showed that the $P^2 - P^1$ element is alright even though energy stability cannot be proved.
Decrease of energy (viscous terms included) and conservation of mass are shown on Figure 2a. Energy is computed using $mesh_1$ and 1, 2, 3 or 4 iterations for the non-linear solver of Section 3.4; the results show that 2 iterations for the nonlinear solver is optimal. Energy is also plotted when $mesh_3$ and $mesh_5$ are used (with 2 iterations). In agreement with the theory, energy decays with time even on a coarse mesh. The decay is due to the numerical dissipation because for the continuous problem it should be constant, and indeed it is less with finer meshes. Note that an impulsive start is not an energy-friendly test case.

Conservation of mass is shown on Figure 2b for $mesh_1$, $mesh_3$ and $mesh_5$. It is within 2%.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a): Energy vs. time computed on $mesh_1$ with 1, 2, 3 and 4 iterations in the solver of Section 3.4. On the same plot energy vs. time is also plotted for $mesh_3$ and $mesh_5$. (b): Mass conservation vs. time when $mesh_1$, (lower curve), $mesh_3$ and $mesh_5$ are used. Curves 3 and 5 are overlapping.</figcaption></figure>

To study convergence we take the results obtained with $mesh_7$ as a reference “almost exact” solution. Then, we observe the point on the torus with maximum radius (i.e., farthest to the right). On Figure 3a the maximum radius is shown versus time for the four meshes; on Figure 3b the pointwise error for $mesh_i$, $i = 1, 3, 5$. Convergence is observed; we make no claim on the speed of convergence versus mesh size because $mesh_7$ is probably too coarse to be reliable.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a): Evolution of the farthest point on the right of the torus versus time computed with $mesh_i$, $i = 1$ (oscillations), 3 and 5. (b): $Mesh_7$ is used as a reference computations. The errors on the farthest point on the right in the torus versus time are shown for $mesh_1$ (highest curve, its mean absolute value is 0.013), $mesh_3$ (lowest curve, its mean absolute value is 0.0023), and $mesh_5$ (flat curve, its mean absolute value is $4.1 \times 10^{-5}$).</figcaption></figure>
## 5.2. Clamped Structure in a Fluid

A set of computations of a beam in vacuum or in fluid has been performed at the request of Larma's PhD advisor, Miguel Fernandez [26]. Experimental measurements are available [20] so it is intended to make them into a well-documented benchmark for the community. The geometry is always a clamped beam in a container; in all but the first case there is a flow in the container but to initialize the computation the position of the beam is computed with a stagnant fluid.

The beam in Hessenthaler's experiment and Larma's simulation is not incompressible, as the Poisson ratio is 0.35 instead of 0.5, so we expect discrepancies. The present method can generalized to compressible material and 2D simulations show similar stability properties but the computer program is not yet ready for 3D cases.

### 5.2.1. Free-Fall of a Clamped Structure in Vacuum

The beam size is 9 × 1 × 1. The beam is clamped on the right to the yz plane and free to fall under its own weight. The density $ρ_s$ is 1, and the gravity is -0.01 along y direction. The Mooney-Rivlin coefficients are $c_1 = 5$ and $c_2 = 2.5$. Computation stops at $T = 100$ after 3200 steps with $δt = 0.03125$. Results are shown at different times in Figure 4.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Free-falling beam in vacuum; snapshots at time (a) $t = 10$, (b) $t = 20$, (c) $t = 40$, (d) $t = 70$, (e) $t = 80$, (f) $t = 90$.</figcaption></figure>

At the end of computation, the error of volume is less than 1%. A discussion about $c_3$ is found later in Section 5.3.3.

### 5.2.2. Free-Fall of a Clamped Beam in a Fluid

Second, we consider a flat beam of size 9 × 1 × 4 with one side attached to a box which is filled with a fluid. The box size is 10 × 7 × 5. Solid density $ρ_s$ is 1. Fluid density $ρ_f$ is 0.5, and viscosity $μ_f$ equals 0.1. The gravity $g$ is -0.5 along the y direction. The Mooney-Rivlin coefficients are $c_1 = 16.67$ and $c_3 = 8.33 (\text{tr}_B^2 - \text{tr}_B^2 - 4)$. Computation stops at $T = 100$ after 200 steps with $δt = 0.5$.

Computing results are shown in Figure 5; volume error is less than 0.5% at the end.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Free-falling beam in fluid; snapshots at time (a) t = 25, (b) t = 50, (c) t = 75, (d) t = 100.</figcaption></figure>

5.3. *The Benchmark*

An experiment was designed by Hessenthaler et al. [20] to calibrate the computer simulations of Larma [26].

The geometry is a cylindrical chamber with a length of 200 mm, diameter 76.2 mm and axis
parallel to x = y = 0. Two inlet circular wholes are on the left z = 0 wall, one with center at
(0, -27.15, 0), the other at (0, 27.15, 0), both of diameter 21.9 mm; the inlet pipes are 20 mm long. A
silicon filament clamped on the z = 0 plane is 2 mm thick, 11 mm wide, and 65 mm long; its center
is located at (0, 0, 0). The computational domain is discretized into 10⁴ vertices and adapted to the
silicone filament to decrease error during its motion. Geometric schematic is shown in Figure 6.

As in [20], the density of the silicone filament is ρ_s = 1.0583 × 10⁻³ g/mm³. The incompressible Mooney–Rivlin hyperelastic coefficients are determined by curve-fitting to uni-axial tensile-load displacement test data. In [26] two coefficients are introduced: Ĉ₁ = 103,533.82/2 Pa and Ĉ₂ = 8891.65/4 Pa. The Poisson ratio is ν = 0.315 in [26] and c₁ = Ĉ₁/(2(1+ν)), c₂ = Ĉ₂/(2(1+ν)). There will be a difference because incompressibility implies ν = 0.5. Gravity g is −9810 mm/s² along the y-direction. Two numerical tests are made, a steady phase I, and a transient phase II.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Geometric schematic. Blue: fluid domain; white: solid domain.</figcaption></figure>

5.3.1. Phase I

In the fluid $\rho_f = 1.1633 \times 10^{-3} \text{ g/mm}^{-3}$ and $\mu_f = 1.25 \times 10^{-2} \text{ Pa}\cdot\text{s}$. A Poiseuille profile flow is applied to the upper and lower pipes with peak velocity of 630 mm/s and 615 mm/s, respectively. The position of center line of the beam along the z direction is recorded as a function of time in Figure 7 and compared with experimental data.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Computational results for phase I (a) position of center line of the beam along z direction (results and experiments overlap), (b) velocity norms about the symmetry plane.</figcaption></figure>

In absence of flow, due to buoyancy forces, the silicone filament reaches an hydrostatic equilibrium with maximum displacement $d_M^0$, at the tip on the right, approximately equal to 29.50 mm in [20] and 31.04 mm in [26]. To compute this displacement we have taken $T = 1$ and 1000 steps with $\delta t = 0.001$.

Then $d^0$ is taken as initial condition to Phase I as the inflows are turned on; the maximum displacements are now approximately 16.41 mm in [20] and 15.99 mm in [26]. To compute this state we have taken $T = 2, 2000$ steps with $\delta t = 0.001$. Our results, shown in Figure 7, overlap the experimental data.
### 5.3.2. Phase II

Finally a transient computation is performed. The fluid properties are $ρ_f = 1.1640 \times 10^{-3}$ g/mm$^{-3}$ and $μ = 1.337 \times 10^{-2}$ Pa·s.

The Poiseuille profile flow applied on the upper pipe has a transient peak velocity $\hat{v}_k(t)$, $k \in [x, y, z]$ shown in Figure 8 and corresponds to a smooth interpolation from Table 1. No flow is applied to the lower inflow pipe.

Computation stops at $T = 5$ after 1600 steps with $\delta t = 0.003125$. The maximum velocity norm of the solid is less than 30 mm/s; the error introduced by the computation is less than 0.1 mm. Deflection of the center line of the silicone filament is recorded versus time and a comparison between numerical results and experimental data is shown in Figure 9 at six different times. These results display a fair—yet not satisfactory—agreement between the simulations and the experimental data. The difference is perhaps due to improper determination of the coefficients of the hyperelastic model, because bending at the end of silicone filament could be observed only when the magnitude of $c_2$ is large enough. Another possible explanation is the wrong Poisson ratio, namely an effect due to incompressibility.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Experimental curves of the peak inflow of inflow velocities used for the boundary condition in the computations.</figcaption></figure>

<table><caption>**Table 1.** For phase II, curve-fitting coefficients of inlet peak velocity for $\hat{v}_k(t) = \sigma_{i=1}^3 n_i t^i / \sigma_{j=0}^4 b_j t^j$ with $\hat{v}_k = 0$ for $t \in I \setminus I_k$. Note that the flow in y direction is applied only in the upper inlet.</caption><thead><tr><th>θ̂</th><th>n<sub>1</sub></th><th>n<sub>2</sub></th><th>n<sub>3</sub></th><th>d<sub>0</sub></th><th>b<sub>1</sub></th><th>b<sub>2</sub></th><th>b<sub>3</sub></th><th>b<sub>4</sub></th><th>Î<sub>k</sub></th></tr></thead><tbody><tr><td>θ̂<sub>x</sub></td><td>-11.37</td><td>-28.99</td><td>7.73</td><td>1.38</td><td>0.24</td><td>3.59</td><td>-3.14</td><td>1</td><td>[0, 4.07]</td></tr><tr><td>θ̂<sub>y</sub></td><td>14.95</td><td>11.88</td><td>-2.17</td><td>2.06</td><td>-2.0</td><td>4.95</td><td>-3.50</td><td>1</td><td>[0, 5.51]</td></tr><tr><td>θ̂<sub>z</sub></td><td>367.10</td><td>363.40</td><td>-62.24</td><td>1.21</td><td>-0.38</td><td>3.76</td><td>-3.19</td><td>1</td><td>[0, 5.27]</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Deflection of the silicone filament in phase II (a) $t = 0.721$, (b) $t = 1.153$, (c) $t = 1.585$, (d) $t = 2.017$, (e) $t = 2.449$, (f) $t = 2.881$.</figcaption></figure>

5.3.3. Variation of Coefficient $c_3$

The coefficient $c_3$ in (10) varies with the structural rheology through $tr_B$ and $tr_{B^2}$. We wish to analyze the non-linearity of $c_3 = \frac{c_1}{2} (tr_B^2 - tr_{B^2} - 4) - 2c_2$ during the simulation.

Maximum and minimum magnitude of $c_3$ and relative volume error with respect to time are plotted in Figure 10a, and the corresponding maximum magnitude of displacement **d** is shown in Figure 10b.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Phase II results (a) maximum and minimum magnitude of $2c_3$ and relative volume error multiplied by 1000, (b) maximum magnitude of displacement **d**.</figcaption></figure>

It is interesting to point out that $c_3$ is not very sensitive to the displacement gradients in this simulation as shown in Figure 11: $c_3$ remains fairly constant during the computation, even though **d** does not. These results seem reasonable since deflection over the whole silicone filament is fairly
uniform. However, we do not rule out that strong deformations could cause a stability problem in the algorithm; this will be investigated in the future.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Deflection **d** = (d₁, d₂, d₃) and a₂ = 2c₃ of the silicone filament in phase II at (a) t = 0.721, (b) t = 1.153, (c) t = 1.585, (d) t = 2.017, (e) t = 2.449, (f) t = 2.881.</figcaption></figure>

**6. Conclusions**

In this study, a fully Eulerian fluid-structure formulation has been presented. An implicit unconditionally stable monolithic finite element scheme based on characteristic-Galerkin discretization has been proposed and studied. At each time step a non-linear system must be solved; however semi-linearization leads to well-posed linear monolithic systems for the velocities and pressures.

What was found in [17,18] for 2D systems could be generalized in 3D with similar stability based on analytical proof that the energy of the discrete system cannot increase.

The numerical method has been implemented with FreeFem++ [25], and validation tests have been presented. For an axisymmetric case, conservation of mass and energy were shown as well as convergence on a sequence of embedded meshes.

The effect of nonlinearity of c₃—the major difference with the 2D case—is studied in Section 5.3.3 with validation tests.

In the near future we intend to parallelize the computer program and use it for larger systems, in particular for hemodynamics.

**Acknowledgments:** This research is partially supported by the Ministry of Education of the Republic of China (R.O.C Taiwan) for the PhD scholarship of Chen-Yu Chiang.

**Author Contributions:** Chen-Yu Chiang contributed to the derivation of the equations in 3D, wrote the 3D-FreeFem++ program and did all the 3D tests. Olivier Pironneau contributed to the derivation of the equations in 3D, to the iterative algorithm and to the energy estimate; he is also responsible for the axisymmetric example and the 2D-axisymmetric FreeFem program. Tony W. H. Sheu and Marc Thiriet are the PhD advisors of Chen-Yu Chiang. They are responsible for quality check of this study and for the future application of this method to hemodynamics which are currently developed in parallel, using the same computer programs.

**Conflicts of Interest:** The authors declare no conflict of interest.
References

1. Martins, J. Wing design via numerical optimization. *SIAG OPT Views News* **2015**, *23*, 2–7.

2. Peskin, C.S. The immersed boundary method. *Acta Numer.* **2002**, *11*, 479–517.

3. Formaggia, L.; Quarteroni, A.; Veneziani, A.; (Eds.). *Cardiovascular Mathematics*; Springer: Milano, Italy, 2009.

4. Patrice Hauret. Numerical Methods for the Dynamic Analysis of Two-Scale Incompressible Nonlinear Structures. Ph.D. Thesis, Ecole Polytechnique, Paris, France, 2004.

5. Hron, J.; Turek, S. A monolithic FEM solver for an ALE formulation of fluid-structure interaction with configuration for numerical benchmarking. In Proceedings of the European Conference on Computational Fluid Dynamics, Egmond aan Zee, The Netherlands, 5–8 September 2006.

6. Basting, S.; Quaini, A.; Čanić, S.; Glowinski, R. Extended ALE Method for fluid-structure interaction problems with large structural displacements. *J. Comput. Phys.* **2017**, *331*, 312–336.

7. Liu, J. A second-order changing-connectivity ALE scheme and its application to FSI with large convection of fluids and near contact of structures. *J. Comput. Phys.* **2016**, *304*, 380–423.

8. Boffi, D.; Cavallini, N.; Gastaldi, L. The finite element immersed boundary method with distributed Lagrange multiplier. *SIAM J. Numer. Anal.* **2015**, *53*, 2584–2604.

9. Wang, Y.; Jimack, P.K.; Walkley, M.A. A one-field monolithic fictitious domain method for fluid-structure interactions. *Comput. Methods Appl. Mech. Eng.* **2017**, *317*, 1146–1168.

10. Maury, B.; Glowinski, R. Fluid-particle flow: A symmetric formulation. *C. R. Acad. Sci. Ser. I Math.* **1997**, *324*, 1079–1084.

11. Glowinski, R.; Pan, T.W.; Hesla, T.I.; Joseph, D.D.; Périaux, J. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: Application to particulate flow. *J. Comput. Phys.* **2001**, *169*, 363–426.

12. Liu, C.; Walkington, N.J. An eulerian description of fluids containing visco-elastic particles. *Arch. Ration. Mech. Anal.* **2001**, *159*, 229–252.

13. Rannacher, R.; Richter, T. An Adaptive Finite Element Method for Fluid-Structure Interaction Problems Based on a Fully Eulerian Formulation. In *Fluid Structure Interaction II*; Springer: Berlin/Heidelberg, Germany, 2011; pp. 159–191.

14. Dunne, T.; Rannacher, R. Adaptive Finite Element Approximation of Fluid-Structure Interaction Based on an Eulerian Variational Formulation. In *Fluid-Structure Interaction*; Springer: Berlin/Heidelberg, Germany, 2006; pp. 110–145.

15. Richter, T. A Fully Eulerian formulation for fluid-structure-interaction problems. *J. Comput. Phys.* **2013**, *233*, 227–240.

16. Wick, T. Fully Eulerian fluid-structure interaction for time-dependent problems. *Comput. Methods Appl. Mech. Eng.* **2013**, *255*, 14–26.

17. Pironneau, O. Numerical Study of a Monolithic Fluid-Structure Formulation. In *Variational Analysis and Aerospace Engineering*; Mohammadi, B., Frediani, A., Cipolla, V., Eds.; Number 116 in Springer Optimization and Its Applications; Springer International Publishing: Cham, Switzerland, 2016; pp. 401–420.

18. Hecht, F.; Pironneau, O. An energy stable monolithic eulerian fluid-structure finite element method. *Int. J. Numer. Methods Fluid* **2017**, doi:10.1002/fld.4388.

19. Pironneau, O. An energy preserving monolithic eulerian fluid-structure numerical scheme. *arXiv* **2016**, arXiv: 1607.08083.

20. Hessenthaler, A.; Gaddum, N.R.; Holub, O.; Sinkus, R.; Röhrle, O.; Nordsletten, D.A. Experiment for validation of fluid-structure interaction models and algorithms. *Int. J. Numer. Methods Biomed. Eng.* **2016**, doi:10.1002/cnm.2848.

21. Pironneau, O. On the transport-diffusion algorithm and its applications to the Navier-Stokes equations. *Numer. Math.* **1982**, *38*, 309–332.

22. Boulakia, M. Existence of weak solutions for the motion of an elastic structure in an incompressible viscous fluid. *C. R. Math.* **2003**, *336*, 985–990.

23. Coutand, D.; Shkoller, S. Motion of an Elastic Solid inside an Incompressible Viscous Fluid. *Arch. Ration. Mech. Anal.* **2005**, *176*, 25–102.

24. Raymond, J.P.; Vanninathan, M. A fluid-structure model coupling the Navier-Stokes equations and the Lamé system. *J. Math. Pures Appl.* **2014**, *3*, 546–596.
25. Hecht, F. New development in freefem++. *J. Numer. Math.* **2013**, *20*, 251–266.

26. Larma, M.L. Coupling Schemes and Unfitted Mesh Methods for Fluid-Structure Interaction. Ph.D. Thesis, Université Pierre et Marie Curie (Paris VI), Paris, France, 2016.

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

© 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).