A TURBULENCE MODEL TAKING INTO ACCOUNT THE LONGITUDINAL FLOW INHOMOGENEITY IN MIXING LAYERS AND JETS

The problem of potential core length overestimation of subsonic free jets by Reynolds-averaged NavierStokes (RANS) based turbulence models is addressed. It is shown that the issue is due to the incorrect velocity pro¦le modeling of the jet mixing layers. An additional source term in ω equation is proposed which takes into account the e ̈ect of longitudinal §ow inhomogeneity on turbulence in mixing layers. Computations con¦rm that the modi¦ed SpezialeSarkarGatski/Launder ReeceRodi-omega (SSG/LRR-ω) turbulence model correctly predicts the mean velocity pro¦les in both initial and far-¦eld regions of subsonic free plane jet as well as the centerline velocity decay rate.


INTRODUCTION
Many RANS-based turbulence models overestimate the potential core length of subsonic free jets [1].Analysis shows that this issue is due to incorrect velocity pro¦le modeling of the jet mixing layers [2,3].To solve this problem, a correction which increases the turbulence di¨usion intensity near the jet axis is proposed in [1].Unfortunately, this correction is not Galilean invariant.In [3], an invariant correction is suggested which enhances the turbulence di¨usion at the edges of turbulent regions and recovers the velocity pro¦le of a single stream mixing layer.
In this paper, a further investigation is conducted.Another way to reduce the jet potential core length in computations is proposed which is not connected to the modi¦cation of turbulence di¨usion coe©cients.Instead, an additional source term in characteristic turbulence frequency equation is proposed.The di¨erential Reynolds stress model (DRSM) SSG/LRR-ω [4] is taken as the basis for improvement.This is a modern DRSM with a great potential in modeling complex problems with nonequilibrium turbulence, separation, secondary §ows, and streamline curvature.The SSG/LRR-ω model uses a blending function to switch between the near-wall and free turbulent coe©cient values.In the present study, only ¤free turbulent part¥ of the model is modi¦ed.
The structure of the paper is as follows.In section 2, the equations of the original SSG/LRR-ω model are quoted as well as the results of self-similar computations of mixing layers and a plane jet.In section 3, a calibration of model coe©cients is conducted using the temporal mixing layer data.In section 4, an additional source term in ω equation is developed which improves the spatial mixing layer velocity pro¦le.Finally, in section 5, the computations using the complete Reynolds equation system closed by the modi¦ed turbulence model are presented and the improvements are shown which result from the changes made.

