The k-ε turbulence model proposed by Harlow and Nakayama in 1968 is by far the most widely-used two-equation eddy-viscosity-using turbulence model. The popularity of the model, and its wide use and testing, have thrown light on both its capabilities and its short-comings, which are well-documented in the literature.

PHOENICS provides the standard high-Reynolds-number form of the k-ε model, as presented by Launder and Spalding [1974], with inclusion of allowance for buoyancy effects.

Also provided are:-

- the modified model of Yap (see Section 3.4.4 below on the Lam-Bremhorst k-ε model);
- a low-Reynolds extension with an option for the so-called 'Yap'
- correction for separated flows (see Section 3.4.4 below);
- the two-scale split-spectrum extension (see Section 3.5 below);
- the modified model of Chen and Kim (see Section 3.4.2 below);
- the RNG-derived model (see Section 3.4.3 below); and
- the Realisable model (see Section 3.4.8 below).

For high turbulent Reynolds numbers, the standard form of the k-ε model may be summarised as follows:

∂/∂t (ρ*k) + **∇.**(ρ***u***k)= **∇.**(ρ*{ν_{l}+ν_{t}/σ_{t,k}}***∇** k )+ ρ*(P_{k} + G_{b} - ε)

∂/∂t (ρ*ε) + **∇.**(ρ***u***ε)=**∇.**(ρ*{ν_{l}+ν_{t}/σ_{t,ε}}***∇** ε ) + {ρ*ε/k}*(C_{1ε}*P_{k} + C_{3ε}*G_{b} -
*C*_{2ε}*ε)

ν_{t} = ℂ_{μ}*k^{2}/ε

P_{k} = ν_{t} * (∂u_{i}/∂x_{j} +∂u_{j}/∂x_{i} ) ∂u_{i}/∂x_{j}

G_{b} = - ν_{t} * **g*****∇**ρ/(ρ*σ_{t,h})

(1)

(2)

(3)

(4)

(5)

where k is the turbulent kinetic energy; ε is the dissipation rate; ρ is the
fluid density; ν_{l} and ν_{t} are the laminar and turbulent
kinematic viscosities, respectively;
P_{k} is the volumetric production rate of k by shear forces;
G_{b} is the volumetric production rate of k by gravitational
forces interacting with density gradients; **g** is the gravitational vector;
and σ_{t,h} is the turbulent Prandtl number.

G_{b} is negative for stably-stratified (dense below light) layers,
so that k is reduced and turbulence damped. For unstably-stratified
(dense above light) layers, G_{b} is positive in which therefore
k increases at the expense of gravitational potential energy. With the
Boussinesq approximation, in which the variations in density are expressed by
way of variations in enthalpy, eqn (5) reduces to:

G_{b} = ν_{t}*{β/(c_{p} * σ_{t,h})}***g*****∇**h

(6)

where *h*_{1} is the enthalpy, c_{p} is the specific heat
at constant pressure, and β is the volumetric coefficient of expansion.

If the energy equation is solved in terms of the temperature *T*, rather
than enthalpy *h*, eqn (6) is replaced by:

G_{b} = ν_{t}*{β/σ_{t,T}}***g*****∇**T

(7)

The following constants are normally used:

ℂ_{μ}=C_{μ}C_{d}=0.09, σ_{t,k}=1.0, σ_{t,ε}=1.314,
C_{1ε}=1.44, C_{2ε}=1.92, C_{3ε}=1.0.

The constant C_{3ε} has been found to depend on the flow situation. It should be close to
zero for stably-stratified flow, and to 1.0 for unstably-stratified flows. The default in
the PHOENICS VR Menu is to compute C_{3ε} from the function proposed by Henkes* et al* [1991]:

C_{3ε} = tanh (*v/u*)

(8)

where *v* and *u *are, respectively, the velocity components parallel
and perpendicular to the gravity vector. This function arranges that C_{3ε}
is unity when the main flow direction is aligned with gravity, and zero when the
main flow direction is perpendicular to gravity. The computed value of
C_{3ε} can be stored by using the command STORE(C3EB) in the Q1
file,as can G_{b} and the density/temperature gradients via
STORE(GENB,DRDX,DRDY,DRDZ).

If PHOENICS VR is not used, then a default value of C_{3ε}=1.0
is employed *via* the setting GCEB=1.0 in the GROUND subroutine GXKEGB (in
the file GXGENK.FOR). A value of C_{3ε}=0.0 may be effected by
simply not using a COVAL statement for ε. The default value of GCEB=1.0
may be overwritten by setting, for example:

SPEDAT(SET,KEBUOY,GCEB,R,0.2)

in the Q1 file. The setting SPEDAT(SET,KEBUOY,GCEB,R,-1.0) arranges that C_{3ε} is computed
from eqn (8) above.

