# An Energy Stable Monolithic Eulerian Fluid-Structure Finite Element Method

Frédéric Hecht, Olivier Pironneau

► To cite this version:

Frédéric Hecht, Olivier Pironneau. An Energy Stable Monolithic Eulerian Fluid-Structure Finite Element Method. International Journal for Numerical Methods in Fluids, 2017, pp.1-17. 10.1002/fld.4388. hal-01326785v2

HAL Id: hal-01326785

https://hal.sorbonne-universite.fr/hal-01326785v2

Submitted on 10 Jul 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
F. Hecht, O. Pironneau, An Energy Stable Monolithic FSI-FEM

Sorbonne Universités, UPMC (Paris VI) Laboratoire Jacques-Louis Lions fred-
eric.hecht@upmc.fr, olivier.pironneau@upmc.fr

Working paper submitted to IJNMF

Abstract

When written in an Eulerian frame, the conservation laws of contin-
uum mechanics are similar for fluids and solids leading to a single set of
variables for a monolithic formulation. Such formulations are well adapted
to large displacement fluid-structure configurations, but stability is a chal-
lenging problem because of moving geometries. In this article the method
is presented; time implicit discretizations are proposed with iterative algo-
rithms well posed at each step, at least for small displacements; stability
is discussed for an implicit in time finite element method in space by
showing that energy decreases with time. The key numerical ingredient
is the Characteristics-Galerkin method coupled with a powerful mesh gen-
erator. A numerical section discusses implementation issues and presents
a few simple tests. It is also shown that contacts are easily handled by
extending the method to variational inequalities. This paper deals only
with incompressible neo-Hookean Mooney-Rivlin hyperelastic material in
two dimensions in a Newtonian fluid; but the method is not limited to
these; compressible and 3D cases will be presented later.

Finite Element, Fluid-Structure, Energy Stable, Monolithic scheme, Characteristics-Galerkin, Eulerian Formulation.
# An Energy Stable Monolithic Eulerian Fluid-Structure Finite Element Method

Frédéric Hecht and Olivier Pironneau

July 10, 2017

## Introduction

For the numerical simulation of fluid-structure interactions, arbitrary Lagrangian Eulerian methods (ALE) are popular and efficient in the case of small displacements (see [32] and the references therein). For large displacements the difficulty is transferred to the mesh [27] [37] and to a lesser extent to the matching conditions at the fluid-solid interface [25]. Furthermore, iterative solvers which rely on alternative solutions in the fluid part and in the structure part are subject to the added mass effect and require special preconditioners [17][8]. Cited references here are certainly not exhaustive as an enormous amount of work has gone into ALE methods.

Immersed boundary methods (IBM) [28] are very efficient for thin structures. For solid shells or membranes in 3D fluids, precision seems to require level sets and for 3D solids mesh adaption [12], yet it is a very flexible method which has been shown stable – for a given deformation map – by Boffi et al [5].

The Eulerian approach makes use of the fact that continuum mechanics does not distinguish between solids and fluids up to the constitutive equations. Classically, one uses the Lagrangian variable, denoted $\mathbf{x} = \mathbf{X}(\mathbf{x}_0, t)$, to restate the equations for the structures in terms of the $x_0$-gradient of its displacement $\mathbf{d}(\mathbf{x}_0, t) := \mathbf{X}(\mathbf{x}_0, t) - \mathbf{x}_0$ defined in the initial domain $\mathbf{x}_0 \in \Omega_0^s$. In an Eulerian formulation one frames the solids by its velocity $\mathbf{u}(\mathbf{x}, t)$ and displacement $\mathbf{d}(\mathbf{x}, t) = \mathbf{x} - \mathbf{X}^{-1}(\mathbf{x}, t)$ in the moving domain $\mathbf{x} \in \Omega_t^s$ with $\mathbf{x}$-gradients of these; here $\mathbf{X}^{-1}$ is the inverse map of $\mathbf{x}_0 \mapsto \mathbf{x} = \mathbf{X}(\mathbf{x}_0, t)$.

These facts have been highlighted particularly in [26] (and perhaps also in [3]) and even more so by Stefan Turek in [22] (yet still using $d(\mathbf{x}_0, t)$ and ALE) who revived the term “monolithic”; we use it here in this sense: one variational equation for the whole system both at the continuous and discrete levels by opposition to a partition or segregated approach.

