Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications

Two new two-equation eddy-viscosity turbulence models will be presented. They combine different elements of existing models that are considered superior to their alternatives. The first model, referred to as the baseline (BSL) model, utilizes the original k-u model of Wilcox in the inner region of the boundary layer and switches to the standard A>e model in the outer region and in free shear flows. It has a performance similar to the Wilcox model, but avoids that model's strong freestream sensitivity. The second model results from a modification to the definition of the eddy-viscosity in the BSL model, which accounts for the effect of the transport of the principal turbulent shear stress. The new model is called the shear-stress transport-model and leads to major improvements in the prediction of adverse pressure gradient flows.


Introduction
T HIS paper is concerned with twoequation eddyviscosity turbulence models with emphasis on an engineering per spective. It is based on the experience of the author in testing a large number of turbulence models against a wide variety of experimental test cases. The test flows cover a significant range of flow situations typically encountered in aerodynamic computations and are believed to allow some conclusions about a model's ability to perform in engineering applications. Two new turbulence models will be presented. They are based on a combination of what the author believes to be the best elements of existing eddyviscosity models.
There is a discrepancy between the large number of publica tions about twoequation models and the slow pace of im provement in accuracy that has been achieved since their intro duction. The basic problem of twoequation models, namely, their failure to correctly predict the onset and amount of separation in adverse pressure gradient flows, is still unre solved. Furthermore, there is no agreement on the standards by which to measure the improvement achieved by proposed new models, or alterations to existing models. Many times new models are based on theoretical concepts, which by themselves involve severe assumptions about the nature of turbulence, not even approximately satisfied in aerodynamic flows (homo geneous turbulence, small pressure gradients, low Reynolds number, flow equilibrium, etc.). It has been the author's experience that small changes (510%) in modeling constants can lead to a significant improvement (or deterioration) of model predictions. None of the available theoretical tools (dimensional analysis, asymptotic expansion theory, use of direct numerical simulations (DNS) data, renormalization group (RNG) theory, rapid distortion theory, etc.) can provide constants to that degree of accuracy. The only way to establish the validity of theoretical arguments under those conditions is to carefully test the resulting model against a number of challenging and welldocumented research flows. Unfortu nately, this is not general practice, and it is often unclear whether the improvements presented for one type of flow (e.g., boundarylayer flows) will not lead to a deterioration for another class of equally important flows (e.g., free shear flows). The author feels that the slow progress in engineering turbulence modeling, and the confusing picture it often pre sents, result to no small extent from an overemphasis of theoretical concepts and a virtual denial of the empirical na ture of the subject. Following an empirical approach, the author has developed two new turbulence models based on elements of existing models which are considered to be superior to their alterna tives. A description of these new models follows as well as an explanation of the rationale behind the choices that have been made in different areas of the flow and an address to antici pated criticism.
The A:co model 1 is the model of choice in the sublayer of the boundary layer. Unlike any other twoequation model, the A:co model does not involve damping functions and, as will be shown, allows simple Dirichlet boundary conditions to be specified. Because of its simplicity, the kco model is superior to other models, especially with regard to numerical stability. Furthermore, it is as accurate as any other model in predicting the mean flow profiles. Wilcox 1 has developed modifications that allow the treatment of rough walls and surface mass injection which can be used in the new model without change.
One point of criticism is that the &co model (like many other models) does not correctly predict the asymptotic behavior of the turbulence as it approaches the wall. However, the Taylor series expansion of the NavierStokes equations that underlies the analysis is only valid in the immediate wall proximity. So close to the surface the eddy viscosity is much smaller than the molecular viscosity and the asymptotic behavior of the mean flow profile is independent of the asymptotic form of the turbulence. Therefore, even if the turbulence model is not asymptotically consistent, the mean flow profile and the wall skin friction are still predicted correctly. A second point of criticism is that the k co model does not accurately represent the k and e distribution in agreement with DNS data. A significant number of damping functions have been developed in the last years for the k-e model which lead to an improved agreement with DNS data. In Refs. 2 and 3, a number of k-e models with different damping functions have been tested for a significant number of flows, with the conclusion that the specific form of the damping functions has little to no effect on the predicted velocity profiles and the skin friction of highReynoldsnumber flows. It should not be forgotten that the main (and often the only) information the mean flow solver gets from the turbulence model is the eddy viscosity. It is not clear why fitting the DNS data for k and e should lead to an improved eddyviscosity distribution. In the end, the agreement with DNS data might only be a matter of interpre tation. In the sublayer, Wilcox equates the quantity k in his model as being proportional to the normal component (with respect to the wall) of the turbulent kinetic energy. This inter pretation leads to a very good agreement with experimental and DNS data. In cases where the agreement with DNS data is considered important, the damping functions developed by Wilcox 4 can be applied to the present model.
The A:co model is also used in the logarithmic part of the boundary layer. It has been shown 1 ' 5 that the behavior of the A:co model in the logarithmic region is superior to that of the k-e model in equilibrium adverse pressure gradient flows and in compressible flows.
In the wake region of the boundary layer, the A:co model has to be abandoned in favor of the k-e model. The reason for this switch is that the A:co model has a very strong sensitivity to the (quiie arbitrary) freestream values co/ specified for co outside the boundary layer. It has been shown in Ref. 6 that the eddy viscosity in boundary and free shear layers can be changed by more than 100% by simply reducing the value of co/. It has also been shown in Ref. 6 that the k-e model does not suffer from this deficiency. There is no mathematical theory to date which distinguishes between twoequation models that suffer from the freestream dependency and those that do not. It is therefore of great importance that the influence of freestream values on the solutions of newly developed models is tested very carefully.
The mathematical analysis of the behavior of twoequation models in adverse pressure gradient flows has been largely restricted to the logarithmic region. 1 ' 7 Although the behavior of the model in the logarithmic region is of importance, espe cially in flows with moderate pressure gradients, it is the level of the eddy viscosity in the wake region that ultimately deter mines the ability of an eddyviscosity model to predict strong adverse pressure gradient flows. This has been clearly demon strated by the improvement that the JohnsonKing model 8 achieved over standard algebraic models by reducing the wake region eddy viscosity in adverse pressure gradient flows. The limited influence of the logarithmic region on the results for strong adverse pressure gradients is also evident in the failure of the original A:co model to accurately predict pressurein duced separation (as will be shown later) despite its superior logregion characteristics. The basic idea behind the Johnson King model is to enforce Bradshaw's observation that the principal turbulent shear stress is proportional to the turbulent kinetic energy in the wake region of the boundary layer. Enforcing this proportionality introduces a lag effect into the equations that accounts for the transport of the principal turbulent shear stress. It will be shown that the classical for mulation of the eddy viscosity in twoequation models violates Bradshaw's relation and thereby misses this important effect. In the new model the eddyviscosity formulation will be mod ified to take the transport effects into account.
Finally, in free shear layers away from surfaces, the stan dard k-e model will be utilized. There does not seem to be a model that accurately predicts all free shear flows (wake, jet, mixing layer) and the k-e seems to be a fair compromise.
To achieve the desired features in the different regions, the standard highReynoldsnumber version of the k-e model will be transformed to a A:co formulation. It will then be multiplied by a blending function (1 -FI) and added to the original A:co model times FI . The blending function F\ will be designed to be one in the sublayer and logarithmic region of the boundary layer and to gradually switch to zero in the wake region. This means that the new model will be based on a A:co formulation, with the original Wilcox model activated in the near wall region and the standard k-e model activated in the outer wake region and in free shear layers. This first step leads to a new model that will be termed the baseline (BSL) model. The BSL model has a performance very similar to that of the original &co model, but without the undesirable freestream dependency.
In a second step, the definition of the eddy viscosity will be modified to account for the transport of the principal turbu lent shear stress. The resulting model will be called the shear stress transport (SST) model. It is this second step that leads to a major improvement in performance over both the original &co and the standard k-e model.

