Improved Dissipation Rate Equation for Swirling and Rotating Pipe Flows

Ronald Ming Cho So

Department of Mechanical Engineering, The Hong Kong Polytechnic University, Hong Kong, China.

DOI: 10.4236/jamp.2022.103048                

Abstract

The nature of turbulent swirling and rotating flow in a straight pipe is investigated using a family of near-wall two-equation models. Specifically, the viability of three different near-wall two-equation models is assessed. These models are asymptotically consistent near the wall. The first two models, one with isotropic and another with anisotropic eddy viscosity invoked, solved a dissipation rate equation with no explicit correction made to account for swirl and flow rotation. The third model assumes an isotropic eddy viscosity but solves an improved dissipation rate equation that has explicit corrections made to account for swirl and flow rotation. Calculations of turbulent flows in the swirl number range 0.25 - 1.3 with and without a central recirculation region reveal that, with the exception of the third model, neither one of the other two models can replicate the mean field at the swirl numbers tested. Furthermore, taking stress anisotropy into account also fails to model swirl effect correctly. Significant improvements can be realized from the third model, which is based on an improved dissipation rate equation and the assumption of isotropic eddy viscosity. The predicted mean flow and turbulence statistics correlate well with measurements at low swirl. At high swirl, the two-equation model with an improved dissipation rate equation can still be used to model swirling and rotating pipe flows with a central recirculation region. However, its simulation of flows without a central recirculation region is not as satisfactory.

Keywords

Dissipation Rate Equation, Swirling and Rotating Pipe Flows

1. Introduction

1.1. Motivation

In order to avoid the necessity of having to incorporate a fuel compartment in ramjet design, a solid-fuel combustor that integrates fuel storage into the combustor is usually adopted. However, ramjet rotation (Figure 1(a)) could complicate the complex reacting flow downstream of the combustor entrance (Figure 1(b)) and could give rise to premature combustion extinction. Since spin is common in ramjet operation, it could impact combustor performance and ramjet trajectory in the duration of its flight. Ahmed et al. ref. [1] have examined the effect of spin on the flow alone inside a ramjet combustor in the laboratory. From this experiment, a large decrease in combustion efficiency could be inferred from the fluid dynamic study. The combustion efficiency decrease could be as large as 50% or more depending on the solid fuel used. Therefore, it could be speculated that spin is a major cause of this decrease.

To shed light on this reduction, it is necessary to understand the combusting flow behaviour inside the combustor. Once the oxidant air enters the combustor, it immediately encounters a sudden expansion, thus, creating a recirculation region and a flow reattachment (Figure 1(a)). Downstream of the reattachment

(a) (b)

Figure 1. (a) Flow behavior downstream of entrance nozzle in a solid fuel ramjet combustor; (b) Details of flow field downstream of recirculation zone. The flow downstream of nozzle entrance is made up of four regions: (1) Air region; (2) Air plus combustion products region; (3) Flame region; (4) Fuel vapor plus combustion products region. Density differences in these four regions give rise to stratification effect in the combusting flow.

point, a wall jet and a boundary-layer type flow will develop (Figure 1(b)) and they will last until the flow reaches the exit expansion nozzle. The ramjet combustor is normally designed in such a way that the recirculation region serves to anchor the flame (Figure 1(a)) and will probably occupy about 20% of the combustor length, while the wall jet, the boundary-layer and the pipe flow that follows dominate the other 80%. The recirculation region is necessary for stable combustion to occur inside the combustor. Therefore, in order to maintain continuous combustion with good efficiency, a thorough understanding of the flow dynamics inside the combustor is necessary.

Besides the backstep flow, there is also a plane wall jet downstream of the reattachment point, and a developing boundary layer as the flow slowly evolves into a pipe flow in its passage to the exit nozzle (Figure 1(b)). These complexities are further complicated by combustion that creates density stratification, and ramjet spin that creates swirl in the combusting flow. In order to achieve optimal design of the combustor, a good understanding of the isothermal flow complexities inside the ramjet is of utmost importance. If the turbulent flow inside the ramjet were to be simulated, a relatively simple turbulence model, such as a two-equation model, is preferred because it could be easily extended to tackle turbulent combusting flow. Furthermore, if the complex flow were to be correctly modelled by relatively simple turbulence models, such as two-equation models, then an accurate modelling of the individual complex flow by a two-equation model needs to be demonstrated and established. Successful modelling of backstep flow, plane wall jets, and flow with density stratification using two-equation models have been demonstrated by So et al. ref. [2], by Gerodimos and So ref. [3], by Sommer et al. ref. [4], and by Yoo and So ref. [5], respectively, while their suitability for rotating pipe flows has been shown to be viable by Yoo et al. ref. [6], who made comparisons of two-equation simulation results with those obtained using Reynolds-stress models. Furthermore, the Yoo et al. ref. [6] study also showed that two-equation models need improvements, especially in their viability and suitability for modelling swirl and rotating flow in ramjet combustor simulation.

1.2. Characteristics of Swirling Flow

Swirl decay in a straight pipe can be classified into two different types depending on the way swirl is imparted into the flow. The swirl in the first type is created by a rotating pipe while the second type is created by a swirler installed at the entrance of the pipe. A forced vortex or a solid-body rotation is generated in the first type, while the vortex could be a free vortex, a forced vortex, or a combined free/forced vortex in the second flow type. Consequently, the resultant swirling flow behaves differently in the two flow types.

In the creation of a rotating pipe flow, the flow could either pass through a continuously rotating pipe, such as that reported by Murakami and Kikuyama ref. [7] and by Kikuyama et al. ref. [8], or through a rotating pipe into a stationary pipe where the swirl velocity decays as the flow moves downstream, such as the Weske and Sturov ref. [9] experiment, and the investigation of Anwer and So ref. [10]. The counterpart to a concentric annulus has been examined by Hirai et al. ref. [11]. If the term “destabilizing” is used to describe the phenomenon of extra turbulence production created by flow rotation, while the term “stabilizing” is used to describe the rotating flow characteristics without the extra turbulence production present, then flow regions with extra turbulence production present (destabilizing) or absent (stabilizing) could be found in different parts of an axially rotating pipe flow. Near the entrance of a rotating pipe, the wall boundary layer is very thin. Consequently, fluid rotation has to decrease rapidly from pipe rotation at the wall to zero outside of the wall boundary layer. The flow near the wall is influenced by a very high mean tangential shear strain and turbulence production is greatly enhanced.

In the entrance region of the rotating pipe, the near-wall flow is destabilized by rotation. Far downstream, the flow becomes fully developed and the fluid rotates as a solid body. Since the solid-body rotation curve has a constant slope, turbulence production due to the tangential mean strain is essentially zero. In this downstream region, rotation gives rise to a stabilizing effect on the flow. Between the entrance and the fully developed region, the flow is influenced by both destabilizing and stabilizing forces created by the rotating pipe. However, the transition from destabilizing to mixed to stabilizing occurs over a finite length and depends on different factors; among them are inlet flow conditions, the flow Reynolds number, Re, and the swirl number, S. In this case, S can be approximated by the ratio of the pipe rotation velocity to the mean bulk velocity of the pipe, or S = ωD/2U_m. This three-dimensional flow has also been investigated and reported by Murakami and Kikuyama ref. [7], and Kikuyama et al. ref. [8].

In the experiments carried out by Weske and Sturov ref. [9] and by Anwer and So ref. [10], the flow very near the wall was immediately affected by a large mean tangential shear strain due to the sudden relaxation of surface rotation. Therefore, the flow in the entrance region of their rotating pipe flow is influenced by both stabilizing and destabilizing effects. Eventually, the solid-body rotation would decay completely, and a fully developed pipe flow would result. Therefore, the flow characteristics downstream of the rotating section would depend on the inlet conditions, on Re and on S. Weske and Sturov ref. [9] measured the Reynolds stresses and the mean flow for two different S and they found that, even at 150D downstream, the turbulence field has not completely recovered its stationary behavior. On the other hand, the experiments of Anwer and So ref. [9] were carried out to study the effects of S on wall shear only and did not examine the decay phenomenon. Their measurements were obtained at 2D downstream of the entrance. From their data, it is clear that the destabilizing effect, although only present in about 10% of the pipe diameter, dominates over the stabilizing effect, thus promoting turbulence production along the pipe.