In a series of papers in the wake of Thomas Dunne’s thesis (2007) (see [14][15]), ALE was replaced by a fully Eulerian algorithm for the structure and the fluid. As more papers have been written on this topic, in particular by Rolf
Rannacher[16] (hierarchical meshes)[33] (adaptive meshes), Thomas Wick[37] (comparison with ALE), Thomas Richter[35][36] ((an XFEM-like strategy to capture the interface and also an immersed boundary-like strategy and fixed meshes) and Stefan Frei[19] (detailed study of contacts) we will refer to them as of the Heidelberg school.

While embarking on Eulerian formulations, in [30], we were not aware of the work of the Heidelberg school and so while rediscovering the formulation which appeared first in Dunne’s thesis we devised different discretization and solution strategies.

More precisely while Dunne works with velocities in the fluid and displacement $d(x, t)$ in the solid, we propose to work with velocities everywhere as in [36]. There are minor differences in the continuous formulation (no inverse matrix of $\nabla d$) and essential differences in the solver at each time step and in the spatial discretizations; while an immerse boundary fixed mesh method is used in [35] and an XFEM-like method in [36], we update the solid region with a Lagrangian mesh; it is a departure from the fully Eulerian paradigm but it is only for the mesh. One reason for this choice is that for moving domains it is hard to find a better mesh adaptation than the one given by the solid deformation itself; another reason is that we can prove energy stability with it.

Several differences with the previous work of the Heidelberg school are highlighted here:

1. In our case every step requires the solution of a well posed problem, at least in the case of small displacements.

2. We use the Characteristics-Galerkin method to solve the equation which connects the velocity to the displacement in Eulerian coordinates:

3. We propose an implicit in time adaptation of the fluid-structure interface, which requires fixed-point type iterations at each time step but which allow us to establish an energy estimate, even in the fully discrete case.

4. Finally we move the vertices of the mesh in the structure with its own velocity $u$ while we remesh at every time step in the fluid.

This Eulerian monolithic approach is presented in two dimensions for Mooney-Rivlin neo-Hookean incompressible hyperelastic materials coupled to incompressible Newtonian flows but it extends to compressible materials [31] and to three dimensions [11].

The paper may be difficult to follow by readers used to work with structures modeled by the equations of elasticity in the initial domain. So a summary of the Eulerian approach is given before embarking on a more detailed derivation of the equations.

In sections 1 and 2 we recall the fundamental principles and the algebra that lead to the Eulerian formulation (5) for a fluid-structure problem. Unsurprisingly this monolithic formulation is in term of velocity-pressure $u, p$ for the
fluid part and displacement-pressure **d**, *p* for the solid part; however it requires
a little rethinking to write the structure part in the moving domain $\Omega_t^s$.

A key point in section 3 for the fully implicit in time discretization of system
(5) is the approximation of the purely convective space-time PDEs by finite-
difference in time approximations of the material derivatives. In variational form
this is known as the Characteristics-Galerkin method [29]. Henceforth $\mathbf{d}^{n+1}$ can
be eliminated in favor of $\mathbf{u}^{n+1}$ so as to obtain a monolithic FSI system which
involves a solver in terms of $\mathbf{u}^{n+1}$ and $p^{n+1}$ only, and in the full domain; this is
(7).

Now there is a delicate problem with PDEs written on time dependent do-
mains: $\mathbf{d}^{n+1}(\mathbf{x})$ (resp. $\mathbf{d}^n(\mathbf{x})$) makes sense only if $\mathbf{x} \in \Omega_{n+1}^s$ (resp. $\mathbf{x} \in \Omega_n^s$).
This forces the following definition for $\Omega_{n+1}^s$: it is the set of $\mathbf{x}$ such that
$\mathbf{x} - \mathbf{u}^{n+1}(\mathbf{x})\delta t \in \Omega_n^s$. This is both a difficulty and an advantage. It requires
an iterative solution of (7) with a readjustment of $\Omega_{n+1}^s$ – we propose in (9) a
simple fixed point iteration scheme– but it is also an asset because, evidently,
$\Omega_{n+1}^s$ is an unknown and an implicit scheme should iterate on it.

Although we cannot prove that the fixed point scheme converges, we can
prove (see section 4) that at each step the problem is well posed for small
displacements and that it is a first order consistent time approximation of the
continuous problem. Moreover we can prove that the scheme is energy stable
(Theorem 2 of section 5).

In section 6, we propose a spatial approximation by the finite element
method. Because of the fluid part – and probably also for the solid part– the
velocity and pressure spaces have to be compatible with the inf-sup condition.
One natural choice is the $P^2 - P^1$ triangular conforming element, yet one has to
allow for discontinuity in the pressure at the fluid-structure interface. Another
well known element is the $P^1$ bubble–$P^1$ element. One great feature of this last
choice is that with a minor modification it is energy stable too. To achieve this
goal the cubic bubble is replaced by a linear bubble constructed by dividing
the pressure triangles into 3 subtriangles using any point in pressure triangles
as fourth vertex. We refer below to it as the $P_3^1 - P^1$ element. Finally, the
simplest is to use a stabilized $P^1$ – $P^1$ element on the same grid for velocities
and pressure; this leads also to an energy stable scheme.

Adapting the mesh to moving domains is difficult and in fact the Heidelberg
school has several of their aforementioned papers devoted to this problem. We
propose a boundary fitted strategy (by opposition to IBM): move the solid mesh
with its own velocity and remesh the fluid part by a Delaunay-Voronoi algorithm
from the knowledge of the new position of the interface boundary. It requires
a projection step in the fluid (most other options we know of do too); here it
is done by linear or quadratic interpolation, depending on the degree of the
element. Such an algorithm is compatible with the energy stability result, at
least in the solid part; in the fluid part it would require an energy diminishing
projection operator in place of interpolation, but our numerical tests do not
point to such a complicate necessity.

Finally, in section 7, we show briefly that contacts modeled by variational
inequalities are fairly easy to handle with this Eulerian monolithic scheme, espe-
cially if a semi-smooth Newton method is used for solving the nonlinear system.

Now, besides energy stability, an important part of this article is devoted to
the numerical implementation, with the following aim:

1. Validate the theoretical propositions: in all our numerical tests we found the method to be quite robust and and capable of giving decent results even with coarse meshes. Hence it is very fast; all tests are done on a Macbook Pro in minutes, except the flag behind a cylinder in subsection 8.2 which requires half an hour.

2. Study the stability of the fixed point iterations: 2 iterations is an optimal choice.

3. Study error estimates, energy and mass conservation.

4. Compare with results available in the literature.

Errors are not as good as expected from the degree of the polynomial approx-
imation. In a way this is expected because at the fluid-structure interface the
gradient of the velocities are discontinuous. This is shown in subsection 8.1 on a
configuration for which we have a semi-analytical solution with a discontinuous
gradient.

The conservation of energy is studied in subsection 8.3. The results are
very good but they do not show that the $P_3^1 - P^1$ element is better that the
$P^1$bubble-$P^1$ element. They show however that two iterations are much better
than one in the fixed point algorithm (it was shown in [30] that 3 iterations or
more make no improvement). For the conservation of mass the results are also
within 1% after 50 time steps.

Comparison with earlier numerical results is made difficult by the fact that
most studies except one by Dunne [33] have been done for compressible solids.
We have performed two cases proposed by Dunne and one by Frei [19] for
contact. Agreement with Dunne’s (subsection 8.2) is not good in one case but
good in the other and excellent with Frei’s (subsection 8.4), even though Frei’s
is for a compressible solid; probably the solid is hard enough and hence almost
incompressible.

# 1 Continuum Mechanics: notations

Consider two non overlapping time dependent bounded sets of $\mathbb{R}^2$, $\Omega_t^f$ and $\Omega_t^s$ for the fluid and solid respectively, $t \in (0, T)$. Let the full domain be $\Omega_t = \Omega_t^f \cup \Omega_t^s$. The fluid-structure interface is denoted $\Sigma_t = \Omega_t^f \cap \Omega_t^s$ and the boundary of $\Omega_t$ is $\partial\Omega_t$. Later we will denote $\Gamma$ the part of $\partial\Omega_t$ where either the structure is clamped or the fluid does not slip (see Figure 1) At initial time $\Omega_0^f$ and $\Omega_0^s$ are prescribed.

The following notations are standard in continuum mechanics (e.g. [24],[3],
[1]):
<figure><img src="image_2.png" /><figcaption>Figure 1: Sketch of the fluid - structure domain and notations</figcaption></figure>

- $\mathbf{X}: \Omega_0 \times (0, T) \mapsto \mathbf{X}(\mathbf{x}_0, t) \in \Omega_t$, the Lagrangian position at $t$ of $\mathbf{x}_0$; in other words, a material point at $\mathbf{x}_0$ at $t=0$ is transported to $\mathbf{X}$ at time $t$ by the motion.

- $\mathbf{u} = \partial_t \mathbf{X}$, the velocity of the motion,

- $\mathbf{F} = \nabla_{\mathbf{x}_0}^T \mathbf{X} := (\nabla_{\mathbf{x}_0} \mathbf{X})^T = ((\partial_{\mathbf{x}_0 j} \mathbf{X}_i))$, the transposed gradient of the motion. We follow [4]; however in some books on continuum mechanics the gradient is the contravariant gradient and the transposed operator is absent.

- $\mathbf{d} = \mathbf{X}(\mathbf{x}_0, t) - \mathbf{x}_0$, the displacement vector of the structure.

We denote by tr$_A$ and det$_A$ the trace and determinant of A. To describe the system, let $\rho(\mathbf{x}, t)$ be the density and $\sigma(\mathbf{x}, t)$ the stress tensor. at position $\mathbf{x}$ and time $t$. Finally and unless specified otherwise, all spatial derivatives are with respect to $\mathbf{x} \in \Omega_t$ and not with respect to $\mathbf{x}_0 \in \Omega_0$.

We recall a few properties:

1. If $\phi$ is a function of $\mathbf{x} = \mathbf{X}(\mathbf{x}_0, t)$, $\mathbf{x}_0 \in \Omega_0$,

$$ \nabla_{\mathbf{x}_0} \phi = ((\partial_{\mathbf{x}_0 i} \phi)) = \sum_j ((\partial_{\mathbf{x}_0 i} \mathbf{X}_j \partial_{\mathbf{x}_j} \phi)) = \mathbf{F}^T \nabla \phi. $$

2. When $\mathbf{X}$ is one-to-one and invertible, $\mathbf{d}$ and $\mathbf{F}$ can be seen as functions of $(\mathbf{x}, t)$ instead of $(\mathbf{x}_0, t)$. They are related by

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

3. Time derivatives are related by (notice the notation for the material derivative)

$$ \mathbb{D}_t \phi := \frac{d}{dt} \phi(\mathbf{X}(\mathbf{x}_0, t), t) = \partial_t \phi(\mathbf{x}, t) + \mathbf{u} \cdot \nabla \phi(\mathbf{x}, t). $$

Finally, it is convenient to introduce the notations:

$$ J = \det_{\mathbf{F}}, \quad \mathbf{D}\mathbf{u} = \nabla\mathbf{u} + \nabla^T\mathbf{u}. $$
Conservation of momentum and conservation of mass take the same form for the fluid and the solid:

$$ \rho \mathbb{D}_t \mathbf{u} = \mathbf{f} + \nabla \cdot \sigma, \quad \frac{d}{dt}(J\rho) = 0, \qquad (2) $$

where $\mathbf{f}$ is the density of volumic forces. So $J\rho = \rho_0$ at all times and, consequently, if the medium is incompressible, $\rho$ is constant in the fluid and constant in the solid if it was constant at $t = 0$. We denote these constants by $\rho^f$ and $\rho^s$ and at any point $\mathbf{x}$ and time $t$, the density is $\rho(\mathbf{x}, t) = \rho^f \mathbf{1}_{\Omega_t^f}(\mathbf{x}, t) + \rho^s \mathbf{1}_{\Omega_t^s}(\mathbf{x}, t)$, where $\mathbf{1}_{\Omega_t^f}(\mathbf{x}) = 1$ if $\mathbf{x} \in \Omega$ and $0$ otherwise.

Continuity of $\mathbf{u}$ and of $\sigma \cdot \mathbf{n}$ at the fluid-structure interface $\Sigma$ (with $\mathbf{n}$ its normal) is part of (2), in absence of external surface force. There are also unwritten constraints pertaining to the realizability of the map $\mathbf{X}$ (see [10],[24]), for instance $\mathbf{X}$ must be invertible.

## 1.1 Constitutive Equations

In this work we shall consider only two dimensional hyperelastic incompressible Mooney-Rivlin materials – also referred as neo-Hookean – and Newtonian incompressible viscous fluids; for those the stress tensors are

$$ \sigma^f := \sigma_{|\Omega_t^f} = -p^f \mathbf{I} + \mu \mathbf{D} \mathbf{u}, \quad \sigma^s := \sigma_{|\Omega_t^s} = -p^s \mathbf{I} + \partial_{\mathbf{F}} \Psi \mathbf{F}^T \qquad (3) $$

The Helmholtz potential $\Psi$ for a Mooney-Rivlin material is given by

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

where $c_1, c_2$ are constant (see [10]). Consequently $\partial_{\mathbf{F}} \Psi \mathbf{F}^T = (2c_1 - 4c_2 \mathrm{tr}_{\mathbf{F}^T \mathbf{F}}) \mathbf{F} \mathbf{F}^T + 4c_2 (\mathbf{F} \mathbf{F}^T)^2$.

In 2d, with any symmetric non singular matrix $B$, one can express $B^{-1}$ and $B^2$ in terms of $B$ and the identity $\mathbf{I}$ by the Cayley-Hamilton theorem. Now by making use of (1) one obtains – for some scalars $\alpha, \alpha'$, functions of $\mathbf{x}$ and $t$,

$$ \partial_{\mathbf{F}} \Psi \mathbf{F}^T = 2c_1 (\mathbf{I} - \nabla \mathbf{d})^{-T} (\mathbf{I} - \nabla \mathbf{d})^{-1} + \alpha \mathbf{I} = 2c_1 (\mathbf{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) + \alpha' \mathbf{I}. \qquad (4) $$

We now have a complete set of equations (2,3,4) in which the initial domain $\Omega_0^s$ does not appear.

# 2 Eulerian Variational Formulation in 2D

Written in variational form system (2,3,4) is:

Given $\Omega_0^f, \Omega_0^s$ and $\mathbf{d}, \mathbf{u}$ at $t = 0$, find $(\mathbf{u}, p, \mathbf{d}, \Omega_t^f, \Omega_t^s)$ with $\mathbf{u}|_{\Gamma} = 0$ and

$$ \int_{\Omega_t^f \cup \Omega_t^s} [\rho \mathbb{D}_t \mathbf{u} \cdot \hat{\mathbf{u}} - p \nabla \cdot \hat{\mathbf{u}} - \hat{p} \nabla \cdot \mathbf{u}] + \int_{\Omega_t^f} \frac{\mu}{2} \mathbf{D} \mathbf{u} : \mathbf{D} \hat{\mathbf{u}} $$
$$ + \int_{\Omega_t^s} c_1(\mathbf{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) : \mathbf{D}\hat{\mathbf{u}} = \int_{\Omega_t^f \cup \Omega_t^s} \mathbf{f} \cdot \hat{\mathbf{u}}, \quad (5) $$

for all $(\hat{\mathbf{u}}, \hat{p})$ with $\hat{u}|_{\Gamma} = 0$, where $\mathbf{d}$, $\Omega_t^s$ and $\Omega_t^f$ are defined forward in time by $\mathbf{D}_t \mathbf{d} = \mathbf{u}$ and

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

We have used the notation $\mathbf{B} : \mathbf{C} = \text{tr}_{\mathbf{B}} \mathbf{r}_{\mathbf{C}}$. Relation (6) defines $\Omega_t^s, \Omega_t^f$ forward in time.

For simplicity we shall consider only the case of homogeneous boundary conditions on $\Gamma \subset \partial\Omega$, i.e. clamped or no-slip, and homogeneous Neumann conditions on $\partial\Omega_t \setminus \Gamma$.

**Remark 1** A similar formulation has been used by Dunne in [15], eq.(9) (see also [14]), but for compressible materials ans using the first form in (4). To solve the problem at each time step they use Newton iterations on the whole system. By contrast we will eliminate $\mathbf{d}$ from the time discretized system at each time step.

# 3 Discretization in Time

## 3.1 A Monolithic Time-Discrete Variational Formulation

**Theorem 1** The following variational problem is a first order in time consistent approximation of (5): find $\mathbf{u}^{n+1}, p^{n+1}, \Omega_{n+1}$ with

$$ \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{p} \nabla \cdot \mathbf{u}^{n+1} \right] + \int_{\Omega_{n+1}^f} \frac{\mu}{2} \mathbf{D}\mathbf{u}^{n+1} : \mathbf{D}\hat{\mathbf{u}} \\ + c_1 \int_{\Omega_{n+1}^s} [D(\tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1}) - \nabla(\tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1}) \nabla^T(\tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1})] : \mathbf{D}\hat{\mathbf{u}} = \int_{\Omega_{n+1}} \mathbf{f} \cdot \hat{\mathbf{u}}, \\ \Omega_{n+1} = \{ \mathbf{x} : \mathbb{Y}^{n+1}(\mathbf{x}) := \mathbf{x} - \delta t \mathbf{u}^{n+1}(\mathbf{x}) \in \Omega_n \} \quad (7) $$

