Numerical simulation of an elementary Vortex-Induced-Vibration problem by using fully-coupled fluid solid system computation

Numerical simulation of Vortex-Induced-Vibrations (VIV) of a rigid circular elastically-mounted cylinder submitted to a fluid cross-flow has been extensively studied over the past decades, both experimentally and numerically, because of its theoretical and practical interest for understanding Flow-Induced-Vibrations (FIV) problems. In this context, the present article aims to expose a numerical study based on fully-coupled fluid-solid computations compared to previously published work [34], [36]. The computational procedure relies on a partitioned method ensuring the coupling between fluid and structure solvers. The fluid solver involves a moving mesh formulation for simulation of the fluid structure interface motion. Energy exchanges between fluid and solid models are ensured through convenient numerical schemes. The present study is devoted to a low Reynolds number configuration. Cylinder motion magnitude, hydrodynamic forces, oscillation frequency and fluid vortex shedding modes are investigated and the “lock-in” phenomenon is reproduced numerically. These numerical results are proposed for code validation purposes before investigating larger industrial applications such as configurations involving tube arrays under cross-flows [4].


INTRODUCTION
Studies on Vortex Induced Vibrations (VIV) encountered a large interest over the past decades because this is a phenomenon appearing in many fields: riser tubes on off-shore oil rigs [28], [29], steam exchangers [23] or bridge cables are VIV-prone.It may cause severe damages when vibrations are repetitive and involve large-magnitude motion.VIV are still a research challenge because of the lack of understanding of certain observable physical processes, and because of the large computational resources and time required for numerical simulations.
In the present work, we are concerned with the numerical simulation of the cross oscillations of an elastically-mounted rigid cylinder at a low Reynolds number without any structural damping as previously studied by [34], [36] providing reference data to be used in the present work.The method to be proposed here consists in using a fully-coupled fluidsolid system computation and applying it to the dynamical analysis of a mechanical system subjected to a fluid flow.The fluid solver to be used is the CFD code: Code_Saturne, developed by EDF R&D, based on a collocated finite volume approach.Momentum equations are solved considering an explicit mass flux.Velocity and pressure coupling is insured by a SIMPLEC prediction/correction method with outer iterations [3].Within the framework of this work, a resolution of structure equations is directly performed in Code_Saturne.Both fluid and solid resolutions are numerically coupled thanks to the modeling of a fluid-structure interface.The main purpose of this paper is to compare our numerical results to those obtained by previous authors with other codes or techniques [14], [20], [31], [24] for simulating the lock-in.Under specific conditions, the vortex shedding frequency can be tuned to the natural frequency of the oscillating cylinder.This phenomenon is an important case in the framework of code validation because of its specificity in the fluid-structure interaction field.The main result of the present work is that "lock-in" * (or "synchronization") patterns and the influence of the mass ratio on the amplitude and frequency response are simulated with a good accuracy.
The paper is organized as follows: main VIV characteristics for rigid circular cylinder are recalled in the first section, then basic principles and choices of numerical approaches are presented.Finally, in the third part numerical results are discussed and compared to [34] and [36] regarding cylinder displacement amplitude, frequency, hydrodynamic forces, vortex shedding modes, and influence of mass ratio on these variables.The capability of the method to reproduce the lock-in phenomenon and its performances are evaluated, giving positive results for further applications involving VIV and more generally Flow-Induced-Vibrations (FIV).