SSG/LRR-ω MODEL PERFORMANCE IN FREE TURBULENT FLOWS
The original SSG/LRR-ω model equations are as follows [4]: Hereinafter, summation over repeated indices is adopted; R ij = u ′ i u ′ j is the Reynolds stress tensor; ω is the characteristic turbulence frequency; k = R ii /2 is the turbulence kinetic energy; and ε = C µ kω is its dissipation rate (C µ = 0.09).Reynolds stress production is computed according to the exact formula: Turbulence kinetic energy production is P = P ii /2.Redistribution term š ij is computed using SSG model [5] in free turbulent regions and LRR model with coe©cients proposed by Wilcox [6] in near-wall regions.Introducing the Reynolds stress anisotropy tensor The coe©cients of the model V are determined as where the blending function F 1 has the form with d w being the distance to the nearest wall.
Coe©cient values for near-wall and free turbulent regions are presented in Table 1.
This turbulence model has been used to compute three free turbulent incompressible self-similar §ows.The ¦rst is temporal mixing layer which forms between two parallel streams with almost equal velocities |u 2 − u 1 | ≪ u 1 + u 2 [7].In the coordinate frame moving with velocity (u 1 + u 2 )/2 and y axis normal to the §ows, the mixing layer is symmetrical relative to its center, longitudinal gradients of every §ow variable, and transverse velocity component are negligible and mixing zone thickness linearly grows with time.To compute the self-similar temporal mixing layer velocity pro¦le, it is su©cient to solve only the longitudinal momentum equation supplied by the shear stress pro¦le: Let us introduce the self-similar variables τ = t and η = y/(u 0 t), where u 0 > 0 is the velocity of either stream in the chosen coordinate frame, and denote f (η) = u(t, y)/u 0 and r ij (η) = R ij (t, y)/u 2 0 .In these variables, Eq. ( 6) takes the form: Equations ( 1) and (2) used to close this equation are converted similarly.Two other §ows are the spatial single-stream mixing layer (Fig. 1a) and the far-¦eld region of a free plane jet (Fig. 1b).Both of them are inhomogeneous in longitudinal direction and contain nonzero transverse (ejection) velocity component v e .
In (ξ, η)-variables, Eqs.(8)(10) read: Equations ( 1) and ( 2) are converted similarly.A program has been implemented which solves Eqs. ( 7), (11), and ( 12) using time marching.An explicit second-order spatial scheme is adopted.At the external boundaries, very weak isotropic turbulence is speci¦ed: and w ∞ ≪ w max where κ max and w max are the maximal values of κ and w in the whole computational domain.It is checked that the solutions obtained within turbulent zones are insensitive to the κ ∞ and w ∞ chosen as described.
A mesh convergence study has been performed, resulting in the decision to use the meshes containing 150 nodes within the turbulent zone.
In Fig. 2, velocity pro¦les obtained with SSG/LRR-ω model are compared with the experimental data (see [9] for temporal mixing layer, the same as in [3] for the other §ows).The temporal mixing layer is plotted using the scaled coordinate η * = f ′ (0)η which makes the velocity slope in the center equal 1.From the data presented, the following conclusions can be made: (1) SSG/LRR-ω model predicts too sharp turbulent zone boundaries.This problem is evident in all the §ows considered; (2) the low-velocity part of the spatial mixing layer is wider than the highvelocity part.Note that this shortcoming is inherent to virtually all di¨erential turbulence models [3,6]; and (3) spatial mixing layer width is underestimated.
Sharp turbulent zone boundaries is a feature of the solution structure at the interface between turbulent and nonturbulent §ow regions (TNT interface) [10].Depending on turbulent di¨usion coe©cients C R , C ω , and α d , di¨erent power-law pro¦les of velocity and turbulence variables can be obtained near a TNT interface.The original SSG/LRR-ω model produces approximately linear velocity pro¦les.They can be smoothed by calibrating the turbulent di¨usion coe©cients.In the following section, this is done using the simplest §ow considered, namely, temporal mixing layer.

COEFFICIENTS CALIBRATION USING TEMPORAL MIXING LAYER DATA
Three series of temporal mixing layer computations have been conducted, in which coe©cients C R , C ω , and α d have been successively varied.The in §uence of these coe©cients is shown in Fig. 3. Coe©cient C R changes the velocity pro¦le very strongly.Its in §uence appears not only near the turbulent zone boundaries, but spreads deeply inside it.As a result, it is impossible to obtain an accurate velocity pro¦le which would fully conform to the experimental data by varying C R only.
Coe©cient C ω determines the velocity pro¦le behavior only near the turbulent zone boundaries and there is ¤saturation¥ at high C ω values, when velocity pro¦le becomes practically straight and does not change any further.To ¦nd the family of coe©cient sets giving the temporal mixing layer velocity pro¦le lying in the experimental data range, a series of parametric computations has been conducted.
In each computation, ¦xed values of C R and C ω have been set, and α d adjusted so that η * 0.99 coordinate (a boundary point where f = 0.99) would be equal, according to experiments, 1.36.As a result, for the values C R = 0.14, 0.18, 0.20, 0.22, and 0.26, the curves α d (C ω ) have been obtained, points on which correspond to temporal mixing layer velocity pro¦les with correct boundary positions.The diagram consisting of these curves is shown in Fig. 4.
Since the condition η * 0.99 = 1.36 does not fully determine the velocity pro¦le, among the coe©cient sets found, there are ones which give distorted velocity pro¦les in the regions between the center of the turbulent zone and its boundaries.The region of coe©cient sets giving velocity pro¦les entirely matching the experimental data is marked on the diagram with black circles.The distribution of the black circles on the diagram indicates that in the 0.18 ≤ C R ≤ 0.22 range, accurate temporal mixing layer modeling is possible.