where the $\sim$ stands for the composition with $\mathbb{Y}^{n+1}$, i.e. $\tilde{\mathbf{d}}^n = \mathbf{d}^n \circ \mathbb{Y}^{n+1}$, and where $\mathbf{d}$ is updated by

$$ \mathbf{d}^{n+1} = \tilde{\mathbf{d}}^n + \delta t \mathbf{u}^{n+1}. \quad (8) $$

We defer the proof to a dedicated section, after showing how the problem can be solved numerically.

## 3.2 Iterative Solution by Fixed Point

In (5), the bilinear part on $\Omega_t^f$ is identical to that of a time dependent Stokes problem. The integral in the solid region, on the other hand, is a trilinear form; as the term of degree 3 is $O(\delta t^2)$, we neglect it and solve the nonlinearity due to $\mathbb{Y}$ and $\Omega^f, \Omega^s$, by an iterative process which is essentially a fixed point algorithm:
1. Set $\rho = \rho_n$, $\Omega = \Omega_n$, $\mathbf{u} = \mathbf{u}^n$, $\mathbb{Y}(\mathbf{x}) = \mathbf{x} - \delta t \mathbf{u}$.

2. Solve, for all $\hat{\mathbf{u}}, \hat{p}$ with $\hat{\mathbf{u}}|_{\Gamma} = 0$,

$$ \int_{\Omega} [\rho \frac{\mathbf{u}^{n+1} - \tilde{\mathbf{u}}^n}{\delta t} \cdot \hat{\mathbf{u}} - p^{n+1} \nabla \cdot \hat{\mathbf{u}} - \hat{p} \nabla \cdot \mathbf{u}^{n+1}] + \int_{\Omega^f} \frac{\mu}{2} \mathrm{D} \mathbf{u}^{n+1} : \mathrm{D} \hat{\mathbf{u}} \\ + \int_{\Omega^s} \delta t c_1 (\mathrm{D} \mathbf{u}^{n+1} - \nabla \tilde{\mathbf{d}}^n \nabla^T \mathbf{u}^{n+1} - \nabla \mathbf{u}^{n+1} \nabla^T \tilde{\mathbf{d}}^n) : \mathrm{D} \hat{\mathbf{u}} \\ = \int_{\Omega} \mathbf{f} \cdot \hat{\mathbf{u}} - \int_{\Omega^s} c_1 (\mathrm{D} \tilde{\mathbf{d}}^n - \nabla \tilde{\mathbf{d}}^n \nabla^T \tilde{\mathbf{d}}^n) : \mathrm{D} \hat{\mathbf{u}} \quad (9) $$

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

4. If not converged return to Step 2.

**Remark 2** It will be seen below that, after discretization in space, setting $\mathbb{Y}(\mathbf{x}) = \mathbf{x} - \delta t \mathbf{u}$, $\Omega^r = \mathbb{Y}^{-1}(\Omega_n^r)$, $r = s, f$ means to move each vertex $\mathbf{q}$ of the triangulation of $\Omega^s$ to $\mathbf{q} + \delta t \mathbf{u}(\mathbf{d})$ and adapt $\Omega^f$ to be compatible with $\Omega^s$, for instance by remeshing entirely from the knowledge of the new position of the interface.

**Remark 3** Notice that (9) contains a bilinear form on $\mathbf{u}^{n+1}$ in $\Omega_{n+1}^s$ which may not be coercive. So we do not know that the problem is well posed, nor can we prove that the iterations converge. However for small displacements $\mathrm{D}\mathbf{u}$ will dominate $\nabla\mathbf{u}\nabla\mathbf{d}$ and so the problem will be well posed for sufficiently small $\nabla\mathbf{d}$.

# 4 Proof of Theorem 1

To prove theorem 1 we need first to recall the Characteristics-Galerkin method and apply it to (5).

## 4.1 Discretization of Material Derivatives

Given $t \in (0, T)$ and $\mathbf{x} \in \Omega$, let $\chi_{\mathbf{u},\mathbf{x}}^t(\tau)$ be the solution at time $\tau$ of

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

If $\mathbf{u}$ is Lipschitz in space and continuous in time the solution exists. The Characteristics-Galerkin method relies on the concept of material derivative:

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

Given a time step $\delta t$, let us approximate $\chi$ by $\mathbb{Y}$:

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

**Remark 4** Note that

$$ \mathbb{J}_{n+1} := \det_{\nabla \mathbb{Y}^{n+1}} = 1 - \nabla \cdot \mathbf{u}^{n+1} \delta t + \det_{\nabla \mathbf{u}^{n+1}} \delta t^2. \quad (10) $$

So, in the incompressible case, one can neglect the Jacobians and still be $O(\delta t)$-consistent.
## 4.2 Variational Form of the Time Discrete Problem

Thus a consistent time discretization of (5) would be that at each time step one must:

Find $\mathbf{u}^{n+1} \in \mathbf{H}_0^1(\Omega_{n+1})$, $p \in L^2(\Omega_{n+1})$, $\Omega_{n+1}^f$, $\Omega_{n+1}^s$, such that with $\Omega_{n+1} = \Omega_{n+1}^f \cup \Omega_{n+1}^s$, $\forall \hat{\mathbf{u}} \in \mathbf{H}_0^1(\Omega_{n+1})$, $\forall \hat{p} \in L^2(\Omega_{n+1})$ the 3 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{p} \nabla \cdot \mathbf{u}^{n+1} \right] + \int_{\Omega_{n+1}^f} \frac{\mu}{2} \mathrm{D}\mathbf{u}^{n+1} : \mathrm{D}\hat{\mathbf{u}} + c_1 \int_{\Omega_{n+1}^s} [\mathrm{D}\mathbf{d}^{n+1} - \nabla \mathbf{d}^{n+1} \nabla^T \mathbf{d}^{n+1}] : \mathrm{D}\hat{\mathbf{u}} = \int_{\Omega_{n+1}} \mathbf{f} \cdot \hat{\mathbf{u}}, \quad (11) $$

$$ \Omega_{n+1} = (\mathbb{Y}^{n+1})^{-1}(\Omega_n) = \{\mathbf{x} : \mathbb{Y}^{n+1}(\mathbf{x}) := \mathbf{x} - \delta t \mathbf{u}^{n+1}(\mathbf{x}) \in \Omega_n\} \quad (12) $$

$$ \mathbf{d}^{n+1} = \mathbf{d}^n \circ \mathbb{Y}^{n+1} + \delta t \mathbf{u}^{n+1}. \quad (13) $$

**Remark 5** Note that $\rho = \rho^f \mathbf{1}_{\Omega_t^f} + \rho^s \mathbf{1}_{\Omega_t^s}$ is convected by the velocity $\mathbf{u}$. Denoting $\rho_n(\mathbf{x}) := \rho(\mathbf{x}, n\delta t)$,

