Simulation of vortex-induced vibrations of a cylinder using ANSYS CFX rigid body solver

This article simulates the vortex-induced oscillations of a rigid circular cylinder with elastic support using the new ANSYS CFX rigid body solver. This solver requires no solid mesh to setup FSI (Fluid Structure Interaction) simulation. The two-way case was setup in CFX only. Specific mass of the cylinder and flow conditions were similar to previous experimental data with mass damping parameter equal to 0.04, specific mass of 1 and Reynolds number of 3800. Two dimensional simulations were setup. Both one-degree-of-freedom and two-degree-of-freedom cases were run and results were obtained for both cases with reasonable accuracy as compared with experimental results. Eight-figure XY trajectory and lock-in behavior were clearly captured. The obtained results were satisfactory.


Introduction
When fluid flows externally over the cylinder in a direction normal to the cylinder, the downstream wake is highly disturbed. Above a critical Reynolds number, the phenomenon of vortex-induced vibrations (VIVs) is initiated. The contrarotating vortices are developed in the immediate downstream wake first on one side of the cylinder and then on the other. Vortices grow from the tube at regular intervals so that a state of regularly disordered flow prevails downstream of the tube. The train of vortices is called the Karman vortex streets. The alternate formation of vortices is accompanied by changes in the pressure distribution of fluid surrounding the tube. This causes a harmonic force to appear perpendicular to the direction of flow acting on the tube and providing the force to induce vibration in the tube.
VIV analysis plays an important role in offshore ocean engineering as well as in industrial and nuclear steam generators and heat exchangers. A large variety of commercial CFD (Computational Fluid Dynamics) packages are emerging and CFD simulation has become more reliable than before. Companies are adopting these packages for efficient production management. CFX, FLUENT, NUMECA, and NASTRAN are a few of these packages.
VIV simulation is normally done by coupling the CFD program to a structural solver like ANSYS Mechanical, and the technique is also known as FSI (Fluid Structure Interaction). Izhar et al. (2014) used CFX and an external mathem-atical solver to obtain one way coupled VIV. FSI is simple to setup by employing these packages. Different people have contributed in this area, Kuehlert et al. (2008) used FLUENT to simulate VIV in steam generator tubes. Izhar et al. (2013) employed CFX-MFX coupling to obtain flow induced displacements and stresses in a full 3D model of a shell and tube heat exchanger.  simulated two-dimensional wake interactions of two cylinders using FLUENT with a focus on harnessing the vortical energy. Lobo et al. (2011) presented VIV simulation of two cylinders and presented a hydrokinetic energy harvesting system from VIVs of submerged bodies by using a CEL (CFX Expression Language) code in CFX for the motion of rigid body. Chen and Li (2009) coupled CFX 10.0 with AN-SYS to find VIVs of an inclined stay cable.  simulated forced oscillations of flexible pipe using MFX multifield coupling scheme provided by ANSYS. MFX can transfer various combinations of variables at either side of fluid-solid interface. Both codes from fluid and structural sides are iteratively coupled. During the outer loop, simulation proceeds; during the stagger loop, the coupled codes have iterated within the given step until a complete implicit solution is obtained. The implicit coupling procedure is critical for obtaining an accurate solution. ANSYS Mechanical and ANSYS CFX can run either simultaneously or sequentially during the stagger loop. Variables are exchanged at each stagger iteration. The above men-tioned researchers have generally coupled the CFD solver to an external structural solver through the MFX solver or by employing custom developed code in either CFX or FLU-ENT to setup a two-way FSI simulation. Coupling with external solver is required to obtain displacement and/or stresses in the solid which is either rigid or constrained. If the solid is considered as rigid body, then ANSYS CFX rigid body solver provides an efficient alternative.
This paper introduces the new CFX rigid body solver which does not require coupling with an external solid mesh to compute the displacements in the solid. The prime advantage of the solver is that actually no solid mesh is needed but only the wall of the solid is defined in CFX-Pre which enables the external fluid to interact with solid wall mesh through the parameters like the mass of the body. The rigid body solver can calculate for three translational and three rotational degrees of freedom. Hence a strongly coupled two way interaction is set up with much less computational load as compared with the coupling of two different solvers if rigid body solution is required. The displacement, velocity and acceleration of the rigid body are therefore calculated.

