High Order Accurate Numerical Simulation of Vortex-Induced Vibrations of a Cooled Circular Cylinder Case using Solution Dependent Weighted Least Square Gradient Calculations

In this paper, the fluid-structure interaction problem: vortex-induced vibration of a cooled circular cylinder involving thermal buoyancy is numerically investigated. The elastically mounted cylinder having a temperature lower than the flowing fluid is modelled using mass-spring-damper hence allowed to vibrate in the transverse direction to the flow direction. The gravity is acting opposite to the flow direction. In-house fluid-structure interaction solver is developed based on Harten Lax and van Leer with contact for artificial compressibility Riemann solver. The arbitrarily Lagrangian-Eulerian formulation is employed here, and the mesh is dynamically moved using radial basis function-based interpolation. The solutiondependent weighted least squares based gradient calculations are developed to achieve higher-order accuracy over unstructured meshes. The laminar incompressible flow at Reynolds number, Re = 200, and Prandtl number, Pr = 0.71, is simulated for the mass ratio of 1 and reduced damping coefficient of 0.001. The flow is investigated for Richardson number (-1, 0) and over a wide range of natural frequencies of the cylinder. The heat transfer characteristics from a cylinder are captured and compared with the existing literature results. From the study, it can be observed that in the presence of the thermal boundary layer, the oscillation of the cylinder increases to its maximum amplitude, particularly for values of natural frequencies (0.063 – 0.3).


