Unsteady RANS and LES Simulation of an Ideal Rankine Vortex Decay

Abstract

The 3D numerical simulation was carried out for an idealized Rankine vortex using nonlinear k-ε model (one kind of RANS model) and large eddy simulation (LES) techniques. In this 3D simulation, the vortex flow field was given to rotate with the vertical axis in a free surface rectangular domain. In order to investigate the predictability of standard (linear) and non-linear k-ε models, the decay of a trailing vortex was simulated and compared with previous DNS data. The governing equations for mean velocities and turbulent flows were discretized with the finite volume method based on a staggered grid system. It was observed that in the growth phase as well as in stabilized phase of turbulence, the decay rate of tangential velocity by RANS model was well comparable with LES simulation as well as previous DNS data. However, in the decay phase of turbulence, RANS model showed slightly faster decay of tangential velocity due to its slower decay of turbulence compared to LES or DNS. The patterns as well as magnitudes of secondary currents predicted by RANS and LES models were well comparable to each other.

1. Introduction

The basic types of plane vortices can be classified into two categories: one with slower velocity at center and maximum at sides and the other with maximum at center and minimum velocity at edges. The rotary fluid motion of the first one is called the solid body rotation, since it is similar to the fluid motion filled in a rotating hollow box. On the other hand, if a long circular rod rotates in a fluid with constant velocity around its axis, the fluid velocity is found highest and equal to the velocity of rod at the rod’s surface (due to adhesion), and with increasing distance from the rod, the velocity is diminished in inverse proportional to the distance. Such a fluid motion is called a potential vortex. Fluid motion composed of a potential vortex and solid body rotation is called Rankine vortex (Figure 1) after the fluid dynamicist Rankine. The radial distance from the center to the maximum tangential velocity is called the radius of vortex core. For a steady circular motion without a velocity component normal to the plane of rotation, the Rankine vortex is the only possible vortex whose velocity is zero at the center as well as far away from it. In addition to these basic vortices, there are other time-dependent rotary motions that have azimuthal velocity component as well as radial and axial components.

(a)

(b)

(c)

Figure 1 

Velocity distribution in (a) solid body rotation, (b) potential vortex, and (c) Rankine vortex.

The existence of such vortices is not only limited to natural phenomena but also in many engineering applications. The Rankine vortex has been used extensively in various studies: for instance, to model the natural phenomena such as hurricanes and tornados, to study the vortex wake hazard caused by flow obstructions, to predict the decay of wing-tip vortices, and so forth. Although some experimental and theoretical studies ([1–3], etc.) as well as DNS simulations [4] were successfully performed to study the structure and decay of such vortex, all of them were performed considering 2D vortex field.

As long as the flow is laminar, the only problem is to solve a system of discretized equations. The accuracy of results depends on both the grid and discretization scheme. The moment that the flow becomes turbulent, a chaotic, random motion is observable. In this case, the instantaneous flow values at a fixed location in the flow field show a fluctuating behavior. Therefore, in order to study such flows numerically, different approaches have been proposed. They have different accuracies, and sometimes their predictability depends on the type of flows. The major approaches to study turbulent flows numerically are Reynolds-averaged Navier-Stokes (RANS), large eddy simulation (LES), and direct numerical simulation (DNS). In this study, the 3D numerical simulations were carried out using non-linear 𝑘-𝜀 model (one kind of RANS model) and large eddy simulation (LES) techniques. The idea of LES comes from the fact that the small scales show an isotropic behavior regardless of the type of flow and boundary conditions. Hence, if the large scales are resolved, it is accurate enough to model the small scales which are called the subgrid scale or SGS. In order to get reasonable accuracy, the grid should be fine to some extent. Important difference between RANS and LES governing equations is the method of filtering. While in RANS the filtering is performed in time, in LES the filtering is performed in space. Another difference compared to RANS is that, as turbulence is a time-dependent three-dimensional phenomenon, the governing equations in LES also need to solve in unsteady three-dimensional format. Another fact regarding LES is that the grid is supposed to be fine enough to capture some part of scales in the inertial subrange. These differences make the LES much more expensive than RANS models.

