Numerical Study of a Monolithic Fluid-Structure Formulation

The conservation laws of continuum mechanic are naturally written in an Eulerian frame where the difference between a ﬂuid and a solid is only in the expression of the stress tensors, usually with Newton’s hypothesis for the ﬂuids and Helmholtz potentials of energy for hyperelastic solids. There are currently two favored approaches to Fluid Structured Interactions (FSI) both working with the equations for the solid in the initial domain; one uses an ALE formulation for the ﬂuid and the other matches the ﬂuid-structure interfaces using Lagrange multipliers and the immersed boundary method. By contrast the proposed formulation works in the frame of physically deformed solids and proposes a discretization where the structures have large displacements computed in the deformed domain together with the ﬂuid in the same; in such a monolithic formulation velocities of solids and ﬂuids are computed all at once in a single variational formulation by a semi-implicit in time and the ﬁnite element method. Besides the simplicity of the formulation the advantage is a single algorithm for a variety of problems including multi-ﬂuids, free boundaries and FSI. The idea is not new but the progress of mesh generators renders this approach feasible and even reasonably robust. In this article the method and its discretization are presented, stability is discussed showing in a loose fashion were are the difﬁculties and why one is able to show convergence of monolithic algorithms on ﬁxed domains for ﬂuids in compliant shell vessels restricted to small displacements. A numerical section discusses implementation issues and presents a few simple tests.


Introduction
In an earlier paper [6] the author and his coauthors proposed a method to compute a fluid in a vessel modeled as a shell with normal displacements as in Nobile and Vergara [23].It was argued that since the model is valid for small displacements only, one may as well use a transpiration approximation for the fluid and do the full computation in a fixed domain.As we were able to prove convergence, an interesting question arose: what is so special about the model that one could prove existence and convergence of the numerical scheme?
This paper answers partially the question: what makes FSI really hard is the moving domain.The same is true of free boundary problems for the Navier-Stokes equations.So it was the transpiration approximation for the moving part which made the analysis possible in [6].
Turning to ALE to work on a fixed domain both for the fluid and the solid is a popular solution [26], but the difficulty is transferred to the mesh [19] and the matching conditions at the fluid-solid interface [17].Even more so with immersed boundary methods (IBM) [24] [11], although the convergence analysis is more advanced [3].
Furthermore, iterative solvers for FSI which rely on alternative solutions of the fluid and the structure parts are subject to the added mass effect and require special solvers [12] [5].
Every so often it is not a bad idea, I guess, to rethink fundamentals and check that what is taken for granted in numerical analysis is still true in the face of hardware and software progress.
So, is there an alternative to ALE and IBM?One old method [1] has resurfaced recently, the so-called actualized Lagrangian methods for computing structures [16] [20] (see also [10] although different from the present study because it deals mostly with membranes).
Continuum mechanics doesn't distinguish between solids and fluids till it comes to the constitutive equations.This has been exploited in many studies but most often in the context of ALE [18] [27].
In the present study we investigate what Stephan Turek [27] calls a monolithic formulation but here in an Eulerian framework, following the displaced geometry of the fluid and the solid.
To the specialist it may appear to be a back to square one idea and it is true: there is nothing new here from the modeling view point; everyone knows that it can be done.What is new is that the mesh generators are now robust and agile at following complex motions of objects, making feasible an Eulerian numerical method.
The first difficulty with Eulerian methods comes from the hyperbolic character of the equations for the displacement of solids while those for the fluid are parabolic in time for the velocity.So let us begin by showing that a wave equation for a displacement can be reformulated as a seemingly parabolic equation for its velocity.

Preliminaries on the Wave Equation
At the core of the numerical scheme proposed here for large displacements is the following rewriting of the first and second order finite difference in time schemes for the wave equation: For ] the following scheme is unconditionally stable on a uniform finite difference grid in 1D(see for instance [21]): where ), these schemes can be rewritten as and initialized by u 0 = d 0 = 0. Evidently this is also unconditionally stable and first order when α = 0, θ = 1 and second order when

General Laws of Continuum Mechanics
Consider a time dependent computational domain Ω t made of a fluid region Ω f t and a solid region Ω s t : At initial time Ω f 0 and Ω s 0 are prescribed.The following notations are standard [8], [22], [1], [27], [18]: • X : Ω 0 × (0, T ) → Ω t : X(x 0 ,t), the Lagrangian position at t of x 0 .• u = ∂ t X, the velocity of the deformation, • F ji = ∂ x 0 i X j , the transposed gradient of the deformation, • J = det F , the jacobian of the deformation.
Remark 1.We use a notation for the gradient which is the transposed of the one found in engineering anglo-saxon books (see [8] for instance).Here the gradient of a scalar being a column vector the gradient of a vector is the row of vectors of the derivative of its components.
We denote by tr A and det A the trace and determinant of A. As usual the following quantities are introduced: B(x,t) the density of volumic and surfacic forces at x,t.
Finally and unless specified all spatial derivatives are with respect to x ∈ Ω t and not with respect to x 0 ∈ Ω 0 .If φ is a function of x = X(x 0 ,t), x 0 ∈ Ω 0 , 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

