Numerical Investigation of Vortex Induced Vibration for Submerged Floating Tunnel under Different Reynolds Numbers

: A 2D numerical model was established to investigate vortex induced vibration (VIV) for submerged floating tunnel (SFT) by solving incompressible viscous Reynolds average Navier-Stokes equations in the frame of Abitrary Lagrangian Eulerian (ALE). The numerical model was closed by solving SST k - ω turbulence model. The present numerical model was firstly validated by comparing with published experimental data, and the comparison shows that good achievement is obtained. Then, the numerical model is used to investigate VIV for SFT under current. In the simulation, the SFT was allowed to oscillate in cross flow direction only under the constraint of spring and damping. The force coefficients and motion of SFT were obtained under different reduced velocity. Further research showed that Reynolds number has not only a great influence on the vibration amplitude and ‘lock-in’ region, but also on the force coefficients on of the SFT. A large Reynolds number results in a relatively small ‘lock-in’ region and force coefficient.


Introduction
Submerged floating tunnel (SFT) is a new type of traffic structure that crosses the strait, bay and lake. It is usually suspended more than 30 m below the water. The SFT has a large internal space, which is sufficient to meet the requirements of roads and even railways. For some fjords with harsh natural conditions, due to environmental conditions and technical constraints, traditional spanning methods (such as: cross-sea bridges, immersed tunnels) are not feasible, and SFT offers the possibility of crossing. Since the SFT is always in a deep-water depth, whose location is more the half wave length of normal wave, the normal wave has less influence on it. In addition, due to severe natural environment in the fjord, the flow velocity is usually fast, which has a greater influence on the SFT.
In view of the coupled analysis of flow and SFT, many scholars have done the relevant researches. Mai [1] considered the fluid-structure interaction effect, and studied the effects of surface velocity, tunnel section form and support form on the dynamic response of the SFT. It is found that the surface velocity would significantly affect the response displacement of the SFT, but it did not affect the stress distribution along the axial direction. Wang [2] analysed the variation law of the load on the SFT structure under the lateral lift force of the flow with the submerged depth, water depth, flow velocity and section size. Long [3] studied the dynamic response of SFT in different Buoyancy Weight Ratios (BWRs), and proposed the optimal range of BWRs for SFT under flow loading. The empirical lift formula based on Morison's formula was used in the above studies when discussing the flow load, but the effect of SFT on the flow field was not considered in detail. In order to find the optimal section of SFT, Luo et al. [4,5] compared the flow field distribution and force acting on the fixed SFT with different cross-section forms by large eddy numerical simulation. It was found that the ear-shaped ( Figure 1) SFT structure had smaller lift coefficient and drag coefficient, which was the most reasonable cross-section shape, followed by circle, ellipse, hexagon and rectangular. Since the SFT is suspended in the sea via the mooring system, it will move under the flow action. The periodically varying lift makes the SFT with elastically support vibrate perpendicular to the inlet flow direction, that is, 'Vortex-Induced Vibration' (VIV). When the vortex shedding frequency is close to the natural frequency of the structure, the phenomenon of resonance or lock-in occurs, which reflects the complex interaction between the fluid and the structure. Therefore, when frequency lockin occurs in the SFT under the flow action, the fatigue damage of the structure will be significantly increased, which will have a negative impact on the safety of the project. Many experimental researches and numerical simulations on the VIV problem were carried out, such as Morse and Williamson [6], Govardhan and Williamson [7], Yan et al. [8], Luo et al. [9], Zheng et al. [10]. VIV experiment of rigid cylinder with elastically support under wind load was successfully studied by Feng [11]. Williamson and Khalak [12,13] and Govardhan et al. [14] performed VIV experiments on rigid columns with low mass ratio and elastic support in the wave tank, which became the verification test for many subsequent numerical simulations. Lu and Dalton [15] numerically studied the cylindrical VIV problem with cross-flow motion in the case of Reynolds number Re = 13,000, and the model used large eddy simulation to close the turbulence equation. Dong and Karniadakis [16] used direct numerical simulation (DNS) to study the force vibration of cylinders with a Reynolds number of 10,000. For the vibration analysis of SFT, Ge et al. [17] and Kang et al. [18] applied the von der Pol equation to simulate the VIV of the SFT in flow, and studied the influence of the spacing of the anchor chain on the vibration amplitude of the tunnel. Su and Sun [19] utilised the wake oscillation model to simulate the VIV of the SFT. It was found that the vortex-induced vibration resonance occurred and the axial stress of the structure increased significantly. Chen et al. [20] proposed a simplified theoretical model for vibration analysis of the coupled SFT tube-cable system under wave and current.
In general, the empirical formulas are employed in most research about the VIV of SFT, and the study based on computational fluid dynamic (CFD) are limited on a specific Reynolds number. For this reason, the investigation of VIV based on CFD in a serious of Reynolds number under a specific engineering background should be proceeded here. The vibration phenomenon of SFT under different Reynolds number is different because the size of the structure and flow velocity vary greatly in different situations. Based on the above background, this study aims to analyse the coupled motion of flow and SFT in the range of Reynolds number from 1000 to 100,000, and study the influence of different Reynolds numbers on the vibration of SFT as well as force coefficient under the flow action, so as to provide reference value for practical engineering.

