NUMERICAL INVESTIGATION OF TURBULENT SWIRLING FLOWS THROUGH AN ABRUPT EXPANSION TUBE

A numerical investigation of turbulent swirling flows through an abrupt expansion tube is reported. The TEFESS code, based on a staggered Finite Volume approach with the standard kε model and first-order numerical schemes built-in, was used to carry out all the computations. The code has been modified in the present work to incorporate the ASM and two second-order numerical schemes. The ASM, which includes the non-gradient convection terms arising from the transformation from Cartesian to cylindrical coordinates, was investigated for isothermal flows by applying it to the flow through an abrupt expansion tube with or without swirl flows. In addition, to investigate the effects of numerical diffusion on the predicted results, two second-order differencing schemes, namely, second-order upwind and the quadratic upstream interpolation, were used to compare with the first-order hybrid scheme. An abrupt expansion tube with non-swirling flow, predicted results using both the k-ε model and the ASM were in good agreement with measurements. For swirling flows, the calculated results suggested that the use of the ASM with a second-order numerical scheme leads to better agreement between the numerical results and experimental data, while the k-ε model is incapable of capturing the stabilizing effect of the swirl.


INTRODUCTION
Swirl flows have been of considerable interest over the past decades because of their occurrence in industrial applications, such as furnaces, utility boilers, internal combustion engines, gas turbine combustors and dust collectors [1,2].Swirl has been used in combustion systems to enhance the flame stability, the mixing and heat transfer besides prolonging the fuel residence time and abating the pollutants.This is because under appropriate conditions, swirl can be employed to induce a central recirculation zone.The recirculating flow generates additional turbulence in the shear layer between the reverse flow and the surrounding forward flow and helps stabilize the flame in combustors.Swirling turbulent flows are physically complex in nature due to the effect of a swirl-turbulence interaction.The turbulence structure in swirling turbulent flows is generally highly non-isotropic and non-homogeneous.
Computation of swirling flows is a difficult and challenging task.Large velocity gradients appear in these flows, so numerical problems and turbulence modelling play a significant role in their analysis.The commonly used, the k-ε model may not be suitable for simulating swirling turbulent flows [1].It is also found that the use of modified k-ε models or even the non-linear kε model [3] leads to no significant improvement of the predictions in swirling flows.The second-order moment closure models, i.e., the Reynolds stress model (RSM) and the algebraic Reynolds stress model provide better methods for the simulation of swirling turbulent flows [4,5,6], but sometime, the original ASM based on Rodi's approximation [7] cannot give satisfactory results for certain aspects of swirling flows [3].The RSM is regarded as a most logical approach to the turbulence closure problem, which does not need any ad hoc modification for extra strain rates.However, in the prediction of swirling flows with the RSM, it is necessary to solve a total of 11 governing differential equations of elliptic type: a continuity equation, three momentum equations, an ε-equation, and six equations for the Reynolds stresses.This leads to much extra computational effort to solve six Reynolds stress transport equations simultaneously [4,5,6] and much attention needs to be paid to numerical stability and inlet boundary conditions.It is for this reason that a simplified algebraic Reynolds stress turbulence model in axisymmetric cylindrical co-ordinates is employed for simulating strongly swirling flows.
This article deals with the numerical simulation of the flow through an abrupt expansion tube with and without swirl by utilising the present ASM and various numerical schemes.The mathematical model including the turbulence models, numerical solution and other computational details is described.Comparisons of the calculated gas tangential and axial velocities with the test data measured in a sudden enlargement duct [8,9] are made to evaluate the turbulence models and the numerical schemes used.

