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); and
- the RNG-derived model (see Section 3.4.3 below).

For high turbulent Reynolds numbers, the standard form of the *k*-ε model may be
summarised as follows, with _{t} denoting differentiation with respect to time and _{i}
denoting differentiation with respect to distance:

(ρ**k*)_{t} + (ρ**u*_{i}**k* - {ρ*ν_{t}/*Pr*_{t}(*k*)}**k*_{i} )_{i} = ρ*(*P*_{k} + *G*_{b} - ε) (1)

(ρ*ε)_{t} + (ρ**u*_{i}*ε - {ρ*ν_{t}/*Pr*_{t}(ε)}*ε_{i} )_{i} = {ρ*ε/*k*}*(*C*_{1}**P*_{k} + *C*_{3}**G*_{b} -
*C*_{2}*ε) (2)

ν_{t} = CMUCD**k*^{2}/ε (3)

Here *k* is the turbulent kinetic energy; ε is the dissipation rate; ρ is the
fluid
density; ν_{t} is the turbulent kinematic viscosity.

*P*_{k} is the volumetric production rate of *k* by shear forces. It is calculated from:

*P*_{k} = ν_{t} * (*u*_{i},j + *u*_{ji}) *u*_{ji} (4)

*G*_{b} is the volumetric production rate of *k* by gravitational forces interacting with
density gradients. It is calculated from:

*G*_{b} = - ν_{t} * gi * {ρ_{i}}/(ρ * *Pr*_{t}(*h*_{1})) (5)

where gi is the gravitational vector and *Pr*_{t}(*h*_{1}) is the turbulent Prandtl number.

*G*_{b} is negative for stably-stratified (dense below light) layers, so that *k* is reduced
and turbulence damped.

*G*_{b} is positive for unstably-stratified (dense above light) layers, 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} * β * gi * {*h*_{1}_{i}}/(*c*_{p} * *Pr*_{t}(*h*_{1})) (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*_{1}, rather than enthalpy
*h*_{1}, equation (6) is replaced by:

*G*_{b} = ν_{t} * β * gi * {*T*_{1}_{i}}/*Pr*_{t}(*T*_{1}) (7)

The following constants are normally used:

*Pr*_{t}(*k*)=1.0, *Pr*_{t}(ε)=1.314, CMUCD=0.09, *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.

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 equation (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 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 eq*u*_{i}valent 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,k,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(*k*MODL-YAP) for the high-Re *k*-ε model; and TURMOD(*k*MODL-LOWRE-YAP) for the
low-Re form of the model.

The extension -YAP is eq*u*_{i}valent 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(*k*CHEN) selects the modified *k*-ε model of Chen and Kim, which is described
below in Section 3.4.2.

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

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

TURMOD(*k*MODL-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.

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).

Press, (1972).

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.