The work presented in this paper can be classified into twofold. In the first part, unsteady numerical simulation with 3D velocity field was carried out for a turbulent trailing vortex using non-linear 𝑘-𝜀 model. To assess the predictability of non-linear 𝑘-𝜀 model, general features of vortex decay are compared with previous DNS data. The simulation was also performed using standard 𝑘-𝜀 model and comparison is shown with non-linear model. In the second part, the 3D numerical simulations were carried out for an idealized Rankine vortex using non-linear 𝑘-𝜀 model and large eddy simulation (LES). 3D staggered grid was used for both RANS and LES simulations. Comparison is shown for the general flow features such as temporal change of vortex decay, radial distribution of tangential velocity, and water surface profile.

2. Numerical Model and Simulation Details

2.1. Nonlinear 𝑘-𝜀 Model

Following 3D flow equations [5–7] were used in unsteady RANS calculation for the chosen flow field.

Continuity equation: 𝜕𝑈_𝑖𝜕𝑥_𝑖=0. (1)

Momentum equation: 𝜕𝑈_𝑖+𝜕𝑡𝜕𝑈_𝑗𝑈_𝑖𝜕𝑥_𝑗=𝑔_𝑖−1𝜌𝜕𝑃𝜕𝑥_𝑖+𝜕𝜕𝑥_𝑗−𝑢_𝑖𝑢_𝑗𝜕+𝜈²𝑈_𝑖𝜕𝑥²_𝑗. (2)

𝑘-equation: 𝜕𝑘+𝜕𝑡𝜕𝑘𝑈_𝑗𝜕𝑥_𝑗=−𝑢_𝑖𝑢_𝑗𝜕𝑈_𝑖𝜕𝑥_𝑗+𝜕𝜕𝑥_𝑗𝜈_𝑡𝜎_𝑘+𝜈𝜕𝑘𝜕𝑥_𝑗−𝜀. (3)

𝜀-equation: 𝜕𝜀+𝜕𝑡𝜕𝜀𝑈_𝑗𝜕𝑥_𝑗=−𝑐_𝜀1𝜀𝑘𝑢_𝑖𝑢_𝑗𝜕𝑈_𝑖𝜕𝑥_𝑗+𝜕𝜕𝑥_𝑗𝜈_𝑡𝜎_𝜀+𝜈𝜕𝜀𝜕𝑥_𝑗−𝑐_𝜀2𝜀²𝑘. (4)

Turbulence model (nonlinear): −𝑢_𝑖𝑢_𝑗=𝜈_𝑡𝑆_𝑖𝑗−23𝑘𝛿_𝑖𝑗−𝑘𝜀𝜈_𝑡³_𝛽=1𝑐_𝛽𝑆_𝛽𝑖𝑗−13𝑆_𝛽𝛼𝛼𝛿_𝑖𝑗. (5)

Here, 𝜈_𝑡=𝑐_𝜇𝑘²𝜀,𝑆_𝑖𝑗=𝜕𝑈_𝑖𝜕𝑥_𝑗+𝜕𝑈_𝑗𝜕𝑥_𝑖,Ω_𝑖𝑗=𝜕𝑈_𝑖𝜕𝑥_𝑗−𝜕𝑈_𝑗𝜕𝑥_𝑖,𝑘𝑆=𝜀12𝑆_𝑖𝑗𝑆_𝑖𝑗𝑘,Ω=𝜀12Ω_𝑖𝑗Ω_𝑖𝑗. (6)