Swirling flow, with either a free or a forced vortex generated by swirl generators, in a straight pipe has not been investigated extensively in the past. In this case, the free or forced vortex evolves in the downstream direction into a combined free and forced vortex before decaying. Among the more relevant studies, the work of Kreith and Sonju ref. [12], Backshall and Landis ref. [13], Yajnik and Subbaiah ref. [14], Murakami et al. ref. [15], Padmanabhan and Janek ref. [16], Ito et al. ref. [17], Kito ref. [18], Kito and Kato ref. [19], Algifri et al. ref. [20], Kitoh ref. [21] and Parchen et al. ref. [22] could be mentioned. The earlier studies up to 1984, including the axisymmetric nature of swirling flow examined by Kito ref. [18], were mainly concerned with the mean field, the pressure drop, and the decay of mean swirl in a pipe. Little attention has been paid to the decay of the turbulence field. On the other hand, some measurements of the turbulence field evolution have been provided by Algifri et al. ref. [20], Kitoh ref. [21] and Parchen et al. ref. [22]. As a result, turbulence structure in the core and the outer region is found to be quite different from those in the near wall region. This difference might be difficult to replicate using conventional turbulence models for swirling flow simulations.

The studies mentioned above cover a wide range of S. Depending on the type of swirler used, the free-to-forced portion of the vortex will vary. The flow characteristics near the pipe centerline also differ significantly depending on S. For flows with S $\ll$ 1, a recirculation region did not exist near the pipe core because the resulting pressure depression is not sufficient to create such a reversed flow. For flows with S ≈ 1, a reversed flow region starts to appear depending on whether other conditions, such as Re, inlet conditions, etc., are favorable to the formation of a recirculation region. For example, the experiments of Kitoh ref. [21] showed a reversed flow region while the measurements of Murakami et al. ref. [15] did not, even though both studies were carried out at an initial S ≥ 1. The decay rate was also found to be dependent on the above conditions. For example, in the experiments of Murakami et al. ref. [15], the mean flow was still evolving even at 190D downstream of the inlet. On the other hand, Weske and Sturov ref. [9] found that the mean flow was quite similar to that of a fully developed turbulent pipe flow at about 150D downstream of the inlet, while the turbulence field was still evolving.

1.3. Modelling Swirling Flow

Due to the complex nature of swirling flow and its decay in a straight pipe, it’s modelling has not been investigated extensively. Among the more notable studies are those of Kreith and Sonju ref. [12], Parchen et al. ref. [22], Kobayashi and Yoda ref. [23], and Khodadadi and Vlachos ref. [24]. Three different approaches to model swirl decay were discussed in these studies. In all but one investigation, wall functions were assumed. An analytical approach was adopted by Kreith and Sonju ref. [12] who assumed a constant eddy viscosity and an axial flow given by the mean axial velocity in a fully developed pipe flow. Thus simplified, the equation governing the tangential velocity, W, can be reduced to an equation with one unknown once the eddy viscosity is specified. Using their own experimental data, they were able to derive an eddy viscosity that was parametric in Re and used it to solve for the evolution of W. They were able to derive an exponential decay law for S. The results thus obtained were in good agreement with their own experimental measurements. Two-equation k-ε models and their variations had been deployed by Kobayashi and Yoda ref. [23] and Khodadadi and Vlachos ref. [24] to simulate swirl decay in a straight pipe. These calculations indicated that the k-ε model was, in general, inadequate and gave poor prediction of the flow in the pipe core region, especially the flow depression along the axis. The coefficients and model constants of the k-ε model could be modified to give better correlations with data as suggested by Kobayashi and Yoda [23]. However, such modifications were rather ad hoc and could not be easily generalized to flows other than those analyzed by Kobayashi and Yoda ref. [23], and Khodadadi and Vlachos ref. [24].

An algebraic stress model plus the standard k and ε equations were adopted by Parchen et al. ref. [22] to simulate swirl decay in a straight pipe. They also compared their results with standard k-ε model calculations and other model predictions obtained by their colleagues. The results showed that, as far as swirl decay was concerned, the standard k-ε model gave a better prediction than the algebraic stress model. However, both models failed to capture the behavior of the mean flow correctly as swirl decayed in the downstream direction. Specifically, the standard k-ε model failed to predict any enhancement of radial momentum exchange as a result of swirl. On the other hand, the algebraic stress model predicted a weak enhancement. Also, for flows with small S, the standard k-ε model performed the best among all models tested. These results were quite puzzling and could be partially explained by attributing the discrepancies to the use of wall functions in the modelling. Consequently, near-wall diffusive transport of all stress components was under-estimated; thus, it could lead to a situation where error cancels out error and give rise to a better prediction by the k-ε model.

These studies led to the realization that not a single k-ε model, without some kind of ad hoc modifications, could correctly track the swirling flow evolution in a straight pipe. Therefore, an alternative approach to model swirling flow is necessary. According to Bardina et al. ref. [25], who studied rotation effects in isotropic turbulence, their findings indicate that the effects of flow rotation and/or swirl should be accounted for in the ε-equation. They suggest modelling the ε-equation to take swirl effects into account directly and explicitly. In addition, Fu et al. ref. [26] and Kim and Rhode ref. [27] showed that using algebraic and second-moment closure models could greatly improve swirl flow simulations. Their model results were in good agreement with rotating flow data. On the other hand, Bernard and Speziale ref. [28] and Speziale and Bernard ref. [29] have found that, if homogeneous shear flow and its asymptotic state were to be calculated correctly, an explicit term representing the effect of vortex stretching should be included in the modified ε-equation.

The two-equation models and their improvements discussed above cannot fulfill the need for a general two-equation model that could simutaneously replicate backstep flow, wall jet, thermal boundary-layer flow, developing pipe flow, and swirling and rotating pipe flow, such as those commonly encounter in spinning ramjet combustors. Recent studies on swirling flow simulations reported by Andrew and Cui ref. [30], Sagol et al. ref. [31], Javadi ref. [32], Li et al. ref. [33] and Alhmadi et al. ref. [34], to mention a few recent studies, also did not emphasize their ability to model the highly complex swirling flow found in spinning ramjet combustors. Similar to other swirling flow models, theirs could also accurately mimic swirling flow without the presence of added flow complexities found in spinning ramjet combustor. Therefore, the need to develop a turbulence model that could properly handle these complexities still has not been fulfilled.

1.4. Present Objectives

Based on the above review on swirling flow modelling, it can be concluded that the time is ripped for the development of a relatively more general swirling flow model based on the two-equation model used in the studies of So et al. ref. [2], Gerodimos and So ref. [3], Sommer et al. ref. [4], Yoo and So ref. [5], and Yoo et al. ref. [6]. Thus configured, the new model could better handle flow complexities found in spinning ramjet combustors. In addition, the swirling flow modelling studies quoted above suggest that the conventional ε-equation also needs modification if the effects of vortex stretching and rotation also were to be simulated. These studies further suggest that the ideas of Bardina et al. ref. [25], Bernard and Speziale ref. [28] and Speziale and Bernard ref. [29] could be used to derive a new near-wall ε-equation. The resultant ε-equation together with a near-wall k-equation, such as that adopted by So et al. ref. [2] to treat backstep flow, by Gerodimos and So ref. [3] to model plane wall jets, by Sommer et al. ref. [4] to treat flow with heat transfer, and by Yoo and So ref. [5] to simulate density stratified flows, could give rise to an improved k-ε model.

The approach used to accomplish this objective is similar to that adopted by So et al. ref. [35] to account for low-Reynolds-number and near-wall turbulence effects in the ε-equation in order to correctly simulate such flows. Therefore, in the present approach, the ε-equation is modified to account for swirl and pipe rotation, while the modified k-equation adopted in previous studies refs. [2] [3] [4] [5] [6] is used without implementing further modelling change to account for swirl effects. Thus formulated, the k-equation and the modified ε-equation together will give an improved k-ε model that is suitable for modelling complex turbulent flows such as those found in spinning solid-fuel ramjet combustor.