Governing Equations
The phenomenon under consideration is governed by the steady two-dimensional axisymmetric form of the continuity and the time-averaged incompressible Navier-Stokes equations.In the Cartesian tensor system these equations can be written in the following form: where ρ, u i , p and x i are the density, mean velocity tensor, mean pressure and coordinate tensor respectively.The mean viscous stress tensor, ij t , is approximated as: where μ is laminar viscosity.The time-averaged Reynolds stress tensor, τ ij = -' ' j i u u ρ , in the above equation is not known and thus, models are needed to express it in terms of the solution variables.In the present study, two turbulence models are used, namely the k-ε model and an algebraic stress model (ASM).The k-ε model has been reviewed in references [10,11] and it will be described briefly.The standard k-ε model version relates the turbulent eddy viscosity to the turbulence kinetic energy k and the dissipation rate ε through Boussinesq's approximation as: ( ) where is the turbulent eddy viscosity and ε is the dissipation rate of turbulence kinetic energy (TKE).The modelled equation of the TKE, k is given by: (5) in which μ e = μ t + μ is effective viscosity.Similarly the dissipation rate of TKE is given by the following equation: where G is the rate of generation of the TKE while ρε is its destruction rate.G is given by: The boundary values for the turbulent quantities near the wall are specified with the wall function method [10].C μ = 0.09, C ε1 = 1.44,C ε2 = 1.92, σ k = 1.0, and σ ε = 1.3 are empirical constants [10,11] in the turbulence transport equations.Reynolds-averaged transport equations can be solved for τ ij , [10,11] the modelled equations for which are: where in which C 1 = 2.5, and C 2 = 0.55 are model constants.

Algebraic Reynolds Stress Model (ASM)
For simplicity in solving the six Reynolds stresses, Rodi's approximation [7] is used in this study and the Reynolds stress transport can be expressed in algebraic form as follows: Substitution of eqs.( 5) and ( 8) into eq.( 9) gives the desired algebraic expression for τ ij : The ASM expressions can thus be rewritten as: (11) where the empirical constant λ ,was found to be 0.135, [12] is defined as:

Common Form for the Equations
All the governing equations can be re-organised and expressed in a standard form that includes the convection, diffusion, and source terms for 2-D axisymmetric flows as follows: where φ may stand for any variable including the velocity components, Γ φx and Γ φr are the exchange coefficients for φ, and S φ is the source term.

Solution Procedure
In the present computation the time-averaged Navier-Stokes equations, equations ( 1) and ( 2); the TKE equation, equation (5); the TKE dissipation rate equation, equation ( 6) are solved numerically by a control-volume finite-difference method [13,14] together with the turbulence model equations, equation ( 4) for the k-ε model or equation for the ASM.The SIMPLE algorithm is utilised for pressure-velocity de-coupling and iteration [13,14].The discretization of the governing equations is accomplished by means of the upwind, the hybrid, the secondorder upwind (SOU) and the quadratic upstream interpolation for convective kinematics (QUICK) schemes and the source term linearisation on a staggered grid cell.The underrelaxation iterative TDMA line-by-line sweeping technique is used for solving the resultant finite-difference equations.Due to the highly non-linear and coupling features of the governing equations for swirling flows, lower underrelaxation factors ranging from 0.005 to 0.1 are chosen for the three velocity components to ensure that the stability and convergence of the iterative calculation.Wall function was used for calculating wall shear stresses at the grid nodes along the walls.The exit boundary was chosen at the center tube outlet where zero gradient conditions were adopted for all variables except the axial velocity u, which is subject to continuity constraints.The computation was carried out using a PC computer.About 20,000 iterations were needed to achieve satisfactory convergence for each calculation case, which requires about 3 hr of computer time.

ABRUPT EXPANSION TUBE FLOW OF YOON (1982)
An experimental study of the turbulent swirling flows through an abrupt expansion tube [8,9] for which measurements are available was selected to validate the model.This flow arrangement shown in Fig. 1 had initial velocity and fluid properties as shown in Table 1.A single air stream entered the test section through a secondary annulus and passed through an adjustable vane swirler as shown in the figure.The inlet swirl vanes of the swirl generator with a central hub that had been blocked-off were set at an angle of 0 or 38 degrees.The secondary annulus expanded into a test chamber whose diameter was twice that of the secondary tube.An expansion block was inserted to position the exit plane of the swirl generator and its hub was 0.032 m upstream of the expansion plane.The hub and the swirl give rise to a significant reverse flow along the axis besides the recirculation zone near the plane of expansion.Experimental mean velocity profiles were provided at the expansion with which to initiate the code predictions.In the present computation, performances of both the standard k-ε model and the ASM with four different discretization schemes, namely, upwind, hybrid, SOU and QUICK, are examined by comparing the predicted velocity profiles at four locations, namely, x = 0.149, 0.298, 0.448, and 0.597 m for both non-swirling and swirling flows.

