Christophe Kassiotis
email: christophe.kassiotis@enpc.fr

Adnan Ibrahimbegovic

Rainer Niekamp

Hermann G Matthies

Nonlinear fluid-structure interaction problem. Part II: space discretization, implementation aspects, nested parallelization and application examples

Keywords: fluid-structure interaction, CTL implementation, finite element, finite volume

The main focus of the present article is the development of a general solution framework for coupled and/or interaction multi-physics problems based upon re-using existing codes into software products. In particular, we discuss how to build this software tool for the case of fluid-structure interaction problem, from finite element code FEAP for structural and finite volume code OpenFOAM for fluid mechanics. This is achieved by using the Component Template Library (CTL) to provide the coupling between the existing codes into a single software product. The present CTL code-coupling procedure accepts not only different discretization schemes, but different languages, with the solid component written in Fortran and fluid component written in C++. Moreover, the resulting CTL-based code also accepts the nested parallelization. The proposed coupling strategy is detailed for explicit and implicit fixed-point iteration solver presented in the Part I of this paper, referred to Direct Force-Motion Transfer/Block-Gauss-Seidel. However, the proposed code-coupling framework can easily accommodate other solution schemes. The selected application examples are chosen to confirm the capability of the code-coupling strategy to provide a quick development of advanced computational tools for demanding practical problems, such as 3D fluid models with free-surface flows interacting with structures.

Introduction

The main thrust of this paper is geared towards the quick developments of computational tools for currently interesting multi-physics problems, pertaining to coupling or interaction of two or more traditional scientific domains. Of main interest for this work is how to provide the most efficient development of the corresponding software for this kind of problems, by coupling the stand alone software products from each traditional domain. The latter are often programmed by experts of a particular domain, with a strong preference to a particular discretization technique, programming language, or other special choices that provide the optimal result. Such an optimal result can in general be produced only by the corresponding domain experts, software designer and programmer, with a deep understanding of both physical phenomena and of numerical analysis issues pertaining to this particular domain. Hardly anybody can gain such level of expertise in many different domains. Therefore, the current trends in engineering developments geared towards multi-physics problems (e.g. fluid-structure interaction, thermomechanical coupling [START_REF] Kassiotis | Multiscale in time and stability analysis of operator split solution procedure applied to thermomechanical problems[END_REF], mechanics and optimization [START_REF] Ibrahimbegovic | Shape optimization of two-phase inelastic material with microstructure[END_REF], mechanics and control [START_REF] Ibrahimbegovic | Optimal design and optimal control of elastic structures undergoing finite rotations[END_REF]) are rather difficult to deal with, not only in terms of theoretical formulations but also in terms of development of corresponding software tools. One way to avoid this difficulty, as proposed and discussed in this work, is by coupling of different software products, each produced by experts for a particular sub-problem.