In (5),𝑆_1𝑖𝑗=𝜕𝑈_𝑖𝜕𝑥_𝛾𝜕𝑈_𝑗𝜕𝑥_𝛾,𝑆_2𝑖𝑗=12𝜕𝑈_𝛾𝜕𝑥_𝑖𝜕𝑈_𝑗𝜕𝑥_𝛾+𝜕𝑈_𝛾𝜕𝑥_𝑗𝜕𝑈_𝑖𝜕𝑥_𝛾,𝑆_3𝑖𝑗=𝜕𝑈_𝛾𝜕𝑥_𝑖𝜕𝑈_𝛾𝜕𝑥_𝑗,𝑐_𝜇=𝑐_𝜇01+𝑐_𝑛𝑆𝑆²+𝑐_𝑛ΩΩ²1+𝑐_𝑑𝑆𝑆²+𝑐_𝑑ΩΩ²+𝑐_𝑑𝑆Ω𝑆Ω+𝑐_𝑑𝑆1𝑆⁴+𝑐_𝑑Ω1Ω⁴+𝑐_𝑑𝑆Ω1𝑆²Ω²,𝑐_𝛽=𝑐_𝛽011+𝑚_𝑑𝑆𝑆²+𝑚_𝑑ΩΩ². (7) The values of the model constants are 𝑐_𝜇0 = 0.09, 𝑐_𝑛𝑆 = 0.005, 𝑐_𝑛Ω = 0.0068, 𝑐_𝑑𝑆 = 0.008, 𝑐_𝑑Ω = 0.004, 𝑐_𝑑𝑆Ω = −0.003, 𝑐_𝑑𝑆1 = 0.00005, 𝑐_𝑑Ω1 = 0.00005, 𝑐_𝑑𝑆Ω1 = 0.00025, 𝑐₁₀ = 0.40, 𝑐₂₀ = 0.0, 𝑐₃₀ = −0.13, 𝑚_𝑑𝑆 = 0.01, 𝑚_𝑑Ω = 0.003, 𝜎_𝑘 = 1.0, 𝜎_𝜀 = 1.3, 𝑐_𝜀1 = 1.44, and 𝑐_𝜀2　= 1.92. In the standard 𝑘-𝜀 model, 𝑐_𝜇 = 0.09 and 𝑐_𝛽 = 0.0.

2.2. Large Eddy Simulation

The filtered equations were developed from the Navier-Stokes equations. Consider that 𝑈_𝑖 is a velocity in 𝑖 direction that contains resolved scale part 𝑈_𝑖 and subgrid-scale part 𝑢_𝑖; then 𝑈_𝑖=𝑈_𝑖+𝑢_𝑖; similarly for pressure𝑃=𝑝+𝑝^′. Substituting these two decompositions into Navier-Stokes equations, the filtered equations of LES can be obtained as follows:

Continuity equation. 𝜕𝑈_𝑖𝜕𝑥_𝑖=0. (8)

Momentum equation: 𝜕𝑈_𝑖+𝜕𝜕𝑡𝑈_𝑗𝑈_𝑖𝜕𝑥_𝑗=𝑔_𝑖−1𝜌𝜕𝑃𝜕𝑥_𝑖𝜕+𝜈²𝑈_𝑖𝜕𝑥²_𝑗−𝜕𝜏_𝑖𝑗𝜕𝑥_𝑗. (9) Here, 𝜏_𝑖𝑗 is the subgrid scale Reynolds stress defined by 𝜏_𝑖𝑗=𝑢_𝑖𝑢_𝑗−𝑢_𝑖𝑢_𝑗. (10) The most commonly used subgrid-scale turbulence model is Smagorisky model proposed by Smagorinsky [8]. In this model, the residual stress takes the Boussinesq eddy viscosity form: 𝜏_𝑖𝑗−13𝜏_𝑘𝑘𝛿_𝑖𝑗=−2𝜈_𝑡𝑆_𝑖𝑗. (11) Here, the resolved strain rate is 𝑆_𝑖𝑗=(1/2)(𝜕𝑈_𝑖/𝜕𝑥_𝑗+𝜕𝑈_𝑗/𝜕𝑥_𝑖).

The eddy viscosity 𝜈_𝑡 is determined by the Smagorisky model: 𝜈_𝑡=𝐶_𝑠Δ²2𝑆_𝑖𝑗𝑆_𝑖𝑗. (12) Δ is the filter width related to grid size calculated as the geometric average of the grid spacing in three directions, Δ=(Δ𝑥₁Δ𝑥₂Δ𝑥₃)^1/3. 𝐶_𝑠 is the Smagorisky model constant, whose value varies from 0.065 to 0.20 [9]; 0.15 was used in the present simulation.

2.3. Computational Schemes and Flow Domain

