Turbulence modeling means the formulation of the mathematical relationships required to obtain solutions of the averaged equations of motion. Averaging is necessary because the time-dependent *Navier-Stokes equation* cannot be solved analytically, and the range of scales occurring in turbulence limits the possibility of numerical solution by supercomputer to simple flow geometries and low Reynolds Numbers. The averaged momentum equations derived by Reynolds (1894) do not form a closed set because they contain new terms in the form of the *Reynolds stresses*,
. Additional equations must be devised to close the system by relating the new terms to known variables. At the lowest level of closure, the
are related to the mean rate of strain through the definition of an *"eddy viscosity"* analogous to the molecular viscosity of laminar flow. In the simplest models the eddy viscosity is determined from a *"mixing length" *
specified at every point in the flow field. One-equation models obtain the intensity of the turbulent velocity fluctuations from a modeled equation for the quantity k ≡
, the *turbulent kinetic energy* per unit mass, but are still limited by their reliance on a specified length scale. This limitation is removed by the introduction of a second modeled equation from which a length scale L can be calculated. Of many such proposals the equation for the rate of *turbulent energy dissipation* ε ≡ k^{3/2}/L has gained wide acceptance and the k-e model may be regarded as the standard method at the present time.

From its inception by *Boussinesq* (1877), it has been argued that the eddy viscosity concept is wrong in principle and its limitations have motivated the development of Reynolds-stress-transport models, or second-order closures, in which the
are obtained from their own modeled equations and the eddy-viscosity idea appears only in relating the triple products
to the gradients of
.

In *large-eddy simulations* (*LES*, or sub-grid-scale modeling), only the small-scale motion is modeled. The simulations are performed on grids coarser than those used for *direct numerical simulations (DNS)*, where the number of nodes needed to discretize all scales increases roughly as the cube of the Reynolds number. Even the low-Re direct simulations currently being attempted require some hundreds of hours of supercomputer time.

In the following pages, turbulence closures are described in increasing order of complexity; from simple mixing-length models for boundary layers to second-order closures based on the equations for . For reasons of space attention is restricted to the incompressible flow of uniform-property fluids, where most current research effort is concentrated. Several general works cover the field. The book by Launder and Spalding (1974) and the volume edited by Bradshaw (1976) are still useful. Surveys by Nallasamy (1987), Hanjalic (1988), Launder (1989,1993), and Speziale (1991) record recent developments.

The instantaneous velocity
and pressure
of a turbulent flow field comprise mean and fluctuating parts:
= U_{i} + u_{i} and
= P + p. (No distinction is made here between time and ensemble averaging.) Substitution in the continuity and Navier-Stokes equations for incompressible flow:

produces Reynolds averaged equations for the mean flow field:

The nonlinear terms of the original Navier-Stokes equations have produced nine components of (six of which are independent) so that, in the general case, there are now ten unknowns but still only four equations. Additional equations for are to be supplied by "modeling".

The Reynolds stresses may be related to the mean rate of strain by the turbulence constitutive equation:

in which ν_{t} is the "eddy viscosity". On dimensional grounds, and by analogy with kinetic theory, it is expressed as the product of a velocity fluctuation and a "mixing length": ν_{t} = υ
, although the mixing length, unlike the mean free path, is not small compared to the distances over which momentum is transported.
corresponds roughly to the scale of the energy-containing turbulent eddies, usually an order of magnitude smaller than a characteristic flow dimension. v_{t} is a property of the turbulence which varies from point to point. The objective is now to determine ν and
.

For certain free shear layers of width δ, it is possible to assume υ α (U_{max} – U_{min}) and
α δ to give

where C is a constant of order 0.1 which varies from flow to flow.

The need to specify in advance of the calculation effectively limits its application to thin shear layers described by the boundary-layer equation:

where U is the mean velocity in the predominant flow direction x, is the turbulent shear stress, and we assume steady flow with zero body forces. Prandtl’s (1925) model relates the eddy viscosity to the mean shear by:

where the modulus is needed to ensure positive values of ν_{1}.

