Sergei Gluchko

Bruno Palpant

Sebastian Volz

Rémy Braive

Thomas Antoni

Thermal Excitation of Broadband and Long-range Surface Waves on SiO 2 Submicron Films

Keywords: numbers: 44.40.+a, 71.36.+c, 78.20.-e Keywords: thermal emission, nearfield thermal radiation, Zenneck modes, surface phonon polaritons

We detect thermally excited surfaces waves on a submicron SiO 2 layer, including Zenneck and guided modes in addition to Surface Phonon Polaritons. The measurements show the existence of these hybrid thermalelectromagnetic waves from near-(2.7 µm) to far-(11.2 µm) infrared. Their propagation distances reach values on the order of the millimeter, several orders of magnitude larger than on semi-infinite systems. These two features; spectral broadness and long range propagation, make these waves good candidates for near-field applications both in optics and thermics due to their dual nature.

Thermal radiation through surface wave diffraction is usually only considered as the result of Surface Phonon Polaritons (SPhPs). SPhPs are hybrid evanescent electromagnetic surface waves generated by the phononphoton coupling, at the interface of polar and dielectric materials (such as SiO 2 and air) [1][2][START_REF] Agranovich | Surface Polaritons: Electromagnetic Waves at Surfaces and Interfaces[END_REF][START_REF] Kim | [END_REF] . The influence of SPhPs on the thermal performance of nanostructured materials has been studied intensively over the last decade, providing an alternative channel of heat conduction when the objects are scaled down 5,6 . Due to this behaviour, they are essential for the improvement of the thermal stability in micro and nanoelectronics [7][8][9] , microscopy 10 , near-field thermophotovoltaics 11 and for thermal radiation [12][13][14] . In addition, SPhPs provide coherent thermal radiation in mid-infrared 13,14 . This feature is now widely used to control thermal radiation but in a frequency range that is limited to the midinfrared because it implies the coupling to transverse optical phonons 1,15 . But this narrow spectrum (typically 8.6 -9.3 µm at a SiO 2 -air interface) in addition to propagation lengths in the range of the wavelength decrease the field of use of SPhPs for many applications such as thermal transport at nanoscale, infrared nanophotonics and coherent thermal emission.

In this letter we demonstrate through experiment that coherent thermal emission, resulting from surface waves, can be extended spectrally. We also prove experimentaly that these surface waves have a long propagation range, when considering isolated submicron layers. Indeed, if the film is thinner than the penetration depth of a) Corresponding author:thomas.antoni@centralesupelec.fr the wave inside the material, the electromagnetic mode can be coupled on both its interfaces allowing for the long-range propagation of two other types of electromagnetic surface waves; Zenneck and subwavelength Transverse Magnetic (TM) guided modes 16 . The propagation length is increased as a consequence of the dramatic decrease in the overlap of the mode with the material, hence its absorption. For example, it is almost two orders of magnitude larger than the wavelength for a 1 µm thick suspended SiO 2 membrane 5 . To prove those predictions, we fabricated a submicron glass layer and characterized its thermal emission by means of Fourier Transform In-fraRed (FTIR) spectroscopy. Our experimental results are then compared with both Finite-Difference Time-Domain (FDTD) simulations and theoretical predictions.

The fabrication process of a suspended submicron thick glass membrane is, however, a challenging task due to the poor selectivity of Si/SiO 2 etching. Nevertheless, as far as electromagnetism is concerned, a sample consisting of a twice thinner layer of SiO 2 deposited over a metallic film -the role of which is to optically isolate the dielectric thin layer from the substrate -is a strictly equivalent structure. To prove this equivalence, Fig. 1(a) reports the calculated dispersion curves of surface waves on a 1 µm thick suspended SiO 2 membrane and a 0.5 µm thick SiO 2 thin film deposited on an aluminum layer obtained by means of FDTD simulations (MEEP code 17 ). This is done by considering real and imaginary parts of glass dielectric function of SiO 2 shown in Fig. 1(b) as extracted from experimental data [START_REF] Palik | Handbook of Optical Constants of Solids[END_REF] . The dispersion curves are superimposed, confirming the electromagnetic equivalence of the two systems. In addition, it can be seen that the curves lie beneath the light line over nearly the full spec- The sample is fabricated by sputtering deposition of an amorphous SiO 2 thin layer on an aluminum layer, itself grown on a polished Si wafer in order to optically mimic the suspended layer of SiO 2 . The thicknesses of SiO 2 and aluminum layers are chosen to be 0.75 µm and 0.25 µm, respectively. A diffraction grating of period Λ = 9.26 µm with a filling factor 0.5 is etched by 0.5 µm in the SiO 2 film using negative UV lithography and anisotropic reactive ion etching. In order to avoid any artificial broadening of the emission spectrum due to finite size effects, the grating has to be much larger than the propagation length of these surface modes. We then choose the lateral size of the grating to be 1 cm in both length and width (as shown on Fig. 2(a)), that is one order of magnitude larger than the maximum value of the estimated propagation length.

We operate a FTIR spectrometer with a spectral resolution of 1 cm -1 , which provides the required sensitivity in the working frequency range from 700 cm -1 to 4000 cm -1 (14.3 -2.5 µm). The sample is heated up to 673 K with a heating stage and the emitted signal is collected for various tilt angles of the sample. KRS-5 holographic wire grid polarizer is used to examine sample emission for TE and TM polarizations providing excellent transmission in the working frequency range. Each experimental spectrum S exp gr is obtained by subtraction of the background radiation S bg and normalization by the emission obtained from a flat region of the sample, without grating S f l , under the same conditions, as follows:

S exp gr (ν, T s , α) = S gr (ν, T s , α) -S bg (ν) S f l (ν, T s , α) -S bg (ν) , (1)
where T s is the sample surface temperature and α is the tilt angle. Note that this emission signal S exp gr (ν, T s , α) is not the same as emissivity as it also features the influence of the topography of the sample on its emission. The normalized spectrum value reduces to unity for any frequency where the grating has no impact on the sample emission. Fig. 2(b) shows this normalized emission spectrum obtained from our sample for T s = 673 K. The sharp peaks, marked with black arrows, are the diffraction orders of the grating which is a clear signature of coherent thermal emission. The peaks without the arrows do not originate from the diffraction by the grating and they are observed due to the different SiO 2 effective thickness in the grating and flat regions of the sample.

We then examine the polarization of the grating emission since the diffraction peaks can only exist for transverse magnetic polarization if they originate from surface waves 16,[START_REF] John Polo | Electromagnetic Surface Waves: A Modern Perspective[END_REF] . Fig. 3(a,c,e) show the grating emission peaks observed by collecting emission signal from the grating on SiO 2 thin film deposited on aluminum substrate heated up to T = 673 K in three different frequency regions. In all these regions the emission peaks disappear for TE polarized signal as expected.

These frequency regions indicate different surface modes according to the values of the relative permittivity of glass shown in Fig. 1(b). These are Zenneck surface modes ( r > 0, i > 0), Surface Phonon Polaritons ( r <air ), and subwavelength TM guided modes ( r > 0, i ≈ 0) according to the classification of surface electromagnetic modes in thin dielectric films 16,[START_REF] John Polo | Electromagnetic Surface Waves: A Modern Perspective[END_REF] . Detection of Zenneck surface modes in a single interface system for these frequencies is usually not possible due to their short propagation length 16 diffraction by the grating. Fig. 3(e) shows the existence of thermally excited subwavelength TM guided modes over a broad frequency range around ν = 2000 cm -1 where the absorption is almost equal to zero. Fig. 3(b,d,f) demonstrate the frequency shift of the emission peaks when the tilt angle is varied. This feature can be used reconstruct the dispersion relation of these surface modes by applying the grating equation:

ω c sin α = k + m 2π Λ , m ∈ Z, (2)
where k is the in-plane wavevector and m is the diffraction order. Fig. 4(a) reports the dispersion relation that has been obtained by following this procedure. Experimental data lie beneath the light line indicating the evanescent nature of the surface waves and is in reasonably good agreement with the numerical predictions. The difference between experimental results and FDTD simulation data can be understood by assuming that the effective thickness of the SiO 2 film with the diffraction grating is not exactly 0.5 µm due to fabrication discrepancies and that the permittivities of the SiO 2 sample and of the FDTD computation differ. Note that we did not observe any diffraction peaks from ν = 1300 cm -1 to ν = 1900 cm -1 due to the presence of water vapor. Lifetime estimation of these surface modes can be obtained by analyzing the width of the grating emission peaks. Considering the experimental dispersion to be photon-like, the coherence length can be deduced as fol-lows:

L = 1 2∆ν , (3)
where ∆ν is the full wave number width at half maximum of the emission peak in cm -1 . Fig. 4(b) shows that the typical coherence length is in the order of 100 µm, that is, ten times larger than the typical coherence length of these surface waves in a semi-infinite SiO 2 -air interface. These values are in agreement with theoretical predictions for a 1 µm thick SiO 2 suspended film reported by Ordonez-Miranda et al. 5 . The obtained length range is also ten times larger than the typical coherence length of the surface waves of a semi-infinite SiO 2 -air interface 1,5 . The coherence length L reaches 700 µm for Zenneck modes which is almost two orders of magnitude larger than their wavelength. Such a large coherence length is achieved due to the fact that most of the electromagnetic energy is propagating in the air close to the interface of the dielectric rather than in the material, decreasing the absorbed power and enhancing the coherence length. Note, that we underestimate the coherence length of the modes since the grating obviously introduces radiative losses and that the values are expected to be even larger for smooth thin films.

In this work we have observed through experiment thermally excited surface waves at the surface of a thin SiO 2 film deposited on aluminum from 882 cm -1 to 3725 cm -1 , whereas an interface between two semiinfinite materials only supports surface waves from 1072 cm -1 to 1156 cm -1 . This spectral broadening is the result of Zenneck and subwavelength TM guided surface waves excitation from 882 cm -1 to 1072 cm -1 and from 1979 cm -1 to 3725 cm -1 in addition to Surface Phonon-Polaritons. From their emission spectra, we were able to reconstruct their dispersion relation and to measure the coherence length of these waves. For Zenneck surface waves, it reaches almost 700 µm which is two orders of magnitude larger than their wavelength. We believe that because of their large coherence length as well as their existence in a very large spectrum, these surface waves can be considered for a wide range of applications in the infrared for both the thermal management of submicron structures and photonics engineering due to their dual nature. Here we focused our experimental study on SiO 2 since it is a very common material in microelectronics but the same phenomena will exist for any dielectric material supporting resonant surface waves.

We wish to thank Jose Ordonez-Miranda, Stéphane Collin, Laurent Tranchant and Mikyung Lim for their fruitful discussions. We also want to acknowledge Thuy-Anh Nguyen, Paul Debue and Romaric de Lépinau for their contribution to the early development of the FDTD code. This work was supported by Renatech project "Phonons Polaritons de Surface large bande".

FIG. 1 .

 1 FIG. 1. (a) Comparison of the FDTD-calculated dispersion relations of a 1 µm thick suspended SiO2 membrane (red line) and of a 0.5 µm thick SiO2 film deposited on aluminum (dashed blue line). (b) Real (magenta line) and imaginary (green line) parts of the relative permittivity of SiO2. Grey region indicates the frequency range of SPhPs.

FIG. 2 .

 2 FIG. 2. Sample design (a) and its emission signal (b) at Ts = 673 K and tilt angle α = 2.6 • . Black arrows indicate the emission peaks due to the diffraction grating.

FIG. 4 .

 4 FIG. 4. Dispersion relation (a) and coherence length (b) of the surface waves for the thin SiO2 film deposited on aluminum obtained from FDTD direct computations (blue line) and experimental measurements (red circles).


