ON THE RELIABILITY OF REYNOLDS-AVERAGED NAVIER STOKES PREDICTION OF MEAN FLOW CHARACTERISTICS OF A RECTANGULAR TURBULENT JET IN CROSSFLOW

A jet of fluid discharging into a cross stream, als o known as Jet In Crossflow (JICF), has received ma ny experimental and numerical investigations. In addit ion to the fundamental understanding of threedimensional mixing and shear flow characteristics, the fluid dynamics research community often regarded it as a benchmark test case for validating turbulence models. Although many authors considere d the canonical case of a jet issuing from a circular orifice, the rectangular shape has received less numerical investigations. The present study deals w ith a jet issuing from a rectangular duct into a confined crossflow domain in which five jet-to-cros sflow velocity ratios ranging from 3.3 to 10 are considered. The analysis focuses on the reliability of three two-equation turbulence models, namely k-ε, k-ω and SST in predicting this type of complex flow phenomena. Comparisons with previous large-eddy simulation results and available test data for the same problem have revealed good agreement in predicting ‘mean’ flow properties, but relative poo r agreement in predicting the second-order statisti cs. It indicates that this type of flow exhibits significa nt non-equilibrium behavior for which the commonly used two-equation turbulence models are unable to d eal with. Thus it is necessary to apply anisotropic turbulence model such as Reynolds stress model or h igh-fidelity large-eddy simulation method.


INTRODUCTION
Turbulent flow mixing of jet flow discharging into a crossflow (shown schematically in Fig. 1) arises in many situations of technologically important fields. For example, a non-reacting transverse jet is a configuration applicable for chimney stacks, Vertical and/or Short Take-Off and Landing (V/STOL) aircrafts, dilution of combustion gases in gas turbines and film cooling of turbine blades. A reacting transverse jet could affect the flame stabilization of a fuel jet issuing into a crossflow as a model of stack flares and also the secondary combustion zones in a gas turbine combustion chamber. On aspect of flow physics, Jet In Crossflow (JICF) has been regarded as a challenging test case for fundamental understanding of three-dimensional turbulent mixing and shear layer flows, thus it is often used for validating turbulent models by the fluid dynamics research community (Chochua et al., 2000;Acharya et al., 2001;Chassaing et al., 1974;Schluter and Schonfeld, 2000;Muppidi and Mahesh, 2007).
So far, various researches have been conducted for circular jet flows issuing into unconfined crossflow domain, because of its application in aerospace engineering such as vertical and/or short take-off and landing aircrafts, steering of rockets, film cooling of turbine blades, fuel injection into combustors, etc. More recently, ecological aspects such as plumes of smoke stacks, volcanoes, or (tunnel) fires have also gained much attention for research.  (Fric and Roshko, 1994) The main efforts were focused on detailed flow field features (such as velocity, pressure, or temperature), dimensionless variables, streamwise trajectory behavior and scaling properties (Andreopoulos and Rodi, 1984;Moussa et al., 1977). In addition, vorticity dynamics such as the Counter-Rotating Vortex Pair (CRVP) have also played an important role in flow mixing process, thus many studies have emphasized on the origin and dynamic evolution of CRVP (Smith and Mungal, 1998;Yuan et al., 1999;Broadwell and Breidenthal, 1984;Cortelezzi and Karagozian, 2001;Sykes et al., 1986;Kelso et al., 1996;Khali and Benmansour, 2009). Other vortex systems presented in the round transverse jet in crossflow include the horseshoe vortex formed upstream of the jet, the jet shear layer structures and the wake vortices that occur downstream of the jet for sufficiently higher jet-to-crossflow velocity ratios. For experimental JICF studies, flow visualization of this highly threedimensional flow has been an invaluable tool (Shan and Dimotakis, 2001;Krothapalli et al., 1990;New et al., 2004). Alternative jet nozzle shapes have also been received attentions to explore their impacts on the behavior of a passive scalar quantity such as temperature for mixing enhancement. In their experimental study, (McMahon and Mosher, 1969) observed from pressure field measurements that a jet issued from a long rectangular shape orifice placed in parallel to a crossflow stream (i.e., streamwise bus) will generate a jet in crossflow penetration distance longer than that of the same jet placed perpendicularly to the crossflow (i.e., blunt bus). Their measurements also suggested that the jet penetration of the rectangular orifice case is more important than that of the circular orifice. In contrary, the confined transverse slot jet of a rectangular duct has received less attention, despite that some studies were conducted (Jones and Wille, 1996;Kalita, 2002;Chen and Hwang, 1991;Haniu and Ramaprian, 1989;Ramaprian and Haniu, 1989;Holdeman, 1993). The studies considered slot jet slot that spans the entire dimension of the confinement (i.e., channel) and results showed that the generated flow field is nominally two-dimensional with respect to mean and turbulence statistics. It is not clear whether there is any CRVP-type structure observed in these confined transverse slot jet cases. A recent study by (Plesniak and Cusano, 2008) explored the flow field for a rectangular jet issuing into a confined crossflow field, with the jet width less than the depth of the crossflow duct. The study showed that the CRVP was established for most conditions with a slot jet spanned up to 80% of the duct depth.
To further advance JICF research, the primary motivation of the present work is to study the flow characteristics of a rectangular jet flow issuing into a main crossflow within a confined rectangular channel domain. The parameters used to characterize this configuration are the jet to crossflow velocity ratio r = U j /U cf and Reynolds number based on the jet bulk velocity and the rectangle width. Results of steady Reynolds-Averaged Navier Stokes (RANS) solution using three turbulence models will be analyzed and compared with previous Large-Eddy Simulation (LES) calculations of (Khali and Benmansour, 2009)

