Christophe Kassiotis
email: christophe.kassiotis@enpc.fr

Adnan Ibrahimbegovic

Hermann G Matthies

July

Partitioned solution to fluid-structure interaction problems in application to free-surface flows

In this work we discuss a way to compute the impact of free-surface flow on nonlinear structures. The approach chosen rely on a partitioned strategy that allows to solve strongly coupled fluid-structure interaction problem. It is then possible to re-use existing and validated strategy for each sub-problem. The structure is formulated in a Lagrangian way and solved by the finite element method. The free-surface flow approach considers a Volume-Of-Fluid (VOF) strategy formulated in an Arbitrary Lagrangian-Eulerian (ALE) framework, and the finite volume are used to discrete and solve this problem. The software coupling is ensured in an efficient way using the Communication Template Library (CTL). Numerical examples presented herein concern 2D validations case but also 3D problems with a large number of equations to be solved.

Introduction

In this work we focus on fluid-structure interaction problems, with free-surface flows. The kind of application motivating this development stems from the coastal engineering problems, chief among them dam-breaking problem and tsunami (e.g. see [START_REF] Fochesato | Numerical modeling of extreme rogue waves generated by directional energy focusing[END_REF]). The main difficulty for this class of problems is not only the usual one for fluids pertaining to the free-surface flow (e.g. see [START_REF] Dutykh | On the relevance of the dam break problem in the context of nonlinear shallow water equations[END_REF][START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF][START_REF] Ghidaglia | On the numerical solution to two-fluid models via a cell centered finite volume method[END_REF][START_REF] Behzadi | Modelling of dispersed bubble and droplet flow at high phase fractions[END_REF][START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF][START_REF] Sobey | Linear and Nonlinear Wave Theory[END_REF][START_REF] Stoker | Water Waves: The Mathematical Theory and Applications[END_REF]) but also the fluid-structure interface. Providing reliable description at the interface for a fluid-structure interaction problem is quite a challenge, since the descriptions generally used for each of the sub-problems are different: for the structure part, it is natural to follow material point motion in a Lagrangian formulation, while an Eulerian formulation is often preferred for the fluid part. Naturally, there are bridges between these two formulations, which have to be derived carefully for fully nonlinear setting.

Moreover, the optimal discretization technics for each sub-problem are in general not the same. For the solid part, the method widely used since the sixties is the Finite Element Method [START_REF] Zienkiewicz | The Finite Element Method[END_REF][START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF], even if it is possible to solve the Partial Differential Equations that describe the continuum problem either by other discretization techniques such as Finite Volume [START_REF] Slone | Dynamic solid mechanics using finite volume methods[END_REF]. For the fluid part, two different methods are mainly used: either the Finite Element Method [START_REF] Zienkiewicz | The Finite Element Method[END_REF][START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF], or the Finite Volume Method [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. The first one has the advantages of a complete mathematical framework and well established convergence properties, but requires special care to stabilize the solution when incompressible flow are aimed at [START_REF] Franca | Stabilized finite element methods[END_REF]. For the second one, the mathematical framework is slightly more difficult to complete but the more physical foundation allows to naturally get good conservative properties. There are bridges that combine advantages of both formulations [START_REF] Pascal | Footbridges between finite volumes and finite elements with applications to CFD[END_REF]. In this work, we choose the Finite Element Method for structure discretization and the Finite Volume Method for fluid discretization. This choice provides the largest selection of existing software products to be used for the software coupling strategy we suggest herein.

For the fluid-structure interaction framework, it is required to consider flow problem in a moving domain, as the fluid domain shares at least one of its boundaries with the structure that is undergoing deformation. It is possible to solve this kind of problem with for instance the fictitious domain approach [START_REF] Wang | The immersed-fictitious element method for fluid-structure interaction: Volumetric consistency, compressibility and thin members[END_REF] or with an Arbitrary Lagrangian-Eulerian approach [START_REF] Demirdžić | Space conservation law in finite volume calculations of fluid flow[END_REF][START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF]. The latter requires to solve the deformations of the fluid underlying grid, that is often not able to support very large deformation without re-meshing; nevertheless it is often the one chosen in existing tools to solve fluid problems in moving domain, and was thus also used herein. The impossibility to represent too large deformations without re-meshing makes the choice of another method necessary when complex fluid domain deformations such as the ones observed for the sloshing of a wave are represented. Here a two-phase approach is chosen where both air and water are represented, and a characteristic function is used to distinguish the two phases. To solve this problem, a Volume-Of-Fluid strategy is considered herein.

For coupled fluid-structure interaction problems, the monolithic [START_REF] Hübner | A monolithic approach to fluidstructure interaction using space-time finite elements[END_REF][START_REF] Perić | Modelling fluidinduced structural vibrations: reducing the structural risk for stormywinds[END_REF][START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF]] and partitioned strategy [START_REF] Farhat | Two efficient staggered algorithms for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems[END_REF][START_REF] Tallec | Fluid structure interaction with large structural displacements[END_REF][START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Fernández | A Newton method using exact Jacobians for solving fluid-structure coupling[END_REF][START_REF] Matthies | Algorithms for strong coupling procedures[END_REF][START_REF] Deparis | Fluidstructure algorithms based on Steklov-Poincaré operators[END_REF][START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF][START_REF] Fernández | Domain decomposition based Newton methods for fluid-structure interaction problems[END_REF][START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF]] can be used. The monolithic approach is abandoned in favor of the partitioned approach. The latter is preferred for its modularity and the possibility of re-using existing software. The partitioned approach used here is based on the Direct Force-Motion Transfer. Both explicit and implicit coupling algorithms for multi-physics problems are detailed. An explicit strategy leads to the so-called "Added Mass effect", and for that justifies the use of more costly implicit solvers.

Implicit solvers can use either a fixed-point strategy that is known to be rather slow to converge, or the Newton strategy that requires building up and evaluating the costly Jacobian. In this work, the fixed-point strategy based on the Block Gauss-Seidel algorithm (DFMT-BGS) with an adaptive relaxation parameter [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF] shows sufficient performances for the example proposed. The properties as well as stability of the implicit coupling DFMT-BGS algorithms used herein are presented in detail in [START_REF] Kassiotis | Partitioned procedure for strongly coupled fluid-structure problems. Part I: implicit algorithm, stability analysis and validation examples[END_REF].

In this work, a general fluid-structure interaction framework based on existing software was used. This framework was built using the middleware Communication Template Library (CTL) [START_REF] Niekamp | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior: Part II-software coupling implementation aspects[END_REF] which offers good performances, and can therefore be used for scientific computing of large systems. An important feature is the possibility to couple types of software product that were initially not programmed to be coupled (here FEAP for the structure and OpenFOAM for the fluid), even if they are based on different discretization techniques (respectively FV for the fluid and FE for the structure) and were programmed in different languages, C++ and Fortran. For more details on the implementation, see [START_REF] Kassiotis | Partitioned procedure for strongly coupled fluid-structure problems. Part II: implementation aspects, re-using existing software as components and nested parallel computations by using ctl[END_REF] The outline of this paper is as follow: in the subsequent section we present the chosen formulations for the structure and fluid sub-problem. In particular, the fluid formulation is given in Arbitrary-Lagrangian Eulerian framework, along with the modifications to account for the free-surface flow. In Section 4, we describe the coupling between the fluid and the structure sub-problems. In Section 5, we give and comment the results of illustrative numerical examples dealing with free-surface flow impacting a structure in two and three-dimensions as well as comparison with existing works. The concluding remarks are given in the last section.