TAKING INTO ACCOUNT THE LONGITUDINAL FLOW INHOMOGENEITY
The ¦rst step to calibrate the coe©cients of the SSG/LRR-ω model using the spatial single-stream mixing layer data is to determine the α d value provided C ω1 = 0.48 and C ω2 = 0.86 so that the mixing layer thickness D 0.1 computed by the points where F = 0.1 and 0.9 would be equal to the experimental value 0.165.The C ω1 and C ω2 values are taken from the middle of the ¤generally accepted¥ ranges 0.44 ≤ C ω1 ≤ 0.52 and 0.80 ≤ C ω2 ≤ 0.92 [6].Coe©cient α d appeared to be equal to 0.566.From the diagram (see Fig. 4), it follows that C R = 0.20 and C ω = 0.0593.The velocity pro¦le of single-stream mixing layer computed with these coe©cient values is presented in Fig. 5a (curve 1).It has the drawback which has already been noticed: the turbulent zone spreads too strongly towards the low-velocity region of the mixing layer and too weakly towards the highvelocity region.The potential core length estimation using this velocity pro¦le is L ini = 8.46 which is about 1.5 times higher than the experimental values 5 ≤ L ini ≤ 6.
It turns out that a turbulence model, which correctly describes the temporal mixing layer, at the same time strongly distorts the single-stream mixing layer.Hence, this distortion is due to the e¨ects absent in the temporal mixing layer: longitudinal gradients of velocity and turbulence variables and transverse velocity component connected to the gas ejection into the turbulent zone.To author£s knowledge, the only attempt to explicitly take into account these e¨ects in a turbulence model is paper [11] in which it is suggested to include an extra source term in the ε equation: where γ is the intermittency factor de¦ned as the probability of the point in space to be in the turbulent region at a given time [12].It increases the production of ε at the low-velocity boundary of a mixing layer and decreases it at the high-velocity boundary which reduces the turbulent zone spread into η < 0 region and enhances ¡ into η > 0 region.
In this paper, another source term is proposed which does not require the intermittency factor and is strictly Galilean invariant.In accordance with the observations made above, it depends on the longitudinal gradients of mean velocity and turbulence variables and on ejection velocity.Like the extra source term in Eq. ( 13), it should behave as an ¤odd function¥ (i.e., to be of di¨erent signs in the regions η < 0 and η > 0).This source term is included in the ω equa-tion which is empirical in nature and allows modi¦cations.On the other hand, the R ij equation ( 1) is left unchanged because it has an exact form derived from the NavierStokes eqation system.
To obtain the most general form of this term, let us restrict ourselves by using only mean velocity, k and ω gradients, but not individual R ij components.If necessary, these restrictions would permit us to adapt the source term to the eddy viscosity models in the future.The second derivatives are not involved for simplicity.Velocity vector itself is not used as it is not the Galilean invariant.
Using the above mentioned parameters, it is possible to form the following dimensionless aggregates: Velocity gradient norms are not used in the denominators of these formulas since they can vanish in the uniform §ow outside the turbulent zones.
In homogeneous turbulence and in the temporal mixing layer, all the expressions N kk , N kω , N ωk , and N ωω are identically zero.It can be shown that in self-similar spatial mixing layer, N kk ≡ 0 and N kω ≡ 0, but N ωk = 0 and N ωω = 0.Moreover, both N ωk and N ωω behave as an ¤odd function¥ relative to the mixing layer center.
Two extra source terms have been successively added to the ω equation: Negative signs in these formulas allow to use positive C ω3 and C ′ ω3 values to get additional ω production at the low-velocity boundary of the mixing layer and ω dissipation at the high-velocity boundary.It is found that I Analysis has shown that the dimensionless aggregate N ωk which enters the I (0) ω does not exceed 0.01 in the core of the mixing layer, but near its boundary, N ωk increases by several orders of magnitude.It is this uncontrollable growth of N ωk which leads to the instability; consequently, N ωk needs to be limited.The second reason for limiting is reduction of I (0) ω impact on the boundaries of turbulent zones, accurate description of which has been obtained earlier using the temporal mixing layer data.
Several modi¦cations of N ωk have been examined, which resulted in the following form of limiting: With this limiting, the restriction C ω3 21 vanishes.To compensate for some loss of e©ciency due to the limiting, C ω3 has been increased to 22.With the presented modi¦cation, the width of the mixing layer practically leaves unchanged (D 0.1 = 0.164) and L ini = 5.88.Velocity pro¦le in the scale of Fig. 5a is undistinguishable from the earlier obtained pro¦le with I ω .Now, one can say that the goal of correct modeling of the spatial single stream mixing layer is reached.
The computation of the plane jet has showed that without I ω , jet halfwidth η 0.5 computed by the point where f = 0.5 equals to 0.139 whereas mean experimental value is 0.105.With I ω , it reduces to 0.098, i. e., I ω increases the dissipation in the plane jet too much.
The following solution of this issue is proposed.It has been noticed that in self-similar mixing layer, N kω ≡ 0 but it is not zero in the plane jet.Let us use this additional term to independently calibrate the model against the plane jet data: The dimensionless aggregate N kω is added with the negative sign to compensate the action of N ωk .A series of computations has been conducted with the use of I ω with di¨erent β values.It is found that η 0.5 = 0.105 is achieved when β = 1.0; so, the ¦nal formulation of the extra source term is as follows: Plane jet velocity pro¦les obtained with (2) and without I ω (1) are depicted in Fig. 5b.With I ω , the pro¦le is slightly too sharp at the outer boundary; however, it falls within the experimental data range.