In the inner region of the boundary layer on an impermeable plate, where V ≈ 0 and ∂U/∂x = –∂V/∂y ≈ 0, Equation (7) reduces to

Integration gives the total (laminar + turbulent) shear stress τ as

where τ_{w} is the shear stress at the wall. When dP/dx = 0, τ = τ_{w} in this inner region. Outside the viscous sublayer, where
> > > ν∂U/∂y, Equation (8) is

where u_{t} is a velocity scale defined by this equation. The scale of the energy-containing eddies is found to increase with increasing distance from the wall. If the variation is assumed to be linear and
= κy, Equation (11) becomes ∂U/∂y = u_{t}/κy to give, on integration, the *"law of the wall"*:

Extensive measurements in boundary layers show that κ ≈ 0.41 and C ≈ 5.2. In the outer region, y > 0.2δ, the influence of the wall is weak and, as in free shear flow, ≈ 0.09δ. Both regions are covered adequately by

with n = 5 for a boundary layer and n = 3 for fully-developed flow in a pipe of radius δ.

As y → 0,
= O(y^{3}),
= O(y^{3/2}) and t_{w} = μ∂U/∂y. The effects of viscosity on
are invariably accounted for by an exponential damping function. In *van Driest's (1956) original model*:

where y^{} ≡ yu_{t}/ν and A^{+} is a constant, ≈26 in constant-pressure flow. The mean-velocity distribution is obtained by integrating:

where u^{+} = U/u_{t},
^{+} =
u_{t}/ν, τ^{+} = τ/τ_{w} ≈ 1 for constant-pressure flow. u^{+} = y^{+} in the viscous sublayer, y^{+} → 0. Equation (15) reduces to du^{+}/dy^{+} = 1/κy^{+} for large y^{+} from which the law of the wall, Eq. (12), is recovered. For κ = 0.41, C ≈ 0.204 A^{+} ≈ 5.3.

The *effects* of *pressure gradient* are at least partially accounted for by writing A^{+} = 26/(τ^{+})^{n} where τ^{+} = τ/τ_{w}, and 0 < n < 2. For large |p^{+}|, and particularly for strongly accelerated flows, the Cebeci-Smith (1974) formula may be used:

McEligot (1985) has made a collection of several such formulas. Typical mixing-length methods are described by Baldwin and Lomax (1978), Cebeci et al. (1984) and Granville (1987).

In the one-equation models attributed to *Kolmogorov* (1942) and *Prandtl* (1945), the velocity scale ν is taken more realistically as
√k, and k ≡
is obtained from an equation derived by multiplying Equation (2) by u_{i}, averaging, and dividing by 2. The time-averaged result is:

where the terms represent

mean-flow convection,

mean-shear production, P,

viscous dissipation, ε,

diffusion by (a) turbulence, (b) viscosity, (c) pressure.

*Turbulent diffusion* is represented by:

with σ_{k} ≈ 1.0. Viscous diffusion is exact. Pressure diffusion has long been assumed to be negligibly small or, alternatively, to be incorporated in Eq. (18). Instead of
a new length scale is defined, on dimensional grounds, as

and the eddy viscosity is written as

where C_{μ} is a constant to be determined. In the constant-stress wall layer ν_{t} =
u_{g}t,
and ε ≈ P ≈
. Then

Shear-flow data suggest
to give C_{μ} ≈ 0.09 and in the constant-stress layer, consistency with the law of the wall requires:

The effects of viscous damping are expressed by

where y* =
, α_{μ} ≈ 2.5 and, in the Hassid and Poreh (1975) model for example, A_{μ} = 84 [see also Gibson et al. (1978)]. The limiting behavior at the wall is

As y → 0, viscous diffusion is balanced by the dissipation and the k-equation, Equation (17), reduces to ε = v∂^{2}k/∂y^{2}. Then, by virtue of Equation (25) the value of the dissipation at the wall may be written as

In the Hassid and Poreh (1975) model:

One-equation models are also described by [Wolfshtein (1967) and by Norris and Reynolds (1975)].

