** PHOENICS VALIDATION CASE
    GROUP 1. Run title and other preliminaries
TEXT(K-OMEGA_1D DEVELOPED PIPE FLOW  :T207
TITLE
  DISPLAY
  The case considered is fully-developed turbulent flow and heat
  transfer in a circular pipe at Re=1.E5 and Pr=3.0. The tube wall
  is held at a constant temperature, and the calculation integrates
  down to the wall through the viscous sublayer. A non-uniform grid
  is employed so as to concentrate cells very close to the wall.
  For this purpose a grid is generated which is a geometric
  progression with the property that the ratio of any two adjacent
  cell lengths is a constant. The turbulence is simulated by one of
  4 low-Re turbulence models, namely: the Lam-Bremhorst k-e model;
  the Chen-Kim k-e model; the Wilcox k-omega model; or the 2-layer
  k-e model which employs the Norris-Reynolds 1-eqn model in the
  near-wall layer.
  ENDDIS
 
  The PHOENICS predictions reveal the following major results:
 
           Lam-Brem   K-Omega  Chen-Kim     2-Layer      Data
   f        0.0197     0.017    0.0180       0.0189     0.018
   Nu        393        352       365          410        392
 
  It should be noted that grid effects have not been investigated.
  The PIL variable WALPRN has been set to T, thereby activating
  printout of the local friction factor (sloc) and Stanton number
  (Stloc) in the RESULT file. In order to convert these values to
  f and Nu above, the following relations should be used:
 
    f = 8.*sloc*(W1(NY)/WIN)**2
 
    Nu = REY*PRNDTL(H1)*Stloc*W1(NY)*(TW-H1(NY))/[WIN*(TW-TB)]
 
  where TB is the bulk temperature printed in the RESULT file.
  ENDDIS
 
   The following AUTOPLOT use file produces four plots;
   the first is the axial velocity profile; the second
   is the temperature profile; the third is the turbulence
   energy profile; and the fourth is the dissipation rate
   profile.
 
   AUTOPLOT USE
   file
   phi 5
 
 
   da 1 w1;col9 1
   msg Velocity (W1) profile
   msg Press RETURN to continue
   pause
   clear
   da 1 tem1;col9 1
   msg temperature profile
   msg Press RETURN to continue
   pause
   clear
   da 1 ke;col9 1
   msg KE profile
   msg Press RETURN to continue
   pause
   clear
   da 1 ep;col9 1
   msg EP profile
   msg Press e to END
   ENDUSE
 
INTEGER(TMODEL);CHAR(CTURB,TLSC)
REAL(DIAM,WIN,REY,MIXL,FRIC,DPDZ,MASIN,DTF,TKEIN,EPSIN)
REAL(DELT1,DELT2,DELT3,US,EPSPLS,KFAC,DELY,AA,GR,GYPLUS)
REAL(GWPLUS,GWR,QIN,DTDZ,COND,CP,TW,AIN,AWAL,NUSS,XR,RIN)
DIAM=0.1;WIN=1.0; REY=1.E5;FRIC=1./(1.82*LOG10(REY)-1.64)**2
FRIC
DPDZ=FRIC*RHO1*WIN*WIN/(2.*DIAM);US=WIN*(FRIC/8.)**0.5
  ** estimate initial values from K+=2 & EP+=1./(ak*Y+)
EPSPLS=1./(0.41*100.); RIN=0.5*DIAM
TKEIN=2.*US*US;ENULA=WIN*DIAM/REY;EPSIN=EPSPLS*US**4/ENULA
  ** the grid-expansion factor Kfac defines a constant ratio of
     lengths of two adjacent cells.
MESG( Enter the required turbulence model:
MESG(  LAMB - Low-Re Lam-Brem. k-e model
MESG(  CHEN - Low-Re Chen-Kim  k-e model
MESG(  TWOL - Low-Re 2-layer   k-e model
MESG(  KO   - Low-Re k-omega model (default)
MESG(
READVDU(CTURB,CHAR,KO)
CASE :CTURB: OF
WHEN LAMB,4
+ TEXT(LAM-BRE KE_1D DEVELOPED PIPE FLOW :T207
+ MESG(Low-Re Lam-Bem. k-e turbulence model
+ TMODEL=1;KFAC=1.1;TLSC=EP
WHEN CHEN,4
+ TEXT(CHEN-KIM KE_1D DEVELOPED PIPE FLOW :T207
+ MESG(Low-Re Chen-Kim k-e turbulence model
+ TMODEL=2;KFAC=1.1;TLSC=EP
WHEN TWOL,4
+ TEXT(2-LAYER KE_1D DEVELOPED PIPE FLOW :T207
+ MESG(Low-Re 2-layer k-e turbulence model
+ SELREF=F;TMODEL=3;KFAC=1.2;TLSC=EP
WHEN KO,2
+ MESG(Low-Re k-omega turbulence model
+ TMODEL=4;KFAC=1.1;TLSC=OMEG
+ STORE(EP);EPSIN=EPSIN/(0.09*TKEIN)
ENDCASE
    GROUP 3. X-direction grid specification
CARTES=F;XULAST=0.1;AIN=0.5*RIN*RIN*XULAST
    GROUP 4. Y-direction grid specification
ENULA=WIN*DIAM/REY
  ** define first dely from wall
DELT1=1.*ENULA/US;DELY=DELT1/(0.5*DIAM)
  ** calculate NY from dely & Kfac
AA=(YVLAST/DELY)*(KFAC-1.0)+1.0;AA=LOG(AA)/LOG(KFAC)+1.0001
NY=AA
  ** define uniform grid initially
IREGY=1;GRDPWR(Y,NY,YVLAST,1.0)
  ** compute expanding grid from north boundary
YFRAC(NY)=1.0;INTEGER(JJM,JJM1)
DO JJ=NY,2,-1
+ JJM=JJ-1
+ YFRAC(JJM)=YFRAC(JJ)-DELY
+ DELY=KFAC*DELY
ENDDO
YVLAST=0.5*DIAM
    GROUP 5. Z-direction grid specification
ZWLAST=0.1*DIAM
    GROUP 7. Variables stored, solved & named
SOLVE(W1,TEM1);STORE(ENUT,LEN1,FMU,REYT,FTWO)
SOLUTN(W1,P,P,P,P,P,N)
CASE (TMODEL) OF
WHEN 1
+ TURMOD(KEMODL-LOWRE);KELIN=3;STORE(REYN,FONE)
WHEN 2
+ TURMOD(KECHEN-LOWRE);KELIN=3;STORE(REYN,FONE)
WHEN 3
+ TURMOD(KEMODL-2L);KELIN=1
WHEN 4
+ TURMOD(KOMODL-LOWRE)
ENDCASE
    GROUP 8. Terms (in differential equations) & devices
TERMS(W1,N,N,P,P,P,P);TERMS(KE,Y,N,P,P,P,P)
TERMS(:TLSC:,Y,N,P,P,P,P);TERMS(TEM1,N,N,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
ENUL=ENULA;PRNDTL(H1)=3.0;MASIN=RHO1*WIN*AIN
  ** prescribe energy flow from slab and fluid specific heat
     estimated from Dittus-Boelter Nu=0.023*Re**0.8*Pr**0.4
     with (Tw-Tb)=5.0
NUSS=0.023*REY**0.8*PRNDTL(H1)**0.4;CP=1.0;AWAL=RIN*XULAST
COND=RHO1*CP*ENUL/PRNDTL(H1);QIN=NUSS*5.0*COND/DIAM
NUSS
  ** compute d(tbulk)/dz for input to single-slab
     thermal solver and prescribe wall temperature
DTDZ=QIN*AWAL/MASIN;TW=10.
    GROUP 11. Initialization of variable or porosity fields
FIINIT(:TLSC:)=EPSIN;FIINIT(KE)=TKEIN;FIINIT(TEM1)=0.5*TW
  ** use log-law for initial W profile
DO JJ=1,NY
+PATCH(IN:JJ:,INIVAL,1,NX,JJ,JJ,1,NZ,1,1)
+GR=0.5*YFRAC(JJ)
IF(JJ.NE.1) THEN
+ JJM1=JJ-1
+ GR=YFRAC(JJM1)+0.5*(YFRAC(JJ)-YFRAC(JJM1))
ENDIF
+GYPLUS=YVLAST*(1.-GR)*US/ENULA
+GWPLUS=LOG(GYPLUS)/0.41+5.25
IF(GYPLUS.LE.11.5) THEN
+ GWPLUS=GYPLUS
ENDIF
+INIT(IN:JJ:,W1,ZERO,GWPLUS*US)
ENDDO
    GROUP 13. Boundary conditions and special sources
WALL(WALLN,NORTH,1,1,NY,NY,1,NZ,1,1);COVAL(WALLN,W1,LOGLAW,0.0)
COVAL(WALLN,TEM1,LOGLAW,TW)
  ** activate pressure-drop calculation in single-slab solver
FDFSOL=T;USOURC=T
PATCH(FDFW1DP,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(FDFW1DP,W1,MASIN,GRND1)
  ** temperature source/sink term for fully-developed flow
PATCH(FDFCWT,PHASEM,1,NX,1,NY,1,NZ,1,1)
COVAL(FDFCWT,TEM1,DTDZ,TW)
    GROUP 15. Termination of sweeps
TSTSWP=-1;LITHYD=6
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
DTF=5.0*ZWLAST/WIN; RELAX(W1,FALSDT,DTF)
CASE (TMODEL) OF
WHEN 1
+ LSWEEP=70; RELAX(KE,LINRLX,0.3); RELAX(EP,LINRLX,0.5)
WHEN 2
+ LSWEEP=70; RELAX(KE,LINRLX,0.3); RELAX(EP,LINRLX,0.5)
WHEN 3
+ LSWEEP=40; RELAX(KE,FALSDT,DTF); RELAX(EP,FALSDT,DTF)
WHEN 4
+ LSWEEP=40; RELAX(KE,FALSDT,DTF); RELAX(OMEG,FALSDT,DTF)
ENDCASE
    GROUP 18. Limits on variables or increments to them
VARMIN(W1)=1.E-10
    GROUP 22. Spot-value print-out
IYMON=NY-2;NZPRIN=1;NYPRIN=2;IYPRF=1;NUMCLS=5
    GROUP 24. Dumps for restarts
  ** store(stan) returns incorrect stanton number in RESULT file
WALPRN=T;STORE(YPLS,SKIN,STAN,STRS)
  ** compute expected Nusselt number from Petukhov
XR=1.07+12.7*(PRNDTL(H1)**.666-1.)*(FRIC/8.)**0.5
NUSS=REY*PRNDTL(H1)*FRIC/(8.*XR)
NUSS
STORE(PRPS);EX(PRPS)=33;FIINIT(PRPS)=33;PRNDTL(TEM1)=CONDFILE

  ** mat no. rho enul cp kond expan
  **   1        air
CSG10='q1'
  MATFLG=T;NMAT=1
  33 1. 1.E-6  1.0 3.333E-7 0