The improved k-ε model is validated against swirling and rotating pipe flow measurements and other two-equation k-ε modelling results. In particular, the new model with isotropic eddy viscosity invoked is assessed against the simulations of two other conventional models that adopt the conventional k-ε equations. One model retains the isotropic eddy viscosity assumption, while another invokes anisotropic eddy viscosity. This way, the effects of stress anisotropy on swirling flow simulations can also be assessed; thus allowing the merits or lack thereof of the explicit modelling of swirl effects in the ε-equation to be thoroughly studied in this investigation.

2. Two-Equation Models

In order to avoid the use of wall function approximations, the present study chooses to apply near-wall models to calculate swirling flows in a straight pipe. The advantage of this approach has been fully demonstrated in the modelling studies of back-step flow, plane wall jet, and boundary-layer flow, refs. [2] [3]. Therefore, three different near-wall two-equation models are investigated in this study. Two of the models share the same set of k and ε equations while a third has a totally different ε-equation, which is derived to specifically account for the effects of rotation and vortex stretching. The first two models are designated as SAYS and ASAYS and they differ only in the eddy viscosity assumed. The SAYS model proposed by So et al. ref. [36] assumes an isotropic eddy viscosity and is based on the SSG model formulated by Speziale et al. ref. [37] without near-wall correction, while the ASAYS model adopts the anisotropic eddy viscosity put forward by Speziale ref. [38]. The SAYS and ASAYS models have been used by So et al. refs. [2] and Gerodimos and So ref. [3], respectfully, to tackle backstep flow and plane wall jet. Therefore, both SAYS and ASAYS models can account for near-wall effects, which consists of both low-Reynolds-number and wall blocking effects but cannot predict the effects of low-Reynolds-number turbulence in a free shear flow. A third model formulated to also take this low-Re effect into account so that it is valid for both free and wall-bounded shear flows has been put forward by So et al. ref. [36]. This new model is designated as the modified SSGZ model. In order to accomplish the present objective, it is necessary to directly model the effects created by swirl and rotating pipe flow in the already improved ε-equation adopted in the SSGZ model. It should be pointed out that this SSGZ model also invokes the assumption of isotropic eddy viscosity. These three models, SAYS, ASAYS, and the modified SSGZ, are briefly described below.

2.1. Isotropic Eddy Viscosity Model—SAYS

Consider an incompressible turbulent swirling flow in a straight pipe. The mean flow and modeled k-ε equations with gradient transport assumed can be written in Cartesian tensor form as

$\frac{\partial U_i}{\partial x_i}=0$,(1)

$\frac{D{U}_i}{Dt}=-\frac{1}{\rho }\frac{\partial P}{\partial x_i}+\nu \frac{{\partial }^2{U}_i}{\partial x_i\partial x_j}-\frac{\partial \overline{u_i{u}_j}}{\partial x_j}$, (2)

$\frac{Dk}{Dt}=\frac{\partial }{\partial x_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]+{\Pi }_k-\epsilon$, (3)

$\frac{D\epsilon }{Dt}=\frac{\partial }{\partial x_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_i}\right]+C_{\epsilon 1}\frac{\epsilon }{k}{\Pi }_k-C_{\epsilon 2}\frac{\epsilon \tilde{\epsilon }}{k}+\xi$, (4)

where D/Dt is the material derivative, U_i is the i^th component of the mean velocity, x_i is the i^th component of the coordinate, ${\Pi }_k=\overline{u_i{u}_j}\left(\partial U_i/\partial x_j\right)$ is the production of k, P is the mean pressure, ρ is fluid density, t is time, $\tilde{\epsilon }=\epsilon -2\nu {\left(\partial \sqrt{k}/\partial x_i\right)}^2$ is a reduced ε, ν_t is the turbulent eddy viscosity, $-\overline{u_i{u}_j}$ is the Reynolds stress tensor, ξ is a near-wall correction function for the ε-equation and σ_k, σ_ε, C_ε1 and C_ε2 are model constants. Here, the overbar is used to denote time average. The near-wall correction function derived for the SASY model is given by So et al. ref. [35] as

$\xi =f_{w1}\left(-L\frac{\epsilon }{k}{\Pi }_k+M\frac{{\tilde{\epsilon }}^2}{k}-N\frac{\epsilon \tilde{\epsilon }}{k}\right)$, (5)

where $f_{w1}=\exp \left[-{\left(R{e}_t/40\right)}^2\right]$ is a damping function used to ensure the disappearance of the effects of ξ away from the wall, $R{e}_t=k^2/\nu \epsilon$ is the turbulence Reynolds number and the constants are given as L = 1.5, M = 0.5 and N = 0.57. Another reduced ε is introduced and this is given by $\bar{\epsilon }=\epsilon -2\nu k/x_2^2$, where x₂ is the wall normal coordinate. The other model constants are specified as σ_k = 1, σ_ε = 1.45, C_ε1 = 1.50 and C_ε2 = 1.83.

If an isotropic eddy viscosity is assumed, the definitions for $-\overline{u_i{u}_j}$ and the eddy viscosity, ν_t, are given by

$-\overline{u_i{u}_j}=2{\nu }_t{S}_{ij}-\frac{2}{3}k{\delta }_{ij}$, (6a)

${\nu }_t=C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon }$, (6b)

where $S_{ij}=\left(\partial U_i/\partial x_j+\partial U_j/\partial x_i\right)/2$ is the mean strain rate, C_μ = 0.096 and f_μ is a damping function introduced to render proper asymptotic behavior for the Reynolds shear stress. The specification for f_μ is

$f_{\mu }=\left[1+\frac{3.45}{\sqrt{R{e}_t}}\right]\tanh \left(\frac{x_2^{+}}{115}\right)$, (7)

where $x_2^{+}=u_{\tau }{x}_2/\nu$ and $u_{\tau }$ is the friction velocity.

The boundary conditions for swirling flows can be stated as U = V = 0, W = 0 for a stationary pipe and W = ωD/2 for a rotating pipe, k = 0 and

${\epsilon }_w=2\nu {\left(\partial \sqrt{k}/\partial x_2\right)}_w^2$ at x₂ = 0. Here, x₂ is taken to measure positive away from

the wall. Its relation to r, the radial coordinate, is x₂ = D/2 − r (Figure 2). At the symmetry axis, the boundary conditions are zero radial gradient for all quantities except the Reynolds shear stress, which is specified as zero.

Figure 2. Co-ordinate system for swirling and rotating pipe flow.

2.2. Anisotropic Eddy Viscosity Model—ASAYS

The above model assumes gradient transport and isotropic eddy viscosity. Therefore, it cannot account for anisotropy effects. In order to capture this particular characteristic without resorting to a full second-order modelling of the flow, it is proposed to couple an anisotropic eddy viscosity with the above two-equation model to calculate swirling flows. The anisotropic model adopted is that proposed by Speziale ref. [38] and can be written as

$-\overline{u_i{u}_j}=2{\nu }_t{S}_{ij}-\frac{2}{3}k{\delta }_{ij}-4\frac{C_D{\nu }_t^2}{k}\left[S_{ik}{S}_{kj}-\frac{1}{3}{S}_{kl}{S}_{kl}{\delta }_{ij}+{\tilde{S}}_{ij}-\frac{1}{3}{\tilde{S}}_{kk}{\delta }_{ij}\right]$, (8a)

where

${\tilde{S}}_{ij}=\frac{\partial S_{ij}}{\partial t}+U_k\frac{\partial S_{ij}}{\partial x_k}-\frac{\partial U_i}{\partial x_k}{S}_{kj}-\frac{\partial U_j}{\partial x_k}{S}_{ki}$, (8b)

is the frame indifferent Oldroyd derivative of S_ij and C_D = 1.68 is a dimensionless constant. It should be noted that the isotropic eddy viscosity is recovered in the limit of CD goes to zero. This anisotropic eddy viscosity along with Equations (3) to (5) is designated as the ASAYS model. The objectives here are first to demonstrate the ease with which the two-equation model can be used with an anisotropic eddy viscosity, where no other modifications are required of the modeled equations, and second to assess the anisotropy effects on the calculation of turbulent flows with high swirl. Once the validity of the concept has been demonstrated, the anisotropic eddy viscosity can be used with other two-equation models. Of course, alternative anisotropic eddy viscosities can be used. An example is the proposal of Gatski and Speziale ref. [39].