2 Description of the structure 2.1 Structure equations of motion

Ω s t = t 0 ∂Ω s,D ∂Ω s,N Ω s ∂Ω s,D ∂Ω s,N t t b u t b u Figure 1: Solid problem
The structure motion is based on a Lagrangian description. Namely, we consider a structure domain Ω s , with imposed displacements u on its Dirichlet boundary ∂Ω s,D , moving under the traction forces t imposed on its Neumann boundary ∂Ω s,N and a body load b applied in the whole domain Ω s (see Fig. 1). The motion is dynamic and evolves in time on the segment [0, T ]. The governing equation for a structure describes the momentum conservation and is also known as the Cauchy equation. The strong form of the structural problem written with respect to the deformed configuration is as follows:

Given: u on ∂Ω s,D × [0, T ], t on ∂Ω s,N × [0, T ] and b in Ω s × [0, T ]. Find: u ∈ Ω s × [0, T ] so that: ∇ • σ + ρ s b -∂ 2 t u = 0 in Ω s × [0, T ] (1)
where ρ s denotes the material density of the solid domain, u its displacement field and ∂ 2 t u the accelerations. The Cauchy stress tensor σ in the deformed configuration can be linked to the second Piola-Kirchhoff stress tensor P formulated in the initial configuration through the gradient F of the deformation and its Jacobian (see [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF]).

To close this Partial Differential Equations system one needs to link the displacements (or one of its derived field) and the stresses together with constitutive behavior law. An elastic material model based on St.-Venant-Kirchhoff constitutive equation or neo-Hookean constitutive equation is assumed in the two presented example to link the stress Cauchy tensor σ and the Green-Lagrange strain tensor E:

F -1 JσF -T = C : E (2)
where C denotes the constitutive fourth-order elasticity tensor1 . It is also possible to use more complex behavior, such as non-linear plasticity, but this is not explored herein. The non-linearity of the problem here comes from the large displacement described by the following relation between Green-Lagrange tensor E and the deformation gradient F = I + ∇u:

E = 1 2 F T F -I (3)

Structure discretization using Finite Element Method

It is a priori impossible to find directly the exact solution to the problem defined above. The idea is to find the best approximation of the solution in a space where the solution can be found numerically. The FEM formulations [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF][START_REF] Zienkiewicz | The Finite Element Method, The Basis[END_REF][START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF] derived from the equilibrium equation ( 1) rely on the associated weak forms of this problem:

Given: t on ∂Ω s,N × [0, T ] and b in Ω s × [0, T ]. Find u ∈ U such that, for all δu ∈ U 0 : G s (u; δu) := Ωs ρ s ∂ 2 t u • δu + Ωs σ : ∇δu - Ωs b • δu - ∂Ωs t • δu = 0
where U and U 0 are functional spaces for the solution and its variation.

The solid domain Ω s is then discretized in a finite number of elements T h = (κ e ) e=1,...,n el so that the whole space is covered by the finite elements that do not intersect. The research space associated with the solution is restrained to the space of continuous piecewise polynomial functions defined in each finite element domain. This space is denoted:

U h = U ∩ u ∈ C 0 (Ω s ) u| κ ∈ P p (κ) , ∀κ ∈ T h (4)
where P p (κ) is the space of polynomials of order p on κ. The same restriction holds on the associated vector space. The FE problem is defined as:

Given: t on ∂Ω s,N × [0, T ] and b in Ω s × [0, T ]. Find u ∈ U h such that, for all δu ∈ U h 0 : G s (u; δu) = 0 (5)
This semi-discrete problem can be written in a matrix form using the real valued

vectors u ∈ | n d-o-f : R s (u s ; λ) := M s üs + f int s (u s ) -f ext s (λ) = 0 (6)
with M the mass matrix, K the stiffness matrix associated with a potentially nonlinear problem, and f the projected loading forces. Here the λ represents the boundary forces computed from the fluid flow problem and imposed on the fluidstructure interaction interface, and therefore a restriction of t. Each matrix and vector of this semi-discrete equation are properly defined by assembling locally computed array on each element with the polynomial basis N e of P (κ e ):

M s,e = κe ρ s N T e N e
The time integration of the structure problem can be carried out by method such as standard time stepping schemes [START_REF] Belytschko | An overview of semidiscretization and time integration procedures[END_REF][START_REF] Hughes | Implicit-explicit finite elements in nonlinear transient analysis[END_REF] or by the so called Finite Element in time as introduced in [START_REF] Argyris | Finite elements in time and space[END_REF]. In [START_REF] Dettmer | An analysis of the time integration algorithms for the Finite Element solution of incompressible Navier-Stokes equations based on a stabilised formulation[END_REF], the use of Finite Element in time for the structure part in a fluid-structure interaction context is shown to increase the overall computational with no noticeable advantages. The Generalized HHT-α method [START_REF] Chung | A family of single-step houbolt time integration algorithms for structural dynamics[END_REF] is used herein. The time interval [0, T ] is discretized into a finite number of time-steps t N such as t 0 = 0 and t Nmax = T . In a typical time step size ∆t = t N +1t N , and time derivatives are approximated with:

u N +1 = u N + ∆t uN + ∆t 2 1 2 -β üN + βü N +1 uN+1 = uN + ∆t [(1 -γ)ü N + γü N +1 ] u N +α f = (1 -α f )u N + α f u N +1 üN+αm = (1 -α m )ü N + α m üN+1 (8)
In the semi-discrete form of the solid equation of motion in Eq. ( 6), the acceleration ü and the displacement u are evaluated at t N +α f and t N +αm For the elastic linear case, it is shown [START_REF] Chung | A family of single-step houbolt time integration algorithms for structural dynamics[END_REF] that there are optimum values for the parameters β, γ, α and α for a given spectral radius ρ ∞ ∈ [0, 1].

β = (1 + α m -α f ) 2 4 , γ = 1 2 + α m -α f , α f = 1 1 + ρ ∞ and α m = 2 -ρ ∞ 1 + ρ ∞ (9)
The spectral radius controls the numerical damping of the time integration scheme. The damping decreases with smaller values of ρ ∞ which is maximum for ρ ∞ = 0. For ρ ∞ = 1 the method is the classic trapezoidal rule. Other time integration schemes can be easily derived from this general formulation [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF].

3 Description of the flow with free-surface

Free-surface flow equations

To model free surface flows, a traditional approach is to consider simplified model such as the Nonlinear Shallow Water Equations first introduced by Saint-Venant in 1837. Extensive literature on how to solve those equations analytically or numerically can be found in a number of previous works (e.g. see [START_REF] Stoker | Water Waves: The Mathematical Theory and Applications[END_REF][START_REF] Sobey | Linear and Nonlinear Wave Theory[END_REF]). If this approach gives good results for the propagation phase of the wave, in [START_REF] Dutykh | On the relevance of the dam break problem in the context of nonlinear shallow water equations[END_REF] these models are shown to be unable to represent accurately the impact process of the wave, as well as its sloshing process (see Fig. 2(a)). For this reason, another appropriate model has to be chosen near a structure hit by water. To represent the complex flow an idea could be to consider a moving fluid domain where the shape of the wave is described by the motion of the domain. However, the domain motion is so complicated that this strategy is not possible to use the traditional ALE formulation given above.

