Alban Leroyer
email: alban.leroyer@ec-nantes.fr

Camille Yvin
email: camille.yvin@sirehna.com

Emmanuel Guilmineau

Michel Visonneau

P Queutey

Efficient FSI codes coupling with possible large added mass effects: applications to rigid and elongated flexible bodies in the maritime field

Keywords: Fluid-Structure Interaction, added-mass, hydrodynamics, codes coupling

INTRODUCTION

Fluid Structure Interaction (FSI) problems are commonly encountered in naval architecture. Even if this can be done through a monolithic approach within a single solver, a partitioned approach is most commonly used through codes coupling, especially well fitted to the resolution of complex FSI problems. Indeed, with this technique, complex models for both the fluid and the structure can be used because each solver can be numerically adapted and dedicated to its own physics. In this work, two different general solvers are used: ISIS-CFD and MBDyn, for the fluid and the structural part, respectively. 1 1 ISIS-CFD is a Navier-Stokes solver developed in the LHEEA Laboratory at Ecole Centrale de Nantes. MBDyn is an open-source solver intended to solve multi-disciplinary problems including non-linear dynamics of rigid and flexible bodies subjected to kinematic constraints, along with active controls. The combination of these two solvers makes possible the study of a broad spectrum of applications in the marine field which cannot be solved with a unique solver. After describing how the code coupling is handled and how the added-mass effects are tackled to reach a robust and efficient algorithm, some applications among those studied are shown. Especially, a test-cases of a ship with active control of appendages (roll damping) is described, where good agreement with the experimental data was obtained. Other applications with flexible slender bodies are also presented to demonstrate the capabilities of such a coupling.

DESCRIPTION OF THE TWO SOLVERS

2.1 The fluid part: ISIS-CFD ISIS-CFD is available as a part of the FINE TM /Marine computing suite which is dedicated to marine applications. This is an incompressible unsteady Reynolds-averaged Navier-Stokes (RANS) solver developed by the DSPM group of the LHEEA Lab. of Ecole Centrale Nantes, UMR-CNRS 6598 . This solver is based on a fully unstructured finite-volume method to build the spatial discretisation of the conservation equations. Pressure-velocity coupling is obtained through a Rhie & Chow SIMPLE-type method : in each time step, the velocity update comes from the momentum equations and the pressure is given by the mass conservation, transformed into a pressure equation. An Arbitrary Lagrangian Eulerian (ALE) formulation is used to take into account modification of the fluid spatial domain [START_REF] Leroyer | Experimental and numerical investigations of the flow around an oar blade[END_REF]. It is associated with robust and fast grid deformation techniques [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF][START_REF] Durand | Interaction fluide-structure souple et légère, application aux voiliers[END_REF]. The temporal discretisation scheme is the Backward Difference Formula of order 2 (BDF2) when dealing with unsteady configurations. For each time step, an inner loop (denoted by non-linear loop) associated to a Picard linearisation is used to solve the non-linearities of the system and to converge all the sequential coupled equations.

Free-surface flow is addressed with an interface capturing method, by solving a convection equation for the volume fraction of water, which is discretised with specific compressive discretisation schemes [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF]. The code is fully parallel using the MPI (Message Passing Interface) protocol. An automatic adaptive grid refinement technique [START_REF] Wackers | Adaptive grid refinement for hydrodynamic flows[END_REF] and a sliding grid method are also included .

The structure part: MBDyn

MBDyn (Multi-Body Dynamics), is an open-source solver under the GNU GPL license developed at the Dipartimento di Ingegneria Aerospaziale of the Politecnico di Milano. It is aimed at the modelling of complex multi-bodies systems and muti-discipline problems including non-linear dynamics, aero-servo-elasticity, smart piezo-structural components and electric and hydraulic components [START_REF] Masarati | MBDyn Theory and Developer's Manual Version 1[END_REF] To solve the kinematic laws of a multi-body mechanical system, the Redundant Coordinate Set (RCS) formulation is used. This means that every inertial body has six rigid body Degrees of Freedom (DOF) even if they are constrained by joints for instance. Additional holonomic or nonholonomic constraint equations are added which introduce algebraic unknowns that are analogous to the Lagrange multipliers and directly represent the reaction forces and couples [START_REF] Masarati | Comprehensive multibody AeroServoElastic analysis of integrated rotorcraft active controls[END_REF]. All these equations are written in the form of a set of first order Algebraic Differential Equations (ADE). A direct resolution of these multi-body ADE is made possible through a combination of a variable substitution of the algebraic unknowns, the scale of constraint equations by the time step and a A-stable multi-step time integration scheme to damp numerical oscillations. In this work, the BDF2 scheme, similar to the fluid solver, is used.