BACKGROUND ON VIV OF A RIGID ELASTICALLY-MOUNTED CYLINDER UNDER CROSS FLOW
Figure 1 shows flow past a cylinder in a (x, y) plane with a uniform flow velocity.Flow upstream velocity is denoted U ∞ .For three-dimensional configurations, it is necessary to consider the third z-direction, in which a lot of complex phenomena interact with twodimensional ones, especially when the Reynolds number increases.
*Definition of "lock-in" as a simple matching of the vortex shedding frequency and the cylinder natural frequency refers to a rather crude description of the process.A thorough discussion on the "true" definition of lock-in is developed in [13] or [21].Fluid and solid parameters are described in Table 1: From a physical point of view, it is interesting to recall dimensionless parameters commonly used to describe the flow and its interactions with a structure (Table 2): the Reynolds number R e which represents the ratio between inertia and viscous forces, the Strouhal number S t which characterizes the reduced frequency of the system, describing the vortex shedding frequency (f s ), and drag and lift coefficients (resp.C D and C L ), respectively calculated from the evaluation of in-line and transverse forces (F x and F y ) on the cylinder induced by the flow.
Considering the cylinder oscillating in the cross direction because of vortex shedding phenomenon and assuming the cylinder oscillates around its mean position in the presence of flow, the cylinder motion equation is: (1) where F Y the transverse force induced by the flow is deduced from the Navier-Stokes equations.The variable Y designates the cylinder displacement due to vortex shedding.A dotted letter represents a time derivation.This dimensional equation can be re-written thanks to dimensionless parameters, as follow: (2) Reduced variables are gathered in Table 3 f N is the proper frequency of the cylinder in quiescent fluid [13], and are respectively the cylinder mass and the added mass (equal to the fluid displaced mass in the case of a simple Dimensionless parameters of a fluid-structure system. cylinder).Maximum dimensionless amplitude A* is deduced from an average of dimensionless displacement maxima.Equation (2) shows that the amplitude of cylinder's vibration a priori depends on mass and damping ratios.The problem can be described by using only one coupled mass-damping parameter (m − ζ) if , (m + C a ) ζ >0.006where the added-mass coefficient Ca represents the ratio m/m d [21].Amplitude response depends of the mass-damping parameter: For high m* −ζ: amplitude response drawn versus U* is a two-branch response, respectively called "initial excitation branch" and "lower branch".For low m* −ζ: there are three distinct branches: initial and lower branches are still present, and there is another one, the so-called "upper branch".Details on this cylinder response (Figure 2) show a hysteretic phenomenon ("H", Figure 2) between initial and lower branch in both cases: depending on the evolution of U*, increasing or decreasing, the obtained response is not the same.Moreover, in case of low m* − ζ, during the transition from upper to lower branch, an intermittent switching ("I", Figure 2) in amplitude values is observed.Maximal amplitude is governed by the coupled mass-damping parameter for a large range of mass and damping ratios [21], while when m* −ζ remains constant, the mass ratio variation does mainly control the width of U* range where the system is in "synchronization" regime.Therefore, the higher the mass ratio m*, the more significant the tendency to develop a "lock-in" zone [36].
There is an m*-dependency of the U*-zone where the "lock-in" is observed: results of [36] for two different mass ratio at low m* −ζ and low Reynolds number are presented on Figure 3.
Many studies are made presenting variables as function of the reduced velocity.A more pertinent parameter is proposed in [36], which provides a unified scaling for system behavior, and which is defined when k* or ζ falls to zero.
More precisely, induced vibrations are observed in [36] even under extreme conditions, where mass, stiffness and damping ratios tend to zero.In case of low damping (under 0.01), it is possible to completely define the system thanks to a single dimensionless parameter

Symbol
Dimensionless ratio Definition taking into account both mass m* and stiffness κ * ratios.This parameter is called "effective stiffness" κ *.Considering a negligible damping and a purely sinusoidal supposed cylinder response and corresponding hydrodynamic forces, the governing equation of motion can be re-written in the following form: Figure 3. Influence of mass ratio m* on "lock-in" zone [36].problem by using fully-coupled fluid solid system computation And since: (4) one gets: (5) Thus, effective stiffness actually combines m * and κ * and for any choice (κ * , m * ) keeping κ * constant, a single response is consistent with the governing equation of motion accordingly to [36].Therefore, this combined parameter is used in [36] to describe a zero-or low damping system at a given value of Reynolds number.
Finally, it is interesting to briefly describe the wake behavior behind the cylinder.Oscillation modes clearly depend on the response regime (initial, upper or lower branch) [21].On the initial branch, vortex shedding occurs alternately from one side of the cylinder to the opposite: it is the so-called "2S" mode, for 2 single vortices shed per cycle (Figure 4a).The regime in the lower branch is characterized by a shedding of two pairs or vortices during a cycle ("2P-mode", Figure 4b).Concerning upper branch, conclusions are uncertain according to [21], while it is suggested in [32] that vortex shedding mode on upper branch corresponds to an intermittent switching between initial and lower branch, with a possible observation of a mixed regime 2S-2P.In [30], thanks to a numerical simulation in good agreement with experiments of the literature, it is noticed that when the "lock-in" phenomenon is possible, a second vortex is shed from one side or the other, in addition to the 2 single vortices.This mode is called "P+S" (Figure 4c).
The Reynolds number has although an influence on the cylinder response, especially at low m * − ζ.Simulations from [38] show that for a fixed value of reduced velocity, change vortex shedding modes may be observed.Table 4 gathers results obtained in [38] for a cylinder free to oscillate in both x and y directions, at U * = 4.92.C(2S) stands for "2S Coalescent": at the end of the wake, vortices merge.
The main purpose of the present paper is to investigate the propensity of the fully-coupled fluid solid system computation to reproduce these patterns, namely the "lock-in" behavior already observed in [36], and the effect of mass ratio on magnitude and frequency response.In the next section principles of fluid-structure coupling procedure are recalled, with a description of various coupling approaches, including time and space discretization choices.
The fluid moving mesh formulation allowing the management of mobile meshes is also introduced.Then, in the last part, results obtained thanks to these numerical investigations are discussed.The question of using quasi-implicit or explicit coupling scheme is asked, as well as the existing limits of the use of such a numerical partitioned procedure.