and
Science Publications AJAS experimental measurements of (Fougairolle, 2009). Finally a conclusion will be made on important findings obtained from this study.

Problem Background and Simulation Details
The flow configuration used in present study follows a previous large-eddy simulation performed by (Khali and Benmansour, 2009), in which a fullydeveloped turbulent jet flow issues perpendicularly from a rectangular duct into a crossflow channel domain. The flow conditions and geometry are the same as those from the experimental work conducted by (Fougairolle, 2009), in which a rectangular jet in crossflow in a closed tunnel has been studied. In their experimental study, the jet was heated up slightly for mixing study purposes, whilst various parameters including the velocity ratio were considered. The main investigations of this experiment focused on the interaction or not of the jet flow with both the opposite wall boundary layer and the adjacent wall where it issues, depending on values of this parameter. Its influence on the mixing properties and some other features of the dynamic behavior of the jet stream was also studied in their experimental work. Figure 2 shows a computational domain. A Cartesian coordinate system (x, y, z), representing streamwise, spanwise and transverse directions respectively, is adopted with its origin located at the centre of the rectangular jet exit plane on the bottom wall of the crossflow domain. The rectangular jet exit has dimensions of l j = 8 cm and h j = 5 cm, in x and y directions, respectively and is centred laterally on the bottom wall with a streamiwse distance of l x = 3l j from the crossflow inlet plane to the jet centre location. The dimensions of crossflow channel domain are L x = 13l j , L y = 6.25l j and L z = 7.5l j in x, y, z directions, respectively. Reynolds number based on the jet bulk velocity, the streamwise length of jet exit l j and viscosity of jet fluid is 26,000, same as that used in a previous LES study (Khali and Benmansour, 2009). A total of five velocity ratios are considered to evaluate its influence on flow characteristics; i.e., r = 3.3, 5, 7.3, 8 and 10. These values also follow a previous experimental investigation (Fougairolle, 2009), excepting for r = 5, 10 they are taken from the LES study (Khali and Benmansour, 2009). In order to study the mixing process of jet and crossflow, a passive scalar equation of temperature field is included in the computation.