Turbulence Model New Baseline Model
The idea behind the BSL model is to retain the robust and accurate formulation of the Wilcox k-u model in the near wall region, and to take advantage of the freestream independence of the k-e model in the outer part of the boundary layer. To achieve this, the k-e model is transformed into a &co formula tion. The difference between this formulation and the original A:co model is that an additional crossdiffusion term appears in the co equation and that the modeling constants are different (A small additional diffusion term is neglected in the transfor mation. It is shown in Ref. 9 that the term has virtually no effect on the solutions). The original model is then multiplied by a function F\ and the transformed model by a function (1 FI), and both are added together. The function F t will be designed to be one in the near wall region (activating the original model) and zero away from the surface. The blending will take place in the wake region of the boundary layer. The lefthand side of the following equations is the Lagrangian derivative: D/Dt: = d/dt + Original A:co model: Transformed k-e model: All constants, as well as the function F l5 are given in the Appendix.

Shear-Stress Transport Model
One of the major differences between eddyviscosity and full Reynoldsstress models, with respect to aerodynamic applica tions, is that the latter accounts for the important effect of the transport of the principal turbulent shear stress r -: -pu'v' (obvious notation) by the inclusion of the term DT dr D? = : Tt + l dx k (8) The importance of this term has clearly been demonstrated by the success of the JohnsonKing (JK) model. 8 Note that the main difference between the JK model and the CebeciSmith model lies in the inclusion of this term in the former, leading to significantly improved results for adverse pressure gradient flows. The JK model features a transport equation for the turbulent shear stress T that is based on Bradshaw's assump tion that the shear stress in a boundary layer is proportional to the turbulent kinetic energy k: with ai being a constant. On the other hand, in twoequation models, the shear stress is computed from: with Q = (du/dy). For conventional twoequation models, Eq. (10) can be rewritten to give: as noted in Ref. 10. In adverse pressure gradient flows the ratio of production to dissipation can be significantly larger than one, as found from the experimental data of Driver, 11 and therefore Eq. (11) leads to an overprediction of r. To satisfy Eq. (9) within the framework of an eddyviscosity model, the eddy viscosity is redefined in the following way: where F 2 is a function that is one for boundarylayer flows and zero for free shear layers. In an adverse pressure gradient boundary layer, production of k is larger than its dissipation (or Q>tf!co) and Eq. (12) therefore guarantees that Eq. (9) is satisfied whereas the original formulation v t -k/u is used for the rest of the flow.
To recover the original formulation of the eddy viscosity for free shear layers [where Bradshaw's assumption, expressed in Eq. (9) does not necessarily hold] the modification to the shearstress transport (SST) model is limited to wall bounded flows. This is achieved in the same way as it is for the BSL model by applying a blending function F 2 (also defined in the appendix). For general flows Q is taken to be the absolute value of the vorticity.

Model Versatility and Generality
The price for avoiding the freestream dependency and achieving the improved performance due to the modification in the eddy viscosity lies mainly in the necessary computation of the blending functions FI, F 2 and the additional crossdiffu sion term. The blending functions involve the distance from the surface which, however, has to be computed only once (as long as there is no grid deformation). Note that the distance from the surface is uniquely defined as being the shortest distance between the present point and all noslip boundaries (distance does not have to be measured normal to a surfacee.g., backwardfacing step). In most application codes, the boundary points have a marker and the computation of the distance function can therefore be automated. The increase in complexity from the Wilcox model to the present model is mainly in terms of coding. The overall computing time, as well as the stability of the code are not affected.
Once the model is implemented, it offers a wide variety of options. An example is a twolayer k-e model 12 with the origi nal k-u model in the sublayer and the k-e model in the high Reynoldsnumber region. This can be achieved by changing the argument of F\ for the BSL model from Eq. (A9) (see Appendix) to: The modification ensures that FI is zero for y + > 70. This twolayer k-e model utilizes the superior sublayer characteris tics of the £co model in much the same way that the model in Ref. 12 introduces an algebraic expression into the e equation. However, in the present approach the blending between the two regions is performed automatically and without user input.
The versatility of the model makes it possible to give the user a number of options, without making it necessary to program various models.

Numerical Method
The mean flow equations are solved by the INS3D code of Rogers and Kwak 13 which is based on a pseudocompressibility method. Important details about the discretization of the tur bulence model are given in Ref. 9. All computations have been performed on different grids to ensure that the presented solutions are grid independent. The airfoil computations were performed on a standard grid kindly provided by Rogers. 14 The standard k-e model is coded as given in Ref. 15.

Flat Plate Boundary Layer
To demonstrate the freestream dependency of the original k-u model, flat plate zero pressure gradient boundarylayer computations with different freestream values for co have been performed. For the first set of computations, the correct freestream values as given in Ref. 6 have been specified at the inflow boundary freestream for both the original and the BSL k-u model. Then, the preceding value was reduced by four orders of magnitude and the computations were repeated with both models. Note that the freestream value of k was also reduced to keep the freestream value of the eddyviscosity constant (the freestream value of the eddy viscosity has no influence, as long as it is small compared to its values inside the boundary layer). Figure 1 shows eddyviscosity profiles for the original and the BSL A:co model. The eddy viscosity of the original model changes by almost 100% due to the changes in co/, whereas the BSL model gives the same results for both cases. The strong sensitivity of the original model to co/ is clearly unacceptable and can lead to a severe deterioration of the results for complex flows, as will be shown later. Results of the SST model are also independent of co/. A more detailed study of the freestream dependency can be found in Refs. 6 and 9.
In each of the following comparisons between the different models, co/was always chosen according to the formula given in Ref. 6.
Zero pressure gradient flat plate boundary layer computa tions are given in Ref. 9. All models give good agreement with the experimental correlations for u + versus y + and c/. The A:co models can be run with the first gridpoint as far out as y+ = 3 without a deterioration of the results.

Free Shear Layers
For free shear layers the SST and the BSL models reduce to the same model (F\ = 0; F 2 = 0), and are virtually identical to the standard k-e model. Because the behavior of the k-e model for free shear layers is well known, and because of space limitations, results are not shown here, but can be found in Ref. 9. Reference 9 also shows the ambiguity of the results of the original A:co model 1 with respect to the freestream values.

Adverse Pressure Gradient Flows
One of the most important aspects of a turbulence model for aerodynamic applications is its ability to accurately predict adverse pressure gradient boundarylayer flows. It is especially important that a model be able to predict the location of flow separation and the displacement effect associated with it.
The test case most widely used to measure the performance of turbulence models under adverse pressure gradient condi tions is the flow reported by Samuel and Joubert. 16 Results for this flow are shown in Ref. 9 and are not reproduced here due to space limitations. It was found in Ref. 9 that all three A:co models reproduce the experimental data very well, whereas the JL k-e model gives values that are too high for c/.
The small differences between the solutions reported in Ref. 9, especially between the different A:co models, do not allow final conclusions about the abilities of the models to predict adverse pressure gradient flows. It appears that the Samuel Joubert flow does not pose a sufficiently strong challenge to the models to assess their potentials for these types of flows. The author has reached a similar conclusion in Ref. 10. It is therefore important to test models under more demanding conditions, with stronger adverse pressure gradients and possi bly separation. The following flowfield, reported by Driver, 11 has proven to be a highly selfconsistent and demanding test case.
In Driver's flow, a turbulent boundary layer develops in.the axial direction of a circular cylinder. An adverse pressure gradient is imposed by diverging wind tunnel walls and suction applied at these walls. The pressure gradient is strong enough to cause the flowfield to separate. The inflow Reynolds num ber is 2.8 • 10 5 based on the diameter D of the cylinder (140 mm). A 60 x 3 x 60 grid 10 was used for the present computa tions. A computation on a 100 x 3 x 100 grid gave almost identical results. Figure 2 shows the wall pressure distribution for Driver's flow as computed by the different models. The SST model gives superior results to the other models due to its ability to account for the transport of the principal turbulent shear stress. As expected, the JL k-e model produces the least accu rate results, with the BSL and the original £co model being close to each other in the middle. Figure 3, depicting the wall shearstress distribution for Driver's flow, shows that the SST model predicts the largest amount of separation, whereas the JL model stays firmly attached. Again, the BSL and the original A:co model produce very similar results.
The differences between the models can be seen in Fig. 4, which shows the velocity profiles. The SST model clearly produces the best agreement with the experiments. The larger displacement effect predicted by this model is reflected in the flattening of the c p distribution as was observed in Fig. 2. The original k co model predicts slightly better results than the BSL model, and the JL k-e model shows very little sensitivity to the pressure gradient, as was already reflected in Figs. 2 and 3.
The reasons for the different behavior of the models can be seen in the following two pictures.  /D = -0.091, 0.363, 1.088, 1.633, and 2.177. obviously predicts significantly higher shearstress levels than the other models in the region where separation is approached. This in turn leads to the firmly attached velocity profiles of Fig. 4. The differences between the models can also be seen from the eddyviscosity distributions. Figure 6 shows the max imum value of the kinematic eddyviscosity profiles for all streamwise (x) stations, nondimensionalized by u e d*. The SST model predicts the reduction of this quantity due to the ad verse pressure gradient in very good agreement with the exper iments. The BSL and the original A:co model are very close to each other up to separation (around x/D = 0), whereas the original model is closer to the experiments in the recovery region. Both models give consistently too high values for the maximum eddy viscosity in the adverse pressure gradient re gion. The k-e model falls only barely below the value of 0.0168 recommended by Clauser for equilibrium boundary layers (and used in the CebeciSmith model) and thereby fails to account for the nonequilibrium effects altogether.

Backward-Facing Step Flow
Results for the flow over a backwardfacing step as reported by Driver and Seegmiller 17 will be discussed next. This flow field was a test case in the 1981 Stanford conference for the evaluation of turbulence models. However, most of the com putations at the time were performed on comparatively coarse grids and there is substantial evidence that significantly finer grids have to be used to obtain gridindependent results. 18 The present computations have been performed on a 120 x 120 grid, with substantial grid refinement near the step. As with the other flowfields, a grid refinement study was made. The present results are virtually identical to those performed on a 90 X 90 and on a 240 x 240 grid. Figure 7 shows a comparison of computed and experimental skin friction distributions. The k-u models all perform signif icantly better than the k-e model. The reattachment length of  /D = -0.091, 0.363, 1.088, 1.633, and 2.177.   Figure 8 shows a comparison of the velocity profiles. All models fail to capture the relaxation downstream of reattach ment correctly. The results of Ref. 19 show that this is also true for a more complex model which accounts for anisotropy effects.

NACA 4412 Airfoil Flow
The following set of computations is for the flow around a NACA 4412 airfoil at 13.87 deg angle of attack. The Reynolds number with respect to the chord length is Re = 1.5210 6 . Experimental data for this flow have been reported by Coles and Wadcock. 20 The grid for the computations consists of 241 x 61 points and was made available by Rogers. 14 It is similar to the one used in Ref. 21. Figure 9 shows a comparison of the computed and the experimental velocity profiles at different streamwise stations. The results are similar to those for the separated case of Driver, Fig. 4 The final test case is the axisymmetric transonic shock wave/turbulent boundarylayer experiment of Bachalo and Johnson. 22 In this experiment, an axisymmetric boundary layer interacts with a shock wave created by a circular arc. It is beyond the scope of this paper to present a detailed study of transonic flows and only the highest Mach number case (M = 0.925) will be shown. The number of gridpoints used was 150 x 3 x 80. Grid independence was established by using different grids (129 x 3 x 60 and 180 x 3 x 100). Figure 10 shows the wall pressure distribution computed by the different models, compared with the experiment. The SST model pre dicts significantly better results than the other models, due to its improved transport features. Detailed comparisons for transonic flows will be presented in the future.

Conclusions
Two new turbulence models have been developed on a strictly empirical basis. They are combinations of what the author considers to be the best available elements of existing eddyviscosity models to date. Both models are based on a A:co formulation which is superior to other formulations with re gard to numerical stability. In a first step, a new baseline (BSL) model has been derived. It utilizes the original A:co model in the sub and loglayer and gradually switches to the standard k-e model in the wake region of the boundary layer.
The k-e model is also used in free shear layers. For boundary layer flows, the BSL model is very similar to the original &co model, but it avoids the strong freestream sensitivity of that model.
In the second step, a modification to the eddy viscosity has been introduced. It is based on the philosophy underlying the JohnsonKing model, which holds that the transport of the principal turbulent shear stress is of vital importance in the prediction of severe adverse pressure gradient flows. The re sulting model is termed the shearstress transport (SST) model.
Both models have been carefully fine tuned and tested for a large number of challenging research flows. The original £co, as well as the standard k-e model are included in the compar ison. As expected, the BSL model gives results very close to the original A:co model of Wilcox but avoids its freestream dependency. The SST model leads to a significant improve ment for all flows involving adverse pressure gradients and should be the model of choice for aerodynamic applications. It is the only available twoequation model that has demon strated the ability to accurately predict pressureinduced sepa ration and the resulting viscousinviscid interaction.
The new models require an increased amount of program ming effort compared to the original A:co model. However, once programmed, the new models consume only insignifi cantly more computing time and more importantly, they have proven to be very stable even in complex applications. 23 The concept underlying the new models is very flexible and lends itself to a multitude of different combinations. An example given in the text is a twolayer k-e model.
It is the author's conviction that a turbulence model has to be tested rigorously for a large number of flows, to establish the boundaries of its usefulness. Because of the limitations of the available theoretical tools and the severe assumptions in volved, this is also true for models based on more theoretical arguments. The new models are presently tested for transonic flows with very encouraging results. An early version of the SST model has been tested for complex threedimensional flows in Ref. 24. The results compare very favorably with the results of a full Reynoldsstress model, but significantly more testing in threedimensional flows will be necessary. (A16) Important detail! In applying this model, it is important that the reader be aware of the following ambiguity in the formulation of the production term of co for the SST model. The definition of the production term of co is sometimes written as: co du; P* = 7 T r u --(All)

Appendix: Baseline and Shear-Stress Transport Models
which introduces the nondimensional group v t (co//:) in front of the strain rate tensor. In the original and in the BSL model this group is equal to one and the two formulations for P w are therefore identical. This is not the case for the SST model because of Eq. (A14). The SST model has been calibrated with respect to Eq. (A2) and Eq. (A17) should therefore not be used.
The term arg! obviously goes to zero far enough away from solid surfaces because of the l/y or l/y 2 dependency in all three terms in Eq. (A9). Inside a boundary layer the three arguments in Eq. (A9) have the following purpose: the first argument is the turbulent length scale divided by y. It is equal to 2.5 in the log layer and goes to zero towards the boundary layer edge. The second argument ensures that FI is equal to one in the sublayer [note that co goes like l/y 2 in the near wall region and is proportional to l/y in the log region, so that l/(co.y 2 ) is constant near the surface and goes to zero in the log region]. The third argument is an additional safeguard against the freestreamdependent solution. It can be shown that the last argument ensures that argj goes to zero near the bound arylayer edge in case the "degenerate" solution given in Ref.
6 is approached. As arg! goes to zero near the boundarylayer edge, so does FI so that the standard k-e is used in that region. The following choice of freestream values is recommended: where L is the approximate length of the computational domain. The boundary condition for co at a solid surface is: where Ayi is the distance to the next point away from the wall. Equation (A 12) simulates the smooth wall boundary condition of Ref. 1 as long as 4y, + <3.

Shear-Stress Transport Model
The SST model is identical to the preceding formulation, except that the constants, 0i, have to be changed to: Set 1 (SST inner): where 0 is the absolute value of the vorticity. F 2 is given by: