Turbulence Modelling & Simulation

Overview

This document is meant to introduce the basic routes employed within TransAT for the modelling and simulation of turbulent fluid flow and scalar convection, without exploring the foundations and model derivations. We briefly present the transport equations and models for RANS and Scale-resolving strategies, including LES and V-LES. The document presents in addition a thorough description of near wall treatment on both approaches. Only the models that have been validated are presented in this document. The extension of these models for multiphase turbulent flows is not treated here. Examples of applications of the models and strategies can be found in internal reports available on the webpage.

Turbulence Modelling & Simulation

Contents
1 Statistical Turbulence Modelling (RANS)
1.1 Reynolds Averaging Procedure
1.2 The Reynolds Averaged Navier-Stokes Equations
1.3 The Reynolds Averaged Energy Equation
1.4 The Reynolds Averaged Scalar Equation
1.5 The Reynolds Stress Equations
1.6 The Dissipation Term
1.7 The Turbulent Kinetic Energy
1.8 Chemical species equations
2 Two-Equation Turbulence Models
2.1 The Standard k - ε Turbulence Model
2.2 The RNG k - ε Turbulence Model
2.3 Extensions to k - ε Turbulence Model
  2.3.1 Yap correction
  2.3.2 Swirl correction
  2.3.3 Compressible correction
3 Near–Wall Modelling
3.1 Near–Wall Treatment in High–Re Flows
  3.1.1 Standard Wall–Function Approach
  3.1.2 Two–layer turbulence models
   k^1∕2 velocity scale based model Rodi (1991): TLK
   (v′²)^1∕2 velocity scale based model Rodi & Mansour (1993): TLV
3.2 The Low-Re k - ε Turbulence Model
  3.2.1 The Jones & Launder and Launder & Sharma Models
  3.2.2 The Abe-Kondoh-Nagano Model
  3.2.3 The Lam-Bremhorst Model
4 Turbulent Heat flux Modelling
4.1 Turbulent Heat Flux Model (THFM)
4.2 Full Algebraic Turbulent Heat Flux Model: ____θ′² - ε_θ′
  4.2.1 Equivalence between THFM and ATHFM heat flux modelling
4.3 Algebraic Turbulent Heat Flux Model: ____θ′²
4.4 WET Model
4.5 Generalized Gradient Diffusion Hypothesis (GGDH)
4.6 Simple Gradient Diffusion Hypothesis (SGDH)
4.7 Variable Turbulent Prandtl Number Model
5 Buoyancy–Driven Turbulent Flows
5.1 The Equations of Motion
6 Scale Resolving Strategies: LES and V-LES
6.1 The filtered Navier-Stokes equations (bases of LES)
6.2 Sub-grid Scale (SGS) Modelling
  6.2.1 The Smagorinsky Model
  6.2.2 The Dynamic Approach (DSM)
  6.2.3 Turbulent Prandtl Number
  6.2.4 The WALE Model
6.3 Variable Turbulence Prandtl Number
6.4 Wall treatment in LES
  6.4.1 The Werner-Wengle model
6.5 Very-Large Eddy Simulation Concept (V-LES)
References

Chapter 1
Statistical Turbulence Modelling (RANS)

In this Chapter, the Reynolds-averaged equations for all quantities used in TransAT are presented. The initial equations, on which Reynolds averaging is applied, are presented in the Equations and Algorithms manual, Chapter Incompressible Navier-Stokes Equations .

1.1 Reynolds Averaging Procedure

According to the Reynolds flow-decomposition procedure, a flow variable ϕ can be expressed as the sum of an average and a fluctuating quantity:

-- ϕ = ϕ + ϕ′. ϕ = vi, p, T ,...; (1.1) v = v-+ u′; p = p+ p′; T = T-+ θ′ i i i

where ϕ is further taken to be the ensemble average of ϕ and ϕ′ is the fluctuation around this mean (ϕ′ = 0). Before proceeding further, it is necessary to recall the statistical rules of Reynolds averaging, often referred to as Reynolds conditions, which any two flow variables ϕ and ψ must obey:

------ -- -- --- -- ψ + ϕ = ψ +-ϕ; -aψ-= a ψ--(a = const.), (1.2) ∂ϕ ∂ϕ ∂ϕ ∂ϕ ---- ---- ∂t-= ∂t; ∂x--= ∂x--; ϕψ = ψ ϕ, j j

and

ϕ′ = ψ-′ = 0, (1.3) ---- ---- -′--′ ϕψ-=-ϕ ψ-+-ϕ ψ , ϕ′ψ = ψ′ϕ = 0.

1.2 The Reynolds Averaged Navier-Stokes Equations

For constant density and viscosity, Reynolds averaging, applied to the system of instantaneous equations of conservation yields the Reynolds averaged mass and momentum conservation equations, known as the Reynolds Averaged Navier-Stokes Equations (RANS):

-- ∂-vi = 0 (1.4) ∂xi- -- D-vi = - 1-∂p-+ ν ∇2v- - ∂τij+ ρg ; i = 1,2,3 (1.5) Dt ρ∂xi i ∂xj i

where τ_ij ≡u′_iu′_j stands for the Reynolds-stress tensor (which is symmetric, i.e. τ_ij = τ_ji ) representing the contribution of the turbulent motion generated by velocity fluctuations to the mean stress tensor in momentum equations. This turbulent motion will in turn result in an increase of momentum exchange and mixing. The off-diagonal components of τ_ij , or shear stresses (u′v′,v′w′, and u′w′), prevail in theory in the transport of mean momentum by turbulent motion, while the diagonal terms, or normal stresses (u′² ,v′² , and w′² ), only play a minor role.

1.3 The Reynolds Averaged Energy Equation

The equation of energy conservation can be averaged with respect to time and solved together with the equations describing the mean flow and turbulence. If this equation is written for the temperature, the Reynolds-averaging procedure applied to T yields

-- ( -- ) ρ c D-T- = -∂-- λ∂-T-- ρ c u′θ′ (1.6) p Dt ∂xj ∂xj p j ◟-------◝◜′′-------◞ - qtotal

where q′′_total = q′′_j + Q′′_j is the total rate of heat transfer due to both molecular and turbulent motions, and Q′′_j = -ρc_p u′_jθ′ represents the turbulent heat-flux.

To solve the system of RANS equations 1.4, 1.5 and 1.6, where each barred quantity represents a Reynolds average, we can take two avenues: Either the turbulent stresses and heat fluxes are determined individually by solving a set of transport equations for each component (a total of six + three), or one introduces proper closure relations for τ_ij and the turbulent heat-fluxQ′′_j . The closure relations must provide a physically coherent representation of turbulence mechanisms.

1.4 The Reynolds Averaged Scalar Equation

In the same way as Navier-Stokes equations, the scalar transport equations can be averaged using the Reynolds-averaging procedure. This gives the following equation

D C- ∂ ( ∂C--) ∂ (---) -- ---- = ---- D ---- - ---- c′u′ + F C (1.7) Dt ∂xj ∂xj ∂xj

where

is the molecular diffusivity, u′c′ is the turbulent scalar transport and F_C is the averaged source term applied to the scalar.

1.5 The Reynolds Stress Equations

The Reynolds stress, τ_ij can also be determined by solving its own transport equation. A special notation is introduced for this purpose, namely the Navier-Stokes operator (u′_i) producing the fluctuating–momentum equations

′ Du-′i 1∂p-′ 2 ′ L(ui) = Dt + ρ∂xi - ν ∇ ui = 0 (1.8)

This operator helps to directly derive an equation for the Reynolds stress tensor through construction of the following second order moment:

-′---′-----′----′- uj L(ui) + uiL (uj) = 0 (1.9)

After some calculations, one arrives at the Reynolds-stress transport equations:

D τij ∂ [ ∂ τij ] ----- = Pij + Πij - εij +---- ν---- - Cijk (1.10) Dt ∂xk ∂xk

where

-- -- Pij = - τik-∂vj- τjk ∂-vi (1.11) ∂xk ∂xk

--(----′--∂u-′)- Πij ≡ p′ ∂u-i+ ---j (1.12) ∂xj ∂xi

------′--′- εij ≡ 2μ ∂u-i∂uj (1.13) ∂xk ∂xk

---′-′-′ -′-′ -′-′ Cijk ≡ ρ uiujuk + p uiδjk + p uj δik (1.14)

designate, respectively, the mechanical production of turbulence P_ij , or the source of Reynolds stress, the pressure-strain correlation or inter-component transfer term Π_ij , the dissipation (by the action of viscosity) rate tensor ε_ij , and the turbulent transport term (C_ijk), acting in addition to the molecular diffusion of τ_ij represented by ν ∇² τ_ij .