$$ \rho_n \circ \mathbb{Y}^{n+1}(\mathbf{x}) = \rho_{n+1}(\mathbf{x}), \quad \mathbf{x} \in \Omega_{n+1} $$

So discretizing the material derivative of $\mathbf{u}$ or the material derivative $\rho \mathbf{u}$ gives the same scheme because in both cases

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

◇

Now to solve (11,12,13), let us plug (13) into (11) to obtain the “monolithic” formulation of Theorem 1, i.e. using,

$$ \nabla \mathbf{d}^{n+1} = \nabla \tilde{\mathbf{d}}^n + \delta t \nabla \mathbf{u}^{n+1} $$

Thus (7) is proved.

# 5 Stability by Energy estimation

## 5.1 Conservation of Energy in the Continuous Case

The following is standard; it is given here as an aid for the proof of the same for the discretized system.

**Proposition 1**

$$ \frac{d}{dt} \int_{\Omega_t} \frac{\rho}{2} |\mathbf{u}|^2 + \frac{\mu}{2} \int_{\Omega_t^f} |\mathrm{D}\mathbf{u}|^2 + \frac{d}{dt} \int_{\Omega_0^s} \Psi(\mathbf{I} + \nabla_{\mathbf{x}_0}^T \mathbf{d}) = \int_{\Omega_t} \mathbf{f} \cdot \mathbf{u} $$
Proof. Choosing $\hat{\mathbf{u}} = \mathbf{u}$, $\hat{p} = -p$ will give the proposition provided

$$2c_1 \int_{\Omega_t^s} (\mathrm{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) : \mathrm{D}\mathbb{D}_t \mathbf{d} = \frac{d}{dt} \int_{\Omega_0^s} \Psi(\nabla_{\mathbf{x}_0} \mathbf{X}).$$

By construction

$$\int_{\Omega_t^s} c_1(\mathrm{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) : \mathrm{D}\hat{\mathbf{u}} = \int_{\Omega_t^s} (\partial_{\mathbf{F}} \Psi(\mathbf{F}) \mathbf{F}^T - \alpha \mathbf{I}) : \mathrm{D}\hat{\mathbf{u}} = \int_{\Omega_0^s} \partial_{\mathbf{F}} \Psi(\mathbf{F}) : \mathrm{D}_{\mathbf{x}_0} \hat{\mathbf{u}}$$

Now as $\frac{d}{dt} \Psi(\mathbf{F}) = \partial_{\mathbf{F}} \Psi(\mathbf{F}) : \partial_t \mathbf{F}$ and $\nabla_{\mathbf{x}_0} \mathbf{u}(\mathbf{x}_0) = \partial_t \nabla_{\mathbf{x}_0} \mathbf{d}(\mathbf{x}_0) = \partial_t \mathbf{F}^T(\mathbf{x}_0)$,

$$\int_{\Omega_t^s} 2c_1(\mathrm{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) : \nabla \mathbf{u} = \int_{\Omega_0^s} \frac{d}{dt} \Psi(\mathbf{F}) = \frac{d}{dt} \int_{\Omega_0^s} \Psi(\mathbf{I} + \nabla_{\mathbf{x}_0}^T \mathbf{d}).$$

The other terms are standard, in particular, with enough regularity on $\Omega_t^r$,

$$\int_{\Omega_t^r} (\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{u} = \frac{d}{dt} \int_{\Omega_t^r} \frac{|\mathbf{u}|^2}{2}, \quad r = f, s$$

**Remark 6** When $\Psi$ is convex, some regularity can be gained from this equality which can lead to existence of solution under severe regularity assumption; one of which (loosely speaking) is that $T$ cannot exceed the first time two pieces of the boundary which were not in contact initially come into contact. The first such result is by Coutand-Shkoller [13], later improved by Boulakia [7] and extended in a special case in [34].

## 5.2 Stability of the Scheme Discretized in Time

**Theorem 2** When $f = 0$ and $\rho$ is constant in each domain $\Omega_n^r$, $r = s, f$,

$$\int_{\Omega_n} \frac{\rho_n}{2} |\mathbf{u}^n|^2 + \delta t \sum_{k=1}^{n} \int_{\Omega_k^f} \frac{\mu}{2} |D\mathbf{u}^k|^2 + \int_{\Omega_0^s} \Psi^n \le \int_{\Omega_0} \frac{\rho_0}{2} |\mathbf{u}^0|^2 + \int_{\Omega_0^s} \Psi^0 \quad (14)$$

We begin the proof with the following lemmas,

**Lemma 1** The Lagrangian map $\mathbf{X}^n : \Omega_0 \mapsto \Omega_n$, implicit from (12), satisfies $\mathbf{X}^{n+1} = (\mathbb{Y}^{n+1})^{-1} \circ \mathbf{X}^n$, $n \ge 1$ and $\mathbf{F}^n := \nabla_{\mathbf{x}_0}^T \mathbf{X}^n = (\mathbf{I} - \nabla \mathbf{d}^n)^{-T}$.

Proof

Notice that $\mathbb{Y}^1(\mathbb{Y}^2(\dots\mathbb{Y}^{n-1}(\mathbb{Y}^n(\Omega_n))\dots)) = \Omega_0$ Hence

$$\mathbf{X}^{n+1} = [\mathbb{Y}^1(\mathbb{Y}^2(\dots\mathbb{Y}^n(\mathbb{Y}^{n+1})))]^{-1} = (\mathbb{Y}^{n+1})^{-1} \circ \mathbf{X}^n.$$

By definition of $\mathbf{d}^{n+1}$ in (13),

$$\mathbf{d}^{n+1}(\mathbf{X}^{n+1}(\mathbf{x}_0)) = \mathbf{d}^n(\mathbb{Y}^{n+1}(\mathbf{X}^{n+1}(\mathbf{x}_0))) + \delta t \mathbf{u}^{n+1}(\mathbf{X}^{n+1}(\mathbf{x}_0))$$
$$ = \mathbf{d}^n(\mathbf{X}^n(\mathbf{x}_0)) + \delta t \mathbf{u}^{n+1}(\mathbf{X}^{n+1}(\mathbf{x}_0)), \quad (15) $$

so $\mathbf{X}^{n+1}(\mathbf{x}_0) = \mathbf{d}^{n+1}(\mathbf{X}^{n+1}(\mathbf{x}_0)) + \mathbf{x}_0$ and therefore

$$ \begin{aligned} \mathbf{F}^{n+1} &= \nabla_{\mathbf{x}_0}^T (\mathbf{d}^{n+1}((\mathbf{X}^{n+1}(\mathbf{x}_0))) + \mathbf{x}_0), \\ &= \nabla \mathbf{d}^{n+1} {}^T \mathbf{F}^{n+1} + \mathbf{I} \Rightarrow \mathbf{F}^{n+1} = (\mathbf{I} - \nabla \mathbf{d}^{n+1})^{-T} \quad (16) \end{aligned} $$

**Lemma 2** For some $\alpha(\mathbf{x}, t)$,

$$ D\mathbf{d}^{n+1} - \nabla\mathbf{d}^{n+1}\nabla^T\mathbf{d}^{n+1} = J_{n+1}^{-2}[(\mathbf{I} - \nabla\mathbf{d}^{n+1})^{-T}(\mathbf{I} - \nabla\mathbf{d}^{n+1})^{-1} + \alpha\mathbf{I} $$

*Proof* With $\beta = \det_{\mathbf{I}-\nabla\mathbf{d}^{n+1}}^2$, and $\gamma$ the trace of $(\mathbf{I} - \nabla\mathbf{d}^{n+1})(\mathbf{I} - \nabla\mathbf{d}^{n+1})^T$, by the Cayley-Hamilton theorem, we have

$$ \begin{aligned} D\mathbf{d}^{n+1} - \nabla\mathbf{d}^{n+1}\nabla^T\mathbf{d}^{n+1} - \mathbf{I} &= -(\mathbf{I} - \nabla\mathbf{d}^{n+1})(\mathbf{I} - \nabla\mathbf{d}^{n+1})^T \\ &= \beta ((\mathbf{I} - \nabla\mathbf{d}^{n+1})(\mathbf{I} - \nabla\mathbf{d}^{n+1})^T)^{-1} - \gamma\mathbf{I} \end{aligned} $$

**Lemma 3**

$$ \int_{\Omega_{n+1}^s} c_1[J_{n+1}^{-1}(\mathbf{I} - \nabla\mathbf{d}^{n+1})^{-T}(\mathbf{I} - \nabla\mathbf{d}^{n+1})^{-1}] : D\hat{\mathbf{u}} = \int_{\Omega_0^s} \partial_{\mathbf{F}}\Psi^{n+1} : \nabla_{\mathbf{x}_0}\hat{\mathbf{u}} $$

*Proof*

From Lemma 1,

$$ \begin{aligned} & \int_{\Omega_{n+1}^s} c_1[J_{n+1}^{-1}[(\mathbf{I} - \nabla\mathbf{d}^{n+1})^{-T}(\mathbf{I} - \nabla\mathbf{d}^{n+1})^{-1}] : D\hat{\mathbf{u}} \\ &= \int_{\Omega_{n+1}^s} c_1 J_{n+1}^{-1}[\mathbf{F}^{n+1}\mathbf{F}^{n+1}{}^T] : D\hat{\mathbf{u}} = \int_{\Omega_0^s} c_1 \mathbf{F}^{n+1} : D_{\mathbf{x}_0}\hat{\mathbf{u}} \\ &= \frac{1}{2} \int_{\Omega_0^s} \partial_{\mathbf{F}}\Psi^{n+1} : D_{\mathbf{x}_0}\hat{\mathbf{u}} = \int_{\Omega_0^s} \partial_{\mathbf{F}}\Psi^{n+1} : \nabla_{\mathbf{x}_0}\hat{\mathbf{u}} \end{aligned} \quad (17) $$

## 5.3 Proof of Theorem 2

Let $r = s$ or $f$. Let us choose $\hat{\mathbf{u}} = \mathbf{u}^{n+1}$ in (11,12,13). By Schwartz inequality

$$ \int_{\Omega_{n+1}^r} (\rho_n^r \mathbf{u}^n) \circ \mathbb{Y}^{n+1} \cdot \mathbf{u}^{n+1} \le \left( \int_{\Omega_{n+1}^r} (\sqrt{\rho_n^r} \mathbf{u}^n)^2 \circ \mathbb{Y}^{n+1} \right)^{\frac{1}{2}} \left( \int_{\Omega_{n+1}^r} (\sqrt{\rho_{n+1}^r} \mathbf{u}^{n+1})^2 \right)^{\frac{1}{2}} $$

because $\rho_n^r \circ \mathbb{Y}^{n+1}(\mathbf{x}) = \rho_{n+1}^r(\mathbf{x})$, $\mathbf{x} \in \Omega_{n+1}^r$. So by a change of variable

$$ \int_{\Omega_{n+1}^r} \rho_{n+1}^r (\mathbf{u}^n \circ \mathbb{Y}^{n+1})^2 = \int_{\Omega_{n+1}^r} (\sqrt{\rho_n^r} \mathbf{u}^n)^2 \circ \mathbb{Y}^{n+1} = \int_{\Omega_n^r} \rho_n^r \mathbf{u}^{n2}. $$
Consequently, using $ab \le \frac{1}{2}a^2 + \frac{1}{2}b^2$,

$$ \int_{\Omega_{n+1}^r} \rho_{n+1}^r \mathbf{u}^n \circ \mathbb{Y}^{n+1} \cdot \mathbf{u}^{n+1} \le \frac{1}{2} \int_{\Omega_n^r} \rho_n^r \mathbf{u}^{n+2} + \frac{1}{2} \int_{\Omega_{n+1}^r} \rho_{n+1}^r \mathbf{u}^{n+1} $$

Finally,

$$ \int_{\Omega_{n+1}} \frac{\rho_{n+1}}{2} |\mathbf{u}^{n+1}|^2 + \delta t \int_{\Omega_{n+1}^f} \frac{\mu}{2} |\mathbf{D}\mathbf{u}^{n+1}|^2 + \int_{\Omega_0} \Psi^{n+1} \le \int_{\Omega_n} \frac{\rho_n}{2} |\mathbf{u}^n|^2 + \int_{\Omega_0^g} \Psi^g $$

**Remark 7** Inequality (18) shows that energy decays at each time step due to fluid viscous effect. Notice however that to reach this result we had to reintroduce in the variation formulation the second order terms in $\delta t$ which we intend to neglect in the numerical implementation; in particular $J_{n+1}$ will be replaced by 1. No significant influence of these simplifications on the numerical tests below has been observed.

# 6 Spatial Discretization with Finite Elements

Let $\mathcal{T}_h^0$ be a triangulation of the initial domain. Spatial discretization can be done with 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 be 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 triangulation, continuous in $\Omega_{n+1}^r$, $r = s, f$; $V_{0h}$ is the subspace of $V_h$ of functions which are 0 on $\Gamma$. A small penalization with parameter $\epsilon$ must be added to impose uniqueness of the pressure.

Discretization in space by the Finite Element Method leads to find $\mathbf{u}_h^{n+1}, p_h^{n+1} \in V_{0h} \times Q_h$ such that for all $\hat{\mathbf{u}}_h, \hat{p}_h \in V_{0h} \times Q_h$,

$$ \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} \right. \\ & \qquad \left. + \mathbf{1}_{\Omega_{n+1}^f} \frac{\mu}{2} \mathbf{D}\mathbf{u}_h^{n+1} : \mathbf{D}\hat{\mathbf{u}}_h \right. \\ & \qquad \left. + c_1 \mathbf{1}_{\Omega_{n+1}^s} [\mathbf{D}(\tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1}) - \nabla(\tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1}) \nabla^T(\tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1})] : \mathbf{D}\hat{\mathbf{u}}_h \right] \\ & = \int_{\Omega_{n+1}} \mathbf{f} \cdot \hat{\mathbf{u}}_h, \quad \Omega_{n+1} = (\mathbb{Y}^{n+1})^{-1}(\Omega_n) = \{\mathbf{x} : \mathbb{Y}^{n+1}(\mathbf{x}) \in \Omega_n\} \quad (19) \end{aligned} $$

with $\mathbf{d}_h^n$ updated by $\mathbf{d}_h^{n+1} = \tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1}$ where $\tilde{\mathbf{d}}_h^n = \mathbf{d}_h^n \circ \mathbb{Y}^{n+1}$ and where

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