Introduction
In engineering applications, flow over bluff bodies like circular, square, etc., are found almost in all applications. To mention few examples: flow over heat exchanger tubes, tubes in power plants, subsea cables for electronic communication and power transmission, submarine oil pipelines, offshore risers, etc. These structures experience the vibration mainly caused due to the surrounding fluid flow hence known as flowinduced vibration (FIV). When periodic lift and drag forces are generated over the cylinder due to separation and vortex shedding leading to structural vibrations are known explicitly as vortex-induced vibrations (VIV). This fluid-structure interaction (FSI) may result in serious engineering problems such as fatigue or fretting failure of the offshore riser and subsea pipelines, vibrations of a tube of the heat exchanger, vibration in structures such as a chimney, high rise building, etc. It has also been found that FIVs influence the convection heat transfer process, where the boundary layer separation, vortex shedding, and wake pattern are affected due to the FIV's. Hence, the prediction of fluid forces and heat transfer characteristics of a cylindrical structure and providing a reliable solution either to suppress or to utilize these vibrations have great importance. * Corresponding author: crsonawane@gmail.com VIV of a circular cylinder has been extensively studied over the last few decades, but the physics behind this complex FSI problem is still elusive and remains the core of the FIV modelling. A detailed review of FIV can be found in [1 -3]. Recently, many researchers have investigated the FIV of circular/cylindrical structures involving heat transfer. Pottebaum and Gharib [4] experimentally investigated the heat transfer characteristics of a circular cylinder due to the transverse oscillations in cross-flow. They observe that the heat transfer enhancement is dependent on cylinder transverse velocity, the cylinder wake mode, and the synchronization with harmonics of the natural shedding frequency. Fu and Tong [5] to investigated the flow structures and heat transfer characteristics of a heated transversely oscillating cylinder in a cross-flow and found that the heat transfer is increased remarkably, particularly in the lock-in regime. Zhang et al. [6] numerically investigated the effect of Reynolds number at various temperature boundary conditions using the immersed-boundary method over a forced oscillating cylinder involving heat transfer. Su et al. [7] used CFX-Ansys software to investigate the FIV and heat transfer of an elastically supported circular cylinder. They observed that cylinder vibration amplitude in the transverse direction is much larger than that of the streamwise direction. They also observed that heat transfer is not significant when the flow velocity is comparable to the cylinder oscillation velocity. Zhou et al. [8] used the discrete vortex method to study one and two degrees of freedom FIV of an elastically mounted circular cylinder. They observed a qualitative agreement between the simulated results obtained from the two degrees of freedom simulation with one degree of freedom. They obtained the maximum amplitude of cylinder vibration when the natural frequency of the fluid-structure system matches with cylinder oscillation frequency. Baratchi et al. [9] investigated the FIV of a heated elastically supported cylinder using the moving overset grids method. They studied the force coefficients, the amplitude of oscillations, vortex shedding pattern, and Nusselt number at Re=200 and Pr=0.7. Wang et al. [10] investigated the forced convection heat transfer characteristics of a heated and elastically supported circular cylinder. They observed that the Nusselt number increases significantly in the lock-in region with increased cylinder vibrating amplitude. Wan and Patnaik [11] investigated the mixed convection heat transfer from elastically supported and heated circular cylinder under transverse VIV. They observed that the VIV of a cylinder could be entirely suppressed for critical values of Richardson number. Izadpanah et al. [12] also investigated the forced convection heat transfer characteristics of a heated and elastically supported circular cylinder under transverse VIV. They observed that heat transfer is significantly varied due to reduced velocity and damping ratio. Chatterjee [13] has investigated the effect of thermal buoyancy on a stationary cylinder at a low Reynolds number. Sharma and Eswaran [14] also investigated the effect of buoyancy on the flow structure and heat transfer characteristics of an isolated square cylinder. Singh and Chandar [15] investigated the thermal effects on the vortex shedding produced by an unsteady flow over a circular cylinder at a Reynolds number of 150. The flow direction and thermal buoyancy force were acting perpendicular to each other. To suppress the vortex shedding, they divided the cylinder into four quadrants and heated and cooled the different quadrants simultaneously. They also concluded that only heating the cylinder surface cannot suppress the vortex shedding completely. Recently, Garg et al. [16] investigated the VIV of an elastically mounted circular cylinder in the presence of thermal buoyancy. The thermal buoyancy was induced by two parallel plates (bottom plate heated, and top plate was cooled). They studied the effect of Richardson number, Prandtl number, mass ratio, and reduced velocity on the cylinder vibrations, oscillation frequency, lift forces, and wake structure. They observe that the thermal buoyancy affect the cylinder vibrations are increases significantly for a wide range of reduced velocities. Garg et al. [17] also investigated the effect of thermal buoyancy over a heated and elastically supported circular cylinder for cross-flow conditions. At higher Richardson numbers, they observed the galloping responses due to the transversely acting thermal buoyancy. Garg et al. [18] further studied the effect of thermal buoyancy on a VIV of a cooled cylinder. They observed that when the flow is aligned in the opposite direction of the thermal buoyancy, a wider lock-in region along with a larger cylinder vibrating amplitude appears. From the above literature review, it can be seen that very few investigations are available on the VIV of a circular cylinder involving thermal buoyancy. Hence, in this paper, we developed an in-house fluid-structure interaction solver for simulating the VIV of a cooled cylinder (cylinder surface is kept at a lower temperature than fluid flow temperature) involving opposing thermal buoyancy. An accurate Harten Lax and van Leer with contact for artificial compressibility (HLLC-AC) Riemann solver [21 -24] are developed for solving incompressible flows in artificial compressibility formulation. The Riemann solver is modified to incorporate arbitrarily Lagrangian-Eulerian (ALE) [25, 19 -20] formulation to take care of moving mesh scenarios where radial basis function [26] based interpolation is used for dynamically moving the mesh. Higher-order accuracy is achieved using solutiondependent weighted least squares (SDWLS) [27] based gradient calculations over unstructured meshes

Arbitrarily Lagrangian-Eulerian (ALE) Formulation
The dimensionless form of the governing equations for the unsteady incompressible flows involving the thermal buoyancy effect can be written as = 0 (1) where i = 1 and 2 correspond to the X and Y direction; respectively, all symbols are defined in the nomenclature. Equation (1) - (3) is modified using artificial compressibility [28] methods with a dual-time stepping approach. To take care of moving boundaries, an ALE formulation [25] is employed. Hence, the modified integral form of the two-dimensional incompressible flow can be written as The convective fluxes in equation (6) can be evaluated using HLLC-AC [21][22][23][24][25]30] Riemann solver developed for stationary boundary problems. The ale flux vector is nothing but the volumetric increment along the face and can be evaluated by considering the Geometric Conservations Law (GCL) [26]. The radial basis function-based interpolation [29] is used here for dynamic mesh movement. The viscous fluxes are evaluated by the Green-Gauss approach based central differencing method. Figure 1 shows the typical stencil used for unstructured meshes. Utilizing the Taylor series of expansion, the left and right state values of a flow variable, W at the interface of cells i and j can be given as

