GROUP 1. Run title and other preliminaries
TEXT(REALISABLE KE_3D FLOW PAST A CUBE IN A CHANNEL: T308
TITLE
  DISPLAY
  The case considered is 3D, steady, incompressible, turbulent flow
  past a surface-mounted cube in a channel. The flow separates in
  front of the cube to form a primary and secondary vortex, and
  the main vortex wraps as a horse-shoe vortex around the cube into
  the wake. The flow separates at the front corners of the cube on
  the roof and the side walls; the reattaches on the side walls but
  not on the roof. A large separation region develops behind the
  cube which interacts with the horseshoe vortex. In the experiments
  vortex shedding as observed from the side walls, and due to
  momentum exchange with the wake, this will lead to a shorter
  separation length than is reported here for a steady simulation.
 
  The height of the cube is 50% of that of the channel. The flow
  Reynolds number based on channel bulk velocity and cube height
  H is 40,000. The inlet plane is located 7H upstream of the cube,
  and the outlet plane 10H downstream of the cube. Because of
  symmetry conditions, only half of the width of the flow is
  calculated. 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 any one of 6 variants of the k-e model
  with scalable wall functions, namely the standard model and the 
  MMK, Kato-Launder, RNG, Chen-Kim and realisable variants. An
  option is provided to also run the standard k-w model. The
  case has been studied experimentally by Martinuzzi & Tropea
  [J.Fluids Engng, 115, p85-92,1993] and numerically by Lakehal &
  Rodi [J.Wind Eng. Ind.Aerodyn, 67 & 68, p65-78, 1997].
 
  For this case, the main parameters that characterise separation
  are the frontal stagnation point Ys/H, the primary upstream
  separation point Zf/H, the roof reattachment point Zr/H and the
  length of the separation zone behind the cube Zb/H. The 
  experimental and computed results for Zb are given below:
 
               K-E   KL   MMK   RKE  CHEN  RNG    KO    EXPT
       Zb/H = 2.08  2.72  2.81  2.46  3.1  2.92  1.70   1.61
 
  These results are not grid independent, and the mesh is not fine
  enough to resolve the expected separation on the roof nor to
  capture adequately the upstream and downstream separation regions. 
  For this rather coarse mesh, all the k-e models overpredict Zb, 
  and the standard k-w model gives close agreement with the data. 
  The standard k-e model is known to produce too small a separation 
  on the roof with unrealistic roof reattachment. The modified k-e 
  models produce longer separation regions and no reattachment, 
  which is in agreement with the data. However, the present 
  computations employ insufficient mesh resolution to exemplify 
  these benefits. However, it is likely that more mesh and the 
  inclusion of unsteady effects are required for a much improved 
  prediction of the separation length behind the cube.
  ENDDIS
 
    This AUTOPLOT sequence provides a plot of the axial
    velocity W1 at the symmetry plane and along the bottom surface 
	of the solution domain versus normalised axial distance X. The 
	axial coordinate 0.0 corresponds to the rear surface of the 
	cube. The reattachment point behind the cube corresponds the 
	axial location where W1 changes from negative to positive.
 
  AUTOPLOT USE
  file                                                                            
  phida 3                                                                         
                                                                                
  d 1 w1 y 1 x 1                                                                  
  plot                                                                            
  redr                                                                            
  shift x -8                                                                    
  1                                                                               
  scale     
  level y 0 
  scale x 0 5
  ENDUSE
 
CHAR(CTURB)
REAL(HCUBE,CLUP,CLDOWN,CHIGHT,CWIDTH)
REAL(REYNO,UIN,TKEIN,EPSIN,MIXL,FRIC,OMIN)
INTEGER(NYC,NZC,NZUP,NZDOWN,NYUP,NXC,NXUP,JKO)
JKO=0
     ** Calculation of domain specifications
HCUBE=1.0;UIN=1.0
CHIGHT=2.*HCUBE;CLUP=7.0*HCUBE;CLDOWN=10.*HCUBE
CWIDTH=4.5*HCUBE
REYNO=4.E4
 
FRIC=0.018;TKEIN=0.25*UIN*UIN*FRIC
MIXL=0.09*CHIGHT;EPSIN=0.1643*TKEIN**1.5/MIXL

NXC=12;NXUP=26
NYC=18;NYUP=18
NZUP=38;NZC=12;NZDOWN=34
    GROUP 3. X-direction grid specification
NREGX=2
IREGX=1;GRDPWR(X,NXC,0.5*HCUBE,1.0)
IREGX=2;GRDPWR(X,NXUP,-(CWIDTH-0.5*HCUBE),1.08)
    GROUP 4. Y-direction grid specification
NREGY=2
IREGY=1;GRDPWR(Y,NYC,-HCUBE,-1.06)
IREGY=2;GRDPWR(Y,NYUP,-(CHIGHT-HCUBE),1.06)
    GROUP 5. Z-direction grid specification
NREGZ=3
IREGZ=1;GRDPWR(Z,NZUP,-CLUP,-1.05)
IREGZ=2;GRDPWR(Z,NZC,HCUBE,1.0)
IREGZ=3;GRDPWR(Z,NZDOWN,-CLDOWN,1.07)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1,W1);SOLUTN(P1,Y,Y,Y,N,N,N);STORE(ENUT)
SOLUTN(U1,P,P,P,P,P,N);SOLUTN(V1,P,P,P,P,P,N)
SOLUTN(W1,P,P,P,P,P,N)
MESG( Enter the required turbulence model:
MESG(  KEM  -  Standard k-e model
MESG(  CHEN -  Chen-Kim k-e model
MESG(  RNG  -  RNG  k-e model
MESG(  MMK  -  MMK  k-e model
MESG(  KLM  -  KL   k-e model
MESG(  KO   -  k-omega model
MESG(  RKE  -  Realisable k-e model (default)
MESG(
READVDU(CTURB,CHAR,RKE)
CASE :CTURB: OF
WHEN KEM,3
+ TEXT(K-E SURFACE-MOUNTED CUBE FLOW :T308
+ MESG(Standard k-e model
+ TURMOD(KEMODL)
WHEN CHEN,4
+ TEXT(CHEN KE SURFACE-MOUNTED CUBE FLOW :T308
+ MESG(Chen k-e model
+ TURMOD(KECHEN)
WHEN RNG,3
+ TEXT(RNG KE SURFACE-MOUNTED CUBE FLOW :T308
+ MESG(RNG k-e model
+ TURMOD(KERNG)
WHEN MMK,3
+ TEXT(MMK K-E SURFACE-MOUNTED CUBE FLOW :T308
+ MESG(MMK k-e model
+ TURMOD(KEMMK)
WHEN KLM,3
+ TEXT(KL K-E SURFACE-MOUNTED CUBE FLOW  :T308
+ MESG(KL k-e model
+ TURMOD(KEKL)
WHEN RKE,3
+ TEXT(RK K-E SURFACE-MOUNTED CUBE FLOW  :T308
+ MESG(RK k-e model
+ TURMOD(KEREAL);STORE(C1E)
WHEN KO,2
+ TEXT(K-O SURFACE-MOUNTED CUBE FLOW :T308
+ MESG(Standard k-w model
+ TURMOD(KOMODL);STORE(EP)
+ JKO=1;OMIN=EPSIN/(0.09*TKEIN)
ENDCASE
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*HCUBE/REYNO
    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(JKO.EQ.1) THEN
+ FIINIT(OMEG)=OMIN
ENDIF
     ** Initialization of variables in blocked region
CONPOR(CUBE,0.0,CELL,-#1,-#1,-#1,-#1,-#2,-#2)
    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)
 
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(JKO.EQ.1) THEN
+ VALUE(INLET,OMEG,OMIN)
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP=1000
    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)
RELAX(U1,FALSDT,DTF)
RELAX(KE,FALSDT,DTF); RELAX(EP,FALSDT,DTF)
IF(JKO.EQ.1) THEN
+ RELAX(OMEG,FALSDT,DTF)
ENDIF 
IYMON=NY-4;IXMON=1;IZMON=NZ-4;NPRMON=100
    GROUP 23. Field print-out and plot control
ITABL=3;NPLT=10;IPLTL=LSWEEP;NZPRIN=2;NYPRIN=2
TSTSWP=-1

DISTIL=T
STORE(PRPS); EX(PRPS)=9.774E-01; EX(VPOR)=9.774E-01
CASE :CTURB: OF
WHEN CHEN,4
+ EX(P1  )=6.623E-02;EX(U1  )=3.166E-02
+ EX(V1  )=2.526E-02;EX(W1  )=9.512E-01
+ EX(KE  )=5.345E-03;EX(EP  )=1.880E-03
+ EX(YPLS)=4.930E+00;EX(ENUT)=4.175E-03
WHEN MMK,3
+ EX(P1  )=6.158E-02;EX(U1  )=2.959E-02
+ EX(V1  )=2.309E-02;EX(W1  )=9.529E-01
+ EX(KE  )=5.234E-03;EX(EP  )=1.703E-03
+ EX(YPLS)=4.916E+00;EX(DWDY)=3.226E-01
+ EX(DWDX)=1.439E-01;EX(DVDZ)=3.854E-02
+ EX(DVDX)=3.064E-02;EX(DUDZ)=3.872E-02
+ EX(DUDY)=3.584E-02;EX(FOMG)=4.360E-01
+ EX(ENUT)=1.407E-03
WHEN KLM,3
+ EX(P1  )=6.087E-02;EX(U1  )=2.910E-02
+ EX(V1  )=2.268E-02;EX(W1  )=9.531E-01
+ EX(KE  )=5.307E-03;EX(EP  )=1.759E-03
+ EX(YPLS)=4.917E+00;EX(DWDY)=3.190E-01
+ EX(DWDX)=1.432E-01;EX(DVDZ)=3.824E-02
+ EX(DVDX)=2.999E-02;EX(DUDZ)=3.867E-02
+ EX(DUDY)=3.519E-02;EX(FOMG)=4.458E-01
+ EX(ENUT)=4.466E-03
WHEN KEM,3
+EX(P1  )=5.786E-02;EX(U1  )=2.727E-02 
+EX(V1  )=2.217E-02;EX(W1  )=9.572E-01 
+EX(KE  )=8.058E-03;EX(EP  )=3.191E-03 
+EX(YPLS)=4.973E+00;EX(ENUT)=5.034E-03 
WHEN KO,2
+EX(P1  )=5.701E-02;EX(U1  )=2.619E-02  
+EX(V1  )=2.080E-02;EX(W1  )=9.569E-01  
+EX(KE  )=1.372E-02;EX(EP  )=4.705E-03  
+EX(YPLS)=5.111E+00;EX(OMEG)=1.820E+00  
+EX(ENUT)=7.999E-03 
WHEN RNG,3
+EX(P1  )=6.545E-02;EX(U1  )=3.136E-02 
+EX(V1  )=2.533E-02;EX(W1  )=9.524E-01 
+EX(KE  )=5.400E-03;EX(EP  )=1.857E-03 
+EX(YPLS)=4.888E+00;EX(ENUT)=3.650E-03 
WHEN RKE,3
+EX(P1  )=6.155E-02;EX(U1  )=2.941E-02
+EX(V1  )=2.435E-02;EX(W1  )=9.567E-01
+EX(KE  )=6.595E-03;EX(EP  )=2.144E-03
+EX(YPLS)=4.925E+00;EX(C1E )=4.491E-01
+EX(DWDZ)=7.234E-02;EX(DWDY)=3.142E-01
+EX(DWDX)=1.289E-01;EX(DVDZ)=4.032E-02
+EX(DVDY)=5.053E-02;EX(DVDX)=3.186E-02
+EX(DUDZ)=3.929E-02;EX(DUDY)=3.484E-02
+EX(DUDX)=5.655E-02;EX(EPKE)=1.823E-01
+EX(CMU )=1.496E-01;EX(ENUT)=7.406E-03
ENDCASE