The basic equations used in LES were three-dimensional, time-dependent, filtered Navier stokes equations. In nonlinear 𝑘-𝜀 model, the turbulent stress terms were solved using two additional equations for turbulent kinetic energy (𝑘) and its dissipation rate (𝜀). The governing equations for mean velocities and turbulent flows were discretized with the finite volume method based on a staggered grid system. For the momentum equation, convective and diffusive fluxes were approximated with quick and central difference schemes, respectively. The hybrid central upwind scheme was used for the 𝑘 and 𝜀 equations. Time advancement was achieved by Adam-Bashforth scheme of second-order accuracy, in each equation. The basic equations were discretized as fully explicit forms and solved successively with the time increment in step by step. The pressure field was solved using iterative procedure at each time step. Three cases have been studied. The initial conditions and simulation details are explained as follows, and the simulation parameters are shown in Table 1.

Table 1 

Simulation parameters for different cases.

2.3.1. Case  1: Simulation of Trailing Vortex by Unsteady RANS

For the 𝑞-vortex, the direct numerical simulation results were reported by Qin [4]. The simulation for case 1 was performed with the same initial condition of Qin. Since the axial velocity of this vortex is not a function of axial distance but of radial distance, 2D numerical grid was used to simulate 3D flow field. The details of the initial flow field are given as follows.

Tangential velocity:𝑉_𝜃𝑉₀=𝑞₀𝑟1−𝑒^−𝑟2. (13)

Radial velocity: 𝑉_𝑟𝑉₀=0. (14)

Axial velocity: 𝑉_𝑧𝑉₀=1𝑞₀1−𝑒^−𝑟2. (15) Here, 𝑉₀ is the scaling velocity, related to the initial mean tangential flow, defined as 𝑉₀=𝑀_𝑝𝛾𝑞₀1−𝑒^−𝑟2. (16) Here, 𝑀_𝑝 = 0.009 and 𝛾 = 1.1209 were used in Qin’s [4] DNS simulation. As reported by Lessen et al. [10], the stability of 𝑞-vortex can be related to the value of swirl number 𝑞. In DNS simulation, the initial swirl number 𝑞₀ was chosen to be unity. The time is nondimensionalized by 𝑇=2𝜋𝑟₀/𝑉_𝜃𝑚 at 𝑡=0. Here, 𝑟₀ is the radial distance where the tangential velocity contains the peak value in the initial flow field (at 𝑡=0).

Variable grid spacing was used with dense grid at centre and coarser towards the boundary. The variable mesh size with radial distance for case 1 is shown in Figure 2. The computational domain consists of 100 grids in each (𝑥 and 𝑦) direction. The size of domain was taken sufficiently large (12𝑟₀×12𝑟₀) to overcome the interference of boundary in the vortex decay process. Cartesian grid was used for computation, and the results are presented in cylindrical coordinate using the geometric conversion. The value of Δ𝑡 was taken as 0.01 sec. The truncation error was allowed as 0.000001 and the maximum number of iterations per time step was set to be 200.

Figure 2 

Variable grid size with radial distance for case 1.

2.3.2. Case  2 and Case  3: Simulation of an Ideal Rankine Vortex by Unsteady RANS and LES

These two cases were performed for the same vortex flow field using different simulation methods: case 2 by RANS and case 3 by LES. The initial condition of the vortex is given as follows.

Tangential velocity: 𝑉_𝜃𝑉₀1=−𝑟1−𝑒^−𝑟2. (17)

Radial velocity: 𝑉_𝑟𝑉₀=0. (18)

Axial velocity: 𝑉_𝑧𝑉₀=0. (19) Here, 𝑉₀ is the scaling velocity, related to the initial mean tangential flow as defined in (16). Like case 1, in these two cases, the initial swirl number 𝑞₀ was also chosen to be unity.

Simulations were performed using 3D Cartesian grid. The size of flow domain for both the cases was 12𝑟₀×12𝑟₀ (in two lateral, 𝑥 and 𝑦 directions) ×1𝑟₀ (in axial/depth-wise, 𝑧 direction). For LES 140×140×15 grids were used in 𝑥, 𝑦, and 𝑧 directions, respectively. For RANS simulation, the grid numbers were reduced to 70×70×10. 3D staggered variable grid system was used for both RANS and LES. The other parameters for the simulation are shown in Table 1.

3. Results and Discussions

