12.10. Turbulence Modeling

The closure problem arising from Reynolds-averaging of the equations of fluid motion has lead to the development of approximate models to close systems of RANS equations. Because of the practical importance of such models for weather forecasting and performance prediction for engineered devices, RANS-closure modeling efforts have existed for more than a century and continue to this day. This section presents a terse development of the so-called k-ε closure model for the RANS mean-flow momentum equation (12.30) and the foundational elements of second-order closures. The many details associated with these and other RANS closure schemes, and large-eddy simulations are described in Pope (2000) and Wilcox (2006). The earlier review article by Speziale (1991) is also recommended.

The primary purpose of a turbulent-mean-flow closure model is to relate the Reynolds stress correlations, uiuj¯ to the mean-velocity field U_i. Prandtl and von Karman developed certain semi-empirical theories that attempted to provide this relationship. These theories are based on drawing an analogy between molecular-motion-based laminar momentum and scalar transport, and eddy-motion-based turbulent momentum and scalar transport. The outcome of such modeling efforts is typically an eddy viscosity ν_T (first introduced by Boussinesq in 1877) and eddy diffusivities κ_T and κ_mT for the closure-model equations:

uiuj¯=23e¯δij−νT(∂Ui∂xj+∂Uj∂xi),uiT′¯=−κT∂T¯∂xi,anduiY′¯=−κmT∂Y¯∂xi.

(12.115, 12.116, 12.117)

Equation (12.115) is mathematically analogous to the stress-rate-of-strain relationship for a Newtonian fluid (4.37) with the term that includes the turbulent kinetic energy e¯ playing the role of a turbulent pressure. It represents the turbulent viscosity hypothesis. Similarly, (12.116) and (12.117) are mathematically analogous to Fourier's law and Fick's law for molecular diffusion of heat and species, respectively, and these equations represent the gradient diffusion hypothesis for turbulent transport of heat and a passive scalar.

To illustrate the implications of such hypotheses, substitute (12.115) into (12.30) to find:

∂Ui∂t+Uj∂Ui∂xj=−1ρ∂P∂xj+∂∂xj([ν+νT](∂Ui∂xj+∂Uj∂xi)−23e¯δij)

(12.118)

for constant-density flow. The factor in [,]-brackets is commonly known as the effective viscosity, and the correspondence between this mean-flow equation and its unaveraged counterpart, (4.86) simplified for constant density, is clear and compelling. Mean-flow equations for T¯ and Y¯ similar to (12.118) are readily obtained by substituting (12.116) and (12.117) into (12.32) and (12.34), respectively. Unfortunately, the molecular-dynamics-to-eddy-dynamics analogy is imperfect. Molecular sizes are typically much less than fluid-flow gradient length scales while turbulent eddy sizes are typically comparable to fluid-flow gradient length scales. For ordinary molecule sizes, averages taken over small volumes include many molecules and these averages converge adequately for macroscopic transport predictions. Equivalent averages over eddies may be unsuccessful because turbulent eddies are so much larger than molecules. Thus, ν_T, κ_T, and κ_mT are not properties of the fluid or fluid mixture, as ν, κ, and κ_m are. Instead, ν_T, κ_T, and κ_mT are properties of the flow, and this transport-flow relationship must be modeled. Hence, (12.118) and its counterparts for T¯ and Y¯ must be regarded as approximate because (12.115) through (12.117) have inherent limitations. Nevertheless, RANS closure models involving (12.115) through (12.117) are sufficiently accurate for many tasks involving computational fluid dynamics.

From dimensional considerations, ν_T, κ_T, and κ_mT should all be proportional to the product of a characteristic turbulent length scale l_T and a characteristic turbulent velocity u_T:

νT,kT,orκmT∼lTuT.

(12.119)

For simplicity, consider fully-developed, temporally-stationary unidirectional shear flow U(y) where y is a cross-stream wall-normal coordinate (Figure 12.26). The mean-flow momentum equation in this case is:

0=−1ρdPdx+∂∂y(ν∂U∂y−uv¯)=−1ρdPdx+∂∂y([ν+νT](∂U∂y)).

