Domenico Panunzio

Maria Adela Puscas

Romain Lagrange
email: romain.lagrange@cea.fr

FSI -Vibrations of immersed cylinders. Simulations with the engineering open-source code TrioCFD. Test cases and experimental comparisons

Keywords: Vibration, Fluid-structure interaction, Added mass, Added damping, Stokes number, ALE method, TrioCFD j ], [C j ] Added mass and damping matrices m

In this paper, we assess the capabilities of the Arbitrary Lagrangian-Eulerian method implemented in the open-source code TrioCFD to tackle down two fluid-structure interaction problems involving moving boundaries. To test the code, we first consider the bi-dimensional case of two coaxial cylinders moving in a viscous fluid. We show that the two fluid forces acting on the cylinders are in phase opposition, with amplitude and phase that only depend on the Stokes number, the dimensionless separation distance and the Keulegan-Carpenter number. Throughout a detailed parametric study, we show that the self (resp. cross) added mass and damping coefficients decrease (resp. increase) with the Stokes number and the separation distance. Our numerical results are in perfect agreement with the theoretical predictions of the literature, thereby validating the robustness of the ALE method implemented in TrioCFD. Then, we challenge the code by considering the case of a vibrating cylinder located in the central position of a square tube bundle. In parallel to the numerical investigations, we also present a new experimental setup for the measurement of the added coefficient, using the direct method introduced by Tanaka. The numerical predictions for the self-added coefficients are shown to be in very good agreement with a theoretical estimation used as a reference by engineers. A good agreement with the experimental results is also obtained for moderate and large Stokes numbers, whereas an important deviation due to parasitic frequencies in the experimental setup appears for low Stokes number. Still, this study clearly confirms that the ALE method implemented in TrioCFD is particularly efficient in solving fluid-structure interaction problems. As an open-source code, and given its ease of use and its flexibility, we believe that TrioCFD is thus perfectly adapted to engineers who need simple numerical tools to tackle down complex industrial problems.

Introduction

The accurate estimation of the force acting on a body moving in a viscous fluid is of crucial importance for engineers, in particular for those working in fields such as turbomachinery [START_REF] Furber | Is the Weis-Fogh principle exploitable in turbomachinery?[END_REF], heat exchangers tube banks [START_REF] Chen | Vibration of nuclear fuel bundles[END_REF][START_REF] Chen | Dynamics of heat exchanger tube banks[END_REF] or energy harvesting of flexible structures [START_REF] Doare | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF][START_REF] Singh | Energy harvesting from axial fluid-elastic instabilities of a cylinder[END_REF][START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF][START_REF] Virot | Coupling between a flag and a spring-mass oscillator[END_REF][START_REF] Eloy | Aeroelastic instability of cantilevered flexible plates in uniform flow[END_REF]. When the immersed body is subjected to a small amplitude of motion, Stokes [START_REF] Stokes | On the effect of the internal friction of fluids on pendulums[END_REF] showed that the fluid force is the sum of two terms: an added mass term related to the body's acceleration and a damping term related to the body's velocity. The concept of added mass and damping terms can be generalized to multiple immersed bodies, introducing some self and cross-added coefficients. In such a case, the self-added coefficients relate the fluid force on a body to its motion. On their hand, the cross-added coefficients relate the forces on a body to the other bodies' motion. The determination of the added coefficients has been the topic of considerable experimental, [START_REF] Clough | Effects of eathquakes on underwater structures[END_REF][START_REF] Chandrasekaran | Virtual mass of submerged structures[END_REF][START_REF] Stelson | Virtual mass and acceleration in fluids[END_REF][START_REF] Keulegan | Forces on cylinders and plates in an oscillating fluid[END_REF][START_REF] Garrison | Hydrodynamic loads induced by earthquakes[END_REF][START_REF] Sarpkaya | Separated flow about lifting bodies and impulsive flow about cylinders[END_REF][START_REF] Chandrasekaran | Vibration of submerged structures[END_REF][START_REF] Skop | Added Mass and Damping Forces on Circular Cylinders[END_REF][START_REF] Sarpkaya | Forces on cylinders and spheres in a sinusoidally oscillating fluid[END_REF] and theoretical studies, mostly based on a Helmholtz decomposition associated with a method of images [START_REF] Carpenter | On the motion of two cylinders in an ideal fluid[END_REF][START_REF] Landweber | Added masses and forces on two bodies approaching central impact in an inviscid fluid[END_REF][START_REF] Hicks | On the motion of two cylinders in a fluid[END_REF][START_REF] Greenhill | Functional images in cartesians[END_REF][START_REF] Basset | A Treatise on Hydrodynamics[END_REF][START_REF] Birkhoff | Hydrodynamics[END_REF][START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF] or a conformal transformation [START_REF] Wang | Interaction of two circular cylinders in inviscid fluid[END_REF][START_REF] Burton | Hydrodynamic forces on two moving discs[END_REF][START_REF] Tchieu | Fluid-structure interaction of two bodies in an inviscid fluid[END_REF][START_REF] Scolan | On the use of conformal mapping for the computation of hydrodynamic forces acting on bodies of arbitrary shape in viscous flow. part 2: multi-body configuration[END_REF][START_REF] Crowdy | Analytical solutions for uniform potential flow past multiple cylinders[END_REF][START_REF] Crowdy | A new calculus for two-dimensional vortex dynamics[END_REF][START_REF] Lagrange | A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant fluid[END_REF][START_REF] Lagrange | New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches[END_REF].

Self-added mass and damping coefficients m cross , c cross

Cross-added mass and damping coefficient m Relative deviation between theoretical/experimental and numerical predictions These theoretical approaches have shown their efficiency in predicting the added-coefficients in some simple configurations, but extending to the industrial context remains delicate. As a consequence, engineers developed fast and robust numerical approaches, such as the immersed body method [START_REF] Gronski | A simple and efficient direct forcing immersed boundary method combined with a high order compact scheme for simulating flows with moving rigid boundaries[END_REF], the cut-cell method [START_REF] Cheny | Ls-stag method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF] or the penalization method [START_REF] Lagrange | New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches[END_REF][START_REF] Mittal | Immersed boundary methods[END_REF][START_REF] Kadoch | A volume penalization method for incompressible flows and scalar advection diffusion with moving obstacles[END_REF][START_REF] Kolomenskiy | Two-and three-dimensional numerical simulations of the clap-fling-sweep of hovering insects[END_REF][START_REF] Schneider | Immersed boundary methods for numerical simulation of confined fluid and plasma turbulence in complex geometries: a review[END_REF][START_REF] Minguez | High-order large-eddy simulation of flow over the Ahmed body car model[END_REF][START_REF] Nore | Numerical simulation of the von karman sodium dynamo experiment[END_REF]. In the present work, we aim to analysis the efficiency of the Arbitrary Lagrangian-Eulerian method [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF][START_REF] Fourestey | A second-order time-accurate ALE Lagrange-Galerkin method applied to wind engineering and control of bridge profiles[END_REF][START_REF] Koobus | Computation of unsteady viscous flows around moving bodies using the k-ε turbulence model on unstructured dynamic grids[END_REF] implemented in the open-source code TrioCFD, developed by the CEA for the nuclear industry. TrioCFD is a C++ object-oriented parallel software designed for calculations of unsteady laminar or turbulent fluid flows [START_REF] Calvin | An object-oriented approach to the design of fluid mechanics software[END_REF]. The calculations are performed on structured (parallelepipeds) and non-structured (tetrahedrons) meshes of several millions of control volumes (hybrid finite volume element, see [START_REF] Bieder | A parallel CFD 3-dimensional code for industrial large eddy simulations[END_REF][START_REF] Bieder | An object oriented code for industrial LES[END_REF][START_REF] Ducros | Verification and validation considerations regarding the qualification of numerical schemes for LES for dilution problems[END_REF]). The code structure is flexible, allowing the user to choose the discretization method, the convection and time schemes, as well as the turbulence model. The reader is referred to the TrioCFD webpage [START_REF]TrioCFD: Webpage[END_REF] and [START_REF] Angeli | Overview of the trioCFD code: Main features, V&V procedures and typical applications to nuclear engineering[END_REF] for a detailed presentation of the code.