Numerical Model
For underwater SFT, the flow around the structure is usually turbulent due to the large scale of the structure itself and the relatively fast flow velocity. At present, the simulation of turbulence can be approximated by direct numerical simulation (DNS) or by using a suitable turbulence model, and the turbulence model is used to simulate the problem in this paper. At the same time, since the length of the SFT structure is much larger than the section size, we can simulate this problem as a twodimensional (2D) flow-structure interaction problem approximately. Although it is well known that the flow is indeed three-dimensional (3D) effect for large Reynolds number, however, 3D simulation cost a lot of computational resources. Therefore, a 2D numerical model was adopted in this study. Although 2D numerical model will overpredict the numerical results, it still can reveal the relationship between reduced velocity, vibration amplitude and force coefficient. In addition, many other researchers have also adopted 2D numerical model to solve similar problems (Lu and Dalton [15], Dong and Karniadakis [16]).

Governing Equation and Turbulence Model
The two-dimensional incompressible Reynolds-Averaged Navier-Stokes equations are adopted to describe the turbulence flow of incompressible viscous fluid. The governing equations in the Arbitrary Lagrangian-Eulerian (ALE) frame can be written as (Liu et al. [21]) where x1 = x, x2 = y are the horizontal and vertical coordinates, respectively, ui is the fluid velocity in the xi-direction, t is the time, m j u is the velocity of moving grid in the xj-direction, p is the pressure, ρ is the fluid density, υ is the kinematic viscosity of the fluid, Sij is the mean strain rate tensor with ( ) The Reynolds stress term in Equation (1) reads can be expressed as where υt is the turbulent eddy viscosity, k is the turbulence kinetic energy and δij is the Kronecker operator.
In order to close the governing equations, the Shear Stress Transport (SST) k-ω turbulence model (Menter [22]; Menter et al. [23]) is adopted. The parameters in the equation have been widely accepted and successfully applied, which has shown good performance in simulating the boundary layer flows with significant adverse pressure gradient. The governing equation of the SST k-ω turbulence model can be written as follows: where Pk is the production of turbulent kinetic energy and the related parameters in Equations (4) and (5) are calculated as follows where Ω is the absolute value of vorticity, y* is the distance to the nearest solid wall, and the parameters F2 and Dkω are By using the blending function F1, the following parameters can be calculated, i.e., The model constant in the SST k-ω model are listed in Table 1. When the flow field and the pressure field are obtained, the fluid force acting on the structure can be obtained by integrating the surface pressure and the viscous shear force over the body surface. The dimensionless drag coefficient CD and the lift coefficient CL are respectively

