Improved Dissipation Rate Equation for Swirling and Rotating Pipe Flows

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.

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 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 twoequation model needs to be demonstrated and established.Successful modelling of backstep flow, plane wall jets, and flow with density stratification using twoequation 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.

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 statio-Journal of Applied Mathematics and Physics nary 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 "stabiliz- 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.
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  1, a recirculation region did not exist near the pipe core be- cause 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.

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 simpli-fied, 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.
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 calcu-R.M. C. So Journal of Applied Mathematics and Physics lated 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.

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 complexi- ties 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 simula-tions 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.

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.

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 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, ( ) is the production of k, P is the mean pressure, ρ is fluid density, t is time, ( ) 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 where ( ) is a damping function used to ensure the disappearance of the effects of ξ away from the wall, 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 , where x 2 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 i j u u − and the eddy viscosity, ν t , are given by 2 2 3 ) 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 where 2 2 x u x τ ν + = and u τ 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 ( ) Here, x 2 is taken to measure positive away from the wall.Its relation to r, the radial coordinate, is x 2 = 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.

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 where 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].

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 con- 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 .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 where [ ] and, as before, f is given by ( ) while the constants are determined by So et al. ref. [43] to be  3), ( 6), ( 9) and (10).However, the damping function f μ is given by So et al.
where ( ) 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, 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 Here, the mean rotation rate tensor is given by In the present analysis, it is obvious that Ω should be adopted rather than ω.

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 ( ) ( ) . 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 .With these suggested improvements, the modified ε-equation for swirling flow can now be written as 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 2 2 x near the wall.Therefore, its intro- duction does not affect the near-wall balance of Equation ( 9), which is carried out to order x 2 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 ( )  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 Re in the model term, therefore, this im- plies 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 ( )  should be of order 1.Note that the additional swirling flow model term does not depend on Ω alone; it also depends on o S and The introduction of Ret allows the local turbulence behavior to be accounted for in a more general way, while the dependence on o S 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 ε .The new ε is given by .Therefore, Equation ( 9) is also applicable to swirling flows provided ε 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 ε , Jones and Launder [40] proposed to solve a transport equa- tion for the reduced dissipation rate, i.e.
. Consequently, their proposed turbulence model becomes more complicated because one more ε -equation has to be solved.

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 i x coordinate system adopted in this paper will take on the following identity 1 x x = , 2 x y = and 3 x 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 5 y + < with the first grid at ap- proximately 1 y + = , and 15 -25 grids are located in the near-wall region 5 65 y + < < , depending on the total number of grids used.The calculation do- main 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 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.R. M. C. So

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 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 Journal of Applied Mathematics and Physics 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 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 Journal of Applied Mathematics and Physics 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 uv − are compared in Figure 9 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 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.

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].

Figure 1 .
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.
ing" 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.

Figure 2 .
Figure 2. Co-ordinate system for swirling and rotating pipe flow.
) are changed to the following set: σ k = 1, σ ε = 1.4,C ε1 = 1.50 and C ε2 = 1.9.Further, the modelled destruction term k εε should be replaced by 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 The ε boundary condition would lead to 0 ε =  at the wall and thus rendering (9) valid as the wall is approached.Therefore, the SSGZ model consists of solving Equations ( 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.

)
Journal of Applied Mathematics and PhysicsThis expression reduces to and Figure5, 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 (Figure3).Besides the calculations from the modified SSGZ model, predictions from SASY and ASASY are also plotted in Figures3-5for comparison.

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.
R. M. C. So sidered.The model thus derived is the SSGZ model reported in So et al. ref. [36].
byBardina et al. ref. [25]is obtained and can be written as R. M. C. So DOI: 10.4236/jamp.2022.103048673 Journal of Applied Mathematics and Physics that proposed