At first, the decay of a trailing vortex simulated by standard and non-linear 𝑘-𝜀 model is compared with previous DNS result. Then the decay of an ideal Rankine vortex predicted by nonlinear 𝑘-𝜀 model and large eddy simulation (LES) is compared. The comparative discussion of simulated results is presented here in after.

3.1. Comparison of RANS Prediction with Qin’s DNS Data

3.1.1. Comparison of Temporal Decay

To investigate the predictability of standard (linear) and non-linear 𝑘-𝜀 model, the decay of trailing vortex as defined in case 1 was simulated with the same initial conditions of Qin [4] using both models. In this section, the simulated results are compared with Qin’s DNS data. Comparison of temporal decay of maximum tangential velocity is shown in Figure 3. It is clearly seen that the standard 𝑘-𝜀 model predicts much faster decay than that of non-linear model. Although the non-linear model shows better comparison with DNS calculation, the decay rate is still faster than that of DNS prediction, especially at about 𝑡 > 4.0𝑇.

Figure 3 

Comparison of temporal decay of maximum tangential velocity calculated by RANS with Qin’s DNS result.

The radial distribution of circumferential velocity, circulation, and axial velocity is compared with the DNS calculation at time 𝑡 = 3.72𝑇, as shown in Figures 4, 5, and 6, respectively. At this age the vortex has already overcome the effect of initial conditions and the turbulent flow field became saturated and stabilized gaining a peak value of the turbulent kinetic energy after passing its growth periods. In these figures, the non-linear model shows much better agreement with DNS data compared to standard model. Even for non-linear model, the magnitude and position of maximum tangential velocity show slightly faster decay than DNS data. Such discrepancy is also seen in the comparison of circulation profile. The cause of deviation is explained by comparing the temporal decay of vortex with the growth/decay of turbulence in the next section.

Figure 4 

𝑘-𝜀 model prediction of tangential velocity profile compared with Qin’s DNS results (𝑡 = 3.72𝑇).

Figure 5 

𝑘-𝜀 model prediction of circulation profile compared with Qin’s DNS results (𝑡 = 3.72𝑇).

Figure 6 

𝑘-𝜀 model prediction of axial velocity profile compared with Qin’s DNS results (𝑡 = 3.72𝑇).

3.1.2. Temporal Growth/Decay of Turbulence

Production of turbulence in the vortex field is due to the tangential and axial velocity gradients. In the vortex lifetime, five different growth/decay rates of turbulent kinetic energy as reported in DNS calculation [4] were also reflected in the simulation of non-linear 𝑘-𝜀 model. Figure 7 presents the growth/decay of turbulence kinetic energy with approximated five different time periods. Initially (zone (i)), the vortex changes very slowly, as it requires adjustment of any unphysical nature of initial conditions. In time zone (ii), an exponential growth of turbulent kinetic energy is seen. The growth slows down in the next time period of zone (iii). It is found that the growth of axial velocity is significantly higher in time zone (ii) and (iii). The turbulent kinetic energy reaches its peak value and remains almost constant throughout the short period of zone (iv). It reveals that in this stabilization period the flow field becomes saturated and cannot support additional turbulence. Finally, the decay of turbulence is started as predicted by most of the previous researches (such as [2, 11], etc.). The logarithmic plot depicts that the decay rate in zone (v) is much slower than the growth rate in zone (ii). The decay of velocity field slows down as the turbulence decay period starts.

Figure 7 

Growth/decay of turbulent kinetic energy with different time periods (nonlinear 𝑘-𝜀 model).

Comparing Figure 3 with Figure 7, we can conclude that in the growth phase as well as in the stabilized phase of turbulence, the decay rate of tangential velocity by nonlinear 𝑘-𝜀 model is well comparable with DNS result. However, in the decay phase of turbulence, RANS model shows faster decay of tangential velocity due to its slower decay of turbulence compared to DNS.

Figure 8 depicts the comparison of the radial distribution of turbulent kinetic energy using three models. It is seen from the figure that the 𝑘-𝜀 model overpredicts the kinetic energy at the peak compared to DNS. This is the cause for the faster decay of tangential velocity in the 𝑘-𝜀 model (as observed in Figure 4).

Figure 8 

𝑘-𝜀 model prediction of turbulent kinetic energy compared with Qin’s DNS results (𝑡 = 3.72𝑇).

