Turbulence Modeling

Seven – Turbulence Modeling Introduction Although there is a strong confidence in our ability to solve the Navier Stokes equations on modern computers, that confidence is limited to laminar – not turbulent – flows. Unfortunately, most flows of engineering interest are turbulent. Although there have been some solutions of the Navier Stokes equations directly for turbulent flows (a process known as DNS for direct numerical simulation), these computations are not practical for engineering analysis. They require extensive computer resources, are applicable only to simple geometries, and do not directly provide the average quantities of interest in engineering design. Instead of direct simulation, CFD applications to turbulent flow use models. These models range from simple algebraic models to those which require the solution of one or more partial differential equations. In general, the more complex the flow, the more complex of a turbulence model is required. Particular problems occur in modeling turbulence with large amounts of swirl as occurs in combustors and turbomachinery. These notes begin with a simple discussion of turbulence and move on to the formal statistical analyses that are used in turbulence. We will see that the formal statistical analysis leads to a closure problem that requires additional assumptions to model physical processes. The most important idea in these notes is that the choice of wall boundary conditions that you use for turbulent flows determines limits on the grid spacing. The choice to (a) model the laminar sublayer at a wall or to (b) use wall functions sets conditions on the grid size near the wall that must be followed to ensure a correct solution. What is turbulence? The transition between laminar and turbulent flow, in any geometry, is characterized by the Reynolds number, Re, based on the appropriate characteristic velocity, U, and the appropriate length parameter, L, for the flow. Re = UL [7-1] ν As usual, ν is the kinematic velocity. When the Reynolds number is below a critical level, the flow is laminar. Increases in the Reynolds number above this critical level lead to a transition region, and further increases lead to a fully turbulent region.* * For free convection flows, the borders between laminar, transition and turbulent flow are determined by the Rayleigh number, Ra = βgΔTL3 να , where β =− 1 ⎛ ∂ρ ⎞ ⎜ ⎟ is the coefficient of ρ ⎝ ∂T ⎠ P thermal expansion, which equals 1/T for an ideal gas, g is the acceleration of gravity, ΔT is a characteristic temperature difference, and α is the thermal diffusivity. Turbulence is characterized by random behavior. The value of the flow properties (velocity components, temperature, species concentrations) at any point in the flow vary in a random manner over time. This random variation may be about a constant (steady-state) average, or it may be about some overall average trend in time. Turbulence is a three-dimensional phenomenon. A fully developed laminar flow in a cylindrical pipe varies only in the radial direction. (There is no variation in the angular coordinate or down the length of the pipe after the flow becomes fully developed.) However, if such a flow were turbulent, there would be variation in the flow in all three coordinate directions. Turbulent Turbulent flows are characterized by turbulent eddies. In pictures of turbulent flows we see underlying structures that move (translate) and rotate in the mean flow. These underlying structures are called turbulent eddies. When we examine a visualization of a turbulent flow we see that these eddies have many different length scales.* The largest eddies have a characteristic velocity and length scale that are of the same order of magnitude as the characteristic velocity and Laminar length scale (U and L) for the main flow. In these large eddies (with large Reynolds numbers) the inertia effects dominate and viscous effects are negligible. [Recall that the Reynolds number is a measure of the ratio of inertia forces (ρU2) to viscous forces (μU/L).] The main flow (at some mean velocity like U) supplies energy to the large eddies. This increases their rotation rate and decreases their size. The larger eddies supply smaller eddies with energy. This creates the following energy cascade. 1. The kinetic energy from the main flow is transferred into kinetic energy of the larger eddies. 2. The kinetic energy from the larger eddies is transferred into the kinetic energy of the smaller eddies. 3. The kinetic energy of the smallest eddies is dissipated by viscous effects. In this way, the turbulence produces increased energy dissipation (and hence increased pressure loss) for the flow. The turbulence structures provide a mechanism by which energy is transformed from the main flow kinetic energy into viscous dissipation.** This mechanism is not present in laminar flows. Although turbulence increases pressure drop because of this energy transfer to viscous forces, it also produces increased mixing rates that are important in heat transfer, mixing processes, and combustion. * Although smoking is not permitted in campus buildings, all of us have probably noticed the smoke from a burning cigarette, which provides a good visualization of the development and structure of turbulent flows. This flow usually starts as a laminar flow next to the cigarette, then goes through a transition region to a turbulent flow. In the turbulent flow there are many eddies with varying characteristic sizes. These eddy structures persist for a while then break up. This is illustrated in the photograph shown here taken from the web site shown below: Humphrey_Bogart_by_Karsh_(Library_and_Archives_Canada).jpg ** Lewis F. Richardson (1881-1953), who wrote an Article entitled “Numerical Prediction of Weather Process: in 1922, wrote a poem that characterizes this process. Big whorls have little whorls That feed on their velocity And little whorls have lesser whorls And so on to viscosity. The typical lengths of the smallest turbulent eddies are about 10-4 to 10-5 meters in turbulent flows that are typical of engineering applications. A direct numerical simulation of turbulent flows would require a computational grid small enough to resolve these length scales while being able to span the entire item being designed. Even for a small item with a length of 0.01 m, 100 to 1000 grid nodes would be required in each coordinate direction. For a three-dimensional flow this means a total of 106 to 109 grid points. Such a large number of grid nodes, even for a very small scale device, makes DNS impractical for applications to typical engineering devices. Time-averaged Navier-Stokes equations (RANS) In order to analyze the random nature of turbulent flow, we write the typical flow property, φ, as the sum of a time-average property, φ̄, and a fluctuating property, φ’. We will first define what we mean by these quantities and then derive some general results for them. We will next use those results to derive the time-averaged Navier Stokes equations for turbulent flows. The resulting equations are often called the Reynolds-Averaged, Navier-Stokes (RANS) equations in honor or Osborne Reynolds who first proposed this approach. The formal definition of the division into a time average and a fluctuating quantity is shown below. ϕ = ϕ +ϕ' [7-2] The time average property, φ̄ , is defined by the following equation. ϕ= Δt 1 ϕdt Δt ∫0 [7-3] In this equation the averaging period, Δt, is assumed to be large enough to obtain a meaningful average. In transient flows (where the mean quantities are changing with time) the averaging is done by the use of an ensemble approach. However, the result is the same; there is an average property and a fluctuation due to the turbulence. In the transient flows the ensemble average is regarded as the average that would be obtained if the same transient flow were studied several times and the measurements of each study were averaged. From the definitions of the average and the fluctuating quantity we can see that the average of a fluctuating quantity is zero. 1 ϕ' = Δt Δt ∫( 0 ) 1 ϕ − ϕ dt = Δt Δt ∫ 0 1 ϕdt − Δt Δt ∫ 0 1 ϕdt = ϕ − ϕ Δt Δt ∫ dt = ϕ − ϕ = 0 [7-4] 0 In equation [7-4] we used the important result that the average of an average flow quantity is simply the original average. However, if we examine the average of two flow quantities, we do not get that result. This is shown below where φ and ψ are used to denote two different flow quantities. ψϕ = 1 Δt Δt Δt Δt Δt Δt 0 0 0 1 1 1 1 ∫ (ϕ − ϕ ')(ψ − ψ ')dt = Δt ∫ ϕψ dt − Δt ∫ ϕψ ' dt − Δt ∫ ϕψ ' dt + Δt ∫ ϕψ ' dt 0 Δt 0 Δt Δt Δt Δt 1 1 1 1 1 dt − ϕ ψ ' dt − ϕ ∫ψ ' dt + ∫ ϕ 'ψ ' dt = ϕψ − 0 − 0 + ∫ ϕ 'ψ ' dt = ϕψ ∫ ∫ Δt 0 Δt 0 Δt 0 Δt 0 Δt 0 [7-5] The overall result of equation [7-5] is usually summarized as follows. ψϕ = ϕψ + ϕ 'ψ ' [7-6] Here we see that the average of the product consists of two terms. The first is the product of the averages. The second is the average of the product of the fluctuations. Only if these fluctuations were perfectly anticorrelated would this average of the product be zero. It is this basic result that provides an introduction to turbulence quantities in the time-averaged Navier-Stokes equations. Equation [7-6] tells us that the average of the product of two fluctuating quantities, which may be the same, will typically not vanish in a turbulent flow. We can use this result to define the root mean square average as our typical measure of the magnitude of a turbulent fluctuation. We define this quantity as follows. Δt 1 (ϕ ' ) 2 dt ∫ t Δ 0 ϕ rms = (ϕ ' ) 2 = [7-7] The velocity fluctuations in a turbulent flow, using Cartesian coordinates and the usual velocity components, are u’, v’, and w’. The kinetic energy of the turbulence, k, is given by the sum of the squares of these velocity fluctuations. k= [ 1 (u ' ) 2 + (v' ) 2 + ( w' ) 2 2 ] [7-8] We see that this is a typical fluid mechanics “kinetic energy” term, which is actually a kinetic energy per unit mass. The time average of the spatial derivative of a flow quantity, φ, is found as follows. Δt ∂ ∂ϕ 1 ∂ϕ dt = = ∫ ∂xi Δt 0 ∂xi ∂xi ⎡ 1 Δt ⎤ ∂ϕ ⎢ ∫ ϕdt ⎥ = ⎣ Δt 0 ⎦ ∂xi [7-9] Thus the average of the derivative (with respect to the coordinate xi) is simply the derivative of the average. This result applies to derivatives of any order. We are now ready to derive the time-average Navier-Stokes equations. We start with equation [1-81]; this general transport equation is copied below. ⎡ ∂ρφ ∂ρ ui ϕ ⎤ + c⎢ ⎥= ∂xi ⎦ ⎣ ∂t ∂ (ϕ ) ∂ϕ Γ + S (ϕ ) ∂xi ∂xi [1-81] To simplify the derivation, we will assume that we have constant properties and a zero source term. Further assume that we have a steady problem so that the time derivative is zero. We will divide the equation by the product, ρc, and write our initial equation as follows. ∂ui ϕ ∂ ∂ϕ = γ (ϕ ) ∂xi ∂xi ∂xi [7-10] In this equation we have defined γ(φ) as follows γ (ϕ ) Γ (ϕ ) = ρc [7-11] Recall that in the general equation c = 1 unless φ is the temperature in the energy equation. In most cases, then, the denominator in equation [7-11] is simply the density, ρ. If φ represents the temperature in the energy equation then c = cp or cv, depending on the equation. In the momentum equation, γ(φ) would be the kinematic viscosity, ν = μ/ρ; in the energy equation with temperature, it would be the thermal diffusivity, α = k/ρcp. (Note that for the momentum equation, Γ(φ) is the dynamic viscosity, μ, and for the energy equation, with temperature as the variable, Γ(φ) is the thermal conductivity.) We want to obtain the time average of the entire transport equation in [7-10]. We do this by substituting both sides of this equation into the definition of the time average from equation [7-3]. Δt Δt 1 ⎡ ∂ui ϕ ⎤ 1 ⎡ (ϕ ) ∂ ∂ϕ ⎤ ⎢ ⎥ dt = ⎢γ ⎥ dt Δt ∫0 ⎣ ∂xi ⎦ Δt ∫0 ⎣ ∂xi ∂xi ⎦ [7-12] Both sides of equation [7-12] are averages of spatial derivatives. We can use the result of equation [7-9] that the average of (spatial) derivatives are simply the derivatives of averages to rewrite equation [7-12] as follows. ∂ ui ϕ ∂ ∂ϕ = γ (ϕ ) ∂xi ∂xi ∂xi [7-13] We can use the result of equation [7-6] for the average of a product to rewrite equation [7-13] as follows. ( ) ∂ ui ϕ + ui 'ϕ ' ∂ ∂ϕ = γ (ϕ ) ∂xi ∂xi ∂xi [7-14] If we compare equations [7-10] and [7-14], we see that the time averaged equation in [7-14] has two differences. The first is that the average quantities appear in [7-14] in place of the instantaneous quantities in [7-10]. The second change is the addition of the correlation term involving the average of two fluctuating quantities. This is the mathematical term that leads to problems in the computations of turbulent flows using the time-averaged Navier Stokes equations. When we use the general balance equation to represent the continuity equation, we take φ = 1 and Γ = 0. In this case the fluctuation term vanishes and we have the constant-property, steadystate equation for turbulent flows as follows. ∂ui = 0 ∂xi [7-15] In principle, we could take the time average of equation [7-14] to get a differential equation for the correlation terms. However, this would introduce new, higher order correlations. This creates what is known as the closure problem. At some point we have to give up using the time averages and try to use models for the correlation terms. This is the next topic of these notes. Introduction to models of turbulent flows The general approach for modeling the fluctuation terms is based on the Boussinesq assumption that this additional term can be modeled in the same was as a viscous or diffusive flux term. That is we assume that we can use an equation like the following. ui 'ϕ ' = − γ t(ϕ ) ∂ϕ ∂xi [7-16] Here, the “t” subscript on the generalized transport coefficient, γ, indicates that we are defining a transport coefficient for turbulent flow. With this definition, we can rewrite equation [7-14] as follows. ( ) (γ ∂ ui ϕ = ∂xi (ϕ ) + γ t(ϕ ) ) ∂ ∂ϕ ∂xi ∂xi [7-17] With this formulation, the basic balance equations for turbulent flow have the same form as the equations for laminar flow. Thus we can use the same computational algorithms for solving the turbulent flows. We have moved the problem from one of calculating the flows to one of calculating the turbulent transport coefficient. In cases where φ is a velocity component, uj, we use a slightly different form of equation [7-16], that is based on the general equations for the viscous stresses and the momentum equations in equations [3-4] and [3-5]. Those equations are copied below. ⎡ ∂ ui τ ij = μ ⎢ ⎢⎣ ∂x j ∂ρu j ∂t + ∂ρ uiu j ∂xi =− + ∂uj ⎤ 2 ⎥ + (κ − μ )Δδ ij ∂xi ⎥⎦ 3 ∂p ∂τ ij + + ρB j ∂x j ∂xi [3-4] j = 1,K3 [3-5] We want to substitute equation [3-4] into equation [3-5] and apply the steady-state, constantproperty assumptions that we are using here. The steady-state assumption allows us to drop the transient term. Our assumption of constant properties, allows us to bring properties outside of derivative operators. For constant density, the dilatation, Δ, is zero. ∂ ui u j ∂xi =− 1 ∂p μ ∂ ⎛⎜ ∂ ui ∂ u j ⎞⎟ + + + Bj ρ ∂x j ρ ∂xi ⎜⎝ ∂x j ∂xi ⎟⎠ j = 1,K3 [7-18] Although flux terms for heat conduction and diffusion involve only a single gradient, we see that the viscous stress terms involve two separate gradients. Thus, we define the general velocity component correlation term as follows. ui ' u j ' = − ⎛ ∂u ∂uj ⎞ 2 μt ⎛⎜ ∂ ui ∂ u j ⎞⎟ 2 ⎟+ δ k + + δ ij k = −ν t ⎜ i + ⎜ ∂x j ∂xi ⎟ 3 ij ρ ⎜⎝ ∂x j ∂xi ⎟⎠ 3 ⎝ ⎠ i = 1,K3; j = 1,K3 [7-19] Here we have defined μt and νt as the turbulent dynamic and kinematic viscosities. Again, we are hiding our ignorance of turbulent flows in these terms. We are left with equations that are similar to the laminar flow equations, but we have to find a way to predict μt and νt. We see that the term − ρ ui ' u j ' follows the same form of the equation as the viscous stress term, τij. These terms are called the Reynolds stress terms. We see that they are symmetric, so there are only six different Reynolds stress terms. We will show below that advanced turbulence models seek to predict each Reynolds stress. The mixing-length model of simple turbulent flows Considerations of turbulent flows started with Reynolds’s work in the late 1800s and the contemporary assumption by Boussinesq that the Reynolds stresses could be related to the gradients of the mean flow velocity components. Subsequently in the early 1900s Prandtl developed the model of the boundary layer in fluid flows. Von Karman later expanded this work. This work developed Prandtl’s mixing-length theory as a useful approach for modeling simple turbulent flows. The mixing-length model remains a useful model for simple two-dimensional flows in jets and near-wall regions. In addition, this basic analysis is used to develop the wall functions that are used as boundary conditions in CFD turbulent flow models. In the mixing-length model, there is a viscous (laminar) sublayer close to the wall. Beyond this sublayer, there is a transition region into a fully turbulent boundary layer. The turbulent boundary layer near the wall is characterized by the following quantities in a two dimensional flow where the x direction is the predominant flow direction, with flow velocity u, and the y direction is the direction perpendicular to the flow. For flows near walls, the point y = 0 is the wall. uτ = τw ρ u+ = u uτ y+ = y ρuτ y = yτ μ [7-20] Here, τw is the wall shear stress and the quantity uτ is called the friction velocity. The law of the wall postulates that there is a relationship between the dimensionless velocity and distance variables defined in equation [7-20]. u+ = f ( y+ ) [7-21] For flows near a wall, this relationship has different forms. The region closest to the wall (y+ < 5) is the laminar sublayer in which the normal equations for viscous stresses are valid. Beyond this region the values of u+ are given by the following equations. ⎧ ⎪ y+ ⎪ ln( y + ) ln( Ey + ) ⎪ u+ = ⎨ +B= κ ⎪ κ 1 y ⎛ ⎞ ⎪ u + = ln⎜ ⎟ − A ⎪⎩ max κ ⎝ δ ⎠ 5 < y + < 30 30 < y + < 500 y + > 500 [7-22] In equation [7-22], we have introduced the constants κ, known as the von Karman constant, which has a value of 0.4, B = 5.5, E = 9.8, and A. The values of these constants are for smooth walls. We have also introduced the boundary-layer thickness, δ, for the velocity profile away from the boundary layer. The results above have not made explicit use of the mixing-length; instead, it has shown the final results. The mixing-length model is a way of determining a characteristic length, l, for the turbulence. If we take a simple dimensional analysis of the turbulent kinematic viscosity, νt, with dimensions of length squared over time, we see that this variable has the same dimensions as the product of a velocity and a length. Thus we can write the turbulent kinematic viscosity as a product of a velocity and a length. If we take the velocity as U and the length scale as l, we can write: νt = Cμ U l [7-23] where = Cμ is an empirical constant. The k-ε model of turbulent flow The most popular model for modeling turbulence in common flows of engineering interest is called the k-ε model. This model uses the kinetic energy of turbulence, defined in equation [7-8] and the turbulence dissipation rate, ε. This dissipation rate measures the rate of energy transfer from kinetic energy in the smallest eddies when they do work against the viscous forces. Formally the dissipation is defined in terms of the deformation (or strain) rates, eij which are defined below. eij = ∂ ui ∂ u j + ∂x j ∂xi so that τ ij = μeij + (κ − 2 μ )Δδ ij 3 [7-24] The formal definition of the dissipation is shown in equation [7-25]. The dimensions of dissipation are energy divided by (mass times time). (As noted in the paragraph following equation [7-8], we conventionally define the energy per unit mass as the “kinetic energy”; that is being done here.) The SI units for ε are m2/s3. ε = 2ν e'ij e'ij [7-25] We are using the summation convention here. The implied summation over the two repeated indices, i and j, gives nine terms in this equation. Also, ν in this equation is the actually kinematic viscosity, not the turbulent kinematic viscosity. The basic approach of equation [7-23] in which the turbulent velocity is treated as the product of a length scale times a velocity scale is also used here. In the k-ε model the velocity scale is the square root of the kinetic energy of turbulence; i.e., U = k1/2. The length scale is equal to k3/2/ε. Substituting these terms into equation [7-23] gives the following equation for the turbulent viscosity. ν t = Cμ k 1 2 k 3 ε 2 = Cμ k2 ε μt = Cμ ρ k2 ε [7-26] The values of k and ε used to compute the turbulent viscosity are found by solving partial differential equations that have the same form as the usual balance equations. These equations are described below. A balance equation for the turbulent kinetic energy can be derived by combining the time averaged Navier-Stokes for all components after multiplying each by the fluctuating velocity components. This is similar to the approach used to obtain equation [1-45] for the kinetic energy of the main flow. The resulting balance has several terms that cannot be evaluated and must be modeled. This balance equation is shown below and the various terms in the equation are discussed following the equation. The summation convention is used in this equation. In some cases a double summation, over repeated i and j indices is implied. ∂u ∂ρk ∂ρ ui k ∂ ∂k ∂ ⎛ρ ⎞ + = μ − ⎜ ui ' u j ' u j ' + p' ui ' ⎟ + ρ ui ' u j ' i − ρε ∂x j ∂t ∂xi ∂xi ∂xi ∂xi ⎝ 2 ⎠ [7-27] The two terms on the left hand side represent the usual transient and convection terms. Similarly, the first term on the right hand side represents the diffusive transport of turbulent kinetic energy. This term is similar to the usual gradient term that we have in our general balance equation. The remaining terms in equation [7-27] are regarded as source terms in the general balance equation. However, these terms are not found directly. Instead, they must be modeled. We will consider each of these terms below. The term ∂ ⎛ρ ⎞ ⎜ ui ' u j ' u j ' + p ' ui ' ⎟ represents the turbulent diffusion of kinetic energy. This is ∂xi ⎝ 2 ⎠ the transport of turbulent kinetic energy by the fluctuations themselves. It is modeled as a gradient term. ⎛ρ ⎞ μ ∂k − ⎜ ui ' u j ' u j ' + p ' ui ' ⎟ ≈ t ⎝2 ⎠ σ k ∂xi [7-28] The quantity σk is called the turbulent Prandtl number for kinetic energy, which is an empirical constant for turbulent flows. See the discussion of the turbulent Prandtl number on page 7.12. The term (really nine similar terms with the summation convention) ρ ui ' u j ' ∂ui is the product of ∂x j a Reynolds stress times a gradient of the mean flow. This term represents the production of turbulent kinetic energy by transfer of kinetic energy from the mean flow, and is usually given the symbol Pk. If we use equation [7-19] for the Reynolds stress terms, we can write this term as follows. Pk = ρ ui ' u j ' ⎛∂u ∂uj ∂ui = −μt ⎜ i + ⎜ ∂x j ∂x j ∂xi ⎝ ⎞ ∂ui ⎟ ⎟ ∂x j ⎠ [7-29] If we substitute equations [7-28] and [7-29] into equation [7-27], we obtain the following result. ∂ρk ∂ρ ui k ∂ ∂k ∂ μ t ∂k + = μ + + Pk − ρε ∂t ∂xi ∂xi ∂xi ∂xi σ k ∂xi [7-30] The first two terms on the right-hand side can be combined to give the equation for the kinetic energy of turbulence. μ ∂ρk ∂ρ ui k ∂ ⎛ ⎜⎜ μ + t + = ∂t ∂xi σk ∂xi ⎝ ⎞ ∂k ⎟⎟ + Pk − ρε ⎠ ∂xi [7-31] The laminar viscosity in the diffusion term is much smaller than the turbulent viscosity (divided by the Prandtl number for k) and is usually not included in the computations. Note that this equation has the same form as the general CFD equation (transient – convection – diffusion – source) so we perform the numerical analysis of this equation by the usual CFD algorithms for other variables. Although we defined the turbulent dissipation, ε, in equation [7-25], we do not use this equation to compute ε. Instead, we have to solve another partial differential equation to compute this quantity. The final equation differential equation for the dissipation, with all the modeling assumptions already made, is shown below. ∂ρε ∂ρ ui ε ∂ μ t ∂ε ε ε2 + = + Cε 1 Pk − ρCε 2 k k ∂t ∂xi ∂xi σ ε ∂xi [7-32] The first three terms in this equation are the usual terms in the general balance equation, representing the transient term, the convection term and the diffusion term. For the turbulent flow, we use the mean velocities in the convection term and the turbulent transport coefficient in the diffusion term. The effective transport coefficient is expressed as the ratio of the viscosity to the Prandtl number for dissipation, σε. The last two terms in equation [7-32], with the empirical constants Cε1 and Cε2, represent the production and destruction of dissipation. This equation is very similar to the corresponding equation [7-31] for the kinetic energy of turbulence. The various empirical constants used in the standard k-ε model are summarized below. C μ = 0.09 Cε 1 = 1.44 Cε 2 = 1.92 σ k = 1.0 σ ε = 1.3 [7-33] These notes have described the original k-ε model which still has wide use in engineering applications. Other versions of this model known as the renormalizable group (RNG) and the realizable k-ε models use alternative derivations of the model to be more consistent with the basic physics of turbulence. These models have different differential equations for the production and dissipation terms and different equations for computing the turbulent viscosity from k and ε. The empirical constants used in these models are also different. These models are recommended in flows that high strain rates (i.e. large velocity gradients). Boundary conditions and wall functions In applying turbulence models like the k-ε model it is necessary to specify boundary conditions for the dependent variables, such as k and ε, at a variety of physical boundaries encountered in the flow. The most complex application of boundary conditions is at solid walls. At large Reynolds numbers, the laminar sublayer is so small that it is not possible to use a grid that is fine enough to resolve this region. As a consequence, the boundary conditions at solid walls are usually handled by the mixing length relations from equation [7-22]. The starting point for the derivation of wall functions is the part of equation [7-22] which is valid for y+ values between 30 and 500. Here we use the second form of that equation. u+ = u ln( Ey + ) 1 ⎛ Eyuτ ⎞ = = ln⎜ ⎟ κ κ ⎝ ν ⎠ uτ [7-34] In the general consideration of wall functions, u+ is regarded as a normalized velocity component parallel to the wall and y+ is regarded as a normalized distance perpendicular to the wall. To use this equation we need a value for the wall shear stress, τw, which is used in the normalizing variable, uτ = τw / ρ . (Equation [7-20] shows how this parameter is used in the definitions of u+ and y+.) This parameter is usually calculated from the equilibrium assumption that the production and dissipation of turbulence are equal. In this case one can derive the following relationship between uτ and the turbulent kinetic energy. 1 uτ = Cμ 4 k [7-35] We can combine equations [7-34] and [7-35] to get an expression for the velocity at the first node in from the wall. ⎛ EyC 14 k ⎞ 1 ⎛ Eyuτ ⎞ μ ⎜ ⎟ u = uτ ln⎜ ⎟ = Cμ 4 k ln⎜ ⎟ ν ν ⎠ ⎝ ⎝ ⎠ [7-36] In this equation Cμ, and E are known dimensionless constants,1 ν is the (laminar) kinematic viscosity, and k and y are the values, respectively, of turbulent kinetic energy and distance from the wall at the first node. Other wall functions are used for k, ε, and other transported variables such as temperature and species concentrations. At outlets, the gradients of turbulent kinetic energy and dissipation are assumed to be zero in the direction of the flow. At inlets it is necessary to specify profiles of k and ε. Ideally these should come from measurements on flows similar to the one being modeled. If such data are not available, the inlet values of k and ε can be estimated form the following equations. 2 k = Cuinlet l = 0.07 L 3 ε = Cμ 4 k l [7-37] Here L is a characteristic length, which may be taken as the hydraulic diameter of an inlet. (Recall that the hydraulic diameter is 4A/P, where A is the area and P is the perimeter. For circular flow passages this is the diameter.) The constant C is a small number in the range of 0.01 to 0.05. Another way to specify k is to specify the turbulence intensity. This is defined as the ratio of the kinetic energy of turbulence to the mean velocity squared. This is equivalent to allowing the user to specify C in the initialization equation for k above The wall functions proposed above are limited to “large” Reynolds numbers. Separate wall functions are required for lower Reynolds numbers. The renormalizable group (RNG) and the realizable k-ε models, mentioned above as alternative k-ε models have different initial condition specifications and different wall functions. 1 E = 9.8 and for the k-ε model, Cμ = 0.09. Wall boundary conditions and grid sizing in CFD codes Perhaps the most important idea of these notes is the relationship between the grid size and the turbulence model. As noted in the previous section there are two possible approaches to wall boundary conditions: (1) use a find grid that computes the laminar sublayer (y+ < 5) or (2) use wall functions based on boundary layer theory of turbulent flow. The choice of which approach to use is coupled to the choice of grid size. • If you want to model the laminar sublayer, it is important to have some number of grid nodes (about ten or more) in the region for which y+ < 5. • If you plan to use wall functions, the first grid node away from the wall must have a y+ value greater than 30 (and less than 500). You will generally not know what the y+ values will be until the computations are completed. However, you must check these values in the completed computations to ensure that they are in the proper range throughout the flow. Fluent provides an option to produce plots of y+ values that you can check. Special wall functions are used for the energy equation in cases where there is a high Mach number that leads to significant heating by viscous dissipation. The usual wall-function approach is based on the assumption of equilibrium turbulence, viz. the production of turbulence equals the dissipation of turbulence. This is a good assumption for high Reynolds numbers. However in regions of low (but turbulent) Reynolds numbers this assumption is not valid and modified wall functions (sometimes called enhanced or non-equilibrium wall functions) are used in this case. Chapter 12 in the Fluent Users’ Manual has complete discussion of the various turbulence models and the different approaches to handling the wall boundary conditions. Other transport equations: the turbulent “Prandtl” number All the discussion to this point has considered a model for the turbulent viscosity that is used in the momentum equations. How do we compute the effect of turbulence on the energy equation or the species balance equation? This is done by making an empirical correction to the turbulent viscosity to obtain the appropriate turbulent diffusion coefficient for the other properties. The empirical parameter used is called the turbulent Prandtl number. In heat transfer the Prandtl number is a property of the fluid defined by the equation Pr = μcp/k. We rewrite the definition of the kinematic viscosity, ν = μ/ρ, as μ = ρν; similarly, we can rewrite the definition of the thermal diffusivity, α = k/ρcp as k = ρcpα. Substituting these results into the definition of Prandtl number gives an alternative expression for the Prandtl number as Pr = μcp/k = (ρν)μcp/(ρcpα) = ν/α. Thus, if we know the laminar kinematic viscosity and Prandtl number for a fluid, we can find the thermal diffusivity for the fluid from the equation α = ν/Pr. This is the basic kind of relationship used in turbulent flows. If we have an empirical model for the turbulent Prandtl number, Prt, we can find the turbulent thermal diffusivity, αt, from the turbulent kinematic viscosity by the equation αt = νt,/Prt. Recall the basic equation [7-17] for turbulent transport copied below. ( ) (γ ∂ ui ϕ = ∂xi (ϕ ) + γ t(ϕ ) ) ∂∂x i ∂ϕ ∂xi [7-17] Recall that we defined the general coefficient γ(φ) = Γ(φ)/(ρc) in equation [7-11]. For momentum transport (where c = 1), γ is the kinematic viscosity, ν, and Γ is the dynamic viscosity, μ. Although the Prandtl number is formally defined as μcp/k = ν/α, in turbulence models, the turbulent “Prandtl” number for a particular flow variable, φ, is defined as the ratio ν/γ(φ) and is given the symbol σφ. (Here I put quotation marks around Prandtl to indicate that the name Prandtl number is applied to ratios other than the original definition of the Prandtl number as ν/α. With this definition of the general turbulent Prandtl number we can rewrite equation [7-17] as follows. ( ) ⎛ ν ⎞ ∂ ∂ϕ ∂ ui ϕ ν = ⎜ + turb ⎟ ⎜ σϕ,lam σϕ,turb ⎟ ∂xi ∂xi ∂xi ⎝ ⎠ [7-38] As noted above, the turbulent transport property is usually much larger than the laminar transport property and is often ignored in computations. The laminar “Prandtl” number is a fluid property that can be found from data supplied with the CFD codes or entered as input by the user from property tables. The turbulent “Prandtl” number is usually set as a default value in the CFD code. These default values are usually satisfactory for most applications. Users can override these default values if they have other data available to them to justify this. Other turbulence models Several other models that solve two partial differential equations have been proposed, and each has their proponents. However none of these seems to have received better acceptance than the k-ε model. For complex flows, typically those with large amounts of swirl, the Reynolds Stress model is used. This model solves six partial differential equations, one for each independent Reynolds stress term. In addition, the dissipation equation [7-32] must be solved to obtain terms used in the differential equations for the Reynolds stresses. The use of the Reynolds stress model has not been completely successful. For some flows this model gives much better agreement with experimental data as compared to the k-ε model. For other flows, little improvement is obtained by going from the k-ε model to the Reynolds stress model. The general equation for the Reynolds stress model is shown below. ∂ ui ' u j ' ∂t + ∂ul ui ' u j ' ∂u ∂ μt ∂ ui ' u j ' ⎛ ∂u ⎞ 2εδ ij = − ⎜⎜ ui ' ul ' j + ul ' u j ' i ⎟⎟ − + Π ij + Ωij [7-39] 3 ∂xl ∂xl σ k ∂xl ∂xl ∂xl ⎠ ⎝ Transient Convection Diffusion Production Dissipation Pressure Strain Rotation The first five terms appearing in full in the above equation, reading from left to right, are the usual transient, convection, diffusion, production and dissipation terms. The remaining terms, the pressure-strain term, Πij, and the rotational term Ωij are described below. The δij term is the Kronecker delta. All other terms in the above equation have been defined previously. Note that we have a repeated sum over the repeated indexl. The i and j indices represent a particular Reynolds stress and we have to solve equation [7-39] for the six unique Reynolds stresses. The pressure strain interaction terms are important ones in this equation, but they are also the most difficult ones to model accurately. They come about from pressure fluctuations that occur when two eddies interact with each other and when an eddy in one flow region interacts with the main flow in another region. The simplest form of this term is ∂u ⎞ 2εδ ij ⎛ 2 ε⎛ ∂u 2 ⎞ Π ij = −C1 ⎜ ui ' u j ' − kδ ij ⎟ + C2 ⎜⎜ ui ' ul ' j + ul ' u j ' i − Pkδ ij ⎟⎟ − k⎝ 3 ∂xl 3 ∂xl 3 ⎠ ⎠ ⎝ [7-40] In this equation, C1 = 1.8, C2 = 0.6, and the Pk term is the total production of kinetic energy defined in equation [7-29]. The rotational term uses the local rotation vector, ω, whose components are ωk. It also uses the symbol εijk which is zero if any two indices are the same and +1 if i = 1, j = 2, and k = 3. For all other combinations of ijk, εijk = -1 if the indices ijk are an odd permutation of ijk or +1 if the indices are an even permutation of ijk. (In an odd or even permutation an odd or even number of moves are required to get 123 into the given set of indices. For example 321 is an odd permutation since it only takes one move, swapping the 3 with the 1. The order 231 is an even permutation since we need two moves. We first swap the 1 and the 3 (giving 321) then we swap the 2 and the 3 to get the 231 result. Other orders of swaps are possible, but it will always take an even number of swaps to permute 123 into 231. With these definitions, the rotation term is ( Ωij = −2ω k u j ' ul 'ε ikl + ui ' ul 'ε jkl ) [7-41] This term has an implied summation over the two repeated indices k and l. In the Reynolds stress model, the kinetic energy of turbulence is found from the sum of the three Reynolds stresses with the same index. k= u j ' u j ' (u1 ' )2 + (u2 ' )2 + (u3 ' )2 = 2 2 [7-42] Because k can be found from this equation, it is not necessary to solve a separate equation for k. However, it is still necessary to solve an equation for the dissipation, ε. Thus, the Reynolds stress model requires the solution of seven partial differential equations to compute the turbulence properties. An approximate form of the Reynolds stress model, known as the algebraic Reynolds stress model, has been used. In this model the source terms for production and destruction of the Reynolds stress terms are assumed to be equal in magnitude. This means that it is not necessary to solve differential equations. Instead, a system of six simultaneous linear equations is solved at each node for the Reynolds stresses at that node. It is still necessary to solve a differential equation for the kinetic energy and the dissipation. A new model called the v2f (or v2-f) model has been used to allow improved results (as compared to the k-ε model) for low Reynolds number turbulent flows. The Spalart-Allmaras model is a one-equation model that solves a transport equation for the turbulent viscosity. It was designed specifically for aerospace applications involving flows along a surface. The detached eddy simulation (DES) model is a modified version of the SpalartAllmaras model that is designed to be a less expensive alternative to LES. These notes have not discussed turbulence models for compressible flows. In such models it is necessary to consider density fluctuations. One approach for doing this is the use of Favre averages in which the instantaneous flow quantity is split into two parts as follows. Δt 1 ρϕdt Δt ∫0 ~ ϕ = ϕ + ϕ"= + ϕ" Δt 1 ρdt Δt ∫0 Models of reacting flows sometimes use a differential equation for the evolution of the probability distribution function (pdf) of the flows. These models have not been extended to general turbulence modeling. The concept of direct numerical simulation was mentioned in the introduction to these notes. A related concept is the use of large eddy simulation (LES). This is a case where the computational grid is fine enough to resolve the large scale structures in the turbulent flow. A model is still required to characterize the smaller scales. These models are called sub-grid-scale models. At present, LES does not provide a significant advantage for routine computations of flows in simple engineering design applications. However, it is the best way to get good results for complex transient turbulent flows. Conclusions on turbulence models Because of the widespread use of turbulent flows in practical systems, it is imperative that a method be developed for computation of turbulence in computational fluid dynamics. Unfortunately, there is no one model that will provide assurance that the properties of the turbulent flows will be computed accurately. The k-ε model and Reynolds stress models are the ones most commonly used in commercial codes and practical applications. The k-ε model is usually adequate for simple turbulent flows in which there is not a large degree of swirl. For swirling flows, the Reynolds stress model usually produces better results than the k-ε model, but this improved performance is not ensured. The Spalart-Allmaras model is the model of choice for external flows, both in aerospace applications and for ground vehicles. An alternative model to the k-ε model is the k-ω model where ω = ε/k; this model is used mainly for free shear flows. Advanced models such as LES have the potential for improving the accuracy of turbulent CFD results, but are presently too costly to be used for routine engineering design problems. DNS is still regarded mainly as a research tool. We have not discussed the overall solution procedures. The turbulence equations may require a finer time and grid scale than the equations for other flow quantities. We also have to be careful to maintain the linkages in calculating k and ε, which are used to compute μt. Once μt.is computed, it is used to compute the mean velocities. The mean velocities are used in the equations for k and ε, and so on. Guidance on turbulence models The choice of the turbulence model to use and the associated choice of how to handle the boundary conditions is a complex question. The Fluent user’s guide states that “no single turbulence model is universally accepted as being superior for all classes of problems.” The simplest approach for economical solutions of internal flow problems is the use of the k-ε model with wall functions. (Enhanced wall functions should be used if the Reynolds number is low.) For external flows the Spalart-Allmaras model should be used. This model may be used with wall functions or resolution of the viscous sublayer. Chapter 12 of the Fluent Users’ manual gives a full discussion of turbulence models and guidance for when to apply a particular model and the types of boundary conditions that can be used with a particular model.