Simulation of Swirling Flow with a Vortex Breakdown Using Modified Shear Stress Transport Model

The well-known Shear Stress Transport (SST k–ω) turbulence model was modified and examined. Two industrially relevant problems with curved and rotating channels have been selected to assess the mod...


■ INTRODUCTION
One of the most remarkable drawbacks of the eddy viscosity models (EVMs) is their inability to capture the flow behavior in a rotating system and streamline curvature. For example, for complex cases of fluidic devices without moving parts where the flow demonstrates a high level of rotational motion, EVMs fail to predict these effects accurately. The main reason for this deficiency is due to the use of the Boussinesq hypothesis, which assumes the eddy viscosity as an isotropic scalar, which is incorrect for complex flows with a rotating system. Since 1980, various attempts of modifications and improvements on different EVMs have been proposed to enhance these models (see refs 1−3). However, these models are still inadequate to address three-dimensional (3D) flows. Later, in 1997, Spalart and Shur 4 proposed an empirical approach on the sensitization of the Spalart−Allmaras turbulence model to rotation and curvature effects. This model demonstrated its valuable ability to capture a high level of turbulence. Unlike the approach proposed by Knight and Saffman, 5 Spalart and Shur's model is more efficient, because it measures the extra influence of the invariant contributor to the turbulence. By 2009, Smirnov and Menter 6 had adapted the rotation−curvature correction function proposed earlier in Spalart and Shur 4 to the shear stress transport k−ω model (SST k−ω). The correction function was applied to the production term in both the k and ω transport equations. As a result, the corrected model, denoted as the SSTCC, is more accurate, computationally efficient, and robust than its predecessor.
In this work, the Richardson number (Ri) has been implemented to avoid the need to calculate the Lagrangian derivatives term (DS ij /Dt), which appears in the nondimensional argument r. This leads to a simpler version of the SSTCC, which is realized by implementing the Ri number. As a result, the obtained numerical code requires lower computational cost and is also competitive, in terms of accuracy, when compared to the conventional EVMs. The new formula for r̃is applied into the rotation function developed in Smirnov and Menter. 6 The initial attempt to investigate the performance and accuracy of the developed SSTCCM model has been discussed in the report by Alahmadi and Nowakowski. 7 The work of Alahmadi and Nowakowski was predominantly focused on a specific case of cyclone operation and therefore was restricted, with regard to thorough analyses of the model. The objective of this contribution is to investigate the performance of the SSTCCM model in a range of different physical scenarios that include a vortex breakdown and to compare it with other available models. This paper is organized as follows. The formulation of the proposed SSTCCM model is presented first. Then, it is applied to the simulation of two cases: (1) a rotating lid in a confined cylinder and (2) turbulent swirling flow through an abrupt expansion. For each test case, a description of the problem, numerical methods, and the simulation results are presented and discussed in details. Finally, major findings and conclusions of the paper are summarized.

■ MODEL FORMULATION
The Navier−Stokes equations describe the physics of fluid motion. The continuity and momentum equation respectively represent the mass conservation equation and incompressible form of Navier−Stokes equations as follows: where u and p represent the velocity and pressure. The term s ij is the strain-rate tensor, which is written as The flow variables in eqs 1 and 2 are decomposed into mean and fluctuating quantities. The time average process then is applied, which yields the final form of the Reynolds-averaged Navier−Stokes equations, which can be written as follows: The term S ij is the time-averaged quantity of the strain-rate tensor. The appearance of the additional term τ ij (which is defined as uu ij i j τ ρ =− ′′ ) is known as the Reynolds stress tensor.
In order to find a solution to the mathematical problem, the system of equations must be closed by evaluating a specific term of τ ij . This is the main principle that must be followed when designing any turbulence model.
Menter SST Turbulence Model. In 1994, Menter 8 introduced the shear stress transport (SST) model, which is an improved version of the two-equation k−ω model. The SST model combines the two turbulence models (k−ω and k−ϵ). The k−ω model is used in the inner part of the boundary layer and switches to the k−ϵ model in the free streamflow. The SST model employs the Boussinesq hypothesis to relate the Reynolds stresses to the mean rate of deformation as follows: The conservation form of the transport equations for both the turbulence kinetic energy (k) and the turbulent frequency (ω) can be written as The production term P k is defined as where the turbulent viscosity is defined as follows: The blinding function F 1 in the free stream is zero (k−ϵ model) and, in the boundary layer, is equal to one (k−ω model), given by F tanh(arg ) 1 1 The constant closure coefficients of the SST model are given in Table 1.
Spalart−Allmaras Rotation Curvature (SARC) Turbulence Model. The standard Spalart−Allmaras (SA) turbulence model does not account for the effects of system rotation and streamline curvature. Spalart and Shur 4 proposed an empirical alteration to the original standard SA model. As a result, the production term of turbulence viscosity equation is multiplied by a rotation function f r1 , which is given by The formula contains the dimensionless quantities r* and r, which are defined as The term DS ij /Dt in eq 18 represents the components of the Lagrangian derivatives of the strain tensor, and ε imn is the Levi− Civita symbol. The empirical constants c r1 , c r2 , and c r3 that appear in eq 16 are given the values 1.0, 2.0, and 1.0, respectively.
The function ( f r1 ) enhances the performance of the SA model. 9 However, it is rather complex and increases the computational cost. This is due to the factor accounting for the rotation-curvature correction that contains Lagrangian derivative DS ij /Dt and a higher-order derivative of the strain-rate tensor term D 4 in eq 18.
SSTCCM Model. The model presented in this paper builds on the shear stress transport model (SST k−ω) and also on the rotation function that was proposed in the work of Smirnov and Menter. 6 The Richardson number Ri, which was defined by Hellsten, 10 is used to avoid calculating the terms D 4 and the complex Lagrangian derivatives in eq 18. To control the production term P k that appears in the transport equations, the modified rotation function is used. The incompressible form of the transport equations of the turbulent kinetic energy and its specific dissipation for the SSTCCM model can be cast in differential form: The production term P k and the blinding function F 1 are given in eqs 9 and 11, respectively. The rotation function (f rot )i s defined as The term f r1 is the rotation function presented in eq 16 and has been defined by Spalart and Shur. 4 The production term that appears in the SA model is based on the vorticity tensor (Ω), which characterizes the rotation, while the production term P k is based on the strain-tensor rate S, which characterizes the total deformation. Therefore, P k is greater than S, which justifies the use of the limiter in the rotating function f rot . 6 The dimensionless quantities r, which exactly represents the Richardson number Ri,a sd e fined by Hellsten, 10 and r* are then given as follows: The Ri term includes the effects of rotation and the significance of the streamline curvature effect by accounting for the mean flow deformation. The direct reference to Ri avoids computing thre Lagrangian derivative DS ij /Dt and a higherorder derivative of the strain-rate tensor term D 4 in eq 18. Here, the strain-tensor rate S and the vorticity tensor Ω represent the following expressions: Numerical Implementation and Simulation Results. The described model has been implemented using the open source code OpenFOAM-2.4.x with flexible and extendable C++ libraries. 11 Finite volume grids have been used to calculate the results for all near-wall flows, resolving the viscous sublayer with y + < 1. To avoid any grid resolution uncertainties, highly refined grids have been created for all two-dimensional (2D) and threedimensional (3D) simulations. The model was first tested for a confined cylinder with a rotating lid and verified with respect to experimental data using different numerical schemes. After that, the 3D swirling flow through a sudden expansion pipe was examined using the proposed model and compared with the classical SST model 8

and various conventional EVMs.
Rotating Lid in a Confined Cylinder. Case Geometry and Numerical Settings. The first test case that was chosen to validate the proposed model and compared to experimental data is a confined cylinder with a rotating lid, which has been studied by Fujimura et al. 12 This case constituted a criterion for the choice of the discretization scheme for the momentum equation convective term. The computational domain is shown in Figure  1a, which can be described as a confined cylinder, where the swirling flow is generated by a rotating lid at the top end-wall of the pipe. Because of the shear forces, the flow in contact with the rotating lid acquires a spinning motion. Consequently, the fluid particles are propelled radially outward, while near the sidewall, the fluid spirals down. As it reaches the bottom end-wall, the fluid reverses its direction and moves from the sidewall toward the central axis before spiraling upward. The characteristics of the flow pattern of this swirling motion are mainly dependent on two dimensionless parameters, i.e., the Reynolds number Re (which is defined as Re = R 2 Ω/ν) and the cylinder aspect ratio Ar (which is defined as Ar = H/R), where R is the radius of the cylinder, Ω is the angular velocity of the lid, and H is the height of the cylinder. An axisymmetric 2D simulation is performed in order to facilitate the computational effort. The 2D axisymmetric computational domain is shown in Figure 1b. The numerical results are presented for aspect ratios of 1.5 and 2.5 with the corresponding Re values of 1010, 1290, and 2200, respectively. The relevant meshes are presented in Figure 2. The mesh is refined according to the Re value with a y + value less than unity for all cases: where Δy 1 is the height of the first cell adjacent to the wall and u τ is the frictional velocity. The relationship between the parameters (Re and Ar) and vortex breakdown observations were reported by Escudier. 13 In this case study, three different regimes are considered. These area, which have been labeled Cases A, B, and C, correspond to the following values: Re = 1290 and Ar = 1.5, Re = 2200 and Ar = 2.5, and Re = 1010 and Ar = 2.5, respectively. The equations of the turbulence quantity variables (k and ω) are both discretized using the second-order linear upwind scheme (LU). The equations residuals for the investigated firstand second-order schemes were set to fall below a predetermined residual tolerance value over a successive number of iterations. For the pressure equation, the residual value was 10 −8 , whereas, for all other variables, a value of 10 −6 was considered sufficient. Figure 3 shows that the upwind scheme provides an erroneous streamlines for cases A and B. On the other hand, the secondorder schemes (Self Filtered Central Differencing Scheme (SFCD) and LU) are both qualitatively in good agreement with the experimental measurements. Unlike the upwind scheme, the location and the number of the vortex breakdowns are accurately predicted by the second-order schemes. Figure 4 shows that the SSTCCM model provides good predictions of the axial velocity profile along the longitudinal centerline, when compared to experimental data. In terms of accuracy and reliability, the second-order schemes are superior to the simple first-order scheme. For example, Table 2 shows that the upwind scheme for the fine mesh underestimates the maximum velocity by ∼10% for case A and by <1% for cases B and C. Although the upwind scheme seems to provide accurate predictions of the velocity profile, it cannot capture the formation of the vortex breakdown, as shown in Figure 3. The results of the second-order schemes demonstrate that the SFCD scheme provides insignificant improvement, in terms of the grid refinement. In contrast, the accuracy of the LU scheme is closely related to the grid refinement (see Tables 3 and 4). Therefore, the second-order LU scheme will be adopted for the sudden expansion case.