We claim that the proof of Theorem 2 used for the spatially continuous case works also for the discrete case if the discrete Lagrangian maps verify

$$ \mathbf{X}^n = \mathbf{X}^{n+1} \circ \mathbb{Y}^{n+1}. \quad (20) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Sketch to understand if $X^n = \mathbb{Y}^{n+1} \circ X^{n+1}$; on the left the case of $P^2$-isoparametric element for the velocities and on the right the case of the $P_3^1$-$P^1$ element where each triangle is divided into four subtriangles on which the velocities are $P^1$ and continuous. A triangle $T_0^k$ in the reference domain (chosen here to be its initial position at time zero) becomes the isoparametric triangles $T_n^k$ and $T_{n+1}^k$: $T_n^k = X^n(T_0^k)$ and $T_{n+1}^k = X^{n+1}(T_0^k)$ in the case of quadratic velocities. But as the vertices and mid-edges of $T_n^k$ are obtained by moving those of $T_{n+1}^k$ by $-\delta\mathbf{u}^{n+1}$ we cannot have $\mathbb{Y}^{n+1}(T_{n+1}^k) = T_n^k$ because the composition of two quadratic maps is of degree 4 in general. However in the case of affine velocities, triangles are transformed into triangles and the mapping composition property holds.</figcaption></figure>

First let us consider the $P^2 - P^1$ element for $V_{0h} \times Q_h$. Then the mapping $\mathbf{x} \mapsto \mathbb{Y}^{n+1}(\mathbf{x})$ is quadratic and may preserve an isoparametric mesh of degree 2 if the vertices and the mid-edge nodes of the mesh at time $t^{n+1}$ are mapped by $\mathbb{Y}^{n+1}$ to the corresponding vertices and mid-edge nodes of the mesh at time $t^n$. Unfortunately such a construction does not satisfy (20) for the obvious reason that $f \circ g$ is of degree 4 when $f$ and $g$ are quadratic.

Now consider the $P_3^1 - P^1$ element: the fluid pressure and the solid pressure are continuous and piecewise linear on the triangulation of $\Omega_t^f$ and $\Omega_t^s$, separately; then each “pressure-triangle” is divided into 3 smaller triangles by a fourth vertex anywhere inside the triangle; on this subdivided triangulation the velocity is chosen continuous and piecewise linear.

Then (20) holds and the proof for the continuous case can be adapted, at least in the solid region. The inner vertex used to construct the fluid mesh will be moved by $\mathbb{Y}^{n+1}$ but $X^{n+1} \circ \mathbb{Y}^{n+1}$ remains linear and for each triangle $T_n^k = \mathbb{Y}^{n+1}(T_{n+1}^k)$ as long as the inner vertex is not transported outside the pressure triangle.

In the fluid region, an additional difficulty occurs when the mesh is not moved by $\mathbb{Y}^{n+1}$: a projection step from the old mesh to the new mesh must be added. Here it is done by polynomial interpolation according to the degree
of the element. A careful analysis of the algorithm shows that we never have to interpolate functions outside their domain of definition; for instance $\mathbf{u}^n$ is defined on $\Omega_n$ and needs to be interpolated when $\mathbf{u}^n(\mathbf{x} - \delta t \mathbf{u}^{n+1})$ is computed with $\mathbf{x} \in \Omega_{n+1}$. As $P^1$-interpolations from one mesh to another reduces the $L^2$ norms we conjecture that the energy decreases from one iteration to the next.

**Remark 8** A stabilized $P^1 - P^1$ element for pressures and velocities could be studied in place of $P_3^1 - P^1$, as in [14, 35] and others. Then a stabilization term $\omega \nabla \cdot \mathbf{u}^{n+1} \nabla \cdot \mathbf{v}$ is added to the first integral in the variational formulation (19) where $\omega$ is a mesh dependent parameter. The proof of Theorem 2 carries over except that $\int_{\Omega_{n+1}^f} \frac{\mu}{2} |D\mathbf{u}^{n+1}|^2$ becomes $\int_{\Omega_{n+1}^f} \frac{\mu}{2} |D\mathbf{u}^{n+1}|^2 + \int_{\Omega_{n+1}} \omega(\nabla \cdot \mathbf{u}^{n+1})^2$.

**Conclusion** In the case of solid only, Theorem 2 holds for the fully discrete problem (19) when linear elements are used. In the case of Fluid-Structure, Theorem 2 is likely to be true too but one would have to analyze mathematically the projection step from the mesh at time $n+1$ to the mesh at time $n$. This conjecture is based on the fact that interpolation of $P^1$ to $P^1$ on different meshes reduces the norm.

# 7 Contacts

