A numerical and experimental study of coaxial jets

An algebraic stress model and the standard k-e{lunate} model is applied to predict the mean and turbulence quantities for axisymmetric, nonswirling coaxial jets without confinement. To investigate discretization schemes, namely, hybrid, power-law, and the flux-spline, are employed. In addition, an experimental study is conducted to provide data of good quality, especially near the inlet, for model assessment. The results show that the use of the algebraic stress model leads to better agreement between the numerical results and experimental data. © 1989.


Introduction
In recent years, considerable research has been directed towards the evaluation of various turbulence models for complex flows (e.g., Refs.1-3).However, a definitive evaluation has been hampered by the presence of excessive numerical (false) diffusion in the computed solutions and the lack of benchmark quality experimental data.Many of the prior computations have been performed using the hybrid (central/upwind) differencing scheme for the convective terms in the transport equations.Such a practice leads to excessive numerical diffusion which may be comparable to the physical diffusion.Further, the hybrid scheme responds very slowly to grid refinement, and an extremely large number of grid points may be required to obtain a grid-independent solution.A solution to the false diffusion problem is the use of higher order discrctization schemes for the convective terms.These schemes have been shown to be more accurate than the hybrid scheme for the same number of grid points and hence have the potential of providing a grid-independent solution without requiring an excessively large number of grid points.Examples of these schemes are QUICK (Ref.4), the skew upwind differencing sceheme (Ref.5), and the flux-splint scheme (Refs.6, 7).
In addition to numerical accuracy, another important consideration in the assessment of a turbulence model is the availability of reliable experimental data.The lack of "correct" boundary conditions may result in predictions which would not compare favorably with the experimental data and may cause erroneous inferences to be drawn about the turbulence model.The errors associated with the discrctization scheme, experimental uncertainty, and the turbulence model occur simultaneously and cannot be separated.It is, therefore, imperative that the errors from the first two sources be minimized so that the discrepancies between the experimental and computed results can definitely be attributed to the failure of the turbulence model.
In this paper, computations for unconfined axisymmetric nonswirling coaxial jets are reported.To assess the effect of false diffusion, solutions have been obtained using various differencing schemes on a fine grid.The schemes used are hybrid, power-law differencing scheme (Ref.8), and the fluxspine scheme (Refs.6,7).Both the standard k-e model (Ref.9) and the algebraic stress model (Ref.10) were applied in this study.These computations have been compared with detailed experimental data obtained using a two-component phase-Doppler technique.Prior studies which deal with some of the aforementioned issues related to the evaluation of turbulence models include the work of Leschziner and Rodi (Refs. 11,12), Hackman, Raithby, and Strong (Ref.13), and Claus (Ref. 14).
In the next section experimental procedure is described.This is followed by the mathematical model, which includes the turbulence models, solution algorithm, discretization schemes, boundary conditions, and other computational details.In the last section, the numerical results are compared with the experimental data, and the turbulence models are evaluated.

Test facility
A testing facility was designed to characterize a wide variety of flows under isothermal conditions (Figure l(a)).For the present study, an unconfined flow configuration was selected, operating with an axial jet injector surrounded by a nonswifling annular jet, as shown in Figure l(b).In this configuration, the injector was directed vertically downward within a 457-mm 2 wire mesh screen.The entire test assembly is surrounded by a flexible plastic enclosure which serves two purposes.First, the enclosure helps damp out extraneous room drafts.Second, and more important, the enclosure allows uniform seeding of the entrained air, thereby permitting unbiased measurements in the jet outer region.The test section air flows into a sealed collection drum and then into a suction vent connected to an exhaust blower.A slide valve on the vent allows for a variable duct back pressure.The support cage is mounted on an optics table from below via a two-axis traverse system in order to allow

Velocity measurement
A two-color, two-component laser anemometer system was used to measure the velocity components, At each spatial point the laser simultaneously measured two orthogonal components of velocity.In order to get all three components, two scans were taken.One was used to measure U, V, ~, v ~, uv components and the other one was able to measure U, W, ~, w ~I, ~-ff components.Thus all three components were measured, with the U velocity and its fluctuation measured twice.

Error estimates
The measurement errors can be broken down into four categories: (1) errors associated with the instrumentation and hardware, (2) uncertainty due to finite number of samples taken at each point, (3) repeatability limitations, and (4) validity of axisymmetric-flow assumption.
Category 1.The instrument accuracy is associated with the fringe spacing in the sampling volume.The error in this value is + 0.1 microns with the optical setup used.This translates into a _+ 1% error in the measurement of mean velocity.The fringes were orthogonal to within 1 °.
Category 2. At each point, 35,000 samples were taken, giving a maximum error of 0.04m/s (to a 95% confidence) in the regions with the highest fluctuating velocities (1.3 m/s).
Category 3.During the course of one measurement sequence (i.e., a single hardware setup and alignment), the mean and fluctuating velocities repeat to within 2% or 3% at a given point.The shear stress values generally repeat to within 5%.
Category 4. The results of mean and rms axial velocities measured along two orthogonal profiles show that the differences are within 2%.