The pressure-strain term Π_ij , gives rise to isotropy of the turbulence field by redistributing the turbulent kinetic energy among its three components and by reducing the shear stresses (Hinze, 1975).

1.6 The Dissipation Term

A typical practice in modelling the energy dissipation term ε_ij (Eq. 1.13) is to assume that the dissipative eddies are isotropic (Hinze, 1975) and thus

---′---′ εij = 2-εδij ; ε ≡ 1εii = ν ∂ui-∂ui (1.15) 3 2 ∂xk∂xk

This hypothesis, known as the Kolmogorov assumption of local isotropy, supposes that small turbulent structures associated with the rate of dissipation do not have a preferential direction. First, a transport equation for the dissipation rate tensor ε_ij has to be established. For the sake of brevity we restrict ourselves to the exact transport equation for ε itself (not for ε_ij ), which can formally be obtained by constructing the moment:

----------- ∂u-′i∂L-(u-′i)- 2 ν∂xj ∂xj = 0 (1.16)

The equation takes the form

D ε --- = P ε1 + Pε2 + Φ ε + D ε + ν∇2 ε (1.17) Dt

where

-- ( --′---′------′---′) Pε1 = - 2ν ∂vi- ∂uk-∂uk-+ ∂ui-∂uj- (1.18) ∂xj ∂xi ∂xj ∂xk ∂xk ------- ′∂u-′i--∂2vi- - 2νu k∂xj ∂xj ∂x ------------k ∂u′ ∂u ′ ∂u′ Pε2 = - 2ν --k----k---i (1.19) ∂xi ∂xj ∂xj (---2-′--)2- Φ ε = - 2ν2 -∂-ui-- (1.20) ∂xj ∂xk ( -------- -----------) -∂-- ∂p′-∂u′j- u′j-∂u′i ∂u′i- D ε = - 2ν ∂xj ∂xi ∂xi + 2 ∂xk ∂xk (1.21)

correspond to the following physical processes: The production by the mean motion (i.e by the mean velocity gradients) P_ε1 , the generation rate of vorticity fluctuations through the self-stretching of vortex tubes P_ε2 , the decay or destruction of the dissipation Φ_ε through the action of viscosity, and the turbulent diffusion of dissipation

_ε (diffusion by molecular effects is represented by ν∇²ε).

1.7 The Turbulent Kinetic Energy

A modelled version of the exact transport equation (1.10) can be derived, by letting k ≡ τ_ii∕2:

-- ( ------ ---) Dk- = - τij ∂vi-- ε - -∂-- 1u ′u ′u′ + p′u ′ + ν∇2k (1.22) Dt --- ∂xj ∂xj 2 i ij j ◟ ◝P◜ ◞

In this equation,

≡ P_ii∕2 represents the mechanical production of turbulence due to the interaction between turbulent stresses and mean velocity gradients. Note also that the contraction of τ_ij yields Π_ii = 0 by virtue of u′_i,i = 0.

1.8 Chemical species equations

Chemical species transport equations are modelled using their respective mass fraction Y _k, where k represents the species index. When the species have different densities, the usual Reynolds average gives new unclosed terms ρ′u_j′ in the mass balance equation corresponding to the correlation between density and velocity fluctuations. To avoid this difficulty, mass-weighted averages (called Favre averages) are usually preferred

ρf- f^= --- ρ

(1.23)

where f is a generic quantity. Splitting Y _i into a Favre mean Y^i and a fluctuation Y _i′′ the transport equation can be written

^ ( ---------) ( ) ρD-Yk = -∂-- (ρDk∇Yk ) - -∂-- ρv ′′Y′′ + ρ^Sk St ∂xj ∂xj k

(1.24)

where D_k is the molecular diffusion coefficient for species k and ^S
k is the Favre average of the source term.

Chapter 2
Two-Equation Turbulence Models

There are different ways to model the Reynolds–stress tensor τ_ij (when the transport equations of τ_ij are not solved explicitly). The best-known approach and the most practical one is based on the eddy viscosity concept or Boussinesq Hypothesis, which assumes that the Reynolds stress is decomposed into an isotropic and a deviatoric part,

2 1 -- -- τij = 3 δijk - 2 νtSij; Sij = 2-(vi,j + vj,i) (2.1)

where δ_ij = 1 for i = j, and zero otherwise. The deviatoric part (2nd term in the RHS of Eq. (2.1)) is a symmetric, traceless tensor and links τ_ij linearly to the rate–of–strain tensor S_ij . The coefficient of proportionality, ν_t , designates the eddy viscosity (or turbulent viscosity), which is a characteristic of the flow.

2.1 The Standard k - ε Turbulence Model

In the k -ε model (Launder & Spalding, 1974), the local state of turbulence is characterised by kinetic energy k and its dissipation rate ε. The combination of ℓ₀ ≈ k^3∕2∕ε and τ₀ ≈ k∕ε (or ℓ₀ and v_I ≡ √ --
k ) leads to the generalised form of the isotropic eddy viscosity ν_t :

νt ≡ ℓ0vI ≡ ℓ20∕τ0 = Cμk2∕ε (2.2)

where C_μ is a constant. The distribution of turbulent kinetic energy k and its rate of dissipation ε appearing in this relation are determined from the following model transport equations:

Dk ∂ ( νt ∂k ) Dt- = ∂x-- σ--∂x-- + P - ε (2.3) j ( k j) 2 D-ε = -∂-- -νt∂-ε- + C P ε-- C ε- (2.4) Dt ∂xj σ ε∂xj ε1 k ε2 k

The reader should refer to Launder & Spalding (1974) for a better understanding of the derivation of the model from the exact transport equations. The empirical constants are assigned the following standard values: C_μ = 0.09; C_ε1 = 1.44; C_ε2 = 1.92; σ_k = 1. and σ_ε = 1.3. The standard model (Eqs. 2.2 - 2.4) with its original coefficients is valid only in flow regions away from solid boundaries (high-Re-number flow regions), and depending on whether it is to be employed in low- or high-Re forms, special near-wall treatments are needed..

2.2 The RNG k - ε Turbulence Model

The RNG k - ε model, is more advanced than standard model and was derived from renormalization group theory. The distribution of turbulent kinetic energy k and its rate of dissipation ε appearing in this relation are determined from the following model transport equations:

( ) Dk- = -∂-- νeff-∂k- + P - ε (2.5) Dt ∂xj ( σk ∂xj ) D ε ∂ νeff ∂ε ε ε2 3 3 ε2 Dt- = ∂x-- -σ--∂x-- + C ε1P k-- C ε2k-- Cμη (1 - η∕η0)∕(1 + βη ) k- (2.6) j e j

where η = Sk∕ε; η₀ = 4.38; β = 0.012 and S denotes strain rate S = ∘ -------
2SijSij

. Effective viscosity instead of turbulent viscosity is used here:

∘ ---- ν = ν (1 + k Cμ)2 (2.7) eff νε

The empirical constants are assigned the following standard values: C_μ = 0.0845; C_ε1 = 1.42; C_ε2 = 1.68; σ_k = 0.71942 and σ_ε = 0.71942. The RNG model (Eqs. 2.5 - 2.7) with its original coefficients is valid only in flow regions away from solid boundaries (high-Re-number flow regions), and depending on whether it is to be employed in low- or high-Re forms, special near-wall treatments are needed. In the energy/temperature equation a new formula for heat diffusion coefficient is employed which is also valid only for high-Re-number flow:

keff = 1.3929 cpνeffρ (2.8)

In the mass transfer equations/concentration equations a new formula for mass diffusivity is employed:

Deff = 1.3929 νeff (2.9)

2.3 Extensions to k - ε Turbulence Model

Depending on the flow conditions, different extensions to the basic k - ε model have been proposed.

2.3.1 Yap correction

The Yap correction (Yap, 1987) consists of an additional source term of the form ρS_ε to the right hand side of the epsilon equation. The source term can be written as,

ε2 (k1.5 ) ( k1.5)2 ρS ϵ = 0.83ρ-- ----- 1 ---- k εle εle

(2.10)

where,

-0.75 le = cμ κyn

(2.11)

with y_n being the normal distance to the nearest wall.

The Yap correction is active in non-equilibrium flows and tends to reduce the departure of the turbulence length scale from its local equilibrium level. Yap (1987) showed improved results with the k - ε model in separated flows when using this extra source term.

2.3.2 Swirl correction

The turbulent eddy viscosity is corrected for strongly swirling flows because the standard k -ε models tend to overestimate turbulent diffusion for swirling flows. The eddy viscosity is modified using a swirl correction factor, f_sw, as follows (Shih et al., 1997),

k2- νt = fswCμ ε

(2.12)

where,

-----1------- fswC μ = kU * 4.0+ AS -ε--

(2.13)

and

