Encyclopaedia Index
Back to start of turbulence-model article

3.4.4 Lam Bremhorst k-ε turbulence model

Contents

  1. Description
  2. Activation of the model
  3. Advice on the use of the model
  4. Sources of further information

1. Description

The Lam-Bremhorst ( hereafter denoted LB ) low-Reynolds-number extension to the k-ε model employs a transport equation for the total dissipation rate, with the advantage that the model requires no additional source terms.

However, a disadvantage of the model is that one of the damping functions requires the calculation of the local distance to the nearest wall.

The form of the LB model implemented in PHOENICS is that described by Patel et al [1984], although it should be noted that since then various slightly-modified versions have been proposed by a number of workers ( see for example Davidson [1990] and Herrero et al [1991]).

The LB low-Reynolds-number k-ε model differs from the standard high-Reynolds-number model in that the empirical coefficients CμCD, C and C are multiplied respectively by the functions:

Fμ = [ 1.-exp(-0.0165*Ren) ]2 (1.+20.5/Ret)

F1 = 1.+(0.05/Fμ)3

F2 = 1.-exp(-Ret**2)

Ren = √k*δnl

Ret = k2/ε/νl

(1)

(2)

(3)

(4)

(5)

 

where δn is the distance to the nearest wall. For high-turbulence Reynolds numbers Ren or Ret, the functions Fμ, F1 and F2 multiplying the three constants tend to unity.

The boundary conditions k=0 and ∂ε/∂y=0 are applied at the wall.

The present implementation of the model has no restriction on its functionality within PHOENICS, although it does not allow simulation of the final-stage decay of isotropic grid turbulence. However, the reader is referred to Davidson [1991] for a fairly simple modification which permits it to do so.

2. Activation of the model

(a) Basic activation


The LB extension to the k-ε model is selected by:

TURMOD(KEMODL-LOWRE)

which is equivalent to TURMOD(KEMODL) plus the PIL commands IENUTA=3 and DISWAL. The DISWAL command activates the solution of a scalar variable LTLS, from which is deduced the minimum distance to the nearest wall δn ( see below, the PHENC entries 'DISWAL' and Section 3.1.2 above on the LVEL turbulence model).

Subsequent WALL (or CONPOR) commands will set boundary conditions for the appropriate velocities, LTLS and turbulent kinetic energy.

No COVAL is required for ε as the default boundary condition is zero normal flux.

Where needed, COVALs for wall PATCHes should take the following form:

    PATCH(WALLS,SWALL,1,1,1,1,1,NZ,1,1)
    COVAL(WALLS,W1,1.0,0.0); COVAL(WALLS,KE,1.0,0.0)
    COVAL(WALLS,LTLS,1.0,0.0)
    COVAL(WALLS,H1,1.0/PRNDTL(H1),HWALL)  or
    COVAL(WALLS,TEM1,1.0,TWALL)

so that laminar boundary conditions are set for the mean-flow variables, and a zero value of k is set at the wall.

When STORE(FMU,FONE,FTWO,REYT,REYN) appears in the Q1 file, the various damping functions and Reynolds numbers defined by equations (2.1) to (2.5) may be printed in the RESULT file or viewed via PHOTON and AUTOPLOT.

(b) Determination of the minimum wall distances

The model requires that the minimum distance to the nearest wall be determined for each cell in the flow field. These distances are calculated automatically via the solution of a differential equation for a generalised length-scale variable LTLS.

This calculation is activated by the TURMOD command, which activates also the command DISWAL.

See PHENC entry: distance from the wall.

(c) The 'Yap' correction for separated flows


In separated flows the standard ε equation gives far too large near-wall values of the turbulence length scale, especially if the equation is integrated to the wall with a low-Reynolds-number extension to the k-ε model. Consequences of this are that boundary-layer separation will tend to be predicted too late and, in separating and reattaching flow regions, excessive heat-transfer coefficients are predicted.