Data sets
The flow conditions used for the case considered in this paper are given in Table 1.The center jet (diameter, D = 24.1 mm) is surrounded by an annular jet with the inner and outer diameters of 29 mm and 36.7 mm, respectively.The effective area ratio and axial velocity ratio of the annular jet to the center jet are 0.87 and 1.8, respectively.The experiment is documented in Ref. 15 and are tabulated following the format outlined in Ref. 16.

Mean flow equations
In this section, the equations which govern the distribution of the mean quantities are summarized.These equations are derived from the conservation laws of mass and momentum using time averaging and are expressed in tensor notation for steady and constant density flow as: where the bar is used to denote time-averaged quantities.
As a consequence of the nonlinearity of Equation 2, the averaging process used introduces unknown correlations which are modeled through a "turbulence model."

Turbulence models
In this paper, two turbulence closure models are considered, namely the standard k-5 model and an algebraic stress model.In the k-5 model, the turbulent fluxes are related to the mean fields through the assumption of an isotropic eddy viscosity and a turbulent Prandtl number as:

\OX~ OX,/
The eddy viscosity (#,) is obtained from the turbulent kinetic energy (k) and its dissipation rate (5) using the relation:

(p ,uA -gy,,-o5
(5) oxj ) The constants used in this model have been taken from Ref. 9 and are given in Table 2.
The k-5 model has been used with success in the calculation of various free shear flows and recirculating flows with and without swirl (e.g., Ref. 17), However, in flows with significant streamline curvature, the isotropic eddy viscosity assumption may not be able to describe the turbulent diffusion effects adequately.
The second turbulence model considered in this study is an algebraic stress model (ASM).The algebraic stress model is a special case of the Reynolds stress transport equation which The order-of-magnitude argument is performed on an equation for (2a~fl), which is obtained by subtracting the product of ~u/3 and the transport equation for 2k from the transport equation for u~j.The resulting equation becomes where the differential operator .~ is used to denote the combined convective and diffusive transport operators.Terms in Equation 8 are evaluated in powers of a and terms of the order of a 2 and higher are neglected.The result of this simplification is: where Pu is the production of the Reynolds stresses, ~b u is the pressure redistribution term, and Pk is the production of the kinetic energy.In the algebraic stress closure, a model for the pressure-strain term (~bq) is required.Here, the model of Launder et al. (Ref.19), which includes both the symmetric and antisymmetric mean-strain effects on redistribution modeling, is selected.These quantities are obtained from the following equations: OXk Oxd where In the above equations, cl, ~t, fl, and ?are model constants.
According to Launder et al. (Ref. 19), ct, fl, and ?can be related to a single constant C 2. Therefore, only two model constants are introduced in Equation 12 rather than four as shown.The constants used in the present study are listed in Table 3.

Solution algorithm
The computations for coaxial jets can be made using a parabolic marching procedure if the radial pressure gradients are small.Such a situation occurs if velocities in the two streams are comparable or the inner stream is faster and if the swirl is weak.However, if the swirl is strong and/or the outer stream is significantly faster, the radial pressure gradients become significant and a region of reverse flow develops.The ultimate goal is to extend this study to swirling flow analysis.Therefore, a calculation procedure based on elliptic flows was used.
The discretization equations are obtained using a controlvolume approach (Ref.8).The details of the differencing schemes for the convection and diffusion terms are given in the next section.The coupling between the continuity and momentum equations is handled using the SIMPLER algorithm (Ref.8).The algebraic equations are solved using a line-by-line tridiagonal matrix algorithm (TDMA).

Discretization schemes
Numerical solutions have been obtained using three schemes for the convection and diffusion terms in the transport equations.These schemes are briefly described below.
Hybrid scheme.In this scheme (e.g., Ref. 8) both the convection and diffusion terms are approximated using the central differencing scheme if the mesh Peeler number (Pc) is less than two.Outside this range, the upwind scheme is used for the convective terms and physical diffusion is neglected.
Power-law scheme.This scheme (Ref.8) is based on a curve-fit to the exact solution of the one-dimensional convectiondiffusion equation without source terms.This scheme becomes identical to the hybrid scheme for Pc> 10.
Flux-spline scheme.The hybrid and power-law differencing schemes can be considered as approximations to the exponential scheme (Ref.20) which results from the exact solution to the one-dimensional convection-diffusion equation without a source.In the derivation of these schemes, the total flux (convection + diffusion) is assumed to be uniform between two grid points.These schemes work well only in problems in which either the flow is closely aligned with the grid lines or there are no strong cross-flow gradients.If such idealized conditions are not encountered, the locally one-dimensional assumption used in these schemes gives rise to numerical (false) diffusion.
In the flux-spline scheme (Refs.6, 7), the total flux is assumed to vary in a piecewise linear manner within a control volume.This assumption leads to a scheme in which the discretization coefficients are identical to those from the exponential scheme but there is an additional source term which involves the differences in fluxes at adjacent faces of a control volume.The presence of this source term enables the flux-spline scheme to respond to the presence of sources and/or multidimensionality of the flow.

Boundary conditions
A calculation procedure for elliptic flow requires boundary conditions on all boundaries of the computational domain.Four kinds of boundaries need consideration, namely, inlet, axis of symmetry, outlet, and the entrainment boundary.At the inlet boundary, which was located at the first measurement plane, the measured profiles of U and V were prescribed.The k-profile was obtained from the measured Reynolds stresses.These profiles are shown in Figure 2.This kinetic energy distribution and the measured shear stress profile were used to derive the e values at the inlet plane through the following relationship: At the axis of symmetry, the radial velocity and the radial gradients of other variables are set to zero.At the outlet, axial diffusuion is neglected for all variables.Along the entrainment boundary, which was placed sufficiently far from the axis of symmetry, the quantity (rV) was assumed constant.In addition, the axial velocity U was assumed zero and k and ~ were assigned arbitrarily low values yielding an eddy viscosity,/~t = 10/~.

Computational details
The computational mesh used for all calculations consisted of 76 x 69 nonuniformly distributed grid points in the axial and radial directions.A finer grid spacing was used near the inlet, centerline, and in the shear layer.The computational domain extended from the first measurement plane, located downstream of the nozzle exit at a distance of 2.0 mm to 40 inner jet diameters downstream of the nozzle exit.In the radial direction, the entrainment boundary was placed at a distance of six jet diameters from the axis of symmetry.
The convergence criterion used to terminate the iterations was that the absolute sums of the mass and momentum residuals at all internal grid points, normalized by inlet mass and momentum fluxes, be less than 5 x 10-3.

Results and discussion
In this section, the results for nonswirling coaxial jets are presented.The numerical results were obtained using two turbulence models with various differencing schemes.The calculated mean and turbulence quantities are also compared with the measurement at selected stations.

The k-8 turbulence mode/
The effect of different discretization schemes is shown by comparing the predicted axial velocity profiles at three streamwise locations, namely, x=15, 35, and 75mm.The velocity profiles are presented in Figure 3.It is noted that, except for some minor differences, all three schemes for the convective terms yield nearly identical results.In earlier studies (Refs.6, 7), it was shown that in the regions of high Peclet number the fiux-spline results are more accurate than those from the power-law scheme• The fact that for the present situation there   are no significant differences between results from these schemes indicates that the results are grid-independent.The differences between the hybrid and the power-law schemes are attributed to the different treatments of the diffusion terms.The computed results at the selected axial stations compare reasonably well with the experimental data.The computations consistently show sharper gradients than the experiment at the points of the maximum and minimum velocity.Figure 4 shows the kinetic energy profiles at three axial locations• The experimental kinetic energy profiles were derived from the measured Reynolds stresses.In these figures, results   from the power-law scheme and flux-spline scheme have been shown.Again, the two sets of computations are in close agreement with each other.Most of the differences are seen in the regions of steep gradients where the flux-spline results are expected to be more accurate.The agreement between the predicted and experimental values of kinetic energy is not as good as that for the axial velocity.Even though the trends are similar, the predicted kinetic energy levels are smaller than those derived from the measurements.
Since the present calculations are essentially free of numerical diffusion, the discrepancies between the experimental data and the predictions can be attributed to two sources--improper boundary conditions at the inlet plane and the deficiencies of the turbulence model.As regards the inlet conditions, all quantities except the dissipation rate were prescribed from the experiment.The e values, however, were derived from the measured shear stresses and the mean velocity gradients.The uncertainties in the measurements and in the evaluation of the velocity gradients may lead to errors in the e values which would adversely affect the calculations at downstream locations.
Numerical experiments indicate that the inlet e profile is the single most important factor in predicting the maximum values of mean and turbulence quantities, provided a reasonable inlet kinetic energy distribution is available.To study the sensitivity to inlet e profile, calculations were also made using an alternative distribution, which were derived from the turbulence kinetic energy and an assigned length scale distribution (3% of the radius).The inlet ~ profiles for both cases are shown in Figure 5.The major differences between these two conditions are near the centerline region; however, the peak values are about the same.The predicted results of mean axial velocity and turbulent kinetic energy at an axial location of 15 mm are shown in Figure 6.The results show that the turbulent kinetic energy has decreased due to excessive inlet dissipation rate in the inner region.On the other hand, the mean velocity is not affected significantly.This can be attributed to the fact that the maximum value of the inlet e in the annular region has not been changed considerably.

The algebraic stress model
In this section, the predictions using the ASM have been compared with those from the k-e model.Similar to the trends observed in the k-~ model calculations, the effect of various discretization schemes on the ASM results was found to be rather insignificant.Consequently, results from different schemes will be shown only for some cases.
The predicted mean axial velocity profiles from the ASM and k-e models have been compared in Figure 7.These results were obtained using the flux-spline discretization scheme.The use of ASM improves the overall agreement between the predictions and the experimental data.The major differences between the two turbulence models are seen in the regions where a maxima or a minima occurs in the velocity profiles.
The predicted turbulent shear stress from the ASM has been compared with the experimental data in Figure 8.Here results from the power-law differencing scheme have also been included to assess the numerical accuracy of the results.Both discretization schemes give nearly identical results, indicating that the solution is grid-independent.The positive peak in the shear stress profile corresponds to the shear layer between the two streams and the negative peak corresponds to the shear layer associated with the expansion.The agreement between the calculation and the experimental data is good, although the peak values are not well predicted.
The normal stresses at different axial locations are shown in Figure 9. Again, results from both the power-law and flux- for ~, underpredicted for ~, and closely predicted for w -'z.This clearly indicates the lack of performance of the pressure-strain model.One reason for this could be that the constant C2 used is not suitable for complex turbulent flows.Since C2 is determined from simple turbulent flows in local equilibrium, it would be more appropriate for equilibrium ASM than for nonequilibrium ASM.Another reason could be the incorrect modeling of the mean strain part of the pressure-strain term.It -

Concluding remarks
Based upon the preceding discussion, the main conclusions can be summarized as follows: (1) The mean flowfield predicted using the k-8 model agrees reasonably well with the experimental data.However, the o Lo. (3) The shear stress predictions using the algebraic stress model are in good agreement with the experimental data.For normal stresses, there are considerable discrepancies between the experimental and numerical results.(4) The discrepancies between the data and the algebraic stress model solution may be related to the pressure-strain correlation model.(5) The effectiveness of the turbulence models, to some extent, can be obscured by boundary conditions.Inlet conditions for e, especially the peak values, are found to be an important factor in determining the maximum velocity values.

Figure I
Figure I Experimental setupthe two degrees of freedom in the horizontal plane.The injector is mounted on a vertical spar to provide the third degree of freedom in the vertical direction.The control volume spatial location is monitored via a three-axis digital indicator which permits positioning to within 0.01 mm.Data were obtained at seven axial stations: 15, 25, 35, 50, 75, 150, and 300mm from the exit plane of the injector.At each axial station, between 10 and 20 radial points were scanned as determined by the desired level of profile resolution.

( 4 )
Two additional partial differential equations are solved to obtain k and 5.These are Figure 2

Figure 3
Figure 3 Comparison of measurements with mean axial velocity calculations

Figure 4
Figure 4 Comparison of measurements with turbulent kinetic energy calculations

uFigure 6 Figure 7
Figure 6 Profiles of mean axial velocity and turbulent kinetic energy using different inlet turbulent dissipation rate Coaxial jets: M. Nikjooy et al.

U V m2/s 2 Figure8
Figure8 Comparison of the algebraic stress model (ASM) prediction of turbulent shear stress with measurements

7 Figure 10 Figure I 1
Figure 10 Comparison of the algebraic stress model prediction of radial normal stress with measurements

Table 2
Values of constants in the k-e model Coaxial jets: 1t4.Nikjooy et al.relates the individual stresses to mean velocity gradient, turbulent kinetic energy, and its dissipation rate by way of algebraic expressions.Algebraic stress models can be classified into two categories.The first is based on a local equilibrium assumption for the turbulence field, whereby the turbulence transport terms are neglected compared to the local production and dissipation of turbulence.A second class of ASM is based on the local nonequilibrium assumption.Approaches of this kind, where the convection and diffusion transport of turbulent stresses are approximated, have been developed byMellor and Yamada  (Ref.10)andRodi(Ref.18).Following the recommendation in Ref. 3, the model proposed byMellor and Yamada (Ref.10)hasbeenadopted in this study.In the Mellor and Yamada model (Ref.10), the Reynolds transport equations are simplified through an order-ofmagnitude argument based on a2=O(a~), where % is the nondimensional measure of anisotropy and is given by