This flexibility in problem formulation and tools development should not be penalized with a severe lack of solution efficiency. Namely, we are interested in development of software tools for solving the problems relevant to industrial applications, which implies using fine discretization and mesh with many d-o-f; the examples of this kind are the problems reaching the half a billion d-o-f from solid mechanics applications to bio-mechanics [START_REF] Adams | Ultrascalable implicit finite element analyses in solid mechanics with over a half billion degrees of freedom[END_REF] or the large CFD computations pertinent to meteorological prediction or climate change [START_REF] Fosdick | An introduction to high-performance scientific computing[END_REF]. For the lack of efficiency, most symbolic interpreter software products, such as Matlab or Octave, are not suitable for this class of problems. In fact, in order to ensure the required level of computational efficiency in the context of fluid-structure interaction [3, [START_REF] Bathe | A mesh adaptivity procedure for CFD and fluid-structure interactions[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Degroote | Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction[END_REF][START_REF] Deparis | Fluid-structure algorithms based on Steklov-Poincaré operators[END_REF][START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF][START_REF] Farhat | Two efficient staggered algorithms for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems[END_REF][START_REF] Charbel Farhat | Mixed explicit/implicit time integration of coupled aeroelastic problems: three-field formulation, geometric conservation and distributed solution[END_REF][START_REF] Felippa | Synthesis tools for structural dynamics and partitioned analysis of coupled systems[END_REF][START_REF] Felippa | Stabilization of staggered solution procedures for fluid-structure interaction analysis[END_REF][START_REF] Fernández | A Newton method using exact Jacobians for solving fluid-structure coupling[END_REF][START_REF] Förster | On the geometric conservation law in transient ow calculations on deforming domains[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] Heil | An efficient solver for the fully coupled solution of large-displacement fluid-structure interaction problems[END_REF][START_REF] Hübner | A monolithic approach to fluid-structure interaction using space-time finite elements[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] Legay | An eulerian-lagrangian method for fluid-structure interaction based on level sets[END_REF][START_REF] Küttler | A solution for the incompressibility dilemma in partitioned fluid-structure interaction with pure Dirichlet fluid domains[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[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] Hermann | Algorithms for strong coupling procedures[END_REF][START_REF] Mehl | An Eulerian approach for partitioned fluid-structure simulations on Cartesian grids[END_REF][START_REF] Perić | Modelling fluid-induced structural vibrations: reducing the structural risk for stormywinds[END_REF][START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems-Part II: energy transfer analysis and three-dimensional applications[END_REF][START_REF] Michael | Treatment of acoustic fluidstructure interaction by localized Lagrange multipliers and comparison to alternative interface-coupling methods[END_REF][START_REF] Tezduyar | Space-time finite element techniques for computation of fluid-structure interactions[END_REF][START_REF] Tezduyar | Interface projection techniques for fluid-structure interaction modeling with moving-mesh methods[END_REF][START_REF] Tezduyar | Multiscale sequentially-coupled arterial FSI technique[END_REF][START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF][START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Wang | Fluid-structure interaction by the discontinuous-Galerkin method for large deformations[END_REF], large size problems are often tackled using dedicated software developments [START_REF] Tezduyar | Parallel finite element simulation of 3D incompressible flows: fluid-structure interactions[END_REF][START_REF] Tezduyar | Space-time finite element techniques for computation of fluid-structure interactions[END_REF][START_REF] Tezduyar | Modelling of fluid-structure interactions with the space-time finite elements: solution techniques[END_REF][START_REF] Tezduyar | Interface projection techniques for fluid-structure interaction modeling with moving-mesh methods[END_REF]. In this work, the large problems from fluid-structure interactionapplications are solved by using existing fluid and structure solvers, along with the corresponding partitioned strategy. The nested parallel computations, exploiting parallelization for both the fluid flow and interface interaction, are used to ensure the computational efficiency.

Granted sufficient efficiency, the strategy of re-using the stand-alone software for particular sub-problems is likely to become the most efficient way in development of software products for multi-physics applications. Indeed, in order to produce a reliable new software (and a bug-free code), the testing and validation is often the most time consuming phase that concerns not only the programmers and code developers but also the end-users. Namely, only after extensive use of a particular software product (by end-users), can we count with sufficient reliability of a particular software product. It is tacitly assumed that a reliable software product will keep its reliability in a more general framework of multi-physics problems.

Therefore, the proposed strategy of re-using software products (or components) for solving multi-physics problems, is currently becoming one of the major trends in scientific computing. The component-coupling framework described in this work, which allows to re-use existing codes for the fluid-structure interaction, is based upon the middleware CTL (Communication Template Library, see [START_REF] Niekamp | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior: Part IIsoftware coupling implementation aspects[END_REF][START_REF] Damijan Markovič | Multi-scale modeling of heterogeneous structures with inelastic constitutive behavior : Mathematical and physical aspects[END_REF]) handling communication between components. The main novelty for CTL implementation presented in this work concerns the use of nested parallelization for fluid-structure interactioncomputations, where the parallel communication is mastered by the CTL, and the call to a parallel component is transparent for the client.

The outline of the paper is as follows. In the next section we give a description of the space discretization methods used for structure and fluid; we combine herein geometrically nonlinear models for structure discretized by Finite Elements [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF] and a free-surface fluid flow formulation (based upon the Volume-Of-Fluid strategy [START_REF] Behzadi | Modelling of dispersed bubble and droplet flow at high phase fractions[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]) discretized by the Finite Volume [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. In Section 3 we describe how to provide the compatible space interpolations between two different discretizations on fluid-structure interface, and thus enforce the corresponding interaction. It is also shown how the chosen method to ensure compatibility of fields between the structure and fluid solvers can be implemented as a separate CTL-component that can be coupled with a number of different codes (see also [START_REF] Jürgens | Survey on software engineering for scientific applications[END_REF]). The performance of parallel computations with fluid component and the CTL-based nested parallelization of fluid-structure interactionproblems is studied in Section 4. The numerical simulations with a couple of large three-dimensional numerical examples, including large number of d-o-f and application to free surface flows, are given in Section 4. The concluding remarks are stated the last section.

Structure and fluid formulations

We start by giving the field discretization details typically used for solids and fluids, which are important to understand.

Structure equations of motion

The structure motion is based on the Lagrangian description. Namely, we consider a structure domain Ω s with imposed displacements u on the Dirichlet boundary ∂ Ω s,D and moving under the loading of traction forces t on the Neumann boundary ∂ Ω s,N and a body force b applied in the whole domain Ω s . This dynamic motion has to be computed in the time interval [0, T ]. The Cauchy governing equation for the structure describes the momentum conservation. The strong form of this equation can be written with respect to the deformed configuration 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: ∇ • Jσ F -T P +ρ 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. We indicated above that the Cauchy stress tensor σ in the deformed configuration can be linked to the first Piola-Kirchhoff stress tensor P formulated in the initial configuration through the gradient F of the deformation and its Jacobian J (see [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF]).

To close this Partial Differential Equations system we will link the displacements (or rather its derivatives) with the stresses through constitutive law. For instance, an elastic material model based on St.-Venant-Kirchhoff constitutive equation is assumed that links the Cauchy stress tensor σ and the Green-Lagrange strain tensor E through:

F -1 Jσ F -T = C : E, E = 1 2 F T F -I and F = I + ∇u (2)
where C denotes the constitutive fourth-order elasticity tensor 1 . Thus, it is a priori impossible to find directly an exact solution to the problem defined above. The idea is to find the best approximation of the solution in a finite-dimensional space where the solution can be found numerically. The FEM approximation [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF][START_REF] Olgierd | 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 weak forms of the equilibrium equation ( 1) and can be written for this problem as 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 sub-domains or elements T h = (κ e ) e=1,...,n el so that the whole space is covered with the finite elements that do not intersect. The solution space associated with this approximation solution is restrained to the space of continuous element-wise polynomial functions, which is denoted:

U h = U ∩ u ∈ C 0 (Ω s ) u| κ ∈ P p (κ) , ∀κ ∈ T h ( 3
)
where P p (κ) is the space of polynomials of order p on κ. The same restriction U h 0 holds on the associated vector space. The semi-discretized FE problem is defined as given t on ∂ Ω s,N and b in Ω s , find u ∈ U h such that, for all δ u ∈ U h 0 :

G s (u; δ u) = 0 (4)
This semi-discrete problem can also be written in a matrix notation by 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 ( 5
)
where M is the mass matrix, f int s with a geometrically nonlinear problem, and f ext s the consistent nodal forces. Here the λ represents the boundary forces computed from the fluid flow problem and imposed on the fluid-structure interface. Each matrix and vector of this semi-discrete equation are properly defined by assembling locally computed array in each element with the polynomial basis N e of P (κ e ):

In order to complete the discretization process, the time integration of the structure problem can be carried out by using 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]. In particular, 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 N max = T . In a typical time step ∆t = t N+1t N , the time derivatives of nodal displacement 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 (7)
In the semi-discrete form of the solid equation of motion in Eq. ( 5), 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 + ρ ∞ (8)
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].

Fluid equations in a moving domain taking into account free-surface

We use here 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] Ferziger | Computational Methods for Fluid Dynamics[END_REF]). In this Volume-of-Fluid (V.O.F.) method, an indicator function (volume fraction, level set or phasefield) is used to represent the interface; the main remaining issue is how to convect the interface without diffusing, dispersing or wrinkling it.This problem has been addressed in [START_REF] Behzadi | Modelling of dispersed bubble and droplet flow at high phase fractions[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]: Volume of Fluid (V.O.F.) methods use precise convection schemes that reconstruct the interface from the volume fraction distribution before advecting it.

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 ] (9)
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 ] (10)
where the characteristic function or fluid volume fraction ι is defined as:

ι(x,t) = 1, for x ∈ Ω 1 (t) 0, for x ∈ Ω 2 (t) (11)
Let us note that when the whole domain is filled with one fluid (ι = 1 ∈ Ω f ), the classical Navier-Stokes equations in ALE framework appear2 . The fluid volume fraction ι and the mass distribution δ Γ are linked by the following relation:

∇ι = δ Γ n (12)
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 ( 13
)
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 ι are the solution of equations ( 9) and ( 10) are verified. Traditional Computational Fluid Dynamic programs solve the fluid equations on a fixed (Eulerian) grid. This present a difficulty in fluid-structure interaction problems because of moving domains at the fluid interface which follows the deformations of the structure. The classical approach to overcome this difficulty is to consider the socalled Arbitrary Lagrangian Eulerian (ALE) method where the whole grid is moved inside the fluid domain, in trying to follow the movement of the boundary. However, this leads to the difficulty of pertaining the quality and the validity of the fluid inner mesh for different new shapes of the boundary This is solved by building a suitable map for the domain motion given its interface displacement. The fluid displacement u m is arbitrary inside the domain Ω f , but on the boundary it has to fulfill the condition:

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

u m = Ext u m ∂ Ω f ( 15
)
The latter is possible to construct by solving the Laplace smoothing equation:

∇ • (γ∇u m ) = 0 in Ω f ( 16
)
This kind of Laplace smoothing equation is known to have some limitation when the deformation of the fluid domain is governed by large rotations (e.g. see []). In our case it shows sufficient quality performance when the diffusion coefficient is made spatially dependent upon the distance to the fluid-structure interface.

The second major difference from solids concerns the favorite discretization technique for fluid in terms of the Finite Volume Method. The latter will transform the weak form of the continuous equations in (9) to ( 16) into a set of algebraic equations that can be solved numerically. One of the goal of this work is to couple different codes, with different dicretization methods (FEM for solids and FVM for fluids). this keeping the main trends of computational scientific software. Namely, 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, Open-FOAM [START_REF] Jasak | OpenFOAM: Open source CFD in research and industry[END_REF] are all FVM based). Moreover, there is a great number of reference books on the subject, among them [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF][START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. We are however aware that work is made in order to solve fluid with Finite Elements (see for instance [START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF][START_REF] Franca | Stabilized finite element methods[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] Legay | An eulerian-lagrangian method for fluid-structure interaction based on level sets[END_REF][START_REF] Hübner | A monolithic approach to fluid-structure interaction using space-time finite elements[END_REF][START_REF] Perić | Modelling fluid-induced structural vibrations: reducing the structural risk for stormywinds[END_REF][START_REF] Tezduyar | Space-time finite element techniques for computation of fluid-structure interactions[END_REF][START_REF] Tezduyar | Interface projection techniques for fluid-structure interaction modeling with moving-mesh methods[END_REF][START_REF] Takizawa | Wall shear stress calculations in spacetime finite element computation of arterial fluidstructure interactions[END_REF][START_REF] Wang | Fluid-structure interaction by the discontinuous-Galerkin method for large deformations[END_REF]).

The FV formulation can be written directly using an integrated form of the conservation equation written in [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. Another possibility is to consider the restriction of the solution space of weak form problems. For consistency with the structure sub-problem formulation we chose to describe the FV strategy in latter framework. The weak form of the Navier-Stokes equation can be written as [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 (17)
where U , I , V and P are suitable functional spaces for for the solution fields.

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 ) with- out overlapping (∩ n el e=1 κ f ,e = / 0
). 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 ) (18)
where ι κ is the characteristic function of the element defined as:

ι κ : Ω f -→ | x -→ 1 if x ∈ κ 0 if x ∈ Ω /κ (19)
The same kind of restriction holds for the solution spaces of equation [START_REF] Deparis | Fluid-structure 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 ( 17) 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 need not to be 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 semi-discrete form of the discretized fluid problem can be written in a matrix form as follows. The fluid mesh motion considers that u m is imposed by the motion of the interface

u: R m (u m ; u) := K m u m -D m u = 0 ( 20
)
where D m is a projection/restriction operator and K m governs the extension of the boundary displacement (see Eq. ( 16)). The volume fraction ι, the d components of velocity v and pressure p are coupled through 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 (ι)v f + N f (ι, v f -um )v + K f (ι)v f + B f p f -f f (ι) B T f v f   = 0 (21)
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 unsymmetric advection matrix, K f is the conduction matrix describing the diffusion terms, and B f is 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].

One way to solve the flow problem is to consider a monolithic solver handling all equations simultaneously. Another way is to consider a split between the mesh motion, the volume fraction advection, the momentum and the continuity equations, and to use an operator split-like procedure often referred to 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 scheme. 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 segregated approach will be used because of its efficiency.

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 ensured explicitly.