CFX rigid body solver
The rigid body solver in ANSYS CFX allows two-way FSI without coupling the fluid solver to external structural solver like ANSYS Mechanical. The rigid cylinder supported by elastic spring and damper representing a mass spring damper system was modeled in an efficient way. The solution algorithm uses a linear momentum solver for the translational motion and a separate angular momentum solver for rotational motion. The translational equation of motion for rigid body is written in Eq. (1). This equation states that for moving rigid body the rate of change of the linear momentum P of the rigid body is equal to the applied force F acting on the body. dP dt = F: In our case translational motions are considered ignoring the rotations of the cylinder. The rigid body algorithm allows for three degrees of freedom for translation motion i.e. X, Y and Z. Eq. (1) becomes where m is the mass of the body, and F is the sum of all forces including hydrodynamic forces, weight of rigid body, and spring/external force. Therefore F is represented as: (3) where F Hydro is the hydrodynamic force, mg is the weight of the cylinder, F ext is the external damping force, and the third term represents the spring force, where K spring is spring stiffness, X n+1 is the position of the cylinder at the current time step which is to be evaluated, and X so is the neutral position of cylinder.
Using the Newmark Integration scheme given in Hughes (1987) and ANSYS (2012), we can write where h is the spatial variable, subscript n represents the known value at the current step and subscript n+1 represents the unknown value at the next step. This scheme depends on two real parameters α and µ. These parameters are directly related to the accuracy and stability of the scheme.
Typically, α and µ are taken as and , respectively. This choice of parameters corresponds to trape-zoidal rule which results in a second-order accurate scheme which is unconditionally stable. Eq.
(2) is represented as: m Ä X n+1 = F = F Hydro +mg ¡ k spring (X n+1 ¡X so )+F ext : (6) Ä X n+1 By substituting into Eq. (4), replacing h by X, and rearranging, X n+1 is obtained as follows: where X n is the position at previous time step. Ä X n +1 From Eq. (4) replacing h by X and arranging for , we have Similarly, from Eq. (5) replacing h by X, we can have the value of :

Simulation setup
The CFX-Pre case was setup by ANSYS 14.5. The unstructured mesh was formed in CFX MESH which is available for ANSYS 12. The mesh file (gtm format) was imported into CFX 14.5. The boundaries defined for the case are shown in Fig. 1. Inflation layers were used near the cylinder to accurately simulate the boundary layer. The mesh was refined around the cylinder and in the wake of the cylinder to accurately represent the vortex street. Point controls with the radius of influence option were used to refine the mesh. The extruded 2D (two-dimensional) mesh option was used as the simulation which was quasi two-dimensional. The extruded distance was only one element thick in Z-direction, with two mesh layers in Z-direction. As the simulation was two-way fluid-structure interaction, a short region around the cylinder was defined as subdomain as shown in Fig. 2. The cylinder wall was setup as the rigid body. This is the main advantage of this solver that solid structure mesh is not required. The mass of the structure is input and the solver can evaluate the state variables i.e. position, velocity and acceleration of the structure as the simulation proceeds. The subdomain was defined to follow the motion from the rigid body solution to preserve the boundary layer mesh around the cylinder. Moreover, the cylinder wall inside the fluid subdomain was also set to follow the rigid body solution whereas the fluid domain outside the subdomain formed a fluid-fluid interface with one side 'Side 2' in the subdomain and other side, 'Side 1' in the outer fluid domain. The mesh motion on both sides was set to follow the rigid body solution but as the subdomain is also set to follow the rigid body solution, the whole subdomain moves rigidly with the cylinder to preserve the boundary layer mesh. The outer fluid domain absorbs the motion of the cylinder as only 'Side 1' is set to follow the rigid body solution. The solver facilitates the use of the spring, and the spring stiffness was input in rigid body parameters. The damper was included by using the external force option available for the rigid body setup. The spring damper representation for the rigid cylinder for two degrees of freedom is shown in Fig. 1. The height of the computational domain was kept sufficiently large and equal to 110 m in order to minimize the effects of blockage ratio on the results. The blockage ratio BR is defined by Eq. (10) where D is the diameter of the cylinder and H is the height of domain in the cross flow direction. For our case, .
The blockage ratio affects the downstream wake as well as the displacement behavior of the cylinder. Prasanth et al. (2006) analyzed the response of the cylinder for blockage ratios of 5% and 1% and concluded that the blockage ratio of 1% was required to minimize the effect of hysteresis. Hysteresis causes a jump in the displacement during the transition from the initial to upper branch, hence the vertical height of the domain was set to 110 m. The length is divided into upstream (L u ) and downstream (L d ) lengths, respectively. L u was set to 25 m and L d was set to 100 m. The step size was based on dividing the period by 40 steps. A general coupling control option was used as it is the most robust scheme available. The maximum number of 19 stagger iterations per time step was used to achieve stability and accuracy. The chosen turbulence model was SAS-SST. This is a class of URANS (Unsteady Reynolds Averaged Navier Stokes) method and it can provide LES (Large Eddy Simulation) like behavior in the detached flow regions. SAS is based on including the Von Karman length scale into the turbulence scale equation. Von Karman length scale allows SAS technique to dynamically adapt to the resolved structures in a URANS simulation which provides an LES like behavior in unsteady regions of flow and simultaneously gives standard RANS capabilities in steady flow field. SAS-SST is particularly suited to flows with large separation zones and vortex formations as described in ANSYS (2012). High resolution advection and second-order backword Euler schemes were used.