Time derivatives are related by
It is convenient to introduce the notation Conservation of momentum and conservation of mass take the same form for the fluid and the solid: So Jρ = ρ 0 at all times and with continuity of u and of σ • n at the fluid-structure interface Σ when B = 0.There are also unwritten constraints pertaining to the realizability of the map X (see [8], [22]).Finally incompressibility implies J = 1 and so ρ = ρ 0 constant.
• For a Newtonian incompressible fluid : where Ψ is the Helmholtz potential which, in the case of a Mooney-Rivlin two dimensional material, is [8] For a compressible Mooney-Rivlin material the same holds but without p s .Here we will only consider incompressible material, but everything said below can be adapted easily.

Computation of the Mooney-Rivlin 2D Stress Tensor
It is readily seen that Now by the Cayley-Hamilton theorem B 2 = cB − bI so By the Cayley-Hamilton theorem again Hence an incompressible two dimensional Mooney-Rivlin material will have, for some α, ∂ F Ψ F T = 2c 1 (Dd − ∇d∇d T ) + αI.

Monolithic Variational Formulation in 2D
Let Γ ⊂ ∂ Ω the part of the boundary on which the solid is clamped or the fluid has a no-slip condition.For incompressible material the final fluid-structure formulation in two dimensions is: for all ( û, p) with û|Γ = 0, where Ω s t and Ω f t are defined incrementally by Initial conditions are: u given, d = 0, Ω r 0 given, r = s, f .
We have used the notation B : C = tr B T C and c1 := ρ s c 1 .

Conservation of Energy
Proposition 1.
When Ψ is convex, an existence of solution result can be gained from this equality (see [14] for example).
Proof.Choosing û = u, p = −p will give the proposition provided

Discretization in Time
It is natural to use the following discretization So if X n is a first order approximation of X(t n+1 − δt) defined by , a consistent first order scheme is to find u n+1 , p n+1 such that u n+1 = 0 on Γ and for all û, p with û|Γ = 0,

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 .A small penalization with parameter ε must be added to impose uniqueness of the pressure when a direct linear solver is used.
At each time step one must find u n+1 h , p n+1 h Then the triangulation must be updated by u n+1 h by moving each vertex from q n i to

Map Preserving Scheme
It is important to understand on which triangulation each variable is defined.One possibilities is to assume that Ω n+1 is constructed by successive iterations such that Hence As above, the Cayley-Hamilton theorem and incompressibility imply that This formula can be used in the variational formulation instead of ∇d.Its numerical performance is similar to the one that uses d but it preserves energy... at the cost of an iterative adjustment of Ω n+1 .

Remark on the Construction of Second order accurate schemes
To build a second order scheme for (3) one would have to • use a second order characteristic method [4]  .
The last item is difficult; the construction of second order schemes is discussed in particular in the thesis of P. Hauret [14] with the Newmark mid-point scheme [15].It seems rather difficult to prove the same here, so we postpone it to a future study.

Perturbation about an Equilibrium
An equilibrium is reached when p, d and Ω s = {x 0 + d(x 0 ) : Assuming for clarity that ρ s = ρ f and Ω is independent of t, if we prime all variations about that state, they much be such that ∂ t d = u and because, provided Σ is smooth, Ω s t φ = Ω s φ + Σ φ u • n.On this formulation we believe that it is not difficult to prove existence, uniqueness and convergence of the time scheme as in [6].

On the Stability of the Scheme
For clarity we drop the subscript h.
Consider the scheme based on (8) and iterated so as to replace Ω n by Ω n+1 Then Now suppose for clarity that ρ f = ρ s and f = 0; the general case will be analyzed later.Choosing p = −p n+1 , û = u n+1 and assuming Ω independent of t, Lemma 1.
Now let us analyze the last term in (12).
By the convexity of Consider Ω r t , r = s or f and 4 Formulation in the initial domain for the Solid We wish to compare the formulation in the moving domain with the formulation in the fixed domain for a single hyperelastic incompressible structure.Recall that the deformation map satisfies Integrated on the initial domain and assuming no additional surface constraint on ∂ Ω 0 , the variational formulation is: for all d zero on Γ .A fully implicit numerical scheme is [27][14]: ) and θ = 1 or 2 and with d n+1 = d n + δtu n+1 .Neglecting the terms in δt 2 , leads to Remark 2. To be consistent with the formulation in the moving domain one ought to have written leading to the formulation Remark 3. The following identity is true for all vector valued functions in H 1 0 (Ω ) 2 (see Costabel et al [9]): It changes the boundary condition on ∂ Ω \Γ to use it in (16), ii.e.
The Mooney-Rivlin model may be degenerate in 2D for incompressible material or the boundary condition implied by this formulation may be more appropriate.This point needs to be studied further.
5 Numerical Tests