Governing Equations
Due to low-speed conditions, the present calculations are based on steady state incompressible turbulent JICF. To allow study of flow and thermal mixing, the jet fluid is slightly heated, the difference of temperatures of the two fluid is ∆T max = 10°C. In these conditions the temperature of the jet is regarded as a passive scalar. With this assumption, the effects of density difference are neglected and without source terms (body forces). However, the equations of this JICF configuration are Equation 1 to 3: where, these equations are continuity, momentum and advection-diffusion of a passive scalar (temperature) respectively. In these equations, U j stand for mean velocity components, ρ for density and T is the temperature. The modified pressure P and the effective viscosity and diffusivity are defined as follows Equation 4 to 6: where, p is pressure, k is the turbulent kinetic energy, µ and Γ are the molecular viscosity and diffusivity, respectively. µ t is the eddy viscosity or turbulent viscosity, which must be modeled. Γ t is the eddy diffusivity written as: where, Pr t is the turbulent Prandtl number taken equal to 0.9 in the present simulation. The above equations can express turbulent fluctuations in terms of function of the mean variables only if the turbulent viscosity, µ t , is known. All the k-∈, k-ω and SST two-equation turbulence models use this variable.
The k-∈ model assumes that the turbulence viscosity is linked to the turbulent kinetic energy and dissipation rate via the relation as Equation 7: where, C µ = 0.09. The values of k and ∈ come directly from the differential transport equations for the turbulence kinetic energy and turbulence dissipation rate, respectively.
In the k-ω model, the turbulent viscosity is assumed to be linked to the turbulence kinetic energy and turbulent frequency via the relation as Equation 8: Finally, in the SST model, the turbulent viscosity is computed by the following relation as Equation 9: With a 1 being a constant, S is an invariant measure of the strain rate and F 2 is a function that is one for boundary-layer flow and zero for free shear flow.

Boundary Conditions
To ensure same 'mean' flow properties at the jet exit plane, a method of interpolating time-averaged LES results from a previous LES mesh (Khali and Benmansour, 2009) onto present RANS meshes is used for both 'mean' velocity and temperature profiles. For the crossflow inlet plane, the definition of flow conditions requires extra treatments and in the present simulation, the velocity profile is provided from a precursor simulation of an incompressible turbulent flatplate boundary layer matching the boundary layer thickness measured in the experiment at the same location, while a prescribed back pressure is imposed at crossflow domain outlet. No-slip and adiabatic wall conditions are used for velocity and temperature, respectively for top and bottom walls in the z-direction. A symmetry condition is imposed for two sidewalls in the y-direction.

Grids and Numerical Details
In RANS computation, three two-equation turbulence models are considered to assess their applicability to this flow configuration; they are k-ε model of Launder and Spalding (1972), k-ω model of Wilcox (2006) and SST model of Menter (1994), while a constant turbulent Prandtl number (0.9) used for computing turbulent diffusivity. In addition, three computational meshes are employed to achieve near grid independent solution. A simple progression functions is used for controlling gridpoint distributions. Three meshes are named coarse mesh (M 1 ), medium mesh (M 2 ) and fine mesh (M 3 ), respectively. All meshes are constructed in structured grid manner with grid refinement enforced around the jet exit and in nearwall regions. The meshes (M 1 , M 2 and M 3 ) have same topology and only differ in mesh density. Figure 3 shows a typical mesh, while Table 1 below summarizes all test parameters used in present simulations.

CFD Solver
The governing equations are discretized by a finitevolume method. The continuity, momentum and a passive scalar equation are solved in the fixed Cartesian directions on a co-located (non-staggered) grid layout such that the control volumes are identical for all transport equations. All the variables are thus stored at the center of the control volume. The velocity components at the control volume faces are computed by the Rhie-Chow interpolation method (Rhie and Chow, 1983) and the pressure-velocity coupling is handled by SIMPLEC method. The convective terms are treated by the hybrid scheme. The Algebraic Multigrid base algorithm is employed for solving the algebraic equations. The solution procedure is iterative and the computations are terminated when the sums of absolute residuals normalized by the inflow fluxes were below 10 −4 for all variables.

Grid Convergence Study
The effect of grid resolution on RANS predictions was studied by performing computations on three successive meshes, M 1 , M 2 , M 3 (see Table 1). Figure 4 displays predicted temperature and velocity profiles for r = 5 in a central XZ plane at y=0 with three turbulence models. The velocity magnitude is normalized by the jet bulk velocity U j and the normalized temperature T d is defined as T d = (T-T cf )/(T j -T cf ) where T j and T cf are the jet and the crossflow temperatures, respectively. Results were plotted at a streamwise location x/l j = 1 downstream of the jet centre. RANS predictions from different mesh resolutions were included for comparison, along with results of previous LES studies ( Khali and Benmansour, 2009;Rhie and Chow, 1983).