Motion Response of SFT
For the vibration problem of the SFT under the flow action, since it is necessary to ensure the anchor cable is always in the elastic range, the whole system can be simplified as a mass-dampingspring system. In this paper, only the vibration response of the SFT in cross flow direction is considered, and its motion equation can be expressed as follows.
where ξ is the damping ratio of structure, fn is the natural frequency and m* is mass ratio.
The following dimensionless relationship can be defined further At the same time according to the definition of lift coefficient, Fy = 0.5ρU 2 DCL, Equation (9) can be transformed in dimensionless form ( ) By introducing the definition of the reduced velocity Ur = U/fnD, the motion equation of the SFT can also be expressed as a dimensionless equation in the form of reduced velocity.
The lift coefficient at the right end of the above formula has been given before, and the dynamic response of the SFT can be calculated by the above formula.

Calculation Model and Boundary Conditions
The calculation model and boundary conditions are shown in Figure 2. Let the origin of the coordinate be at the initial centre of the cylinder, and the dimensionless cylinder diameter D = 1. The dimensionless velocity u = 1, v = 0 are set up at inlet. Symmetrical boundary conditions are applied at side wall ∂u/∂y = 0, v = 0. The outlet velocity boundary condition is ∂ui/∂t + c∂ui/∂xi = 0, where c is local average flow velocity. non-slip boundary conditions applied to the cylindrical surface u = dx/dt, v = dy/dt. In the calculation, in order to ensure the first layer of the grids is in the viscous boundary layer, the distance between the surface of the circular cylinder and the first layer of the grids is less than 0.05% D. In addition, the choice of boundary layer can be employed by the method of Palm et al. [24]. In the calculation, the pressure at outlet p = 0 and the pressure boundary condition ∂p/∂n = 0 is applied at other boundary conditions, n is unit normal vector pointing out the fluid domain. At initial time, the velocity and pressure in the fluid domain are set up to zero, i.e., the initial velocity field satisfies the continuous equation.

Numerical Dispersion and Grid Update
The convection-diffusion equation is solved by streamline upwind/petrov-galerkin finite element method in this paper (Brooks and Hughes [25]), and this method has been applied in the solution of impressive flow problem successfully (Mochida and Murakami [26], Kim et al. [27], Guilmineau and Queutey [28]). The distribution method is utilised in the time integration of the momentum equation. First, the pressure term is neglected, and the intermediate velocity of convection and diffusion term is considered; then the pressure equation is solved to calculate the pressure of the next time step; finally, the pressure gradient term is considered to correct the flow field. Streamline upwind method is used to predict the velocity.
The Newmark-β method is used to solve the motion equation of the structure. Given the displacement, velocity and acceleration at an initial moment, the time step Δt, the parameters β and γ are selected, and then the equivalent stiffness is formed. The effective load at time t + Δt is obtained, and the displacement at the time t + Δt can be solved. For time advance, according to CFL (Courant-Friedrichs-Lewy) conditions, the following dynamic time steps are adopted: where Sc is the mesh area, ue is the flow velocity at grid centre, min indicates the minimum in the computational domain, Cs is the safe coefficient, Cs = 0.2. Due to the reciprocating motion of the SFT under the flow action, the dynamic grid method based on ALE is used to simulate the fluid-structure coupling problem. In this paper, the mesh in the computational domain is assumed to be an elastic one, as shown in Figure 3, in order to achieve the purpose of adapting to the movement of the grid boundary nodes and internal nodes. The motion and deformation of the mesh can be obtained by solving the governing equation of linear elastodynamics (Johnson and Tezduyar [29]). The mesh updating method make the displacement of the mesh nodes more uniform and improve the stability of numerical calculation. In addition, the possibility of mesh distortion can be reduced by controlling the elastic modulus of the computing element. Specifically, the balance length of the mesh is equal to the length of the mesh itself at the initial moment. When the two ends of the joint move relative to each other, the mesh will be stretched or compressed, correspondingly. The mesh still satisfies Hooke's law, so the total force vector of any node i is ( ) where Fi is the total force vector on node i, αij is the mesh stiffness between the nodes, υi represents the number of nodes connected to node i, j = [i, υi], δi and δj are displacement vector of node i and j.
In order to avoid collisions of mesh nodes, the following expression is usually used to calculate the mesh stiffness where xi, xj are the position vectors of nodes i and j, that is, the value of αij is considered to be the reciprocal of the side length. This mesh updating method has also been widely used in the VIV research (Tang et al. [30], Lu et al. [31]).