The model presented above is applicable only in regions where the turbulence Reynolds number is high. Near walls, where the Reynolds number tends to zero, the model requires the application of so-called 'wall functions' (see Section 8 below ), or alternatively, the introduction of a low-Reynolds-number extension (see Sections 3.3.1 and 3.4.4. The standard model employs the wall-function approach.

It should be mentioned that the standard form of the k-ε model has been found to perform less than satisfactorily in a number of flow situations, including:

- separated flows,
- buoyancy,
- streamline curvature,
- swirl,
- turbulence-driven secondary motions,
- rotation,
- compressibility,
- adverse pressure gradients and
- axi-symmetrical jets.

Nevertheless because the model is so widely used, variants and/or *ad hoc* modifications
aimed at improving its performance abound in the literature. The most well-known include:

(a) the Realisable, RNG, Chen-Kim and Yap variants for use in separated flows; and

(b) the *ad hoc* Richardson-number modification of Launder et al [1977] for curvature,
swirl and rotation.

The standard k-ε model is activated by inserting the PIL command TURMOD(KEMODL) in the Q1 file, which is equivalent to:

SOLVE(KE,EP);ENUT=GRND3;EL1=GRND4;KELIN=0 PATCH(KESOURCE,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP) COVAL(KESOURCE,KE,GRND4,GRND4) COVAL(KESOURCE,EP,GRND4,GRND4);GENK=T TERMS(KE,N,Y,Y,Y,Y,N);TERMS(EP,N,Y,Y,Y,Y,N)

The sources for k and ε are calculated and inserted in subroutine GXKESO called from GROUP 13 of GREX. Different linearizations of these sources are selected by the KELIN parameter. The generation rate used in the source terms can be stored by the command STORE(GENK).

The WALL and CONPOR commands create the required wall-function COVAL settings automatically,

COVAL(WALLN,KE,GRND2,GRND2); COVAL(WALLN,EP,GRND2,GRND2)

or

COVAL(WALLN,KE,GRND3,GRND3); COVAL(WALLN,EP,GRND3,GRND3)

if the user sets WALLCO=GRND3 in the Q1 file.

Thus, the PHOENICS default is equilibrium (GRND2) wall functions. However, in separated
flows the use of non-eq*u*_{i}librium (GRND3) wall functions is recommended, especially if wall
heat transfer is present.

The buoyancy source terms are activated in the k and ε equations when the user introduces the following PATCH and COVAL statements in the Q1 file:

PATCH(KEBUOY,PHASEM,1,NX,1,NY,1,NZ,1,1) COVAL(KEBUOY,KE,GRND4,GRND4) COVAL(KEBUOY,EP,GRND4,GRND4)

The Cartesian components of the gravitational vector, which appears in equation (5), are defined by setting BUOYA, BUOYB and BUOYC in the Q1 file.

These parameters must be set in connexion with the buoyancy source term in the momentum equations. For example, with gravity acting in the negative y-direction in Cartesian coordinates, the following PIL sequence introduces a buoyancy source in the V1 equation:

BUOYA=0.;BUOYB=-9.81;BUOYC=0.0;BUOYD=0.5*(ρB+ρT) PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,1) COVAL(BUOY,V1,FIXFLU,GRND2)

Here BUOYD is the reference density defined by the user.

With the Boussinesq approximation, and solution of the energy equation *via* *h*_{1}, the
buoyancy force would be activated as follows:

BUOYA=0.;BUOYB=-9.81;BUOYC=0.0;BUOYD=....;BUOYE=.... PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,1) COVAL(BUOY,V1,FIXFLU,GRND3)

Here BUOYD = -β/*c*_{p} and BUOYE = -BUOYD**h*_{1},ref and *h*_{1},ref
is the reference enthalpy.

If the energy equation is solved *via* the *T*_{1} variable,
then

BUOYD = -β and BUOYE = *T*_{1}_{ref}

where *T*_{1}_{ref}
is the reference temperature; if the volumetric coefficient of expansion is
variable, then BUOYD=GRND10 is appropriate.

The Boussinesq form of the k-ε buoyancy source terms is used only when the density is
set as a constant *via* the Q1 or the menu, and not when a constant value is set *via* GREX or
GROUND.

For separated flows, the so-called Yap-correction to the ε eqn can be introduced into the k-ε model, as follows:

TURMOD(KEMODL-YAP) for the high-Re k-ε model; and TURMOD(KEMODL-LOWRE-YAP) for the low-Re form of the model.

The extension -YAP is equivalent to the following PIL commands:

DISWAL;PATCH(EYAP,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP) COVAL(EYAP,EP,FIXFLU,GRND5) FIINIT(WDIS)=0.1*AMIN1(XULAST,YVLAST,ZWLAST)

TURMOD(KECHEN) selects the modified k-ε model of Chen and Kim, which is described below in Section 3.4.2.

TURMOD(KERNG) selects the RNG-derived k-ε model of Yakhot and Orszag, which is described below in Section 3.4.3.

TURMOD(KEMODL-LOWRE) selects the Lam-Bremhorst low-Re extension to k-ε model, which is described below in Section 3.4.4.

TURMOD(KEMODL-2L) activates the two-layer low-Re model, as an alternative to the conventional low-Re k-ε models. The two-layer model employs the high-Re k-e model away from the wall, and a one- equation in the near-wall region.

TURMOD(KEREAL) selects the realisable k-ε model, which is described below in Section 3.4.8.

- F.H.Harlow and P.I.Nakayama, 'Transport of turbulence energy decay rate', LA-3854, Los Alamos Science Lab., U. California, USA, (1968).
- R.A.Henkes,F.van der Flugt and C.Hoogendoorn, 'Natural convection in a square cavity with low-Re turbulent fluids', Int.J.Heat Mass Transfer, Vol.34, p 1543-1557, (1991).
- B.E.Launder and D.B.Spalding, 'The numerical computation of turbulent flows', Comp. Meth. in Appl. Mech. & Eng., Vol.3, pp 269, (1974).
- B.E.Launder, C.H.Priddin and B.R.Sharma, 'The calculation of turbulent boundary layers on spinning and curved surfaces', ASME J Fluids Engng., Vol.99, p 321, (1977).
- W.Rodi, 'Examples of calculation methods for flow and mixing in stratified fluids',J.Geo.Res., Vol.92, No.C5, p 5305, (1987).

See also the PHOENICS-Instruction-Course lectures on Turbulence Modelling.