Comparisons for an Incompressible Beam
The monolithic method set in the moving domain is compared with the more classical method ( 16) set in the initial domain in the case of a single rectangular elastic beam of length to width ratio equal to 10 and bent by its own weight from rest and clamped either n both side or on one side.The spatial discretization is performed by using Lagrangian Triangular Finite Elements with polynomials of degree 2 for the velocities and displacements and degree 1 for the pressures.FreeFem++ [13] has been used to implement the algorithms.All linear systems are solved by direct solvers of the libraryMUMPS.
The method in the fixed initial domain did not seem to work with (16) but it did with (17); we have used the second order approximation, θ = 1.
The penalization parameter is ε = 0.01 for the formulation in the initial domain and 10 −6 for the one in the moving domain.The other parameters are The gravity differs : 1. when the beam is clamped at both sides, it is set to f = 2. and when the beam is clamped on the left only f = −0.002.
The beam clamped at both side went through one and a half oscillation cycle during the 45 time steps while the one clamped on the left only went through a half cycle.With the method using the initial domain there is a slight change of surface area, going from 10 initially to 10.57 after 45 time iterations while the change is less that 1% for the other one.

Fall of an Hyperelastic Beam Clamped on One Side
The same beam clamped on one side only is allowed to fall under its own weight for 200 iterations.The results obtained by the method set in the moving domain are shown on figure 2. All parameters have the same value except the time step which is 0.5.auxiliary Laplace equation is solved to move the inner vertices at each time step: The vertices of the mesh in the solid part are moved by u n+1 δt by calling movemesh(), the mesh moving function of FreeFem++.
The mesh is moved by calling the FreeFem++ function movemesh().If it flips over a triangle the function adaptmesh() is called.
The second example is the free motion of an hyperelastic solid submitted to its own weight and to the force due to the rotation of the fluid induced by the sliding of the lower horizontal boundary at unit speed.Initially the solid is a disk.The fluid domain is a rectangle of size 10 × 7 and ν = 0.1, ρ f = 0.5 while ρ s = 1.The gravity is −0.2.With the same data two layers of fluids with different densities, the top one having ρ 2 = 2, rotate under the action of the sliding lower boundary.No-slip conditions is applied to the vertical walls.Here when movemesh() flip over a triangle, it is detected with checkmovemesh() and then a call to adaptmesh(th,h,IsMetric=1) rebuild the mesh th with mesh size average h.Once more we stress the fact that all cases have been computed with the same computer program -given in appendix for the case of figure 4 -without modification of the core algorithm.

Conclusion
A fully Eulerian fluid-structure formulation has been presented and an attempt at deriving an implicit unconditionally stable monolithic finite element discretization has been proposed and studied.The method has been implemented with FreeFem++.It is reasonably robust but needs to be made unconditionally stable by implicit iterations on the moving domain, a path currently under investigation.
n+1 h it suffices to copy the array of values of d n h in the array of values of d n+1 h and add δtu n+1 h (q n i ) to the array.The vertices in the fluid are moved by ũ solution of −∆ ũ = 0 in the fluid and ũ = u on Σ and zero on other boundaries.Moving the vertices of T n h gives a new triangulation T n+1 h .

Fig. 1
Fig. 1 Top: Maximum bent under its own weight for the formulation in the initial domain (left) and for the formulation in the moving domain (right).Bottom: Free fall under its own weight of the same solid clamped on the left only; position at time 50 computed in the initial domain (left) and the moving domain (right).For the full swing see figure 4.

Fig. 2
Fig. 2 Free fall under its own weight of the same hyperelastic solid clamped on the left only; position at iteration 50, 100, 150, 200.

Fig. 3
Fig. 3 Free fall of the same beam, clamped on the right, in a fluid initially at rest.Every 10 th steps the geometry and the norm of velocities are shown.At the 110 th time steps the mesh is unusable, so adaptmesh() is called.

Fig. 4 Fig. 5 A
Fig.4Free motion of an hyperelastic incompressible solid submitted to its own weight and to the force due to the rotation of the fluid induced by the sliding lower horizontal boundary.