Model Validation
In order to obtain reliable numerical results, this paper firstly uses the free vibration cylinder problem under Re = 30,000 and Ur = 6.00 as an example to verify the mesh convergence of the numerical model. Four different meshes are considered in Table 2. From Table 2, it can be seen that the numerical results under the four meshes are very close, indicating that the numerical results have converged under the current grid density. Considering the computational efficiency, the latter numerical calculation takes mesh 3 (Figure 4) as the benchmark. In the table, Ymax denotes the maximum vibration amplitude of the cylinder, Dismin denotes the minimum distance between the circular cylinder and the first layer gird, M D C is the mean drag coefficient, RMS L C is root mean square (RMS) of lift coefficient.   [32] and the numerical results of Dong and Karmiadakis [16], Zhao et al. [33], Song [34], which are shown in Table 3.  Table 3  Reynolds number. The problem of flow around a fixed cylinder is actually only a unilateral hydrodynamic calculation problem. It does not involve the motion response calculation of the circular cylinder itself, so it is not a true fluid-structure coupling problem. The vibration of cylinder with elastic support under the flow action involves fluid-structure coupling problem. There have been many experimental results on the VIV of rigid cylinders with elastic supports under the flow action, such as Khalak and Williamson [12], which have done systematic experimental studies. In order to facilitate comparison with the experimental results, the same calculation parameters as in the tests of Khalak and Williamson [12] were used. The Reynolds number Re = 12,000, the mass ratio m* = 2.4, and the mass damping ratio m*ξ = 0.013, and the maximum dimensionless vibration amplitudes versus reduced velocity are calculated in the paper.
From the comparison results in Figure 5, it can be seen that the maximum dimensionless vibration amplitude is close to 1 and the 'lock-in' region is from Ur = 4.0~10.0. In addition, the upper branch and lower branch can also be described in the numerical model clearly as shown in the experiment. Therefore, the results in the numerical model agree well with the experimental results of Khalak and Williamson [12]. It is illustrated that the model established in this paper can be used to investigate the fluid-structure coupling problem with high Reynolds number.

