TALK=T;RUN( 1, 1)
 ** LOAD(x307) from the x Input Library
    GROUP 1. Run title and other preliminaries
TEXT(REALISABLE KE_2D FLOW PAST A SQUARE RIB: T307
TITLE
  DISPLAY
  The case considered is 2D, steady, incompressible, turbulent flow
  past a surface-mounted square rib in a channel. This case has
  been studied experimentally by D.Crabb et al, Proc. 4th Brazilian
  Congress on Mech. Engng., Florianopolis, Brazil, p415,(1997).
 
  The height H of the rib is 8.5% of that of the channel. The flow
  Reynolds number based on channel bulk velocity and rib height
  H is 300,000. The inlet plane is located 6H upstream of the rib,
  and the outlet plane 20H downstream of the rib. A fixed-pressure
  boundary condition is applied at the outlet, and uniform flow
  profiles are specified at the inlet.
 
  The case is set up to run one of six versions of the k-e model
  with scalable wall functions, and four different versions of the 
  k-w model.
 
  For this case, the main parameter characterising separation is
  the length of the separation zone behind the rib. The experimental
  and computed results are:
 
          KE  MMK   KL   CK   RKE  RNG  KW  KWR  MKW  SST  EXPT
 
   Lr/H = 8.3 11.4 11.2 15.1 12.9 13.7  8.0 11.7 8.3  12.3 12.3  
 
  where the separation length Lr is measured from the front of the
  rib. These results are not grid independent, and the mesh is not
  fine enough to resolve the separation regions around the rib.
  The standard k-e & k-w models seriously underpredict the length
  of the separation behind the rib, whilst the k-w-SST, the Wilcox
  2008 k-w and realisable k-e models give closest agreement with 
  the data.
  ENDDIS
 
    The AUTOPLOT sequence below provides a plot of the axial
    velocity W1 along the bottom surface of the solution domain
    versus normalised axial distance. The axial coordinate 0.0
    corresponds to the rear surface of the rib. The reattachment
    point behind the cube corresponds the axial location where W1
    changes from negative to positive.
   AUTOPLOT USE
    AUTOPLOT
    FILE
    PHIDA 3
 
    D 1 W1 Y 1
    PLOT
    LEVEL Y 0
    SHIFT X -7 1
    REDR
    SCALE X 0 15
    msg Press e to END
  ENDUSE
CHAR(CTURB)
REAL(HRIB,CLUP,CLDOWN,HCHAN,WCHAN,RECHAN)
REAL(RERIB,UIN,TKEIN,EPSIN,MIXL,FRIC,OMIN)
INTEGER(NYC,NZC,NZUP,NZDOWN,NYUP)
     ** Calculation of domain specifications
HRIB=1.0;UIN=1.0
HCHAN=11.75*HRIB
CLUP=6.*HRIB;CLDOWN=20.*HRIB

RERIB=3.E5
RECHAN=RERIB*HCHAN/HRIB

  ** NB: Channel hydraulic diameter=2.*hchan
RECHAN=2.*RECHAN
  ** Estimate inlet conditions assuming fully-developed
	 duct flow at the inlet
FRIC=1./(1.82*LOG10(RECHAN)-1.64)**2
FRIC 
TKEIN=0.25*UIN*UIN*FRIC
MIXL=0.045*HCHAN;EPSIN=0.1643*TKEIN**1.5/MIXL
 
NYC=28;NYUP=62
NZUP=34;NZC=12;NZDOWN=64
    GROUP 3. X-direction grid specification
    GROUP 4. Y-direction grid specification
NREGY=2
IREGY=1;GRDPWR(Y,NYC,HRIB,1.0)
IREGY=2;GRDPWR(Y,-NYUP,-(HCHAN-HRIB),1.08)
    GROUP 5. Z-direction grid specification
NREGZ=3
IREGZ=1;GRDPWR(Z,NZUP,-CLUP,-1.05)
IREGZ=2;GRDPWR(Z,NZC,HRIB,1.0)
IREGZ=3;GRDPWR(Z,NZDOWN,-CLDOWN,1.04)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1);SOLUTN(P1,Y,Y,Y,N,N,N);STORE(ENUT)
SOLUTN(V1,Y,Y,Y,P,P,N);SOLUTN(W1,Y,Y,Y,P,P,N)
MESG( Enter the required turbulence model:
MESG(  KE   -  Standard k-e model
MESG(  MMK  -  MMK  k-e model
MESG(  KL   -  KL   k-e model
MESG(  CK   -  Chen-Kim  k-e model
MESG(  RKE  -  Realisable k-e model (default)
MESG(  RNG  -  RNG k-e model
MESG(  KW   -  Wilcox 1988 k-w model
MESG(  KWR  -  Wilcox 2008 k-w model
MESG(  KWM  -  Menter 1992 k-w model
MESG(  KWS  -  k-w SST model
MESG(
READVDU(CTURB,CHAR,RKE)
CASE :CTURB: OF
WHEN KE,2
+ TEXT(K-E 2D SQUARE RIB FLOW  :T307
+ MESG(Standard k-e model
+ TURMOD(KEMODL)
WHEN MMK,3
+ TEXT(MMK K-E SQUARE RIB FLOW :T307
+ MESG(MMK k-e model
+ TURMOD(KEMMK)
WHEN KL,2
+ TEXT(KL  K-E SQUARE RIB FLOW :T307
+ MESG(KL  k-e model
+ TURMOD(KEKL)
WHEN CK,2
+ TEXT(Chen-Kim K-E SQUARE RIB FLOW :T307
+ MESG(Chen Kim k-e model
+ TURMOD(KECHEN)
WHEN RNG,3
+ TEXT(RNG K-E SQUARE RIB FLOW :T307
+ MESG(RNG k-e model
+ TURMOD(KERNG)
WHEN RKE,3
+ TEXT(RK K-E SQUARE RIB FLOW :T307
+ MESG(RK k-e model
+ TURMOD(KEREAL);STORE(C1E)
WHEN KW,2
+ TEXT(Wilcox 1988 k-w SQUARE RIB FLOW :T307
+ MESG(Wilcox 1988 k-w model
+ TURMOD(KWMODL)
+ OMIN=EPSIN/(0.09*TKEIN)
WHEN KWR,3
+ TEXT(Wilcox 2008 k-w SQUARE RIB FLOW :T307
+ MESG(Wilcox 2008 k-w model
+ TURMOD(KWMODLR)
+ OMIN=EPSIN/(0.09*TKEIN)
+ FIINIT(FBP)=1.0;STORE(CDWS)
WHEN KWM,3
TEXT(Menter k-w SQUARE RIB FLOW :T307
+ MESG(Menter 1992 k-w model
+ TURMOD(KWMENTER);;FIINIT(BF1)=1.0
+ OMIN=EPSIN/(0.09*TKEIN);STORE(CDWS)
WHEN KWS,3
TEXT(SST k-w SQUARE RIB FLOW :T307
+ MESG(Menter 1992 k-w SST model
+ TURMOD(KWSST);FIINIT(BF1)=1.0;FIINIT(BF2)=1.0
+ STORE(CDWS);OMIN=EPSIN/(0.09*TKEIN)
ENDCASE

BOOLEAN(KWMOD)
KWMOD=IENUTA.EQ.10.OR.IENUTA.EQ.15.OR.IENUTA.EQ.17.OR.IENUTA.EQ.19
IF(KWMOD) THEN
+ SOLUTN(KE,Y,Y,Y,P,P,P);SOLUTN(OMEG,Y,Y,Y,P,P,P)
ELSE
+ SOLUTN(KE,Y,Y,Y,P,P,P);SOLUTN(EP,Y,Y,Y,P,P,P)
ENDIF
STORE(YPLS)
SCALWF=T  ! Scalable wall functions
    GROUP 8. Terms (in differential equations) & devices
    GROUP 9. Properties of the medium (or media)
RHO1=1.0;ENUL=UIN*HRIB/RERIB
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=UIN;FIINIT(P1)=1.3E-4
FIINIT(KE)=TKEIN;FIINIT(EP)=EPSIN;FIINIT(V1)=0.001*UIN
IF(KWMOD) THEN
+ FIINIT(OMEG)=OMIN
ENDIF
     ** Initialization of variables in blocked region
STORE(PRPS)	 
PATCH(RIB,INIVAL,#1,#1,#1,#1,#2,#2,1,1)
INIT(RIB,PRPS,0.,198)
EGWF=T

    GROUP 13. Boundary conditions and special sources
INLET(INLET,LOW,#1,#NREGX,#1,#NREGY,#1,#1,1,1)
VALUE(INLET,P1,UIN);VALUE(INLET,W1,UIN)
VALUE(INLET,KE,TKEIN);VALUE(INLET,EP,EPSIN)

IF(KWMOD) THEN
+VALUE(INLET,OMEG,OMIN)
ENDIF
 
PATCH(OUTL,HIGH,#1,#NREGX,#1,#NREGY,#NREGZ,#NREGZ,1,1)
COVAL(OUTL,P1,1.0E3,0.0)
COVAL(OUTL,W1,ONLYMS,0.0);COVAL(OUTL,V1,ONLYMS,0.0)
COVAL(OUTL,KE,ONLYMS,0.0);COVAL(OUTL,EP,ONLYMS,0.0)
 
WALL(WALLN,NORTH,#1,#NREGX,#NREGY,#NREGY,#1,#NREGZ,1,1)
WALL(WALLS,SOUTH,#1,#NREGX,#1,#1,#1,#NREGZ,1,1)
IF(KWMOD) THEN
+ COVAL(WALLN,OMEG,GRND2,GRND2)
+ COVAL(WALLS,OMEG,GRND2,GRND2)
ENDIF
 
    GROUP 15. Termination of sweeps
LSWEEP=1800
    GROUP 16. Termination of iterations
SELREF=T
LITER(P1)=50;LITER(KE)=5;LITER(EP)=5
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=ZWLAST/UIN/NZ/2
RELAX(W1,FALSDT,DTF);RELAX(V1,FALSDT,DTF)
IF(KWMOD) THEN
+ LITER(OMEG)=5
+ RELAX(KE,FALSDT,DTF);RELAX(OMEG,FALSDT,DTF)
ELSE
+ KELIN=3
+ RELAX(KE,LINRLX,0.3);RELAX(EP,LINRLX,0.3)
ENDIF
IYMON=NYC-4;IXMON=1;IZMON=NZUP+NZC+10;NPRMON=100
IYMON=30;IXMON=67
    GROUP 23. Field print-out and plot control
ITABL=3;NPLT=25;IPLTL=LSWEEP;NZPRIN=4;NYPRIN=4
TSTSWP=-1
SPEDAT(SET,GXMONI,PLOTALL,L,T)
SPEDAT(SET,OUTPUT,NOFIELD,L,T)
  ** Revised Wilcox moddel
IF(IENUTA.EQ.15) THEN	 
+ SCALWF=F   ! deactivate scalable wall functions
+ RELAX(FBP,LINRLX,0.05)
+ OUTPUT(FBP,Y,N,Y,N,Y,Y);OUTPUT(XWP,Y,N,Y,N,Y,Y)
+ OUTPUT(CDWS,Y,N,Y,N,Y,Y)
+ LSWEEP=3000
ENDIF
  ** k-w-SST model
IF(IENUTA.EQ.19) THEN
+ LSWEEP=2000
ENDIF

OUTPUT(ENUT,Y,N,Y,N,Y,Y)
 
DISTIL=T
CASE :CTURB: OF
WHEN KE,2
+EX(EPKE)=1.136E-01;EX(PRPS)=9.661E-01
+EX(P1  )=1.279E-01;EX(V1  )=6.438E-02 
+EX(W1  )=8.001E-01;EX(KE  )=1.590E-02 
+EX(EP  )=3.310E-03;EX(YPLS)=4.941E+00 
+EX(ENUT)=1.601E-02 
WHEN MMK,3
+EX(P1  )=1.464E-01;EX(V1  )=6.987E-02 
+EX(W1  )=7.850E-01;EX(KE  )=1.200E-02 
+EX(EP  )=2.174E-03;EX(EPKE)=1.096E-01 
+EX(PRPS)=9.661E-01;EX(YPLS)=4.969E+00 
+EX(DWDY)=2.871E-01;EX(DVDZ)=4.714E-02 
+EX(FOMG)=5.301E-01;EX(ENUT)=9.002E-03 
WHEN KL,2
+EX(P1  )=1.454E-01;EX(V1  )=6.957E-02
+EX(W1  )=7.856E-01;EX(KE  )=1.200E-02
+EX(EP  )=2.182E-03;EX(EPKE)=1.111E-01
+EX(PRPS)=9.661E-01;EX(YPLS)=4.943E+00
+EX(DWDY)=2.858E-01;EX(DVDZ)=4.671E-02
+EX(FOMG)=5.751E-01;EX(ENUT)=1.397E-02
WHEN CK,2
+EX(P1  )=1.590E-01;EX(V1  )=7.342E-02
+EX(W1  )=7.788E-01;EX(KE  )=8.821E-02
+EX(EP  )=1.647E-03;EX(EPKE)=1.196E-01
+EX(PRPS)=9.661E-01;EX(YPLS)=4.981E+00
+EX(ENUT)=1.129E-02
WHEN RKE,3
+EX(P1  )=1.512E-01;EX(V1  )=7.132E-02 
+EX(W1  )=7.815E-01;EX(KE  )=1.017E-02 
+EX(EP  )=1.776E-03;EX(PRPS)=9.661E-01 
+EX(YPLS)=4.931E+00;EX(C1E )=4.664E-01 
+EX(DWDZ)=7.146E-02;EX(DWDY)=2.955E-01 
+EX(DVDZ)=4.759E-02;EX(DVDY)=7.292E-02 
+EX(EPKE)=1.117E-01;EX(CMU )=1.096E-01
+EX(ENUT)=1.484E-02
WHEN RNG,3
+EX(P1  )=1.561E-01;EX(V1  )=7.308E-02
+EX(W1  )=7.816E-01;EX(KE  )=9.889E-02
+EX(EP  )=1.682E-03;EX(EPKE)=1.181E-01
+EX(PRPS)=9.661E-01;EX(YPLS)=4.934E+00
+EX(ENUT)=1.149E-02
WHEN KW,2
+EX(P1  )=1.271E-01;EX(V1  )=6.348E-02
+EX(W1  )=8.010E-01;EX(KE  )=1.747E-02
+EX(EP  )=3.564E-03;EX(EPKE)=9.661E-11
+EX(PRPS)=9.661E-01;EX(YPLS)=5.138E+00
+EX(OMEG)=1.182E+00;EX(ENUT)=1.770E-02
WHEN KWR,3
+EX(P1  )=1.450E-01;EX(V1  )=6.907E-02 
+EX(W1  )=7.812E-01;EX(KE  )=1.266E-02 
+EX(EP  )=2.346E-03;
+EX(PRPS)=9.661E-01;EX(YPLS)=5.091E+00 
+EX(CDWS)=1.075E-01;EX(DWDZ)=6.958E-02 
+EX(DWDY)=2.747E-01;EX(DVDZ)=4.752E-02 
+EX(DVDY)=7.104E-02;EX(GEN1)=4.946E+00 
+EX(FBP )=9.607E-01;EX(XWP )=1.437E+00 
+EX(OMEG)=1.109E+00;EX(ENUT)=1.331E-02 
WHEN KWM,3
+EX(P1  )=1.279E-01;EX(V1  )=6.430E-02
+EX(W1  )=8.090E-01;EX(KE  )=1.633E-02
+EX(EP  )=3.330E-03;
+EX(PRPS)=9.661E-01;EX(YPLS)=4.999E+00
+EX(CDWS)=2.394E-01;EX(LTLS)=6.637E+00
+EX(WDIS)=1.560E+00;EX(BF1 )=7.176E-01
+EX(OMEG)=1.219E+00;EX(ENUT)=1.649E-02
WHEN KWS,3
+EX(P1  )=1.467E-01;EX(V1  )=7.121E-02 
+EX(W1  )=7.809E-01;EX(KE  )=1.151E-02 
+EX(EP  )=1.893E-03;
+EX(PRPS)=9.661E-01;EX(YPLS)=4.845E+00 
+EX(CDWS)=6.563E-02;EX(LTLS)=6.637E+00 
+EX(WDIS)=1.560E+00;EX(GEN1)=5.684E+00 
+EX(BF2 )=9.350E-01;EX(BF1 )=7.172E-01 
+EX(OMEG)=1.183E+00;EX(ENUT)=1.255E-02 
ENDCASE
STOP