Frédéric Hecht

O Pironneauan

Olivier Pironneau
email: olivier.pironneau@upmc.fr

An Energy Stable Monolithic Eulerian Fluid-Structure Finite Element Method

Keywords: Finite Element, Fluid-Structure, Energy Stable, Monolithic scheme, Characteristics-Galerkin, Eulerian Formulation

When written in an Eulerian frame, the conservation laws of continuum 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 challenging problem because of moving geometries. In this article the method is presented; time implicit discretizations are proposed with iterative algorithms 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 Characterics-Galerkin method coupled with a powerful mesh generator. 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.

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 [START_REF] Formaggia | Cardiovasuclar Mathematics[END_REF] and the references therein). For large displacements the difficulty is transferred to the mesh [START_REF] Liu | A second-order changing-connectivity ALE scheme and its application to FSI with large convection of fluids and near-contact of structures[END_REF] [START_REF] Th | Fully Eulerian fluidstructure interaction for time-dependent problems[END_REF] and to a lesser extent to the matching conditions at the fluid-solid interface [START_REF] Tallec | Energy conservation in fluid-structure interactions[END_REF]. 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 [START_REF] Fernandez | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF] [START_REF] Bukaca | Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement[END_REF]. Cited references here are certainly not exhaustive as an enormous amount of work has gone into ALE methods.

Immersed boundary methods (IBM) [START_REF] Peskin | The immersed boundary method[END_REF] 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 [START_REF] Th | Implicit Boundary and Adaptive Anisotropic Meshes[END_REF], yet it is a very flexible method which has been shown stable -for a given deformation map -by Boffi et al [START_REF] Boffi | The finite element immersed boundary method with distributed Lagrange multiplier[END_REF].

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 x = X(x 0 , t), to restate the equations for the structures in terms of the x 0 -gradient of its displacement d(x 0 , t) := X(x 0 , t) -x 0 defined in the initial domain x 0 ∈ Ω s 0 . In an Eulerian formulation one frames the solids by its velocity u(x, t) and displacement d(x, t) = x -X -1 (x, t) in the moving domain x ∈ Ω s t with x-gradients of these; here X -1 is the inverse map of x 0 → x = X(x 0 , t).