The only remaining question is the choice of velocity variation in the time step ∆t = t N+1 -t 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 proved [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 integration schemes (like implicit Euler or Crank-Nicholson):

um = u m,N+1 -u m,N ∆t (22)
The volume fraction function is supposed to be sharp at the interface between water and gazs, and therefore, standard FV discretization that can be strongly diffusive cannot be applied, as they would smear the interface. One way to guarantee a sharp and bounded solution is by using 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 requires time sub-cycles, since an explicit treatment referred to as MULES (Multidimensional Universal Limiter with Explicit Solution) is used (see [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 equlibrium, 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] Sv Patankar | A calculation procedure for heat, mass and momentum transfer in threedimensional 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 pressurecorrection techniques is given and in [4], a comparison of two pressure-correction algorithms showed the overall better performances of PISO algorithms over SIMPLE ones.

In this work, PISO algorithms are used to solve our CFD problem. The semi-discrete form of the Navier-Stokes equation is discretized in time using implicit integration schemes (the Euler implicit, or the second order Crank-Nicholson scheme). The discretized momentum Eq. ( 21) is split into 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 ( 23
)
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 (24)
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. It 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 (meaning that the transport due to the flow over a cell is not well captured) the PISO algorithm fails to converge.

Space interpolation between solvers

The use of different discretization and time-integration schemes, for the fluid and the structure part, do not provide in general a matching mesh at the interface. This feature is taken into account in the nonlinear stability analysis performed in the Part I of this paper by interpolation matrix in the contraction parameter α. 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 example 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 sub-problems 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]. For the sake of simplicity, the use of Dirac functions has often been proposed [START_REF] Ibrahimbegovic | Strong coupling methods in multi-phase and multi-scale modeling of inelastic behavior of heterogeneous structures[END_REF][START_REF] Hautefeuille | Numerical Modeling strategy for Heterogeneous Materials: A FE Multi-scale and Component-based Approach[END_REF].

In our approach, it was decided not to favor any particular mesh-based representation at the interface, and even allow that the fluid problem can also be solved by a meshfreebased method [START_REF] Jj Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Dalrymple | Numerical modeling of water waves with the SPH method[END_REF]. Namely, an interpolation strategy relying on radial basis function is chosen herein. As suggested in [START_REF] Beckert | Multivariate interpolation for fluid-structure-interaction problems using radial basis functions[END_REF], this method fits well with the current work where, both sides of the problem are based on different approximation techniques. First, for the structure we employ the FE approximation that requires values imposed at the nodes and gives in return nodal values as answers. Second, in the FV discretization for the fluid part, the fields at the interface are defined at the center of the cells, and not at the nodes as it is done for FE. However, as the client has no reason to know the connectivity table of the coupled component nor the faces centers coordinates, it was decided to interpolate the fields internally in ofoam-2 and to set or get them at the points of the interface.

The Steklov-Poincaré operators for the fluid and the structure (described in Part I of this paper and based on the FE and FV solvers) request imposed boundary values at the nodes placed interface, and provide the values, for potentially different nodal points x f and x s . The interpolation step from fluid to structure consist in solving a linear system as explained subsequently.

Let us consider that the displacement vector at the N s solid nodes (x s,i ) i=1...N s is given as (u s,i ) i=1...N s . We want to build the interpolation of this displacement field at the N f fluid nodes in order to impose the fluid mesh motion as needed in ALE fluid computations. An interpolation for each scalar field of the projected displacements in one of the Cartesian directions has then to be built. This scalar field is denoted as u i = u(x s,i ) at each points, and the interpolation has to be performed d times for each vector field of dimension d.

Then the interpolation at a point x takes the form:

u(x) = N s ∑ i=1 c i Φ(x -x s,i ) (25)
where Φ is a fixed basis function which is radial with respect to the Euclidian distance, i.e.:

Φ(x) = Φ( x 2 ) = Φ( d ∑ i=1 x 2 i ) (26)
The coefficients c i are determined by the interpolation condition:

u j = N s ∑ i=1 c i Φ(x s, j -x s,i ) (27)
Thus, for interpolation of any field to interpolate, the coefficients c i have first to be computed as the inverse of the following relation:

M i j c j = u j ( 28
)
where the interpolation matrix entries are computed as:

M i j = Φ(x s, j -x s,i )
The last step is to compute the interpolated field by using the expression given in Eq. ( 25).

The choice of the radial basis function is governed by the requirement that the influence of a center has to be smaller when the distance to an evaluation node increases. Elsewhere, the global character of the radial basis function tends to smooth out all local effects. Another desirable property is to use the functions with compact support which results with a sparse interpolation matrix. In [START_REF] Beckert | Multivariate interpolation for fluid-structure-interaction problems using radial basis functions[END_REF], a comparison between different smooth compact radial basis functions for field transfer at the interface is given. In this work, global radial basis function are used in the exponential form:

Φ(x) = exp - x h ( 29
)
where h is a characteristic distance between points. The Interpolator implementation is described in [START_REF] Jürgens | Survey on software engineering for scientific applications[END_REF]. The advantage of the component use is the possibility to plug-in any interpolation software to the fluid-structure interaction framework that matches the Component Interface. The inverse problem of the full matrix M i j is here relying on the lapack or the Intel Math Kernel Library; the cost of such an operation is around N 3 , but as the interface d-o-f are by far less than number of d-o-f for each subproblem, in general this cost cannot be the bottleneck of the computations. One of the advantages of the proposed approach is that the interpolation is exact for matching nodes. It can be used without any lack of precision in a generic way for matching meshes where the ordering is not the same for both subproblems.

4 COupling the software COmponents by a Partitioned strategy (cops) and parallel computations 4.1 Software coupling applied to fluid-structure interaction Coupling a fluid code and a structure code in order to produce the code for solving the fluid-structure interaction problem is carried by master coupling approach, where a master code sends request and receive data from the coupled software. Traditionally this approach is quite intrusive, and imposes new developments inside the software. However, the use of component technology with the middleware CTL (see [START_REF] Niekamp | Multi-scale modelling of heterogeneous structures with inelastic constitutive behavior: Part IIsoftware coupling implementation aspects[END_REF]) eliminates the difficulties in the master coupling approach. We will here insist on two points: the development of a component based upon an existing code requires a very good understanding of its architecture, and a deep knowledge of the physics behind. However, once the component is developed, it can actually be used as a black box.

The COupling COmponents by a Partitioned Strategy (cops) provides a generic implementation of the explicit and implicit DFMT-BGS algorithms for fluid-structure interaction. The key idea is to re-use existing codes, here FEAP for structures and OpenFOAMfor fluids, as elementary bricks to build upon the CTL components.

We note that cops itself is a component and has been successfully re-use in [2] for sequential multi-scale approach to solving fluid-structure interaction problem. Dependency graph of cops component is presented by a general overview of its architecture given in UML syntax in Fig. 1.

The interface, defined in a .ci file, is realized by a C++ class named cops. The main routine calls the coupling component, providing the path where the control file associated to cops lies. The coupling component cops instantiates two sub-solvers, the one of coFeap for the structure and another of ofoam-2 for the fluid problem and with a partitionedSolver binding the two subproblems. PartitionedSolver is a template class that couples any component matching a SimuCI or a CFDSimuCI where method is necessary to build the Steklov-Poincaré operators, and optionally goback in time if sub-iterations of one subproblem in implicit computations is asked for.

The PartitionedSolver instantiates a Picard method that is in charge of one Picard iteration upon the coupled subproblem. One Picard iteration allows to determine the new residual needed to check the convergence of the DFMT-BGS solver. The Picard class itself relies on a template SteklovPoincare class.

Performance Comparison between system calls and file reading/writing for the fluid component

In its first implementation (see [START_REF] Krosche | Ofoam's manual. Informatikbericht[END_REF]), the ofoam component architecture was based on file reading and system calls. Behind the CFDSimu.ci interface is implemented a class that reads output files from OpenFOAM and uses system calls in order to perform computation. This implementation has the advantage of being not intrusive at all, since the component depends neither on header files nor on OpenFOAM libraries. However, the slow execution speed for this kind of implementation makes it prohibitively expensive for any large scale computations.

The second implementation described herein above concerns a component directly linked to OpenFOAM, and it has the major advantage of working than the first. In order to evaluate the performance gain obtained with the new implementation, we compare the CPU time required by the call of methods of each component (ofoam-1 -based on file reading -and ofoam-2 -the wrapper class and component build around OpenFOAM) on the same test case.

The latter concerns a Newtonian flow in a cylinder for fully 3D case. The chosen dimensions of the cylinder: diameter 1m and length 10m. The velocity on the inlet is imposed as being uniform (10m.s -1 ), whereas at the outlet, the pressure gradient is set to 0. Perfect wall conditions (zero velocity) are applied to tube walls. We consider following values for material properties: density ρ = 1.0kg.m -3 and dynamic viscosity µ = 0.01m 2 .s -1 . The numerical simulations spans time interval from t = 0s to t = 0.01s. The discretization in space is carried out by second order Finite Volume approximation. For time discretization, the implicit Euler scheme with dt = 0.001s is applied. The algorithm chosen to handle the computations with incompressibility constraint is based on PISO [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. At each time step, two outer corrections are performed to ensure the pressure velocity coupling. The solvers for the fluid velocity and pressure fields are based on PBiCG and PCG respectively.

The computations are performed for meshes with 9 × 10 3 , 72 × 10 3 and 576 × 10 3 cells (the number of d-o-f equals roughly 4 times the number of cells). The velocity field magni- 1 Performance comparison between ofoam and ofoam-2 in terms of CPU time for each method given in seconds; the latter has been measured using the CTL profiler (CTL Profile).

The overall performance of ofoam-1 and ofoam-2 components as expected (see Table 1). The ofoam-1 component based on file reading will be intrinsically less efficient for methods that require intrinsically a lot of data exchange compared to the computation cost (methods get nodes, pressure and velocity and set velocity and mesh displacements), since the speed of access to data by opening and reading files written on the hard drive cannot compete with a direct access to the memory. However, for all the methods that require a lot of computations, the performances of ofoam-1 and ofoam-2 implementations are comparable and mainly depends on the implementation of OpenFOAM itself, as the time spent in data exchange is very small compared to the one spent in the computation, especially for a huge number of d-o-f. In Fig. 3, the bottleneck of a fluid computation pertains, as expected, to the problem solving (matrix inversion for any method used). Note that the solve method is a bit faster for ofoam-2 implementation since there is no need for the software to reload data before performing computation. Thus, the slight advantages offered in terms of CPU time does not alone justify the choice of building a new component ofoam-2. Rather the possibility of reusing any class developed in the OpenFOAM project provides a better argument in favor of the new implementation

Parallel CFD Component Performance and nested parallelization

For most fluid-structure interaction problems, the bottleneck of performances are flow computations. For instance, in the problem presented in this article, where the fluid domain surrounds the structure, more than 95% of the computational time is spend in the fluid computation. Therefore, even if the parallelization of the coupling algorithm is of interest, this can result in sufficient increase of efficiency only in the case if the fluid and the structure computations require the same CPU time. For all other cases we will used the nested parallelization, where the parallel coupling algorithm can accommodate parallel fluid computation. The implementation of nested parallel coupling algorithms can then be accomplished with CTL in a quite straightforward manner by using a certain number of instance of the ofoam-2 component controlled by the CTL.

The performance of the parallel version of the ofoam-2 component is tested on the same problem of flow in a cylinder. For space discretization and the chosen mesh is handled by METIS [START_REF] Karypis | METIS, A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices[END_REF] in each sub-domain with the same weight (Fig. 4). In each computation, the bottleneck is computation of the flow evolution between two time steps. The computational time reported in this section is the average time T CPU required to solve one time step of the given problem and to write the associated results. The initialization time is therefore not taken into account.

We consider the main indicator to measure the performance of the nested parallelization to be the speed-up χ defined as:

χ = T CPU 1 T CPU N ( 30
)
where T CPU 1 is the computational time for the problem solved on one processor and T CPU N for N processors. A linear speed-up: it means that taking N processors, one expects to divide the computational time by N, cannot be obtained with ofoam-2. Namely, the communication between processes, that is highly linked to size of interfaces, has to be taken into account. For this reason, the notion of efficiency can be introduced in hoping to have it close to 1:

ξ = T CPU 1 N × T CPU N ( 31
)
We tested the nested parallel computations on the same multi-processor machine or on a cluster architecture where computers communicate throughout a network, and we report both set of results. We underline that the given results are empha priori different from the one obtained comparing the parallelization of OpenFOAM, as the communication is made between CTL components.

Parallelization on the same multi-processor machine

For this case we consider a rather coarse mesh with 72 × 10 3 cells (≃ 216 × 10 3 d-o-f). Computations are run with the same parameters as the one used for the comparison between ofoam-1 and ofoam-2 (see Sec. 4.2). The CPU time results are given in Table 2.

number of method solve method init processor N T CPU N in (s) T CPU N in (s) 1 5
Processor number (N) Fig. 5 Parallel speed-up on a multi-processor machine

The graphic illustration of the results is given in Fig. 6, showing that efficiency of nested parallel computing is maintained around 0.7, even with 7 processors. Similarly, in parallel computation we 7 processors one can observe a Speed-up (Fig. 5) of around 5 times. Processor number (N)

Fig. 6 Parallel efficiency on a multi-processor machine However, the question of the representative value of this computation can be raised, as the mesh size is maybe not big enough and compared to the time spend in the solving, the communication cost between processes may not be negligible. For that reason, we compute the flow in a cylinder problem for different mesh refinements -from 27 × 10 3 to roughly 18 × 10 6 d-o-f -on one processor, and then repeat parallel computations on 6 processors. Here is represented the computational time spent for the solve method with direct and compressed transfers between the processes. In compressed transfer, each double is converted into a float, leading to a small loss of accuracy.

These results are illustrated in Fig. 7 in terms of efficiency of the parallel computation on 6 processors as a function of number of cells for different mesh refinement. One observes a decrease in the performance for coarser grids, as for these cases the communication between processes is not negligible. This is even more noticed for standard transfer when one chooses not to compress the double into float.

Moreover, we also observe a decrease in the performance for finer meshes. This is the consequence of imposing the incompressibility conditions, which requires more iterations in a parallel than in a serial run. Namely, it is known that state solver convergence is affected by running in parallel as the preconditioning and smoothing operations in the cells adjacent to the boundaries are less effective, and lead to a small increase in the number of iterations. Moreover, in any such case, decreasing the accuracy when the values between processes are transmitted with compact transfer is not a good choice, since the gain obtained on the com- munication process (especially on the same machine) does not justify the loss of precision and result instead in an increase of the number of iterations. Finally, for the finest meshes, the performance of the parallel computation drops significantly. This phenomenon seems to be due to the architecture of one cluster node where not only processor velocity but also the bus communication speed for exchanging data between processors limit the overall performance, especially when huge transfer of data is required. This difficulty can be overcome with parallel computations on a cluster with several machines, as shown in the next section.

Parallelization on a cluster and network communication

For the parallel computations on several nodes (machines) of a cluster, we only consider the finer grid with 4.608 × 10 6 cells (roughly 18 × 10 6 d-o-f). Each node is linked with the other via a fast network, and possesses 8 Intel processors with a 3.0GHz frequency. Computations are run on eleven nodes of the cluster. The load of each node is maintained at the same level as much as possible3 . The results obtained for the solve method are presented in Table 4. The time spent for other methods like initialization is not presented, as the result in Table 2 has shown that the parallelization is of little influence for these methods. Fig. 8 represents the real speed-up observed to solve one time-step against the theoretical linear speed-up. The latter is shown to confirm that almost linear speed-up is observed until more than 4 processors of each node are loaded. The same value was observed in the previous section, when the parallelization was done on the processors of the same machine, thus indicating some limitation of the architecture of the cluster nodes. The efficiency of the parallel computing, as shown in Fig. 9 remains around optimum value of 1 before it starts decreasing. Note that the efficiency for parallel computation in a small number of cluster nodes is above the theoretical optimum. This is due to the fact that less data need to be handled by the memory on each nodes. It is also of interest to compute the CPU time for the same range of mesh size as the one given in Table 5. We also computed the efficiency obtained for different mesh sizes and present the results in Fig. 10. For this computation, each grid is split into 6 sub-domains, then solved on 6 different machines. Contrary to the phenomenon observed for parallel runs on the same machine (Fig. 7), the expected gain in efficiency is observed with increasing number of d-of and parallel computation on different machines. In other words, the bigger is the mesh, the more efficient is the parallel computation, since the communication between components then remains small compared to the time spent for solving the problem.

Handling data transfer with compacting the exchanged values from double to float for is not advantageous for all computations of the studied flow case. Namely, the speed-up in communication for compact transfer is diminished by the loss in accuracy that leads to more iterations to smooth the values and precondition the solvers, especially near the sub-domain boundaries. The fluid problem is discretized with FV method, and then split in 6 sub-domains by METIS software tool in order to performs parallel computations (see Fig. 12). The chosen boundary conditions, specified in Fig. 11. are as follows: for lateral walls of the fluid domain, the velocity boundary condition allows for slipping; at the inflow, a constant velocity is imposed with v = (100cm.s -1 , 0, 0); Zero gradient pressure is specified at the outflow.

In order to smooth the first steps of the fluid-structure interaction computations, we do not start from the rest, but rather perform the fluid only computation from t = -2s to t = 0s, with the inflow velocity increasing in a smooth way the velocity according to: 1 2 sin π t+1 2 + 1 The fluid-structure interaction will start at t = 0s. All the points of the fluid mesh will move in the ALE strategy, with their motion governed by a smoothing process based on the Laplacian operator and a diffusivity coefficient whose value depends on the distance to the appendix. The fluid discretization techniques and solvers are equivalent to the one used in the two-dimensional example.

In the previous computation of this problem achieved in [START_REF] Malte Von Scheven | Effiziente Algorithmen für die Fluid-Struktur-Wechselwirkung[END_REF], the discretization of the structure problem is performed by shell finite elements. Three-dimensional elements with quadratic shape functions are used herein, with each elements therefore containing 27 nodes. -Venant-Kirchhoff constitutive model (see Eq. 2 and [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF]). The time integration is handled by a generalized-α scheme with: ρ ∞ = 1 2 ; β = 4 9 ; γ = 5 3 ; and α = 2 3 . At each iteration, the linearized system of solid equations of motion is solved by a direct solver for real value non-symmetric system matrix. The total number of d-o-f for the coupled problem is given in Table 6 for both coarse and fine mesh. The computation is carried out with a coupling time step of 1 × 10 -3 s for each mesh. The coupling scheme used is DFMT-BGS with Aitken's relaxation. The initial relaxation parameter is ω = 1.0. The chosen value of the residual tolerance for coupled simulation is: r The displacement of the points at the free-end represented in Fig. 14, show that the motion of the flag corresponds to the first flexural mode. It is hard to predict the exact solution for such a complex three-dimensions flow in with a relatively high Reynolds number. It is interesting to note that the results in [START_REF] Malte Von Scheven | Effiziente Algorithmen für die Fluid-Struktur-Wechselwirkung[END_REF] indicate that, after a certain time, some torsion modes occur in the flag, that could not be confirmed here at the same Reynolds number. Consequently, the motion amplitude presented herein is around 4 times bigger than the one obtained in [START_REF] Malte Von Scheven | Effiziente Algorithmen für die Fluid-Struktur-Wechselwirkung[END_REF]. Our results are more in agreement with the ones provided for 2D version of the same problem (see [START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF][START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF][START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF]) with a flexible appendix.

Three-dimensional sloshing wave impacting a flexible structure

In [START_REF] Kassiotis | Partitioned solution to fluidstructure interaction problems in application to free-surface flows[END_REF], we have studied and validated our coupling framework in application free-surface flow in interaction with non-linear structures. The problem presented herein emphasize the 3D capacities of our approach as well as the possibility to deal with complex flow. This example is a simplified representation of a dam-breaking event that brings about the sloshing wave impact on a flexible structure standing in its way, as presented in Fig. 15. At initial time t = 0s, a three-dimensional water column starts falling down under the gravity loading. Spreading further at a later time, it hits the obstacle which is a slender plate-like solid body made of elastic material that can undergo large deformations. The dimension of the problem and the imposed boundary conditions are given in Fig. 15.

In order to prevent that the water will bounce-back and again hit the structure after breaking on the walls, only the left and bottom planes of the fluid domain are defined as non-slipping walls, while the others are defined with atmospheric boundary condition for the pressure.

The material properties are chosen as follows: for the high density fluid (the water) the density and the kinematic viscosity are ρ f ,1 = 1 × 10 3 kg.m -3 and ν f ,1 = 1 × 10 5 m.s -1 , indicator equation, as well as the half-implicit nature of the coupling between the momentum predictor and the pressure corrector. At this scale of modelling it is not required to consider surface tension between the two fluids.

The structure model is constructed by using three-dimensional elements with quadratic shape functions, where each element has 27 nodes. The material properties used herein correspond to a neo-Hookean elastic material with Young's modulus E s = 1 × 10 6 Pa, Poisson's ratio ν s = 0, and a density ρ s = 2500kg • m -3 , The chosen model an represent finite deformation. The time integration is carried out by a Generalized-α scheme with the chosen parameters as ρ ∞ = 1 2 β = 4 9 , γ = 5 3 and α = 2 3 . The total number of d-o-f given in Table 7 for both the coarse and fine discretization. The computation of the coupled problem is carried out by an iterative scheme. The results of fluid and structure 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 is DFMT-BGS with Aitken's relaxation with the initial parameter ω = 0.25 and the predictor of order 1. The absolute tolerance for coupled computation is equals to: The number of iterations required to reach the convergence criteria is given in Fig. 16. Note that there is no coupling iteration before the water hits the structure (the effect of air flow can almost be deemed negligible with respect to the structure). During the waterstructure contact the number of iteration depends on chosen the discretization density. In reaching the opposite wall, the water does not rebound on the wall but simply flows away, which again does not require any iteration.

r (k) N ≤ 1 × 10 -6 (32)
In Fig. 18, the water or 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 obtain at t = 0.25s, before the solid comes back to free-vibration phase. In Fig. 17 the motion of the free-end of the obstacle is plotted. Contrary to the twodimensional version of this example presented in [START_REF] Kassiotis | Partitioned solution to fluidstructure interaction problems in application to free-surface flows[END_REF], small drops of high density fluid are not interacting with the obstacle after the main shock. Therefore, the motion of the solid part remains rather well described even with the coarsest grid.

Conclusion

We propose in this work a solution approach to fluid-structure interactionproblems that allows coupling different space discretization methods, such as FE for structures and FV for fluids. As illustrated by numerical results, the proposed strategy is applicable to demanding problems of this kind that deal with complex free-surface flows interacting with geometrically non-linear structures. The proposed strategy can also employ different time integration schemes, and interface matching conditions between fluid and structure in the spirit of implicit analysis and chosen interface field representation. The stability in this nonlinear context is demonstrated by numerical examples, confirming the proof from mathematical analysis provided in Part I of this paper).

The key idea of this work pertains to producing the final software tools for complex FSI problems by coupling the existing software products that were developed (and tested) previously for either structure or fluid. The code-coupling strategy is capable of dealing with full 3D models, where the computational efficiency is of paramount interest. The latter is ensured by the nested parallelization, where the parallelization is carried out not only for coupling of fluid and structure through interface matching condition, but also for flow computations that are the most expensive task. Both code-coupling and nested parallelization are handled by the CTL, with the latter as its new feature.

  (u e ) = κ e ∇N : P (u e N e )

Fig. 1 A

 1 Fig. 1 A simplified UML class diagram for cops

Fig. 2 1 cells

 21 Fig. 2 Velocity field magnitude in the cylinder in m.s -1

Fig. 3

 3 Fig.3Performance comparison between some ofoam and ofoam-2 methods for different mesh refinements.

Fig. 4

 4 Fig. 4 Cylinder meshed with 9 × 10 3 cells split via METIS in 8 sub-domains for parallel runs

Fig. 7

 7 Fig. 7 Parallel computation performance for different mesh sizes

Fig. 8

 8 Fig. 8 Parallel computation speed-up for runs on different cluster nodes

Fig. 9

 9 Fig. 9 Parallel computation efficiency for runs on different cluster nodes

Fig. 10

 10 Fig. 10 Parallel computation performance on six different cluster nodes for increasing mesh sizes from 216 × 10 3 to 18 × 10 6 d-o-f

Fig. 11

 11 Fig. 11 Flag in the wind: geometry and boundary conditions

Fig. 13

 13 Fig.13Flag in the wind: motion of the structure and stream-tube snapshots for some time steps.

  Oscillating 3D Appendix: extremity displacements for points A = (10.0, 5.5, 3.0), B = (10.0, 5.5, 5.0) and C = (10.0, 5.5, 7.0). +SFEM f DFMT-BGS (b) Oscillating 3D Appendix: extremity displacements for points C comparison with[START_REF] Malte Von Scheven | Effiziente Algorithmen für die Fluid-Struktur-Wechselwirkung[END_REF] (DFMT-BGS coupling FEM for structure and stabilized FEM for fluid).

Fig. 14

 14 Fig. 14 Displacement of structure extremity

Fig. 16

 16 Fig.[START_REF] Demirdžić | Space conservation law in finite volume calculations of fluid flow[END_REF] Number of iterations in order to make the DFMT-BGS algorithm converge for the three-dimensional dam-breaking problem

Fig. 17

 17 Fig.17Three-dimensional dam break example: obstacle displacement measured at the center of the top face

Fig. 18

 18 Fig. 18 Tri-dimensional dam break problem. Evolution of the free surface and motion of the structure.

Table 3

 3 CPU Time for the solve for different meshes 3.0Ghz with standard and compact transfer

Table 4

 4 CPU Time for the solve method parallelized on 11 cluster nodes.

Table 5

 5 CPU Time for the solve for different meshes with standard and compact transfer

		number of	reference	standard transfer	compact transfer
		cells	T CPU 1	in (s)	T CPU 6	in (s)	T CPU 6	in (s)
		9 × 10 3	3.83 × 10 -1	2.74 × 10 -1	3.33 × 10 -1
		30 × 10 3	1.55 × 10 0	5.44 × 10 -1	6.38 × 10 -1
		72 × 10 3	4.38 × 10 0	1.16 × 10 0	1.29 × 10 0
		243 × 10 3	2.45 × 10 1	4.02 × 10 0	4.58 × 10 0
		576 × 10 3	8.11 × 10 1	1.28 × 10 1	1.57 × 10 1
		4608 × 10 3	1.55 × 10 3	2.59 × 10 2	2.90 × 10 2
		1	compact transfer standard transfer
	(ξ )	0.8
	Efficiency	0.4 0.6
		0.2
		7 • 10 3		7 • 10 4	Number of cells	7 • 10 5	7 • 10 6

Table 6

 6 Number of d-o-f for coarse and fine discretization of the three-dimensional oscillating appendix

Table 7

 7 Number of d-o-f for coarse and fine discretization of the three-dimensional dam-breaking problem

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])