2.3. The Modified SSGZ Model for Swirling Flow—Modified SSGZ

The improved ε-equation for swirling flows is derived in two steps. For the sake of completeness, a brief description of this model is given below. First, the abilities to model free shear flows and their decay as well as near-wall flows are considered. The model thus derived is the SSGZ model reported in So et al. ref. [36]. This is followed by a second modification where the effects of swirl and their modelling are examined. These requirements, therefore, necessitate a complete re-examination of the e-equation and its modelling. However, the form of the modeled equation should remain simple and should not deviate substantially from the established form of a modeled production and a modeled destruction term. In other words, the e-equation should have the same form as that given in Equation (4), but with different constants and damping functions assumed for the production and destruction terms. Furthermore, the near-wall correction function ξ should again be derived to render the improved equation asymptotically correct near a wall.

According to Launder et al. ref. [40], the spreading rate of free shear layers can be predicted correctly if the constants in Equations (3) and (4) are changed to the following set: σ_k = 1, σ_ε = 1.4, C_ε1 = 1.50 and C_ε2 = 1.9. Further, the modelled destruction term $\epsilon \tilde{\epsilon }/k$ should be replaced by ${\epsilon }^2/k$. In addition, as pointed out by Launder et al. ref. [41], the set of constants chosen should yield a correct von Karman constant κ for the log-linear region of wall-bounded turbulent flows. However, such an ε-equation could not yield correct results for the simulation of homogeneous turbulence decay. On the other hand, Hanjalic and Launder ref. [42] found that the decay of homogeneous free turbulence can be correctly simulated if and only if the destruction term in the dissipation rate equation is modified by a damping function that is dependent on Re_t. They proposed the following damping function $f_{\epsilon }=1-0.22\exp \left[-{\left(6R{e}_t\right)}^2\right]$. This suggests that $C_{\epsilon 1}\left(\epsilon /k\right){\Pi }_k$ in Equation (4) should remain the same while $C_{\epsilon 2}\left(\epsilon \tilde{\epsilon }/k\right)$ should be replaced by $C_{\epsilon 2}{f}_{\epsilon }\left({\epsilon }^2/k\right)$. Thus modified, the near-wall correction function ξ can be derived following the procedure outlined in So et al. ref. [43]. The resultant function that yields asymptotically consistent results for k and ε can be incorporated into the ε-equation to render it valid for both low-Reynolds-number turbulence and near-wall flows. The modified ε-equation is given by

$\frac{D\epsilon }{Dt}=\frac{\partial }{\partial x_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_i}\right]+C_{\epsilon 1}{f}_1\frac{\epsilon }{k}{\Pi }_k-C_{\epsilon 2}{f}_2\frac{\epsilon \tilde{\epsilon }}{k}$ (9)

where $f_1=\left[1-\alpha f_{w1}\right]$ and, as before, $f_{w1}$ is given by $f_{w1}=\exp \left[-{\left(R{e}_t/40\right)}^2\right]$ and $f_2=f_{\epsilon }$ while $\tilde{\epsilon }$ is given by

$\tilde{\epsilon }=\epsilon -\frac{C_{\epsilon 3}}{C_{\epsilon 2}{f}_2}\nu {\left[\frac{\partial \sqrt{k}}{\partial x_i}\right]}^2$, (10)

while the constants are determined by So et al. ref. [43] to be $C_{\epsilon 3}=2.955$ and $\alpha =10$. The ε boundary condition would lead to $\tilde{\epsilon }=0$ at the wall and thus rendering (9) valid as the wall is approached. Therefore, the SSGZ model consists of solving Equations (3), (6), (9) and (10). However, the damping function f_μ is given by So et al. ref. [44] as

$f_{\mu }=\left[1+80\exp \left(-R{e}_t\right)\right]{\left[1-\exp \left(-\frac{R{e}_{\epsilon }}{43}-\frac{R{e}_{\epsilon }^2}{330}\right)\right]}^2\left(1+3R{e}_t^{-3/4}\right),$ (11)

where $R{e}_{\epsilon }={\left(\nu \epsilon \right)}^{1/4}{x}_2/\nu$ is a Reynolds number based on the Kolmogorov velocity scale. This model has been validated against a wide variety of flows by So et al. refs. [43] [44] and good correlation with data is obtained for each of the flow case studied.

As will be seen later, the ε-equation given by (9) fails to improve the predictions of swirling flows. The reason is that it does not account for the effects of swirl in an explicit manner. In a study on isotropic turbulence under the influence of a constant angular velocity w, Bardina et al. ref. [25] found that the effects of rotation should appear explicitly in the ε-equation. Furthermore, their simulation results show that the effect on ε is approximately linear in ω. As a result, they proposed a reduced ε-equation, which can be written as,

$\frac{\text{d}\epsilon }{\text{d}t}=-C_{\epsilon 2}\frac{{\epsilon }^2}{k}-C_{\epsilon 3}\omega \epsilon$, (12)

to analyze rotating isotropic turbulence. The last term in this equation was proposed by Bardina et al. ref. [25] to account for constant rotation effects in isotropic turbulence and the calculated results from Equation (12) are in good agreement with data. However, for more general flows, where the angular velocity Ω is a function of position and/or where mean strain might be present, a more appropriate angular velocity or rotation rate such as $\Omega ={\left[{\Omega }_{ij}{\Omega }_{ij}/2\right]}^{1/2}$ should be used instead. Here, the mean rotation rate tensor is given by ${\Omega }_{ij}=\left[\partial U_i/\partial x_j-\partial U_j/\partial x_i\right]/2$. In the present analysis, it is obvious that Ω should be adopted rather than ω.

3. Modified ε-Equation for Swirling and Rotating Pipe Flow

It can be inferred from the proposal given in Equation (12) that the effect of swirl could not appear explicitly as a coefficient in either the production or destruction term in the ε-equation. For isotropic turbulence, the proposed effect is linearly proportional to ω. The rotation effect may not be linear for non-isotropic turbulence. However, as a first attempt, a linear assumption is invoked. Since swirl affects the vorticity field and could contribute to vortex stretching, this suggests that the model term to account for swirl effects explicitly should take on a form similar to that introduced by Bernard and Speziale ref. [28], and Speziale and Bernard ref. [29] to model the effects of vortex stretching in the ε-equation. The model vortex stretching term was introduced to allow the asymptotic limit of homogeneous shear flow turbulence to be approached correctly. According to Bernard and Speziale ref. [28], the term should be proportional to $\left[{\left(R{e}_t\right)}^{1/2}\left({\epsilon }^2/k\right)\right]$. The present analysis adopts this proposal as the base and modifies it to account for swirl effects. Therefore, combining this proposal with the idea of Bardina et al. ref. [25] and the condition that the term should be dimensionally consistent with other terms in the ε-equation, a model similar to that proposed by Bardina et al. ref. [25] is obtained and can be written as $\Omega \left(k/\epsilon \right)R{e}_t^{1/2}\left({\epsilon }^2/k\right)$. With these suggested improvements, the modified ε-equation for swirling flow can now be written as

$\frac{D\epsilon }{Dt}=\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{\epsilon 1}{f}_1\frac{\epsilon }{k}{\Pi }_k-C_{\epsilon 2}{f}_2\frac{\epsilon \tilde{\epsilon }}{k}+C_{\epsilon 4}\Omega R{e}_t^{1/2}$ (13)

where C_ε4 is a coefficient associated with the swirling flow term. The exact form of C_ε4 will be determined later. It should be noted that the additional swirling flow model has a leading term of order $x_2^2$ near the wall. Therefore, its introduction does not affect the near-wall balance of Equation (9), which is carried out to order x₂ only.