Input parameters
The specific mass of the cylinder (m s ), i.e. the mass per unit length of the cylinder is represented as follows: where m is the mass of the cylinder, ρ is the density of the fluid, and L is the length of the cylinder. The mass m of the cylinder was calculated from Eq. (11) for the specific mass of cylinder, m a is the added mass i.e. the mass of fluid displaced by the cylinder. This can be defined as the effectively additional inertia of an accelerating structure in fluid caused by the acceleration of fluid in addition to the structure as mentioned by Kuehlert et al. (2008). The added mass is calculated from Eq. (12).
(12) where r is the radius of the cylinder.
The natural frequency (f n ) and the spring constant (k sp ) of the system are related, as given by Eq. (13): By rearranging Eq. (13) and substituting for f n , the spring constant is calculated by Eq. (14) This is mentioned here that the natural frequency of the system was assumed to be 14 Hz and the corresponding parameters like the damping coefficient adopted the assumed natural frequency of the cylinder, the diameter of the cylinder, and the added mass.  IZHAR Abubakar et al. China Ocean Eng., 2017, Vol. 31, No. 1, P. 79-90 81 The mass damping parameter m s ζ was taken as 0.04. With the specific mass, m s =1, the damping constant c is calculated by Eq. (15) as follows: where ζ is damping ratio. Therefore we can obtain the values of the damping constant and the spring stiffness of the spring damper system for the assumed natural frequency of 14 Hz.
3.2 1DOF (one degree of freedom) simulation The cylinder was setup to move in the lift (Y) direction only for a one-degree-of-freedom solution. The translational degree of freedom was set to Y-only in the rigid body dynamics tab in CFX-Pre. The rotational degrees of freedom were ignored. The linear spring constant as well as the external damping force was considered in the Y-direction only. The equilibrium position was at the center of the rigid body coordinate frame. Twenty cases were run gradually to increase the velocity at the inlet boundary thus give 20 reduced velocities corresponding to the assumed natural frequency and diameter of the cylinder. Each case took approximately 2200 steps with a time step of 0.001 s. Within this number of iterations, approximately 14 periods of oscillations were obtained to reasonably estimate the frequency of oscillation. The rigid body displacement was monitored in solver. Fig. 3 shows the vorticity contours obtained through simulation at different reduced velocities (V r ). This figure shows both the 2S mode and 2P mode. 2S mode has two vortices generated in each cycle and 2P mode has two pairs of vortices generated in each cycle with each pair consisting of two vortices. 2S and 2P modes comprise the initial and upper branch of response, respectively as were observed by Williamson (1997, 1999).
The displacement signals for several reduced velocities obtained are shown in Figs. 4-7. These signals show that displacement is well converged and steady state oscillation is achieved. Fig. 8 shows the amplitude response as the function of the reduced velocity. The average of upper 10% amplitude values were calculated from the data at each reduced velocity to obtain this amplitude response. The frequency of oscillation for each case was obtained using Discrete Fourier Transform (DFT) analysis of the displacement signal in CFD-Post. The frequency spectrums at several reduced velocities are shown in Figs. 9-11. Fig. 12 shows the ratio of oscillation frequency to the      natural frequency of the cylinder. It is obvious from the graph that the oscillation frequency remains close to the natural frequency of the cylinder for the reduced velocities ranging from about 4 to 7.9 and hence it somewhat represents a straight line in this region. This portion is the lock-in region in which the vortex shedding frequency coincides and locks with the natural frequency of the system and the tube oscillates with larger amplitude. This region of reduced velocities has contained the upper branch and lower branch with 2P vortex shedding mode until the desynchronization region starts at larger reduced velocities from about 9 to 12 in which the frequency of oscillation gets higher than the tube natural frequency, the tube comes out from the lock in the region and the amplitude of oscillation decreases. The initial branch corresponds to the reduced velocities ranging from 2 to 4 in which the frequency of oscillation is much lower than the natural frequency of the tube and is characterized by 2S vortex shedding mode. Fig. 13 shows the evaluated values of the lift coefficient obtained from the simulation. The force_y function of CFX was used to monitor the lift force (F y ) at the rigid body wall during the run and the coefficient can be calculated from Eq. (16) where A, L and D are the projected area, length and diameter of the cylinder, respectively. The values of C L for all time steps at each reduced velocity (V r ) were evaluated using the above equation and plotted. One value of C L at each reduced velocity was obtained from each respective plot, which was equal to the amplitude of oscillation of the lift coefficient plot. It is evident from Fig. 13 that the lift coeffi- Fig. 9. Frequency magnitude spectrums for 1DOF simulation at V r =3 and 4.5. Fig. 10. Frequency magnitude spectrums for 1DOF simulation at V r =5.25 and 7.9. Fig. 11. Frequency magnitude spectrums for 1DOF simulation at V r =11. Fig. 12 where f s is the vortex shedding frequency and U is the incoming free stream velocity. Khalak and Williamson (1999) mentioned that during resonance the cylinder oscillation frequency (f o ) is equal to the natural frequency (f n ) of the cylinder and is also equal to the vortex shedding frequency from the static cylinder (f ss ). f ss is known as the Strouhal frequency. These frequencies are equal at resonance i.e. f o =f s =f n =f ss .The lift coefficient (C L ) in this range has larger values due to 2S vortex shedding mode which imparts greater energy from fluid to the cylinder compared with 2P mode.

