GROUP 1. Run title and other preliminaries
TEXT(LAM-BRE KE_1D DEVEL CHANNL FLOW   :T206
TITLE
  DISPLAY
  The case considered is 2d fully-developed turbulent flow in a
  plane closed channel at a Re=1.E5. The turbulence is simulated by
  use of the Lam-Bremhorst low-Re k-e model. The calculation
  integrates down to the wall and the 1d solution is obtained by
  use of the single-slab solver with a specified mass flow rate. 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 calculation may also be
  performed with the low-Re Chen-kim k-e model and the Wilcox low-Re
  k-omega model. The Lam-Bremhorst model produces a friction factor
  f=8*tauw/(rho*wb**2)=0.0178, whereas the Chen-Kim k-e and k-omega
  models produce values of f=0.0163, which are in closer agreement
  with the measured value of 0.0156. No grid-dependency studies have
  been carried out. All 3 models produce similar velocity and
  turbulence profiles, but the k-omega model has the advantage of
  very smooth and rapid convergence behaviour and not requiring the
  additional calculation of a wall distance.
  ENDDIS
 
   The following AUTOPLOT use file produces three plots;
   the first is the axial velocity profile; the second
   is the turbulence energy profile; and the third is the
   eddy-viscosity profile.
 
   AUTOPLOT USE
   file
   phi 5
 
   da 1 w1;col9 1
   msg Velocity (W1) 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 ENUT;col9 1
   msg Press e to END
   ENDUSE
 
CHAR(CTURB,TLSC);BOOLEAN(VARLAM);VARLAM=T
REAL(HGHT,WIN,REY,REYH,DHYDR,MIXL,FRIC,DPDZ,MASIN,DTF)
REAL(TKEIN,EPSIN,DELT1,DELT2,DELT3,US,EPSPLS)
REAL(KFAC,DELY,GR,GYPLUS,GWPLUS,GWR);INTEGER(JJM,JJM1)
  ** NB: The the hydraulic diameter is equal to 2.*duct height,
         so that pipe-flow correlations still apply
         with diameter replaced by 2.*height
HGHT=0.1;WIN=1.0; REY=1.E5;DHYDR=2.*HGHT; REYH=2.*REY
FRIC=1./(1.82*LOG10(REYH)-1.64)**2
US=WIN*(FRIC/8.)**0.5;DPDZ=0.5*RHO1*WIN*WIN*FRIC/DHYDR
  ** estimate initial values from K+=2 & EP+=1./(ak*Y+)
EPSPLS=1./(0.41*100.)
TKEIN=2.*US*US;ENULA=WIN*HGHT/REY;EPSIN=EPSPLS*US**4/ENULA
    GROUP 4. Y-direction grid specification
ENULA=WIN*HGHT/REY
  ** define first dely from wall and the grid-expansion
     factor Kfac which defines a constant ratio of lengths of
     two adjacent cells.
DELT1=1.*ENULA/US;KFAC=1.1;DELY=DELT1/(0.5*HGHT)
  ** calculate NY from dely & Kfac
REAL(AA)
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
DO JJ=NY,2,-1
+ JJM=JJ-1
+ YFRAC(JJM)=YFRAC(JJ)-DELY
+ DELY=KFAC*DELY
ENDDO
YVLAST=0.5*HGHT
    GROUP 5. Z-direction grid specification
ZWLAST=0.1*HGHT
    GROUP 7. Variables stored, solved & named
SOLVE(W1);STORE(ENUT,LEN1);SOLUTN(W1,P,P,P,P,P,N)
MESG( Enter the required turbulence model:
MESG(  CK   -  Chen-Kim low-Re k-e model
MESG(  LB   -  Lam-Bemhorst low-Re k-e model (default)
MESG(  KO   -  Wilcox low-re k-omega model
MESG(
READVDU(CTURB,CHAR,LB)
CASE :CTURB: OF
WHEN CK,2
+ TEXT(CHEN-KIM KE_1D DEVEL CHANNL FLOW   :T206
+ MESG(Chen-Kim low-Re k-e model
+ TURMOD(KECHEN-LOWRE);STORE(REYN);TLSC=EP
WHEN LB,2
+ MESG(Lam-Bremhorst low-Re k-e model
+ TURMOD(KEMODL-LOWRE);STORE(REYN);TLSC=EP
WHEN KO,2
+ TEXT(K-OMEGA_1D DEVEL CHANNL FLOW    :T206
+ MESG(k-omega low-Re model
+ TURMOD(KOMODL-LOWRE);TLSC=OMEG
+ STORE(EP);EPSIN=EPSIN/(0.09*TKEIN)
ENDCASE
    GROUP 8. Terms (in differential equations) & devices
  ** Deactivate convection
TERMS(W1,N,N,P,P,P,P);TERMS(KE,Y,N,P,P,P,P)
TERMS(:TLSC:,Y,N,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
ENUL=ENULA
   ** test for ground-set enul
IF(VARLAM) THEN
+ TMP1=CONST;TMP1A=0.0;ENUL=LINTEM
ENDIF
    GROUP 11. Initialization of variable or porosity fields
FIINIT(:TLSC:)=EPSIN;FIINIT(KE)=TKEIN
  ** 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
WALLCO=GRND3
WALL(WALLN,NORTH,1,1,NY,NY,1,NZ,1,1);COVAL(WALLN,W1,LOGLAW,0.0)
MASIN=RHO1*WIN*0.5*HGHT
  ** Activate single-slab solver with specified mass flow rate
FDSOLV(FLOW,MASIN)
    GROUP 15. Termination of sweeps
LSWEEP=20;LITHYD=10
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
DTF=5.0*ZWLAST/WIN; RELAX(W1,FALSDT,DTF)
IF(:TLSC:.EQ.OMEG) THEN
+ DTF=DTF*10.; RELAX(KE,FALSDT,DTF); RELAX(:TLSC:,FALSDT,DTF)
ELSE
+ KELIN=3; RELAX(KE,LINRLX,0.6); RELAX(EP,LINRLX,1.0)
ENDIF
    GROUP 18. Limits on variables or increments to them
VARMIN(W1)=1.E-10
    GROUP 22. Spot-value print-out
ITABL=3;IYMON=NY-2;NZPRIN=1;NYPRIN=2;IYPRF=1;NUMCLS=5
TSTSWP=-1
    GROUP 24. Dumps for restarts
STORE(FMU,REYT,FONE,FTWO,STRS,YPLS,SKIN)