(12.120)

A Mixing Length Model

An eddy viscosity for this equation can be constructed by interpreting l_T as a mixing length, defined as the cross-stream distance traveled by a fluid particle before it gives up its momentum and loses its identity. In this situation, an eddy of size l_T driven by a local shear rate of dU/dy produces a velocity fluctuation of u_T ∼ l_T (dU/dy) as it turns over, so

Figure 12.26 Schematic drawing of an eddy of size l_T in a shear flow with mean-velocity profile U(y). A velocity fluctuation, u or v, that might be produced by this eddy must be of order l_T (dU/dy). Therefore, we expect that the Reynolds shear stress will scale like uv¯∼lT2(dU/dy)2.

−uv¯=νTdUdy∼lTuTdUdy∼lT(lTdUdy)dUdy=lT2(dUdy)2.

The mixing-length concept was first introduced by Taylor (1915), but the approach was fully developed by Prandtl and his coworkers. For a wall-bounded flow it makes sense to assume that l_T is proportional to y when y = 0 defines the wall. Thus, setting l_T = κy, where κ is presumed to be constant, completes a simple mixing-length turbulence model, and (12.120) becomes:

0=−1ρdPdx+∂∂y(νdUdy+κ2y2(dUdy)2).

(12.121)

When the pressure gradient is zero or small enough to be ignored, (12.121) can be integrated once to find:

νdUdy+κ2y2(dUdy)2=const.=τwρ,

where the final equality comes from evaluating the expression on the left at y = 0. For points outside the viscous sub-layer, where the turbulence term dominates, the last equation reduces to a simple ordinary differential equation that is readily integrated to reach:

dUdy≅τwρ1κy,orUu∗≅1κlny+const.,

(12.122)

which replicates the log-law (12.88). This simplest-level turbulence model is known as an algebraic or zero-equation model. Such mixing length models can be generalized to a certain extent by using a contracted form of the mean strain-rate tensor or the mean rotation-rate tensor in place of (dU/dy)². However, there is no rational approach for relating l_T to the mean-flow field in general.

Since the development of modern computational techniques for solving partial differential equations, the need for simple intuitive approaches like the mixing length theory has essentially vanished, and Prandtl’s derivation of the empirically known logarithmic velocity distribution has only historical value. However, the relationship (12.119) remains useful for estimating the order of magnitude of the eddy diffusivity in turbulent flows, and for development of more sophisticated RANS closure models (see below). Consider the estimation task first via the specific example of thermal convection between two horizontal plates in air. If the plates are separated by a distance L = 3 m, and the lower plate is warmer by ΔT = 1°C. The equation for the vertical velocity fluctuation gives the vertical acceleration as:

Dw/Dt∼gαT′∼gΔT/T,

(12.123)

since T′ is expected to be of order ΔT and α = 1/T for a perfect gas. The time t_r to rise through a height L will be proportional to L/w, so (12.123) gives a characteristic vertical velocity acceleration of:

w/tr=w2/L∼Dw/Dt∼gΔT/T→w∼gLΔT/T≅0.3m/s.

The largest eddies will scale with the plate separation L, so the thermal eddy diffusivity, κ_T, is:

κT∼wL∼0.9m2/s,

which is significantly larger than the molecular thermal diffusivity, 2 × 10⁻⁵ m²/s.

One-Equation Models

Independently, Kolmogorov and Prandtl suggested that the velocity scale in (12.98) should be determined from the turbulent kinetic energy:

uT=ce¯,

where c is a model constant. The turbulent viscosity is then obtained from an algebraic specification of the turbulent length scale l_T, and the solution of a transport equation for e¯ that is based on its exact transport equation (12.47). In this case, the dissipation rate ε¯ and the transport terms must be modeled. For high Reynolds number turbulence, the scaling relationship (12.48) and the gradient diffusion hypothesis lead to the following model equations for the dissipation and the transport of turbulent kinetic energy:

ε¯=Cεe¯3/2/lTand−1ρ0puj¯+2νujS′ij¯−12ui2uj¯=νTσe∂e¯∂xj,