■ RESULTS AND DISCUSSION
The 3D computations should resolve other flow features as geometry induces more intricate dynamics and the flow field is influenced by the stretching of vortical structures. The mechanism of vortex stretching is absent in 2D flows; therefore, the results could be qualitatively different from 3D results. There are also restrictions related to k−ω turbulence models. Simpson and Ranade 14 highlighted excessive predictions of turbulence kinetic energy in the vicinity of stagnation points, as well as the differences between 2D axisymmetric models and fully 3D approaches.
Swirling Flow through a 3D Sudden Expansion Pipe. Turbulent swirling flow through a 3D abrupt expansion can be found in various fluidic devices without moving parts. Examples are hydrocylones, 15−17 cyclone separators, 18 jet ejectors, 19 or spray dryers. 20 The flow field is complex. It includes various dynamic phenomena that are manifested by extremely high levels of turbulence, recirculation, separation and reattachment as well as vortex breakdown. The initiation and determination of the location of the vortex breakdown and the flow instabilities can be construed because of the diverging nature of the expansion flow. Therefore, in order to predict this flow behavior correctly, particular attention must be paid to the choice of the turbulence closure. Despite of the numerous studies that have investigated vortex breakdown through the sudden expansion pipe, including refs 21−25, this type of flow is still not fully understood, particularly for high swirling flow with swirl numbers near unity.
The present work evaluates the capability of the new formulation by comparing its performance with different EVMs: standard k−ϵ; RNG k−ϵ, and SST k−ω. It then validates the results, using the experimental measurements of Dellenback et al. 21 Case Geometry and Numerical Settings. The numerical configuration was based on the experiment of Dellenback et al. 21 The experiment was designed to investigate the turbulent swirling flow through an abrupt expansion with an aspect ratio of  Industrial & Engineering Chemistry Research pubs.acs.org/IECR Article 1:1.94, using the laser Doppler anemometer LDA method. Figure 5 and Table 5 depict the computational domains used in the present study. The structured grid consists of hexahedral elements. A high-resolution discretization is applied close to the wall with a y + value of less than unity. The simulation was performed first on a course mesh, then on refined mesh that was obtained by doubling the number of grid points in each direction. This resulted in a negligible difference for all of the primitive variables. The final mesh that was implemented, which is depicted in Figure 6, consisted of 1 567 954 cells.
The Reynolds number, based on the upstream diameter D,is Re =3× 10 4 . The inlet section of the flow is placed at x = −2D, a n da na x i s y m m e t r i cs u d d e ne x p a n s i o nd o w n s t r e a mi s encountered at x = 0. Two cases for the inlet velocity has been set: either a pure axial velocity of (0.473), as stated in Table 6,or a swirling inlet boundary condition that is adjusted as an axisymmetric profile. The axial and tangential components of the inlet velocity profile are adapted from the experimental work of Dellenback et al., 21 while the radial velocity is set to zero. The For the swirling case, the swirl number is ∼0.6, based on the inlet radius (R = D/2); υ t and υ z indicate the time-averaged tangential and axial velocities, respectively. Figure 7 shows the applied velocity profiles at the inlet for the swirling case (S = 0.6).
The transient PisoFOAM solver is employed to handle the pressure−velocity coupling. For the temporal differentiation, the second-order backward scheme is applied. All simulations were run with a constant time step of (2 × 10 −4 ). The minimum required time step was calculated using the Courant−Friedrichs-Lewy (CFL) stability condition. The considered EVMs were   For the pressure equation, the residual value was 10 −8 , whereas, for all other variables, a value of 10 −6 was considered to be sufficient. Characteristic profiles representing primitive variables values were also monitored, to ensure that these are converged.

