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:-
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 )+ ρ*(Pk + Gb - ε)
∂/∂t (ρ*ε) + ∇.(ρ*u*ε)=∇.(ρ*{νl+νt/σt,ε}*∇ ε ) + {ρ*ε/k}*(C1ε*Pk + C3ε*Gb - C2ε*ε)
νt = ℂμ*k2/ε
Pk = νt * (∂ui/∂xj +∂uj/∂xi ) ∂ui/∂xj
Gb = - ν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; Pk is the volumetric production rate of k by shear forces; Gb 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.
Gb is negative for stably-stratified (dense below light) layers, so that k is reduced and turbulence damped. For unstably-stratified (dense above light) layers, Gb 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:
Gb = νt*{β/(cp * σt,h)}*g*∇h
(6)
where h1 is the enthalpy, cp 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:
Gb = νt*{β/σt,T}*g*∇T
(7)
The following constants are normally used:
ℂμ=CμCd=0.09, σt,k=1.0, σt,ε=1.314, 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 C3ε 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,as can Gb and the density/temperature gradients via STORE(GENB,DRDX,DRDY,DRDZ).
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 ε. 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 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:
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-equilibrium (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 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 = -β/cp and BUOYE = -BUOYD*h1,ref and h1,ref is the reference enthalpy.
If the energy equation is solved via the T1 variable, then
BUOYD = -β and BUOYE = T1ref
where T1ref 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.
See also the PHOENICS-Instruction-Course lectures on Turbulence Modelling.