Solution reconstruction -evaluation of interface values of primitive variables (WL and WR):
The quadratic order of accuracy can be obtained by keeping three terms on the right-hand side (RHS) of the equation (7) at the cell centre, need to be evaluated for two-dimensional simulation. Figure 1 (b) shows the typical stencil for vertex-based neighbours used for the current weighted least squares formulation. For any cell i which has n number of neighbouring cells (j =1, ..., n), using Taylor's series expansion applied over these neighbouring cells, the solution variable at the cell centre i can be given as

Solution dependent weighted least square (SDWLS)
The required spatial derivatives can be evaluated by solving the least square formulation with weights as where Wt = diag[Wt,1, ... , Wt, n] is the diagonal weight matrix, which is dependent on the solution itself [27,31] and is solved using the QR decomposition [32] algorithm.

Structural modelling of flow-induced vibration
As shown in figure 2, the mass-spring-damper system was used to model the elastically supported circular cylinder. The dynamic equation of this system can be written as Here , ̇, and ̈ are the non-dimensional displacement, velocity, and acceleration of the cylinder, respectively.   Figure 2 shows the domain as well as boundary conditions used for VIV of a circular cylinder, where the outer boundary (inout) is located at 15 times diameter from the cylinder wall. The various boundary conditions are, at the inlet, vertical velocity (keeping x component of velocity as zero) is supplied corresponding to Reynolds number of 200, and pressure is interpolated from the interior domain. At the outlet, the atmospheric pressure is specified, and velocities are extrapolated from the interior domain. At the cylinder wall, the noslip boundary condition is maintained for the stationary cylinder, whereas for the VIV of a cylinder, the vibrating cylinder velocity is evaluated based on the structural dynamic equation and provided at the cylinder wall. Also, the cylinder wall is maintained at a constant temperature, which is at a temperature lower than that fluid flow temperature; hence the Richardson number is negative. The gravity is acting opposite to the flow direction. Figure 3 shows the hybrid mesh used for the current study having 30410 mesh elements and 26216 nodes. The mesh is selected based on mesh convergence study (details are not presented here). The first quad mesh is located at a distance of 0.005 from the cylinder wall distance.   Figure 4 shows the lift and drag coefficient, whereas Table 1 shows the variation of Nusselt number for stationary cylinder compared with the numerical results of Baratchi et al. [9].

Result and Discussion
To investigate the VIV of the circular cylinder under the influence of thermal buoyancy, now the cylinder is allowed to vibrate in a transverse direction. The gravity is acting opposite to the flow direction. The cylinder is maintained at a lower temperature than the fluid temperature; hence the Richardson no, is negative. Firstly results for Ri = 0, i.e. in the absence of thermal buoyancy, are simulated. Figure 5 shows the comparison of present results, cylinder amplitude achieved for a range of Strouhal number achieved due to the VIV of a circular cylinder without thermal buoyancy (Ri =0), with that of Baratchi et al. [9]. Figure 5 also shows the VIV characteristics with thermal buoyancy effect for Ri = -1 condition. It can be observed that due to thermal buoyancy, the cylinder vibration increases significantly, particularly for the smaller strouhal number values. It can also be seen that the cylinder vibrates with a larger amplitude for a broader range of strouhal numbers

Conclusions
The laminar incompressible flow at Reynolds number, Re = 200, and Prandtl number, Pr = 0.71, is simulated for Richardson number (-1, 0) and over a wide range of natural frequencies of the cylinder. It can be seen that that in the presence of the thermal boundary layer, the oscillation of the cylinder increases to its maximum amplitude, particularly for values of natural frequencies (0.063 -0.3).