Fig. 2a: Effect of grid distributions on the velocity profiles for the k-ε model
First, the effect of grid sizes on the solution was investigated with two grid densities of 50 × 30 and 90 × 50 and is shown in Figs.2a and 2b for the k-ε model and the ASM, respectively.It is seen that both grids give nearly similar predictions for each turbulence model, indicating that a grid of 50 × 30 or finer led to results that were sufficiently grid-independent.The effects of numerical schemes on the predicted results were also examined with a 50 × 30 grid and are presented in Figs.3a and 3b for both the turbulence models.It is noted that, except for some minor differences, all schemes yield almost identical results.Further examination reveals that the QUICK and the SOU results are more accurate than those from the other schemes when compared with the experimental data.The predicted radial profiles of axial and radial velocities for the ASM and the k-ε model using the QUICK scheme are compared in Fig. 4. The ASM led to better overall agreement between predicted and measured axial velocity profiles although the slight under-prediction is visible in the near inlet region.Radial velocity predictions by both the turbulence models fail to match the experimental data.This flow has also been calculated by Sloan et al. [15] with only the k-ε model and their results closely match those obtained here.They also reported that the integrated experimental mass flow rates differed from the input mass flow rate by as much as 25 -30%.

Fig. 2b:
Effect of grid distributions on the velocity profiles for the ASM

Predictions with the standard k-ε turbulence model
The effect of grid sizes on the calculated results was examined and the results are compared with measurements at four downstream stations.Figure 5 shows the sensitivity of the QUICK scheme and different grid densities to the predicted axial and tangential velocity profiles.The axial velocity variation along the centre-line is shown in Fig. 6.Grids finer than 40 × 20 give almost the same results, especially at locations far from the inlet.Close to the inlet, however, finer the grid spacing led to a more accurate solution.The tangential velocity predictions with the k-ε model display a rapid decay to a forced vortex profile (solid-body rotation) while the experimental data correspond to a combined free and forced vortex profile as evident from Fig. 5.The results of Fig. 6 indicate that the predicted axial velocity recovers and progresses to uniformity at a faster rate than in the measurements.The effect of different discretization schemes on the predicted axial and tangential velocity distributions for the four schemes used is presented in Fig. 7 with a 60 × 30 grid.It is found that the QUICK and the SOU generally give better results than the other two, upwind or hybrid.respectively.Two recirculation zones are identified; one is at the corner and the other, an internal recirculation zone (IRZ), near the inlet.It is clear that the tangential velocity predictions do not agree with the measured data even for trends.Since the calculations were fairly free of numerical diffusion, discrepancies between the data and the predictions can be attributed to two sources: -improper boundary conditions at the inlet plane and deficiencies of the turbulence model.All inlet values except the ε and the k, were obtained indirectly from the experiment.The values of ε were derived from a constant length-scale assumption and those of k from an approximation of total kinetic energy.To study the sensitivity of the inlet ε and k profiles to the flow field, calculations were made using different ε and k distributions, by increasing or decreasing their estimated values by a factor of 10.It was found that the calculated velocity was not affected significantly.The poor agreement between the predictions and measurements is, therefore, more likely to be due to the deficiencies of the turbulence model, as observed by many other past investigator.

Fig. 8:
Streamlines predicted by the k-ε model

Fig. 9:
Vector plots of velocity predicted by the k-ε model