Turbulence Model Influence
Turbulence model influence at different levels is carried out, following a procedure proposed by (Sagaut and Deck, 2009). In the present case, two levels with the increased profoundness are investigated. They are 'mean' flow properties represented by 'mean' velocity which constitutes the first-order temporal statistics, proceeding with the second-order temporal statistics, i.e., the Turbulent Kinetic Energy (TKE). The results presented thereafter are from the fine mesh M 3 which is the same as the one used in LES. The velocity ratio is r = 10 and velocity and TKE profiles will be compared to those predicted by LES calculation for the same condition (Khali, 2010). All profiles are plotted along segment in the centre-plane XZ at x/l j = 0 which corresponds to the centre of the jet exit. Figure 5 depicts RANS predicted 'mean' velocity and turbulent kinetic energy profiles from three turbulence models in comparison with LES predictions. Figure 5 (lower graphs) gives turbulent kinetic energy profiles normalized by U j 2 at same streamwise locations as those 'mean' velocity profiles (Fig. 5 upper graphs).

Jet Trajectory and Decay of Passive-Scalar
The jet trajectory path presented here corresponds to a z-coordinate position z max where the maximum peak temperature is reached at each x-location in the central XZ plane (y = 0). Thus the trajectory presented is merely based upon a local maximum 'mean' passive-scalar (i.e., temperature in present study), rather than the velocity field as commonly used in the literature (Smith and Mungal, 1998)). Figure 6 displays simulated trajectories at five different velocity ratios while Fig. 7 shows the jet trajectory and the envelope of the jet for two velocity ratio of r = 5 and r = 7.3, respectively. The jet envelope is characterized by a parameter b, which is defined as b(x) = |z 1/2 (x)-z max (x)| for each vertical profile (Strzelecki et al., 2009). In this expression, z 1/2 (x) is z-coordinate on both sides of the jet stream 'tube' where temperature is equal to half a maximum value identified by z max (x).
To complete the analysis, results obtained by RANS approach are compared to LES results previously obtained by (Khali, 2010) for the same configuration. Figure 8 depicts RANS and LES (Khali, 2010) predicted trajectory at two velocity ratios r = 5, 10. Please note that in the case r = 10 the jet impinges onto the opposite wall.

Decay
Rate of Passive-Scalar and Comparison with Experiment Figure 9 shows maximum normalized temperature T dmax Vs curvilinear abscissa, denoted as 's', for four velocity ratios (r = 3.3, 5, 8 and 10). The curvilinear abscissa is calculated along the jet trajectory, i.e., a distance measured along jet trajectory path from the injection start point.

Fig. 9. Maximum temperature plotted with distance s normalized by lj
Re-plotting of the results in the log-log coordinates with s/l j parameter as suggested by (Smith and Mungal, 1998) is shown by Fig. 9b. Figure 10 shows further comparison of temperature decay rate for velocity ratios r = 3.3 and 10 with available experimental measurement of (Fougairolle, 2009) for the same configuration.

Mean Velocity and Temperature Profiles
Quantitative comparisons of steady RANS prediction with 'mean' LES results from reference paper (Khali, 2010) for velocity ratio r = 5 are carried out at three streamwise locations in the centre plane (y = 0) for velocity and temperature profiles, respectively (Fig. 11). Figure 12a and b give contours of passive-scalar (temperature) and velocity magnitude, respectively in a vertical plane (y-z) at a location of x/l j = 5.