■ RESULTS AND DISCUSSION
The first case study is the no-swirling flow at S =0 . Figure 8 compares the axial velocity contours predicted by the standard k−ϵ, RNG k−ϵ, and the SST k−ω models with the new SSTCCM model. It can be seen that the standard k−ϵ model predicts a very weak recirculation zone, with a percentage error of 63.9% in predicting the reattachment position, compared to the experimental value (X r /h = 9.3). This could be attributed to the isotropic nature of this model. The RNG k−ϵ model performs better than the former, generating 39.3% error, but still not as required. The congruence between the results of both the SST k−ω and SSTCCM models and the experimental data was manifested by a smaller percentage error of 25.8%. Figure 9 shows that the SSTCCM model performs well, with regard to predicting axial mean velocity profiles along the longitudinal axis when validating it against the experimental results of Dellenback et al. 21 At a moderate to a high level of swirling strengths, the flow field starts exhibiting an indistinct and complex phenomenon of an unsteady asymmetry that is usually observed in swirling flows.
This feature is initiated with the on-axis recirculation and the vortex breakdown phenomenon. Figures 10 and 11 show the extent of the recirculation zone at (S = 0.6) through contours of axial velocity and static pressure. As the swirling level increases to moderate values, the flow has a tendency toward asymmetrical form with central recirculation, signaling the start of vortex breakdown. This is due to the viscous dissipation of the tangential velocity component, which produces an adverse pressure gradient on the centerline of the tube. Figure 12 reassures that the SSTCCM model performs well with regard to predicting axial and tangential mean velocity profiles along the longitudinal axis, when validating it against the experimental results of Dellenback et al. 21 There is a significant challenge in precisely predicting the unsteady characteristics of the swirling flow with vortex breakdown utilizing the conventional EVMs. This weakness is clearly observed in Figure 10, where the standard k−ϵ, RNG k−ϵ, and the SST k−ω models were unable to capture the vortex breakdown phenomenon. The SSTCCM model is better-suited to capture the location and extent of the central recirculation zone in the swirling flow immediately after expansion.  Figure 5. Schematic diagram of the sudden expansion geometry.  The Reynolds Stress Model (RSM) could be treated as an alternative approach for the simulations presented in this study. The model accounts for anisotropic Reynolds stresses in the flow, contrary to the EVMs models, which do not provide the directional dependence of the turbulent stresses. The model explicitly accounts for the effects of streamline curvature and rotation. Although the RSM approach is superior to any EVM approaches and has been successfully used for computing swirling flows 26,27 (for example, to resolve secondary flows in gas vortex units), 28 it has been found that the model is incapable of resolving all of the deficiencies of the two-equation models for simulating turbulent swirling flows. Tsai et al. 29 compared the RSM simulations with measurements and reported that, for weakly swirling flows (S = 0.3), the strength of the decay of swirlinduced deceleration of the axial velocity is not reproduced correctly. Lu and Semiao 30 reported that, for this flow regime, the RSM model does not resolve the intensity of turbulence along the center line correctly. The RSM model also carries significant computations overhead, because the additional equations for the Reynolds stresses must be solved in three dimensions and, therefore, the model was not considered in the present study. On the other hand, the common deficiency of the closure model modifications is related to the fact that a scalar formula does not distinguish directional components of the covariance tensor of Reynolds stresses. To elevate this problem, various other modifications of a closure model have been applied. The turbulent production limiters meant to avoid excessive turbulent kinetic energy prediction in stagnation regions were proposed by Menter 8 and by Kato and Launder. 31 The production limiters replace the quadratic strain contribution with the strain rate multiplied by the vorticity rate. In other methods, the original Kato−Launder approach was replaced with a linear combination of the unmodified production term and the Kato−Launder modified production term. Reboud et al. 32 have proposed the introduction of a production limiting term, which is a simple arbitrary multiplier acting on the turbulence viscosity term.

■ CONCLUSIONS
The numerical simulation have been performed using the developed shear stress transport model with curvature correction modification. The results were compared with other popular turbulence closure models in flow regimes that are subject to directional forces, such as those due to system rotation. The model was successfully validated against experimental data representing the flow in a confined cylinder with a rotating lid and the 3D swirling flow through a sudden expansion pipe. The model that was designed for practical simulations of complex cases requiring a quick turnaround time to complete calculations proved to reproduce the flow profiles with better accuracy than the conventional EVMs. The proposed approach allowed capturing the location of the vortex breakdown phenomenon in the swirling flow.