The length scale need not be determined when it can be evaluated at every point as the solution of a second equation. The length scale itself is not a suitable dependent variable, not least because it actually increases as turbulence decays. Several other variables have been suggested for this duty of which the most widely used is the dissipation rate, ε = k^{3/2/L}. In the k-ε model, Equation (20) for the eddy viscosity is rewritten as

where f_{μ} is the low-Re damping function multiplying L.

The exact equation for ξ contains numerous intractable terms [see Rodi and Mansour (1993)] but it may be condensed to the basic form containing the four essential processes of advection, diffusion (viscous, turbulent and pressure), production, dissipation:

If ε is interpreted as the rate at which turbulence kinetic energy is transferred from large scales to small in the spectral cascade, its rate of transfer can be expressed in terms of the time scale t = k/e:

where C_{ε2} is considered to be a constant, at least for turbulence at high Reynolds numbers. The production of ε by vortex stretching ought arguably to be modeled in terms of turbulence quantities, but this is not feasible in two-equation modeling. Instead, for high-Re turbulence, it is assumed that

A gradient-diffusion process is assumed for turbulent transport:

and the modeled ε-equation is finally written as

where the functions f_{1} and f_{2} are introduced to deal with low-Re effects on production and dissipation, and symbol E stands for additional terms which become important in the immediate vicinity of the wall.

The decay constant C_{ε2} is determined by reference to the decay of homogeneous grid turbulence for which the k and e Equations (17) and (34) reduce to:

Measurements of k downstream of the grid are fitted by the power law k = k_{0}x^{−α}, with α ≈ 1.25 in the initial period. Substitution in Equations (35) and (36) gives C_{ε2} = (α + 1)/α ≈ 1.8. However Launder and Spalding's (1974) C_{ε2} = 1.92 has stood the test of time and endures in the "standard" k-ε model.

A relationship between the three constants of the ε-equation is found from conditions in the constant-stress wall layer where P ≈ ε ≈ . Convection by the mean flow is negligibly small and Equation (34) reduces to:

from which, with dk/dy ≈ 0

Launder and Spalding (1974) found that the predicted spreading rates of free shear layers were highly sensitive to (C_{ξ2} – C_{ξ1}) and that "best-fit" values were: C_{ε1} = 1.44, C_{ε2} = 1.92, σ_{ε} = 1.3.

Because very fine grids are needed to resolve the small-scale turbulence in the viscous layers close to a wall, and because wholly satisfactory low-Re models have yet to be developed, these regions are often bridged by "wall functions" based on the universal laws of the wall. Alternatively, the problem of adapting the e-equation may be avoided by using a one-equation model in which the length-scale distribution is prescribed through the viscous layers, as advocated by Rodi (1991) and Rodi et al. (1993).

As the wall is approached the k-equation (17) reduces to ε = ν∂^{2}k/∂y^{2}. This result causes problems which Jones and Launder (1972) avoid by using an equation for the "isotropic" part of the dissipation, which is zero at the wall:

A tendency in more recent work has been to retain e as the dependent variable while defining the time and length scales with , thus: ; . As y → 0, from Equation (25)

When this is compared with Equation (26), ε/ν = 2α + 4βy + O (y^{2}), it is seen that
= 0(y^{2}).
is proportional to y at the wall, as it is in the logarithmic layer where
→ ε, and f_{μ} in Equation (29) ought strictly to be proportional to y as y → 0, although this condition is seldom satisfied in practice. Rodi and Mansour (1993) compare various f_{μ} functions with the results of numerical simulations of channel flow. Their best-fit function is

The drawback to this and many of the other functions in the literature is that y^{+} is unsuitable for general flow calculations. The Launder and Sharma (1974) formula is better adapted to general purpose calculations in that it is based on the local turbulence Reynolds number, R_{t} = k^{2}/νε, although it compares poorly with the DNS channel-flow results:

In nearly all current models the ε-production rate is left unchanged by setting f_{1} = 1.0. The dissipation-rate function f_{2}, must fit measurements of the decay of grid turbulence in the final period, k α t^{−α}, with α ≈ 2.5. Then C_{ε2}f_{2} = (α + 1)/α ≈ 1.4, whence f_{2} ≈ 1.4/1.8 ≈ 0.78 as R_{t} → 0. According to Hanjalic and Launder (1976) the data of Batchelor and Townsend (1948) are fitted by:

where R_{t} ≡ k^{2}/νε. Viscous effects are confined to the layers very close to the wall, e.g. f_{2} ≈ 0.99 for R_{t} = 10, and R_{t} is to be preferred to y^{+} as the wall-distance parameter.

Different forms of the term E in Equation 34 are listed in the survey by Rodi and Mansour (1993). The most widely used has been Jones and Launder's (1972):

which, according to Rodi and Mansour (1993), compares unfavourably with near-wall values of the exact production term:

in the low-Re channel-flow simulations of Kim et al. (1987).

A passive scalar (temperature) is expressed as the sum of mean and fluctuating parts, T + , which is substituted in the -conservation equation. The result upon averaging is

where α is the molecular diffusivity and the
are turbulent scalar fluxes. An eddy diffusivity α_{t} is defined by

and is related to ν_{t} via a "turbulent Prandtl number" Pr_{t} ≡ ν_{t}/α_{t}.

Near a wall, for small ∂T/∂x, the counterpart of Eq. (10) is

where is the heat flux through the wall when scalar T is the temperature. The damped distribution is like that of , Eq. (14):

and the T^{+} (≡ –ρu_{τ}c_{p}(T – T_{w})/
) profile is obtained by integrating:

to give using Eq. (15), for q^{+} = 1 and large y^{+}:

Typically [Gibson (1990)]
= 0.45,
= 0.91. With
= A^{+} = 26, C
= 3.0 when Pr = 0.7 (for air). A temperature wall function

is used in two-equation modeling with

obtained from pipe-flow data by [Jayatillaka (1969)].

The equation for k ≡ is obtained by multiplying the (T + ) equation by , dividing by two, and averaging. Then

This is like the k-equation, Equation (17), without the pressure term. Turbulent transport is modeled by gradient diffusion. The dissipation rate may be represented by ε = k ε/kR, where R is the ratio of time scales which is found to range from about 0.5 to 1.0 in different flows. It is difficult to devise an equation for ε because there are two relevant time scales. The problems are identified by Launder (1976) and Lumley (1978). Modeled ε -equations are described by Jones and Musonge (1988), Malin and Younis (1990), Craft and Launder (1991), Youssef et al. (1992), and Sommer et al. (1993), mainly in connexion with second-moment closures.

Although two-equation models are capable of giving good results for thin shear layers and some recirculating flows, they fail to account adequately for the effects of complex strain fields and body forces. These limitations, combined with the acknowledged defects of eddy-viscosity modeling generally and the availability of ever-increasing computer power, have motivated the development of second-order closures based on modeled transport equations for the Reynolds stresses themselves. The equations for
are obtained by multiplying the equations for (U_{i} + u_{i}) and (U_{j} + u_{j}) by u_{j} and the u_{i}, respectively, adding and averaging. The result is

Models are required for viscous dissipation, turbulent and pressure diffusion, and energy redistribution by pressure-strain. The closure ought strictly to satisfy "realiseability" constraints which do not permit negative normal stresses, or correlations that fail to satisfy Schwartz's inequality [Lumley (1978)]. These conditions have often been disregarded in practice.

Rotta's (1951) ideas were developed by Daly and Harlow (1970), Hanjalic and Launder (1972), Launder et al. (1975) and Lumley (1978). Recent advances are described by Shih and Lumley (1986), Speziale et al. (1991), Launder and Tselepidakis (1993) and Durbin (1993). Launder (1989) and So et al. (1991) survey the field.

When the turbulence Reynolds number is sufficiently high the small-scale motion is nearly isotropic and the dissipation may be assumed to be equally divided between components. The form

also allows for viscous energy redistribution when these conditions do not obtain and the term in square brackets is nonzero. ξ is to be obtained as the solution of a modeled transport equation as in two-equation modeling.

with C_{s} = 0.11. A shorter form attributed to Daly and Harlow (1970), with only the last term, is also used with C_{s} = 0.22.