Yap [1987] proposed an additional source term to the ε equation which goes some way towards removing this deficiency. The volumetric source term takes the form:

S = max < 0.83*ρ*(L/Le-1)*{(L/Le)2}*(ε2)/k,0.>

(6)

where L = k3/2/ε and Le = CLn, with CL=κ/(CμCD)3/4 (=2.495) and δn is the distance from the wall.

The term vanishes in local-equilibrium wall turbulence because then L=Le; it also becomes small at large distances from the wall, since then L/Le << 1.

However, if L > Le the term is positive leading to increased values of ε and reduced values of L.

The Yap correction may be selected by using:

    TURMOD(KEMODL-LOWRE-YAP)

which is equivalent to:

    TURMOD(KEMODL-LOWRE)

plus:

    PATCH(EYAP,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
    COVAL(EYAP,EP,FIXFLU,GRND5)
    FIINIT(IWDIS)=0.1*AMIN1(XULAST,YVLAST,ZWLAST)
The Yap correction may also be applied to the high-Reynolds-number form of the k-ε model by using:
    TURMOD(KEMODEL-YAP).

(d) Other matters

The FORTRAN coding sequences for the LB model may be found in: Group 1 Sections 1 and 2, and in Group 19 Section 3 of subroutine GREX3; in subroutine GXENUT which resides in the file GXPROP.FOR; and in subroutines GXKESO, GXREYN, GXREYT and GXLRDF which reside in the file GXTURB.FOR.

The coding sequences for the Yap correction may be found in: Group 13 of subroutine GREX3 and in subroutine GXEYAP.

The convergence of the turbulence equations can be more problematic than with the standard high-Reynolds-number form of the k-ε model. In cases where convergence proves difficult the user is advised to set KELIN=1 which invokes an implicit increase in the relaxation on k and ε through a different linearisation of the turbulence-model source terms.

Finally, the low-Reynolds-number extension may be applied to the modified k-ε model of Chen and Kim by selecting TURMOD(KECHEN- LOWRE), which is equivalent to TURMOD(KEMODL) plus the following PIL commands:

          IENUTA=4;PRT(KE)=0.75;PRT(EP)=1.15;STORE(WDIS)
          PATCH(KECHEN,PHASEM,1,NX,1,NY,1,NZ,1,1)
          COVAL(KECHEN,EP,FIXFLU,GRND4)
For more details, the reader is referred to the Encyclopaedia entry provided under 'CHEN-Kim modified k-ε turbulence model'.

(e) Activation with "conjugate heat transfer"

The PHOENICS default boundary condition for KE, EP, C1, C2, etc at the fluid-solid interface is finite diffusion flux normal to the wall.

For consistency with the standard PHOENICS default boundary condition this ought to be zero normal flux. In order to achieve the desired boundary conditions at the interface, the diffusive links through the wall are set to zero by use of the so-called Group 12 Q1 facility. Thus for example the PIL commands:

          PATCH(GP12DFNT,CELL,1,NX,NYGL,NYGL,1,NZ,1,1)
          COVAL(GP12DFNT,KE,0.0,0.0)
          COVAL(GP12DFNT,EP,0.0,0.0)
          COVAL(GP12DFNT,C1,0.0,0.0)

will zero the north-face diffusion fluxes for KE, EP and C1 at the locations defined by the PATCH limits.

(f) Examples

A number of Q1 files may be found in the advanced-turbulence-models library of low-Reynolds-number turbulence models which demonstrate the use of the model.

3. Advice on the use of the model

Since the low-Reynolds-number extension does not employ wall functions, and the flow field needs to be meshed into the laminar sublayer and down to the wall, the computer storage and run-time requirements for this approach are much greater than those of the wall-function approach.

In general, the grid normal to the main flow direction needs to be distributed so as to give a high concentration of grid cells near the wall, with the wall-adjacent node positioned at y+=1.0 or even less.

In any case, the user is advised that the location of the first grid point normal to the wall should not exceed Y+=4.0, and for reasonable accuracy it is recommended that at least five points are located in the region Y+ < 11.5. Here Y+ is defined by:

Y+ = Uτnl

(7)

 

where Uτ is the friction velocity at the wall (= √(τw/ρ)).

If the user wishes to know the near-wall values of y+ he may activate printout of their values in the PHOENICS RESULT file by setting YPLS=T.

Please note, however, that this facility requires that the user define the boundary condition on the velocity parallel to the wall via

COVAL(WALLS,W1,GRND2,0.0)

which is equivalent to setting a COefficient of unity provided that the near-wall value of y+ is less than 11.5. For "conjugate-heat- transfer" calculations one must set

COVAL(WALLS,W1,GRND2,SAME)

to elicit printout of the near-wall y+ values when YPLS=T. CONPOR- generated COVALs for wall patches produce CO=GRND2 by default.

While there are a number of choices for generating a non-uniform mesh across the flow ( such as, for example, the power-law option provided in PHOENICS ), one of the most economical options is to use a geometric progression with the property that the ratio of length of any two adjacent cross-stream intervals is a constant; i.e.

DY(J)=K*DY(J-1).

(8)

 

If, for example, we consider the case of a flow bounded by a single wall, the distance of the Jth cell face is given by the formula:

Y(J)=DY(1)*(K**J-1)/(K-1), K > 1.

(9)

 

There are two parameters in the above equation: DY(1), the length of the first grid cell, and K, the ratio of two successive cell sizes. The total number of cells NY can be calculated from the following formula:

NY = LN [1+(K-1)(Y(NY)/DY(1)]/LN[K].

(10)

 

DY(1) may be defined so that y+=1.0 by way of an estimated friction velocity Uτ.

Thus, for turbulent flow through a smooth pipe at a Reynolds number of 1.E5, DY(1) may be calculated as follows:

          REY=1.E5
          FRIC=1./(1.82*LOG10(REY)-1.64)**2
          US=UIN*(FRIC/8.)**0.5
          ENUL=WIN*DIAM/REY
          YPLS=1.0
          DY(1)=YPLS*ENUL/US
The value of k may be chosen typically as 1.1, and so approximately 60 grid cells across the flow are sufficient to represent the laminar and turbulent regions of the flow. The foregoing example can readily be extended to accommodate two walls. Geometric-progression grids created by PIL commands in the Q1 file are exemplified in several of the test cases provided in the advanced turbulence model library.

Some other useful formulae for the generation of non-uniform grid distributions are given by Mansole and Lage [1993].

4. Sources of further information

  1. Y.S.Chen and S.W.Kim, 'Computation of turbulent flows using an extended k-ε turbulence closure model', NASA CR-179204, (1987).
  2. L.Davidson, 'Calculation of the turbulent buoyancy-driven flow in a rectangular cavity using an efficient solver and two different low Reynolds number k-ε turbulence models', Num. Heat Transfer, Part A, Vol.18, p129, (1990).
  3. J.Herrero, F.X.Grau, J.Grifoll and F.Girault, 'A near-wall k-ε formulation for high Prandtl number heat transfer', Int.J.Heat Transfer, Vol.34, No.3, p711, (1991).
  4. C.K.G.Lam and K.Bremhorst, 'A modified form of the k-ε model for predicting wall turbulence', ASME J. Fluids Engng., Vol.103, p456,(1981).
  5. D.M.Mansole and J.L.Lage, ' Nonuniform grid accuracy test applied to natural-convection flow within a porous medium cavity.', Numerical Heat Transfer, Part B, Vol.23, p351, (1993).
  6. V.C.Patel, W.Rodi & G.Scheurer, 'Turbulence models for near-wall and low-Reynolds-number flows: A review', AIAA J, Vol.23, No.9, p1308, (1984)
  7. C.Yap,'Turbulent heat and momentum transfer in recriculating and impinging flows', PhD Thesis, Faculty of Technology, University of Manchester, (1987).