FLUID AND SOLID MODELS
Here governing fluid and structure equations are recalled in order to grasp their numerical resolution (see also [23], [25] and [37]).For incompressible viscous fluid flow, conservation equations may be written as follows: (6) where u is the fluid velocity, F represents volume forces and the stress tensor.Suffix f indicates fluid variables.The structure dynamic behavior in the cross direction without damping is governed by: (7) Boundary conditions are the following for the time, and the fluid velocity is imposed on Γ v part of the fluid domain boundary (Figure 5): Finally, two continuity conditions are to be satisfied at the interface: ( ) Table 4. Vortex shedding modes depending on the Reynolds number [38].

Reynolds number
Figure 5. Spatial boundary conditions on the studied domain.problem by using fully-coupled fluid solid system computation Fluid equations added to the pressure law yield to Navier-Stokes equations: (9) where v is the kinematics viscosity.Here, we consider a rigid cylinder: consequently, we do not use stress-displacement laws.
Exchanges through the interface between fluid and structure domains have to be precise, that is why the space discretization is a real challenge for the information conservation.Moreover, diffusion or overestimations all the way through the time must be avoided as far as possible, so the time coupling represents another problem, still under investigations [2], [23] and [25].

TIME MARCHING OF FLUID/STRUCTURE COUPLING
The main difficulties lie in the coupling scheme, which yields to a time-lag between fluid and structure equation resolutions.Moreover equation formulations are fundamentally different: fluid formulation is Eulerian while the structure one is Langrangian.Interaction between both fluid and structure domains is sketched up on Figure 6: A partitioned procedure has been chosen for the present study, because it is strategically convenient for using existing solvers for computation of both fluid and structure problems.A monolithic approach [6], [17] would imply a full new code  Figure 7. Time-lag implied by the partitioned approach (according to [18]).development.A partitioned scheme consists in the successive resolution of fluid and solid equations.
Figure 7 describes the time-lag previously mentioned between both resolutions.By using structure displacement at time t n+1 , calculations can be performed to obtain a fluid state on date t n+1 and fluid forces are thus deduced.Then, from these forces at t n+1 structure displacement at t n+1 is deduced thanks to the state at time t n .Thus, because of this time-lag, energy exchanges are not exactly conserved through the interface; that is why a good displacement predictor choice is crucial.
Among the partitioned schemes, we have to distinguish the weakly coupled from the strongly coupled ones.A scheme is said to be strongly coupled when equation ( 7) is exactly satisfied at each time-step, while a weakly coupled scheme does not enforce exactly equation (7).In the present work, we are interested in strongly coupled schemes, as far as possible.
Within the framework of our investigations, two partitioned coupled schemes are considered: explicit asynchronous [10], [33] and semi-implicit [1], [11], [16].The first one consists in considering an intermediary time step t n+1/2 .Effort prediction is simply identified to calculated fluid force: (10) where the superscript P indicates a predicted value.Interface position X f is given at time t n+1/2 thanks to the variables calculated at t n : (11) Force balance is respected with this first order scheme, as well as geometry.It is also possible to use a second order prediction for fluid forces of the form : (12) The second scheme, the so-called "semi-implicit" one, is built in order to improve explicit predictions: fluid sub-cycles are performed at each time step, so that the position prediction should converge.The first part of the algorithm consists in an initialisation step, and then predictions on the structure displacement are carried out thanks to the values of the previous sub-cycle [19] (Figure 8).