The fluctuating pressure appears in "diffusion" and "redistribution" terms on account of the manipulation

where the last term has the useful property of contracting to zero in incompressible flow. It therefore serves only to redistribute energy between components. may be expressed as the sum of two terms in a volume integral. One of these terms contains only averaged fluctuating quantities; the second involves the mean rate of strain.

Models for the turbulence, or "slow" part of rest on the proposition that strain-free turbulence tends to the isotropic state. For the decay of turbulence downstream of a wind-tunnel contraction the energy and stress Equations (17) and (57) combine to give

where there is no summation on repeated Greek indices. In Rotta's (1951) widely used linear model

where C_{1} is a constant and k/ε is the time scale. C_{1} values in use range from 1.5 to 3.0. There is no return to isotropy when C_{1} = 1. The evidence is that the process is nonlinear and that the tendency is weak for small anisotropy. The nonlinear model described by Lumley (1978) has two coefficients which are generally functions of the turbulence Reynolds number, k^{2}/νε, and the invariants b^{2} = b_{ij}b_{ij} and b^{3} = b_{ij}b_{jk}b_{ki} of the anisotropy tensor

Sliih and Lumley (1985) use Equation (62) with:

where

is a parameter which has the useful properties (for near-wall Reynolds-stress modeling) of being unity in isotropic turbulence and zero in the two-dimensional limit [Lumley (1978)].

A simple but effective model for the mean-strain, or "rapid", part of is

where the right-hand side is proportional to the anisotropy of production. This term is actually one of three components formally obtained by Launder et al. (1975), but good results can be obtained by this abbreviated form alone. On combining the two components, generalized to cover the shear stresses, there results

For two-dimensional thin shear layers close to local equilibrium (P ≈ ε) the model gives the following algebraic expressions for the stresses:

where φ ≡ (1 – C_{2})/C_{2}. Evidently φ must be approximately 0.2 to give
≈ 0.5, C_{μ} =
≈ 0.1. Gibson and Launder (1978) chose C_{2} = 0.6 in order to satisfy an exact result for isotropic turbulence, setting C_{1} = 1.8 to obtain φ = 0.22. Gibson & Younis (1986), got better results for curved and swirling flows with C_{1} = 3.0, C_{2} = 0.3, disregarding the results of rapid-distortion theory and achieving arguably a better fit to the return-to-isotropy data [but see Taulbee (1989) for a contrary assessment].

A rigid boundary reflects pressure fluctuations and thereby inhibits the transfer of energy into the component normal to the wall. Launder et al. (1975), followed by Gibson and Launder (1978), add wall-effect terms to both parts of their pressure-strain models, in each case expressed as linear functions of distance from the wall. Craft and Launder (1991) retain only the wall correction to the rapid component, and that in more complicated form. Speziale et al. (1991) and Durbin (1993) have devised models that require no explicit wall-reflection terms.

In the near-wall viscous region where the Reynolds number is low, the energy-containing and dissipation ranges of the energy spectrum overlap, the fine-scale motion is no longer isotropic, and the dissipation rate is no longer divided equally between components. Rotta (1951) initially assumed that ε_{ij}/ε ≈
but it is now known that this approximation seriously underestimates ε_{22}, and it has been corrected in later work. The ε-equation differs only slightly from the form used in two-equation modeling. In the immediate vicinity of the wall
becomes so attenuated that the turbulence becomes almost two-dimensional. The models devised by Launder and Shima (1989), Shima (1993), Launder and Tselepidakis (1993) involve Lumley's (1978) "flatness factor". A, of Eq. (65) which is unity in isotropic turbulence but → 0 as the turbulence approaches two-dimensionality in the near-wall layers.

The equations have received less attention than the stress equations, and the closure methods are less highly developed, though naturally similar. The survey by Launder (1976) remains a useful introduction, recent advances are described by Lai and So (1991) and Durbin (1993).

The equations are obtained by multiplying the equation for (T +
) by u_{j} and adding it to the (U_{i} + u_{j}) equation multiplied by
. The result after averaging is