where C_ε, and σ_e are model constants. So, for constant density, the turbulent kinetic energy model equation is:

∂e¯∂t+Uj∂e¯∂xj=∂∂xj(νTσe∂e¯∂xj)−ε¯−uiuj¯∂Ui∂xj,

(12.124)

and this represents one additional nonlinear second-order partial-differential equation that must be solved, hence the name one-equation model. As mentioned in Pope, one-equation models provide a modest accuracy improvement over the simpler algebraic models.

Two-Equation Models

These models eliminate the need for a specified turbulent length scale by generating l_T from the solutions of transport equations for e¯ and ε¯. The popular k-ε closure model of Jones and Launder (1972) is described here. A k-ω closure model also exists. (Throughout much of the turbulence modeling literature “k” used for the turbulent kinetic energy, so the model name “k-ε“ is merely a specification of the dependent-field variables in the two extra partial differential equations.) The k-ε model is based on the turbulent viscosity hypothesis (12.115) with ν_T specified by (12.119), lT=e¯3/2/ε¯, and u_T = e¯^1/2:

νT=Cμ[e¯3/2/ε¯]e¯=Cμe¯2/ε¯,

(12.125)

where C_μ is one of five model constants. The first additional partial-differential equation is (12.124) for e¯. The second additional partial-differential equation is an empirical construction for the dissipation:

∂ε¯∂t+Uj∂ε¯∂xj=∂∂xj(νTσε∂ε¯∂xj)−Cε1(uiuj¯∂Ui∂xj)ε¯e¯−Cε2ε¯2e¯.

(12.126)

The standard model constants are from Launder and Sharma (1974):

Cμ=0.09,Cε1=1.44,Cε2=1.92,σe=1.0,andσε=1.3,

and these have been set so the model's predictions reasonably conform to experimentally determined mean-velocity profiles, fluctuation profiles, and energy budgets of the type shown in Figures 12.15 and 12.16 for a variety of simple turbulent flows. More recently renormalization group theory has been used to justify (12.126) with slightly modified constants (Yakhot & Orszag, 1986; Lam 1992; see also Smith & Reynolds, 1992).

When the density is constant, (12.27), (12.30), (12.115), and (12.124) through (12.126) represent a closed set of equations. Ideally, the usual viscous boundary conditions would be applied to U_i. However, steep near-wall gradients of the dependent field variables pose a significant computational challenge. Thus, boundary conditions on solid surfaces are commonly applied slightly above the surface using empirical wall functions intended to mimic the inner layer of a wall-bounded turbulent flow. Wall functions allow the mean-flow momentum equation (12.30) and the turbulence model equations, (12.124) and (12.126), to be efficiently, but approximately, evaluated near a solid surface. Unfortunately, wall functions that perform well with attached turbulent boundary layers are of questionable validity for separating, impinging, and adverse-pressure-gradient flows. Furthermore, the use of wall functions introduces an additional model parameter, the distance above the wall where boundary conditions are applied.

Overall, the k-ε turbulence model is complete and versatile. It is commonly used to rank the performance of fluid dynamic system designs before experimental tests are undertaken. Limitations on its accuracy arise from the turbulent viscosity hypothesis, the ε¯ equation, and wall functions when they are used. In addition, variations in inlet boundary conditions for e¯ and ε¯, which may not be known precisely, can produce changes in predicted results. In recent decades, two equation turbulence models based on the eddy viscosity hypothesis have begun to be eclipsed by Reynolds stress models or second-order closures that directly compute the Reynolds stress tensor from a modeled version of its exact transport equation (12.35).

Reynolds Stress Models

Reynolds stress closures for RANS equations are attractive because they eliminate the need for an eddy viscosity, but the resulting equations are considerably more complicated than those of two-equation closures. Reynolds stress models are generally superior to two-equation closures because they incorporate flow-history effects since uiuj¯ is not directly linked to ∂U_i/∂x_j as it is in (12.115); they include streamline curvature and rotation effects through their direct use of Duiuj¯/Dt; and (3) they do not require the normal stresses to be equal when ∂U_i/∂x_j vanishes, as is the case for (12.115). For Reynolds stress models, every dependent field variable in the Reynolds-averaged momentum equation (12.30) is computed, so the pallet of scalars, vectors, and tensors for the creation of closure models includes: P, U_i, and uiuj¯. Additionally, a model equation like (12.126) for the kinetic energy dissipation rate is typically solved as well, so ε¯ is also used in Reynolds stress closure models. The summary provided here is drawn largely from Speziale (1991), Pope (2000), and Wilcox (2006).

