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

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

A popular approach to numerically match the velocities and stresses at the fluid-solid interfaces is the arbitrary Lagrangian Eulerian formulation (ALE) whereby at every time step the fluid equations are mapped back into the solid domain [5].Difficulties with the mesh used in the fluid part arise in the case of large displacements of the solids [6,7].
The immersed boundary method (IBM) proposed by Charles Peskin [2] and analyzed in [8] is efficient for shells in fluid but is more difficult to implement for thick structures [9].Particles in suspension in a fluid have been simulated by IBM-like methods but mostly for hard particles (see for instance [10,11]).
A third way is the fully Eulerian approach, advocated by Liu et al. [12] and Dunne-Rannacher [13,14] (see also [15,16]).This approach is well suited to problems with large displacements.
In [17,18] a similar monolithic numerical method was proposed for the fluid-structure equations in two dimensions.It was also shown to be energy-stable.For more details on the motivations behind such a radically different approach to FSI, see [18].The basic idea is to write the equations in terms of velocity, pressure, and displacement, as usual but in an Eulerian frame, and then eliminate the displacement in the time-discretized equations; the result is a monolithic formulation which is well-posed at each time step.Furthermore, to obtain energy stability it was suggested in [18] to update the structure with its own velocity and remesh the fluid domain at every time step.
In this paper we keep the same approach but in 3D.The Cayley-Hamilton theorem was used in 2D and it is different in 3D.Hence, the equations for a structure in Eulerian variables need to be derived again.The analysis is performed for hyper-elastic Mooney-Rivlin incompressible material.The result is not very different except that one coefficient of the constitutive relation is now a nonlinear function of the deformation tensor.
The article has four sections.
• The fully Eulerian formulation is derived in Section 1.
• In Section 2 the equations are discretized implicitly in time as in [18,19]; then the displacements are eliminated and a non-linear system for the velocities and pressures remains.Spatial discretization is obtained by using stable finite element spaces like linear pressures and quadratic velocities on a tetrahedral mesh.• In Section 3 it is shown that the energy decays at each time step, which is an indication that the method is robust.• Finally, in Section 4, the method is implemented in 2D for axisymmetric systems and in 3D for general systems.Energy decrease, mass conservation and convergence are analyzed on an axisymmetric case: an elastic torus in a cylindrical canister filled with a fluid at a Reynolds number of a few hundred.Then, the method is carefully evaluated on the test case proposed in [20] and comparisons with previous publications are made.

Notations
The computational domain is denoted Ω t ; it is time-dependent and is the union of the fluid part Ω t f and solid part Ω t s .Technically Ω t must be an open set and so we denote by Ω t its closure.The fluid and solid domains must not overlap: Ω t f ∩ Ω t s = ∅.Initially, Ω 0 f and Ω 0 s are prescribed.Let: be the fluid-structure interface, and ∂Ω t be the boundary of Ω t , • Γ be the part of ∂Ω t where the structure is clamped or the fluid does not slip.It is assumed to be independent of t.
The following standard notations are used: is the Eulerian velocity of the deformation at t and x = X(x 0 , t), is the transposed gradient of the deformation, • J = det F is the Jacobian of the deformation.
Let ρ f , σ f be the density and stress tensor in the fluid, as for the solid with ρ s , σ s , respectively.It is convenient to define a unique density and stress tensor by using the set function indicators , the stress tensor.
For readability, vectors, tensors and matrices are noted in bold, except x and x 0 .Unless specified otherwise, all spatial derivatives are with respect to x ∈ Ω t and not with respect to When X is 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 a notation,

Conservation Laws
Conservation of momentum and conservation of mass take the same form for the fluid and the solid.With f (x, t), the volumic force in the system is: with continuity of u and of σ • n at the fluid-structure interface σ t in the absence of external surface force, like surface tension.There are also unwritten constraints pertaining to the realizability of the map X.Finally, incompressibility implies J = 1 and so ρ = ρ 0 along the Lagrangian trajectories.
In particular, if ρ 0 is piecewise constant and equal to ρ f , ρ s in the fluid and the solid, at initial time, then it remains so at later times.
• For a Newtonian incompressible fluid, where Ψ is an Helmholtz potential, like,