In this form, Equation (13) does not reduce to Equation (9) when swirl vanishes because Ω is not zero as long as mean strain is present in the flow. Since the C_ε4 term is introduced to account for swirl effects explicitly, it should vanish as swirl goes to zero. In view of this, C_ε4 should be made parametric in S so that for a given swirling flow C_ε4 is a constant. This suggests that C_ε4 should be parametric in So, the inlet swirl number, which is constant for a given flow. A simple form proposed for C_ε4 is $C_{\epsilon 4}=C_{S1}\left[1-\exp \left(-C_{S2}{S}_o\right)\right]$ where C_S1 and C_S2 are constants. The value of C_S1 can be estimated based on that proposed by Bardina et al. ref. [25] for Equation (12). They suggested a value of 0.15 for C_ε3 in Equation (12). The present proposal has $R{e}_t^{1/2}$ in the model term, therefore, this implies that C_ε4 should be at least one order of magnitude less than that adopted by Bardina et al. ref. [25]. In other words, C_ε1 can at most take on a value of 0.015 while $\left[1-\exp \left(-C_{S2}{S}_o\right)\right]$ should be of order 1. Note that the additional swirling flow model term does not depend on Ω alone; it also depends on $S_o$ and $R{e}_t^{1/2}$. The introduction of Ret allows the local turbulence behavior to be accounted for in a more general way, while the dependence on $S_o$ permits the initial swirl to be taken into account. Later comparison of model calculations with experimental measurements shows that C_S2 = 5.0 is most appropriate, because it allows the effects of very small initial swirl to be accounted for properly in swirling flow modelling. These small swirl effects are known to influence the flow substantially according to Yajnik and Subbaiah ref. [14]. Thus improved, the ε-equation is more general than Equation (4) and Equation (9) because it is valid for both low-Reynolds-number flows with and without the presence of a wall and could account for swirl effects explicitly. The two-equation model based on Equations (3), (6), (10), (11) and (13) is the modified SSGZ or MSSGZ model for swirling flows. As before, the model can also be used in conjunction with an anisotropic eddy viscosity model.

Finally, it should be noted that Equation (13) can again be cast in a form similar to Equation (9) by modifying the definition for $\tilde{\epsilon }$. The new $\tilde{\epsilon }$ is given by

$\tilde{\epsilon }=\epsilon -\frac{C_{\epsilon 3}}{C_{\epsilon 2}{f}_2}\left[\nu {\left(\frac{\partial \sqrt{k}}{\partial x_i}\right)}^2+\frac{C_{S1}\left\{1-\exp \left(-C_{S2}{S}_o\right)\right\}}{C_{\epsilon 3}}\Omega R{e}_t^{1/2}k\right]$ (14)

This expression reduces to $\tilde{\epsilon }=\epsilon -2\nu {\left[\partial \sqrt{k}/\partial x_2\right]}_w^2$ at the wall and thus allows the wall boundary condition ${\epsilon }_w=2\nu {\left(\partial \sqrt{k}/\partial x_2\right)}_w^2$ to be satisfied exactly with

$\tilde{\epsilon }=0$. Therefore, Equation (9) is also applicable to swirling flows provided $\tilde{\epsilon }$ is now defined by Equation (14) instead of Equation (10). Casting the equation in this form provides a reasonable interpretation for the additional model terms in the ε-equation. Essentially, Equation (14) shows that the effects of viscosity and of swirl could be accounted for by introducing a reduced ε for the destruction term of the ε-equation; an idea first broached by Jones and Launder ref. [45]. However, instead of replacing one of the ε in the destruction term by a reduced dissipation rate $\tilde{\epsilon }$, Jones and Launder [40] proposed to solve a transport equation for the reduced dissipation rate, i.e. $\tilde{\epsilon }=\epsilon -2\nu {\left[\partial \sqrt{k}/\partial x_2\right]}_w^2$. Consequently, their proposed turbulence model becomes more complicated because one more $\tilde{\epsilon }$ -equation has to be solved.

4. Numerical Simulation

Swirling flows are calculated using the elliptic TEACH code developed by Gosman and Ideriah ref. [46]. Although it was originally developed to incorporate a two-equation model, it provides the basic structure for the implementation of more advanced turbulence models. The code uses a finite volume approach where the transport equations are integrated around their respective control volumes. The SIMPLE method of Patankar ref. [47] is adopted to solve the nonlinear system of equations together with an alternating direction, line-by-line, tridiagonal matrix solution scheme. The convective and diffusive terms are discretized using the hybrid scheme which is a combination of the central difference and upwind schemes depending on the cell Peclet number. The steady form of the equations is solved by introducing an under-relaxation factor. The iterations are continued until the sum of the normalized residuals reduces to 10^-6 for the dependent variables. In all the calculations, the $x_i$ coordinate system adopted in this paper will take on the following identity $x_1=x$, $x_2=y$ and $x_3=z$ in the TEACH code. Special attention is given to the placement of grids near the wall. A minimum of 5 grids are placed within $y^{+}<5$ with the first grid at approximately $y^{+}=1$, and 15 - 25 grids are located in the near-wall region $5<y^{+}<65$, depending on the total number of grids used. The calculation domain varies depending on the flow case considered. It is chosen so that the parabolic condition can be assumed at the outlet of the domain. A grid independent study has been carried out using 51 × 56, 102 × 56 and 133 × 56 grids along the axial (x) and radial direction (r). The grids are closely spaced near the pipe entrance and gradually become coarser towards the outlet. Usually, the calculations are started at a location where reasonably reliable experimental data are available; e is either interpreted from the measured $-\overline{uv}$ and k distributions or from the mean velocity, U, and the equilibrium turbulence assumption. Grid independent results are obtained with a 102 × 56 grid.

5. Results and Discussion

The merit of an anisotropic eddy viscosity assumption versus an explicit modelling of swirl in the ε-equation is assessed by using SASY, ASASY and the modified SSGZ to simulate three different cases of swirling flows. In subsequent plots, all results labelled SSGZ are to be interpreted as represented by the modified SSGZ; the adjective “modified” is ommitted in order to avoid clutter in the model prediction labelling. The first case is the rotating pipe experiments of Kikuyama et al. ref. [8]. These experiments were carried out with three rotation numbers, R_o = 0, 0.5 and 0.83 where R_o = W_o/U_o, W_o = wD/2 is the rotating pipe surface velocity and U_o is the maximum mean velocity. This gives rise to S_o = R_o/2 because in the experiments of Kikuyama et al. ref. [8] a uniform U was specified at the inlet. The present calculations are compared with the case where R_o = 0.5 or S_o = 0.25. Two other cases with a fairly high S_o are also considered. One case is the experiment of Kitoh ref. [21] with S_o = 0.97. The swirl in this case is generated by a swirler located at the entrance of the pipe. As a result of high swirl, a strong recirculation region was observed along the central core. An empirical decay law was proposed by Kitoh ref. [21] to describe swirl decay in this particular experiment. Another case is the experiment of Weske and Sturov ref. [9] with S_o = 1.3. The inlet swirl was generated by a rotating pipe; therefore, the rotating velocity field was represented by a solid-body rotation, and the inlet U was maintained uniform upstream of the rotating pipe. Even though the swirl number was larger than that specified in the experiment of Kitoh ref. [21], no central recirculation was observed in the experiment of Weske and Sturov ref. [9]. The three cases therefore bracket a swirl number range of 0.25 to 1.3 and a situation where a central recirculation is present. In the following, the case of Kitoh ref. [21] is used to assess the choice of C_S1 and C_S2 since it is the most complicated among those selected. Consequently, the two high swirl cases are presented first; this is then followed by a detailed comparison of the low swirl case.

In the experiment of Kitoh ref. [21], measurements at six axial locations are available. The mean velocity and Reynolds stress measurements at x/D = 0.57 are used as inputs to carry out the calculations, while ε is estimated from the measured k profile. Swirl decay is compared in Figure 3. Reference to the modified SSGZ in all figures is simply stated as SSGZ. In Figure 4 and Figure 5, only sample comparisons of the mean velocity distributions with measurements at x/D = 12.3 and 39.0 are shown because other stations show similar behavior. Different values for C_S2 have been attempted. This ranges from 1 to 5. Since the appropriate value is found to be 5; therefore, C_S2 = 5 is used from this point on. Four different values of C_S1 are used in the modified SSGZ model to perform the calculations; they are 0, 0.005, 0.01 and 0.015, respectively. Therefore, the U and W distributions at each location shown consist of results from four values of C_S1. On the other hand, only the C_S1 = 0.015 result is compared in the swirl decay plot (Figure 3). Besides the calculations from the modified SSGZ model, predictions from SASY and ASASY are also plotted in Figures 3-5 for comparison.