DATA PROJECTION AND SPATIAL DISCRETIZATION OF THE INTERFACE
Sometimes geometries of fluid and structure models do not coincide at the interface and it is necessary to use specific techniques to keep the information conserved through the interface.However, in the present work, the cylinder is modelled as a single point mass associated with a stiffness and information on forces and displacements is interpolated at this point.
The present chosen interfacial conservation method consists of convenient projections ensuring data transfers through the interface.Fluid forces are estimated at the centre of each cell face located on the interface.They are averaged and condensed at the mass point.Conversely the structure displacement is imposed at each node of the fluid mesh located on the moving boundary.In the present case we are dealing with a projection method combined with a condensation technique.This procedure is also possible with a beam modelling for the structure.
From a general point of view, the main difficulty of data projection techniques is limited to the possible incompatibility between fluid and structure interface modelling.It is necessary to distinguish two cases: 1) cases where interfaces are jointed, when the physical interface exactly corresponds to fluid and structure interfaces at each time step; 2) cases where they are disconnected or incompatible as explained in [11], [18] and [26].In the presence of incompatibilities several techniques are possible but master-slave techniques are suitable as far as fluid/structure interface problems are considered.An example is illustrated in [18].

FLUID MOVING MESH FORMULATION
The fluid domain discretization presents difficulties because structure displacement, and hence interface, have to be followed up.As previously mentioned, the fluid computational process is based on an Eulerian formulation, where a boundary displacement cannot be followed, while the structure resolution relies on a Langrangian formulation, which requires large computational times since the system encounters large-scale displacements.To get round this problem, a so-called "Arbitrary Lagrange Euler" (ALE) formulation is involved.It consists in giving a combined Eulerian and Lagrangian description to the system [9].Thus, a new calculation domain is created, following the system physical boundaries, even for high displacement problems.
In [39], [40] ALE implementation in fluid equations is studied for fluid-structure interaction problems.We now present the analytic formulation of fluid equations expressed on moving mesh.
Let Ω m be a material domain where fluid particles X are followed in their motions according to the Lagrangian description, Ω s being a spatial fixed domain with fixed spatial positions x and Ω ale a given domain whose arbitrary motion is different from the Ω m 's one (Figure 9).In Ω ale , the particular derivative is defined by: (13) represents the arbitrary referential velocity in Ω s .Practically, it is interesting to choose Ω ale in such a way that a part of its boundary is connected with the moving interface.
Thus, ξ and w respectively correspond to the mesh position and velocity: from now on, Navier-Stokes equations are solved in the referential Ω ale .So if the domain Ω s is moving, Eulerian formulation of Navier-Stokes equations is modified with the particular derivative definition in the Ω ale domain: (14) After variable change [7] [35], for incompressible flow, the new system becomes: (15) It is necessary to add the Geometric Conservation Law (GCL) [22] to this equation system.It ensures numerical conservation of physical quantities: this law states that the variation of an elementary volume during a time step ∆t must be equal to the volume covered by its boundary faces during the same period.This law is commonly used as a necessary condition to provide mesh stability.
The integral form of the GCL is ( 16): (16) In order to work with regular cells, the mesh velocity choice can either be based on the node displacement or on the velocity.In the case where the mesh displacement (or the mesh velocity) governing equation is linear, a first order approximation gives: where and are respectively the displacement at (n+1) th and n th time step and w i n the mesh velocity at time t = n∆t.
x x t i , ) Furthermore, the mesh actualization must be adapted to various cell forms and must avoid distortion of the mesh topology.Moreover, as it is called at each time step, the algorithm has to spare calculation time and to be as robust as possible.The grid updating is performed as follows: considering a domain with moving boundaries, each node of the mesh is linearly linked to two nodes of the boundary (Figure 10), thanks to the relation between positions noted X: (17) where α represents the parametric coordinates.
Thanks to this method, the boundary conditions are preserved.Finally, a general algorithm for a finite volume Navier-Stokes resolution with ALE method is presented.We first have to pay attention to the fact that physical magnitudes are calculated at cell centres, but the mesh deformation requires the information at nodes.The result is that an interpolation from mesh deformation to determine node displacements, respecting the GCL, must be performed.
Thus, the general ALE algorithm is given by: -Time loop ALE displacement boundary condition defined at the fluid-structure interface Transformation of these boundary conditions into velocity boundary conditions, GCL enforced Mesh deformation velocity calculation Mesh deformation thanks to velocity data and GCL Navier-Stokes equation resolution