COMPLETE REYNOLDS EQUATION SYSTEM COMPUTATIONS
The SSG/LRR-ω turbulence model with the ability to enable the modi¦cations proposed in the paper has been implemented within TsAGI in-house code de- A plane cold subsonic free jet has been computed.Nozzle width based Reynolds number Re = u 0 h/ν is 6.7 • 10 6 and Mach number is 0.30.Computational domain, the mesh, and far-¦eld boundary conditions are shown in Fig. 6.The region of ¦ne mesh (of characteristic size 200h) is surrounded by bu¨er blocks to move the computational domain boundaries 500h1000h far from the nozzle exit section which is positioned in the coordinate frame origin.The lower computational domain boundary is the jet symmetry plane.
Turbulent boundary layer is modeled on the internal nozzle surface and contains 2530 cells across it.To avoid the KelvinHelmholtz type instability, longitudinal mesh re¦nement near the nozzle lip is not made.The jet mixing layer contains 60 cells across it.In the jet far ¦eld, the mesh is re¦ned in the transverse direction so that there are 70 cells across the jet half-width.Mesh convergence study has con¦rmed this mesh to be enough to accurately resolve the self-similar regions of the jet.
In Fig. 7, the longitudinal velocity §ow ¦eld obtained with the modi¦ed SSG/LRR-ω model is shown.The velocity pro¦le of the jet mixing layer in x/h = 3 cross section is compared to the self-similar pro¦le obtained earlier in Fig. 8a.The far-¦eld velocity pro¦les (x/h = 60 and self-similar) are depicted in Fig. 8b.The di¨erences between the pro¦les are insigni¦cant, which con¦rms the applicability of self-similar approximation to the analysis conducted in the paper.The most pronounced discrepancies are on the low-velocity boundary of the turbulent zones and are due to the absence of the vertical wall above the nozzle exit section in the complete Reynolds equation system computations.
Axial velocity u c distributions in the initial and part of the transitional region of the jet are shown in Fig. 9a.The modi¦ed model is compared to the original In Fig. 9b, axial velocity distributions in the far ¦eld of the jet are presented.First, it can be observed that (u 0 /u c ) 2 grows almost linearly obeying the selfsimilar relation, growth rate being in accordance with the experimental data.Second, the modi¦ed model corresponds to the experiments in x/h > 15 region while both the original SSG/LRR-ω and the SST models shift the velocity distribution downstream.
To ensure that the proposed modi¦cations do not negatively in §uence the description of boundary layers, the §at plate boundary layer computations have been conducted.In §ow Mach number is equal 0.8 and the plate is adiabatic.The length of the plate L is 1.5 m which corresponds to Re L = 1.7 • 10 7 .The outer boundaries of the computational domain are about 3L far from the plate.To analyze the boundary layer velocity pro¦le, a cross section is taken where momentum thickness Reynolds number Re θ is equal 16 • 10 3 .There are approximately 100 cells across the boundary layer, the height of the ¦rst cell above the wall y + being equal 0.1.
Skin friction distribution along the plate is shown in Fig. 10a.Again, the original SSG/LRR-ω model (1), its modi¦ed version (2), and the SST model (3) are used.In Fig. 10b, boundary-layer velocity pro¦le is presented in the cross section mentioned above.The experimental correlations are taken from [14,15].First, all the considered models give satisfactory results.The solutions obtained with the original SSG/LRR-ω and SST models are hardly distinguishable.Second, the modi¦ed SSG/LRR-ω predicts slightly higher skin friction and alters the outer part of the boundary layer.The near-wall region is left unchanged and mostly corresponds to the experimental data.