Introduction
There are few studies available in literature that deal with cylinder response with two degrees of freedom (Jauvtis and Williamson, 2003. Sanchis et al. (2008) carried out experimentation with XY degrees of freedom. They have observed that the freedom to oscillate in inline direction does not affect the lift response of the cylinder much as compared with Y-only case. They have clarified that previous Y-only studies still hold their validity when compared with 2DOF. The amplitude response and frequency ratios from Jauvtis and Williamson (2003) are shown in Figs. 14 and 15, respectively. It was found that the peak amplitude response for 2DOF case was 10% larger than that for Y-only case. Furthermore, the response branches in the Y direction were the same in two cases i.e. initial, upper and lower branches were easily observed. Moreover, the griffin plot for Y motion did not show much difference for XY and Yonly cases. The mass damping parameter (m s +C A ) ζ for the griffin plot was between 0.001 and 0.1, where C A is the added mass coefficient. The vortex modes corresponding to different response branches for Y-only motion were the same as for 2DOF motion and this similarity was attributed to the low amplitude inline response which did not have much effect on wake. Sanchis et al. (2008) studied 2DOF response of a rigid low specific mass (m s =1.04) cylinder. The mass damping parameter was high enough which resulted in a two-branch type response and super upper branch (SUB) not being observed. The peak lift response was comparable to Y-only case. It was concluded that, with the combination of low specific mass m s and high mass damping (m s +C A ) ζ, the system exhibited the characteristics of high mass damping systems like that of Feng (1968) where the upper branch was absent and the peak amplitude was obtained in the initial branch and the frequency response was much similar to the small mass damping. Sanchis et al. (2008) found that the response was increased by a combination of low specific mass and two degrees of freedom as their data were in the ranges investigated by Feng (1968), Khalak and Williamson (1999) at moderate mass ratios and the same mass damping.

CFX simulation setup for 2DOF
The two degrees freedom setup was similar to that of 1DOF except the translational motion was setup in X and Y directions. The spring stiffness was specified in both X and Y directions, and the damper was also included in the X direction in addition to the Y direction. Fig. 16 shows vorticity    Kuehlert et al. (2008) and Dahl et al. (2006) with Crescent and Figure  Eight leaning downstream. It can be seen from these figures that patterns in the initial branch like at V r equal to 2 and 3.2 as well as in desynchronization branch (V r =9.2 and 11) do not show regular behavior as compared with regular patterns in the upper and lower branch regions of V r =4 to 7.5 which is also consistent with previous observations. Fig. 18 shows the temporal profiles of lift and inline force coefficients on the cylinder at different reduced velocities. These figures show that the forces are well converged to steady oscillations state. It is apparent that inline forces have higher frequency than that of lift forces which have approximately half the frequency of inline forces. This frequency difference is analogous to previous results.
Figs. 19 and 20 show the nondimensional lift (Y/D) and inline (X/D) displacements of the cylinder from 2DOF simulations, respectively. It is observed that the Y displacement is not altered much as compared with 1DOF lift displacement. The main reason is that the low X displacement around 0.04 does not significantly affect the Y displacement. The frequency ratio from 2DOF simulation is shown in Fig.  21, and it is obvious from Fig. 21 that the frequency of oscillation increases with reduced velocity V r and as a result the frequency ratio increases, the cylinder comes out of lock-in region at V r ≈7.9 and ratio reaches a maximum value of 1.59 at V r ≈11.
A total of 20 cases were run with the varied reduced velocity. The displacement data for 2DOF simulation were imported into CFDPost and frequency spectrum was obtained by using FFT at each reduced velocity. Fig. 22 shows the obtained X and Y oscillation frequencies. The red line indicates the frequency of drag oscillation and the blue line indicates the frequency of lift oscillation. It can be seen from the figure that the red line has the double frequency value of the blue line except for some initial and final values of reduced velocities. It is because at initial and final values of V r , the IZHAR Abubakar et al. China Ocean Eng., 2017, Vol. 31, No. 1, P. 79-90 cylinder does not follow a regular figure of Eight patterns and the response is random without clear path. However at these values of V r , the secondary frequency peaks do show up at the double value of lift frequency peaks.

Convergence
It is mentioned here that CFX allows the flexibility to set the under-relaxation factor and convergence targets for force and mesh motion data transfer but this option was kept off as the default settings to give reasonable accuracy with an affordable number of computational iterations. These solver control settings provide reasonably accurate results with force data transfer residual convergence criterion equal to 10 -2 and motion data transfer residual criterion equal to 10 -4 . Forcing the force data transfer residual to lower value like 10 -4 increased the number of iterations two folds and there was a little difference in the simulated results as compared with that obtained with convergence target equal to 10 -2 . The residuals convergence behavior is given in Fig.  23. It can be seen in the figure that residuals are well converged. Hover et al. (1998) studied the free oscillation response of cylinders with the span of 60 cm experimentally at Re=3800. A hybrid system was developed to study free motions employing a pair of force transducers to measure the lift forces at both edges of the cylinder. Mass damping (m s ζ) of 0.04 was used with a specific mass of 1. Fig. 24 shows the comparison of the amplitude response of the two simulations and the experimental results of Hover et al. (1998). It is observed that non-dimensional response Y/D from Hover et al. (1998) is around 0.79 while that from current simulations is around 0.64. This discrepancy can be attributed to the two-dimensional approximation for the current simulation. The experimental setup of Hover et al. (1998) comprised three-dimensional cylinders with the span of 60 cm. Williamson and Govardhan (2004) stated that numerical simulations can predict the peak amplitude is around 0.6D   when experiments produce peak response well over 1.0D. This statement is seen to be agreed with these results. Xie et al. (2012) carried out two-dimensional numerical simulations for a flexible cylinder for Re=1000. They found a peak response of 0.57D for 2D (two-dimensional) and a peak response of 1.02D for 3D (three-dimensional) simulation. Wang et al. (2010) presented high resolution of 2D CFD simulation for VIV of risers. They obtained peak amplitude of 0.6D as compared with the 0.9D experimentally observed by Khalak and Williamson (1997). Blackburn et al. (2001) carried out a comparative study of experimental and 2D/3D simulation results and concluded that 3D simulation is actually required to produce the response envelope observed experimentally. The response curve from Blackburn et al. (2001) is shown in Fig. 25. The X-axis in Fig. 25 is the normalized and reduced velocity in order to account for the disparity of the Strouhal numbers in 2D and 3D cases. The maximum amplitude in this figure from experiments is larger than that of 2D and 3D numerical simulations. The 3D simulation results are closer to but still lower than experi-   IZHAR Abubakar et al. China Ocean Eng., 2017, Vol. 31, No. 1, P. 79-90 87 mental results. Moreover, the range of the reduced velocity for which the upper branch extends is also small in 2D case. Blackburn et al. (2001) described that the reason for lower amplitude in simulation could be due to the small domain length in axial direction. Somehow the results obtained through current simulations agree well with the results of Sanchis et al. (2008) whose experimental data are almost similar to our simulation data with the exception of Reynolds number hence keeping in view the results from the above mentioned references. It can be stated that the trend of the response amplitude from the current simulation is similar to that of experiments and is much satisfactory as far as 2D simulation is concerned. The two degrees of freedom results from current simulation are also shown in Fig. 25 and it is observed that 2DOF results are not much different in the initial and lower branches while the amplitude in the de-synchronization region is larger than that of one degree freedom simulation. Sanchis et al. (2008) verified that for systems characterized by high mass damping, the ability to move in-line and cross-flow with a very small mass ratio m s =1 does not affect dramatically the lift response and this is the case with our simulations. Fig. 26 gives comparison of frequency ratios for the three cases and the ratios are rather close to each other, showing little effect of inline oscillations on the lift oscillations. Fig. 27 shows the lift coefficients from the simulations and experiments. The coefficients agree well however some difference is observed in the upper and lower branch regions. The difference in the lift coefficients is acceptable as the coefficients were obtained experimentally for cylinders with the span of 60 cm with two force signals measured at each edge of the cylinder. Fig. 28 shows the graph of phase difference (in degrees) between the lift force and the displacement signals. It is observed that general trend of the phase difference from the simulation is similar to that of experiments; however, some scattering of the data points is observed from the simulation. It is because the frequency spectrum was obtained by applying the Discrete Fourier Transform (DFT) to the obtained numerical data and DFT has the limitation such as the leakage or spread to adjacent frequencies as well as the effects of the discrete character of the numerical method which causes an inaccurate calculation of phases. Moreover, the lack of periodicity in the force or displacement signals obtained through simulation adds inaccuracy to the phase difference calculation. Fig. 29 shows the griffin plot from different previous authors representing the peak amplitude as the function of the mass damping parameter α, and it can be seen that our data points for two simulations adjust reasonably among those from these authors. As there is much scatter in the griffin plot, Govardhan and Williamson (2006) introduced a modified griffin plot taking into consideration the effect of Reynolds number on peak amplitudes. They developed a functional relationship among A*, A M * and Re as follows:    where A* is non-dimensional amplitude and A M * is modified non-dimensional amplitude. Fig. 30 shows the modified griffin plot from Govardhan and Williamson (2006) along with the current results of our simulation for both 1DOD and 2DOF plotted on the modified griffin plot and shows reasonable agreement with the data from previous authors (Izhar et al., 2014).

Conclusion
We have introduced and demonstrated the application of ANSYS CFX 14.5 Rigid Body Solver to simulate vortex-induced oscillations of the cylinder with two way fluid structure interaction. FSI was set up in CFX standing alone without coupling to any external solver which is the basic advantage of this solver and it is more efficient as compared with CFX-MFX, ANSYS Flotran-MFS or Fluent-MP-CCI (Multi Physics Code Coupling Interface) coupling methods. The obtained results were satisfactory for 2D cases, and the lock in behaviour and the resonance phenomenon were captured with much accuracy. Eight figure and Crescent shaped motions were clearly observed for 2DOF (Two-Degree-of-Freedom) simulation. The discrepancies in the results are due to 2D approximation as the 2D case neglects the complex wake behaviour along the span of the cylinder which prevails in 3D cases. The authors are trying to apply this solver for an accurate 3D simulation to reveal the complex behaviour of VIV phenomenon.  IZHAR Abubakar et al. China Ocean Eng., 2017, Vol. 31, No. 1, P. 79-90 89