NUMERICAL RESULTS AND DISCUSSION
By default, the time coupling scheme is the explicit asynchronous one and an unsteady method is adopted, as we are interested in the setting up of an established regime.Pressure calculation is performed thanks to a second order scheme.
In order to avoid 3D effects, investigations are performed at low Reynolds number (R e = 100).On Figure 11  In what follows several configurations are discussed.A large range of dimensionless parameter values is considered enabling both to improve the physical understanding of phenomena and to identify the borderlines of numerical scheme validity domain, especially for time integrators involved in the partitioned coupling procedure.

X X X
node boundary 1

VORTEX SHEDDING IN THE WAKE OF A FIXED RIGID CYLINDER UNDER CROSS FLOW
A time-step convergence is set up, so that the Courant-Friedrichs-Levy condition should be satisfied.Our aim is to be able to observe the vortex shedding phenomenon, so a flow perturbation is imposed at the beginning of the calculation with the fixed cylinder.Then, when the steady state has been reached, calculations with moving mesh are performed.Numerical results are presented in Table 5 and compared to reference data [14] [31] [34].
Here, respectively represent a mean of the drag performed at each time step, the mean on local maxima of the fluctuating lift signal, and the Root Mean Square (RMS) value of the fluctuating lift signal.Results are reasonably encouraging for the rest of the investigations.

VORTEX INDUCED VIBRATIONS AT MASS RATIO M* = 3.3
The cylinder is allowed to vibrate in the cross direction (y-direction), while the inlet flow velocity is uniform in the x-direction.Mechanical characteristics are chosen to be the same as those of [34] and [36].Thus, mass ratio is equal to 3.3, and the range of U* (or κ*, notation chosen to indicate k* eff in [36]) is chosen between U* = 3.00 and U* = 11.32 (or κ* from 1.47 to 17.64).Various dimensionless variables are plotted against effective stiffness κ*: the maximal amplitude A*, reduced frequency f* (where f* = fd/U ϱ with f the actual vibration frequency of the cylinder), and hydrodynamic coefficients C D,mean and C L,max .These results are gathered in Figures 11a) to 11d).
Good qualitative and quantitative reproduction of the "lock-in" phenomenon is obtained.It is encouraging for the use of coupled schemes in fluid-structure interaction, even with an explicit partitioned scheme.In a next part, the advantages of an implicit scheme will be discussed.Comparison between present simulations and [34], reference computation is [36], R e = 100.

C
A frequency analysis on maximal dimensionless amplitude reveals that the Khalak et al. [21] dimensionless frequency f* K&W (with f* K&W = f/f N ) according as κ* (or U*) be high or low is not the same.For low κ*, f* K&W is close to 2, and for high κ*, f* K&W is almost equal to 0.5.
For κ* around 2.44 (U* = 4.91), f* K&W is very close to unity: it confirms a "resonance" phenomenon, where the cylinder frequency is close to its proper frequency in water, all the more so since for each variable (A*, C D,mean and C L,max ), the biggest value is obtained for κ* 2.44.On the other hand, when U* 4.40, we observe a beating phenomenon (several frequency existence is revealed by the FFT of signals), also observed in [21], [34] and [36].This aspect is considered in [21]  At R e = 100, it is difficult to give a precise conclusion, but observing vorticity contours on Figure 13, it is possible to distinguish two different mode emission, to be compared to [34].Thus, when κ* is low or high (limits of the "lock-in" zone), the shedding mode is close to the ≅ ≅ Strouhal's one, (observed for a fixed cylinder) i.e.Von Kármán vortex street, while if κ* has a midway value, i.e. κ* = 2.44 ("lock-in"), vortices are not so stretched and the wake is clearly separated in two parts on both sides of the centerline.In spite of the fact that this study is led at low Reynolds number, the change in vortex shedding modes usually observed in the literature for higher Reynolds number is here perceptible.
Figure 14 provides a description of the link between the frequency behavior and the "synchronization": straight lines represent the Strouhal frequency versus reduced velocity, for R e = 100 and R e situated in [10 3 ,10 4 ], respectively calculated from the present Strouhal number and from empirical formula proposed in [31].The reduced frequency is f* K&W, plotted against reduced velocity.Data points that are situated close to Strouhal lines mean that for this range of reduced velocity, the Strouhal frequency (i.e.vortex shedding frequency) is governing the motion.On the other hand, data points situated just above f* = 1 depict the range of reduced velocity where the "lock-in" phenomenon is present.For higher Reynolds numbers (reference experimental data come from [21]), we can clearly distinguish three response branches.