Predictions with the ASM
Predictions using the ASM are compared with those from the k-ε model and data from the measurement.The effects of various discretization schemes on the ASM results are as significant as in the case of the k-ε model calculations.Results with the QUICK or the SOU were closest to the measured data, and hence only results with the QUICK scheme will be discussed here.The radial distributions of axial and tangential velocities predicted by the ASM are compared with the measurements in Fig. 10.These results were obtained using the QUICK scheme, β = 0.4, and different grid size distributions.Grids finer than 40 × 20 appear to yield solutions which were sufficiently grid -independent.The use of ASM results in better overall agreement with the measurements than possible with the k-ε model.The major difference between calculations with the two turbulence models is most clear in the recirculation region and in the tangential velocity profiles.The agreement between the calculation with the ASM and the experimental data is fairly good.The peak values and a combined forced and free (Rankine) vortex motion for the tangential velocity profiles are well predicted.The recirculation zone (IRZ) is, however, longer and wider than that with the k-ε model.
Figure 11 shows the variation of centre-line axial velocity for the QUICK scheme and different grid densities.It is found that grids finer than 40 × 20 give nearly the same results.The negative and positive peaks in the profiles are well predicted.In comparison with the measurements, the ASM performs better than the k-ε model.The effect of β on the ASM results is shown in Fig. 12 that presents the radial profiles of axial and tangential velocities.It is found that a value of β between 0.3 and 0.4 appears to give the best agreement between prediction and measurement.
The ASM with added convection (β > 0) under-predict the axial velocity profiles in the core region if β is larger than 0.4 but predict well their maximum values.In addition, the peak values of the tangential velocity profiles are slightly under-predicted for all β values at downstream locations.Streamlines and velocity vectors predicted with the ASM are illustrated in Fig. 13 and Fig. 14 respectively.It is observed that the size of the recirculation zone calculated by the ASM is larger than that for the k-ε model.

Fig. 10:
Comparison of measurements with the ASM predictions using various grids The flows of Yoon [8] were also predicted by Sloan et al. [15] using a 38 × 35 grid with a modified hybrid/power-law differencing scheme and the same inlet conditions as the present study.Four different turbulence models, namely, the k-ε model, LPS Richardson number (a modified k-ε model for swirling flows), the algebraic stress model, and the algebraic stress model with added convection, were used.The results with the k-ε model were similar to the present calculations with the k-ε model and an upwind scheme.Their predictions by the algebraic stress model and the LPS Richardson number differ only marginally from those with the standard k-ε model while the algebraic stress model with added convection led to slightly better predictions but not as good as the present ones, although their model was based on the same assumptions.This can be attributed to (1) a coarse grid, (2) numerical diffusion and (3) improper boundary conditions at the inlet plane.Numerical experimentation has shown that, for sufficiently accurate results with the QUICK scheme, a grid of at least 20 × 25 was needed between the inlet and one diameter from it in order to cover the central toroidal recirculation zone.With the duct length of around 10 diameters, the 38 × 35 grid used in their predictions was probably not sufficiently fine near the inlet region.Also, Nikjooy and Mongia [16] have reported the use of first-order differencing schemes for the transport equations introduces significant numerical diffusion for a coarse grid.The measured inlet stress profiles were not provided and Reynolds stress estimation for the ASM at the inlet can be difficult.In their predictions, Sloan et al. reported that dominant stresses from a plausible approximation were used as input and the remainder was calculated theoretically from the model.Erroneous inlet conditions can have a substantial effect on the predicted flow field.Abujelala and Lilley [17] also calculated the flows of Yoon.Conditions used in their calculations are the same as the Sloan et al. cases except that they used the experimental inlet profile shapes measured by Sander and Lilley [18], at flow rates different from those of Yoon.Their predictions are similar to those of Sloan et al., and show the same rate of decay of swirl velocity profiles, and manifest the same lack of agreement.

CONCLUSION
In the present work, a numerical study has been conducted to investigate flow-field characteristics in an abrupt expansion tube.From the computation results, it can be concluded as follows: • The major difference between calculations with the two turbulence models is most clear in the recirculation region and in the tangential velocity profiles.The calculated results of gas tangential and axial velocities using both the k-ε model and the ASM are in generally good agreement with the measurements.

•
Influences of the numerical schemes for convection transport are found to be significant in simulating annular swirling turbulent flows.The computations show that the use of the QUICK and the SOU leads to better agreement between the numerical results and experimental data than the upwind and hybrid schemes.
• It is found that a value of β between 0.3 and 0.4 appears to give the best agreement between prediction and measurement.

Table 1 :
Data for flow through an abrupt expansion tube