This article is organized as follows. Section 2 presents the ALE method and the main equations solved by TrioCFD. In Section 3, we test the code, considering the vibration of two coaxial cylinders separated by a fluid layer. The added coefficients are extracted from the numerical predictions of the fluid forces, and the results are compared with those of the literature. In Section 4, we challenge the code, considering the vibration of a cylindrical tube located in the central position of a square tube array. Our numerical results for the fluid coefficients are compared with a phenomenological estimation used by engineers of the nuclear industry and with some experimental measurements that we have performed on a new experimental set-up built at CEA. Finally, Section 5 summarizes our findings.

Numerics

The numerical simulations of the Navier-Stokes equations have been conducted on the 1.8.0 version of TrioCFD, a programmable CFD code based on the TRUST platform and a C++ language architecture, see [START_REF] Angeli | Overview of the trioCFD code: Main features, V&V procedures and typical applications to nuclear engineering[END_REF][START_REF] Angeli | FVCA8 benchmark for the Stokes and Navier-Stokes equations with the trioCFD code -benchmark session[END_REF]. This open source code allows the user to choose a wide range of options and parameters, among which the discretization schemes in space and time, the turbulence models (RANS or LES approaches) or the boundary conditions. Given its flexibility, the code is widely used in the nuclear industry for massive parallel and high performance calculations.

In what follows, we briefly present the ALE method implemented in TrioCFD to solve fluid problems with moving boundaries.

Principle of the ALE numerical method

To determine the flow of a fluid, it is necessary to describe the kinematics of all its material particles throughout time. To do so, one can adopt either an Euler description of motion, in which a fluid particle is identified by its instantaneous position, or a Lagrange description of motion, in which a fluid particle is identified by its initial position. Both descriptions are totally equivalent, leading to different forms of the Navier-Stokes equations that can be discretized on a stationary mesh grid (Euler) or a mesh grid that follows the motion of the fluid particles (Lagrange). In both cases, the mesh grids do not account for the motion of the boundaries, which makes the numerical simulations of the related Navier-Stokes equations delicate.

To overcome this problem, several approaches, such as the immersed boundary methods [START_REF] Puscas | A three-dimensional conservative coupling method between an inviscid compressible flow and a moving rigid solid[END_REF][START_REF] Puscas | A time semi-implicit scheme for the energy-balanced coupling of a shocked fluid flow with a deformable structure[END_REF][START_REF] Puscas | A conservative embedded boundary method for an inviscid compressible flow coupled with a fragmenting structure[END_REF], or the ALE method [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF][START_REF] Fourestey | A second-order time-accurate ALE Lagrange-Galerkin method applied to wind engineering and control of bridge profiles[END_REF][START_REF] Koobus | Computation of unsteady viscous flows around moving bodies using the k-ε turbulence model on unstructured dynamic grids[END_REF] have been developed.

In the ALE approach, a fluid particle is identified by its position relative to a frame moving with a nonuniform velocity v ALE . In this new frame of reference, the Navier-Stokes equations write

     ∇ • v = 0, ∂Jv ∂t = J ν∆v -∇ • ((v -v ALE ) ⊗ v) - 1 ρ ∇p , ( 1a
) (1b)
with J the Jacobian of the transformation between the ALE and the Lagrange descriptions. The ALE method is actually a hybrid description between the Euler and the Lagrange descriptions, both of them corresponding to the particular cases v ALE = 0 and v ALE = v, respectively.

In the ALE framework, the choice of v ALE is arbitrary as long as the deformation of the mesh grid remains under control. For moderate deformations, v ALE is usually defined as the solution of an auxiliary Laplace problem, see [START_REF] Duarte | Arbitrary Lagrangian-Eulerian method for Navier-Stokes equations with moving boundaries[END_REF]:

     ∆v ALE = 0 in the fluid domain, v ALE = v solid at a solid interface, v ALE = 0 at a free surface, (2a) (2b) (2c)
from which the kinematics of the mesh grid is updated, i.e. x new = x old + ∆tv ALE .

Space and time discretizations

In case of unstructured meshes (triangles in 2D or tetrahedrons in 3D), TrioCFD uses Finite Volume-Elements (FVE) approach to solving discretized Navier-Stokes equations. These methods combine the Finite Elements (FE) method with the Finite Volumes (FV) method, gathering the advantages of each approach for incompressible Navier-Stokes problems [START_REF] Angeli | FVCA8 benchmark for the Stokes and Navier-Stokes equations with the trioCFD code -benchmark session[END_REF]. The governing equations (1) are solved through a staggered approaches: a primary grid for the discrete pressure p h and a face-based staggered dual grid for the discrete velocity v h = (v x,h , v y,h ) T . As shown in Fig. 1, the velocity is evaluated on the centre of the 2D edges (or on the center of the 3D faces in 3D) while pressure has more degrees of freedom (DoF) and can be evaluated in the element gravity center and nodes.

As in Finite Volume approach, the local equations are integrated over the control volumes. The control volumes for the mass equation are the primal mesh cells whereas the dual mesh cells (denoted by w hereafter) are the control volumes of momentum. The control volume w associated to each face is obtained by joining the gravity centers of the two adjacent cells G i and G j with the vertices S 1 and S 2 of the shared face (see 1(b)). The fluxes and the differential operators are computed by means of a FE formulation. Hereafter, we recall the main ideas of the FVE method [START_REF] Angeli | Overview of the trioCFD code: Main features, V&V procedures and typical applications to nuclear engineering[END_REF][START_REF] Fiorini | Sensitivity equation method for the navier-stokes equations applied to uncertainty propagation[END_REF] and for this purpose, in this section we only consider the spatial part of the systems [START_REF] Furber | Is the Weis-Fogh principle exploitable in turbomachinery?[END_REF]. Also, in this paragraph, in order to simplify the presentation, we limit ourselves to the two-dimensional case. Let T h be a triangulation of the domain and K j ∈ T h a triangle (j = 1, . . . , N T ). We denote with x i the nodes, which are the middle points of the edges of the triangles and by w i (i = 1, . . . , N N ) their associated control volumes, see Fig. 1(b). We introduce the following spaces:

v h v h v h p h • p h p h p h G i G j • • x i x j • • S 1 S 2 • • K i K j w i (a) (b)
Q h = {q h : ∀K j ∈ T h , q h | Kj ∈ P 0 (K j )}, ( 3a
)
W h = {w h continuous at x i : ∀K j ∈ T h , w h | Kj ∈ P 1 (K j )}, (3b)
W h = {w h = (w x,h , w y,h ) T : w x,h , w y,h ∈ W h } = W 2 h , ( 3c
)
with P n a polynomial function of the n-th order. The space Q h is spanned by the indicator functions of the triangles, χ Kj , and W h is spanned by φ i (x), with φ i ∈ W h and φ i (x j ) = δ ij . We seek an approximate solution for the systems (1):

(v h , p h ) ∈ W h × Q h .
In order to have a discrete formulation, we integrate the mass equation over the triangles K j and the momentum equation over the control volumes ω i as:

         - ∂ωi (ν∇v h - 1 ρ p h [I]) • ndσ i + ∂ωi ((v h -v h,ALE ) ⊗ v h ) • ndσ i = 0 ∀i ∈ [1, . . . , N N ] , ∂Kj v h • ndσ j = 0 ∀K j ∈ T h , (4a) (4b)
with dσ i and dσ j a Lebesgue measure of ∂ω i and ∂K j , respectively, and [I] the identity matrix. Then, one can expand (v h , p h ) in the bases of the correspondent spaces as follows:

v h (x) = NN i=1 v h (x i )φ i (x), p h (x) = NT j=1 p h (K j )χ Kj .
(5)

Plugging ( 5) into (4) for the linear terms and using the identity (a ⊗ b) • c = a(b • c) for the nonlinear terms, results in:

                       - NN j=1 ∂ωi (νv h (x j ) ⊗ ∇φ j ) • ndσ i + NT j=1 1 ρ p h (K j ) ∂ωi χ Kj ndσ i + + ∂ωi v h (v h • n)dσ i - ∂ωi v h (v h,ALE • n)dσ i = 0 ∀i ∈ [1, . . . , N N ] , NN i=1 v h (x i ) • ∂Kj φ i ndσ j = 0 ∀K j ∈ T h . (6a) (6b)
The unknown vectors are defined as follows:

V h = U h = [v x,h (x i )] i=1,...,NN V h = [v y,h (x i )] i=1,...,NN , P h = 1 ρ [p h (K j )] j=1,...,NT , (7)
and the following matrix elements:

[ Ã] ij := - ∂ωi ν∇φ j • ndσ i , [D ℓ ] ij := ∂Kj φ i n ℓ dσ j , [G ℓ ] ij := ∂ωi χ Kj n ℓ dσ i ,
where the subscript ℓ indicates the components of the normal vector n = (n x , n y ).

The convection term can be written as:

∂ωi v h (v h • n)dσ i = ∂ωi v h   NN j=1 v h (x j )φ j • n   dσ i . ( 8
)
The x-projection of (8) yields an equation for the first component v x,h of v h :

∂ωi v x,h   NN j=1 v h (x j )φ j • n   dσ i = NN j=1 v x,h (x j ) ∂ωi v x,h φ j n x dσ i + NN j=1 v y,h (x j ) ∂ωi v x,h φ j n y dσ i .
We can define the following matrix elements:

[L x (U h )] ij = ∂ωi v x,h φ j n x dσ i , [L y (U h )] ij = ∂ωi v x,h φ j n y dσ i .
The second component would give as a result [L x (V h )] ij and [L y (V h )] ij . Finally, by introducing the following notation:

[A] = [ Ã] 0 0 [ Ã] , [D] = [D x ] [D y ] , [G] = [G x ] [G y ] , [L(V h )] = [L x (U h )] [L y (U h )] [L x (V h )] [L y (V h )] , (9)
the discrete system can be written in the compact form:

[A]V h -[L(V h )]V h + [L(V h )]V h,ALE -[G]P h = 0, [D]V h = 0. (10a) (10b)
In the present work, the integral ( 8) is computed numerically by using MUSCL, i.e. a Monotone Upstream-Centred Scheme for Convective flows [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF], see Appendix A. Applying the numerical FVE discretization on Eq. ( 1), along with the Forward Euler scheme for the time discretization, results in the following discrete system:

     [M ] J n+1 V n+1 h -J n V n h ∆t = J n+1 [A]V n+1 h -[L(V n h )]V n+1 h + [L(V n )]V n+1 h,ALE -[G]P n+1 h , [D]V n+1 h = 0, (11a) (11b)
where ∆t is the time step (see Appendix B), [M ] is the mass matrix operator, [A] is the discrete diffusion operator, [L(V h )] is the non-linear discrete convection operator, [G] is the discrete gradient operator and [D] is the discrete divergence operator, defined above. In order to solve the velocity-pressure coupling, a multi-step (projection-correction) technique [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF] is employed, where an intermediate (predicted) velocity V * h is computed:

J n+1 1 ∆t [M ] -[A] + [L(V n h )] V * h = 1 ∆t [M ]J n V n h + J n+1 [L(V n h )]V n h,ALE -J n+1 [G]P n h , (12)
and the mass conservation is then enforced by solving a Poisson equation for pressure:

[D][M ] -1 [G]P ′ h = 1 ∆t [D]V * h . ( 13
)
Finally, the velocity is updated using the predicted velocity V * h and the pressure increment P ′ h :

V n+1 h = V * h -∆t[M ] -1 [G]P ′ h , P n+1 h = P n h + P ′ h . ( 14
)
In the present work, the linear systems ( 12) is solved by the iterative solver GMRES from the PETSc library [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF].

Study 1. Vibrations of two coaxial cylinders in a viscous fluid

In this section, we test the capabilities of the code TrioCFD along with the ALE approach, considering the 2D case of two coaxial cylinders vibrating in a viscous fluid. The goal of this preliminary work is to introduce all the important concepts that will be used in the study of the vibrations of a tube in a square bundle, considered in §4. In what follows, we introduce the problem, present the numerical setup and finally discuss our numerical results in comparison with the theoretical estimations of [START_REF] Chen | Added mass and damping of a vibrating rod in confined viscous fluids[END_REF][START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF].

Presentation of the problem and governing equations

Let C 1 and C 2 be two concentric cylinders with diameters D j and boundaries ∂C j , j = {1, 2}, see Fig. 2. One cylinder, either C 1 or C 2 , oscillates in the (x, y) plane with a simple harmonic motion of angular frequency Ω and a displacement amplitude U . The fluid in the annulus region is Newtonian, homogeneous, with mass density ρ and kinematic viscosity ν. The fluid flow generated by the oscillation of one cylinder is assumed as incompressible and two-dimensional. The Navier-Stokes equations and the boundary conditions for the fluid flow (v, p) are:

             ∇ • v = 0, ∂v ∂t + (v • ∇)v + 1 ρ ∇p -ν∆v = 0, v - dU dt = 0 on ∂C j , j = {1, 2}. (15a) (15b) (15c)
The equation (15c) ensures the continuity of velocities at the cylinder boundaries ∂C j . The fluid force acting on C j is the sum of a pressure and a viscous term

F j = - ∂Cj pn j dL j + ρν ∂Cj [∇v + (∇v) T ] • n j dL j , ( 16
)
where n j is the outward normal unit vector to ∂C j , see Fig. 2, and dL j is an infinitesimal line element of integration. Picking D 1 /2 and Ω -1 as a characteristic length and time, the dimensionless quantities t * , u * , v * , p * and f * j are defined as

t = Ω -1 t * , U = U u * , v = U Ωv * , p = ρU D 1 2 Ω 2 p * , F j = ρU D 1 2 Ω 2 f * j . ( 17
)
To reduce the number of parameters of the problem, the following rescaled quantities are also defined

ε = D 2 D 1 , KC = U D 1 , Sk = D 2 1 (Ω/2π) ν , ( 18
)
as the dimensionless separation distance, the Keulegan-Carpenter number and the Stokes number, respectively. Introducing ( 17) in ( 15) yields the dimensionless Navier-Stokes equations

           ∇ • v * = 0, ∂v * ∂t * + 2KC(v * • ∇)v * + ∇p * - 2 πSk ∆v * = 0, v * - du * dt * = 0 on ∂C j , j = {1, 2}. ( 19a
) (19b) (19c)
The dimensionless fluid force acting on C j is

f * j = - ∂Cj p * n j dl j + 2 πSk ∂Cj [∇v * + (∇v * ) T ] • n j dl j , ( 20
)
with dl j = 2dL j /D 1 .

For KC << 1, the dimensionless Navier-Stokes equations are linear. It follows that the fluid forces are linear combinations of the cylinder velocity du * /dt * and acceleration

d 2 u * /dt * 2 f * j = -π [M j ] d 2 u * dt * 2 + [C j ] du * dt * , (21)
with

[M j ] = m (j) self m (j) self and [C j ] = c (j) self c (j) self if C j is moving, (22a) [M j ] = m cross -m cross and [C j ] = c cross -c cross if C j is stationary. ( 22b
)
The self-added mass and damping coefficients m For a harmonic displacement u * = sin (t * )e x , the dimensionless fluid forces (21) reduce to

f * j = π m (j) self sin (t * ) -c (j) self cos (t * ) e x , if C j is moving, (23a)
f * j = π (m cross sin (t * ) -c cross cos (t * )) e x , if C j is stationary, (23b)
or similarly f * j = h j sin (t * + ϕ j ) e x , with

h j = π m (j) self 2 + c (j) self 2 and ϕ j = -arctan c (j) self m (j) self if C j is moving, (24a)
h j = π m 2 cross + c 2 cross and ϕ j = π -arctan c cross m cross if C j is stationary. (24b)
Yeh et al. [START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF] derived some analytical expressions of the fluid added coefficients. Introducing ℜ and ℑ the real and imaginary part operators, these expressions write

m (1) self = ℜ(a 11 ), m (2) self = ε 2 ℜ(a 22 ), m cross = εℜ(a 12 ), (25a) c (1)
self = ℑ(-a 11 ), c (2) self = ε 2 ℑ(-a 22 ), c cross = εℑ(-a 12 ), (25b)
with

a 11 = -(1 + 2a), a 22 = 2a/ε, a 12 = 1 -2a/ε 2 and a = 1 1 J 1 (β 1 ) Y 1 (β 1 ) 0 1 ε -1 J 1 (β 2 ) ε -1 Y 1 (β 2 ) 2 2 β 1 J 0 (β 1 ) β 1 Y 0 (β 1 ) 0 2 β 1 J 0 (β 2 ) β 1 Y 0 (β 2 ) 1 1 J 1 (β 1 ) Y 1 (β 1 ) ε -2 1 ε -1 J 1 (β 2 ) ε -1 Y 1 (β 2 ) 0 2 β 1 J 0 (β 1 ) β 1 Y 0 (β 1 ) 0 2 β 1 J 0 (β 2 ) β 1 Y 0 (β 2 ) -1 . ( 26
)
In [START_REF] Wang | Interaction of two circular cylinders in inviscid fluid[END_REF], J n and Y n are the Bessel functions of the first and second kind, with arguments