Note that the stability proof for coupled fluid-structure interaction problem given in Part I is for this classical case

For instance, for the parallel computing on 32 processors, the first ten nodes were loaded on 3 of their processors whereas those on the last one two processors were used.

Fig. 15 Three-dimensional wave impacting an obstacle: geometry and boundary conditions whereas for the low density fluid (the air in the remaining part of the domain) ρ f ,2 = 1kg.m -3 and ν f ,2 = 1 × 10 6 m.s -1 . The mesh motion problem is solved using a Laplacian smoothing material where the diffusion coefficient is a quadratic inverse function of the distance to the interface between structure and fluid.

The fluid domain is discretized with Finite Volume cells covering always the complete domain either by one phase or by the other. The computations are performed for two different meshes with the chosen discretization and the number of cells given in Table 7. An explicit-implicit algorithm is used to compute the two phase flow evolution from the corresponding the Navier-Stokes equations. In the Volume of Fluid (V.O.F.) method used herein (see [START_REF] Ghidaglia | On the numerical solution to two-fluid models via a cell centered finite volume method[END_REF]), an indicator function (volume fraction, level set or phase-field) is used to represent the phases; leaving only the remaining issue on how to convect the interface without diffusing, dispersing or wrinkling it. This is particularly troublesome when the volume fraction is chosen as an indicator function because the convection scheme has to guarantee that the volume fraction stays bounded, with the values that remain within its physical bounds of 0 and 1. We here follow the idea proposed in [START_REF] Behzadi | Modelling of dispersed bubble and droplet flow at high phase fractions[END_REF]: Volume of Fluid (V.O.F.) method using convection schemes that reconstructs the interface from the volume fraction distribution before advecting it. The equation associated with the characteristic function is solved with an explicit time integration scheme whereas the remaining terms are solved with implicit time integration schemes (see [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]). The fluid is handled by second order space discretization with a Van Leer limiter used for the advection terms, and the implicit Euler time-integration scheme. Note that small time steps are required for the explicit solution of the phase function