The ability to take into account multidisciplinary complex systems and the simplicity of implementation are the main advantages of this formulation.

CODES COUPLING AND FSI

General algorithm

FSI problems can simply be expressed as two continuum materials (structural and fluid, denoted by an index s and f , respectively) sharing a common interface Γ s = Γ f = Γ, where the two following conditions operate: kinematic condition: δ s = δ f on Γ , (δ refers to the position)

(1) dynamic condition: σ s + σ f = 0 on Γ , (σ refers to the stress vector)

(

Using a domain decomposition point of view, let's introduce the Steklov-Poincaré operator S and its inverse S -1 , which can be seen as the transfer function at the interface of the fluid and structure solver (subscript d represents the considered domain):

S d (δ d ) = σ d S -1 d (σ d ) = δ d , (3)
The conditions (1) and ( 2) can be expressed with the operators previously defined in different formulations which provide different resolution strategies ([8]). In particular, the classical fixed-point algorithm, leading to an implicit Block Gauss-Seidel (BGS) approach (also called Dirichlet-Neumann algorithm) can be simply represented with the following expression:

σ f = -S f • S -1 s (σ f ) (4)
Equation ( 4) is a non-linear time dependant equation and has to be fully solved at each time step. Explicit coupling schemes in time, which perform only one iteration to solve equation (4) using predictor and corrector, can produce or dissipate energy at the interface and leads to unphysical results if the problem is highly coupled. Moreover, they are very sensitive to the destabilizing added-mass effect [START_REF] Badia | On some fluid-structure iterative algorithms using pressure segregation methods. Application to aeroelasticity[END_REF]. Hence, in this work, only implicit time integration is considered. It means that equation (4) needs to be solved through an iterative procedure, indexed by i. The simplest implicit BGS algorithm to solve time iteration n + 1 can be represented by equation ( 5):

3 σ f | i+1 n+1 = -S f • S -1 s σ f | i n+1 δ| i+1 n+1 (5)
At each coupling iteration, the structure problem is solved using the current fluid loads σ f | i n+1 at time step n + 1 coming from the last iteration. This provides a new position δ| i+1 n+1 , which is used as an input data for the fluid solver to compute new fluid variables at the same time step n + 1, leading to an updated fluid stress field at the interface σ f | i+1 n+1 . This algorithm is quite easy to implement since it only operates at the temporal loop level, but needs a lot of coupling iterations before reaching convergence. In order to reduce the simulation time, Aitken ∆ 2 relaxation technique is often used [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF]. This technique is based on a geometrical approximation (tangent method) of the coupled problem. To be efficient, the solvers have to produce physical results, and consequently, each operator has to be solved accurately. This means that for an implicit coupling algorithm, the fluid problem has to be solved several times at each time step which is not acceptable in term of CPU time. In this work, the classical Steklov-Poincaré operator for the fluid problem is then modified to reduce the CPU time of this part. This new operator, denoted as S * f does not represent a global fluid resolution any more but only one Picard iteration of the fluid solver as already considered in [START_REF] Badia | On some fluid-structure iterative algorithms using pressure segregation methods. Application to aeroelasticity[END_REF]. It can be seen as a fluid linearised version of the Steklov-Poincaré operator. In a few words, the fluid problem converges at the same time as the coupled problem. In term of implementation, it means integrating the resolution of the structure solver inside the non-linear iteration of the fluid which becomes mixed up with the FSI coupling loop. To be efficient, the CPU time related to the structural solver has to be weak compared with that requested for the fluid part. A second requirement is to have available fast grid deformation techniques since such procedure is now called at each non-linear iteration to update the fluid mesh with the new boundary nodal positions related to the solved structures. The final algorithm (including a relaxation operator which is going to be introduced in section 3.2) is displayed on figure 1.

Added-mass effect and stabilization

When dealing with large fluid density or light structure, added-mass effects become an important physical phenomenum for the coupled problem. At the discrete level, this inertia effect which leads to fluid loads dependant on the acceleration of the structure makes all the code coupling algorithms previously described unstable. This numerical instability has been clearly highlighted in [START_REF] Söding | How to integrate free motions of solids in fluids[END_REF][START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF]. When implicit algorithm is used, a popular technique to tackle this diverging behaviour is to modify the structure equation to decrease the acceleration dependancy of the source term. This modification can be interpreted as an approximated (and then iterative) resolution of a block-LU factorization of the monolithic system. Let's represent the formal monolithic system as:

F C SF C F S S x F x S = s F s S (6) 4
x F and x S refer to the fluid and structure variables, respectively. s F and s S include the source term of the fluid and structure operators (F and S, respectively). C SF and C F S define the data transfer operator at the interface between the two domains. A block-LU factorization leads to the following system:

(S -C F S F -1 C SF )x S = s S -C F S F -1 s F (7) Fx F = s F -C SF x S (8)
This resolution exhibits operator composed with F -1 , which is not acceptable to build in the context of co-simulation. However, by reducing the Jacobian of C F S F -1 C SF only to its inertial effect (i.e. by the added mass operator -M a , approximation all the more valid as the time step is small [START_REF] Badia | Splitting methods based on algebraic factorization for fluid-structure interaction, added-mass operator (dt -¿ 0)[END_REF][START_REF] Joosten | Analysis of the block Gauss-Seidel solution procedure for a strongly coupled model problem with reference to fluid-structure interaction[END_REF][START_REF] Durand | Interaction fluide-structure souple et légère, application aux voiliers[END_REF]), it is possible to transform the direct resolution of in an iterative procedure without any operator composed with F -1 . It comes indeed:

(S + M a ) x k+1 S = s k S -C F S x k F + M a x k S ( 9
)
Fx k+1 F = s k F -C SF x k+1 S ( 10
)
Whereas equation ( 10) refers to a classical resolution of the fluid solver, equation ( 9) is a modified version of the structure solver in which the right-hand side term is far less dependent on acceleration. In a code coupling context, this procedure is not suitable since it requires to modify the structural solver. To avoid this, it can be shown that the equation ( 9) is almost equivalent to a classical resolution of the structure solver added to a second step of relaxation with an operator R, depending on the mass matrix of the structure and the added-mass operator (or an approximation of it). In the context of code coupling, the two approaches slightly differ. As a matter of fact, the evaluation of s S and the non-linearities due to the rotation motion are not solved in the same manner, since for the second case, the relaxation occurs after the global resolution of the structure. The final algorithm is then simply described by equation [START_REF] Söding | How to integrate free motions of solids in fluids[END_REF] and is illustrated with a graph in figure 1.

σ f | i+1 n+1 = -S * f • R • S -1 s σ f | i n+1 ( 11
)
For a body with one degree of freedom (DOF), it can be shown that the relaxation value for the acceleration is equal to 1/(1 + ma m ), where ma is an approximation of the added mass. In case of a six DOF rigid body, the interface displacement can be replaced by the generalized position (position and orientation) of the body denoted as δ and the added mass effect can be represented by a symmetric matrix of rank six denoted as M a . The corresponding approximated added mass matrix is denoted as MA. Thus, the acceleration relaxation step is defined by equation [START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF]. δ is the direct result coming from the unmodified structure solver and R a is the acceleration relaxation operator.

δ i+1 n+1 = δ i n+1 + R a δ i+1 n+1 -δ i n+1 , with R a = I d + M -1 MA -1 (12) 5
After the acceleration relaxation step, velocities and generalized positions are reconstructed according to the time integration scheme used. When using the implicit algorithm described in figure 1, the operator MA does not need to accurately represent the physical added-mass operator MA to reach the converged solution. It is also worth noting that this converged solution does not depend on MA. However, it was shown in [START_REF] Yvin | Interaction fluide-structure pour des configurations multi-corps. Applications aux liaisons complexes, lois de commande d'actionneur et systèmes souples dans le domaine maritime[END_REF] that the number of iterations to reach a given gain is noticeably smaller when this operator is close to the physical added-mass operator. In this case, it was found that the number of iteration is even similar to an unsteady configuration without coupled FSI (for example, a simulation with an imposed body motion), which can be seen as an optimal efficiency. Tests have been carried out not only with classical hydrodynamics cases but also with configurations with extreme added mass effects. Even if the relaxation operator can lead to very small value, it does not compromise the efficiency of the coupling since this low relaxation has a physical origin. Let's now describe how the operator MA is computed for a rigid body and how it is generalized for elongated beam bodies.

Computation and approximation of the added-mass operator

This physical added-mass operator, which is based on the instantaneous response to acceleration, generalizes the infinite-frequency added-mass operator deduced from classical linear potential flow approaches. But to avoid any recourse to such an external tool, we developed inside the CFD solver an integrated computation of the added-mass operator, which contributes to the originality of the present work. It is based on the resolution of the pressure field due to a brutal variation of the rigid body velocity (i.e. an acceleration step). Due to the different time scales of the Navier-Stokes equation, this perturbation denoted ũ is guided by the following equation, where p refers to the additional pressure field due to this perturbation [START_REF] Söding | How to integrate free motions of solids in fluids[END_REF][START_REF] Yvin | Interaction fluide-structure pour des configurations multi-corps. Applications aux liaisons complexes, lois de commande d'actionneur et systèmes souples dans le domaine maritime[END_REF].

∂ ũ ∂t == - 1 ρ ∇p (13)
By taking the divergence of equation ( 13) and using the Green-Ostrogradski theorem, it comes the following finite-volume form since ũ is divergence free:

" S - 1 ρ ∇p • n dS = 0 (14)
This equation has to be solved with the Neumann boundary conditions for p with a unit acceleration γ b imposed at the surface of the considered body by a given DOF (one equation to solve for each one) and with zero for the possible other bodies or walls.

-

1 ρ ∇p • n = ∂ ũb ∂t • n = γ b • n at the wall of the considered body (15)
Equation ( 14) is finally similar to the pressure equation used to solve the real flow. The free surface shape is naturally taken into account by the spatial variation of the density ρ 6 6 when a multi-phase flow is considered. No additional equation has to be solved to take into account the free-surface position. This is not the case when a classical potential flow solver is used because the kinematic and the dynamic boundary conditions at the free surface have also to be respected. In addition to the easy-to-use feature due to its integration within the CFD solver (no need to generate another mesh, possible update during the computation), this method does not suffer from any limitations of linear potential solvers. For example, contrary to the latter, fine meshes can be used so all the details of the geometry can be entirely respected, and a larger amount of space configurations, like complex free surface position (wave breaking) or important confinement (shallow water, interactions with close bodies) is naturally taken into account and can be easily updated during the computation, if the conditions are strongly modified.

Case of rigid body

The added-mass operator is here restricted to a 6 × 6 matrix, which can be computed as follows: after resolution of p for a given DOF index i, the integration of p on the body surface provides the six coefficients of the column i of the added-mass matrix from the three forces and three moments. Some validations and comparisons with results coming from potential solver are available in [START_REF] Yvin | Interaction fluide-structure pour des configurations multi-corps. Applications aux liaisons complexes, lois de commande d'actionneur et systèmes souples dans le domaine maritime[END_REF].

Case of flexible elongated body

Assuming the structure of such a body is described with a beam structure, the number of DOF may be quite large compared to a rigid body, since each node of a classical Euler-Bernoulli beam gets at least three DOF of interest for the added-mass effects (the two transversal translation and the rotation of the section along the neutral line). As a consequence, the computation of each pressure field for each DOF becomes too CPU time consuming. Taking advantage of the elongated property of the body, the different added-mass coefficients associated to a given section (and then a given node) can be deduced considering the 2D potential flow in a plane normal to the neutral line as soon as the curvature effects are small ( [START_REF] Candelier | Three-dimensional extension of Lighthill's large-amplitude elongated-body theory of fish locomotion[END_REF]). This reasonable approximation, resulting of the so-called elongated body theory is well suited for a body described with analytical definitions. However, for the sake of simplicity, this approach, which requires to extract 2D sections and to build a 2D mesh around them from a possible complex 3D body shape was not retained in this way in the context of a 3D CFD solver. Alternatively, a slightly similar approximation was done to avoid computing this series of 2D computations section by section : it consists here in the application of a unit acceleration corresponding to the same DOF with respect to the neutral line on all nodes of the beam simultaneously. The pressure field obtained within a single 3D computation is then concatenated towards the different nodes of the beam with the same loads transfer used during the FSI computation. The resulting nodal forces and moments which corresponds to the loads on the section surrounding the beam node, give access to the added-mass coefficients for each node. This 7 procedure can be applied for the three main DOF of interest and for the longitudinal translation along the neutral line too if it is needed. As a result, only a maximum of four resolution of Poisson operator is required to deduce a sufficient approximation of the added-mass operator for the entire DOF of each flexible beam. It can be noticed that this approach removes the possible coupling between the different nodes, by assuming the acceleration between the nearest neighbours is similar. This greatly simplifies the computation of the relaxation operator [START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF], since I d + M -1 MA remains block diagonal.

Roll damping with active control of appendages

This application is dedicated to the study of the roll decay of a frigate with active control of fins (see figure 3). Passive fins are also studied as a reference point. The numerical results are compared to experimental tests. The roll motion of the ship and the relative rotations of the fins are taken into account through a sliding grid technique. A combination of rigid motion and deformation algorithms are also used to take over the other degrees of freedom. The mesh is made up of 31 millions cells. The scenario is divided into three parts. First, the ship velocity is imposed (others degrees of freedom are released) until a steady state is obtained. Then, the roll angle is imposed to perturb the ship in a close way compared with the experimental procedure. In figures 4 and 5, this ramp is applied between t 2 and t 3 . Finally, the roll motion is released and its decay is analysed. Comparisons between experimental tests and simulations are shown in figures 8 4, 5 and 6. One can see that both global quantities (roll angle) and local effects (drag and lift forces on fins) are well captured in spite of the difficulty to simulate this kind of complex experiments. It also can be seen in figure 4 that the roll damping is much more important when the fins are active than when they are passive.

fins bilge keels

Figure 3: DTMB with active fins -Geometry This application is a typical validation case for fluid structure interaction problems dealing with flexible beams. A flexible beam is attached to the upstream face of a fixed rigid square. Due to the vortex shedding created by the square, the flexible beam is exited and presents strong oscillations which are taken into account with a deformable mesh technique. The settings of this numerical two-dimensional case and others references can be found in [START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF]. The results obtained are shown in figure 7. They are between the classical values which can be found in literature. Visualisations of the velocity fields are also proposed in figure 8.

t 2 t 3 t 3+1 t 3+2 t 3+3 t 3+4 t 3+5 t 3+6
t 2 t 3 t 3+1 t 3+2 t 3+3 t 3+4 t 3+5 t 3+6 -1.0 -0.5 0.0 0.5 1.0 φ/φ 0 [-]
2 t 3 t 3+1 t 3+2 t 3+3 -1.0 -0.5 0.0 0.5 1.0 f y /B4t [-] (a) Starboard fin -Lift force t 2 t 3 t 3+1 t 3+2 t 3+3 -0.2 -0.1 0.0 0.1 0.2 f x /B4t [-] (b) Starboard fin -Drag force

3D flexible barge

The current work is dedicated to the study of a very flexible barge in waves. The results could be compared to experimental tests [START_REF] Malenica | Hydroelastic response of a barge to impulsive and non-impulsive wave loads[END_REF] and numerical simulations using a simplified model [START_REF] Senjanović | Investigation of ship hydroelasticity[END_REF]. The deformations of the barge for the first simulation with ISIS-CFD and MBDyn can be seen at two different instants in figure 9. This paper describes an efficient and robust algorithm, able to handle configuration with severe added mass effects, while remaining efficient. The coupling iteration occurs during the non-linear iterations of the fluid solver because it is the most costly part in our applications. To tackle the destabilising added-mass effects, a relaxation technique of the structural kinematics is used. The latter can be viewed as a modification of the artificial added mass method and does not need any intrusive modification of the structure solver. The relaxation operator is naturally related to the artificial added mass operator. Efficiency of the coupling is optimal when the artificial added mass is close to the physical one. The computation of the approximated or exact added mass operator for rigid and beam-like bodies is carried out through an on-line resolution of an original pressure-like equation integrated in the fluid solver ( [START_REF] Söding | How to integrate free motions of solids in fluids[END_REF][START_REF] Yvin | Interaction fluide-structure pour des configurations multi-corps. Applications aux liaisons complexes, lois de commande d'actionneur et systèmes souples dans le domaine maritime[END_REF]). Compared to a non-FSI simulation (unsteady simulation with imposed motion), the number of non-linear iterations to converge the fluid part is similar, even when large added mass effects occur. The additional cost of the FSI cases is then reduced to the resolution of the structural part and to the mesh deformation technique. This computation chain consisting here of ISIS-CFD and MBDyn has been validated numerically for 6 DOF rigid bodies with strong added mass effects [START_REF] Yvin | Interaction fluide-structure pour des configurations multi-corps. Applications aux liaisons complexes, lois de commande d'actionneur et systèmes souples dans le domaine maritime[END_REF] as well as for flexible bodies. Among these studies, some of them are briefly described here and comparisons with experimental data or other numerical results are shown.

Figure 1 :

 1 Figure 1: coupling algorithm

Figure 2 :

 2 Figure 2: computation of the added-mass operator for an elongated body

  CO POSSIBLE LARGE ADDED M APPLICATIONS TO RIGID AN FLEXIBLE BODIES IN THE MC. YVIN * , A. LEROYER † , E. GUILMINEAUFebruary 18, 2015

Figure 4 :

 4 Figure 4: Roll angle -Passive fins (left) and active fins (right) :

  t

Figure 6 :

 6 Figure 6: Hydrodynamic forces on fins -Active fins :

f

  [s -1 ] y [m] y = 1.05 m f = 2.90 s -1

Figure 7 :Figure 8 :

 78 Figure 7: 2D flexible membrane -Vertical displacement of the beam end (left) and FFT (right)

Figure 9 :

 9 Figure 9: 3D flexible barge -Barge at two different times