Visualization of Flow and Mean Temperature Field
The passive scalar (temperature) field was further visualized in successive (y-z) planes of x/l j = 1, 3, 6. shows clearly that near grid independent solution has been achieved for medium mesh M 2 based on the fact that profiles were well-collapsed with those of fine mesh M 3 . For 'mean' velocity profiles (Fig. 5  upper graphs), influence of turbulence model is clearly seen where the k-ω model results give poorer 'mean' velocity profile in comparison to that of LES (Khali, 2010), than those from k-ε and SST models. This may be attributed to the known deficiency of k-ω model by its strong dependency on the freestream ω value (Wilcox, 1991). In fact, low level of inlet turbulence intensity of 1% for the crossflow is imposed in present simulations and using a ratio of turbulent and molecular viscosity rather than the ω value. Therefore, an overestimation of turbulent viscosity could be produced, which may result in an unphysical damping of spatial and temporal turbulent fluctuations. Comparing to the k-ω model, other two turbulence models, k-ε and SST, have shown overall good agreement with LES data. This means that a proper prediction of 'mean' quantities for a jet in crossflow also requires good spatial and temporal resolution of inherent dynamics. The second-order statistics of TKE are also compared with LES data (Khali, 2010) as displayed by Fig. 5 lower graphs. It is worth noting that all TKE profiles obtained by three turbulence models are overestimated in the vicinity of the jet exit, i.e., z/l j = 0, comparing to LES predictions. This trend could be due to the fact that unlike LES, RANS approach can only capture large-scale flow motions, rather than small-scale ones which are important in contribution towards turbulent kinetic energy distributions. Again, the TKE profile obtained by using the k-ω model shows a significant overestimation compared to that of LES by approximately a factor of two, while two other turbulence models over-predict the LES data by max 50%. The influence of the jet-tocrossflow velocity ratio parameter on the jet penetration shown in Fig. 6 in term of jet trajectories illustrate that a common feature of two distinct parts is clearly visible. From the origin of jet exit up to a certain abscissa (depending on velocity ratio), the lines are fairly grouped, except for a low velocity ratio of r = 3.3, indicating that the jet is penetrating in near straight manner. Beyond that position, the lines are bended in slopes and each follows its own evolution path depending on the velocity ratio. This behavior agrees experimental observations made by (Humber e al., 1993;Strzelecki et al., 2009) for a rectangular JICF, in which they noticed a unique penetration zone for different velocity ratio studied, followed by a change in the law of evolution along the centerline. The trajectories obtained in the present study also indicate the influence of the confinement (i.e., the upper wall), which imposes an earlier bending of the jet flow paths. Figure 7 shows that the jet spreads more quickly toward the bottom wall rather than the upper wall, comparing to the jet issuing into a crossflow of open space domain. It can be seen that jet trajectories obtained Science Publications AJAS by present RANS modelling are in good agreement with LES data in the near field region. After that, RANS predicted jet penetration depth into the crossflow is somewhat shorter than that of LES. This discrepancy increases with downstream distance and it may be due to the fact that in downstream region, the mixed jet/crossflow field is predominated by the newly formed CRVP, a well-known flow feature for this type of configuration. However, in this region RANS produced higher streamline curvature indicates anisotropic flow features, suggesting that Reynolds stresses may not be predicted accurately by two-equation turbulence models. Same issue was also discussed by (Demuren, 1993). Hence, an inadequate eddy-viscosity value is most likely computed by two-equation turbulence model that leads to an excess of turbulent diffusion, while LES approach uses a more universal sub-grid-scale model.
The passive scalar mixing properties is performed in this study by the decay of its maximum value along the jet center line (trajectory) as illustrated by Fig. 9. In Figure 9a, two regions are clearly shown; i.e., while close to the jet exit, a potential zone where temperature remains broadly constant and further downstream a diffusion zone characterized by a rapid decrease of temperature value. The length of the potential zone is equal to 1l j in case of r = 5 and 2l j in cases of r = 8 and 10. When re-plotting results in the log-log coordinates with s/l j , a clear slope in diffusion zone can be obtained that is dependent on velocity ratio. For velocity ratio r = 3.3 temperature decreases initially in s −1 law in the nearfield region, then later in s −2/3 law in the far-field region. Despite similar trends of decay in far field, present decay rate in s −1 law produced by RANS modelling is slightly different from those obtained by (Smith and Mungal, 1998) who derived an initial decay rate of s −1.3 law which is faster than present RANS prediction. This difference could be due to the fact that in experimental study of (Smith and Mungal, 1998), a top-hat velocity profile was used in defining the jet inlet conditions, compared to a fully-developed velocity profile used in present study. In fact, s −1 power-law decay was later obtained experimentally by Su and Mungal (2004), who used a fully-developed flow conditions at the jet exit. For velocity ratios r = 5, 8 and 10, temperature decreases in s −1/3 law in both near-field and far-field. The behavior observed at these velocity ratios is related to features observed in Fig. 6 and 7. As highlighted above, for r>3.3, the confinement imposes a greater bending constrain on the jets development, which limits their spreading towards the upper wall and induces an overall faster decay of the passive-scalar in the far-field for higher velocity ratio as shown in Fig. 9b. Figure 10 shows further comparison of temperature decay rate for velocity ratios r = 3.3 and 10 with available experimental measurement of Fougairolle (2009) for the same configuration. In general, RANS predicted results are in good agreement with these test data. Figure 11a shows that the LES predicted jet velocity remains almost constant in the vicinity of the jet exit up to a vertical position of about z = 1l j , after that it decreases rapidly to the crossflow velocity magnitude at z = 2.5l j . Similar behaviour is observed for temperature profile at this location (Fig. 11d). The jet evolutions at two downstream locations, x = 1l j , 3l j , are shown in Fig.  11 and 11c. For RANS results, both figures have shown clear peaks, inherited from higher jet exit velocity. Considering the velocity profiles (Fig. 11b), the first peak occur at z≈2l j near the exit (where the jet issues), indicating velocity increase from an imposed non-slip wall condition to a wake velocity at this position, whereas the second peak occurs at the same position as that of a temperature scalar profile, i.e., z≈3.5l j (Fig.  11e). This will be used to define the jet trajectory path explained in section 4 above. All these peaks appear inside the jet potential core. However for Fig. 11c and  11f, it can be seen that the second peak of the velocity profile appears at a position higher than that of temperature profile. To understand this observation, by comparison to the contours of passive-scalar (temperature) and velocity magnitude given in Fig. 12a and b, respectively, it can be seen that the maximum velocity occurs at a position z/l j = 0.3, whereas the passive-scalar maxima lies below this position (cf. Fig.  12a and b). This explains the reason why the jet trajectory computed with the maximum velocity magnitude lies above the one defined by the maximum passive scalar (temperature) concentration, as commonly used. Nevertheless, overall agreement of RANS with LES 'mean' data is generally acceptable.

