Numerical simulation of vortex induced vibrations of a circular cylinder: isothermal and heat transfer cases

In this paper, the Fluid Structure Interaction (FSI) problem - vortex induced vibration of a circular cylinder are studied using a numerical method. An accurate Harten Lax and van Leer with contact for artificial compressibility (HLLC-AC) Riemann solver have been used for flow computation. The Riemann solver is modified to incorporate arbitrarily Lagrangian-Eulerian (ALE) formulation in order to take care of mesh movement in the computation, where radial basis function is used for dynamically moving the mesh. Higher order accuracy over unstructured meshes is achieved using quadratic solution reconstruction based on solution dependent weighted least squares (SDWLS). Two flow scenario of 1 Degree of Freedom (DOF) oscillation cases: isothermal and with heat transfer are presented here. The results obtained by the present method are compared with the literature.


Introduction
Over the years, flow past circular cylinder has been studied [1][2][3][4][5][6][7][8][9][10][11] extensively due to its academic and practical significance. Periodic lift and drag forces are generated over the cylinder due to separation and vortex shedding leading to structural vibration also known as vortex induced vibrations (VIV). This fluid-structure interaction (FSI) may result in serious engineering problems such as fatigue failure of the off-shore riser and sub-sea pipelines, vibrations of tube of heat exchanger, vibration in structures such as chimney, high rise building, etc. Therefore, prediction of fluid forces on the cylinder has great importance from the point of view of structural design.
In the case of an oscillating cylinder, there exists phenomenon referred to 'lock-in' which was first observed experimentally by Bishop et al. [1]. Griffin and Ramberg [2] pointed out that lock-in augments the intensity of vortices. They observed experimentally that the horizontal spacing between the vortex streets is reduced as the amplitude of cylinder vibration is increased and the vertical spacing is independent of the changes in cylinder vibration amplitude. Khalak and Williamson [3] in their experiment found out the three different branches (initial, upper, and lower branch) for cylinder with very low mass damping ratio. Raghavan and Bernistas [4], in their experiment studied the transition from laminar to turbulent flow and observed that for high values of Reynolds number the vortex shedding mode is 2P. Soft lock-in range was investigated by Mittal and Kumar [5] at Re=325 using stabilized space-time finite formulation for various values of natural frequencies for the oscillator. They also carried out parametric study by varying oscillator frequency to observe its effect on in-line and transverse oscillator amplitudes. Zhou et al. [6] conducted a numerical study on one and two degree of freedom flow induced vibration of an elastic circular cylinder using discrete vortex method. They observed that the results obtained from one degree of freedom study are only in qualitative agreement with the two degree of freedom simulation. They also showed that maximum amplitude of vibration occurs when cylinder oscillation frequency matches with the natural frequency of fluid-structure system.
Studies on heat transfer enhancement through vortex induced vibration are very few in numbers. Fu and Tong [7] in their study numerically investigated the flow structures and heat transfer characteristics of heated transversely oscillating cylinder in a cross flow. They also studied the effect of Reynolds number, oscillating amplitude, oscillating speed on the flow structures and heat transfer from the cylinder. Zhang et al. [8] used immersed-boundary method to study heat transfer problems of flow over a circular cylinder. In this they also studied the effects of Reynolds number and different temperature boundary conditions on heat transfer through forced oscillating cylinder. The numerical work of Su et al. [9] is among the rare studies that have investigated flow-induced oscillations and the heat transfer process of a circular, which was installed in elastic supports. They used Ansys Workbench and CFX software and observed that vibration amplitude in the transverse direction is much higher than that in the stream wise direction. They inferred that when the maximum value of the cylinder velocity is much smaller than flow velocity, cylinder vibration does not enhance the heat transfer process. Wang et al. [10] used pressurebased high-order upwind scheme using PISO algorithm to numerically investigate the heat transfer from a vibrating cylinder under natural coupling. They concluded that the heat transfer was greater at high Reynolds number values. Baratchi et al. [11] used moving overset grids method to numerically investigate flow-induced vibrations of a heated elastically supported cylinder in a laminar flow at Re=200 and Pr=0.7. From the above literature study it can be seen that there is a scope as well as need to capture the flow physics of vortex induced vibration with heat transfer problem in a better manner. In this paper, the FSI problems: vortex induced vibration generated in the downstream of a circular cylinder for two different cases i.e. isothermal and with heat transfer are studied using a numerical method. An accurate Harten Lax and van Leer with contact for artificial compressibility (HLLC-AC) Riemann solver [12 -14] developed for solving incompressible flows in artificial compressibility formulation have been used for flow computation. The Riemann solver is modified to incorporate arbitrarily Lagrangian-Eulerian (ALE) [15] formulation in order to take care of mesh movement in the computation, where radial basis function [16] is used for dynamically moving the mesh. Higher order accuracy is achieved using quadratic solution reconstruction based on solution dependent weighted least squares (SDWLS) [17]. The results obtained by the present method are validated with those reported in the literature.