Several alternative approaches are thus proposed:

Meshed Lagrangian approach: (see Fig. 2(a)) are used to represent the free surface fluid flows by using constant re-meshing; therefore, efficient mesh generators are needed in order to compute the new topology at each iteration. For instance, the PFEM method [START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF] uses one such fast meshing algorithm and solve this kind of problems.

Meshless Lagrangian approach: (see Fig. 2(b)) such as the SPH has the advantage not to require the re-meshing step of the previous strategies [START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF][START_REF] Dalrymple | Numerical modeling of water waves with the SPH method[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF][START_REF] Violeau | Numerical modelling of complex turbulent freesurface flows with the SPH method: an overview[END_REF].

Surface tracking methods (ALE): (see Fig. 2(d)) has the advantage that the representation of the fluid domain and its free surface is independent from the representation of the flow field. Hence, the resolution of the surface and that of the flow field may be chosen independently from one another. Of course, the level of details should be comparable in order to resolve the fluid motion properly. However, this freedom is helpful, for example, in order to improve the accuracy of the evaluation of the surface tension force (e.g. see [START_REF] Ma | QALE-FEM for numerical modelling of non-linear interaction between 3D moored floating bodies and steep waves[END_REF]).

The relative efficiency of this method will not be discussed herein. As the result of the initial choice made to couple a FV code, non of them was used. We use an Arbitrary Lagrangian-Eulerian description of a two-phase flow (taking into account both water and surrounding air) discretized by Finite Volume (see [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF][START_REF] Ghidaglia | On the numerical solution to two-fluid models via a cell centered finite volume method[END_REF][START_REF] Behzadi | Modelling of dispersed bubble and droplet flow at high phase fractions[END_REF]). In this Volume-Of-Fluid (V.O.F.) method, an indicator function (volume fraction, level set or phase-field) is used to represent the interface (see Fig. 2(e)); the main remaining issue is how to convect the interface without diffusing, dispersing or wrinkling it [START_REF] Dias | A two-fluid model for violent aerated flows[END_REF][START_REF] Ghidaglia | On the numerical solution to two-fluid models via a cell centered finite volume method[END_REF][START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF][START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF].

Remark: The incompressible hypothesis can be discussed, especially for highly aerated waves. For more details on compressible multiphasic simulations, the reader is invited to consults works of [START_REF] Dias | A two-fluid model for violent aerated flows[END_REF].

For complex flows (with jets, cavitation and aeration in the sloshing wave) it is natural to consider the Navier-Stokes equations for two immiscible and incompressible flows (water and air for instance) occupying transient domains Ω i (t) so that the whole fluid domain considered

Ω f (t) = Ω 1 (t) ∪ Ω 2 (t).
The interface between the both domains Ω 1 and Ω 2 is denoted Γ . In space-time domain Ω f (t) × [0, T ], the Navier-Stokes equations formulated in an ALE framework apply. For an arbitrary motion of the total fluid domain Ω f described by a displacement field u m , it can be written as:

ρ∂ t v + ρ (v -um ) ∇ • v -∇ • 2µD(v ) = -∇p + f Γ + ρg in Ω f (t) × [0, T ] ∇ • v = 0 in Ω f (t) × [0, T ] (10
) where p denotes the pressure field, v the velocity. We also introduce g that depicts the gravity field, f Γ the surface tension forces. It can be expressed as f Γ = σκδ Γ , where σ is the surface tension, κ the curvature of the free-surface (i.e.interface between Ω 1 and Ω 2 ) and δ Γ the mass distribution concentrated at the surface (equivalent to a Dirac distribution). Fluid material properties are the dynamic viscosity µ and the density ρ. To write a unique formulation in the whole domain Ω f (t) we express the local material property values as the function of ι:

ρ = ιρ 1 + (1 -ι)ρ 2 and µ = ιµ 1 + (1 -ι)µ 2 in Ω f (t) × [0, T ] (11)
where the characteristic function or fluid volume fraction ι is defined as:

ι(x , t) = 1, for x ∈ Ω 1 (t) 0, for x ∈ Ω 2 (t) (12)
Thus, the fluid volume fraction ι and the mass distribution δ Γ are linked by the relation:

∇ι = δ Γ n (13)
To close the set of equations we have to write the conservation of ι. When no reaction between phases occurs, the fluid volume fraction evolves only by advection:

∂ t ι + (v -um ) • ∇ι = 0 (14)
The conservation equation system on the whole domain Ω f in d dimensions can be written as a function of the 2d + 2 unknowns: u m , v , p and ι which are sought in their respective spaces U , V, P and I that denote the restriction of Hilbert and Sobolev spaces with suitable boundary condition. Thus, we aim to find: (u m , v , p, ι) ∈ U × V × P × I, such that equations ( 10) and ( 11) are verified.

Arbitrary moving domain description

A difficulty in fluid-Structure interaction problems pertains to unsteadily moving domains for the fluid part, as the fluid-structure interface follows the deformations of the structure. Traditional Computational fluid Dynamic programs solve the fluid equations on a fixed (Eulerian) grid. The classical approach to overcome this difficulty is to consider the so-called Arbitrary Lagrangian Eulerian (ALE) method where the whole grid is moved inside the fluid domain following the movement of the boundary [START_REF] Demirdžić | Space conservation law in finite volume calculations of fluid flow[END_REF][START_REF] Förster | On the geometric conservation law in transient ow calculations on deforming domains[END_REF][START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF][START_REF] Hughes | Lagrangian-Eulerian finite element formulation from incompressible viscous flows[END_REF][START_REF] Lesoinne | Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations[END_REF][START_REF] Wang | The immersed-fictitious element method for fluid-structure interaction: Volumetric consistency, compressibility and thin members[END_REF][START_REF] Ma | QALE-FEM for numerical modelling of non-linear interaction between 3D moored floating bodies and steep waves[END_REF]. However, this leads to a new difficulty: knowing the new shape of the boundary, how can one keep the quality and the validity of the fluid inner mesh? In more mathematical terms, this means building a suitable map for the domain motion given its interface displacement. The fluid displacement u m is arbitrary inside the domain Ω f , but has to fulfill the condition:

u m = u m on ∂Ω f (t) × [0, T ] (15)
Thus, the fluid displacement is an arbitrary extension of u m | ∂Ω f inside the fluid domain Ω f :

u m = Ext u m ∂Ω f (16)
Let us consider a domain T h (Ω f ,t N ) which represent the mesh configuration at a time t N . This domain is moving according to the imposed condition at the boundaries, in terms of either imposed displacement or velocity. The issue is to build a new valid mesh T h Ω f ,t N +1 at t N +1 = t N + ∆t knowing the initial valid mesh T h (Ω f ,t N ) and imposed boundary conditions.

It can be proved that the mesh motion problem is formally equivalent to a solid body under large deformations. The cost to solve the resulting set of non-linear equations is tremendous and this strategy cannot be applied within an efficient 3D CFD code. A traditional way to overcome this difficulty is to consider simplified solid equations for the mesh motion:

Spring analogy aims at linking each point of the mesh by fictitious spring.

In [START_REF] Jasak | Automatic mesh motion for the unstructured finite volume method[END_REF], authors show this yields failure modes. The failure modes of this kind can be eliminated by introducing non-linear and torsional springs. However, the cost induced by the improvement of this imperfect system by nature can be deemed as too large.

Pseudo-solid equations can be considered as a simplification of the 3D nonlinear problem to the linear one; the analogy is a mechanical problem with small deformations.

Laplacian smoothing operator is different strategy where we consider the Laplace smoothing equation:

∇ • (γ∇u m ) = 0 in Ω f (17)
A study of some mesh motion algorithms for the examples presented herein is given in [START_REF] Kassiotis | Which strategy to move the mesh in the Computational Fluid Dynamic code OpenFOAM?[END_REF] 2 . The main conclusion of this work can be summarized as:

Laplace smoothing operator with a non constant diffusion coefficient can handle large deformation of the mesh while the rotations remains small.

Pseudo-solid formulation is of interest when the motion of the mesh is undergoing large deformation and is mainly governed by rotation.

When the mesh motion are too large (this is not the case herein), strategies such as Cartesian Cut Cells [START_REF] Causon | Calculation of shallow water flows using a Cartesian cut cell approach[END_REF], Immersed Fictitious Elements [START_REF] Wang | The immersed-fictitious element method for fluid-structure interaction: Volumetric consistency, compressibility and thin members[END_REF], or even meshless methods [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] can be used.

Finite Volume discretization

The continuous equations are transformed into a set of algebraic equations that can be solved numerically. Indeed, the goal of this work is first to emphasize the possibility of coupling different codes, with different methods (FEM for solid, FVM for fluid) and second to follow the trends of the computational scientific software market. And, as FVM, contrary to FEM, leads to intrinsically conservative methods at the local stage they are often preferred in commercial and non-commercial softwares (Phoenics, fluent, flow3D, Star-CD, Code Saturne, OpenFOAM [START_REF] Jasak | OpenFOAM: Open source CFD in research and industry[END_REF] are FVM based, Adina is FEM based). There are a great number of books of reference on the topic, among them [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF][START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF].

The FV formulation can be written directly using an integrated form of the conservation equation written in [START_REF] Bullock | Violent breaking wave impacts. Part I: Results from large-scale regular wave tests on vertical and sloping walls[END_REF]. Another possibility is to consider the restriction of the solution space of weak form problems. For homogeneity with the solid method we chose to describe the FV strategy in this framework. The weak form of the Navier-Stokes equation can be written as follows [START_REF] Glowinski | Numerical methods for fluids (Part III)[END_REF]. Find (u m , ι, v , p) ∈ U ×I×V×P, such that, for all (δu m , δι, δv , δp) ∈ U 0 ×I 0 ×V 0 ×P 0 :

G f := Ω f ∇ • (γ∇u m )δu m + Ω f (∂ t ι + (v -um ) ∇ι) δι Ω f ρ∂ t v • δv + Ω f ρ∇ (v -um ) ⊗ v • δv - Ω f ∇ • µ f D(v ) • δv + Ω f p∇ • δv + Ω f ∇ • v δp + [B.C. in a weak form] = 0
(18) For this method, the whole volume Ω f is divided into a set of of discrete elements, here called discrete volumes (κ f ,e ) e=1,n el covering the whole domain (Ω f = ∪ n el e=1 κ f ,e ) without overlapping (∩ n el e=1 κ f ,e = ∅). For a Finite Element discretization, the solution spaces are restricted to suitable spaces of piecewise polynomial functions over the set of discrete elements. For a Finite Volume discretization used herein, the test functions are chosen in the space of characteristic discrete volume functions. For instance, the velocity can be approximated as:

V h = V ∩ v v | κ ∈ span(ι κ ), ∀κ ∈ T h (Ω f ) ( 19
)
where ι κ is the characteristic function of the element defined as:

ι κ : Ω f -→ | x -→ 1 if x ∈ κ 0 if x ∈ Ω/κ (20)
The same kind of restriction holds for the solution spaces of equation [START_REF] Deparis | Fluidstructure algorithms based on Steklov-Poincaré operators[END_REF]. Therefore, the function are piecewise constant by elements, and with the restriction of the weak formulation reduces to: find

(u h m , ι h , v h , p h ) ∈ U h × I h × V h × V h , such that,for all (δu h m , δι h , δv h , δp h ) ∈ U h 0 × I h 0 × V h 0 × V h 0 : G f ((u h m , ι h , v h , p h ); δu h m , δι h , δv h , δp h ) = 0 .
The divergence terms in equation [START_REF] Deparis | Fluidstructure algorithms based on Steklov-Poincaré operators[END_REF] can be written in terms of flux at the boundary of volume controls using the Gauss theorem. Hence, the weak formulation can be written as:

0 = κ ∂κ dΓ • (γ∇u m ) -[B.C.] 0 = κ κ ∂ t ι + ∂κ dΓ • (v -um ) ι -[B.C.] 0 = κ κ ρ∂ t v + ∂κ ρ dΓ • (v -um ) ⊗ v - ∂κ dΓ • µ f D(v ) + ∂κ p dΓ -[B.C.] 0 = κ ∂κ dΓ • v -[B.C.]
where dΓ is the elementary surface vector. Note that there is no continuity requirement for the solution (contrary to classical FE), and therefore the approximate solutions are not properly defined at the interface. Their flux can be computed without being imposed by the restriction of the solution space to the FV space. The only difficulty is now to build an accurate representation of the fluxes at the boundaries from a piecewise constant field.

On each control volume, three levels of numerical approximations are applied to build the boundary fluxes: interpolation: to express variable values at the control volume surface in terms of nodal values (depending on where the variable is stored).

differentiation: to build convective and diffusive fluxes the value of the gradient of the quantity of interest -or at least its approximation -is required.

integration: to approximate surface and volume integral using quadrature formulae.

The details of constructing an accurate representation of the derived fields (especially of the convection terms that require often special cares) will not be exposed here. Some schemes seek to stabilize the solution at the expense of accuracy, others, such as the MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) exhibits good stability properties for a high-order scheme [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF][START_REF] Piperno | A nondiffusive finite volume scheme for the three-dimensional Maxwell's equations on unstructured meshes[END_REF].

The fluid mesh motion considers that u m is imposed by the motion of the interface u and can be written in a semi discrete form as:

R m (u m ; u) := K m u m -D m u = 0 (21)
where D m is a projection/restriction operator and K m governs the extension of the boundary displacement either by a diffusion process or a pseudo-solid equation (see Section 3.2 above). The volume fraction ι, the d components of velocity v and pressure p are coupled through with a set of non-linear equations.

Written in a matrix forms, it gives the following semi-discrete problem:

R f (ι, v f , p f ; u m ) :=   M ι ι + N ι (v f -um )ι M f (ι) vf + N f (ι, v f -um )v + K f (ι)v f + B f p f -f f (ι) B T f v f   = 0 (22)
where M ι and N ι are the matrices associated to the advection problem of the fluid volume fraction, M f is a positive definite mass matrix, N f is an unsymmetrical advection matrix, K f is the conduction matrix describing the diffusion terms, and B f stands for the gradient matrix, whereas f f is the discretized nodal loads on the flow. This matrix form also takes into account the boundary conditions; special care has to be taken concerning the discretization of boundary conditions -and especially normal flux -when using the Finite Volume Method [START_REF] Ghidaglia | The normal flux method at the boundary for multidimensional finite volume approximations in CFD[END_REF].