β 1 = (1 -i) √ π/2 √ Sk or β 2 = εβ 1 .
From an asymptotic expansion of these functions as Sk → ∞, we show that (m

(j) self , m cross , c (j) self , c cross ) are equivalent to ( m (j) self , m cross , c (j) self , c cross ) with m (j) self = m (j)pot self + 4 √ π 1 √ Sk , m cross = m pot cross - 4 √ π 1 √ Sk , (27a) c (j) self = 4 √ π 1 √ Sk ε 4 + ε (ε 2 -1) 2 , c cross = - 4 √ π 1 √ Sk ε 4 + ε (ε 2 -1) 2 , ( 27b
)
and

m (1)pot self = ε 2 + 1 ε 2 -1 , m (2)pot self = ε 2 (ε 2 + 1) ε 2 -1 , m pot cross = - 2ε 2 ε 2 -1 . ( 28
)
The terms m (j)pot self and m pot cross are the inviscid limits of m (j) self and m cross as Sk → ∞. These terms can also be obtained from a potential theory in which the fluid forces are only due to the pressure field. Note that c (j)pot self = c pot cross = 0 as there are not damping effects in the inviscid framework.

In §3.3 the theoretical estimations ( 25) and ( 27) are compared with the numerical predictions performed with TrioCFD. Numerically, the fluid forces acting on ∂C j are computed by the sum of the pressure and viscous terms in Eq. ( 16) given by TrioCFD. To extract the added coefficients from the numerical simulations of the dimensionless fluid forces, we introduce the Fourier inner product over five periods

f (t * ), g(t * ) = 1 5π 10π 0 f (t * )g(t * )dt * . ( 29
)
From [START_REF] Basset | A Treatise on Hydrodynamics[END_REF], it follows that the added-coefficients are

m (j) self = sin(t * ), f * j (t * ) • e x π and c (j) self = - cos(t * ), f * j (t * ) • e x π if C j is moving, ( 30a
)
m cross = sin(t * ), f * j (t * ) • e x π and c cross = - cos(t * ), f * j (t * ) • e x π if C j is stationary. (30b)

Numerical setup