To each simply connected disjointed part $S_t^i, i = 1..n_s$ of $\Sigma_t$ or $\Gamma_t$ is associated a signed distance function $\mathbf{x} \mapsto d_{S_t^i}(\mathbf{x})$ measuring the Euclidian distance of $\mathbf{x}$ to $S_t^i$ with the sign indicating whether $\mathbf{x}$ is in the structure ($d_{S_t^i}(\mathbf{x}) < 0$) or in the fluid ($d_{S_t^i}(\mathbf{x}) > 0$). When $d_{S_t^i}(\mathbf{x}) = 0$ for some $\mathbf{x} \in \partial\Omega \setminus S_t^i$ there is a contact. Note that contacts between points on the same part $S_t^i$ is not covered by this framework.

At contact points the equality in (5) becomes an inequality. Then the system is modeled by a variational inequality with Lagrange multiplier $\lambda$:

$$ \int_{\Omega_t} \left[ \rho \mathbb{D}_t \mathbf{u} \cdot \dot{\mathbf{u}} - p \nabla \cdot \hat{\mathbf{u}} - \hat{p} \nabla \cdot \mathbf{u} \right] + \int_{\Omega_t^f} \frac{\mu}{2} \mathrm{D}\mathbf{u} : \mathrm{D}\hat{\mathbf{u}} + \int_{\Omega_t^s} c_1 (\mathrm{D}\mathbf{d} - \nabla \mathbf{d} \nabla^T \mathbf{d}) : \mathrm{D}\hat{\mathbf{u}} + \int_{\partial\Omega_t} \lambda \mathbf{n} \cdot \hat{\mathbf{u}} = \int_{\Omega_t} \mathbf{f} \cdot \hat{\mathbf{u}}, \\ \lambda(\mathbf{x}, t) \le 0, \quad d_{S_t^i}(\mathbf{x})\lambda(\mathbf{x}, t) = 0 \quad \forall \mathbf{x} \in S_t^i, \forall t. \qquad (21) $$

To solve the problem we apply the discretization described above and the Semi-Smooth Newton method proposed by Ito and Kunisch [23] which replaces inequality constraints by equality constraints at each time step; here it is:

$$ \lambda^{n+1}(\mathbf{x}) - \min\{0, \lambda^{n+1}(\mathbf{x}) + c_0 d_{S_{n+1}^i}(\mathbf{x} + \delta t \mathbf{u}_n^{n+1}(\mathbf{x})\mathbf{n}(\mathbf{x}))\} = 0 \quad \forall \mathbf{x} \in S_{n+1}^i $$

and similarly on $\Gamma$; $c_0$ is any positive constant. These equality constraints are only left and right differentiable, but it is enough for a Newton-type algorithm
to converge. In the end one needs to solve iteratively in $k$ for each $n$,

$$ \int_{\Omega_{n+1,k}} \left[ \rho_{n+1} \frac{\mathbf{u}_h^{n+1,k+1} - \tilde{\mathbf{u}}_h^n}{\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} + \mathbf{1}_{\Omega_{n+1}^f} \frac{\mu^f}{2} \mathrm{D}\mathbf{u}_h^{n+1,k+1} : \mathrm{D}\hat{\mathbf{u}}_h \right. \\ + c_1 \mathbf{1}_{\Omega_{n+1}^s} [\mathrm{D}(\tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1,k+1}) - \nabla(\tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1,k}) \nabla^T (\tilde{\mathbf{d}}_h^n + \delta t \mathbf{u}_h^{n+1,k+1})] : \mathrm{D}\hat{\mathbf{u}}_h] \\ + \sum_{i=1}^{n_s} \int_{S_{n+1}^i \cap \{\mathbf{x}: \lambda + c_0 d_{S_{n+1}^i}(\mathbf{x} + \delta t \mathbf{u}_h^{n+1,k}(\mathbf{x})) < 0\}} c_3(\mathbf{x} + \delta t(\mathbf{n} \cdot \mathbf{u}_h^{n+1,k+1})\mathbf{n})\hat{\mathbf{u}}_h = \int_{\Omega_{n+1}} \mathbf{f} \cdot \hat{\mathbf{u}}_h, \\ \left. \Omega_{n+1,k+1} = \{\mathbf{x}: \mathbb{Y}^{n+1}(\mathbf{x}) := \mathbf{x} - \delta t \mathbf{u}_h^{n+1,k+1} \in \Omega_n\} \right] \quad (22) $$

where $c_3$ is a very large constant. The purpose of $c_3$ is to impose $\mathbf{x}+\delta t \mathbf{u}_h^{n+1,k+1}\mathbf{n}(\mathbf{x}) = 0$ on the boundary where the constraint is almost active (it is not a penalization of the constraint). Indeed it is standard practice to impose an homogeneous Dirichlet condition at node $i$ by changing the $i^{th}$ diagonal of the discrete system to a very large value so that off-diagonal terms do not contribute; here the same trick is used and a low order quadrature for the integral (mass lumping) will make sure that $c_3$ appears only on the diagonal of the linear system.

This shows that Eulerian formulations are well adapted to an easy treatment of contacts as observed by Richter[35]. The argument above shows also that semi-smooth Newton fits very well the Eulerian framework: it is simply an additional surface integral. Moreover the iterations of index $k$ can be combined with the fixed point iterations required by the algorithm. Notice also that the treatment of contact is symmetric; nowhere is it assumed that it is boundary $S_i$ which hits boundary $S_j$.

# 8 Numerical Tests

In all tests except one, $P^2-P^1$ elements are used. The vertices in the structures are moved by their own velocities. A procedure to extract the modified solid boundary has been implemented in the public domain software freefem++[18] and from this knowledge the mesh in the fluid region is rebuilt every time step by a Delaunay-Voronoi algorithm. The number of fixed point iterations is fixed at 2 in most tests; less results in a risk to lose mass, unless the time step is quite small, and more does not make a difference.

## 8.1 Validation with a Rotating Disk

The purpose of this test is to compare the numerical solution of (19) with a semi-analytical solution which can be computed to any desired accuracy.

A cylinder contains a fixed rigid cylindrical rod in its center, a cylindrical layer of hyperelastic material around the rod; the rest is filled with a fluid (see figure 3). First the system is at rest and then a constant rotation is given to the outer cylinder. This causes the fluid to rotate with an angular velocity which depends on the distance $r$ to the main axis; in turn because the friction of the fluid at the interface, the hyperelastic material will be dragged into a angular
<figure><img src="image_4.png" /><figcaption>Figure 3: A fluid-structure system inside a rotating cylinder (giving a constant angular velocity to the fluid outer boundary) with a fixed rod in its center. Left: sketch of the system. Right: a 2d calculation showing the velocity vectors at time 0.85 for the coarser mesh.</figcaption></figure>

velocity $u_θ$ which is also only a function of $r$ and time. Due to elasticity $u_θ$ will oscillate with time until numerical dissipation and fluid viscosity damps it.

In a two dimensional cut perpendicular to the main axis, the velocities and displacements are two dimensional as well. Hence the geometry is a ring of inner and outer radii, $R_0$ and $R_1$, with hyperelastic material between $R_0$ and $R$ and fluid between $R$ and $R_1$. Because of the incompressibility of the fluid and axial symmetry, $R$ is constant.

In this test $R_0 = 3$, $R = 4$, $R_1 = 5$. The solid is neo-Hookean with $c_1 = 2$ and $ρ^s = 2$. The Newtonian fluid has $μ = 2$, $ρ^f = 1$. The radial velocity of the outer cylinder has magnitude 3.

As everything is axisymmetric the computation can be done in polar coordinates $r, θ$, and the fluid-solid system reduces to

$$ρ∂_t u_θ - \frac{1}{r}∂_r [ξ^f r∂_r u_θ + ξ^s r∂_r d] = 0, \quad ∂_t d = u_θ, \ r ∈ (R_0, R_1), \ v|_{R_0} = 0, \ u_{θ|_{R_1}} = 3, (23)$$
with $\rho = \rho^s \mathbf{1}_{r \le R} + \rho^f \mathbf{1}_{r > R}$, $\xi^s = 2c_1 \mathbf{1}_{r \le R}$, $\xi^f = \mu \mathbf{1}_{r > R}$, and with $d(r, 0) = 0$.

Comparison between this one dimensional system and the Eulerian 2D system (19) is given on Figure 3, on the right, at $T = 0.85$ (which is the first maximum rotation of the solid) after 70 time steps with a coarse mesh with 520 vertices.

Then the same is computed on a finer mesh with 140 time steps and 756 vertices and finally with a mesh with 1773 vertices and 280 time steps. Results are displayed on Figure 4.

<figure><img src="image_5.png" /><figcaption>Figure 4: Rotating cylinder. Left: Evolution of the $L^2$ error on the radial velocity, between the semi-analytical solution and the 2D computation, versus time for the 3 meshes. Right: velocities normal to the ray at $\theta = \pi/4$ versus $r$, computed on the 3 meshes. The solution of the one dimensional equation is shown with green crosses; the curve corresponding to the finer mesh (blue crosses) is hard to see because it is hidden by the green crosses when $r - 3 < 0.3$ and then hidden by the 2 other curves when $r - 3 > 0.5$ (middle mesh in dotted blue and coarse mesh in red). Notice the discontinuity of the gradient of the velocity at the fluid-solid interface.</figcaption></figure>
This test has two qualities: a/ the exact solution is easy to compute to any
precision; b/ the geometry does not change and hence the spatial discretization
of the Eulerian formulation can be assessed independently to the equation which
relates the displacement to the velocity because it can be simplified to $\mathbf{d}^{n+1} =
\mathbf{d}^n + \delta t \mathbf{u}^{n+1}$ instead of (8).