3.2. Decay of an Ideal Rankine Vortex by LES and RANS Simulations

In this 3D simulation, the vortex flow field was given to rotate with the vertical axis in a free surface rectangular domain. The temporal change of velocity at a point is shown in Figure 9 for both LES and RANS simulations. In addition to temporal decay of velocity, the typical velocity fluctuations are observed in LES results; on the other hand, RANS simulation shows a mean decay profile.

Figure 9 

Comparison of temporal fluctuations of velocity at a point by RANS and LES simulations.

3.2.1. Comparison of Temporal Decay

Figure 10 shows the radial distribution of tangential velocity predicted by LES and RANS models at 𝑡 = 3.06𝑇 and 𝑡 = 4.48𝑇. Although the velocity profiles at 𝑡 = 3.06𝑇 show good comparison, RANS model shows faster decay than LES at 𝑡 = 4.48𝑇.

Figure 10 

Radial distribution of depth calculated by RANS and LES (symbols indicates LES prediction, and lines show the corresponding RANS prediction).

The variation of water surface for different times predicted by RANS and LES is shown in Figure 11. The water surface is found to be depressed near the central area with minimum water depth at the centre of vortex. Since the strength of vortex decreases due to the decay of tangential velocity, the amount of depression in water surface is also decreased with time. Similar to Figure 10, although the depression of water surface at 𝑡 = 3.06𝑇 shows good comparison, RANS model predicted less depression in water surface than that of LES at 𝑡 = 4.48𝑇.

Figure 11 

Radial distribution of tangential velocity calculated by RANS and LES.

Figure 12 shows the calculated temporal decay of maximum tangential velocity for both LES and non-linear 𝑘-𝜀 model. It is found that although the prediction of RANS model shows good agreement with LES in initial times, some deviation is observed at about 𝑡 > 4.0𝑇. Such discrepancy in the decay phase of turbulence is also observed in the comparison of non-linear 𝑘-𝜀 model predictions with previous DNS results of trailing vortex decay (case 1). Therefore, similar to previous case, it can be concluded that in the growth phase as well as in stabilized phase of turbulence, the decay rate of tangential velocity by RANS model is well comparable with LES. However, in the decay phase of turbulence, RANS model shows slightly faster decay of tangential velocity due to its slower decay of turbulence compared to LES.

Figure 12 

Time decay of maximum tangential velocity calculated by RANS and LES.

3.2.2. Secondary Currents

The secondary currents calculated by both the models are shown in Figure 13. The flow vectors of secondary current show three circulation cells in each side of the centerline. Secondary current inside the vortex core is found relatively weaker than the outside of core. It is explained that (in Figure 1) the Rankine vortex is a combination of nearly solid body rotation and potential vortex. For both type of vortices, the individual patterns of secondary currents are reported by Laugt [12]. The pattern of secondary currents predicted by LES as well as RANS is found well agreed to Laugt. In the vortex core, the secondary current shows the pattern similar to that of solid body rotation; outside the core, the pattern is similar to that of potential vortex. The patterns as well as magnitudes of secondary currents predicted by two models are well comparable to each other.

(a)

(b)

Figure 13 

Pattern of secondary current at 𝑡 = 40 sec by (a) RANS and (b) LES simulations.

4. Conclusions

3D numerical simulations were carried out for an ideal Rankine vortex using both RANS (non-linear 𝑘-𝜀) and LES models. In the growth phase as well as in stabilized phase of turbulence, the decay rate of tangential velocity by RANS model was found to be well comparable with LES. However, in the decay phase of turbulence, RANS model showed faster decay of tangential velocity compared to LES due to its slower decay of turbulence. The patterns as well as magnitudes of secondary currents predicted by two models are in good agreement with each other. The water surface was found to be depressed near the central area of the vortex with a minimum water depth observed at the centre of vortex.

In order to investigate the predictability of standard (linear) and non-linear 𝑘-𝜀 models, the decay of a trailing vortex was simulated and compared with previous DNS data. The comparison of non-linear 𝑘-𝜀 model prediction with previous DNS result showed similar phenomenon as observed in the comparison with LES as described previously. Due to the inherent deficiency of standard 𝑘-𝜀 model, it predicted much faster decay than that of non-linear model.