An Energy Stable Monolithic Eulerian Fluid-Structure Finite Element Method

When written in an Eulerian frame, the conservation laws of continuum mechanics are similar for ﬂuids and solids leading to a single set of variables for a monolithic formulation. Such formulations are well adapted to large displacement ﬂuid-structure conﬁgurations, 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 ﬁnite 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 ﬂuid; 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 [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 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 [26] (and perhaps also in [3]) and even more so by Stefan Turek in [22] (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 [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 ∇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 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 (5) 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 [29].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 (7).Now there is a delicate problem with PDEs written on time dependent domains: . This forces the following definition for Ω s n+1 : it is the set of x such that x − u 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 [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.

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.[24], [3], [1]): 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, , 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.
-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: 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 3. Time derivatives are related by (notice the notation for the material derivative) Finally, it is convenient to introduce the notations: Conservation of momentum and conservation of mass take the same form for the fluid and the solid: 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 ρ( 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 [10], [24]), 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 The Helmholtz potential Ψ for a Mooney-Rivlin material is given by where c 1 , c 2 are constant (see [10]) .Consequently 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, We now have a complete set of equations (2,3,4) 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 for all (û, p) with û|Γ = 0, where d, Ω s t and Ω f t are defined forward in time by We have used the notation B : 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 [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 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 where the stands for the composition with Y n+1 , i.e. dn = d n • Y n+1 , and where d is updated by 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: 2. Solve, for all û, p with û|Γ = 0, 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).
If u is Lipschitz in space and continuous in time the solution exists.The Characteristics-Galerkin method relies on the concept of material derivative: Given a time step δt, let us approximate χ by Y: Remark 4 Note that 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: So discretizing the material derivative of u or the material derivative ρu gives the same scheme because in both cases Now to solve (11,12,13), let us plug ( 13) into (11) to obtain the "monolithic" formulation of Theorem 1, i.e. using, 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 Proof.Choosing û = u, p = −p will give the proposition provided

By construction
The other terms are standard, in particular, with enough regularity on Ω r t , 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 [13], later improved by Boulakia [7] and extended in a special case in [34].

Stability of the Scheme Discretized in Time
Theorem 2 When f = 0 and ρ is constant in each domain Ω r n , r = s, f , We begin the proof with the following lemmas, Lemma 1 The Lagrangian map X n : Ω 0 → Ω n , implicit from (12), satisfies By definition of d n+1 in (13), so X n+1 (x 0 ) = d n+1 (X n+1 (x 0 )) + x 0 and therefore Lemma 2 For some α(x, t), , and γ the trace of (I − ∇d n+1 )(I − ∇d n+1 ) T , by the Cayley-Hamilton theorem, we have

Proof of Theorem 2
Let r = s or f .Let us choose û = u n+1 in (11,12,13).By Schwartz inequality So by a change of variable Consequently, using ab Finally, 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 δ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 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 ; 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 : 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 (20) 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 (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 Y n+1 but X n+1 • Y n+1 remains linear and for each triangle ) 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 [14,35] and others.Then a stabilization term ω∇ • u n+1 ∇ • v is added to the first integral in the variational formulation (19) where ω is a mesh dependent parameter.The proof of Theorem 2 carries over except that

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.

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 λ: 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: 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, 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 [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 .

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.

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

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 × 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.  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 [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.

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.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: 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 [30][Figure 5] and show no improvement over 2).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 [35] and Stefan Frei in [19].
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 [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.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++ [18] 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 [31] and a 3D implementation is being developed [11].Finding a second order in time extension seems much harder.

Figure 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 :
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.

Figure 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 :
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 :
Figure 7: Final position of the flag at t=0.98 after it has fallen under its own weight.

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

Figure 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 :
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.