MASS RATIO INFLUENCE
Now, we are interested in the propensity of the system to have a larger reduced velocity range where amplitude is the highest for a low mass ratio m*, than for a high mass ratio.At the same time, reduced frequency f* is most clearly "locked" close to unity for high mass ratio.Figure 15 gathers results obtained for m* = 1, m* = 33 and .The code reproduces well this mass-ratio dependency.
It is worth noting that for the case m* = 1, we have to choose the semi-implicit coupled scheme [7]: without sub-cycles, the explicit one is diverging.Comparing both synchronous explicit and semi-implicit schemes for κ* = 2.44 (U* = 4.91), we obtain that, as expected, the implicit scheme is more precise and an error between both amplitude response sets up.But this error remains low, and thanks to the good comparison with literature data, an explicit scheme is sufficient as long as convergence is assured.On the other hand, as soon as we want to work with a very low mass-ratio, neither of explicit or semi-implicit schemes is sufficient, even if a large number of sub-cycles is used.The mass ratio m* = 0.8 seems to be a limit of the use of such a coupled algorithm unless an additional stabilization procedure is involved.

CONCLUSION
Numerical simulations of the VIV phenomenon and particularly the frequency "lock-in" of a cylinder in a uniform cross flow, without damping and at low Reynolds number has been performed and discussed in the present study.Simulations have been performed thanks to a CFD code, using a partitioned coupled scheme for simulation of coupling with the structure motion.An ALE formulation has been chosen to solve fluid equations with a moving mesh.The aim of these investigations is to validate the computational process for this "lock-in" phenomenon, and for the mass ratio influence.Quantitative agreement has been highlighted, respectively for amplitude, frequency, hydrodynamic forces, and qualitative agreement has been observed for vortex shedding patterns.Moreover, mass ratio influence is quite well reproduced numerically.However, in the extreme case of very low mass ratio, a stabilized coupled schemes should be considered.
From a numerical point of view the objective is to show the interest of this numerical method compared to other strategies involved in the literature.Moreover from a physical point of view one of the objectives is to investigate physical mechanisms involved through suitable parametric studies in order to point of the influence of key dimensionless parameters characterizing the fluid structure interaction.This influence can be major from a physical point of view (example of the effective stiffness) and although from a numerical point of view (example of the mass ratio) because numerical properties like stability of numerical schemes like coupling time integrators depend on these parameters.
These results provide positive results confirming the capability of fully-coupled fluid solid system computation methods to build realistic models of fluid-structure interaction problems.

Figure 2 .
Figure 2. Scheme for both high and low m*-ζ types of amplitude response [21].

Figure 6 .
Figure 6.Fluid and structure domain interaction.

Figure 9 .
Figure 9. Coordinate mapping between the three domains. h fluid mesh involved in the present calculations are shown.

Figure 11 .
Figure 11.Fluid domain mesh.a) Global fluid mesh b) Zoom around the cylinder as a vortex shedding mode change.a) Dimensionless amplitude, b) Dimensionless frequency, c) Maximal lift coefficient, all variables are plotted versus κ* d) Mean value of the drag coefficient, all variables are plotted versus κ*.

Figure14.
Figure14.Representation of the frequency lock-in phenomenon.

Figure15.
Figure15.Dimensionless frequency and maximal amplitude for various mass ratios.

Table 1 .
Fluid solid system mechanical parameters

Table 5 .
Hydrodynamic forces and Strouhal number for a fixed cylinder