Figure 3. Comparison of calculated swirling decay with data of Kitoh ref. [21].

(a)(b)

Figure 4. (a) Comparison of calculated U with data of Kitoh ref. [21] at x/D = 12.3; (b) Comparison of calculated W with data of Kitoh ref. [21] at x/D = 12.3.

(a)(b)

Figure 5. (a) Comparison of calculated U with data of Kitoh ref. [21] at x/D = 39.0; (b) Comparison of calculated W with data of Kitoh ref. [21] at x/D = 39.0.

Kitoh ref. [21] suggested an empirical decay law, $S=S_0\exp \left(-\beta x/D\right)$ with β = 0.02362, to correlate the measurements. The calculated decay by SASY and ASASY tracks this behavior well in the near field (Figure 3) but deviates from it in the far field (x/D > 40). It should be pointed out that the swirl decay obtained from SASY and ASASY is essentially identical. The reason is that neither anisotropic nor isotropic eddy viscosity assumption has much effect on the prediction of the mean flow (Figure 4 and Figure 5). The ASAYS model gives additional information on the normal stresses, while the modified SSGZ model under-predicts swirl decay. This could be the result of a fairly correct prediction of U (Figure 4(a) and Figure 5(a)) but a less accurate calculation of W (Figure 4(b) and Figure 5(b)). As for SAYS and ASAYS, their predications of U and W are both incorrect. Consequently, the under prediction of U and W could give rise to a better prediction of swirl decay because S is related to U and W by $S=2{\displaystyle {\int }_0^{D/2}\left(UW{r}^2\text{d}r\right)/\left(D^3{U}_m^2/8\right)}$.

Both SAYS and ASAYS fail to predict a velocity depression in the central core at x/D = 12.3 and 39.0 (Figure 4(a) and Figure 5(a)). On the other hand, the modified SSGZ model with C_S1 = 0 is also incapable of producing a significant velocity depression (Figure 4(a)). It can be seen that explicit swirl correction in the ε-equation have a substantial effect on the prediction of the mean velocities; the effect is more pronounced in U than in W. With C_S1 = 0.015, a very significant axial velocity depression in the central core is predicted at x/D = 12.3. The model even manages to give rise to a substantial region of central recirculation (Figure 4(a)). Further downstream, even though measurements do not show a recovery from flow recirculation at x/D = 39.0, the prediction at the same location reveals that the velocity depression in the central core has recovered substantially (Figure 5(a)). This discrepancy could be partially attributed to measurement errors. It seems that a different value of C_S2 would be more appropriate because, with C_S2 = 5.0, significant discrepancy between predictions and measurements of W still exist at x/D = 12.3 (Figure 4(b)). However, calculations of the decay behavior (Figure 3) and the other two cases reveal that C_S2 = 5.0 is an appropriate choice for swirling flows over the range of S considered. In general, increasing C_S1 improves the agreement with U, but only slightly with W (Figure 4(b) and Figure 5(b)). It is believed that the W discrepancy noted in the central core is related to the W discrepancy shown in the near-wall region (Figure 4(b) and Figure 5(b)). If this discrepancy is to be improved, perhaps, swirl effects should also be accounted for explicitly in the near-wall region.

From this point on, only C_S1 = 0.015 and C_S2 = 5.0 are used in the modified SSGZ model to carry out the calculations of the other two cases. Comparisons for the case of Weske and Sturov ref. [9] with S_o = 1.3 are shown in Figures 6-8. Again, these calculations are carried out with inlet conditions given by the mean flow and normal stress measurements obtained at x/D = 0.35 while ε is estimated from the measured k and the shear stresses are deduced from the measured U

Figure 6. Comparison of calculated swirling decay with data of Weske & Sturov ref. [9].

(a)(b)

Figure 7. (a) Comparison of calculated U with data of Weske & Sturov ref. [9] at x/D = 5.1; (b) Comparison of calculated W with data of Weske & Sturov ref. [9] at x/D = 5.1.

and W assuming gradient transport. The decay behavior is compared in Figure 6 while some sample distributions of U and W at two axial locations are shown in Figure 7 and Figure 8, respectively. It can be seen that the calculated decay by the SAYS and ASAYS models is again identical (Figure 6). However, their decay curve lies below that given by the modified SSGZ model and agrees better with data, much like that shown in Figure 3 for the case studied by Kitoh ref. [21]. These results tend to substantiate the finding of Parchen et al. [22]; namely that conventional k-e models yield a slightly more accurate prediction of swirl decay. On the other hand, the modified SSGZ model calculations of U and W at x/D = 5.1 (Figure 7) and x/D = 20.0 (Figure 8) are in better agreement with data than those given by SASY and ASAYS. It appears that both the shape and magnitude of the velocity profiles are predicted fairly correctly at the two axial locations. The axial velocity depression in the central core is over-predicted by the

(a)(b)

Figure 8. (a) Comparison of calculated U with data of Weske & Sturov ref. [9] at x/D = 20; (b) Comparison of calculated W with data of Weske & Sturov ref. [9] at x/D = 20.

modified SSGZ model with C_S1 = 0.015, while both SAYS and ASAYS under-predict the velocity depression in the central core. This tends to suggest a smaller C_S1 value would be more appropriate for swirling flows. However, the third case, to be presented below, seems to lend credence to C_S1 = 0.015. The turbulence measurements were not reliable in the two high swirl number cases, therefore, the calculated k and $-\overline{uv}$ are not compared with the calculations. This, however, will be carried out for the low swirl number case where the turbulence measurements are much more accurate.

The final case to compare is the rotating pipe flow experiment of Kikuyama et al. ref. [8] with S_o = 0.25. In this case, measurements were available at x/D = 28.5 only. At the inlet, U is known from measurement and W is assumed to be given by W = wr. Again, the turbulence statistics are estimated from the mean U and W distributions assuming gradient transport and, once the normal stresses are known, e is estimated from k. The flow is still developing at x/D = 28.5, therefore, it is dominated by both stabilizing and destabilizing effects. In view of this, the present comparison not only validates the model's ability to replicate swirl effects at small S but also its ability to correctly predict the development of the rotating boundary layer in the pipe entrance region. The calculated U, W, k and $-\overline{uv}$ are compared in Figure 9(a), Figure 9(b), Figure 10(a) and Figure 10(b), respectively. Again, the modified SSGZ model results with C_S1 = 0.015 are in much better agreement with data than those of SAYS and ASAYS. The good agreement is not limited to the mean flow (Figure 9(a) and Figure 9(b)) but extends to k (Figure 10(a)) and $-\overline{uv}$ (Figure 10(b)). The improvements afforded by ASASY are quite evident in this case. In general, ASAYS yields calculations that are closer to the data (Figure 9 and Figure 10). However, the improvements are not as substantial as those given by the modified SSGZ model. In

(a)(b)

Figure 9. (a) Comparison of calculated U with data of Kikuyama et al. ref. [8] at x/D = 28.5; (b) Comparison of calculated W with data of Kikuyama et al. ref. [8] at x/D = 28.5.

(a)(b)

Figure 10. (a) Comparison of calculated k with data of Kikuyama et al. ref. [8] at x/D = 28.5; (b) Comparison of the calculated shear stress with data of Kikuyama et al. ref. [8] at x/D = 28.5.

fact, the modified SSGZ is the only model among the three that predicts the rotating boundary layer development correctly, including the behavior of k and $-\overline{uv}$ in the central core. Both SAYS and ASAYS yield a finite k in the pipe center (Figure 10(a)), which is not consistent with data. Only the modified SSGZ predicts the vanishing k and $-\overline{uv}$ correctly as the pipe center is approached (Figure 10(a) and Figure 10(b)). It should be pointed out that the modified SSGZ model results are very similar to those given by Yoo et al. ref. [6], who calculated the flow using a near-wall Reynolds-stress model (Model A-1 in Yoo et al. ref. [6] ). Their results are also reproduced in Figure 9 and Figure 10 for comparison with the model calculations. It can be seen that their calculations also yield an incorrect prediction of the boundary layer thickness. While the modified SSGZ model predicts the boundary layer edge to occur at about (1 − 2r/D) ≈ 0.8, the Reynolds-stress model of Yoo et al. ref. [6] yields a value closer to 1 and is comparable to that given by SAYS and ASAYS. In fact, the modified SSGZ model gives a much better prediction of k and $-\overline{uv}$ compared to those given by the Reynolds-stress model. This comparison, therefore, lends further support to the choice of C_S1 = 0.015 for swirling flow with a wide range of S.