√-- AS = 6 cos(ϕ )√ -- ϕ = 1arccos( 6W *) 3 W * = SijSjkSki- -S3----- U* = ∘ S2 + Ω2 ∘ ------ S = ∘ SijSij- Ω = ΩijΩij

2.3.3 Compressible correction

Initially turbulence models were designed for low-speed, isothermal flows. The compressibility correction is devised to deal with additional effects seen for higher Mach number flows, specifically, the effects of compressibility on the dissipation rate of the turbulence kinetic energy. For free shear flows, this is exhibited as the decrease in growth rate in the mixing layer with increasing Mach number. Standard turbulence models do not account for this Mach number dependence, and thus a compressibility correction is used. For compressible flows, two extra terms, known as the dilatation dissipation, ε_d, and the pressure-dilatation occur in the turbulence kinetic energy equation. The pressure-dilatation term is usually neglected because its contributions have been shown to be small (Wilcox, 1993).

The dilatation dissipation term is included in addition to the solenoidal, or incompressible, dissipation, ε. Thus, the effect is that the growth rate of turbulence is inhibited when the correction is active. Sarkar et al. (1991) modeled the ratio of the dilatation dissipation to the solenoidal dissipation, ε_d∕ε, as a function of the turbulence Mach number, M_t, defined as

M 2= 2k- t c2

(2.14)

where c is the speed of sound. For the model by Sarkar et al. (1991), the additional source term to the k equation is,

2 ρεd = ρεM t

(2.15)

Chapter 3
Near–Wall Modelling

3.1 Near–Wall Treatment in High–Re Flows

The near-wall treatment to be presented next applies strictly to the standard k -ε model (Eqs. 2.2 - 2.4), Their combination is referred to as the high-Re k - ε model. This near-wall modelling consists in bridging the viscous sub layer, i.e. the first grid point must be located outside this zone.

3.1.1 Standard Wall–Function Approach

At a rigid boundary, the no-slip condition in setting the velocity of the fluid to zero (or to that of the wall to account for moving-boundaries) is used for laminar flows. In high-Re flows, the wall-function approach, whose theoretical basis is discussed next, is applied in place of the no-slip condition. This treatment is based on two assumptions: (i) The immediate near-wall region is in a state of local equilibrium ( = ε), and (ii) the velocity profile merges with the ”log-law” given by

-- U+ = u--= 1-ln(E y+ ) (3.1) uτ κ

However one needs to make the extra assumption , that, the equilibrium layer near the surface, where production balances dissipation, is one-dimensional and the stress across is constant. In these circumstances, it can be demonstrated that

k 1 --p2= --1∕2- or uτ = k1p∕2C1μ∕4 (3.2) u τ Cμ

where the subscript p refers to the first control volume centre relative to the wall. Laboratory data for flow over smooth surfaces indicate that k_p∕u_τ² ≈ 3.3, and hence C_μ ≈ 0.09. In a second step and making use of (3.2) we can relate the induced retarding wall shear-stress

_w to the velocity vector

_p by

⃗τ = - λ ⃗V (3.3) w w p

where

