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, C1ε and C2ε are multiplied respectively by the functions:
F1 = 1.+(0.05/Fμ)3
F2 = 1.-exp(-Ret**2)
Ren = √k*δn/νl
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.
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.
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.
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:
(6)
where L = k3/2/ε and Le = CL*δn, 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).
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'.
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.
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.
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:
(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.
(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:
(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].