6. Conclusions

Swirling flow simulations have been investigated using two-equation turbulence models with and without modelling stress anisotropy, as well as with and without a modified ε-equation. Based on these simulations, the following conclusions can be drawn.

1) An improved ε-equation with an explicit term to simulate swirl effects is formulated and is used in conjunction with the conventional k-equation and an isotropic eddy viscosity to calculate swirling flows in pipes. This is designated as the modified SSGZ model.

2) Application of the modified SSGZ model to calculate swirling flows with large S yields results that are dramatically improved from those obtained using a regular ε-equation. For the first time, a two-equation model capable of giving mean flow results comparable to experimental data is available. For swirling flow with large S and a central recirculation, the modified SSGZ could also replicate qualitatively the depression in U along the central core of the swirling pipe flow.

3) Among the two-equation models tested, good predictions of the mean and turbulence field of the swirling flow case with S_o = 0.25 as reported in Kikuyama et al. ref. [8] are obtained using the modified SSGZ model with C_S1 = 0.015 and C_S2 = 5.0 specified. The modified SSGZ model predictions of this case are comparable to those given by a near-wall Reynolds-stress model reported in Yoo et al. ref. [6].

4) Near-wall two-equation models with isotropic eddy viscosity (SAYS) and anisotropic eddy viscosity (ASAYS) fail to replicate swirling flows in the S range of 0.25 - 1.3. This is especially true if the flow also involves a central recirculation region. The predictions of the mean field by SAYS and ASAYS are essentially identical, except in the rotating boundary layer flow case where ASAYS shows a slight improvement. However, ASAYS yields additional information for the normal stresses. This shows that modelling stress anisotropy in a two-equation closure fails to improve the prediction of flows in the range of S = 0.25 - 1.3. Failure of two-equation models with isotropic or anisotropic eddy viscosity can be traced to the lack of an explicit term to model swirl effects in the ε-equation.

5) In view of these findings, a modified k-ε model is now available for complex flow simulation such as backstep flow, wall jet, and boundary layer flow with spin or swirl superposed on the individual elemental or combined flow.

Acknowledgements

Library, computer and other administrative support granted to RMCS by the Department of Mechanical Engineering, The Hong Kong Polytechnic University, Hong Kong, during the course of this work is gratefully acknowledged. Their help enable him to get hold of some key references cited in this paper.

Nomenclature

$C_{\epsilon 1},\cdots ,C_{\epsilon 4}$ model constants

D pipe diameter

g_i gravitational vector

g gravitational constant

$f_{w1},\cdots ,f_{\mu }$ damping functions introduced to in near-wall region

k turbulent kinetic energy

P mean static pressure

p fluctuating static pressure

r radial distance from origin of co-ordinate system

Re = U_mD/ν Reynolds number based on pipe diameter

$R{e}_t=k^2/\nu \epsilon$ turbulence Reynolds number

$R{e}_{\epsilon }={\left(\nu \epsilon \right)}^{1/4}{x}_2/\nu$, a Reynolds number based on the Kolmogorov velocity scale

R_o = W_o/U_o

S = ωD/2U_m is the ratio of angular momentum to the moment of linear momentum

S_o inlet swirl number

U_j mean velocity vector

U_m mean bulk velocity of pipe flow

U mean velocity along x-axis

U_o maximum mean velocity

V mean velocity along r-axis or y-axis

W circumferential mean velocity along the f direction

W_o = wD/2 is the rotating pipe surface velocity

u, v, w fluctuating velocities along x, y, z or z, r, f directions, respectively

u_τ friction velocity

x, y, z co-ordinates along the flow, normal to the flow, and normal to the x-y plane, respectively

Greek Symbols

ϵ_ijk alternating tensor

ε local turbulent dissipation rate

θ fluctuating temperature

ρ fluid density

ν fluid kinematic viscosity

ν_t eddy viscosity or turbulent kinematic viscosity

σ_k, σ_ε model constants introduced in Equations (2) and (3)

Τ turbulent shear stress

ω pipe rotational speed

Ω rotational speed of co-ordinate system about z-axis