CONCLUDING REMARKS
It is shown that turbulence di¨usion models are not the main contributors to the discrepancies of the modern turbulence models in resolving the velocity pro¦le of the single-stream mixing layers.Obtaining the correct solution is possible with the use of an extra source term in ω equation which accounts for the longitudinal gradients of velocity and turbulence variables and transverse velocity component connected to the gas ejection into the turbulent zone.
The proposed modi¦cations to the SSG/LRR-ω turbulence model allow to obtain correct velocity pro¦les in the temporal and spatial mixing layers and in the free plane jet far ¦eld.It has been demonstrated in the computation of plane cold subsonic free jet that the modi¦ed model is superior to the original one in the whole §ow ¦eld.The description of §at plate boundary layer remains almost unchanged.
The ¦nal formulation of the model is as follows: The production term, redistribution term, and blending function are left unchanged and are given by formulas (3), (4), and (5), respectively.
The coe©cient values are given in Table 2.

Figure 1
Figure 1 Flow ¦eld scheme of single-stream mixing layer (a) and free plane jet far ¦eld (b)

Figure 2
Figure 2 Self-similar solutions obtained with the original SSG/LRR-ω model: (a) temporal mixing layer; (b) single-stream mixing layer; and (c) free plane jet far ¦eld

Figure 5
Figure 5 Velocity pro¦le of single-stream mixing layer (a) and free plane jet far ¦eld (b) before (1) and after (2) the introduction of the extra source term I (0) ω (a) or Iω (b) in the SSG/LRR-ω turbulence model with the modi¦ed coe©cient values

ω
and with C ω3 ≥ 20, it produces the mixing layer velocity pro¦le entirely matching the experimental data.The C ω3 = 21 value is chosen which gives the potential core length estimation L ini = 5.80 and mixing layer width D 0.1 = 0.165 without any extra turbulence model tuning.In Fig. 5a, velocity pro¦les are depicted without I (0) ω (1) and with it (2).A shortcoming of the I (0) ωsource term is that with C ω3 ≥ 22, the solution becomes unphysical: near the high-velocity mixing layer boundary, a running turbulent front appears on which instability develops.Moreover, even with C ω3 = 21 obtaining the stationary solution for the plane jet, far ¦eld becomes a challenge.

Figure 6
Figure 6 Computational domain, the mesh, and far-¦eld boundary conditions in the computation of the plane