To illustrate the form of common Reynolds stress models, consider the following equivalent form of (12.35) for constant density:

∂uiuj¯∂t+Uk∂uiuj¯∂xk=−uiuk¯∂Uj∂xk−ujuk¯∂Ui∂xk−εij+Mij+Nij

(12.127)

(see Exercise 12.45). The first two terms on the right side are the Reynolds-stress production terms and do not need to be modeled. The third term on the right side (12.127) is the Reynolds-stress dissipation rate tensor:,

εij≡2ν∂ui∂xk∂uj∂xk¯≅23ε¯δij.

(12.128)

It is typically modeled as being isotropic via the approximate equality in (12.128). The fourth term on the right side of (12.127) is the Reynolds-stress transport tensor:

Mij≡∂∂xk(ν∂∂xkuiuj¯−uiujuk¯−uip¯ρδjk−ujp¯ρδik),

(12.129)

which includes viscous-, turbulent-, and pressure-transport contributions. The viscous contribution need not be modeled, and is typically negligible except near walls. The others contributions are commonly modeled by an equation representing gradient diffusion that embodies appropriate symmetries, such as:

−uiujuk¯−uip¯ρδjk−ujp¯ρδik≅23Cse¯2ε¯[∂uiuj¯∂xk+∂ujuk¯∂xi+∂ukui¯∂xj],

(12.130)

where C_s ≈ 0.11 (Mellor and Herring 1973, Launder et al. 1975). The final term on the right of (12.127) is the pressure-rate-of-strain tensor:

Nij≡pρ(∂ui∂xj+∂uj∂xi)¯.

(12.131)

It is important for the solution of (12.127) and difficult to model because simultaneous measurements of velocity-fluctuation gradients and pressure fluctuations within a turbulent flow are essentially non-existent.

Insight into the character of the pressure-rate-of-strain tensor can be obtained from the Poission equation for the pressure fluctuation p:

1ρ∂2p∂xk2=−2∂Ui∂xj∂uj∂xi−∂2∂xi∂xj(uiuj−uiuj¯),

(12.132)

(see Exercise 12.46). On the basis of this equation, pressure fluctuations are presumed to fall into three categories: (1) rapid pressure fluctuations that respond immediately to changes in mean-flow gradients via the first term on the right of (12.132), (2) slow pressure fluctuations that occur in response to changes the second term on the right of (12.132), and (3) harmonic pressure fluctuations for which ∂2p/∂xk2=0 that arise to satisfy boundary conditions on p. Overall, (12.132) indicates that pressure fluctuations may be influenced by the entire flow field; thus, an accurate local closure scheme for (12.127) is impossible. However, an assumption of locally-homogeneous turbulence allows a Green's function solution of (12.132) to be combined with (12.131) to reach:

Nij=14π∭V(∂ui∂xj+∂uj∂xi)∂2(ukul)∂yk∂yl¯d3y|x−y|+12π∂Uk∂xl∭V(∂ui∂xj+∂uj∂xi)∂ul∂yk¯d3y|x−y|.

(12.133)

where x is the position vector, y is a position-vector integration variable, and V is a volume that encompasses the entire flow field (see Exercise 12.47). Most Reynolds stress models use a form of (12.133) with modeled terms replacing the volume integrals (see Wilcox 2006).

Thus, a complete Reynolds stress closure model for constant-density turbulent flow consists of (12.27), (12.30), (12.126) or its equivalent, (12.127)–(12.129), (12.130) or its equivalent, and (12.133) or a model equation that replaces it. The performance of Reynolds-stress models is described in Pope (2000) and Wilcox (2006).