However, even with this simplification, Figure 4 on the right shows that the $L^2$ error is closer to $O(h)$ where $h$ is the mesh size. The right side of the same figure shows that the solution is not smooth: there are slope discontinuities at the fluid-structure interface, as expected. So optimal order convergence with $P^2$ elements for the velocities cannot be expected. There may be another reason, similar to the Babuska paradox[2] (which says that in some problems involving the trace of the stress on a curved boundary, approximated by a polygonal line, convergence does not occur when the polygon converges to the curved line: isoparametric elements have to be used).

**Conclusion:** *The method converges but not at as fast as expected from the degree of the polynomial approximation.*

## 8.2 Flag Attached to a Cylinder

*This is a reference test for compressible material, but not yet for incompressible structures.*

The test is called FLUSTUK-FSI-3 in [16]; it is a well established challenge, considered to be hard, for compressible structures; FLUSTUK-FSI-3* is the incompressible variant. A beam, shaped like a flag, made of an hyperelastic neo-Hookean material is clamped at the back of a hard fixed cylinder in a fluid contained in a rectangular pipe. The cylinder is slightly off the symmetry line of the pipe so as to trigger vortex shedding at a relatively low Reynolds number.

**Geometry** The flag is a rectangle of size $l \times h$ at $t = 0$. It is attached behind a cylinder of radius $r$ so that the top and bottom left corners are symmetrically on the cylinder; the left vertical boundary of the rectangle is inside the cylinder. The fluid computational domain is a rectangle of size $L \times H$; the flow enters from the left and is free to leave on the right. The center of the cylinder is at $(c, c)$ with $c = 0.2$; other parameters are $l = 0.35$, $h = 0.02$, $L = 2.5$, $H = 0.41$.

**Fluid** The density is $\rho^f = 10^3 \text{kg/m}^3$ and the reduced viscosity $\nu^f = \frac{\mu}{\rho^f} = 10^{-3} \text{m}^2/\text{s}$; the inflow velocity is horizontal and parabolic in $y := x_2$ with flux $\bar{U}H$, with $\bar{U} = 2$; namely $\mathbf{u}(0, y) = \bar{U} \left( \frac{6}{H^2} y(H-y), 0 \right)^T$. Top and bottom boundaries are no-slip walls.

**Solid** $\rho^s = \rho^f$, $c_1 = 10^6 \text{kgm}^{-1}\text{s}^{-1}$ and no external force.
<figure><img src="image_6.png" /><figcaption>Figure 5: An hyperelastic beam (flag) attached to a cylinder in a Navier-Stokes flow. A snapshot of the velocity norms in the full computational domain at $t = 3.786$.</figcaption></figure>

Initially all velocities and displacements are zero. Around $t \sim 2$, the flow starts to oscillate and creates a Karman vortex street.

Figure 5 shows the position of the flag at $t = 3.786$ and a color map of the velocity norms, while Figure 6 shows the x and y coordinates of a vertex on the right vertical tip of the flag as a function of time. The results displayed are obtained with a mesh with 2200 vertices and a time step equal to 0.005, as in Dunne[14]. The results concur roughly with [16] in the sense that the frequency is around $5.4s^{-1}$ compared to $6.7s^{-1}$ and the amplitude is around 0.015 compared to 0.013 in [16] however we do not consider this agreement reliable.

In truth, mesh and time step refinements bring the frequency and the amplitude closer to $5.4s^{-1}$ and 0.025 respectively. Moreover the amplitude is very sensitive to $c_1$: a 10% increase induces a two-fold change.

Let it be noted, for comparison, that the same approach gives the correct frequency (within 10%), when the structure is a compressible Mooney-Rivlin hyperelastic material (test case FSI-3 (see [31])) and also in absence of flag (Navier-Stokes flow around a cylinder at $Re = 100$). This disagreement, which may be due to unconverged results, will not be solved until more participants do the tests.

## 8.3 Clamped Beam Falling Freely in a Fluid

The purpose of this test is to check how the algorithm preserves the invariants: energy and mass in the case of large displacements and compare with [16].

In [16] the free fall under its own weight of the flag attached to the cylinder, initially horizontal, is studied; the geometry is identical to the previous test.
<figure><img src="image_7.png" /><figcaption>Figure 6: On the left (resp. right) the x (resp. y) position of a vertex on the right tip of the flag versus time. The frequency is around 5.4s⁻¹ and the amplitude around 0.03. The mesh has 2200 vertices and the time step is 0.005.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Final position of the flag at t=0.98 after it has fallen under its own weight.</figcaption></figure>
The density of the solid is now 20, gravity is 1 and $c_1 = 0.2510^6$; all other parameters are identical to the previous case. After 50 time steps the lower right tip of the flag is at height $y = 0.01$ at $t = 0.97$. In the present simulation $y = 0.0174$ at $t = 0.98$. Figure 7 shows the velocity norms at this time.

We have used this test case to study conservation of mass, i.e. the area of the flag versus time, and conservation of energy, i.e. the value versus time of the relative difference between the energy at time $n + 1$ and the initial energy minus the work of the force (gravity) at time $n + 1$:

$$1 - \left( \int_{\Omega_{n+1}} \frac{\rho_{n+1}}{2} |\mathbf{u}^{n+1}|^2 + \sum_{j=1}^{n} \delta t \int_{\Omega_{j+1}} \frac{\mu}{2} |\mathbf{D}\mathbf{u}^{j+1}|^2 + \int_{\Omega_0} \Psi^{n+1} \right) / \left( \int_{\Omega_0} \Psi^0 - \sum_{j=1}^{n+1} \delta t \int_{\Omega_j} g \rho_j \mathbf{u}_j^j \right) \quad (24)$$

Three elements have been tested: $P^2 - P^1$, $P^1 \text{bubble} - P^1$ and $P^1 - P^1$ without stabilization (it seems that using $D$ instead of $\nabla$ for the viscous term has a stabilizing effect).

Figure 8 displays the energy index (24) versus time. It shows that the $P^1 - P^1$ element is the best, as the theory predicts. It also shows that one iteration in the fixed point algorithm is not enough, two is best (similar plots with 3 iterations have been given in [30][Figure 5] and show no improvement over 2).