Arbitrarily Lagrangian-Eulerian (ALE) formulation
An artificial compressibility [18] methods with dual-time stepping the unsteady Navier-Stokes incompressible equations is modified here to take care of moving boundaries using the arbitrarily Lagrangian-Eulerian (ALE) formulation [15]. The integral form of the twodimensional governing equation in arbitrarily Lagrangian-Eulerian form can be written as where, = − and = − are the convective velocities in referential frame with and are the velocities of the moving grid in X and Y directions respectively. Note that, equation (1) Now equation (4) can be discretized in a very similar manner to that for unsteady Navier-Stokes equation for stationary boundary problem [12 -14]. The additional effort need to be added for ale flux vector. This additional term, ale flux, is nothing but the volumetric increment along the face and can be evaluated by considering the Geometric Conservations Law (GCL) [16]. The radial basis function [19]: Thin-Plate Spline (TPS) with global support is used here for mesh movement. The Thin Plate Spline with global support generates meshes of high quality after deformation along with the computational efficiency.
The convective fluxes at cell interface, that is, the stationary reference convective fluxes are evaluated using the Harten Lax and van Leer with contact for artificial compressibility (HLLC-AC) [12][13][14] Riemann solver. For viscous fluxes, a central differencing method based on Green-Gauss approach is used.

Convective flux computation using HLLC-AC upwind method
According to the HLLC-AC upwind method [12][13][14], which is a modified and developed from Harten Lax and van Leer with contact (HLLC) [20] upwind scheme used for the compressible flows, the normal components of the convective flux vector F ij = (E c l x + G c l y + H c l z ) at the interface between two cells 'i' and 'j' is evaluated based on a Riemann solver using three-wave structure (Fig. 1) as The convective flux is thus evaluated for each face of the cell using Eq. (5) once * and * are evaluated as per above Eq. (6).

Solution reconstruction: Evaluation of interface values of primitive variables (W L and W R )
The left and right state values of any solution variable, W at the interface between two adjacent cells i and j (Fig. 2) are found using Taylor series expansion of solution variables about the cell center and given by where, � and � are the values at the cell center i and cell center j, respectively; W L and W R are the interpolated values of the solution variables just at the left and right sides of the interface. The quantities ∇ � and ∇ � represent the gradients of � at left and right cell centers. The vectors r i , r j and r int are respectively the position vectors of the cell center i, cell center j and the interface between the cell i and j. The quantities H i and H j represent the Hessian matrices computed at the left and right cell centers, respectively. The order of accuracy of the scheme depends upon the number of terms kept on the right-hand side (RHS) of the expression (7). The third order scheme is obtained by keeping all the three terms, which requires an evaluation of both first and second order derivatives of with respect to space variable at the cell centers. Thus, quadratic reconstruction in three dimensions requires the calculation of nine derivatives at the cell centers are described in Fig. 2.

Solution dependent weighted least square (SDWLS) gradient calculation
The stencil consisting of vertex based neighbours, shown in Fig. 2 where Δx j , Δy j , Δz j are the distances along the x, y and zaxis between cell centers of vertex neighbour cell j (j =1, 2, ..., n) and the cell i at which the nine (three first order and six second order) spatial derivatives in (8) are being evaluated. This leads to an over determined system of equations as these derivatives were calculated by solving weighted normal equations of the form where W t = diag [W t,1 , ..., W t, n ] is the diagonal weight matrix. The expressions for weights [17,21] used here is Other weights given in [21] can also be used. For the three-dimensional case with a vertex-based stencil, depending on grids and flow conditions, the SDWLS results are found [22] to suffer for solution accuracy when the least squares problem is solved using weighted normal equations. It is well known that formation of normal equations can result in severe loss of accuracy [23]. Hence, SDWLS method above is suitably modified here for the use of QR algorithm. Since, the solution of the Eq. (11) is not always well behaved [22], for the reason mentioned earlier, it is not solved directly. Instead, the Eq. (11) is re-written as which can be considered as the normal equation for a plain least squares problem with rescaled matrix ( 1/2 ) and ( 1/2 ∆ ). The equivalent least squares system is a solved using QR algorithm, with 1/2 = � � ,1 … � , �, where expressions for w i are given in Eq. (12).
Bischot et al. [23] had introduced Lapack routine DGELS is used to solve the equivalent least squares problem. This routine uses the Householder transformation for the QR method to solve the full rank least squares problem.