Example Analysis
Based on the above numerical model, this paper calculates the motion of the SFT under different constraint stiffness and different Reynolds numbers. The Reynolds number is calculated from 1000 to 100,000. Firstly, the time history curves of the cross-flow direction under the conditions of Reynolds number 50,000, mass ratio m* = 2.5, damp ξ = 0.007 and reduced velocity Ur = 2.0, 5.0 and 12.0, respectively, are introduced, as shown in Figure 6. Then, the Fast Fourier Transform (FFT) is utilised in the lift force coefficient time history, and the result is displayed in Figure 7.
From Figure 5, it can be found that when reduced velocity are 2.0 and 12.0, the vibration amplitude in cross flow direction of the SFT is small, and when the reduced velocity is 5.0, the vibration amplitude is large. Subsequently, the flow field is analysed in the case of larger and smaller vibration amplitudes, as shown in Figures 8 and 9.  From Figure 6, it can be seen that the vortex shedding frequency is about 0.23 Hz. According to the definition of Ur = U/fnD, the natural frequencies (fn) of SFT in the case here are 0.5 Hz, 0.2 Hz and 0.083 Hz, respectively. Therefore, it can be seen that when the reduced velocity of the SFT are 2.0 and 12.0, the vortex shedding frequency is far away from the natural frequency, and the VIV of the SFT is not obvious. The wake pattern is shown in Figure 9, and the wake shape is also regular in one vibration period. However, when the reduced velocity is 5.0 by adjusting the spring stiffness, the frequency of vortex shedding is close to the natural frequency. Under the action of flow lift force, a large VIV phenomenon occurs in the SFT. The vibration amplitude is larger, even reaching 0.7 times of the outer diameter of the SFT, whose wake pattern is shown in Figure 8. It can be seen that the wake has a long 'tail' after being separated, and the wake shape is irregular. Next, we compare the coupling effect of the flow and the SFT versus different reduced velocity under this Reynolds number, and the maximum dimensionless vibration amplitude of the SFT has been statistically obtained, as shown in Figure 10.   It can be seen from the calculation results that when the reduced velocity is from 4.0 to 10.0, the structure is 'locked' under the flow action, while at other reduced velocity, the vibration amplitude of the structure is small. Then the VIV of the SFT under different Reynolds numbers are compared, and the calculation results are shown in Figure 11. From the results of VIV of SFT under different Reynolds numbers, it can be seen that Reynolds number has not only a great influence on the vibration amplitude, but also on the 'lock-in' region. In general, the lower the Reynolds number is, the larger the amplitude is out of the 'lock-in' region. The minimum vibration amplitude in the 'lockin' region is about 0.4D, while the maximum vibration amplitude in the 'lock-in' region can reach to 0.8D. It can also be seen that larger Reynolds number leads to narrow 'lock-in' region in Figure 10. Therefore, when the size of the SFT is small or the inlet flow velocity is slow, VIV 'lock-in' phenomenon is more likely to occur for SFT because of lower Reynolds number. Figure 11. Vibration amplitude of SFT versus natural frequency in different Reynolds numbers. Figure 12 is the force coefficient on SFT versus reduced velocity at different Reynolds numbers. It can be seen that reduced velocity has a greater influence on the mean drag coefficient and RMS of lift coefficient. When the Reynolds number is low, the mean drag coefficient and RMS of lift coefficient of the SFT are relatively large from 1000 to 10,000. As the Reynolds number increases, the mean drag coefficient and RMS of lift coefficient become smaller. Therefore, when the size of the SFT is small or the inlet flow velocity action on the structure is slow, the force coefficient is large, and when the size is large or the flow velocity is fast, the mean drag coefficient and lift force coefficient of the structure are small.

Conclusions
Based on the FEM solution of incompressible viscous Reynolds average Navier-Stokes equations, combining the frame of Abitrary Lagrangian Eulerian, through accurate computational fluid dynamic numerical simulation, the vortex-induced vibration problems of submerged floating tunnel with different Reynolds numbers are studied. Main conclusions are as follows: Firstly, the analysis of uniform flow and fixed single cylinder proves the accuracy of the model in the case of high Reynolds number. Then, by simulating the vortex-induced vibration of a cylinder at high Reynolds number and comparing with other scholars' experimental results, it is proved that the model established in this paper can be used to study the fluid-structure coupling problem at high Reynolds number.
Through the research of vortex-induced vibration of SFT under the flow action, the force coefficient and motion of the SFT versus different reduced velocity at different Reynolds numbers are analysed. The results show that the Reynolds number has not only a great influence on the vibration amplitude and 'lock-in' region, but also on the force coefficient on the SFT. When the Reynolds number is low, the 'lock-in' region, the mean drag coefficient and RMS of lift coefficient of the SFT are relatively large. As the Reynolds number increases, the 'lock-in' region, the mean drag coefficient and RMS of lift coefficient become smaller. Therefore, when the size of the SFT is small or the flow velocity action on the structure is slow, the force coefficient and 'lock-in' region are relatively large, while when the size is large or the flow velocity is fast, the force coefficient and 'lock-in' region are relatively small.