AJAS
This near-field initiation of CRVP was suggested by various experimental observations (Smith and Mungal, 1998;Kelso et al., 1996;Khali and Benmansour, 2009). The effect of velocity ratio can also be seen in these Science Publications AJAS figures; i.e., while the velocity ratio increases with the size of the jet 'tube' increases. Furthermore, a more distinct separation of the two lobes in the passive scalar field occurs and the jet penetrates farther away from the wall where it issues into the crossflow field.

CONCLUSION
The interaction and evolution of a rectangular duct turbulent jet issuing perpendicularly into a crossflow has been numerically investigated using a commercial package ANSYS-CFX solver for several jet to crossflow velocity ratios and Reynolds number 26,000. The flow parameters have been chosen based on a previous incompressible Large-Eddy Simulation (LES). Three two-equation turbulence models were used; i.e., k-ε, k-ω and SST. Turbulent inflow data for the jet duct flow were generated by interpolating data from a previous LES simulation. This method was found to work well compared with instantaneous turbulent inlet velocity profiles at the jet plane from a fully-developed duct flow. This demonstrated the applicability and adequacy of the present approach with inflow-generation method for a turbulent jet flow interacting with a crossflow. RANSpredicted flow field results were found in qualitatively good agreement with previous published data, including a kidney shape structure (i.e., CRVP) captured for all velocity ratios considered, similar to those obtained by other researchers. At a given downstream location, the lobes of the kidney shape became more pronounced as the velocity ratio increases. The effect of domain confinement on the jet development was found to be complicated and also largely dependent on the velocity ratio. Despite that the first-order 'mean' quantities (such as velocity and temperature) predicted by RANS were in fairly good agreement with that from LES, all three turbulence models tested exhibited different levels of discrepancies in predicting the second-order temporal statistics such as the turbulent kinetic energy. Overall the RANS predictions using the k-ε model were in better agreement with the LES data than other two models. The RANS predicted passive-scalar decay along the jet centre line was also found in good agreement with available experiment data. This demonstrated the reliability of the RANS computation using two-equation eddy-viscosity models such as k-ε and SST in predicting 'mean' quantities for this jet in crossflow flow configuration. However, it is necessary to apply anisotropic turbulence model such as Reynolds stress model or high-fidelity large-eddy simulation method for better prediction of the second-order statistics in future work.

ACKNOWLEDGMENT
This study was conducted at Kingston University London and the University of the West of England, UK. The first author would also like to thank the Saad Dahleb University at Blida in Algeria for the funding support of this research project.