At high Reynolds numbers the molecular diffusion terms are negligible, and the dissipation is zero by virtue of the local-isotropy hypothesis. Closure approximations for turbulent diffusion and the pressure-scalar-gradient correlation are akin to those of the Reynolds-stress equations, but are further complicated by the need to allow for different time scales in the fluctuating scalar and velocity fields. Turbulent diffusion of is modeled by

with C_{s}
≈ 0.11 or, frequently, with only the last of these terms. The arguments which led to the modeled pressure-strain term of the stress equations suggest for the equivalent term here:

where the first term on the right is the counterpart of Equation (62) and the second involves the mean-strain part of production.

For the idealized local-equilibrium shear layer, with negligible mean-flow or turbulent transport, the equations reduce to [Durbin (1993)]

giving

C
_{1} = 3.0 and C
_{2} = 0.33 are chosen on the basis of the experimental data for free flow examined by Launder (1976): Pr_{t} ≈ 0.7,
≈ –1.3. The wall-reflection correction is small; substitution of near-wall stress values from Eq. (68) gives Pr_{t} ≈ 0.8,
≈ –1.8. Lai and So (1990) and Durbin (1993) describe low-Re scalar-flux models.

#### REFERENCES

Baldwin, B. S. and Lomax, H. (1978) AIAA paper, 78–257.

Batchelor, G. K. and Townsend, A. A. (1948) *Proc. Roy. Soc*, A194, 538.

Boussinesq, J. (1877) *Mem. Pre. par div. Sav.*, 23, Paris.

Bradshaw, P., (ed) (1976) *Topics in Applied Physics*, 12, Springer.

Cebeci, T. and Smith, A. M. O. (1974) *App. Math. Mech.*, 15 Academic Press.

Cebeci, T., Stewartson, K, and Whitelaw, J. H. (1984) *Numerical and Physical Aspects of Aerodynamic Flows II*, T. Cebeci, Ed., 1, Springer.

Craft, T. J. and Launder, B. E. (1991) *Proc. 8th Turbulent Shear Flows Symposium*, Munich.

Daly, B. J. and Harlow, F. H. (1970) *Phys. Fluids.*, 13, 2634.

Durbin, P. A. (1993) *J. Fluid Mech.*, 249, 465–498.

Gibson, M. M. and Launder, B. E. (1978) *J. Fluid Mech.*, 86, 491.

Gibson, M. M. and Younis, B. A. (1986) *Phys. Fluids*, 29, 38.

Gibson, M. M. (1990) *Near-Wall Turbulence*, S. J. Kline and N. H. Afghan, Eds., 157, Hemisphere.

Gibson, M. M. and Rodi, W. (1989) */ Hydraulic Research*, 27, 233.

Gibson, M. M., Spalding, D. B. and Zinser, W. (1978) *Letters in Heat and Mass Transfer*, 5, 73.

Granville, P. S. (1987) *AIAA J.*, 25, 1624.

Hanjalic, K. and Launder, B. E. (1972) *J. Fluid Mech.*, 52, 609.

Hanjalic, K. and Launder, B. E. (1976) *J. Fluid Mech.*, 74, 593.

Hanjalic, K. (1988) *Near-Wall Turbulence*, S. J. Kline and N. H. Afghan, Eds., 762–781, Hemisphere.

Hassid, S. and Poreh, M. (1975) *ASME J. Fluids Eng.*, 100, 107.

Jayatillaka, C. V. L. (1969) *Prog. Heat Mass Transfer*, 1, 193.

Jones, W. P. and Launder, B. E. (1972) *Int. J. Heat Mass Transfer*, 15, 301.

Jones, W. P. and Musonge, P. (1988) *Phys. Fluids*, 31, 3589–3604.

Kim, J., Moin, P. and Moser, R. (1987) *J. Fluid Mech.*, 177, 133.

Kolmogorov, A. N. (1942) *Izv. Akad. Nauk. SSSR Ser. Phys.*, 6, 56.

Lai, Y. G. and So, R. M. C. (1990) *Int. J. Heat Mass Transfer*, 33, 1429.