Vortex induced vibration and its coupling with flow solver
An elastically supported cylinder of mass * and diameter * is kept in either a uniform flow field or flow field generated by mixed convection. The motion of the cylinder is governed by the following dimensionless equation: where i = 1 and 2 correspond to the X-and Y-direction respectively. , ̇, and ̈ are the non-dimensional displacement, velocity, and acceleration of the cylinder, respectively. The non-dimensional frequency is = * * * where the natural frequency of the system in vacuum is * = � * 4 2 * ; = * * is the damping ratio, where * = 2 * * (angular frequency * = 2 * ) is the critical damping of the system. Here, is the mass ratio, defined as = 4 * / * * 2 per unit length of the cylinder * and is the fluid density.
In Eq. (13), 1 ( = 2 * * * 2 * ) and 2 ( = 2 * * * 2 * ) are corresponding to drag and lift coefficient per unit length of the cylinder, respectively, where * and * are forces exerted on the cylinder in the X-and Y-direction, respectively. We implement a one-way coupling in which the fluid dynamic forces from the previous time step is used to find the solution of Eq. (14) in the current time step. The coupling scheme is stable for relatively large to moderate mass ratios and is computationally efficient.

Result and discussion
In the present paper, two cases for a flow induced vortex vibration over a circular cylinder allowed for 1 degree of freedom: isothermal and heat transfer cases are simulated. Figure 3 shows the domain considered, where outer boundary (inout) is located at 30 times diameter from the cylinder wall. Figure 4 shows the unstructured triangular grid having 23548 mesh element and 12037 nodes is selected based on the mesh convergence study. As per boundary conditions are concerned, at inlet, the velocities (x velocity, whereas y velocity is kept zero) are supplied according to the specified Reynolds number and pressure is interpolated from the interior. At outlet the atmospheric pressure (zero gauge) is specified and velocities are extrapolated from the interior. At cylinder wall, the y directional velocity is specified which is based on the cylinder displacement computed from (14). For heat transfer case, The cylinder wall is maintained at high temperature Tw =1 as compared to the fluid supplied.

Case i)Vortex-induced vibration (VIV) with 1 DOF at Reynolds number (Re) 150
In order to provide the initial condition to the isothermal case of VIV problem, firstly unsteady flow over cylinder at Reynolds number of 150 is obtained. Mandal et al [17] had presented the comparison of lift and drag coefficient at Re= 150 with literature. With unsteady flow at Re = 150 as initial conditions, the VIV computation is simulated and cylinder is allowed to move in Y direction as per forces (lift forces) applied by the flow. The lift and drag forces whereas calculated as average over a cylinder wall. The cylinder displacement where evaluated and compared with the literature. Figure 5 shows the comparison of the maximum traverse displacement attained by the cylinder for a various reduced velocities �U r = U (D. f n ) � � at Re =150. The plot also compares the results obtained with the literature [24 -26]. From figure 5, it can be seen that the simulation results produced using present formulation agree well with the literature data [24 -26]. The maximum transverse amplitude occurs at a reduced velocity U r =4. with Ymax= 0.51. Hence we can infer that the lock-in region for Re=150, m=8/π lies at =4.

Case ii) Vortex-induced vibration (VIV) 1 DOF at Reynolds number (Re) 200 with heat transfer
In order to simulate VIV with heat transfer case, initially unsteady flow cylinder at Re = 200 with heat transfer is simulated. The average Nusselt number, lift and drag coefficient is obtained and compare with that of literature [11]. Now as initial condition (obtained as unsteady results above) provided the simulations of vortex induced vibration with heat transfer (with cylinder maintained at constant temperature, T W ) of circular cylinder having diameter D =1 is performed here. The simulations are performed for reduced mass m*=1, Pr = 0.7, damping co-efficient ξ = 0.01 at different natural frequencies, Stn = 0.0633, 0.17, 0.21 and 0.95, which is similar to that of [11]. Fig 6 shows the comparison of the average Nusselt number obtained at various natural frequencies with the literature data [11]. The matching of results is found to be very good.

Conclusion
An accurate Harten Lax and van Leer with contact for artificial compressibility (HLLC-AC) Riemann solver with arbitrarily Lagrangian-Eulerian (ALE) formulation has been developed and used for computing flow past a circular cylinder vibrating due to induced vortex. The solution dependent weighted least square (SDWLS) based gradient method is applied to produce higher order accurate results over unstructured meshes. The numerical results produced by current solver matches well with numerical results reported in the literature.