Semi-implicit non-linear solver for the fluid problem

One way to solve this problem is to consider a monolithic application. Another way is to consider a split between the mesh motion, the volume fraction advection, the momentum and the continuity equations, and to solve it thanks to an operator split-like procedure often termed as the segregated approach [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. This approach is favored for its computational efficiency compared to the monolithic one. Indeed, even with a simple fixed point iteration strategy its cost is less important than that of the monolithic approach for large size problems [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. In the work presented herein, the later will be used.

Let us note that for a given motion of the fluid domain, the coupling between the mesh deformation and the Navier-Stokes equation is weak, in the sense that no variable like velocity v or pressure p influences the fluid domain deformation under imposed boundary displacements. The coupling between the mesh motion problem and the fluid momentum equation can therefore be insured explicitly.

The only remaining question is the choice of velocity in the time step ∆t = t N +1t N . As the mesh motion is arbitrary and does not rely on any physical phenomenon, it is a priori possible to take any velocity evolution on the window [T N , T N +1 ] so that the initial mesh deformation is equal to u m,N and the final mesh deformation is equal to u m,N +1 .

The Geometric Conservation Law demands a numerical scheme to reproduce exactly and independently from the mesh motion a constant solution. This condition can be found in the literature for ALE formulation discretized either by the Finite Volume [START_REF] Demirdžić | Space conservation law in finite volume calculations of fluid flow[END_REF] or the Stabilized Finite Element methods [START_REF] Förster | On the geometric conservation law in transient ow calculations on deforming domains[END_REF]. It is proven [START_REF] Farhat | The discrete geometric conservation law and the nonlinear stability of ale schemes for the solution of flow problems on moving grids[END_REF] that the velocity of the dynamic mesh needs to be computed for all first-and second-order time accurate methods like implicit Euler or Crank-Nicholson that the mesh motion velocity so that:

um = u m,N +1 -u m,N ∆t (23)
The volume fraction function is supposed to be sharp at the interface between water and gaze, and therefore, standard FV discretization that can be strongly diffusive cannot be applied, as they would smear the interface. A way to guarantee a sharp and bounded solution is to solve it with a numerical scheme designed for the multi-dimensional advection equation [START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF]. We will not enter into the details of such a treatment, but let us note that, in our case, the treatment of the volume fraction require time sub-cycles, as an explicit treatment termed as MULES (Multidimensional Universal Limiter with Explicit Solution) is used (see for instance [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF][START_REF] Behzadi | Modelling of dispersed bubble and droplet flow at high phase fractions[END_REF]).

For the coupling between the pressure equation and the momentum equilibrium, strategy based upon ACM (Artificial Compressibility Method) [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF], or pressure correction techniques such as SIMPLE (Semi-Implicit Method for Pressure Linked Equations) [START_REF] Patankar | A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[END_REF][START_REF] Van Doormaal | Enhancements of the SIMPLE method for predicting incompressible fluid flows[END_REF] or PISO (Pressure Implicit with Splitting of Operators) [START_REF] Issa | Solution of the implicitly discretised fluid flow equations by operator-splitting[END_REF][START_REF] Issa | Solution of the implicitly discretised reacting flow equations by operator-splitting[END_REF] are traditionally used in CFD. They often rely on the use of suitable relaxation parameters [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF] in order to reach convergence for the stiff coupled problem. In [START_REF] Rogers | Steady and unsteady solutions of the incompressible Navier-Stokes equations[END_REF], a comparison between ACM and pressure-correction techniques is given and in [START_REF] Barton | Comparison of SIMPLE-and PISO-type algorithms for transient flows[END_REF], a comparison of two pressure-correction algorithms showed the overall better performances of PISO-like algorithms over SIMPLE ones.

In this work, PISO-like algorithms are used to solve our CFD problem. The semi-discrete form of the Navier-Stokes equation is discretized in time using implicit or explicit integration schemes, such as Euler explicit and implicit, or a second order Crank-Nicholson scheme. The discretized momentum Eq. ( 22) is split in the following way when an implicit integration scheme is used:

A f (v N +1 )v N +1 -H f (v N , v N +1 ) = -B f p N +1 (24)
where A f stands for time derivative terms in a cell (and is therefore diagonal) and H f takes into account all neighboring velocity and source terms in elements.

The incompressibility condition can be re-written in a discrete form using the previous split as:

B T f 1 A f (v N +1 ) B f p N +1 -B T f 1 A f (v N +1 ) H f (v N , v N +1 ) = 0 (25)
For the PISO algorithm, the coupling between the incompressibility condition and the momentum equilibrium parts of the Navier-Stokes equation is assured in an iterative way as detailed in Algorithm 1.

The PISO algorithm is not fully implicit, as corrective term of the velocity is introduced explicitly. Hence, in the correction step, it is supposed that the influence of the transported term is negligible compared to the pressure gradient correction terms. Therefore, even with an implicit time integration scheme, the stability of the PISO algorithm remains conditional, and when the Courant Number becomes too large (it means that the transport due to the flow over a cell is not well captured) the PISO algorithm fails to converge.

Remark: in the algorithm given, the non-orthogonal correctors applied to build more accurate flux terms when the mesh is not orthogonal grid are not detailed. For further details, see [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF][START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF].

Algorithm 1 Pressure Implicit Splitting of Operators algorithm 1: Given: initial velocity v 0 and pressure p 0 . 2: while T < T max do

3:

Momentum predictor:

A f (v (0) N +1 )v (0) N +1 -H f (v N , v (0) N +1 ) = -B f p N 4: for (k) = 0 to (k) < (k max ) do 5:
Pressure correction:

B T f • 1 A f (v (k) N +1 ) B T f p (k) N +1 = B T f • H f (v N , v (k) N +1 ) A f (v (k) N +1 ) 6: Explicit velocity correction v (k+1) N +1 = v (k) N +1 - 1 A(v (k) N +1 ) B f p (k) N +1 7:
end for 8: end while 4 Implicit coupling using algorithm Direct Force-Motion Transfer and Block-Gauss-Seidel with Aitken's relaxation

By enforcing the continuity of primal variables at the interface we can eliminate the energy errors that characterize the explicit interface matching. When coupling incompressible flow with structure, the implicit interface matching is required for stability reason, as proved in [START_REF] Kassiotis | Partitioned procedure for strongly coupled fluid-structure problems. Part I: implicit algorithm, stability analysis and validation examples[END_REF]. This ought to be done by iterating on the following residual to reduce its value below the chosen tolerance:

r N +1 := u s,N +1 -u f ,N +1 ≃ 0 ≤ TOL (26)
In this way we obtain an implicit algorithm requiring more than one iteration to enforce the interface matching condition. The chosen order of iterations, corresponds to the Block-Gauß-Seidel algorithm for fluid-structure interaction problem [START_REF] Matthies | Algorithms for strong coupling procedures[END_REF]. Let us note that not only the value at synchronization points T n or T n+1 , but also the interpolated evolution of variables have to be exchanged in the entire time-interval t ∈ [T n , T n+1 ] when the time steps are not matching between fluid and structure sub-problems.

Contrary to explicit algorithms which generate spurious energy at the interface, the present implicit interface matching algorithm enforce the same evolution of the primal variables at the fluid-structure interface. In other words, an iterative solution for primal (displacements) continuity as well as the dual (forces) equilibrium equations at the interface is performed by using the Picard iteration:

u (k+1) N +1 = G u (k) N +1 ; G = S s -1 • -S f ( 27
)
where S f and S s are Steklov-Poincaré operators for fluid and structure defined as define in [START_REF] Deparis | Fluidstructure algorithms based on Steklov-Poincaré operators[END_REF]. These operators can be formulated using transfer operators and equations ( 6), ( 22) and ( 21):

S s = T λ s • R s • T u s ; S f = T λ f • R f • R m • T u s ( 28
)
where the transfer of structure displacement to fluid-structure interface displacement is T u s , the transfer of fluid displacement to interface displacement T λ s , the transfer of structure stresses to the interface T u f and the transfer of fluid stresses to the interface T λ f . The Picard iterations will continue until convergence of interface residual is achieved:

r (k) N +1 = u (k) s,N +1 -u (k) f ,N +1 = G u (k) N +1 -u (k) N +1 (29)
Such a Block-Gauß-Seidel algorithm with implicit interface matching, further denoted as DFTM-BGS, can be presented as a natural generalization of the explicit algorithms. We can thus write: Algorithm 2 Direct Force-Motion Transfer Block-Gauß-Seidel 1: Given: initial time T = T 0 , final time T max , window size ∆t, initial interface displacement u 0 . 2: while T < T max do 3:

(k) = 0 4: Predict displacement: u (0) N +1 = P(u (kmax) N , u(kmax) N , u (kmax)
N -1 , . . . )

5:

repeat 6:

Perform Picard iteration: G u (k) N +1

7:

Compute residual: r

(k) N +1 = G u (k) N +1 -u (k) N +1 8:
Update interface primal variable: u (k+1)

N +1 = u (k) N +1 + r (k) N +1 9:
do (k) ← (k) + 1 N ← N + 1 and T ← T + ∆t 12: end while It is clear that this fixed-point algorithm based on Picard iterations has the main drawback that the search directions for u and λ variables at the interface do not exploit any information from the fixed-point function G nor the Steklov-Poincaré operators S f and S s . Therefore, quite a few iterations may be needed to reach the convergence.

The stability of such a coupling algorithm is studied in [START_REF] Kassiotis | Partitioned procedure for strongly coupled fluid-structure problems. Part I: implicit algorithm, stability analysis and validation examples[END_REF]. We give a formal proof of potential numerical instability due to the added-mass effect also observed in [START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF][START_REF] Tallec | Fluid structure interaction with large structural displacements[END_REF]. In order to improve the convergence of the DFMT-BGS method, we can use a relaxed update:

u (k+1) N +1 = u (k) N +1 + ω (k) r (k) N +1 (30)
Our favorite choice for constructing ω (k) is using a secant methods which can keep the cost of each iteration as low as possible. The Aitken's relaxation strategy has been extensively used in fluid-structure interaction [START_REF] Deparis | Fluidstructure algorithms based on Steklov-Poincaré operators[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF], and shown sufficient performances to be used in the following.

The use of different different solvers, for the fluid and the structure part, do not provide in general a matching mesh at the interface. Furthermore, even for matching meshes, as the geometries of the domains are not the same on both sides of the interface, an optimal numbering of the nodes can lead to different orders for the interface nodes. In the examples proposed herein, only this latter point is of interest. Last but not least, different discretization techniques (Finite Element versus Finite Volume) or different order p of the polynomials can be used for constructing solution to fluid-structure interaction problem. In the domain of FE applied to mechanical engineering, extensive literature can be found on how to build a consistent interpolation for both subproblems at the interface [START_REF] Felippa | Synthesis tools for structural dynamics and partitioned analysis of coupled systems[END_REF]. For the fluid-structure interaction problems, an interesting review can be found in [START_REF] De Boer | Review of coupling methods for non-matching meshes[END_REF]. In our framework, it was decided not to favor any mesh-based representation of the interface, since, in the most general case, the fluid problem can also be solved by a meshfree-based method [START_REF] Dalrymple | Numerical modeling of water waves with the SPH method[END_REF]. Namely, an interpolation strategy relying on radial basis function is here chosen. This method has already been employed for FSI in [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Beckert | Multivariate interpolation for fluid-structureinteraction problems using radial basis functions[END_REF].

The algorithm presented here is simple to implement. We use for this work the Communication Template Library (CTL, see [START_REF] Markovič | Multi-scale modeling of heterogeneous structures with inelastic constitutive behavior : Mathematical and physical aspects[END_REF][START_REF] Niekamp | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior: Part II-software coupling implementation aspects[END_REF]) that allows to re-use existing codes in a generic way, either called as libraries on the same computer, or as remote executables through network. With the CTL, we are able to couple existing stand-alone software, in a quite straightforward way, even if they are programmed in different langages (Fortran for the structure part, C++ for the fluid part), and to conserve the inner parallelism of each component. For more details on the implementation, the reader is invited to see [START_REF] Kassiotis | Partitioned procedure for strongly coupled fluid-structure problems. Part II: implementation aspects, re-using existing software as components and nested parallel computations by using ctl[END_REF].

5 Numerical simulations

Two-dimensional sloshing wave hitting a rigid structure

To validate the proposed V.O.F. fluid strategy, the classical dam break example is computed. It is studied from the experimental [START_REF] Bullock | Violent breaking wave impacts. Part I: Results from large-scale regular wave tests on vertical and sloping walls[END_REF], the analytical [START_REF] Stoker | Water Waves: The Mathematical Theory and Applications[END_REF] as well as the numerically point of view [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF][START_REF] Dutykh | On the relevance of the dam break problem in the context of nonlinear shallow water equations[END_REF], to describe the collapse of a water column hitting a rigid structure. The geometry of the problem is represented in Fig. 3: it mainly consists on a 400mmmm box that contains a water column (initially 146mmmm) collapsing and hitting a 16mmmm rigid obstacle. The high density fluid representing water has the following properties: density ρ 1 = 1 × 10 3 kg.m -3 and kinematics viscosity ν 1 = 1 × 10 -6 m.s -2 . The low density fluid -air -is described with a density ρ 2 = 1kg.m -3 and a kinematics viscosity ν 2 = 1.48 × 10 -5 m.s -2 . The flow in each phase is considered to be Newtonian (no turbulence effect is taken into account). The surface tension between each phase is 0.07kg.s -2 , but regarding the characteristic size of the problem, the surface tension is of little influence and can be neglected.

The computation is run with three meshes, with 2268, 9072 and 36288 cells: around 46 and 50 cells in e x and e y directions for the coarse grid, 92 cells for the medium grid, and 92 cells are considered fine grid. At t = 0, the initial conditions are at rest, and the water column starts to collapse. The computation is run with an time step ∆t = 0.005. The time step can be decreased in order Time t = 0.1s Time t = 0.2s Time t = 0.4s Time t = 0.6s Figure 3: Evolution of the volume fraction ι < 0.3 for the two-phase flow (fine mesh) and comparison with experimental results from [START_REF] Koshizuka | A particle method for incompressible viscous flow with fluid fragmentation[END_REF] to maintain the Courant-Friedrichs-Lewy condition under 0.5 during the whole simulation. For the coarsest grid considered, no decreasing of the time step is required whereas it is required for the finest grid after the water column hits the rigid obstacle.

In Fig. 3, the evolution of the flow for the finest grid, and a comparison with experimental results [START_REF] Koshizuka | A particle method for incompressible viscous flow with fluid fragmentation[END_REF] is given. We give here only this qualitative comparison. To get more quantitative comparison, one can see the same dam-break problem without obstacle that as been largely used to validate the VOF strategy (see e.g. [START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF][START_REF] Dutykh | On the relevance of the dam break problem in the context of nonlinear shallow water equations[END_REF]). As we are interested in fluid-structure interaction, the dambreak with an obstacle seems here more relevant with the examples proposed in the following sections. In Fig. 4 the forces acting upon the obstacle are also represented. Before the water column hit the structure, the three meshes give roughly the same response. After that, the fine mesh is able to capture the separation of the water phase in many drops and to represent more accurately their separated impact on the right side of the obstacle. It explains the difference observed for t ≥ 0.5. The use of a turbulence model should be here probably required to obtain more physical results, but at this is not the topic of this paper (e.g. see [START_REF] Launder | The numerical computation of turbulent flows[END_REF][START_REF] Pope | Turbulent flows[END_REF]), we will use only laminar models with the coarse and medium mesh densities for the fluid-structure interaction examples provided in the following.

Two-dimensional sloshing wave hitting a flexible structure

The problem solved herein is a modification of the dam-break problem presented in Sec. 5.1. At initial time t = 0s, the same water column starts to fall under the gravity loading. Instead of hitting a rigid structure, its geometry is modified in order to get a more slender obstacle (Fig. 5). Furthermore, the material properties of the obstacle are modified in order to obtain a flexible elastic structure. This problem was studied experimentally and numerically [START_REF] Walhorn | Fluid-structure coupling within a monolithic model involving free surface flows[END_REF][START_REF] Baudille | Modelling FSI problems in FLUENT: a dedicated approach by means of UDF programming[END_REF].

Ω f ,1 ρ f,1 , µ f ,1 Ω f ,2 ρ f,2 , µ f ,2 ρ s , E s , ν s Ω s g
p + 1 2 ρ v 2 = p atm (31)
with p atm = 0. The material properties are imposed as follows: the density and the kinematic viscosity are ρ f ,1 = 1 × 10 3 kg.m -3 and ν f ,1 = 1 × 10 5 m.s -1 for the high density fluid, whereas ρ f ,2 = 1kg.m -3 and ν f ,2 = 1 × 10 6 m.s -1 are considered for the low density fluid domain. Due to the scale of the simulation there is no need to consider surface tension.

Mesh motions based on smoothing operator like the Laplacian equation fails for the used meshes since the points around the obstacle have difficulties to follow the large displacement and rotations of the mesh around the structure. For this reason, the mesh motion problem is solved by using a pseudo-elastic material model where the rigidity is a quadratic inverse function of the distance to the interface between solid and fluid.

The fluid problem is discretized with Finite Volume Method. The results are presented for two meshes with respectively 3340 and 13760 cells. The fluid is handled by second order space discretization and a Van Leer limiter is used for the advected terms. The time integration scheme employed herein is implicit Euler. Small time steps are required by the explicit nature of the coupling between the phase function indicator problem, the momentum prediction and the pressure correction step.

For structure mechanics problem, it is proposed to use two dimensional elements with quadratic shape functions. Each element therefore contains 9 nodes. For the coarse grid, 51 nodes are considered while 165 nodes are used to discretize the fine grid. The total number of d-o-f given in Tab. 1. The model can handle finite deformations. The material properties used for the solid are: a neo-Hookean elastic material with Young modulus E s = 1 × 10 6 P a and Poisson's ratio ν s = 0. The chosen value of mass density

ρ s = 2500kg • m -3 .
The time integration is handled by a Generalized-α scheme with: The computation is carried out run with a time step of 1 × 10 -4 for the coarse and 2 × 10 -5 for the fine discretization. The coupling scheme used here is DFMT-BGS with Aitken's relaxation, with the initial parameter value of ω = 0.25. The predictor is of order 1, since computations with second order predictor fail for the finnest grid. The absolute tolerance considered is:

ρ ∞ = 1
r (k) N ≤ 1 × 10 -6
The total number of iteration required to reach the convergence criteria at each time step is given in Fig. 6. Note that no iteration is required before the water hits the structure since the effect of air flow can almost be deemed negligible for this structure. Subsequently, the number of iteration depends on the discretization density. Namely, the finest scale is able to represent fine water sloshing after the breaking of the dam, that will increase the number of iterations needed to ensure the convergence of the implicit coupling. This large number of iteration make the computation costly, and for the finest grid, the total time required to perform the whole coupled simulation on a single 3.0GHz Intel processor is 143 × 10 3 s.

In Fig. 8 we represent the domain occuped by high density fluid as well as the stream-lines in the high and low density fluid domains (compare again the results in Fig. 3 for rigid obstacle). For first 0.1s of the simulations, the water column is falling under the gravity loading. There is no effect whatsoever on the structure until the high density fluid reaches it. The maximum amplitude of the motion is obtained at t = 0.25s, before the solid comes back to its initial position due to flow friction.

After one second of simulation, the fluid is not entirely at rest, but its main effects on the structure are mostly captured. Moreover, we are able to catch the Figure 7: Obstacle extremity displacement in e x direction for different meshes and comparison with results from literature main impact with a good accuracy with both the coarse and the fine mesh. Our results are in accordance with the one obtained using a monolithic stabilized-FEM strategy described in [START_REF] Walhorn | Fluid-structure coupling within a monolithic model involving free surface flows[END_REF]. In [START_REF] Baudille | Modelling FSI problems in FLUENT: a dedicated approach by means of UDF programming[END_REF], a tight coupling strategy between a FV and a FEM solver is used, that probably explained the differences observed.

After, the impact high density fluid is highly fragmented (see Fig. 8 for time greater than 0.3s), and the results mostly depend on the ability of the chosen fluid mesh to represent these fine effects, as well as the use of turbulence fluid models not explored herein.

Three-dimensional sloshing wave impacting a flexible structure

The problem solved in this example is a modified three-dimensional representation of dam-breaking event that brings about a sloshing wave impact on a flexible structure presented in Fig. 9. At initial time t = 0s, a three-dimensional water column starts falling down under the gravity loading and eventually hits the obstacle placed in the way. The flexible obstacle is a slender plate-like body made of elastic material that can undergo large deformation. The chosen dimension of the problem, as well as the boundary conditions are given in Fig. 9.

Let us note that we propose to use open boundary conditions far from the obstacle in order to avoid the water bounces-back and hits again the structure after breaking off the walls. For that reason, only the left and bottom planes of the fluid domain are defined as non-slipping walls, while the others are defined The material properties are chosen as follows: the density and the kinematic viscosity are ρ f ,1 = 1 × 10 3 kg.m -3 and ν f ,1 = 1 × 10 5 m.s -1 for the high density fluid (water in the reservoir), versus ρ f ,2 = 1 × 10 3 kg.m -3 and ν f ,2 = 1 × 10 6 m.s -1 for the low density fluid (air in the remaining part of the domain). The mesh motion problem is solved by using a Laplacian smoothing material where the diffusion coefficient is a quadratic inverse function of the distance to the interface between solid and fluid.

The results are computed for two meshes with the chosen discretization and the number of cells given in Tab. 2. For the finest grid, around 64, 60 and 40 cells are used in e x , e y and e z direction. The mesh is refined gradually refined around the structure, and initially the cell dimensions are between 8.9 × 10 -8 and 3.2 × 10 -5 . For this finest grid, the maximum skewness of the mesh observed is 2.947, that does not generate too large errors. The fluid is handled by second order space discretization and a Van Leer limiter is used for the advection terms. The time integration scheme employed in this computation is implicit Euler. For such a scale of modeling it is not required to consider surface tension between the two fluids. For this problem the fluid computation is parallelized, but reduction of the CPU time is obtained by using a Generalized Algebraic-MultiGrid (GAMG) linear solver.

Note that small time steps are required for the explicit solution of the phase function indicator equation, as well as the half-implicit nature of the coupling between the momentum predictor and the pressure corrector.

For the structure part, we propose here to use three-dimensional elements with quadratic shape functions, where each element has 27 nodes. The material properties used for the solid are: a neo-Hookean elastic material with Young's modulus E s = 1 × 10 6 P a and Poisson's ratio ν s = 0 and a density ρ s = 2500kg • m -3 , which can represent finite deformation. The time integration is carried out by a Generalized-α scheme with the same parameters as the one used for the 2D case. The computation of the coupled problem (with total number of d-o-f given in Tab. 2), is carried out by an implicit iterative scheme. The results of fluid and solid computations are matched for a time step of 1 × 10 -4 for the coarse and 2 × 10 -5 for the fine discretization. The coupling scheme used is DFMT-BGS with Aitken's relaxation. The initial parameter is ω = 0.25 and the predictor is of order 1. The absolute tolerance considered is:

r (k) N ≤ 1 × 10 -6 (33)
The number of iterations required to reach the convergence criteria is given in Fig. 10. Note that there is no iteration required before the water hits the strategy is shown to provide reliable (convergent) results with implicit scheme computation for both structure and fluid component. The Direct Force-Motion Transfers algorithm is used to pass the pertinent information between structure and fluid motion, featuring also the implicit solution strategy based upon the fixed-point iterations with a dynamic relaxation parameter used to accelerate convergence [START_REF] Kassiotis | Partitioned procedure for strongly coupled fluid-structure problems. Part I: implicit algorithm, stability analysis and validation examples[END_REF].

Most importantly, the proposed solution method allows to perform coupled simulations and obtain reliable solution to complex fluid-structure interaction by using the existing codes, that were initially developed to support either fluid or structure motion computation. This is achieved thanks to the use of the component technology [START_REF] Kassiotis | Partitioned procedure for strongly coupled fluid-structure problems. Part II: implementation aspects, re-using existing software as components and nested parallel computations by using ctl[END_REF][START_REF] Niekamp | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior: Part II-software coupling implementation aspects[END_REF] providing the coupling between existing software products. Therefore, the proposed solution method for fluid-structure interaction can utilize very different discretization strategies to obtain the optimal accuracy; The case in point concerns FE for the structure and FV for the fluid. The use of these popular methods for the fluid and solid parts allows to benefit from the advanced features of the two families of methods, each developed by the experts from the corresponding domain. Accordingly, on the fluid side, it is possible to use a very efficient semi-implicit solver for incompressible flow (PISO), inverse techniques (Algebraic Multigrid) or advanced models for free-surface flows. A very good performance of the proposed technology for fluid-structure interaction is illustrated with 2D and 3D models for dam breaking examples, which also involve flexible obstacles.

K

  s,e (u e )u e = κe ∇N : σ (u e N e )

Figure 2 :

 2 Figure 2: Techniques to solve a free-surface flow

Figure 4 :

 4 Figure 4: Horizontal and vertical forces acting upon the obstacle for coarse (2268 cells) and fine (9072 cells) grids

Figure 5 :

 5 Figure 5: Dam break interacting with an obstacle: geometry (given in mm) and boundary conditions All boundaries a non-slip boundary condition is applied, except for the upper boundary z = H where the tank is open and a total pressure condition is used:

Figure 6 :

 6 Figure 6: Number of iteration in order to make the DFMT-BGS algorithm converge for the two-dimensionnal dam-break example

Figure 8 :

 8 Figure 8: Bi-dimensional dam breaking problem. Evolution of the free surface, structure motion and streamlines for water and air.

Figure 9 :

 9 Figure 9: Three-dimensional water column impacting an obstacle: geometry (given in mm) and boundary conditions

5 × 10 5 Table 2 :

 552 Fine 104 × 10 3 520 × 10 3 2205 6.6 × 10 3 1 Number of d-o-f for coarse and fine discretization of the threedimensional dam-breaking problem

Figure 12 :

 12 Figure 12: Tri-dimensional dam break problem. Evolution of the free surface and motion of the structure.

Table 1 :

 1 Number of d-o-f for coarse and fine discretization of the twodimensional dam-break problem

		2	; β =	4 9	; γ =	5 3	and α =	2 3	(32)
	At each iteration, the linear system is solved by a direct solver for real value
	asymmetric matrices.
			fluid			solid		number of
	Discretization	cells		d-o-f nodes d-o-f time steps
	Coarse	3440 17.2 × 10 3	51	102	1 × 10 4
	Fine 13760 68.8 × 10 3	165	330	5 × 10 4

We note in passing that the formulation developed herein for fluid-structure interaction problem would also apply to more elaborate inelastic constitutive models (see[START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF])

See also the reference thread in http://www.cfd-online.com/Forums/openfoam-solving/ 57916-moving-mesh-problem-openfoam-141-a.html OpenFOAM forum

Figure 10: Number of iterations in order to make the DFMT-BGS algorithm converge for the three-dimensional dam-breaking problem structure since the effect of air flow can almost be deemed negligible with respect to the structure. Subsequently, the number of iterations depends on the discretization density. In reaching the opposite wall, the water does not rebound on the wall but simply flows away. For the finest grid, the total time required to perform the whole coupled simulation on a single 3.0GHz Intel processor is 279 × 10 3 s.

In Fig. 12, the high density fluid domain is represented, as well as some part of the fluid mesh and the structure displacement. The first 0.1s of the simulations, the water column falls under the gravity loading. There is no effect whatsoever on the structure until the high density flow reaches its bottom. The maximum amplitude of the motion is obtained at t = 0.25s, before the solid comes back to its initial position and oscillates after the shock. In Fig. 11 the motion of the extremity of the solid obstacle is plotted. Contrary of the two-dimensional example, small drops of high density fluid are not interacting with the obstacle after the main shock. Therefore, the motion of the flexible structure remains fairly smoothed and it is rather well described with the coarsest grid.

Conclusion

The presented work deals with a partitioned strategy for complex fluid-structure interaction problems in presence of free-surface flows. The proposed coupling