<figure><img src="image_9.png" /><figcaption>Figure 8: Free fall of a beam, clamped on the left, in a fluid initially at rest. The energy index (24 is plotted, versus time, for several finite elements and one or two fixed point iterations. Notice that $P^1 - P^1$ with 2 fixed point iterations is best.</figcaption></figure>

Finally Figure 9 displays the area of the flag versus time. Here we see that no element is superior, but no more that 1% is lost in 50 time steps.
<figure><img src="image_10.png" /><figcaption>Figure 9: *Free fall of a beam, clamped on the left, in a fluid initially at rest. Conservation of mass when using various finite elements and one or two fixed point iterations.*</figcaption></figure>

We have also studied the $P_3^1 - P^1$ element with moving mid vertices; small differences are observed but not worth the switch from $P^1 bubble - P^1$ to $P_3^1 - P^1$. Figure 10 shows the triangulation at one instant of time during the computation. It is seen that most inner vertices stay fairly in the center of their pressure triangles. Small differences are observed for energy balance and surface area of the beam but not worth implementing the $P_3^1 - P^1$ element.

<figure><img src="image_11.png" /><figcaption>Figure 10: *The $P_3^1 - P^1$ element is shown in action: zoom on the right tip of the beam showing that most but not all inner vertices are in the center of the pressure triangles. The colors correspond to the norm of the velocity vectors. The edges in the structure are red.*</figcaption></figure>

**Conclusion:** The $P^1 - P^1$ element conserves energy better but all element tested conserve energy and mass within 1% after 50 time steps even with a
rather coarse triangulation. One fixed point iteration is not enough and we
suggest to use 2. We also suggest to use the stabilized P¹ – P¹ element instead
of P³¹ – P¹ element.

8.4 Bouncing Ball

The purpose of this test is to analyze the performance of the algorithm in the
case of contact. The test was proposed by Thomas Richter[35] and Stefan Frei
in [19].

A circular structure of radius $R = 0.4$, centered at $\mathbf{x}_c = 1$, $y_c = 1$, in a
rectangular cavity $[0, 2] \times [0, 1.5]$ filled with a fluid of density $\rho^f = 1000$ and
viscosity $\mu = 1000$. The lateral walls of the cavity are numerical boundaries
where no condition is imposed in the variational formulation, i.e. $\sigma^f \mathbf{n} = 0$.

The Mooney-Rivlin coefficient is $c_1 = 10^5$ and $\rho^s = 1000$. The structure is
subject to a gravity $g = 1$ but the fluid is not. Hence the ball falls till it hits
the bottom of the cavity and then rebounds.

A snapshot shows the velocity norms in the middle of the rebound on Figure 11.

In Figure 12 the position of the lowest point on the disk is shown, versus
time, together with the vertical diameter of the disk, i.e. the distance between
the lowest and highest points on the disk versus time.

The results are identical to those obtained by Frei [19][Figure 12.4 page 138] (see also [35] for another case and [20] for a description of their Eulerian formulation)

**Conclusion:** The present Eulerian formulation is well adapted to contacts and agrees with a simulation by another team with a similar method but using penalty instead of Semi-Smooth Newton iterations.

Acknowledgement

Implementation was done with the PDE solver FreeFem++[18]; the program
scripts can be sent upon request by email to the second author.

**Summary and Conclusion**

An Eulerian fluid-structure formulation has been presented with steps to derive
an implicit in time energy stable scheme. The scheme is monolithic in the sense
that every iteration in a fixed point loop involves a linear solver on the velocities
and pressures only. Finite Element discretizations are proposed and studied
and one of them has been shown to be energy stable; it involves a standard P¹
element for the pressures (one for the fluid part and another for the structural
part) and a fancy P¹ element for the velocities on a mesh obtained by dividing
<figure><img src="image_12.png" /><figcaption>Figure 11: Free fall of an hyperelastic disk in a fluid with rebounds on the bottom. Two snapshots, of the velocity norms show, at $t = 1.65$, the beginning of the contact and at $t = 1.75$, its end.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Free fall of an hyperelastic disk in a fluid with rebounds on the bottom. On the left the height of the lowest point on the disk versus time computed with 3 triangulations each twice finer than the previous one and similarly for the time steps; the finest results are shown in blue. On the right the maximum vertical diameter of the disk versus time computed with the finest mesh.</figcaption></figure>
the pressure elements into 3 subtriangles and an inner vertex which is allowed to
move with the structure velocity field when the mesh is updated. Alternatively
the stabilised $P^1 - P^1$ element can also be used with the same energy stability
property.

While it is interesting to know that there is at least one finite element space
which is energy stable, in practice it is unnecessarily complex and the usual
$P^2 - P^1$ or $P^1 \text{bubble} - P^1$ elements or any other compatible element for the
inf-sup condition are acceptable.

The method has been implemented with the software freefem++[18] and
validation tests have been presented. Validation by comparison with other au-
thors is difficult because not many FSI codes exist for incompressible material
and large displacement.

The scheme is reasonably robust; occasionally it stalls while updating the
mesh in the structural region when an element is turned over by the structural
velocity. In principle such a situation could be fixed by remeshing; however
numerical tests show that when this happens a few time iterations later the
same trouble arises, indicating a deeper flaw, perhaps in the realizability of the
system.

We have also identified a precision issue due to the lack of regularity of the
velocity at the fluid-structure interface.

Nevertheless, the scheme is unusually stable and it passes most of the tests
considered difficult by the community, including contact. The code is also fast
and gives meaningful results even with coarse triangulations.

A generalization to compressible structure is available [31] and a 3D imple-
mentation is being developed [11]. Finding a second order in time extension
seems much harder.

References

[1] S.S. ANTMAN. Nonlinear Problems of Elasticity (2nd ed.). Springer series in Applied Mathematical Sciences, Vol. 107, 2005.

[2] I. BABUSKA, J. PITKARANTA. The Plate Paradox for Hard and Soft Simple Support SIAM J. Math. Anal., 21 (3, 1990): 551-576.

[3] K.J. BATHE. Finite element procedures, Prentice-Hall, Englewood Cliffs, New-Jersey. 1996.

[4] Ted Belytschko, Wing Kam Liu, Brian Moran, Khalil i. Elkhodary Nonlinear Finite elements For Continua and structures, Wiley 2014.

[5] D. BOFFI, N. CAVALLINI, L. GASTALDI, The finite element immersed boundary method with distributed Lagrange multiplier, arXiv:1407.5184v2, 2015 (to appear)
[6] K. BOUKIR, Y. MADAY, B. METIVET. A high order characteristics method for the incompressible Navier-Stokes equations. Comp. Methods in Applied Mathematics and Engineering 116 (1994), 211-218.

[7] M. BOULAKIA, Existence of weak solutions for the motion of an elastic structure in an incompressible viscous fluid, C. R. Math. Acad. Sci. Paris 336 No 12 (2003), 985-990.

[8] M. BUKACA, S. CANIC, R. GLOWINSKI, J. TAMBACAC, A. QUAINIA. Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement. Journal of Computational Physics 235 (2013) 515541.

[9] G. CHECHKIN, D. LUKKASSEN, A. MEIDELLOn the Sapondzhyan-Babuska Paradox, Applicable Analysis 87(12):1443-1460 December 2008.

[10] P.G. CIARLET. Mathematical Elasticity. North Holland, 1988.

[11] C.Y. CHIANG, O. PIRONNEAU, M. THIRIET AND TONY SHEU Numerical Study of a 3D Eulerian Monolithic Formulation for Fluid-Structure-Interaction. Special volume on Mechanics of Fluid-Particles Systems and Fluid-Solid Interactions. To appear in *Fluids*.

[12] TH. COUPEZ, L. SILVA, E. HACHEM, Implicit Boundary and Adaptive Anisotropic Meshes, S. Peretto, L. Formaggia eds. *New challenges in Grid Generation and Adaptivity for Scientific Computing*. SEMA-SIMAI Springer series, vol 5, 2015.

[13] D. COUTAND AND S. SHKOLLER, Motion of an elastic solid inside an incompressible viscous fluid, Arch. Ration. Mech. Anal. 176 No 1 (2005), 25-102.

[14] TH. DUNNE, Adaptive Finite Element Approximation Of Fluid-Structure Interaction Based On An Eulerian Variational Formulation. Dissertation, University of Heidelberg, 2007 and ECCOMAS CFD 2006 P. Wesseling, E. Oñate and J. Périaux (Eds) Elsevier, TU Delft, The Netherlands, 2006.

[15] TH. DUNNE. An Eulerian approach to fluid-structure interaction and goal-oriented mesh adaptation. Int. J. Numer. Meth. Fluids 2006; 51:1017-1039.

[16] TH. DUNNE AND R. RANNACHER Adaptive Finite Element Approximation of Fluid-Structure Interaction Based on an Eulerian Variational Formula- tion. In *Fluid-Structure Interaction: Modelling, Simulation, Optimization*. Lecture Notes in Computational Science and Engineering, vol 53. p110-146, H-J. Bungartz, M. Schaefer eds. Springer 2006.

[17] M. A. FERNANDEZ, J. MULLAERT, M. VIDRASCU. Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures, Comp. Methods in Applied Mech. and Engg. 267, 566593, 2013.
[18] F. HECHT New development in FreeFem++, J. Numer. Math., 20 (2012), pp. 251-265. www.FreeFem.org.

[19] S. FREI. Eulerian finite element methods for interface problems and fluid-structure interactions. Ph.D. dissertation, University of Heidelberg, 2016.

[20] S. FREI, T. RICHTER, AND T. WICK. Long-term simulation of large deformation, mechano-chemical fluid-structure interactions in ALE and fully Eulerian coordinates. Journal of Computational Physics, Volume 321 Issue C, September 2016 Pages 874-891.

[21] P. HAURET. Méthodes numériques pour la dynamique des structures non-linéaires incompressibles à deux échelles. Doctoral thesis, Ecole Polytechnique, 2004.

[22] J. HRON AND S. TUREK. A monolithic fem solver for an ALE formulation of fluid-structure interaction with configuration for numerical benchmarking. European Conference on Computational Fluid Dynamics ECCOMAS CFD 2006 P. Wesseling, E. Onate and J. Periaux (Eds). TU Delft, The Netherlands (2006)

[23] K. ITO AND K. KUNISCH, Semi-Smooth Newton Method for Variational Inequalities of the First Kind, ESAIM:M2AN, Vol. 37, No 1, 2003, pp. 4162.

[24] J. MARSDEN AND T.J.R. HUGHES Mathematical Foundations of Elasticity. Dover publications 1993.

[25] P. LE TALLEC AND P. HAURET. Energy conservation in fluid-structure interactions. In P. Neittanmaki Y. Kuznetsov and O. Pironneau, editors, Numerical methods for scientific computing, variational problems and applications, CIMNE, Barcelona, 2003.

[26] P. LE TALLEC AND J. MOURO. Fluid structure interaction with large structural displacements. Comp. Meth. Appl. Mech. Eng., 190(24-25):3039-3068, 2001.

[27] J. LIU. A second-order changing-connectivity ALE scheme and its application to FSI with large convection of fluids and near-contact of structures. Submitted to Journal of Computational Physics, (2015).

[28] C.S. PESKIN. The immersed boundary method. Acta Numerica, 11:479-517, 2002.

[29] O. PIRONNEAU. *Finite Element Methods for Fluids*. Wiley, 1989.

[30] O. PIRONNEAU. Numerical Study of a Monolithic Fluid-Structure Formulation. in the proceedings of Erices *Variational Analysis and Aerospace Engineering: Mathematical Challenges for the Aerospace of the Future* Aldo Frediani, Bijan Mohammadi, Olivier Pironneau (eds). Springer 2016.
[31] O. PIRONNEAU. An Eulerian Monolithic Fluid-Structure Formulation. *Anniversary volume in honor of P.G. Ciarlet. Li-Ta-Tsien ed.* 2018. Also in *hal.archives-ouvertes.fr/hal-01348648*.

[32] L. FORMAGGIA, A. QUARTERONI, A. VENEZIANI. *Cardiovasuclar Mathematics*. Springer MS&A Series. Springer-Verlag, 2009.

[33] R. RANNACHER AND T. RICHTER, An Adaptive Finite Element Method for Fluid-Structure Interaction Problems Based on a Fully Eulerian Formulation. Springer Lecture Notes in Computational Science and Engineering, Vol 73, 2010.

[34] J-P. RAYMOND, M. VANNINATHAN A fluid-structure model coupling the Navier-Stokes equations and the Lamé system. J. Math. Pures Appl. 102 (2014) 546-596

[35] TH. RICHTER. A Fully Eulerian formulation for fluidstructure interaction problems. Journal of Computational Physics, 233:227-240, 2013.

[36] TH. RICHTER AND TH. WICK. Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates. Comput. Methods Appl. Mech. Engrg. 199 (2010) 2633-2642.

[37] TH. WICK. Fully Eulerian fluidstructure interaction for time-dependent problems. Comput. Methods Appl. Mech. Engrg. 255 (2013) 14-26.