{ + λw = μ∕yp1∕4 1∕2 if yp < 11.6 (3.4) ρC μ kp κ∕ ln (Ey+p ) otherwise

with

y+p = ρ C1μ∕4k1p∕2yp∕ μ, κ = 0.41 (3.5)

Furthermore, the diffusive flux of k is equal to zero at the wall. In a second step, the wall shear-stress (supposed to be uniform and equal to that at the wall, y ≥ y_w ) acting at a distance y_n = 2y_p from the wall induces a rate of turbulence production given by

| 1 ∫ yn ∂u-|| τ2w Pw = - y-- τw ∂y-|| dy = κμy+- (3.6) n yw w p

This quantity (

_w) has to replace the production term

in the transport equation of k (2.3). Likewise, over the fully turbulent region (y_w ≤ y ≤ y_n ) where k was taken to be uniform (in the pace interval y_w ≤ y ≤ d_p ), the mean value of ε_w can be obtained by the expression

∫ yn 3∕2 3∕4 3∕2 εw = 1-- C Δ kp--⋅dy = C-μ--kp- (3.7) yn yw y κyp

It should be noted, however, that since the assumption of ”local equilibrium” near the wall is not justifiable in the viscous sub layer (y_p⁺ < 11.6), the wall-function approach requires the value of y_p⁺ at the cell adjacent to the wall to fall in the range 11.6 ≤ y_p⁺ ≤ 200 - 500, otherwise, the result is systematically biased.

3.1.2 Two–layer turbulence models

The two–layer approach adopted here consists of resolving the viscosity–affected regions close to walls with a one–equation model, while the outer core flow is resolved with the standard k - ϵ model described above. In the one–equation model, the eddy viscosity is made proportional to a velocity scale and a length scale l_μ. The distribution of l_μ is prescribed algebraically while the velocity scale is determined by solving the k–equation (as in Eq. (2)). The dissipation rate ε appearing as sink term in the k–equation is related to k and a dissipation length scale l_ε which is also prescribed algebraically. The different two-layer versions available in the literature differ in the use of the velocity scale and the way l_μ and l_ε are prescribed. It should be mentioned that in the fully turbulent region the length scales l_μ and l_ε vary linearly with distance from the wall. However, in the viscous sublayer l_μ and l_ε deviate from the linear distribution in order to account for the damping of the eddy viscosity and the limiting behaviour of ε at the wall.

k^1∕2 velocity scale based model Rodi (1991): TLK

The approach combines the standard k - ϵ model in the outer region with a one–equation model due to Norris and Reynolds Norris & Reynolds (1975) in the viscous–sublayer employing

νt = Cμk1∕2lμ; ε = k3∕2∕lε (3.8)

In this model, the length scale l_μ is damped in a similar way as the Prandtl mixing length by the Van Driest function, so that it involves an exponential reduction governed by the near–wall Reynolds number R_y = Uy_n∕ν. However, in contrast to the original Van Driest function, R_y uses k^1∕2 as a velocity scale U instead of U_τ which can go to zero for separated flows.

( R ) lμ = Cl ynfμ with: fμ = 1- exp - --y-25+ (3.9) A μA

The constant C_l is set equal to κC_μ^-3∕4 to conform with the logarithmic law of the wall. The empirical constants appearing in the f_μ–function are assigned the values A_μ = 50.5 and A⁺ = 25. The reader is referred to Rodi (1991) for a review and further details on the choice of the constants. For the dissipation scale the following distribution is used near the wall:

l = -----Clyn-----; C = 13.2 (3.10) ε 1 + Cε∕(RyCl ) ε

The outer (k - ϵ) and the near–wall model are matched at a location where the damping function f_μ reaches the value 0.95, i.e. where viscous effects become negligible.

The combination of the Kato–Launder correction with the TLK model is hereafter labeled TLKK.

(v′²)^1∕2 velocity scale based model Rodi & Mansour (1993): TLV

The development of this model was motivated by the fact that the length scales-functions as proposed in Norris & Reynolds (1975) particularly the l_ε–function, are not in agreement with direct numerical simulation (DNS) data, and that the normal fluctuations (v′²)^1∕2 are a more relevant velocity scale for the turbulent momentum transfer near the wall than is k^1∕2. Therefore, the following model using (v′²)^1∕2 as a velocity scale was proposed in Rodi & Mansour (1993)

∘ --- ∘ --- --2 --2 νt = v′lμ,ν, ε = v′ k∕lε,ν (3.11)

[ ∘ --- ] -′2- with lν,μ = 0.33yn, and lε,ν = 1.3yn∕ 1.+ 2.12ν ∕ v yn (3.12)

which is based on DNS data for fully developed channel flow. As an equation for k is solved, v′² needs to be related to k, which is done through the following DNS based empirical relation

∘ --- v′2∕k = 4.65× 10-5(Ry )2 + 4.00× 10- 4Ry, Ry = k1∕2yn∕ν (3.13)

valid only very near the wall. The matching between the outer and near–wall model is performed at a location where R_y = 80.

3.2 The Low-Re k - ε Turbulence Model

The standard k - ε model was developed for high Reynolds number flows and is therefore not valid in flow regions very close to the wall, i.e. within the viscous sub layer. In this approach, the following relation for ν_t was proposed; it includes a damping function, f_μ , varying from almost zero near the wall to ≈ 1 at the outer edge of the viscous sub layer:

νt = Cμfμk2∕ε (3.14)

Experiments and DNS have also shown that the model constants C_ε1 and C_ε2 appearing in the source terms of the transport equation for ε need also to incorporate damping functions f₁ and f₂ in order to reproduce the steep gradient of ε near the wall, so that equation (2.4) takes the form

[( ) ] 2 D-ε = -∂-- ν + νt- ∂ε--+ C ε1f1P ε-- Cε2f2 ε + Ξ (3.15) Dt ∂xi σε ∂xi k k

where the molecular diffusion of ε is now re-introduced, unlike in (2.4) for high-Re-number flows. In this equation Ξ represents an additional term to account for the fact that dissipation processes are not isotropic within the viscous sub layer (Jones & Launder, 1972). The different low-Re models proposed so far differ either in the use of the model functions f_μ , f₁ and f₂ and the term Ξ, or the way the dissipation rate is obtained.

3.2.1 The Jones & Launder and Launder & Sharma Models

The Low-Re k - ε model employed most frequently is due to Jones & Launder (1972). This model is

[ 2] 2 fμ = exp - 3.4∕(1+ Rt∕50 ) ; Rt = k ∕ν ε (3.16) ( 2) -- 2 f1 = 1.0 ; f2 = 1- 0.3exp - R t ; Ξ = 2ννt(v,yy) (3.17)

The model constants C_ε1 and C_ε2 are assigned the same values as those of the standard k -ε model, i.e. C_ε1 = 1.44 and C_ε2 = 1.92 (as in Launder & Sharma (1974)) . This model is found to result in a very rapid growth of f_μ towards 1 (i.e. towards the edge of the viscous sub layer).

3.2.2 The Abe-Kondoh-Nagano Model

Various new low–Re k - ε have been proposed. One of them is proposed by Abe et al. (1994)

* 2[ - 3∕4 ( 2)] fμ = [1 - exp (- y ∕14)] 1+ 5 Rt exp - [Rt∕200] (3.18) [ ( )] f = 1; f = [1 - exp (- y*∕3.1)]2 1 - 0.3exp - [R ∕6.5]2 (3.19) 1 2 t

The wall boundary conditions used together with this type of model require the dissipation rate at the wall ε|_wall to be adjusted so as to reproduce the correct asymptotic behaviour, i.e. ε|_wall = 2νk∕y_p². In contrast to the wall–function approach used with the standard k -ε model for high–Re number flows, in the low-Re schemes the no-slip condition is to be employed for the velocity, along with the condition of zero turbulent kinetic energy at the wall. Additionally, the model should meet a criterion similar to the one required for the WF approach, i.e. y_p⁺ < 0.1 - 0.5. In general, a typical layer of about 30 grid–points lying within the viscous sub layer (where f_μ < 0.95) is necessary to correctly reproduce the steep gradient of the dissipation rate.

An alternative approach to the use of the damping function f₂ to avoid having to face a singularity in the destruction term -C_ε2 ε²∕k consists in replacing this term by

ε- - C ε2T (3.20)

so as to allow the turbulent time scale

to recover τ₀ = k∕ε far from boundaries – at high–Re turbulence –, and makes it proportional to the Kolmogorov time scale ~ ∘ ----
ν∕ ε

for low-Re near wall turbulence (Durbin, 1991, 1993):

[k ( ν)1∕2] T = max -; CK -- (3.21) ε ε

where C_K is a constant of order one.

3.2.3 The Lam-Bremhorst Model

Lam & Bremhorst (1981) extended the high Reynolds number form of the k - ε model and tested by application to fully developed pipe flow. The relationship is of the form:

( 20.5) 2 fμ = 1 + Re-- [1 - exp(- 0.0165Rey )] t

(3.22)

(0.05)2 ( ) f1 = 1 + ---- ; f2 = 1- exp - Re2t fμ

(3.23)

2 0.5 Ret = k-; Rey = k--ymin- νε ν

(3.24)

Chapter 4
Turbulent Heat flux Modelling

In order to model the turbulent heat fluxes, the Reynolds Averaged Navier-Stokes equations (RANS) and the Reynolds Averaged Energy equation introduced in the first chapter need to be solved.

-- ∂-vi ∂xi = 0 (4.1) Dvi- 1 ∂p- -- ∂τij ---- = - -----+ ν ∇2vi - ----+ ρgi; i = 1,2,3 (4.2) Dt- ρ∂x(i -- ∂xj) DT-- --∂- -∂T- -′-′ ρcp Dt = ∂xj λ∂xj - ρ cpujθ (4.3) ◟-------◝◜-------◞ - q′′total

The different modelling approaches for turbulent heat flux modelling are presented in Fig. 4.1 wherein, the simplest approach is called the Simple Gradient Diffusion Hypothesis (SGDH) and the most complex one is the Turbulent Heat Flux Modelling (THFM) approach which requires the solution of 5 additional equations to model the turbulent heat flux.

Figure 4.1: Turbulent heat flux modelling overview, where k_T represent the variance of the temperature fluctuations θ′² and e_T represent the dissipation rate of the heat fluxes ε_θ′.

4.1 Turbulent Heat Flux Model (THFM)

The most accurate approach to account for thermal effects and heat transfer within the time-averaged context is to solve the transport equation for the turbulent heat fluxes u′_iθ′, by incorporating a detailed modelling of the buoyancy term. In this model, transport equations for the variance of the temperature fluctuations θ′² and its dissipation rate ε_θ′ have to be solved in addition to close the set of equations.

The modeled transport equations for turbulent heat fluxes are given by: (for i=1,2,3) Carteciano & Grötzbach (2003)

( ) Du-′θ′ ∂ [ k2 κ+ ν ] ∂u-′θ′ ( ---- ∂T- ---- ∂u-) ---i--= ---- CT D---+ ----- ---i- - u′iu′j----+ u ′jθ′ --i- - Gu-′θ′ + πi + εu′θ′ Dt ∂xj ε 2 ∂xj ∂xj ∂xj i i

(4.4)

where κ is the thermal diffusivity: κ = -λ-
ρcp .

The term π_i represent the modeling of the pressure-temperature gradient correlation with Monin’s and Launder’s correlation.

---- ---- -- --- ----- 3 πi = - CT 1εu′θ′ + CT 2u′θ′∂-ui+ CT 3βgiθ′2 - CT 4εu′ θ′δi[n] k-2- k i j ∂xj k [n] x [n]ε

(4.5)

where δ_i is the Kronecker delta and the index n represent the normal direction to a wall.
The dissipation rate of the heat fluxes ε_{u′_iθ′} can be modelled by:

--- -1-+-Pr--( ε) -′-′ εu′iθ′ = - 2√P-r√R- k exp [- CT 5(Ret + Pet)]uiθ

(4.6)

Where Re_t is the local turbulent Reynolds number, Pe_t = Re_t × Pr is the local Peclet number and R is the turbulent time-scale ratio. ε_{u′_iθ′} is negligible in the transport equation (4.4) of the heat fluxes at high Peclet numbers but its contribution is important at low Peclet numbers.
Another expression for the dissipation rate is available for very low Peclet numbers:

1 ( 1 ) (P r)0.7( ε) ---- εu′iθ′ = - 2 1 + Pr- -R- k- u′iθ′

(4.7)

The buoyancy term G_{u′_iθ′} is expressed in terms of the variance of the temperature fluctuations θ′².

G ′′ = βg θ′2- uiθ i

(4.8)

For a detailed description of the buoyancy effects, a transport equation of the variance of the temperature fluctuations θ′² is solved:

--′2 -′2 [( 2 ) -′2] ---- -- ( ∘ -′2)2 ∂ρθ-- + ∂ρujθ--= -∂-- ρCT T k-+ ρ κ ∂θ-- - 2ρ u′θ′ ∂T-- 2ρ εθ′ - 2 ρκ ∂-θ-- ∂t ∂xj ∂xj ε ∂xj j ∂xj ∂xj

(4.9)

The previous transport equation requires the dissipation rate of the temperature variance ε_θ′ to be modeled using an other transport equation.

[( 2 ) ] ( -′-′ -- ) ∂ρεθ′+ ∂ρujεθ′= -∂-- ρC k--+ ρκ ∂εθ′ - ρε ′ C εθ′+ C ε-+ C ujθ-∂T--- C Pk- ∂t ∂xj ∂xj DD ε ∂xj θ D1 θ′2 D2k P 1θ′2 ∂xj P 2k ( -- )2 + 2ρ κκ -∂2T--- (4.10 ) t ∂xk∂xj

with:

( -- --) -- P = νt ∂ui+ ∂uj- ∂-ui k ∂xj ∂xi ∂xj

(4.11)

4.2 Full Algebraic Turbulent Heat Flux Model: θ′² - ε_θ′

In this model, the turbulent heat flux u′_iθ′ is modelled algebraically instead of solving an additional set of partial differential equations. The model is referred to as the Algebraic Turbulent Heat Flux Model (ATHFM) (Launder, 1988). The turbulent heat flux (u′_iθ′) can be approximated by its production rate times the turbulent time scale. The production terms and turbulence fluctuations from the differential transport equations are locally in balance S. Kenjeres (2000).

---- (---- -- ---- -- --) u′θ′ = - Cθ′ k u′u′ ∂T-+ ξu′θ′ ∂ui-+ ηβ giθ′2 i ε i j ∂xj j ∂xj

(4.12)

In order to close the above system, the values of θ′² and ε_θ′ are needed. In the case where transport equations are solved for both the above quantities we obtain the ATHFM -θ′² - ε_θ′ model.

--- --- ( ---) ∂ρθ′2 ∂-ρujθ′2 -∂-- ∂θ′2 ∂t + ∂xj = ∂xj ρ(κ + κt)∂xj + 2ρP θ′ - 2ρεθ′ (4.13)

and

∂ρεθ′ ∂ρujεθ′ ∂ ( ∂ ˜εθ′) ′ ε˜θ′ ′ ˜εθ′ ----- + -------= ---- ρ (κ + κt)---- + C θε1ρP θ′′2 + Cθε3ρP --- ∂t ∂xj ∂xj ∂xj θ k θ′ ˜ε2θ′- θ′ ′ε ˜θ′˜ε ′ - C ε4ρ θ′2 - Cε5ρfεθ k + Eθ (4.14)

where:

---- ∂u- ----∂T- ---- ( ∂2T )2 P = -u ′iu′j---i, Pθ′ = -u′jθ′---, G = - βgiuiθ′, E θ′ = 2ρκκt ------- , ∂xj ∂xj ∂xj∂xk

( √ -)2 ( ∘ -′2)2 ˜ε = ε - 2ν ∂--k- , ε˜θ′ = εθ′ - κ ∂--θ-- , ∂xk ∂xk

[ ] k2 - 3.4 κt = C Φfμ--, f μ = exp (-------)2 , fεθ′ = 1. ˜ε 1 + Re50t

In Kenjereš et al. (2005), the model for the Reynolds stress was extended with a buoyancy and turbulent heat flux term given as,

---- ( ) ( ---- ---- ---- ) u′u′ = 2kδij - νt ∂ui-+ ∂uj- + C θτβ giu′θ′ + gju′θ′ - 2gku′θ′δij (4.15) i j 3 ∂xj ∂xi j i 3 k

with C_θ = 0.15.

This model was further developed and calibrated by Shams et al. (2014) where the eddy viscosity ν_t was given as,

νt = fμC μ[kτ ] (4.16) k- (ν-) τ = max ε,CT ε , and [ ( ∘ ---- 2)] fμ = 1 - exp - Cd0 Red + Cd1Red + Cd2Re d

Where Re_d = √--
k

d∕ν with d being the distance to the wall, and the constants being, C_μ = 0.09,C_T = 0.6,C_d₀ = 0.091,C_d₁ = 0.0042,C_d₂ = 0.00011.

The algebraic turbulent heat flux model used (Shams et al., 2014; Kenjereš et al., 2005) involves an additional term proportional to the Reynolds stress anisotropy. The algebraic turbulent heat flux is also redefined with different names for the model constants and is given as,

---- ( ---- ---- --) ---- u′θ′ = - Ct0τ Ct1u′u′∂T--+ Ct2u′θ′∂ui-+ Ct3βgiθ′2 + Ct4aiju′θ′ (4.17) i i j∂xj j ∂xj j

with g_i the gravity vector and a_ij the Reynolds stress anisotropy tensor defined as,

---- u′iu′j 2 aij = --k- - 3δij

(4.18)

The transport equation for θ′² is solved and ε_θ′ is evaluated from the thermal to mechanical time-scale ratio (Eq. (4.23)) where an algebraic expression is used by Kenjereš et al. (2005) and a constant value of 0.5 by Shams et al. (2014). The assumption of the constant time-scale ratio R worked reasonably well in a number of flow regimes when compared with the respective reference DNS or experimental data.

In the work of Shams et al. (2014), the model proposed by Kenjereš et al. (2005) is referred to as AHFM-2005 which was calibrated for natural and mixed convection flow regimes close to unity Prandtl number. Shams et al. (2014) proposed a new calibration for forced convection and low Prandtl number fluids (liquid metal) referred to as AHFM-cc. The model was further extended by introducing a logarithmic dependence of coefficient C_t₁ on the Reynolds and Prandtl numbers with a constraint to be positive for RePr > 180. This variant of AHFM-cc was called AHFM-NRG. The model constants for the above models are given in Table 4.1.

C_t₀

C_t₁

C_t₂

C_t₃

C_t₄

AHFM-2005

0.15

0.6

1.5

0.5

AHFM-cc

0.2

0.25

0.6

2.5

0.0

0.5

AHFM-NRG

if RePr > 180

0.2

0.053 x ln (ReP r)

- 0.27

0.6

2.5

0.0

0.5

Table 4.1: Model constants for the algebraic heat flux model

Note that the nomenclature of the constants used by Shams et al. (2014) will be used in this report. In the model presented by Kenjereš et al. (2005), an additional coefficient C_t₁ has been introduced as compared to S. Kenjeres (2000) which was by default set to one in earlier publications.

4.2.1 Equivalence between THFM and ATHFM heat flux modelling

The algebraic expression for the turbulent heat flux can be obtained from the partial differential model (Eq. 4.4) by simply dropping the unsteady, advection, and diffusion terms to obtain,

-- ---- 1 (k ) [----∂ T ----∂ui --] u ′iθ′ = -C--- ε- u′iu′j∂x--+ (1 - CT2)u′jθ′∂x--+ (1- CT 3)βgiθ′2 T1 j j

(4.19)

where the near wall term in π_i and the dissipation rate of the turbulent heat fluxes have been ignored. By comparison with the algebraic model proposed by Kenjereš et al. (2005) in Eq. 4.17, we note that,

1 Ct0 = C--- T1 Ct1 = 1 Ct2 = 1 - CT 2 Ct3 = 1 - CT 3 (4.20)

The work of Shams et al. (2014) has shown that C_t₁ values are much different from 1.0. For example, Kenjereš et al. (2005) propose a value of 0.6 which is further reduced by Shams et al. (2014) for low Prandtl numbers to 0.25 in the so called AHFM-cc model. In the AHFM-NRG model proposed by them, C_t₁ is made into a function of Reynolds and Prandtl numbers with values always below 1.0. With this in mind, the THFM model is extended by introducing a new coefficient and a new nomenclature where,

Ct0 = -1-- CT0 Ct1 = 1 - CT 1 Ct = 1 - CT 2 2 Ct3 = 1 - CT 3 (4.21)

such that a one-to-one equivalence is established between the algebraic model of Kenjereš et al. (2005) and the differential model of Carteciano & Grötzbach (2003). Due to this alteration the expression for π_i in the transport equations for the turbulent heat fluxes changes to,

-- -- 3 π = - C εu′θ′ + C u′u′∂T--+ C u-′θ′∂-ui+ C βg θ′2-- C -εu′-θ′δ -k-2- i T0k i T1 i j∂xj T 2 j ∂xj T3 i T4k [n] i[n]x[n]ε

(4.22)

The final set of model coefficients for the THFM model is given in Table 4.2.

u_iθ′ tr. equation

θ′² tr. equation

ε_θ′ tr. equation

coeff.

value

coeff.

value

coeff.

value

C_TD

0.11

C_TT

0.13

C_DD

0.13

C_T0

3.0

C_D1

2.2

C_T1

0.0

C_D2

0.8

C_T2

0.33

C_P1

1.8

C_T3

0.5

C_P2

0.72

C_T4

0.5

Table 4.2: Set of empirical coefficients for the equivalent THFM

4.3 Algebraic Turbulent Heat Flux Model: θ′²

The ATHFM -θ′² - ε_θ′ model can be further simplified by calculating the dissipation rate of the temperature variance ε_θ′ using the definition of the mechanical to thermal turbulent time-scale ratio R. The time scale ratio R is either set to a constant value or prescribed algebraically. The dissipation rate of temperature fluctuation variance can then be directly obtained as,

1 (ε) --- εθ′ = 2R k- θ′2

(4.23)

This represents a measure of the relative importance of the relaxation timescales of the mechanical and thermal dissipation. With this simplification, the need to solve a transport equation for ε_θ′ doesn’t exist. This model is referred to as the ATHFM -θ′² model. Although the value and variation of R in most situations is not known, such an assumption has displayed remarkable success in a number of thermal flows driven by gravitational effects. Typically, a value of 0.5 is used.

4.4 WET Model

Further simplification is achieved by determining the turbulent heat flux by applying the WET (Wealth ≡ Earnings × Time) theory, a syllogism applied by Launder (1988) to turbulent heat fluxes which lead to the expression: (Value of Second Moment ≡ Production Rate of Second Moment × Turbulent Time Scale). This, together with the turbulent time k∕ε, yields,

( -- --) u′θ′ = - C k- C u′u′ ∂T--+ C u′θ′ ∂ui (4.24) i t0ε t1 i j ∂xj t2 j ∂xj

The WET model is supposed to remedy the drawback of simpler variants, in which the heat flux is only generated by temperature gradients; which is not always the case. For example, the mixed layer formed close to a heated wall featuring a uniform vertical temperature gradient is not necessarily linked to turbulence, so the heat flux is actually over-represented in the relative sense. The same is true when vertical temperature gradients are small: here it is the velocity gradients that cause the wall-to-flow heat transfer.

4.5 Generalized Gradient Diffusion Hypothesis (GGDH)

The WET model on further simplification yields the GGDH model. The turbulent heat-flux is modelled in an manner analogous to the turbulent transport term in the Reynolds stress equation, in particular by reference to Daly & Harlow (1970) model.

---- ( ---- --) u′θ′ = - Ct k Ct u′u′-∂T- i 0ε 1 i j∂xj

(4.25)

This approach is known as the anisotropic eddy-diffusivity model, or the generalized gradient diffusion hypothesis (GGDH), a definition that points to the fact that heat transfer is driven by an anisotropic thermal diffusion (Λ_ij),

k ---- Λij = Ct0ε-Ct1u ′iu′j (4.26)

This approach has the merit to conform to many experimental findings, including the measurements of turbulent heat transfer in pipes and boundary layer flows by Bremhorst & Bullock (1973), and by many others. Indeed, these authors demonstrated that turbulent heat flux in the flow direction are two to three times larger than in the direction normal to the wall while the streamwise temperature gradient is negligible compared to that normal to the surface.

4.6 Simple Gradient Diffusion Hypothesis (SGDH)

In the context of linear RANS modelling for convective heat transfer, the turbulent heat flux is linked to the temperature gradient via the expression below where the turbulent stress u_i^′u_j^′ is replaced by the by its trace u_i^′u_i^′ equal to 2k, the turbulent kinetic energy. This gives rise to the isotropic thermal-diffusivity model, known as the Simple Gradient Diffusion Hypothesis (SGDH).

---- u′iθ′ = --νt-∂T-- P rt∂xi

(4.27)

The turbulent Prandtl number Pr_t introduced in the above expression is typically set to a constant value of 0.9. However, two-equation models for thermal diffusivity have been developed in parallel with the k -ε model of turbulence. This idea was motivated by the fact that turbulent heat diffusion should also be characterized by a scalar (thermal) time scale that varies in space and time, just like the turbulent time scale τ = k∕ε. Such models are also referred to as variable turbulent Prandtl number models.

The model coefficients for the whole chain for models from WET to ATHFM -θ′² -ε_θ′ is given in Table 4.3.


C_t₀	C_Φ	C_ε1^θ′	C_ε3^θ′	C_ε4^θ′	C_ε5^θ′	C_t₁	C_t₂	C_t₃	C_t₄


0.2	0.09	1.3	0.72	2.2	0.8	0.25	0.6	2.5	0.0

Table 4.3: Coefficients for ATHFM

4.7 Variable Turbulent Prandtl Number Model

Within the Reynolds-averaged Navier Stokes model, where the turbulent heat flux is modelled by the SGDH model, the turbulent Prandtl number Pr_t is typically assumed to be a constant. This modelling approach limits the generality of these models because the turbulent Prandtl number can alter the mean flow prediction. One approach to extend the applicability of the SGDH model is to model the variation of the turbulent Prandtl number.

The turbulent thermal diffusivity is defined using a composite time scale using both the mechanical and thermal turbulent time scales given by,

∘ ----- k θ′2 κt = Cλfλk ----- ε εθ′

(4.28)

where C_λ is a model constant and f_λ is a near-wall damping function for the turbulent thermal diffusivity applicable for low Reynolds number turbulence models. Using the standard expression for eddy diffusivity, the turbulent Prandtl number can be expressed as,

∘ ----- P r = C-μfμ k-εθ′ t C λfλ ε θ′2

(4.29)

In order to close Eq. (4.29) the variance of the temperature fluctuations, and its dissipation rate are required. If transport equations of these quantities are solved, we obtain the 2-equation variable Pr_t turbulent heat flux model.

The equations governing the transport of temperature variance and its dissipation rate are given as, N. Chidambaram & Kenzakowski (2001)

--- --- ( [ ] ---) ( )2 ∂ρθ′2 ∂ρujθ′2- -∂-- -κt- ∂θ′2- ∂T-- ′ ∂t + ∂xj = ∂xj ρ κ+ σ θ′2- ∂xj + 2ρκt ∂xj - 2ρεθ

(4.30)

∂ρε ′ ∂ρu ε ′ ∂ ( [ κ ] ∂ε ′) ( ε′ ε ) ( ∂T )2 ε ′ ---θ-+ ---j-θ-= ---- ρ κ+ --t- --θ- + Cp1 -θ-+ Cp2-- ρ κt ---- + Cp3-θ-P- ∂t ∂xj ∂xj σεθ′ ∂xj θ′2 k ∂xj k ˆεθ′ ˆε- - Cd1θ′2ρεθ′ - Cd2k ρεθ′ + ξεθ′ (4.31)

where the operator

denotes the low Reynolds number modifications, ξ_{ε_θ′} is the near-wall correction function, and P is the turbulent kinetic energy production term. According to T.P. Sommer & Lai (1993), different values of C_p₁ can be used: a value of 1.8 for boundary layers and 2.0 for internal flows

( --)2 ∂ui- P-= μt ∂xj

(4.32)

The near-wall correction function is given as,

[ ′ *2 ′ ] ξεθ′ = fwρ (Cd1 - 4) ˆεθεθ′ + Cd2 ˆεεθ′ - εθ′ - (2 - Cp1 - Cp2P r)-εθ-Pθ′ θ′2 k θ′2 ρθ′2

(4.33)

where P_θ′ is the production due to the mean temperature gradient in the streamwise direction to account for a wall heat-flux boundary, the damping function is defined as,

[ ( )2] ( ∘ -′2-)2 ( √--)2 -′2 fw = exp - Ret- , ˆεθ′ = εθ′- κ ∂--θ-- , ˆε = ε- 2ν ∂--k- ε*′ = εθ′- κθ-- 80 ∂xj ∂xj θ d2

(4.34)

The near-wall damping function f_λ is given as,

fwC [ ( y+ ) ]2 fλ = ----11λ+ 1 - exp - --+ Ret4 A

(4.35)


C_p₁	C_p₂	C_p₃	C_d₁	C_d₂	σ_θ′²	σ_{ε_θ′}	C_λ	A⁺	C_1λ

2.0	0.0	0.72	2.2	0.8	1.0	1.0	0.14	45	0.1

Table 4.4: Coefficients for Variable turbulent Prandtl number model

Chapter 5
Buoyancy–Driven Turbulent Flows

5.1 The Equations of Motion

The buoyancy effects give rise to an additional body-force term f_i = -β g_i (T -T₀) proportional to the perturbations in the temperature,

Dvi- 1 ∂p-- ∂ τij ---- = -------+ ν∇2u ′ ----- + fi; i = 1,2,3 (5.1) Dt ρ ∂xi ∂xj

this body force per unit volume, applied to a fluid element, has the potential of modifying the turbulent structure through a production (or destruction) of Reynolds stress at a rate proportional to (f_i′u′_j + f_j′u′_i ), where f_i′ is the turbulent perturbation superimposed on f_i . The effects of this body force will change according to the physical processes in question, i.e. through the coupling with the flow or scalar fields.

Chapter 6
Scale Resolving Strategies: LES and V-LES

6.1 The filtered Navier-Stokes equations (bases of LES)

In Large Eddy Simulation (LES) the motion of the super-grid turbulent eddies is directly captured whereas the effect of the smaller scale eddies is modelled or represented statistically by means of simple models, very much the same way as in Reynolds-averaged models (RANS); i.e. the usual practice is to model the sub-grid stress tensor by an eddy viscosity model. In terms of computational cost, LES (Sagaut, 2005) lies between RANS and DNS and is motivated by the limitations of each of these approaches. Since the large-scale unsteady motions are represented explicitly, LES is more accurate and reliable than RANS. In the present work, use was made to the Dynamic sub-grid scale model of (Moin et al., 1991); both in the single-phase and boiling flow case. LES involves the use of a spatial filtering operation

∫ --- +∞ ′ ′ ′ Fk (x,t) = -∞ F (x,t)G(x - x )dx (6.1)

where the fluctuation of any variable F from its filtered value is denoted by f^′ = F_k - F_k. Filter function G(x-x^′) is invariant in time and space, and is localized, obeying the properties:

G (x) = G (- x ) (6.2)

∫ +∞ -∞ G (x)dx = 1 (6.3)

Applying the filtering operation to the instantaneous Navier-Stokes equations under incompressible flow conditions leads to the system of filtered transport equations for turbulent convective flow:

∇ ⋅u-= 0 (6.4)

∂ -- ---- -- ∂t(ρu )+ ∇ ⋅(ρu u) = ∇ ⋅(π - τ) + Fb + Fc (6.5)

where u is the velocity vector, ρ is the density, π=(-p.I + σ) is the Cauchy stress embedding pressure and viscous terms. The source terms in the RHS of the momentum equation represents the body force, F_b, and the convolution-induced, F_c. Further, filtering introduces the SGS stress tensor defined as τ = ρ(uu- u u), of which the deviatoric part is to be statistically modeled. This way, turbulent scales larger than the filter width imposed by filter function G(x-x^′) are directly solved, whereas the diffusive effects of the SGS scales are modeled.

6.2 Sub-grid Scale (SGS) Modelling

In turbulent flows, small-scale eddies are known to be simpler to model than the entire spectrum, since, in this high-wave number zone turbulence is likely to be homogeneous and isotropic (Kolmogorov, 1942). In other words, the subgrid-scale turbulence is much less problem–dependent than the resolvable one. Use is generally made of the Eddy Viscosity Concept, linking linearly the SGS eddy viscosity and thermal diffusivity to the gradients of the filtered velocity and temperature, respectively:

--- 1 ----2 τij = - 2μsgsSij + -δijτll; μsgs = (Cs Δ )2ρ|S| (6.6) 3

-- q” = - α ∂T-; α = μsgs (6.7) j θ∂xj θ Prt

where Pr_t is the turbulent Prandtl number, the strain rate tensor S_ij and the temperature gradients are determined from the large-scale turbulence. In the above two relations, the eddy viscosity μ_t and thermal diffusivity α_θ parameterise locally the non-resolvable dynamic stresses and the heat flux in terms of the local rate of strain tensor S_ij and the temperature gradients, calculated from the large-scale turbulence. The model constant (Cs) is either fixed or made dependent on the flow; this latter option is precisely the spirit of the dynamic model, known as DSM (Moin et al., 1991). A damping function is often introduced for the model constant Cs to accommodate the asymptotic behavior of near-wall turbulence. Similarly, the same strategy could be used to close the turbulent SGS heat flux, where the thermal diffusivity could be determined either based on the resolved thermal-flow field, or alternatively based on the eddy viscosity (defined dynamically) and a fixed Pr_t. In the present simulations, we have tested two variants based on fixed and variable model coefficient Cs, namely the DSM and the WALE sgs models (Nicoud & Ducros, 1999), the latter has been shown to behave very well in wall-bounded flows, and to be less dissipative, capable to capture the thin-shear layer accurately.

6.2.1 The Smagorinsky Model

The Smagorinsky model (Smagorinsky, 1963) is based on the Boussinesq approach, Eqs. (6.6 and 6.7). It is nothing more than a mixing-length theory

|| -|| -- -- -- -- -- -- νt = ℓ2 ||dv|| = ℓ2 | S | = (Cs Δ )2 |S | with: | S |= (2SijSij)1∕2 (6.8) dy

where the length scale ℓ is based on the cell size Δ³ = (Δ₁Δ₂Δ₃) rather than on the distance to the wall, i.e. on the mixing length ℓ₀ = κ ⋅ y_n . The constant C_s is referred to as the Smagorinsky constant.

Other options for near-wall treatment are also available. The Deardorff model is given as (Deardorff, 1970)

-- ℓ = min(Cs Δ,κdw )

(6.9)

where d_w is the minimum distance to the wall. In addition to this the Schmidt and Schumann model (Schmidt & Schumann, 1989) is given as,

-- ℓ = m(in (CsΔ, ℓmix)) (6.10) -1-- -1-- -1 ℓmix = c Δ-+ κdw (6.11) 2

where c₂ is taken to be c₂ = 0.1.

Assuming a direct analogy between momentum and heat flux diffusion by turbulence, one can determine the thermal diffusivity α_θ by reference to the eddy viscosity given by α_θ = ν_t∕Pr_t , as is usually done in turbulence modelling. This practice is quite robust in the context of the Smagorinsky model.

6.2.2 The Dynamic Approach (DSM)

In this approach developed by Germano et al. (1991), A priori adjustment of the Smagorinsky constant C_s is not needed. Instead, the dynamic model is used to determine this unknown variable C_s from the information contained in the resolved velocity field.

The main idea consists of introducing a test filter with a larger width than the original one, i.e. -^
Δ > Δ. This test filter is then applied to the filtered Navier-Stokes equations (the NS equations are filtered twice), yielding a sub–test–scale stress tensor T_ij similar in form to τ_ij that takes the following form:

---- ----- Tij = ^u′iu′j - ^uiu^j (6.12)

By virtue of the Germano identity, the two SGS stress tensors τ_ij and T_ij are connected through the following relation

----- Tij - ^τij ≡ Lij = u′iu′j - u^i ^uj- (6.13)

Assuming now the Smagorinsky functional form to hold and a variable coefficient C_s to be used to close the deviatoric parts of τ_ij and T_ij , we get

τ = - 2(C Δ)2 | S-| S (6.14) ij s ij

-- -- -- Tij = - 2(Cs ^Δ)2 | ^S | ^Sij (6.15)

where

= (

_i,j +

_j,i)∕2 is the rate–of–strain tensor of the test–filtered velocity field. Rearranging the last three equations results in

⌊--2 ⌋ -- 2 Δ^ ^- ^- -- -- Lij = - 2(CsΔ ) ⌈--2 | S |Sij - |S |Sij⌉ (6.16) Δ

A variant of this model (Lilly, 1992) uses a least–squares approach to obtain values for (C_s Δ)², leading to

-- 1 L M (CsΔ )2 = - ----ij---ij- (6.17) 2 MijMij

where

_ij is the term in brackets in Eq. (6.16).

Following the procedure employed to derive the dynamic parameter C_s , we can arrive at a relation for C_θ , the equivalent parameter for the thermal diffusivity, through

--2 α θ = (C θΔ ) | S | (6.18)

We get

-- LθjM θj (Cθ Δ)2 = - ---θ---θ- (6.19) M jM j

where

⌊--2 -- --⌋ θ -- 2 ^Δ ^- ^∂T -- ∂ T L j = - (C θΔ ) ⌈--2 |S | ∂x--- |S |∂x--⌉ (6.20) Δ j j

6.2.3 Turbulent Prandtl Number

The thermal diffusivity can be determined within the Newtonian closure framework using α_θ = ν_t∕Pr_t , but with the turbulent Prandtl number Pr_t , not α_θ , derived via the dynamic procedure of Germano (Moin et al., 1991). Here as well, information on the resolvable flow and temperature fields serve to determine the dynamic Pr_t . The turbulent Prandtl number takes the following form:

⌊ --2 -- -- ⌋ --2 ⌈ ^Δ-- ^- -^∂T- ∂T-- -- -- ⌉ Prt = (Cs Δ ) Δ2 |S |∂xj ∂xj - |S ||T | (6.21)

where

-- -- -- ∂T ∂T |T |= ∂x--∂x-- (6.22) j j

and (C_s Δ)² is determined from Eq. (6.17).

Alternatively, a constant turbulent Prandtl number can be assumed, typically (Pr_t ≈ 0.85). The turbulent Prandtl number can be also calculated based on the model proposed by Kays (1994) given as,

0.7 ν P rt = 0.85+ P-r ν t

(6.23)

where Pr = μC_p∕λ is the molecular Prandtl number, nu is the kinematic viscosity. As the grid is refined Pr_t tends towards large values. This is acceptable since the eddy thermal diffusivity tends to zero. However, it can be shown that the thermal diffusivity tends to zero much faster than the eddy kinematic viscosity.

6.2.4 The WALE Model

In the WALE concept (Nicoud & Ducros, 1999), the eddy viscosity is given as

(SdSd )3∕2 νt = (Cw Δ )2-¯-¯--5ij∕2--ij--d--d-5∕4- (SijSij) + (SijS ij)

(6.24)

where, C_w = ∘ -------
10.6 C2s and _ij^d is defined as,

d 1 2 2 1 2 Sij = 2(¯gij + ¯gji)- 3δij¯gkk (6.25) ∂u¯ ¯gij = ---i (6.26) ∂xj

where, g_ij² = g_ikg_kj and δ_ij is the Kronecker symbol.

6.3 Variable Turbulence Prandtl Number

Turbulent Prandtl number is given by Kays (Churchill, 2002)

( ) 0.7- ν- P rt = 0.85 + Pr νt (6.27)

where ν_t is limited by

ν = max (ν , ν-) (6.28) t ( t F ) 85- 0.85 F = --0.7---- P r (6.29)

6.4 Wall treatment in LES

At high Reynolds number, resolving eddies in the near-wall region would require the use of a very fine mesh close to the wall. To avoid this, Near-wall treatment is available in TransAT.

6.4.1 The Werner-Wengle model

The Werner-Wengle wall law (Werner & Wengle, 1991) consists in a two-layer approximation, with an assumption that a 1∕7^th power law holds outside the viscous sub-layer. It can then be written

{ + + u+ = y if y ≤ 11.8 8.3(y+ )1∕7 if y+ > 11.8

(6.30)

with

+ yuτ- y = ν (6.31) + -u- u = uτ (6.32)

being the dimensionless wall distance and dimensionless velocity, respectively. ν is the kinematic viscosity and u_τ = ∘ -----
τw∕ρ

is the wall shear stress velocity, with τ the shear at the wall and ρ the density. The wall shear stress is then computed using the instantaneous velocity in the flow.

6.5 Very-Large Eddy Simulation Concept (V-LES)

Very Large Eddy Simulation (V-LES) is based on the concept of filtering a larger part of turbulent fluctuations as compared to LES. This necessitates the use of a more complex sub–grid modelling strategy. The V-LES used in TransAT based on the use of the k - ϵ model as a sub–filter model. The filter width can be chosen by the user; it must be larger than the grid size. Increasing the filter width beyond the largest length scales in the flow, will lead to a standard RANS simulation, whereas in the limit of a small filter width the model will tend towards a large eddy simulation (although with a different sub–grid scale model).

The V-LES concept is based on the k -ε model. A filter is applied to this model, so that turbulent structures smaller than the filter width are not solved. A length-scale limiting function f has been derived (Johansen et al., 2004) and can be written

( ) ( ) Δ-ε- -Δε- f C3 k3∕2 = Min 1,C3 k3∕2

(6.33)

where Δ is the filter width and

C3 = ----γ∘-----= 1.0 4C μ 3∕2

(6.34)

is a new constant defined by Johansen et al. (2004), with γ ≤ 1 the anisotropic factor and C_μ = 0.09 the k -ε constant. Applying this function to the k -ε model gives the following expression for turbulent viscosity

[ Δε ] k2 νt = CμMin 1,C3k3∕2 ε--

(6.35)

Near the wall boundaries, the length-scale limiting function is equal to one, which means that the standard k - ε model is applied. This enables one to use standard wall-function of RANS models for V-LES.

Bibliography

Abe, K., Kondoh, T. & Nagano, Y. 1994 A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows–ii. the thermal field calculations. Int. J. Heat Mass Transfer 37, 139–151.

Bremhorst, K. & Bullock, K. 1973 Spectral measurements of turbulent heat and momentum transfer in fully developed pipe flow. International Journal of Heat and Mass Transfer 16 (12), 2141–2154.

Carteciano, L. N. & Grötzbach, G. 2003 Validation of turbulence models in the computer code flutan for a free hot sodium jet in different buoyancy flow regimes. Tech. Rep. FZKA-6600. Forschungszentrum Karlsruhe GmbH, Karlsruhe, Institut fr Kern- und Energietechnik.

Churchill, S. W. 2002 A reinterpretation of the turbulent prandtl number. Ind. Eng. Chem. Res. 41, 6393–6401.

Daly, B. & Harlow, F. 1970 Transport equations in turbulence. Phys. Fluids, A 13, 2634–2649.

Deardorff, J. 1970 A numerical study of three-dimensional turbulent channel flow at large reynolds numbers. J. Fluid Mech. 41, 453–480.

Durbin, P. 1991 Near-wall turbulence closure without damping function. Theoretical & Comput. Fluid Dynamics 3, 1–13.

Durbin, P. 1993 A reynolds stress model for near-wall turbulence. J. Fluid Mech. 249, 465–498.

Germano, M., Piomelli, U., Moin, P. & Cabot, W. 1991 A dynamic subgrid–scale eddy viscosity model. Phys. Fluids A 3, 1760–1765.

Hinze, J. O. 1975 Turbulence. New York, USA: MacGraw-Hill.

Johansen, S. T., Wu, J. & Shyy, W. 2004 Filtered-based unsteady RANS computations. Int. J. of Heat and Fluid Flow 25, 10–21.

Jones, W. & Launder, B. 1972 Prediction of relaminarization with a two-equation turbulence model. Int. J. Heat Mass Transfer 15, 301–314.

Kays, W. 1994 Turbulent prandtl number - where are we? In The 1992 Max Jacob Memorial Lecture. ASME J. Heat Transfer, , vol. 116, pp. 284–295.

Kenjereš, S., Gunarjo, S. & Hanjalić, K. 2005 Contribution to elliptic relaxation modelling of turbulent natural and mixed convection. International Journal of Heat and Fluid Flow 26 (4), 569–586.

Kolmogorov, A. 1942 The equations of turbulent motion in an incompressible fluid. Ivz. Acad. Sci. USSR, Phys. 6, 56–58.

Lam, C. K. G. & Bremhorst, K. 1981 A modified form of the k - ε model for predicting wall turbulence. Journal of Fluids Engineering 103, 456–460.

Launder, B. 1988 On the computation of convective heat transfer in complex turbulent flows. ASME J. Heat Transfer 110, 1112–1128.

Launder, B. & Spalding, D. 1974 The numerical computation of turbulent flows. Comput. Meths. Appl. Mech. Eng. 3, 269–289.

Launder, B. E. & Sharma, B. I. 1974 Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc. Letters in Heat and Mass Transfer 1, 131–138.

Lilly, D. 1992 A proposed modification of the germano subgrid–scale closure method. Phys. Fluids A 4, 633–635.

Moin, P., Squires, K., Cabot, W. & Lee, S. 1991 A dynamic subgrid-scale model for compressible turbulence and scalar transport. Phys. Fluids A 3, 2746–2757.

N. Chidambaram, S. D. & Kenzakowski, D. 2001 Scalar variance transport in the turbulence modeling of propulsive jets. JPP 17.

Nicoud, F. & Ducros, F. 1999 Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbulence and Combustion 62, 183–200.

Norris, L. & Reynolds, W. 1975 Turbulent channel flow with a moving wavy boundary. In Rept. No. FM-10. Stanford University, Dept. Mech. Eng.

Rodi, W. 1991 Experience with two–layer models combining the k - ε model with a one–equation model near the wall. AIAApaper 91–0216.

Rodi, W. & Mansour, N. 1993 Low reynolds k-ε modelling with the aid of direct numerical simulation. J. Fluid Mech. 250, 509–529.

S. Kenjeres, K. H. 2000 Convective rolls and heat transfer in finite-lenght rayleigh-bnard convection: A two-dimensional numerical study. Tech. Rep.. Department of Applied Physics, Thermo-Fluids Section, Delft University of Technology, The Netherlands.

Sagaut, P. 2005 Large Eddy Simulation for Incompressible Flows: An Introduction. Springer.

Sarkar, S., Erlebacher, G., Hussaini, M. Y. & Kreiss, H. O. 1991 The analysis and modelling of dilatational terms in compressible turbulence. Journal of Fluid Mechanics 227, 473–493.

Schmidt, H. & Schumann, U. 1989 Coherent structure of the convective boundary layer derived from large-eddy simulation. J. Fluid Mech. 200, 511–562.

Shams, A., Roelofs, F., Baglietto, E., Lardeau, S. & Kenjeres, S. 2014 Assessment and calibration of an algebraic turbulent heat flux model for low-prandtl fluids. International Journal of Heat and Mass Transfer 79, 589 – 601.

Shih, T.-H., Zhu, J., Liou, W., Chen, K.-H., Liu, N.-S. & Lumley, J. L. 1997 Modeling of turbulent swirling flows .

Smagorinsky, J. 1963 General circulation experiments with the primitive equations, i, the basic experiment. Mon. Weather Rev. 91, 99–165.

T.P. Sommer, R. S. & Lai, Y. 1993 Near-wall variable-prandtl-number turbulence model for compressible flows. AIAA Journal 31 (1), 27–35.

Werner, H. & Wengle, H. 1991 Large-eddy simulation of turbulent flow over and around a cube in a plate channel. In 8th Symposium on Turbulent Shear Flows. Munich, Germany.

Wilcox, D. 1993 Turbulence Modeling for CFD. DCW Ind., Inc., La Cañada, California.

Yap, C. R. 1987 Turbulent heat and momentum transfer in recirculating and impinging flows. Doctoral thesis, Univ. Manchester, United Kingdom.

Overview