The Mooney-Rivlin 3D Stress Tensor
Note that ∂ F tr F T F = 2F and ∂ F tr (F T F) 2 = 4FF T F. Hence, In order to compute the tensors at Eulerian points (x, t) rather than Lagrangian points x 0 , t , it is necessary to compute x → F (x, t) or x → d (x, t) in terms of Ψ.By (2) By Cayley-Hamilton theorem in 3D, for an invertible symmetric matrix, and by incompressibility det B = 1.Multiplying (7) by B −2 , and the result by B −1 leads to polynomial expressions in terms of B −1 : With the help of ( 8), (5) becomes: Consequently, as a function of d, the 3D Mooney-Rivlin stress tensor is obtained from (6), with a complex expression for α in (10) which is unimportant for the time being because it is absorbed by the pressure.This is compared with the 2D model of [19] obtained with the same Ψ: with c 3 given by (10).

From 3D to 2D
In the 3D model, c 3 varies with tr B and tr B 2 instead of being constant, like c 1 in 2D model.The way to calculate c 3 is to compute directly tr B and tr B 2 by (6) which involves a 3 by 3 inverse matrix at each point of the solid domain.
The 3D model degenerates into a 2D system when the geometry and the data are invariant with respect to one coordinate, like a translation invariance with respect to z for instance, in a Cartesian coordinate system, or a rotation invariance in θ, in a cylindrical coordinate system.
Assume invariance with respect to the third coordinate.Then, Consequently tr Hence a comparison between the 2D model and the 3D model on a 2D configuration requires replacing the c 1 of the 2D Helmholtz potential by c 1 − 2c 2 from the 3D Helmholtz potential.

Variational Formulation
For simplicity, we shall consider only the case of homogeneous boundary conditions on Γ ⊂ ∂Ω t , i.e., clamped or no-slip, and homogeneous Neumann conditions on ∂Ω t \ Γ.
So, given Ω 0 f , Ω 0 s , d, and u at t = 0, we must find u, p, d, Ω t f , Ω t s with u |Γ = 0 and, for all ( û, p) with û|Γ = 0 and where Ω t s and Ω t f are defined by:

Characteristic-Galerkin Derivatives
The characteristic-Galerkin method is applied to (12) to discretize the total derivatives.Let u : Ω t × (0, T) → R 3 be given, Lipschitz in space and continuous in time.Let χ t u,x (τ) be the solution at time τ < t of: So χ is the particle path in a velocity field u and χ t u,x (τ) is the position of a particle at time τ which will be at position x at (future) time t.The method of characteristics relies on the concept of total derivative: for any differentiable function w : Given a time step δt, the position of a particle at time nδt, which will be at position x at time (n + 1)δt, is approximated by Y n+1 (x): Then a first-order consistent scheme for the total derivative of w is: Combined with a Galerkin approximation in space the scheme has been analyzed in depth in [21] in the variational context of the finite element methods.

A Monolithic Time-Discrete Variational Formulation
The following variational problem is a first order in time-consistent approximation of ( 12): where dn stands for d n (Y n+1 ).

Some comments
1. Formula ( 14) has been used to approximate D t u and D t d in (12).
2. One may wonder why the scheme is applied to u and not to ρu.
. This shows that discretizing the total derivative of u or the total derivative of ρu gives the same scheme: Proof.With the characteristic Galerkin method, a consistent time discretization of ( 12) would be such that at each time step one must: 0 Ω n+1 ; the three relations below hold: A fully implicit monolithic formulation with variables u n+1 , p n+1 can be derived by substituting d n+1 in ( 17) with (18): System ( 15) is found by expanding the nonlinear terms in (20), and keeping only the zero order terms and the terms of order one in δt.

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

Solution Algorithm
Equation ( 21) must be solved iteratively because Ω n+1 is updated by (19).The most natural method is to freeze some coefficients, c 3 included, so as to obtain a well-posed linear problem and iterate.To explain the iteration process let us rename, in (19) (21).In this study a direct solver is used for the linear system.
Remark 1. Ω n+1 s we move the vertices q j in the solid part by its velocity u q j .As dn is moved by Y, it is naturally convected on the new mesh.Similarly, ρ being constant on Ω s , it is well defined on the new mesh.Finally, c 3 being defined in terms of u, it is also well defined on the new mesh.
In the fluid region we generate a new mesh with a Delauney-Voronoi automatic mesh generator from the new position of the interface σ.There are many other possibilities like solving a Laplace equation in the fluid part with Dirichlet conditions equal to the velocity of the fluid-solid interface.It works equally well but what we advocate is, in our experience, more robust, even if it is at the cost of a possibly greater interpolation error at every step because ũn needs to be projected on the new mesh.A P 2 polynomial interpolation is used.
Remark 2. An earlier study in [17] showed that 2 iterations are necessary but 3 or more make invisible improvements; this result is confirmed in 3D by the numerical tests below.
By inspection of ( 23) below, notice that this iterative algorithm computes the stationary point As Ψ is convex, it should be a convergent process, but it requires a more rigorous proof.