These facts have been highlighted particularly in [START_REF] Pironneau | Fluid structure interaction with large structural displacements[END_REF] (and perhaps also in [START_REF] Bathe | Finite element procedures[END_REF]) and even more so by Stefan Turek in [START_REF] Hron | A monolithic fem solver for an ALE formulation of fluid-structure interaction with configuration for numerical benchmarking[END_REF] (yet still using d(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 [START_REF] Th | Adaptive Finite Element Approximation Of Fluid-Structure Interaction Based On An Eulerian Variational Formulation[END_REF] [START_REF] Th | An Eulerian approach to fluid-structure interaction and goaloriented mesh adaptation[END_REF]), 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 [START_REF] Th | Rannacher Adaptive Finite Element Approximation of Fluid-Structure Interaction Based on an Eulerian Variational Formulation[END_REF] (hierarchical meshes) [START_REF] Rannacher | An Adaptive Finite Element Method for Fluid-Structure Interaction Problems Based on a Fully Eulerian Formulation[END_REF] (adaptive meshes), Thomas Wick [START_REF] Th | Fully Eulerian fluidstructure interaction for time-dependent problems[END_REF] (comparison with ALE), Thomas Richter [START_REF] Th | A Fully Eulerian formulation for fluidstructure interaction problems[END_REF] [START_REF] Th | Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates[END_REF] ((an XFEM-like strategy to capture the interface and also an immersed boundary-like strategy and fixed meshes) and Stefan Frei [START_REF] Frei | Eulerian finite element methods for interface problems and fluidstructure interactions[END_REF] (detailed study of contacts) we will refer to them as of the Heidelberg school.

While embarking on Eulerian formulations, in [START_REF] Pironneau | Numerical Study of a Monolithic Fluid-Structure Formulation[END_REF], 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 [START_REF] Th | Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates[END_REF]. There are minor differences in the continuous formulation (no inverse matrix of ∇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 [START_REF] Th | A Fully Eulerian formulation for fluidstructure interaction problems[END_REF] and an XFEM-like method in [START_REF] Th | Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates[END_REF], 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 [START_REF] Pironneau | An Eulerian Monolithic Fluid-Structure Formulation[END_REF] and to three dimensions [START_REF] Chiang | Numerical Study of a 3D Eulerian Monolithic Formulation for Fluid-Structure-Interaction[END_REF].

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 velicity-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 Ω s t . A key point in section 3 for the fully implicit in time discretization of system [START_REF] Boffi | The finite element immersed boundary method with distributed Lagrange multiplier[END_REF] is the approximation of the purely convective space-time PDEs by finitedifference in time approximations of the material derivatives. In variational form this is known as the Characteristics-Galerkin method [START_REF] Pironneau | Finite Element Methods for Fluids[END_REF]. Henceforth d n+1 can be eliminated in favor of u n+1 so as to obtain a monolithic FSI system which involves a solver in terms of u n+1 and p n+1 only, and in the full domain; this is [START_REF] Boulakia | Existence of weak solutions for the motion of an elastic structure in an incompressible viscous fluid[END_REF]. Now there is a delicate problem with PDEs written on time dependent domains:

d n+1 (x) (resp. d n (x)) makes sense only if x ∈ Ω s n+1 (resp. x ∈ Ω s n )
. This forces the following definition for Ω s n+1 : it is the set of x such that xu n+1 (x)δt ∈ Ω s n . This is both a difficulty and an advantage. It requires an iterative solution of ( 7) with a readjustment of Ω s n+1 -we propose in (9) a simple fixed point iteration scheme-but it is also an asset because, evidently, Ω s n+1 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 1 3 -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 approximation. 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 1 3 -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 [START_REF] Pironneau | Numerical Study of a Monolithic Fluid-Structure Formulation[END_REF] 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 [START_REF] Rannacher | An Adaptive Finite Element Method for Fluid-Structure Interaction Problems Based on a Fully Eulerian Formulation[END_REF] have been done for compressible solids. We have performed two cases proposed by Dunne and one by Frei [START_REF] Frei | Eulerian finite element methods for interface problems and fluidstructure interactions[END_REF] 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.

Continuum Mechanics: notations

Consider two non overlapping time dependent bounded sets of R 2 , Ω f t and Ω s t for the fluid and solid respectively, t ∈ (0, T ). Let the full domain be Ω t = Ω f t ∪ Ω s t . The fluid-structure interface is denoted Σ t = Ω f t ∩ Ω s t and the boundary of Ω t is ∂Ω t . Later we will denote Γ the part of ∂Ω t where either the structure is clamped or the fluid does not slip (see Figure 1) At initial time Ω f 0 and Ω s 0 are prescribed.

The following notations are standard in continuum mechanics (e;g. [START_REF] Marsden | Hughes Mathematical Foundations of Elasticity[END_REF], [START_REF] Bathe | Finite element procedures[END_REF], [START_REF] Antman | Nonlinear Problems of Elasticity[END_REF]):

Ω f t Γ Ω s t Σ t
Figure 1: Sketch of the fluid -structure domain and notations -X : Ω 0 × (0, T ) → X(x 0 , t) ∈ Ω t , the Lagrangian position at t of x 0 ; in other words, a material point at x 0 at t = 0 is tranported to X at time t by the motion.

-u = ∂ t X, the velocity of the motion,

-F = ∇ T x0 X := (∇ x0 X) T = ((∂ x0 j X i ))
, the transposed gradient of the motion. We follow [START_REF] Belytschko | Elkhodary Nonlinear Finite elements For Continua and structures[END_REF]; however in some books on continuum mechanics the gradient is the contravariant gradient and the transposed operator is absent.

-d = X(x 0 , t) -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 ρ(x, t) be the density and σ(x, t) the stress tensor. at position x and time t. Finally and unless specified otherwise, all spatial derivatives are with respect to x ∈ Ω t and not with respect to x 0 ∈ Ω 0 .

We recall a few properties:

1. If φ is a function of x = X(x 0 , t), x 0 ∈ Ω 0 , ∇ x0 φ = ((∂ x0 i φ)) = j ((∂ x0 i X j ∂ xj φ)) = F T ∇φ.
2. When X is one-to-one and invertible, d and F can be seen as functions of (x, t) instead of (x 0 , t). They are related by

F T = ∇ x0 X = ∇ x0 (d + x 0 ) = ∇ x0 d + I = F T ∇d + I, ⇒ F = (I -∇d) -T (1)
3. Time derivatives are related by (notice the notation for the material derivative)

D t φ := d dt φ(X(x 0 , t), t) = ∂ t φ(x, t) + u • ∇φ(x, t).
Finally, it is convenient to introduce the notations:

J = det F , Du = ∇u + ∇ T u.
Conservation of momentum and conservation of mass take the same form for the fluid and the solid:

ρD t u = f + ∇ • σ, d dt (Jρ) = 0, ( 2
)
where f is the density of volumic forces. So Jρ = ρ 0 at all times and, consequently, if the medium is incompressible, ρ is constant in the fluid and constant in the solid if it was constant at t = 0. We denote these constants by ρ f and ρ s and at any point x and time t, the density is ρ(

x, t) = ρ f 1 Ω f t (x, t) + ρ s 1 Ω s t (x, t), where 1 Ω (x) = 1 if x ∈ Ω and 0 otherwise.
Continuity of u and of σ • n at the fluid-structure interface Σ (with 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 X (see [START_REF] Ciarlet | Mathematical Elasticity[END_REF], [START_REF] Marsden | Hughes Mathematical Foundations of Elasticity[END_REF]), for instance X must be invertible.

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

σ f := σ |Ω f t = -p f I + µDu, σ s := σ |Ω s t = -p s I + ∂ F Ψ F T (3)
The Helmholtz potential Ψ for a Mooney-Rivlin material is given by

Ψ(F) = c 1 (tr F T F -2) + c 2 (tr (F T F) 2 -tr 2 F T F -2).
where c 1 , c 2 are constant (see [START_REF] Ciarlet | Mathematical Elasticity[END_REF]) . Consequently

∂ F Ψ F T = (2c 1 -4c 2 tr F T F )FF T + 4c 2 (FF 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 I by the Cayley-Hamilton theorem. Now by making use of (1) one obtains -for some scalars α, α , functions of x and t,

∂ F Ψ F T = 2c 1 (I -∇d) -T (I -∇d) -1 + αI = 2c 1 (Dd -∇d∇ T d) + α I. (4)
We now have a complete set of equations [START_REF] Babuska | The Plate Paradox for Hard and Soft Simple Support SIAM[END_REF][START_REF] Bathe | Finite element procedures[END_REF][START_REF] Belytschko | Elkhodary Nonlinear Finite elements For Continua and structures[END_REF] in which the initial domain Ω s 0 does not appear.

Eulerian Variational Formulation in 2D

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

Given Ω f 0 , Ω s 0 and d, u at t = 0, find (u, p, d, Ω f t , Ω s t ) with u |Γ = 0 and

Ω f t ∪Ω s t ρD t u • û -p∇ • û -p∇ • u + Ω f t µ 2 Du : Dû + Ω s t c 1 (Dd -∇d∇ T d) : Dû = Ω f t ∪Ω s t f • û, (5)
for all (û, p) with û|Γ = 0, where d, Ω s t and Ω f t are defined forward in time by

D t d = u and dχ dτ = u(χ(τ ), τ ), χ(t) ∈ Ω r t ⇒ χ(τ ) ∈ Ω r τ ∀τ ∈ (0, t), r = s, f (6)
We have used the notation B :

C = tr B T C . Relation (6) defines Ω s t , Ω f t forward in time.
For simplicity we shall consider only the case of homogeneous boundary conditions on Γ ⊂ ∂Ω, i.e. clamped or no-slip, and homogeneous Neumann conditions on ∂Ω t \Γ.

Remark 1 A similar formulation has been used by Dunne in [START_REF] Th | An Eulerian approach to fluid-structure interaction and goaloriented mesh adaptation[END_REF], eq.( 9) (see also [START_REF] Th | Adaptive Finite Element Approximation Of Fluid-Structure Interaction Based On An Eulerian Variational Formulation[END_REF]), 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 d from the time discretized system at each time step.

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 u n+1 , p n+1 , Ω n+1 with

Ωn+1 ρ n+1 u n+1 -u n • Y n+1 δt • û -p n+1 ∇ • û -p∇ • u n+1 + Ω f n+1 µ 2 Du n+1 : Dû +c 1 Ω s n+1 [D( dn + δtu n+1 ) -∇( dn + δtu n+1 )∇ T ( dn + δtu n+1 )] : Dû = Ωn+1 f • û, Ω n+1 = {x : Y n+1 (x) := x -δtu n+1 (x) ∈ Ω n } (7)
where the stands for the composition with Y n+1 , i.e. dn = d n • Y n+1 , and where d is updated by

d n+1 = dn + δtu n+1 . (8)
We defer the proof to a dedicated section, after showing how the problem can be solved numerically.

Iterative Solution by Fixed Point

In ( 5), the bilinear part on Ω f t 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(δt 2 ), we neglect it and solve the nonlinearity due to Y and Ω f , Ω s , by an iterative process which is essentially a fixed point algorithm:

1. Set ρ = ρ n , Ω = Ω n , u = u n , Y(x) = x -δtu.
2. Solve, for all û, p with û|Γ = 0,

Ω [ρ u n+1 -ũn δt • û -p n+1 ∇ • û -p∇ • u n+1 ] + Ω f µ 2 Du n+1 : Dû + Ω s δtc 1 (Du n+1 -∇ dn ∇ T u n+1 -∇u n+1 ∇ T dn ) : Dû = Ω f • û - Ω s c 1 (D dn -∇ dn ∇ T dn ) : Dû (9) 3. Set u = u n+1 , Y(x) = x -δtu, Ω r = Y -1 (Ω r n ), r = s, f , update ρ by Ω. 4. If not converged return to Step 2.
Remark 2 It will be seen below that, after discretization in space, setting Y(x) = x -δtu, Ω r = Y -1 (Ω r n ), r = s, f means to move each vertex q of the triangulation of Ω s to q + δtu(d) and adapt Ω f to be compatible with Ω 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 u n+1 in Ω s n+1 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 Du will dominate ∇u∇d and so the problem will be well posed for sufficiently small ∇d.

Proof of Theorem 1

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

Discretization of Material Derivatives

Given t ∈ (0, T ) and x ∈ Ω, let χ t u,x (τ ) be the solution at time τ of χ(τ ) = u(χ(τ ), τ ) , ∀τ ∈ (0, t), with χ(t) = x.

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

D t v(x, t) := d dτ v(χ(τ ), τ )| τ =t = ∂ t v + u • ∇v.
Given a time step δt, let us approximate χ by Y:

χ (n+1)δt u n+1 ,x (nδt) ≈ Y n+1 (x) := x -u n+1 (x)δt
Remark 4 Note that

J n+1 := det ∇Y n+1 = 1 -∇ • u n+1 δt + det ∇u n+1 δt 2 . ( 10
)
So, in the incompressible case, one can neglect the Jacobians and still be O(δt)consistent.

Variational Form of the Time Discrete Problem

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

Find u n+1 ∈ H 1 0 (Ω n+1 ), p ∈ L 2 (Ω n+1 ), Ω f n+1 , Ω s n+1 , such that with Ω n+1 = Ω f n+1 ∪ Ω s n+1 , ∀û ∈ H 1 0 (Ω n+1 ), ∀p ∈ L 2 (Ω n+1
) the 3 relations below hold:

Ωn+1 ρ n+1 u n+1 -u n • Y n+1 δt • û -p n+1 ∇ • û -p∇ • u n+1 + Ω f n+1 µ 2 Du n+1 : Dû + c 1 Ω s n+1 [Dd n+1 -∇d n+1 ∇ T d n+1 ] : Dû = Ωn+1 f • û, (11)
Ω n+1 = (Y n+1 ) -1 (Ω n ) = {x : Y n+1 (x) := x -δtu n+1 (x) ∈ Ω n } (12) d n+1 = d n •Y n+1 + δtu n+1 . ( 13
) Remark 5 Note that ρ = ρ f 1 Ω f t + ρ s 1 Ω s t is convected by the velocity u. Denot- ing ρ n (x) := ρ(x, nδt), ρ n • Y n+1 (x) = ρ n+1 (x), x ∈ Ω n+1
So discretizing the material derivative of u or the material derivative ρu gives the same scheme because in both cases

1 δt v n+1 (x) -v n (Y n+1 (x)) = (∂ t v+u∇v)| x,t n+1 +O(δt) with v = u or v = ρu
Now to solve [START_REF] Chiang | Numerical Study of a 3D Eulerian Monolithic Formulation for Fluid-Structure-Interaction[END_REF][START_REF] Th | Implicit Boundary and Adaptive Anisotropic Meshes[END_REF][START_REF] Coutand | Motion of an elastic solid inside an incompressible viscous fluid[END_REF], let us plug ( 13) into [START_REF] Chiang | Numerical Study of a 3D Eulerian Monolithic Formulation for Fluid-Structure-Interaction[END_REF] to obtain the "monolithic" formulation of Theorem 1, i.e. using,

∇d n+1 = ∇ dn + δt∇u n+1
Thus (7) is proved.

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

d dt Ωt ρ 2 |u| 2 + µ 2 Ω f t |Du| 2 + d dt Ω s 0 Ψ(I + ∇ T x0 d) = Ωt f • u
Proof. Choosing û = u, p = -p will give the proposition provided

2c 1 Ω s t (Dd -∇d∇ T d) : DD t d = d dt Ω s 0 Ψ(∇ x0 X).

By construction

Ω s t c 1 (Dd -∇d∇ T d) : Dû = Ω s t (∂ F Ψ(F)F T -αI) : Dû = Ω s 0 ∂ F Ψ(F) : D x0 û Now as d dt Ψ(F) = ∂ F Ψ(F) : ∂ t F and ∇ x0 u(x 0 ) = ∂ t ∇ x0 d(x 0 ) = ∂ t F T (x 0 ), Ω s t 2c 1 (Dd -∇d∇ T d) : ∇u = Ω s 0 d dt Ψ(F) = d dt Ω s 0 Ψ(I + ∇ T x0 d).
The other terms are standard, in particular, with enough regularity on Ω r t ,

Ω r t (∂ t u + u • ∇u) • u = d dt Ω r t |u| 2 2 , r = f, s
Remark 6 When Ψ 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 [START_REF] Coutand | Motion of an elastic solid inside an incompressible viscous fluid[END_REF], later improved by Boulakia [START_REF] Boulakia | Existence of weak solutions for the motion of an elastic structure in an incompressible viscous fluid[END_REF] and extended in a special case in [START_REF] Raymond | Vanninathan A fluid-structure model coupling the Navier-Stokes equations and the Lamé system[END_REF].

Stability of the Scheme Discretized in Time

Theorem 2 When f = 0 and ρ is constant in each domain Ω r n , r = s, f ,

Ωn ρ n 2 |u n | 2 + δt n k=1 Ω f k µ 2 |Du k | 2 + Ω s 0 Ψ n ≤ Ω0 ρ 0 2 |u 0 | 2 + Ω s 0 Ψ 0 (14)
We begin the proof with the following lemmas, Lemma 1 The Lagrangian map X n : Ω 0 → Ω n , implicit from [START_REF] Th | Implicit Boundary and Adaptive Anisotropic Meshes[END_REF], satisfies

X n+1 = (Y n+1 ) -1 • X n , n ≥ 1 and F n := ∇ T x0 X n = (I -∇d n ) -T . Proof Notice that Y 1 (Y 2 (..Y n-1 (Y n (Ω n ))..)) = Ω 0 Hence X n+1 = [Y 1 (Y 2 (..Y n (Y n+1 )))] -1 = (Y n+1 ) -1 • X n .
By definition of d n+1 in (13),

d n+1 (X n+1 (x 0 )) = d n (Y n+1 (X n+1 (x 0 ))) + δtu n+1 (X n+1 (x 0 )) = d n (X n (x 0 )) + δtu n+1 (X n+1 (x 0 )), (15)
so X n+1 (x 0 ) = d n+1 (X n+1 (x 0 )) + x 0 and therefore

F n+1 = ∇ T x0 (d n+1 ((X n+1 (x 0 ))) + x 0 ), = ∇d n+1 T F n+1 + I ⇒ F n+1 = (I -∇d n+1 ) -T (16)
Lemma 2 For some α(x, t),

Dd n+1 -∇d n+1 ∇ T d n+1 = J -2 n+1 [(I -∇d n+1 ) -T (I -∇d n+1 ) -1 + αI
Proof With β = det 2 I-∇d n+1 , and γ the trace of (I -∇d n+1 )(I -∇d n+1 ) T , by the Cayley-Hamilton theorem, we have

Dd n+1 -∇d n+1 ∇ T d n+1 -I = -(I -∇d n+1 )(I -∇d n+1 ) T = β (I -∇d n+1 )(I -∇d n+1 ) T -1 -γI Lemma 3 Ω s n+1 c 1 [J -1 n+1 (I -∇d n+1 ) -T (I -∇d n+1 ) -1 ] : Dû = Ω s 0 ∂ F Ψ n+1 : ∇ x0 û Proof From Lemma 1, Ω s n+1 c 1 [J -1 n+1 [(I -∇d n+1 ) -T (I -∇d n+1 ) -1 ] : Dû = Ω s n+1 c 1 J -1 n+1 [F n+1 F n+1 T ] : Dû = Ω s 0 c 1 F n+1 : D x0 û = 1 2 Ω s 0 ∂ F Ψ n+1 : D x0 û = Ω s 0 ∂ F Ψ n+1 : ∇ x0 û (17)

Proof of Theorem 2

Let r = s or f . Let us choose û = u n+1 in [START_REF] Chiang | Numerical Study of a 3D Eulerian Monolithic Formulation for Fluid-Structure-Interaction[END_REF][START_REF] Th | Implicit Boundary and Adaptive Anisotropic Meshes[END_REF][START_REF] Coutand | Motion of an elastic solid inside an incompressible viscous fluid[END_REF]. By Schwartz inequality

Ω r n+1 (ρ r n u n )•Y n+1 •u n+1 ≤ Ω r n+1 ( √ ρ r n u n ) 2 • Y n+1 ) 1 2 Ω r n+1 ( √ ρ r n+1 u n+1 ) 2 1 2 because ρ r n • Y n+1 (x) = ρ r n+1 (x), x ∈ Ω r n+1 .
So by a change of variable

Ω r n+1 ρ r n+1 (u n • Y n+1 ) 2 = Ω r n+1 ( ρ r n u n ) 2 • Y n+1 = Ω r n ρ r n u n2 .
Consequently, using ab

≤ 1 2 a 2 + 1 2 b 2 , Ω r n+1 ρ r n+1 u n • Y n+1 • u n+1 ≤ 1 2 Ω r n ρ r n u n2 + 1 2 Ω r n+1 ρ r n+1 u n+1 2 .
Finally,

Ωn+1 ρ n+1 2 |u n+1 | 2 + δt Ω f n+1 µ 2 |Du n+1 | 2 + Ω0 Ψ n+1 ≤ Ωn ρ n 2 |u n | 2 + Ω0 Ψ n . ( 18
)
Remark 7 Inequality [START_REF] Hecht | New development in FreeFem++[END_REF] 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 δ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.

Spatial Discretization with Finite Elements

Let T 0 h 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 Σ; therefore Q h is the space of piecewise linear functions on the triangulation, continuous in Ω r n+1 , r = s, f ; V 0h is the subspace of V h of functions which are 0 on Γ. A small penalization with parameter must be added to impose uniqueness of the pressure.

Discretization in space by the Finite Element Method leads to find

u n+1 h , p n+1 h ∈ V 0h × Q h such that for all ûh , ph ∈ V 0h × Q h , Ωn+1 ρ n+1 u n+1 h -u n h • Y n+1 δt • ûh -p n+1 h ∇ • ûh -ph ∇ • u n+1 h +1 Ω f n+1 µ 2 Du n+1 h : Dû h +c 1 1 Ω s n+1 [D( dn h + δtu n+1 h ) -∇( dn h + δtu n+1 h )∇ T ( dn h + δtu n+1 h )] : Dû h = Ωn+1 f • ûh , , Ω n+1 = (Y n+1 ) -1 (Ω n ) = {x : Y n+1 (x) ∈ Ω n } (19
Y n+1 (x) = x -δtu n+1 h (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

X n = X n+1 • Y n+1 . ( 20
)
T k n T k 0 X n T k n+1 X n+1 Y n+1 q j n q j n+1 T k 0 T k n T k n+1 X n Y n+1 X n+1 q j n q j n+1 Figure 2: Sketch to understand if X n = Y n+1 • 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 1 3 -P 1 element where each triangle is divided into four subtriangles on which the velocities are P 1 and continuous. A triangle T k 0 in the reference domain (chosen here to be its initial position at time zero) becomes the isoparametric triangles T k n and T k n+1 :

T k n = X n (T k 0 ) and T k n+1 = X n+1 (T k 0 )
in the case of quadratic velocities. But as the vertices and mid-edges of T k n are obtained by moving those of T k n+1 by -δtu n+1 we cannot have Y n+1 (T k n+1 ) = T k n 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.

First let us consider the P 2 -P 1 element for V 0h × Q h . Then the mapping x → Y n+1 (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 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 [START_REF] Frei | Long-term simulation of large deformation, mechano-chemical fluid-structure interactions in ALE and fully Eulerian coordinates[END_REF] for the obvious reason that f • g is of degree 4 when f and g are quadratic.

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

Then [START_REF] Frei | Long-term simulation of large deformation, mechano-chemical fluid-structure interactions in ALE and fully Eulerian coordinates[END_REF] 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 Y n+1 but X n+1 • Y n+1 remains linear and for each triangle

T k n = Y n+1 (T k n+1
) 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 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 u n is defined on Ω n and needs to be interpolated when u n (x -δtu n+1 ) is computes with x ∈ Ω 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 1 3 -P 1 , as in [START_REF] Th | Adaptive Finite Element Approximation Of Fluid-Structure Interaction Based On An Eulerian Variational Formulation[END_REF][START_REF] Th | A Fully Eulerian formulation for fluidstructure interaction problems[END_REF] and others. Then a stabilization term ω∇ • u n+1 ∇ • v is added to the first integral in the variational formulation [START_REF] Frei | Eulerian finite element methods for interface problems and fluidstructure interactions[END_REF] where ω is a mesh dependent parameter. The proof of Theorem 2 carries over except that

Ω f n+1 µ 2 |Du n+1 | 2 becomes Ω f n+1 µ 2 |Du n+1 | 2 + Ωn+1 ω(∇ • u n+1 ) 2 .

Conclusion

In the case of solid only, Theorem 2 holds for the fully discrete problem [START_REF] Frei | Eulerian finite element methods for interface problems and fluidstructure interactions[END_REF] 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.

Contacts

To each simply connected disjointed part S i t , i = 1..n s of Σ t or Γ t is associated a signed distance function x → d S i t (x) measuring the Euclidian distance of x to S i t with the sign indicating whether x is in the structure (d S i t (x) < 0) or in the fluid (d S i t (x) > 0) . When d S i t (x) = 0 for some x ∈ ∂Ω\S i t there is a contact. Note that contacts between points on the same part S 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 λ:

Ωt ρD t u • û -p∇ • û -p∇ • u + Ω f t µ 2 Du : Dû + Ω s t c 1 (Dd -∇d∇ T d) : Dû + ∂Ωt λn • û = Ωt f • û, λ(x, t) ≤ 0, d S i t (x)λ(x, t) = 0 ∀x ∈ S i t , ∀t. (21)
To solve the problem we apply the discretization described above and the Semi-Smooth Newton method proposed by Ito and Kunisch [START_REF] Ito | Semi-Smooth Newton Method for Variational Inequalities of the First Kind[END_REF] which replaces inequality constraints by equality constraints at each time step; here it is:

λ n+1 (x) -min{0, λ n+1 (x) + c 0 d S i n+1 (x + δtu n+1 n (x)n(x))} = 0 ∀x ∈ S i n+1
and similarly on Γ; 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,

Ω n+1,k ρ n+1 u n+1,k+1 h -ũn h δt • ûh -p n+1 h ∇ • ûh -ph ∇ • u n+1 h + 1 Ω f n+1 µ f 2 Du n+1,k+1 h : Dû h +c 1 1 Ω s n+1 [D( dn h + δtu n+1,k+1 h ) -∇( dn h + δtu n+1,k h )∇ T ( dn h + δtu n+1,k+1 h )] : Dû h + ns 1 S i n+1 ∩{x:λ+c0d S i n+1 (x+δtu n+1,k h (x))<0} c 3 (x + δt(n • u h n+1,k+1 )n)û h = Ωn+1 f • ûh , Ω n+1,k+1 = {x : Y n+1 (x) := x -δtu n+1,k+1 h ∈ Ω n } ( 22
)
where c 3 is a very large constant. The purpose of c 3 is to impose x+δtu n n+1 h n(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 [START_REF] Th | A Fully Eulerian formulation for fluidstructure interaction problems[END_REF]. 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 .

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++ [START_REF] Hecht | New development in FreeFem++[END_REF] 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.

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 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 θ - 1 r ∂ r [ξ f r∂ r u θ + ξ s r∂ r d] = 0, ∂ t d = u θ , r ∈ (R 0 , R 1 ), v |R0 = 0, u θ |R1 = 3, (23)
with ρ = ρ s 1 r≤R + ρ f 1 r>R , ξ s = 2c 1 1 r≤R , ξ f = µ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. 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.

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 d n+1 = d n + δtu n+1 instead of [START_REF] Bukaca | Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement[END_REF].

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 [START_REF] Babuska | The Plate Paradox for Hard and Soft Simple Support SIAM[END_REF] (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.

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 [START_REF] Th | Rannacher Adaptive Finite Element Approximation of Fluid-Structure Interaction Based on an Eulerian Variational Formulation[END_REF]; 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 × 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 × 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. Solid ρ s = ρ f , c 1 = 10 6 kgm -1 s -1 and no external force. Initially all velocities and displacements are zero. Around t ∼ 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 [START_REF] Th | Adaptive Finite Element Approximation Of Fluid-Structure Interaction Based On An Eulerian Variational Formulation[END_REF] . The results concur roughly with [START_REF] Th | Rannacher Adaptive Finite Element Approximation of Fluid-Structure Interaction Based on an Eulerian Variational Formulation[END_REF] 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 [START_REF] Th | Rannacher Adaptive Finite Element Approximation of Fluid-Structure Interaction Based on an Eulerian Variational Formulation[END_REF] 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 [START_REF] Pironneau | An Eulerian Monolithic Fluid-Structure Formulation[END_REF]) 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.

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 [START_REF] Th | Rannacher Adaptive Finite Element Approximation of Fluid-Structure Interaction Based on an Eulerian Variational Formulation[END_REF].

In [START_REF] Th | Rannacher Adaptive Finite Element Approximation of Fluid-Structure Interaction Based on an Eulerian Variational Formulation[END_REF] 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. 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 - Ω n+1 ρn+1 2 |u n+1 | 2 + n 1 δt Ω f j+1 µ 2 |Du j+1 | 2 + Ω 0 Ψ n+1 / Ω 0 Ψ 0 - n+1 1 δt Ω j gρju j 2 (24)
Three elements have been tested: P 2 -P 1 , P 1 bubble -P 1 and P 1 -P 1 without stabilization (it seems that using D instead of ∇ 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 [START_REF] Pironneau | Numerical Study of a Monolithic Fluid-Structure Formulation[END_REF][Figure 5] and show no improvement over 2). 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.

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. We have also studied the P 1 3 -P 1 element with moving mid vertices; small differences are observed but not worth the switch from P 1 bubble-P 1 to P 1 3 -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 1 3 -P 1 element.

Figure 10: The P 1 3 -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.

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 1 -P 1 element instead of P 1 3 -P 1 element.

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 [START_REF] Th | A Fully Eulerian formulation for fluidstructure interaction problems[END_REF] and Stefan Frei in [START_REF] Frei | Eulerian finite element methods for interface problems and fluidstructure interactions[END_REF].

A circular structure of radius R = 0.4, centered at x c = 1, y c = 1, in a rectangular cavity [0, 2] × [0, 1.5] filled with a fluid of density ρ f = 1000 and viscosity µ = 1000. The lateral walls of the cavity are numerical boundaries where no condition is imposed in the variational formulation, i.e. σ f n = 0.

The Mooney-Rivlin coefficient is c 1 = 10 5 and ρ 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 [START_REF] Frei | Eulerian finite element methods for interface problems and fluidstructure interactions[END_REF][Figure 12.4 page 138] (see also [START_REF] Th | A Fully Eulerian formulation for fluidstructure interaction problems[END_REF] for another case and [START_REF] Frei | Long-term simulation of large deformation, mechano-chemical fluid-structure interactions in ALE and fully Eulerian coordinates[END_REF] 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. 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 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++ [START_REF] Hecht | New development in FreeFem++[END_REF] and validation tests have been presented. Validation by comparison with other authors 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 [START_REF] Pironneau | An Eulerian Monolithic Fluid-Structure Formulation[END_REF] and a 3D implementation is being developed [START_REF] Chiang | Numerical Study of a 3D Eulerian Monolithic Formulation for Fluid-Structure-Interaction[END_REF]. Finding a second order in time extension seems much harder.

Figure 3 :

 3 Figure3: 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.

5 Figure 4 :

 54 Figure4: 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 θ = π/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.

2 ;

 2 other parameters are l = 0.35, h = 0.02, L = 2.5, H = 0.41. Fluid The density is ρ f = 10 3 kg/m 3 and the reduced viscosity ν f = µ ρ f = 10 -3 m 2 /s; the inflow velocity is horizontal and parabolic in y := x 2 with flux Ū H, with Ū = 2; namely u(0, y) = Ū 6 H 2 y(H -y), 0 T . Top and bottom boundaries are no-slip walls.

Figure 5 :

 5 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.

Figure 6 :

 6 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 -1 and the amplitude around 0.03. The mesh has 2200 vertices and the time step is 0.005.

Figure 7 :

 7 Figure 7: Final position of the flag at t=0.98 after it has fallen under its own weight.

Figure 8 :

 8 Figure8: 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.

Figure 9 :

 9 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.

Figure 11 :Figure 12 :

 1112 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.

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 1 element for the pressures (one for the fluid part and another for the structural part) and a fancy P 1 element for the velocities on a mesh obtained by dividing