Lang, N. J. and Shih, T. H. (1991) *NASA TM* 105237.

Launder, B. E. and Sharma, B. I. (1974) *Letters in Heat and Mass Transfer*, 1, 131–138.

Launder, B. E. (1976) in Bradshaw 1976.

Launder, B. E. (1989) *Int. J. Heat Fluid Flow*, 10, 282.

Launder (1993) *Engineering Turbulence Modeling and Experiments*, 2 W.

Rodi and F. Martelli, Eds., 3-22, Elsevier.

Launder, B. E. and Shima, N. (1989) *AIAA J.*, 27, 1319.

Launder, B. E. and Spalding, D. B. (1974) *Computer Methods in Applied Mechanics and Engineering*, 3, 269.

Launder, B. E. and Tselepidakis, D. P. (1993) *Turbulent Shear Flows*, 8, F. Durst et ai., Ed. 81–96, Springer.

Launder, B. E., Reece, G. J. and Rodi, W. (1975) *J. Fluid Mech.*, 68, 537.

Lumley, J. L. (1978) *Adv. Appl. Mech.*, 18, 123.

Malin, M. R. and Younis, B. A. (1990) *Int. J. Heat and Mass Transfer*, 33, 2247–2264.

McEligot, D. M. (1985) *Lecture Notes in Physics*, 235, 292-303, Springer.

Nagano, Y. and Shimada, M. (1993) *9th Turbulent Shear Flows Symposium*, Kyoto.

Nallasamy, M. (1987) *Computers and Fluids*, 15, 151.

Norris, L. H. and Reynolds, W. C. (1975) *MED Rep. FM-IO*, Stanford.

Prandtl, L. (1925) *ZAMM*, 5, 136.

Prandtl, L. (1945) *Nachrichten Akad. Wissenschaft Gottingen*.

Reynolds, O. (1894) *Phil. Trans. Roy. Soc*, 186 A, 123.

Rodi, W. (1991) AIAA paper 91-0216.

Rodi, W. and Mansour, N. N. (1993) *J. Fluid Mech.*, 250, 509.

Rodi W., Mansour, N. N. and Michelassi, V. (1993) *Trans ASME J. Fluids Eng.*, 115, 196.

Rotta, J. C. (1951) *Z. Phys.*, 129, 547 and 131, 51.

Shih, T.-H. and Hsu, A. T. (1991) AIAA paper 91-0611

Shin, T-H., Lumley, J. L. and Chen, J-Y. (1990) *AIAA J.*, 28, 610.

Shih, T-S. and Lumley, J. H. (1985) *Cornell Rep.* FDA-85-3.

Shih, T-S. and Lumely, J. H. (1986) *Phys. Fluids*, 29, 971.

Shima, N. (1988) *Trans. ASME J. Fluids Eng.*, 110, 38.

Shima, N. (1993) *Proc. 9th Turb. Shear Flows Symp.*, Kyoto.

Shir, C. C. (1973) *J. Atmos. Sci.*, 30, 1327.

So, R. M. C, Lai, Y. G., Zhang, H. S. and Hwang, B. C. (1991) *AIAA J.*, 29, 1819.

Speziale, C. G. (1991) *Ann. Rev. Fluid Mech.*, 23, 107.

Speziale, C. G., Abid, R. and Anderson, E. C. (1992) *AIAA J.*, 30, 2, 324.

Speziale, C. G., Sarkar, S. and Gatski, T. B. (1991). *J. Fluid Mech.*, 227, 245.

Sommer, T. P., So, R. M. C. and Zhang, H. S. (1993) AIAA paper 93-0088.

Taulbee, D. B. (1989) *Recent Advances in Turbulence*, (Ed. W. K. George and R. Arndt, Eds., 75, Hemisphere.

Van Driest, E. R. (1956) *J. Aero. Sci.*, 23, 1007.

Wolfshtein, M. (1967) Ph.D. thesis, Imperial College.

Youssef, M. S., Nagano, Y. and Tagawa, M. (1992) *Int. J. Heat Mass Transfer*, 35, 3095.