Conservation of Energy
Proposition 2. In the continuous case, it holds that: Proof.Choosing û = u, p = −p in Equation ( 12), then one term becomes (9) : ), so conservation of the Helmholtz potential Ψ can be derived as follows: The other terms in (12) are standard, in particular, with enough regularity on Ω t , Remark 3. When Ψ is convex, some regularity can be gained from this energy conservation which can lead to existence of solution up to time T which, loosely speaking, in the most optimistic perspective, is the first time of contact of two boundaries which were not in contact initially (see [22][23][24] etc.).

Stability of the Scheme Discretized in Time
Lemma 1.
Theorem 1.When ( 17) is used and if f = 0 and ρ r is constant in each domain Ω n r , r = s, f , Proof.The proof is similar to that given in [19] for the 2D case.With û = u n+1 , p = −p n+1 , and f = 0, (17) becomes: We note that ; so: Schwartz inequality applied to the right side followed by Young's inequality leads to: Because (Y n+1 ) −1 (Ω n+1 ) = Ω n , and det ∇Y n+1 = 1, for any f , The third integral in ( 25) can be transformed with Equation ( 9): The last equality is due to Lemma 1. Now, by the convexity of Ψ:

Energy Inequality for the Fully Discrete Scheme
It was argued in [18] that the same proof holds for scheme discretized with P 1 triangular elements for velocity and pressure because then Y transforms a tetrahedron into a tetrahedron and so (27) holds also in the discrete case.Compatible linear elements are either the P 1 − P 1 stabilized element or the P 1 element for the pressure with P 1 also for the velocity but on a mesh where each tetrahedron is divided in sub-tetrahedra from an inner additional vertex.In the first case the above stability analysis needs to be modified to include the stabilization term.Numerical tests in 2D showed that the P 2 − P 1 element is alright even though energy stability cannot be proved.

Numerical Tests
The first set of tests is meant to evaluate the precision of the method.Due to limited computer resources (the program is not parallelized), we chose an axisymmetric example so as to work in 2D with four embedded meshes and study convergence, conservation of mass and energy decay.
The second set concerns a well-documented test case in 3D with experimental comparisons.
Computing time on a Mac-Intel core 2 i7 at 3 GHz ranges from 7sec for mesh 1 to 2800sec for mesh 7 on the axisymmetric case with two iterations in the solver.For the benchmark 3D case it takes several hours.The computer program is written in the FreeFem++ language [25] (see also www.freefem.org),which is an extension of C++.

An Axisymmetric Torus
A torus of rubber-like material is immersed in a fluid filling a cylindrical canister.Assuming axial symmetry, the computational domain is a rectangle (0, L) × (0, H) with a disk inside of radius R and center x c , y c .
The disk is the cross-section of the torus which is an incompressible Mooney-Rivlin material characterized by c1 = c 1 − 2c 2 and ρ s , as explained in Section 2.5.The rectangle is filled with a Newtonian fluid characterized by its density ρ f and viscosity µ f .
At time zero a centripetal horizontal uniform velocity U 0 is imposed on the torus.Consequently, the torus contracts, and its cross-section, changes to preserve mass, until a position nearest to the axis of symmetry is reached.Then, it bounces away from its axis to a maximal position and repeats these oscillations a few times until fluid viscosity dampens the motion.
The tests are done with T = 3, corresponding to 2.3 periods; the parameters are: The problem is solved on four triangulations, called mesh i , i = 1, 3, 5, 7; the associated time step size is δt = T/(50 × i).The number of triangles in the meshes are 775, 6522, 17,749 and 35,308.
Figure 1 shows the initial state and some snapshots at times of maximum vertical and horizontal elongations, and farthest travel to the right.Figure 1.Deformation of a rubber-like torus in a fluid due to initial centripetal velocity of the torus; snapshots at t = 0 (a), 0.12 (b), 0.28 (c), 1.45 (d) show also the four meshes, from finest to coarsest.The torus moves towards its axis (left boundary) and then bounces to the right to a maximum position and return to its initial position after a few oscillations.Notice the change of area to preserve mass.
Decrease of energy (viscous terms included) and conservation of mass are shown on Figure 2a.Energy is computed using mesh 1 and 1, 2, 3 or 4 iterations for the non-linear solver of Section 3.4; the results show that 2 iterations for the nonlinear solver is optimal.Energy is also plotted when mesh 3 and mesh 5 are used (with 2 iterations).In agreement with the theory, energy decays with time even on a coarse mesh.The decay is due to the numerical dissipation because for the continuous problem it should be constant, and indeed it is less with finer meshes.Note that an impulsive start is not an energy-friendly test case.
Conservation of mass is shown on Figure 2b for mesh 1 , mesh 3 and mesh 5 .It is within 2%.On the same plot energy vs. time is also plotted for mesh 3 and mesh 5 .(b): Mass conservation vs. time when mesh 1 , (lower curve), mesh 3 and mesh 5 are used.Curves 3 and 5 are overlapping.
To study convergence we take the results obtained with mesh 7 as a reference "almost exact" solution.Then, we observe the point on the torus with maximum radius (i.e., farthest to the right).On Figure 3a the maximum radius is shown versus time for the four meshes; on Figure 3b the pointwise error for mesh i , i = 1, 3, 5. Convergence is observed; we make no claim on the speed of convergence versus mesh size because mesh 7 is probably too coarse to be reliable.The errors on the farthest point on the right in the torus versus time are shown for mesh 1 (highest curve, its mean absolute value is 0.013), mesh 3 (lowest curve, its mean absolute value is 0.0023), and mesh 5 (flat curve, its mean absolute value is 4.1 × 10 −5 ).

Clamped Structure in a Fluid
A set of computations of a beam in vacuum or in fluid has been performed at the request of Larma's PhD advisor, Miguel Fernandez [26].Experimental measurements are available [20] so it is intended to make them into a well-documented benchmark for the community.The geometry is always a clamped beam in a container; in all but the first case there is a flow in the container but to initialize the computation the position of the beam is computed with a stagnant fluid.
The beam in Hessenthaler's experiment and Larma's simulation is not incompressible, as the Poisson ratio is 0.35 instead of 0.5, so we expect discrepencies.The present method can generalized to compressible material and 2D simulations show similar stability properties but the computer program is not yet ready for 3D cases.At the end of computation, the error of volume is less than 1%.A discussion about c 3 is found later in Section 5.3.3.

Free-Fall of a Clamped Beam in a Fluid
Second, we consider a flat beam of size 9 × 1 × 4 with one side attached to a box which is filled with a fluid.The box size is 10 × 7 × 5. Solid density ρ s is 1.Fluid density ρ f is 0.5, and viscosity µ f equals 0.1.The gravity g is −0.

The Benchmark
An experiment was designed by Hessenthaler et al. [20] to calibrate the computer simulations of Larma [26].
The geometry is a cylindrical chamber with a length of 200 mm, diameter 76.2 mm and axis parallel to x = y = 0. Two inlet circular wholes are on the left z = 0 wall, one with center at (0 , −27.15, 0) , the other at (0, 27.15, 0), both of diameter 21.9 mm; the inlet pipes are 20 mm long.A silicon filament clamped on the z = 0 plane is 2 mm thick, 11 mm wide, and 65 mm long; its center is located at (0, 0, 0).The computational domain is discretized into 10 4 vertices and adapted to the silicone filament to decrease error during its motion.Geometric schematic is shown in Figure 6.
As in [20], the density of the silicone filament is ρ s = 1.0583 × 10 −3 g/mm 3 .The incompressible Mooney-Rivlin hyperelastic coefficients are determined by curve-fitting to uni-axial tensile-load displacement test data.In [26] two coefficients are introduced: ĉ1 = 103, 533.82/2Pa and ĉ2 = 8891.65/4Pa.The Poisson ratio is ν = 0.315 in [26] and c 1 = ĉ1 2(1+ν) , c 2 = ĉ2 2(1+ν) .There will be a difference because incompressibility implies ν = 0.5.Gravity g is −9810 mm/s 2 along the y-direction.Two numerical tests are made, a steady phase I, and a transient phase II.In absence of flow, due to buoyancy forces, the silicone filament reaches an hydrostatic equilibrium with maximum displacement d 0 M , at the tip on the right, approximately equal to 29.50 mm in [20] and 31.04 mm in [26].To compute this displacement we have taken T = 1 and 1000 steps with δt = 0.001.
Then d 0 is taken as initial condition to Phase I as the inflows are turned on; the maximum displacements are now approximately 16.41 mm in [20] and 15.99 mm in [26] .To compute this state we have taken T = 2, 2000 steps with δt = 0.001.Our results, shown in Figure 7, overlap the experimental data.The Poiseuille profile flow applied on the upper pipe has a transient peak velocity vk (t), k ∈ [x, y, z] shown in Figure 8 and corresponds to a smooth interpolation from Table 1.No flow is applied to the lower inflow pipe.
Computation stops at T = 5 after 1600 steps with δt = 0.003125.The maximum velocity norm of the solid is less than 30 mm/s; the error introduced by the computation is less than 0.1 mm.Deflection of the center line of the silicone filament is recorded versus time and a comparison between numerical results and experimental data is shown in Figure 9 at six different times.These results display a fair-yet not satisfactory-agreement between the simulations and the experimental data.The difference is perhaps due to improper determination of the coefficients of the hyperelastic model, because bending at the end of silicone filament could be observed only when the magnitude of c 2 is large enough.Another possible explanation is the wrong Poisson ratio, namely an effect due to incompressibility.Table 1.For phase II, curve-fitting coefficients of inlet peak velocity for vk (t) = σ 3 i=1 n i t i /σ 4 j=0 b j t j with vk = 0 for t ∈ I\I k .Note that the flow in y direction is applied only in the upper inlet.It is interesting to point out that c 3 is not very sensitive to the displacement gradients in this simulation as shown in Figure 11: c 3 remains fairly constant during the computation, even though d does not.These results seem reasonable since deflection over the whole silicone filament is fairly uniform.However, we do not rule out that strong deformations could cause a stability problem in the algorithm; this will be investigated in the future.

Conclusions
In this study, a fully Eulerian fluid-structure formulation has been presented.An implicit unconditionally stable monolithic finite element scheme based on characteristic-Galerkin discretization has been proposed and studied.At each time step a non-linear system must be solved; however semi-linearization leads to well-posed linear monolithic systems for the velocities and pressures.
What was found in [17,18] for 2D systems could be generalized in 3D with similar stability based on analytical proof that the energy of the discrete system cannot increase.
The numerical method has been implemented with FreeFem++ [25], and validation tests have been presented.For an axisymmetric case, conservation of mass and energy were shown as well as convergence on a sequence of embedded meshes.
The effect of nonlinearity of c 3 -the major difference with the 2D case-is studied in Section 5.3.3 with validation tests.
In the near future we intend to parallelize the computer program and use it for larger systems, in particular for hemodynamics.

Figure 2 .
Figure 2. (a): Energy vs. time computed on mesh 1 with 1, 2, 3 and 4 iterations in the solver of Section 3.4.On the same plot energy vs. time is also plotted for mesh 3 and mesh 5 .(b): Mass conservation vs. time when mesh 1 , (lower curve), mesh 3 and mesh 5 are used.Curves 3 and 5 are overlapping.

Figure 3 .
Figure 3. (a): Evolution of the farthest point on the right of the torus versus time computed with mesh i , i = 1 (oscillations), 3 and 5. (b): Mesh 7 is used as a reference computations.The errors on the farthest point on the right in the torus versus time are shown for mesh 1 (highest curve, its mean absolute value is 0.013), mesh 3 (lowest curve, its mean absolute value is 0.0023), and mesh 5 (flat curve, its mean absolute value is 4.1 × 10 −5 ).

5. 2 . 1 .
Free-Fall of a Clamped Structure in Vacuum The beam size is 9 × 1 × 1 .The beam is clamped on the right to the yz plane and free to fall under its own weight.The density ρ s is 1, and the gravity is −0.01 along y direction.The Mooney-Rivlin coefficients are c 1 = 5 and c 2 = 2.5.Computation stops at T = 100 after 3200 steps with δt = 0.03125.Results are shown at different times in Figure 4.

5
along the y direction.The Mooney-Rivlin coefficients are c 1 = 16.67 and c 3 = 8.33 tr 2 B − tr B 2 − 4 .Computation stops at T = 100 after 200 steps with δt = 0.5.Computing results are shown in Figure 5; volume error is less than 0.5% at the end.

Figure 7 .
Figure 7. Computational results for phase I (a) position of center line of the beam along z direction (results and experiments overlap), (b) velocity norms about the symmetry plane.

Figure 8 .
Figure 8. Experimental curves of the peak inflow of inflow velocities used for the boundary condition in the computations.

Figure 10 .
Figure 10.Phase II results (a) maximum and minimum magnitude of 2c 3 and relative volume error multiplied by 1000, (b) maximum magnitude of displacement d.