Ω_ij mean rotation rate tensor

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Ahmed, S.A., So, R.M.C. and Mongia, H.C. (1985) Density Effects on Jet Characteristics in Confined Swirling Flow. Experiments in Fluids, 3, 231-238. https://doi.org/10.1007/BF00265106
[2]	So, R.M.C., Lai, Y.G., Hwang, B.C. and Yoo, G.J. (1988) Low-Reynolds-Number Modelling of Flows over a Backward-Facing Step. Zeitschrift für angewandte Mathematik und Physik ZAMP, 39, 13-27. https://doi.org/10.1007/BF00945719
[3]	Gerodimos, G. and So, R.M.C. (1997) Near-Wall Modeling of Plane Turbulent Wall Jets. Journal of Fluids Engineering, 119, 304-313. https://doi.org/10.1115/1.2819135
[4]	Sommer, T.P., So, R.M.C. and Lai, Y.G. (1992) A Near-Wall Two-Equation Model for Turbulent Heat Fluxes. International Journal of Heat and Mass Transfer, 35, 3375-3387. https://doi.org/10.1016/0017-9310(92)90224-G
[5]	Yoo, G.J. and So, R.M.C. (1989) Variable Density Effects on Axisymmetric Sudden-Expansion Flows. International Journal of Heat and Mass Transfer, 32, 105-120. https://doi.org/10.1016/0017-9310(89)90095-1
[6]	Yoo, G.J., So, R.M.C. and Hwang, B.C. (1991) Calculation of Developing Turbulent Flows in a Rotating Pipe. Journal of Turbomachinery, 113, 34-41. https://doi.org/10.1115/1.2927735
[7]	Murakami, M. and Kikuyama, K. (1980) Turbulent Flow in Axially Rotating Pipes. Journal of Fluids Engineering, 102, 97-103. https://doi.org/10.1115/1.3240633
[8]	Kikuyama, K., Murakami, M. and Nishibori, K. (1983) Development of Three-Dimensional Turbulent Boundary Layer in an Axially Rotating Pipe. Journal of Fluids Engineering, 105, 154-160. https://doi.org/10.1115/1.3240955
[9]	Weske, D.R. and Sturov, G.Y. (1974) Experimental Study of Turbulent Swirled Flows in a Cylindrical Tube. Fluid Mechanics-Soviet Research, 3, 77-82.
[10]	Anwer, M. and So, R.M.C. (1989) Rotation Effects on a Fully Developed Turbulent Pipe Flow. Experiments in Fluids, 8, 33-40. https://doi.org/10.1007/BF00203062
[11]	Hirai, S., Takagi, T., Tanaka, K. and Kida, K. (1987) Effect of Swirl on the Turbulent Transport of Momentum in a Concentric Annular with a Rotating Cylinder. Transaction of the Japan Society of Mechanical Engineers, B53, 432-443.
[12]	Kreith, F. and Sonju, O.K. (1965) The Decay of a Turbulent Swirl in a Pipe. Journal of Fluid Mechanics, 22, 257-271. https://doi.org/10.1017/S0022112065000733
[13]	Backshall, R.G. and Landis, F. (1969) The Boundary-Layer Velocity Distribution in Turbulent Swirling Pipe Flow. Journal of Fluids Engineering, 91, 728-733. https://doi.org/10.1115/1.3571214
[14]	Yajnik, K.S. and Subbaiah, M.V. (1973) Experiments on Swirling Turbulent Flows. Part 1: Similarity in Swirling Flows. Journal of Fluid Mechanics, 60, 665-687. https://doi.org/10.1017/S0022112073000406
[15]	Murakami, M., Kito, O., Katayama, Y. and Iida, Y. (1976) An Experimental Study of Swirling Flow in Pipes. Bulletin of the Japan Society of Mechanical Engineers, 19, 118-126.
[16]	Padmanabhan, M. and Janek, C.R. (1980) Swirling Flow and Its Effect on Wall Pressure Drop within Pipes. In: Swift, W.L., Barna, P.S. and Dalton, C., Eds., Vortex Flows, ASME Special Publications, New York, 65-75.
[17]	Ito, S., Ogawa, K. and Kuroda, C. (1980) Turbulent Swirling Flow in a Circular Pipe. Journal of Chemical Engineering of Japan, 13, 6-10.
[18]	Kito, O. (1984) Axisymmetric Character of Turbulent Swirling Flow in a Straight Circular Pipe. Bulletin of the Japan Society of Mechanical Engineers, 27, 683-690.
[19]	Kito, O. and Kato, T. (1984) Near Wall Velocity Distribution of Turbulent Swirling Flow in Circular Pipe. Bulletin of the Japan Society of Mechanical Engineers, 27, 1659-1666.
[20]	Algifri, A.H., Bhardwaj, R.K. and Rao, Y.V.N. (1987) Prediction of the Decay Process in Turbulent Swirl Flow. Proceedings of the Institution of Mechanical Engineers, 201, 279-283. https://doi.org/10.1243/PIME_PROC_1987_201_121_02
[21]	Kitoh, O. (1991) Experimental Study of Turbulent Swirling Flow in a Straight Pipe. Journal of Fluid Mechanics, 225, 445-479. https://doi.org/10.1017/S0022112091002124
[22]	Parchen, R., Steenbergen, W. and Voskamp, J. (1993) A Study of Swirling Flows in Long Straight Pipes. In: Rodi, W. and Martelli, F., Eds., Engineering Turbulence Modeling and Experiments, Elsevier Science Publishers, Amsterdam, 207-216. https://doi.org/10.1016/B978-0-444-89802-9.50024-1
[23]	Kobayashi, T. and Yoda, M. (1987) Modified k-ε Model for Turbulent Swirling Flow in a Straight Pipe. Japan Society of Mechanical Engineers International Journal, 30, 66-76.
[24]	Khodadadi, J.M. and Vlachos, N.S. (1990) Effects of Turbulence Model Constants on Computation of Confined Swirling Flows. AIAA Journal, 28, 750-752. https://doi.org/10.2514/3.10455
[25]	Bardina, J., Ferziger, J.H. and Rogallo, R.S. (1985) Effect of Rotation on Isotropic Turbulence: Computation and Modeling. Journal of Fluid Mechanics, 154, 321-336. https://doi.org/10.1017/S0022112085001550
[26]	Fu, S., Huang, P.G., Launder, B.E. and Leschziner, M.A. (1988) A Comparison of Algebraic and Differential Second-Moment Closures for Axisymmetric Turbulent Shear Flows with and without Swirl. Journal of Fluids Engineering, 110, 216-221. https://doi.org/10.1115/1.3243537
[27]	Kim, N. and Rhode, D.L. (2003) Refined Turbulence Modeling for Swirl Velocity in Turbomachinery Seals. International Journal of Rotating Machinery, 9, 451-459. https://doi.org/10.1155/S1023621X03000447
[28]	Bernard, P.S. and Speziale, C.G. (1992) Bounded Energy States in Homogeneous Turbulent Shear Flow—An Alternative View. Journal of Fluids Engineering, 114, 29-39. https://doi.org/10.1115/1.2909995
[29]	Speziale, C.G. and Bernard, P.S. (1992) The Energy Decay in Self-Preserving Isotropic Turbulence Revisited. Journal of Fluid Mechanics, 241, 645-667. https://doi.org/10.1017/S0022112092002180
[30]	Escue, A. and Cui, J. (2010) Comparison of Turbulence Models in Simulating Swirling Pipe Flow. Applied Mathematical Modeling, 34 2840-2849. https://doi.org/10.1016/j.apm.2009.12.018
[31]	Sagol, E., Reggio, M. and Ilina, A. (2012) Assessment of Two-Equation Turbulence Models and Validation of the Performance Characteristics of an Experimental Wind Turbine by CFD. International Scholarly Research Notices, 2012, Article ID: 428671. https://doi.org/10.5402/2012/428671
[32]	Javadi, A. (2016) Turbulence-Resolving Simulations of Swirling Flows. Ph.D. Thesis in Thermo and Fluid Dynamics, Department of Applied Mechanics, Division of Fluid Dynamics, Chalmers University of Technicology, Gothenburg.
[33]	Li, W.H., Ren, J., Jiang, H.D., Luan, Y.G. and Ligrani, P. (2016) Assessment of Six Turbulence Models for Modeling and Predicting Narrow Passage Flows, Part 2: Pin Fin Arrays. Numerical Heat Transfer, Part A: Applications, 69, 445-463.
[34]	Alhmadi, Y.H., Awadh, S.A. and Nowakowski, A.F. (2021) Simulation of Swirling Flow with a Vortex Breakdown Using Modified Shear Stress Transport Model, I&EC Research. Industrial & Engineering Chemistry Research, 60, 6016-6026.
[35]	So, R.M.C., Sarkar, A., Gerodimos, G. and Zhang, J. (1997) A Dissipation Rate Equation for Low-Reynolds-Number and Near-Wall Turbulence. Theoretical and Computational Fluid Dynamics, 9, 47-63. https://doi.org/10.1007/s001620050031
[36]	So, R.M.C., Aksoy, H., Sommer, T.P. and Yuan, S.P. (1994) Development of a Near-Wall Reynolds-Stress Closure Based on the SSG Model for the Pressure Strain. NASA Contractor Report-4618.
[37]	Speziale, C.G., Sarkar, S. and Gatski, T.B. (1991) Modeling the Pressure-Strain Correlation of Turbulence: An Invariant Dynamical Systems Approach. Journal of Fluid Mechanics, 227, 245-272. https://doi.org/10.1017/S0022112091000101
[38]	Speziale, C.G. (1987) On Nonlinear K-l and K-ε Models of Turbulence. Journal of Fluid Mechanics, 178, 459-475. https://doi.org/10.1017/S0022112087001319
[39]	Gatski, T.B. and Speziale, C.G. (1993) On Explicit Algebraic Stress Models for Complex Turbulent Flows. Journal of Fluid Mechanics, 254, 59-78. https://doi.org/10.1017/S0022112093002034
[40]	Launder, B.E., Morse, A., Rodi, W. and Spalding, D.B. (1972) Prediction of Free Shear Flows—A Comparison of the Performance of Six Turbulence Models. NASA SP-321, 361-426.
[41]	Launder, B.C., Reece, G.J. and Rodi, W. (1975) Progress in the Development of a Reynolds-Stress Turbulence Closure. Journal of Fluid Mechanics, 68, 537-566. https://doi.org/10.1017/S0022112075001814
[42]	Hanjalic, K. and Launder, B.E. (1976) Contribution towards a Reynolds-Stress Closure for Low-Reynolds-Number Turbulence. Journal of Fluid Mechanics, 74, 593-610. https://doi.org/10.1017/S0022112076001961
[43]	So, R.M.C., Lai, Y.G, Zhang, H.S. and Hwang, B.C. (1991) Second-Order Near-Wall Turbulence Closures: A Review. AIAA Journal, 29, 1819-1835. https://doi.org/10.2514/3.10807
[44]	So, R.M.C., Aksoy, H., Yuan, S.P. and Sommer, T.P. (1996) Modeling Reynolds-Number Effects in Wall-Bounded Turbulent Flows. Journal of Fluids Engineering, 118, 260-267. https://doi.org/10.1115/1.2817372
[45]	Jones, W.P. and Launder, B.E. (1972) The Prediction of Laminarization with Two-Equation Model of Turbulence. International Journal of Heat and Mass Transfer, 15, 301-314. https://doi.org/10.1016/0017-9310(72)90076-2
[46]	Gosman, A.D. and Ideriah, F.J.K. (1976) TEACH-T: A General Computer Program for Two-Dimensional, Turbulent, Recirculating Flows. Department of Mechanical Engineering, Imperial College, London.
[47]	Patankar, S.V. (1980) Numerical Heat Transfer and Fluid Flow. McGraw Hill, New York, NY.