The KE-EP turbulence model proposed by Harlow and Nakayama in 1968 is by far the most widely-used two-equation eddy-viscosity turbulence model, mainly because the EP was long believed to require no extra terms near walls. The popularity of the model, and its wide use and testing, has 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 KE-EP 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 KE-EP 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 KE-EP model may be summarised as follows, with ,t denoting differentiation with respect to time and ,i denoting differentiation with respect to distance:

(RHO*KE),t + (RHO*Ui*KE - {RHO*ENUT/PRT(KE)}*KE,i ),i = RHO*(Pk + Gb - EP) (1)

(RHO*EP),t + (RHO*Ui*EP - {RHO*ENUT/PRT(EP)}*EP,i ),i = {RHO*EP/KE}*(C1*Pk + C3*Gb - C2*EP) (2)

ENUT = CMUCD*KE**2/EP (3)

Here KE is the turbulent kinetic energy; EP is the dissipation rate; RHO is the fluid density; ENUT is the turbulent kinematic viscosity.

Pk is the volumetric production rate of KE by shear forces. It is calculated from:

Pk = ENUT * (Ui,j + Uj,i) Ui,j (4)

Gb is the volumetric production rate of KE by gravitational forces interacting with density gradients. It is calculated from:

Gb = - ENUT * gi * {RHO,i}/(RHO * PRT(H1)) (5)

where gi is the gravitational vector and PRT(H1) is the turbulent Prandtl number.

Gb is negative for stably-stratified (dense below light) layers, so that KE is reduced and turbulence damped.

Gb is positive for unstably-stratified (dense above light) layers, in which therefore KE 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:

Gb = ENUT * BETA * gi * {H1,i}/(CP * PRT(H1)) (6)

where H1 is the enthalpy, CP is the specific heat at constant pressure, and BETA is the volumetric coefficient of expansion.

If the energy equation is solved in terms of the temperature TEM1, rather than enthalpy H1, equation (6) is replaced by:

Gb = ENUT * BETA * gi * {TEM1,i}/PRT(TEM1) (7)

The following constants are normally used:

PRT(KE)=1.0, PRT(EP)=1.314, CMUCD=0.09, C1=1.44, C2=1.92, C3=1.0.

The constant C3 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 C3 from the function proposed by Henkes et al [1991]:

C3 = tanh (v/U) (8)

where v and U are, respectively, the velocity components parallel and perpendicular to the gravity vector. This function arranges that C3e 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 C3 can be stored by using the command STORE(C3EB) in the Q1 file.

If PHOENICS VR is not used, then a default value of C3=1.0 is employed via the setting GCEB=1.0 in the GROUND subroutine GXKEGB (in the file GXGENK.FOR). A value of C3=0.0 may be effected by simply not using a COVAL statement for EP. 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 C3 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 KE-EP 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 KE-EP 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 KE and EP 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-equilibrium (GRND3) wall functions is recommended, especially if wall heat transfer is present.

The buoyancy source terms are activated in the KE and EP 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*(RHOB+RHOT) 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 H1, 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 = -BETA/CP and BUOYE = -BUOYD*H1,ref and H1,ref is the reference enthalpy.

If the energy equation is solved via the TEM1 variable, then

BUOYD = -BETA and BUOYE = TEM1,ref

where TEM1 is the reference temperature; if the volumetric coefficient of expansion is variable, then BUOYD=GRND10 is appropriate.

The Boussinesq form of the KE-EP 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 EP eqn can be introduced into the KE-EP model, as follows:

TURMOD(KEMODL-YAP) for the high-Re KE-EP 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 KE-EP model of Chen and Kim, which is described below in Section 3.4.2.

TURMOD(KERNG) selects the RNG-derived KE-EP 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 KE-EP 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-e 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, p1543-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, pp269, (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, p321, (1977).

W.Rodi, 'Examples of calculation methods for flow and mixing in stratified fluids', J.Geo.Res., Vol.92, No.C5, p5305, (1987).

See also the Instruction Course lectures on Turbulence Modelling.