To discretize the fluid domain, we use an unstructured grid of triangles generated by the Gmsh platform, see [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in preand post-processing facilities[END_REF]. Due to its adaptive and automatic algorithm, Gmsh makes it possible to choose between two different local sizes for the 2D mesh: a small local size, lc f ine , for elements close to the moving cylinder and a large local size, lc, for elements close to the stationary cylinder. By this way, a refined mesh is used in the regions with large gradient fields whereas a loose mesh is used in the areas with low gradient fields. In Appendix E, a mesh sensitivity analysis clearly shows a convergence of the mass coefficients as l c and lc f ine are changed. The convergence of the damping coefficients is less obvious, especially for high values of Sk. Physically, this is related to the thickness of the boundary layer, which tends to zero as Sk increases. It follows that a finer mesh is required close to a cylinder boundary to account for the thickness of the boundary layer and obtain an accurate estimation of the damping terms. In this work, a compromise between the time of calculation and the precision (throughout of this work, we indicate the relative deviations of our numerical results compared to some theoretical or experimental references) is made to choose adequate values for lc and lc f ine .

Results and discussion

The inner cylinder is imposed a sinusoidal displacement in the x-direction for five periods of time, i.e. u * = sin(t * )e x with t * ∈ {0, 10π}. The numerical simulations are performed for ε ∈ {1.25, 1.5, 2}, Sk ∈ {10 1 , 10 2 , 10 3 , 10 4 } and KC = 10 -2 .

The time evolution of the dimensionless fluid forces is represented in Figs. 3 and4. A perfect agreement is obtained between theoretical estimations and numerical predictions. The fluid forces are sinusoidal functions in phase opposition. In Figs. 5 and6, we show that the amplitudes h j (resp. phases ϕ j ) decrease (resp. increase) with both Sk and ε, recovering the asymptotic limits as Sk → ∞ (inviscid limit) and as ε → ∞ (isolated cylinder limit).

The evolution of the added coefficients is depicted in Figs. 7 and8. The numerical predictions are in very good agreement with the theoretical estimations, even if a tiny difference is observed in the range of low Stokes numbers in which the theoretical estimations reach their limit of validity. Still, we confirm that the self (resp. cross) added coefficients decrease (resp. increase) with both Sk and ε, recovering the asymptotic limits as Sk → ∞ (inviscid limit) and as ε → ∞ (isolated cylinder limit). . The red and blue solid lines correspond to the theoretical estimations of [START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq. [START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF]. The dashed lines (indistinguishable from the solid lines for Sk > 10 2 ) correspond to the theoretical asymptotic expansion Sk → ∞, see Eq. [START_REF] Burton | Hydrodynamic forces on two moving discs[END_REF]. The symbols correspond to the numerical predictions.

£ uid force, f* j C 1 C 2 C 1 C 2 C 1 C 2 C 1 C 2 -25 -20 -15 -10 -
We introduce the quantity ι, defined as the relative deviation between the numerical and the exact theoretical predictions [START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF] of some quantity Q : ι = |Q th. -Q num. |/|Q th. |. The Fig. 9 and the Tabs. D.2, D.3 in Appendix D show that ι is very low for the added mass coefficients and it becomes greater for the added damping coefficients as Sk increases but seems to be weakly influenced by ε. In fact, for large values of Sk, the damping coefficients become very small and even small differences from the theoretical value produce significant deviations. Physically, this is related to the thickness of the boundary layer, which tends to zero: to reduce ι for the damping coefficients a finer mesh at the cylinder boundaries would be required. We shall note that the theoretical approach [START_REF] Chen | Added mass and damping of a vibrating rod in confined viscous fluids[END_REF][START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF] is fully linear since the convective term 2KC(v * • ∇) v * of the Navier-Stokes equation (19b) is neglected. In the numerical simulations, the nonlinear convective term is retained through a small but nonzero Keulegan-Carpenter number KC = 10 -2 . However, this difference might slightly affect the deviation between the theoretical and numerical results. In any case, the relative deviation for m [START_REF] Furber | Is the Weis-Fogh principle exploitable in turbomachinery?[END_REF] self (resp. m cross ) is smaller than ι ≤ 0.6% (resp. ι ≤ 0.2%) while the relative deviation for c [START_REF] Furber | Is the Weis-Fogh principle exploitable in turbomachinery?[END_REF] self (resp. c cross ) is more pronounced, with ι ≤ 5.5% (resp. ι ≤ 6.7%).

A further analysis is carried out to check the scale invariance of the numerical results obtained for ε = 2, Sk = 10 4 and KC = 10 -2 . In the numerical simulation, the dimensional quantities were : self , c cross = (6.71, -2.71, 0.0460, -0.0460). It follows that the relative deviation between the two approaches is ι ≤ 12.8%. We attribute this deviation to the mesh refinement and the larger extension of the surface on which the ALE module acts, compared to the case in which the inner cylinder moves. In fact, in the ALE approach we need to solve the auxiliary Laplace problem (2) to compute the mesh velocity, with its own approximations and therefore with an additional numerical error. We shall note that a finer mesh would have reduced this error and normally we would not have observed such differences. The red and blue solid lines correspond to the theoretical estimations of [START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq. ( 25). The dashed lines correspond to the theoretical asymptotic expansion Sk → ∞, see Eq. ( 27). The horizontal dashed lines correspond to the theoretical inviscid limits as Sk → ∞. The symbols correspond to the numerical predictions.

D 1 = 1, D 2 = 2, U = 10
Magnitude, h j Phase, φ j (rad) and KC = 10 -2 . The red and blue solid lines correspond to the theoretical estimations of [START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq. ( 25). The dashed lines correspond to the theoretical asymptotic expansion Sk → ∞, see Eq. ( 27). The horizontal dashed lines correspond to the theoretical isolated cylinder limits as ε → ∞. The symbols correspond to the numerical predictions.

As expected, the numerical simulations confirm that the cross added coefficients do not depend on which cylinder is moving. Indeed, for ε = 2, Sk = 10 4 and KC = 10 As a conclusion to study 1, we show in Fig. 10 the pressure distribution in the fluid layer when C 1 is imposed a sinusoidal displacement in the x-direction and C 2 is stationary, for Sk = 10 4 , ε = 2 and KC = 10 -1 . The numerical and theoretical predictions are in very good agreement as the two snapshots look very similar. Note that these snapshots are taken at a time when C 1 has a positive acceleration. As expected from the inertia effect of the added mass, the fluid pushes C 1 in a direction opposite to its acceleration vector, by creating a negative (resp. positive) pressure on the left (resp. right) side of C 1 .

In this section we have thoroughly tested the capabilities of our CFD code in predicting the fluid forces in the case of two coaxial cylinders. In what follows we challenge the code with a more complicated configuration c (1) self m cross m (1) self c cross . The solid lines refer to the theoretical estimations of [START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq. [START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF]. The dashed lines correspond to the theoretical asymptotic expansion Sk → ∞, see Eq. ( 27). The horizontal dashed lines correspond to the theoretical inviscid limits as Sk → ∞. The symbols correspond to the numerical predictions.

consisting of a tube bundle immersed in a viscous fluid.

Study 2. Vibrations of a cylinder in a square tube bundle immersed in a viscous fluid

The importance of fluid-elastic forces in tube bundle vibrations can hardly be over-emphasized, in view of their damaging potential. In the last decades, advanced models for representing the fluid-elastic forces through added-coefficients have therefore been developed by the community of the domain. Those models are nowadays embedded in the methodologies that are used on a regular basis by both steam generators providers and operators, in order to prevent the risk of a tube failure with adequate safety margins. From an R&D point of view however, the need still remains for more advanced models of the fluid added-coefficients, in order to fully decipher the physics underlying the observed phenomena.

In what follows, we aim to determine numerically and experimentally (new measurement setup built at CEA) the effect of the Stokes number Sk on the fluid coefficients, considering the case of a square tube bundle immersed in a viscous fluid at rest. The numerical results obtained with TrioCFD are thoroughly compared with the experimental results and also a theoretical approach proposed by [START_REF] Pettigrew | Flow-induced vibration specifications for steam generators and liquid heat exchangers[END_REF].

Presentation of the problem

Let C C be a cylinder of diameter D, located in the central position of a square tube bundle with pitch P , see Fig. 11 a). The cylinder C C oscillates in the (x, y) plane with a simple harmonic motion of angular c (1) self m cross m (1) self c cross The solid lines refer to the theoretical estimations of [START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq. ( 25). The dashed lines correspond to the theoretical asymptotic expansion Sk → ∞, see Eq. ( 27). The horizontal dashed lines correspond to the theoretical isolated cylinder limits as ε → ∞.

C D E F G H I P Q R S T U V m (1) self W X Y `a b c d e f h i p q r s t c (
The symbols correspond to the numerical predictions.

frequency Ω and a displacement amplitude U . All the other cylinders are stationary. The tube bundle is immersed in a Newtonian and homogeneous fluid, with mass density ρ and kinematic viscosity ν. The fluid flow generated by the oscillation of the central cylinder is assumed as incompressible and two-dimensional. Following the approach presented in §3 with D 1 = D 2 = D, the dimensionless fluid force acting on the central cylinder writes

f * = -π m self d 2 u * dt * 2 + c self du * dt * , ( 31
)
with m self and c self the self-added mass and damping dimensionless coefficients, respectively. For a given configuration, these coefficients are functions of the pitch ratio P/D, the Keulegan-Carpenter number KC = U/D and the Stokes number Sk = D 2 Ω/(2π/ν).

In order to determine the variations of m self and c self with Sk, an experimental setup consisting of 3 × 5 cylinders (plus two columns of 5 half-cylinders at the boundaries) has been built at CEA, see Fig. 11(b). This setup is briefly presented in the following and the reader is referred to [START_REF] Caillaud | The measurement of fluidelastic effects at low reduced velocities using piezoelectric actuators[END_REF][START_REF] Caillaud | Active vibration control for the measurement of fluidelastic effects[END_REF][START_REF] Piteau | Experiments and computations of a loosely supported tube in a rigid bundle subjected to single-phase flow[END_REF] for an extensive description. Following the principle of the direct method developed by [START_REF] Tanaka | Fluidelastic vibration of tube arrays in cross-flow[END_REF], an electro-dynamic shaker (PRODERA 200 N) is used to impose a harmonic motion to the central tube, supported by two parallel flexible steel blades allowing large vibrations in the x direction. The fluid coefficients are directly extracted from the measure of the fluid force, obtained with two piezo-electric transducers (KISTLER 9132 A). Despite its apparent simplicity, the method of [START_REF] Tanaka | Fluidelastic vibration of tube arrays in cross-flow[END_REF] has been progressively abandoned due to its difficulty to implement experimentally. Still, here we reconsider this method, performing some new experiments with a tube bundle with pitch P = 0.045 m, made with cylinders of diameter D = 0.03 m and length L = 0.3 m. The central cylinder is imposed a harmonic displacement in the x-direction with amplitude U = 0.003 m and a frequency F ∈ {5, 10, 15, 20, 25} Hz. In terms of dimensionless numbers, this corresponds to a pitch ratio P/D = 1.5, a Keulegan-Carpenter In what follows, we aim to compare the experimental results with TrioCFD predictions. In the discussion of our results we also refer to the phenomenological estimations of the fluid coefficients provided by [START_REF] Pettigrew | Flow-induced vibration specifications for steam generators and liquid heat exchangers[END_REF] for KC << 1

m self = π 4 (D e /D) 2 + 1 (D e /D) 2 -1 , c self = 4 √ π 1 √ Sk 1 + (D/D e ) 3 [1 -(D/D e ) 2 ] 2 , ( 32
)
where D e = P (1.07 + 0.56P/D) represents an equivalent diameter introduced to model the confinement of the tube bundle.

Numerical setup

The fluid domain is discretized with an unstructured grid of triangles generated by the SALOMÉ platform, see [START_REF] Ribes | Salomé platform component model for numerical simulation[END_REF].

Before analyzing the meshes used in the numerical simulations and their properties, we start discussing on the technique used to design the geometry. In the experimental setup, the geometry is perfectly symmetrical so that there is no fluid force acting on the cylinders as they are stationary. Numerically, to ensure a zero hydrostatic force, the same discretization scheme of the fluid domain close to every single cylinder is required, while maintaining an x and y symmetry of the mesh grid. To do so, we introduce some artificial lines to subdivide the fluid domain into blocks, each of them being divided into four slices. The automatic and optimized mesh algorithm ensures that the arrangement of the generated triangles, on each side of the artificial lines, is perfectly symmetrical, so that all the slices are discretized equally. The MG-CADSurf Module of SALOMÉ is used for this purpose. We define two mesh sizes: a local size, lc f ine , for elements close to the cylinders and a user size, lc for elements far from the cylinders. A gradation parameter is also set, ∆ = 1.1, as the maximum ratio between the lengths of two adjacent edges.

Firstly, a mesh sensitivity analysis is carried out on a truncated fluid domain (with the last rows of cylinders cut in half) in order to reduce the time of calculation and speed up the convergence of the grid. This analysis shows that a mesh grid with lc = 0.0005 and lc f ine = 0.00015 is suitable for an accurate estimation of the fluid forces and the related added coefficients, see Appendix E. Then, the mesh sizes are kept constant whereas the size of the fluid domain is gradually increased by a distance of N D diameters in the y-directions, with N D from 1 to 5.

Results and discussion

The central cylinder is allowed to move in the x-direction with a sinusoidal displacement of five periods of time, i.e. u * = sin(t * )e x with t * ∈ {0, 10π}. The numerical simulations are performed for Sk ∈ {0.45, 0.90, 1.35, 1.80, 2.25} × 10 4 corresponding to imposed frequencies in the range F ∈ {5, 10, 15, 20, 25} Hz. The other parameters are P/D = 1.5 and KC = 10 -1 .

Before analyzing the variations of the fluid coefficients, we start discussing on the choice of the boundary conditions used in the numerical simulations. In our experiments, the fluid is laterally confined by two parallel plates. Numerically, a classical wall condition (v = 0) is thus applied on the two lateral sides. The choice of the boundary condition at the top and bottom sides of the computational domain (y directions) is more delicate as the tube bundle is experimentally surrounded by two columns of water. To overcome this problem, we determine the size of the computational domain which yields numerical results poorly sensitive to this choice. To do so, we test the effect of two types of y boundary conditions: a wall condition and a Neumann zero stress condition (pn -ρν[∇v + (∇v) T ] • n = 0). In Fig. 12, we show the results of this test for Sk = 1.8 × 10 4 as the size of the fluid domain increases by a distance of N D diameters in the y-directions.

As N D increases, we observe that m self converges to a constant value that does not depend on the choice of the y boundary condition. Similarly, the two types of boundary conditions yield very tiny differences on the values of c self . From these observations, we thus decide to perform all our numerical simulations with N D = 4 and a wall boundary condition in the y-directions. The time evolution of the dimensionless fluid force is represented in Fig. 13. We note that our simulations are in good agreement with both the theoretical and the experimental estimations, despite a small amplitude deviation, in particular with the experimental measurements for F ∈ {5, 20} Hz, i.e. Sk ∈ {0.45, 1.8} × 10 4 . This small deviation will be analyzed more precisely in what follows, considering the variations of the fluid coefficients m self and c self .

The evolution of the self-added coefficients with the Stokes number, Sk, is depicted in Fig. 14. Our simulations are in good agreement with the theoretical estimations, in the sense that similar trends are recovered, bringing out the same behavior of the fluid coefficients. Still, we note that the simulations tend to underestimate (resp. overestimate) m self (resp. c self ), with a relative deviation always smaller than ι ≤ 3.1% (resp. ι ≤ 27.5%), see Tab. 1. We attribute these deviations to the fact that the theoretical approach of [START_REF] Pettigrew | Flow-induced vibration specifications for steam generators and liquid heat exchangers[END_REF] is based on a strong approximation in which the tube bundle is modeled as a system of two concentric cylinders, through the definition of a fictive equivalent diameter determined empirically. Also, the theory of [START_REF] Pettigrew | Flow-induced vibration specifications for steam generators and liquid heat exchangers[END_REF] relies on a linear approach, contrarily to our numerical simulations in which the nonlinear convective term of the . The blue solid lines correspond to the theoretical estimations of [START_REF] Pettigrew | Flow-induced vibration specifications for steam generators and liquid heat exchangers[END_REF], see Eq. [START_REF] Lagrange | A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant fluid[END_REF]. The red solid lines correspond to the experimental measurements. The orange symbols correspond to the numerical predictions. (For interpretation of the colors in this figure, the reader is referred to the web version of this article).

dimensionless Navier-Stokes equation is taken into account. The agreement between the numerical predictions and our experimental results is partial, showing a good agreement for the self-added mass coefficient m self , with ι < 7%, but some important deviation for c self , especially for the lowest Stokes number Sk = 0.45 × 10 where ι ≈ 105%. It shall be noted that it is very delicate to measure the self-damping coefficient in our experiments as the determination of c self relies on a tiny phase difference between the fluid force and the imposed displacement. Also, for such a low Stokes number, i.e. low frequency, our experiments become very sensitive to parasitic frequencies in the experimental setup, leading to a bad signal to noise ratio. Still, these first comparisons are very encouraging and should foster further developments of the experimental direct method.

As a conclusion to study 2, we show in Fig. 15 the pressure distribution when C C is imposed a sinusoidal displacement in the x-direction, for Sk = 1.8 × 10 4 , P/D = 1.5 and KC = 10 -1 . Note that this snapshot is taken at a time when C C has a positive acceleration. As expected from the inertia effect of the added mass, the fluid pushes C C in a direction opposite to its acceleration vector, by creating a negative (resp. positive) pressure on the left (resp. right) side of C C .

Conclusion and perspectives

In this study, we have first considered the vibration of two coaxial cylinders separated by a viscous fluid. The ALE numerical approach, based on a mesh update technique that considers the motion of the boundary, has been carried out with TrioCFD to estimate the fluid forces acting on the two cylinders.

We have started studying the case in which the outer cylinder is stationary while the inner one is imposed a harmonic motion. The numerical predictions for the fluid forces and the corresponding added coefficients are in very good agreement with the theoretical estimations, even if a tiny difference is observed in the range of low Stokes numbers in which the theoretical estimations reach their limit of validity. More specifically, we have shown that the two forces are in phase opposition and sensitive to the Stokes number, Sk, and the dimensionless separation distance, ε. The amplitude (resp. phase) of the forces decreases (resp. increases) as Sk or ε increases, recovering the inviscid limit for large Stokes numbers or the limit of an isolated cylinder for large ε. The fluid coefficients variations are also correctly reproduced by the numerical simulations, which are shown to respect the scale invariance expected from the dimensional analysis. Secondly, we have studied the case in which the inner cylinder is stationary while the outer one is imposed a harmonic motion. From this study, we have shown that the force acting on the stationary cylinder does not depend on which cylinder is moving (inner or outer).

After having successfully tested the capabilities of the ALE method implemented in TrioCFD, in the configuration of two coaxial cylinders, we have considered the problem of a vibrating cylinder located in the central position of a square tube bundle immersed in a viscous fluid. We have shown that the numerical predictions for the self-added mass coefficient are in good agreement with the results of our experiments led at CEA and a theoretical estimation used by engineers. The numerical predictions for the self-added damping coefficient are also in good agreement with the theoretical estimation. They partially agree with the experimental results, showing an important divergence for low Stokes numbers. However, such a deviation is mainly due to our experimental setup, which exhibits a strong sensitivity to parasitic frequencies for low Stokes numbers.

In conclusion, the ALE method implemented in TrioCFD is particularly efficient in solving fluid-structure interaction problems with a fluid initially at rest. We are now confident that further developments will make it possible to tackle three-dimensional problems with an incident fluid flow. m (1) self c (1) self m cross c cross m (1) self c (1) self m cross c cross m (1) self c (1)

  self as Sk → ∞ m cross , c cross Asymptotic expansion of m cross and c cross as Sk → ∞ m (j)pot self , c (j)pot self Inviscid limits of m (j) self and c (j) self m pot cross , c pot cross Inviscid limits of m cross and c cross ι

Figure 1 :

 1 Figure 1: (a) DoF of 2D element (black squares for velocity v h , black dots and circles for pressure p h ). (b) Control volume w i between two triangles K i and K j of respective center G i and G j .

Figure 2 :

 2 Figure 2: Configuration with two concentric cylinders. The fluid in the annulus region is homogeneous with mass density ρ and kinematic viscosity ν.

  self relate the fluid force on the moving cylinder to its own motion. The cross-added mass and damping coefficients m cross and c cross relate the fluid force on the stationary cylinder to the motion of the other cylinder. The fluid added coefficients are functions of the dimensionless separation distance ε and the Stokes number Sk.

2 Figure 3 :

 23 Figure 3: Study 1. Time evolution of the dimensionless fluid forces, for Sk = {10 1 , 10 2 , 10 3 , 10 4 }, ε = 2 and KC = 10 -2. The red and blue solid lines correspond to the theoretical estimations of[START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq.[START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF]. The dashed lines (indistinguishable from the solid lines for Sk > 10 2 ) correspond to the theoretical asymptotic expansion Sk → ∞, see Eq.[START_REF] Burton | Hydrodynamic forces on two moving discs[END_REF]. The symbols correspond to the numerical predictions.

Figure 4 : 2 )

 42 Figure 4: Study 1. Time evolution of the dimensionless fluid forces, for ε = {1.25, 1.5, 2}, Sk = 10 4 and KC = 10 -2. The red and blue solid lines correspond to the theoretical estimations of[START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq.[START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF]. The symbols correspond to the numerical predictions.

Figure 5 :

 5 Figure 5: Study 1. Evolution of the magnitudes h j and phases ϕ j with the Stokes number Sk, for ε = 2 and KC = 10 -2 .The red and blue solid lines correspond to the theoretical estimations of[START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq.[START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF]. The dashed lines correspond to the theoretical asymptotic expansion Sk → ∞, see Eq. (27). The horizontal dashed lines correspond to the theoretical inviscid limits as Sk → ∞. The symbols correspond to the numerical predictions.

Figure 6 :

 6 Figure6: Study 1. Evolution of the magnitudes h j and phases ϕ j with the dimensionless separation distance ε, for Sk = 10 4 and KC = 10 -2 . The red and blue solid lines correspond to the theoretical estimations of[START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq.[START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF]. The dashed lines correspond to the theoretical asymptotic expansion Sk → ∞, see Eq.[START_REF] Burton | Hydrodynamic forces on two moving discs[END_REF]. The horizontal dashed lines correspond to the theoretical isolated cylinder limits as ε → ∞. The symbols correspond to the numerical predictions.

  -2 , when C 1 moves we obtained (m cross , c cross ) = (-2.71, -0.0488) whereas when C 2 moves we have (m cross , c cross ) = (-2.70, -0.0488). In other words, the dimensionless fluid force acting on C 1 as C 2 moves is the same as the dimensionless fluid force acting on C 2 as C 1 moves, as shown in Fig. C.18(b).

Figure 7 :

 7 Figure 7: Study 1. Evolution of the added coefficients with the Stokes numbers Sk, for ε = 2 and KC = 10 -2. The solid lines refer to the theoretical estimations of[START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq.[START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF]. The dashed lines correspond to the theoretical asymptotic expansion Sk → ∞, see Eq.[START_REF] Burton | Hydrodynamic forces on two moving discs[END_REF]. The horizontal dashed lines correspond to the theoretical inviscid limits as Sk → ∞. The symbols correspond to the numerical predictions.

  Dimensionless separation distance, εDimensionless separation distance, εDimensionless separation distance, εDimensionless separation distance, ε

Figure 8 :

 8 Figure 8: Study 1. Evolution of the added coefficients with the dimensionless separation distances ε, for Sk = 10 4 and KC = 10 -2 .The solid lines refer to the theoretical estimations of[START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq.[START_REF] Gibert | Vibration of structures in a static fluid medium[END_REF]. The dashed lines correspond to the theoretical asymptotic expansion Sk → ∞, see Eq.[START_REF] Burton | Hydrodynamic forces on two moving discs[END_REF]. The horizontal dashed lines correspond to the theoretical isolated cylinder limits as ε → ∞. The symbols correspond to the numerical predictions.

Figure 9 :Figure 10 :

 910 Figure 9: Study 1. Evolution of the relative deviation, ι with the Stokes number, Sk for different value of ε and KC = 10 -2 . The red dotted lines with points refer to ε = 2, the blue points refer to ε = 1.5 while the orange points refer to ε = 1.25.

Figure 11 :

 11 Figure 11: Study 2. (a) Configuration with 3 × 5 cylinders plus two columns of 5 half-cylinders at the boundaries. The tube bundle is immersed in a homogeneous fluid with mass density ρ and kinematic viscosity ν. (b) Photo of the experimental setup designed to measure the fluid coefficients by the direct method.

Figure 12 :

 12 Figure 12: Study 2. Evolution of the self-added coefficients with the size of the fluid domain (the size of the fluid domain increases by a distance of N D diameters in the y-directions). The blue (resp. red) dashed line corresponds to the numerical predictions obtained with a wall boundary condition (resp. a Neumann boundary condition) imposed on the top and bottom boundaries of the computational domain, see Fig. 11. The dimensionless parameters are Sk = 1.8 × 10 4 , P/D = 1.5 and KC = 10 -1 .

Figure 13 :

 13 Figure 13: Study 2. Time evolution of the dimensionless fluid forces, for Sk ∈ {0.45, 0.90, 1.35, 1.80, 2.25} × 10 4 , P/D = 1.5 and KC = 10 -1. The blue solid lines correspond to the theoretical estimations of[START_REF] Pettigrew | Flow-induced vibration specifications for steam generators and liquid heat exchangers[END_REF], see Eq.[START_REF] Lagrange | A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant fluid[END_REF]. The red solid lines correspond to the experimental measurements. The orange symbols correspond to the numerical predictions. (For interpretation of the colors in this figure, the reader is referred to the web version of this article).

Figure 15 :

 15 Figure 15: Study 2. Dimensionless pressure field, for Sk = 1.8 × 10 4 , P/D = 1.5 and KC = 10 -1 , computed with TrioCFD.

Figure C. 18 :

 18 Figure C.18: Study 1. Testing the symmetry of the fluid added mass and damping matrices. (a) Time evolution of the dimensionless fluid forces, for Sk = 10 4 , ε = 2 and KC = 10 -2. The red and blue solid lines correspond to the theoretical estimations of[START_REF] Yeh | The effect of fluid viscosity on coupled tube/fluid vibrations[END_REF], see Eq. (25). The symbols correspond to the numerical predictions. (b) Time evolution of the dimensionless fluid force on the stationary cylinder, for Sk = 10 4 , ε = 2 and KC = 10 -2 . The solid line corresponds to the numerical predictions as C 1 oscillates. The symbols correspond to the numerical predictions as C 2 oscillates.

Table D. 2 : 6 :

 26 Study 1. Case 1. Table of the fluid added coefficients as a function of the Stokes number, Sk and relative deviation, ι. The notations ET and AT refer to the exact and asymptotic theories. The numerical values refer to symbols shown in Fig. 7. The dimensionless separation distance is ε = 2 and the Keulegan-Carpenter number is KC = 10 -2 . Study 1. Effect of the mesh size on the fluid added coefficients. The dimensionless separation distance is ε = 1.25 and the Keulegan-Carpenter number is KC = 10 -2 .

4 : 5 : 7 : 8 : 9 : 10 :Table E. 11 :

 457891011 Study 1. Testing the scale invariance. Table of the fluid added coefficients and the relative deviation, ι. The notations ET and AT refer to the exact and asymptotic theories. The numerical values are extracted from the dimensionless fluid forces shown in Fig. C.17. The Stokes number is Sk = 10 4 , the dimensionless separation distance is ε = 2 and the Keulegan-Carpenter number is KC = 10 -2 . Study 1. Testing the symmetry of the fluid added mass and damping matrices. Table of the fluid added coefficients and the relative deviation, ι. The notations ET and AT refer to the exact and asymptotic theories. The numerical values are extracted from the dimensionless fluid forces shown in Fig. C.18(b). The Stokes number is Sk = 10 4 , the dimensionless separation distance is ε = 2 and the Keulegan-Carpenter number is KC = 10 -2 . Study 1. Effect of the mesh size on the fluid added coefficients. The dimensionless separation distance is ε = 1.5 and the Keulegan-Carpenter number is KC = 10 -2 . Study 1. Effect of the mesh size on the fluid added coefficients. The dimensionless separation distance is ε = 2 and the Keulegan-Carpenter number is KC = 10 -2 . Study 1. Testing the scale invariance. Effect of the mesh size on the fluid added coefficients. The dimensionless separation distance is ε = 2 and the Keulegan-Carpenter number is KC = 10 -2 . Study 1. Testing the symmetry of the fluid added mass and damping matrices. Effect of the mesh size on the fluid added coefficients. The dimensionless separation distance is ε = 2 and the Keulegan-Carpenter number is KC = 10 -2 . Study 2. Effect of the mesh size on the fluid added coefficients. The pitch ratio is P/D = 1.5 and the Keulegan-Carpenter number is KC = 10 -1 .

Table 1 :

 1 Study 2. Evolution of the self-added coefficients with the Stokes number Sk, for P/D = 1.5 and KC = 10 -1 . The red, blue and orange dashed lines with points refer respectively to the experimental measurements, the theoretical estimations of[START_REF] Pettigrew | Flow-induced vibration specifications for steam generators and liquid heat exchangers[END_REF] and the numerical predictions. Study 2. Table of the self-added coefficients as a function of the Stokes number, Sk and relative deviation, ι. The notations ER and AT refer to the experimental results and the approximate theory. The pitch ratio is P/D = 1.5 and the Keulegan-Carpenter number is KC = 10 -1 .

	self
	Self added mass, m

4

,

  Table D.3: Study 1. Table of the fluid added coefficients as a function of the dimensionless separation distance, ε and relative deviation, ι. The notations ET and AT refer to the exact and asymptotic theories. The numerical values refer to symbols shown in Fig. 8. The Stokes number is Sk = 10 4 and the Keulegan-Carpenter number is KC = 10 -2 .

Acknowledgements

The authors fully acknowledge P. Piteau and T. Valin (CEA Saclay) for implementing the Tanaka's direct method of measurement on the DIVA experimental setup located at CEA Saclay, and for sharing their results as part of the present study. The authors also acknowledge financial support for this work, which was performed in the framework of a joint research program co-funded by FRAMATOME, EDF and CEA (France).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Discretization of the non-linear convection term: MUSCL scheme

In this paragraph, in order to simplify the presentation, we restrict ourselves to the two-dimensional case. We recall that in the VDF formulation with each degrees of freedom x i of the speed, we associate a control volume w i (see Fig. A.16). We have:

where γ i represents the faces of the control volume w i . In the MUSCL scheme (Monotone Upstream-Centred Scheme for Convective flows) [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF], the flow across the face γ ij is approximated based on the Simpson interpolation formula:

where we have used the notations shown in Fig.

Different slope limiters are used in order to calculate the ∇v i from the gradients associated with the triangles K i and K j sharing the face S 1 S 2 :

Appendix B. Calculation of the time step

Flexibility of the code lets us to choose the most appropriate time integration schemes, with implicit or explicit temporal time marching.

The explicit schemes require an additional stability criteria (Courant-Friedrichs-Lewy condition) over the time step of the simulation, ∆t:

∆t

where CF L = O(1) and the time step ∆t stab is calculated as the harmonic mean of the convection time step ∆t conv and diffusion one ∆t dif f as:

and ∆x the minimum mesh size of the fluid domain. Thus, the time step is calculated as:

where ∆t stab is calculated by TrioCFD before each time step and ∆t max is set by the user. We can see that the restriction on the stability time step is due to the diffusion time step, especially in the case of a fine mesh (quadratic term in ∆x). Thus, an implicit treatment of the diffusion term can considerably reduce the time of calculation. To overcome this problem, one can use an implicit scheme, as in our case, or can implicitly treat the diffusion term, at each time step using a conjugate gradient approach. The time step is then calculated only from the convection time step (when an implicit scheme is used, the time step is calculated as if the convection is completely explicit even if it is semi-implicit).

Appendix C. Study 1. Further analyses

In this Appendix, we report further results obtained for the case study 1.

Firstly, we check the scale invariance of the numerical results, considering two self similar configurations in terms of dimensionless numbers: ε = 2, Sk = 10 4 and KC = 10 -2 . The time evolution of the dimensionless fluid forces is represented in Fig. C.17, for both cases. As expected, the numerical predictions are strictly identical for both configurations, confirming that the dimensionless added coefficients are functions of Sk, ε and KC.

Appendix D. Study 1. Tables of comparison numerics versus theory

In this Appendix, we report the tables of comparison between the theoretical and numerical values of the fluid added coefficients for the case study 1, with two coaxial vibrating cylinders. The relative deviation ι is also reported in the tables. The notations ET and AT refer to the exact and asymptotic theories.

Appendix E. Effect of the mesh size and time step on the fluid added coefficients

In this Appendix, we report the mesh converge analysis performed for both cases study. For the case study 1, in Tabs. E.6, E.7, E.8, E.9 and E.10 we clearly show a convergence of the mass coefficients as l c and lc f ine are changed. The convergence of the damping coefficients is less obvious, especially for high values of Sk. Physically, this is related to the thickness of the boundary layer, which tends to zero as Sk increases. It follows that a finer mesh is required close to a cylinder boundary to account for the thickness of the boundary layer and obtain an accurate estimation of the damping terms. Moreover, all the simulations are performed with CF L = 1 to ensure maximum precision, except for the simulations with a low value of the Stokes number, Sk ∈ {10 1 , 10 2 } and ε = 2 where the CF L is progressively increased from 1 to 10, without influencing the quality of the results. Following the numerical approach of the case study 1, we perform a mesh sensitivity analysis on the case study 2, see Tab. E.11.


Let (X 1 , . . . , X n ) be i.i.d. r.v.

Theorem 1 (S. [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF] Under mild conditions and for a sufficiently large u, N n can be approximated by a non-homogeneous Poisson process N of intensity measure Λ with parameters θ = (µ, σ, ξ) ∈ R × R + × R, such that for all x > u, Under mild conditions and for a sufficiently large u, N n can be approximated by a non-homogeneous Poisson process N of intensity measure Λ with parameters θ = (µ, σ, ξ) ∈ R × R + × R, such that for all x > u,

Λ(I x ) = +∞ x λ(t)dt = 1 + ξ x-µ σ -1 ξ + if ξ = 0 , exp(-x-µ σ ) if ξ = 0.
Λ(I x ) = +∞ x λ(t)dt = 1 + ξ x-µ σ -1 ξ + if ξ = 0 , exp(-x-µ σ ) if ξ = 0.
Likelihood: p PP (x | µ, σ, ξ)= exp -m 1 + ξ u -µ σ -1 ξ σ -n n i=1 1 + ξ x i -µ σ -1 ξ -1
The Bayesian Point of view

Inference on θ = (µ, σ, ξ) using Bayesian algorithms such as MCMC.

The main advantages are (S. G. [START_REF] Coles | Bayesian Methods in Extreme Value Modelling: A Review and New Developments[END_REF]):

• Consideration of expert information with informative prior,

• Can be used in any case, even when the likelihood is irregular (for ξ < -1),

• Access to the posterior predictive distribution:

p(x | x) = p(x | θ) randomness of future obs. p(θ | x)
parameter uncertainty dθ.

The Bayesian Point of view

Inference on θ = (µ, σ, ξ) using Bayesian algorithms such as MCMC.

The main advantages are (S. G. [START_REF] Coles | Bayesian Methods in Extreme Value Modelling: A Review and New Developments[END_REF]):

• Consideration of expert information with informative prior,

• Can be used in any case, even when the likelihood is irregular (for ξ < -1),

• Access to the posterior predictive distribution:

p(x | x) = p(x | θ)
randomness of future obs.

p(θ | x)

parameter uncertainty dθ.

Computational Bayesian issues: prior elicitation, reparameterization, etc.

The uniform prior on (µ, log σ, ξ) leads to a proper posterior for the Poisson process model for extremes as soon as n ≥ 3.

Theorem 3 (Moins et al.)

Jeffreys prior for a Poisson process for extremes with parameters θ = (µ, σ, ξ) exists provided ξ > -1/2, and can be written as

π J,PP (µ, σ, ξ) ∝ 1 + ξ u-µ σ -3 2ξ -1 σ 2 (1 + ξ)(1 + 2ξ) 1 2
.

It is improper and leads to a proper posterior as soon as ξ > -1/2.

Alternative parameterisation: instead of (µ, σ, ξ), run MCMC on an orthogonal parameterisation suggested for another purpose by [START_REF] Chavez-Demoulin | Generalized additive modelling of sample extremes[END_REF]:

(r , ν, ξ) = m 1 + ξ u -µ σ -1/ξ , (1 + ξ)(σ + ξ(u -µ)), ξ .
r : intensity of PP, expected number of exceedances. This ensures full orthogonality, which means that I (r , ν, ξ) is diagonal.

Jeffreys prior:

π J,PP (r , ν, ξ) ∝ r 1/2 ν(1 + ξ)(1 + 2ξ) 1/2 .
Alternative parameterisation: instead of (µ, σ, ξ), run MCMC on an orthogonal parameterisation suggested for another purpose by [START_REF] Chavez-Demoulin | Generalized additive modelling of sample extremes[END_REF]

: (r , ν, ξ) = m 1 + ξ u -µ σ -1/ξ , (1 + ξ)(σ + ξ(u -µ)), ξ .
r : intensity of PP, expected number of exceedances. This ensures full orthogonality, which means that I (r , ν, ξ) is diagonal.

Jeffreys prior:

π J,PP (r , ν, ξ) ∝ r 1/2 ν(1 + ξ)(1 + 2ξ) 1/2 .

Our experiments show better convergence with orthogonal parameters.

Theorem 2 (Moins et al.)

The uniform prior on (µ, log σ, ξ) leads to a proper posterior for the Poisson process model for extremes as soon as n ≥ 3.

Alternative parameterisation: instead of (µ, σ, ξ), run MCMC on an orthogonal parameterisation suggested for another purpose by [START_REF] Chavez-Demoulin | Generalized additive modelling of sample extremes[END_REF]

r : intensity of PP, expected number of exceedances. This ensures full orthogonality, which means that I (r , ν, ξ) is diagonal.