TALK=T;RUN( 1, 1)
 ** LOAD(x100) from the x Input Library
    GROUP 1. Run title and other preliminaries
TEXT(CHEN-KIM K-E_1D PLANE COUETTE FLOW :T100
TITLE
  DISPLAY
  The problem considered is turbulent couette flow in a plane
  channel with one moving wall. The shear stress is uniform
  across the flow and hence the axial pressure-gradient is zero.
  The velocity profile is S-shaped and symmetrical about the 
  central plane, and so the average velocity is one half of
  the velocity of the moving wall. In the calculation made here
  the Reynolds number is 1.E5 based on the channel height and 
  the average velocity.
  ENDDIS
 
  Calculations with the standard k-e model plus wall functions
  show that the turbulent kinetic energy k is essentially
  uniform across the flow. This result is in accordance with the
  exact solution of the k-e model equations reported by Henry
  and Reynolds [1984], but not with the measurements reported
  in the literature, which show that k is about 70% larger near
  the walls than it is at the centre plane. Gibson [1988] and
  Schneider [1988] have demonstrated that the inclusion of a
  counter-gradient diffusion term in the k-equation leads to very
  good agreement with the measured k profiles.
 
  Calculations may be made with the high-Re forms of either the 
  standard k-e model, the Chen-Kim k-e model, the RNG k-e model,
  the Wilcox 1988 and 1992 k-w model, the Menter k-w model, and 
  the k-w SST model. For this case all models produce similar 
  results.
 
  ENDDIS
 
  Comparisons between measured and computed skin-friction
  coefficients Cf are made below:
                                          Cf
 
    Data - Telbany & Reynolds [1982]    3.070E-3
    Standard k-e model                  3.378E-3
    Chen-Kim k-e model                  3.300E-3
    RNG k-e model                       3.243E-3
    Wilcox 1988 k-w model               3.280E-3
    Wilcox 2008 k-w model               3.220E-3	
    Menter low-Re k-w model             3.286E-3
    k-w SST low-Re model                3.167E-3
 
   where Cf =2.*tauw/(rho*uav)**2 and Uav is the average velocity.
 
   The following AUTOPLOT use file produces two plots;
   the first is the axial velocity profile; and the
   second is the turbulence energy profile.
 
   AUTOPLOT USE
   file; phida 3
 
   cl; da 1 w1; col9 1
   msg        Velocity (W1) profile; Press RETURN to continue
   pause; cl; da 1 ke; col9 1; scale y 0. 2.e-3; redraw
   msg        KE profile. Press RETURN, then e to END
   pause
   ENDUSE
 
CHAR(CTURB,TLSC)
REAL(HEIGHT,WTOP,REY,TKEIN,EPSIN,MIXL,DTF,WAV,WSTAR,MASIN)
HEIGHT=0.1;WTOP=1.0; REY=1.E5;WAV=0.5*WTOP
  ** wstar from data of El Telbany & Reynolds [1982]
WSTAR=WAV*0.196/LOG10(REY);TKEIN=WSTAR*WSTAR/.3
MIXL=0.045*HEIGHT;EPSIN=TKEIN**1.5/MIXL*0.1643
    GROUP 4. Y-direction grid specification
ENUL=WAV*HEIGHT/REY;NY=30;YVLAST=HEIGHT;GRDPWR(Y,NY,YVLAST,1.0)
DTF=20.*ZWLAST/WAV
    GROUP 7. Variables stored, solved & named
SOLVE(W1);STORE(ENUT,LEN1);SOLUTN(W1,P,P,P,P,P,N)
STORE(STRS)
(stored of CF is 2.*STRS/(RHO1*:WAV:*:WAV:))
MESG( Enter the required turbulence model:
MESG(  CK   -  Chen-Kim k-e model (default)
MESG(  KE   -  Standard k-e model
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,CHEN)
CASE :CTURB: OF
WHEN CK,2
+ MESG(Chen-Kim k-e model
+ TURMOD(KECHEN);KELIN=1;TLSC=EP
WHEN KE,2
+ TEXT(K-E_1D PLANE COUETTE FLOW :T100
+ MESG(Standard k-e model
+ TURMOD(KEMODL);KELIN=1;TLSC=EP
WHEN RNG,3
+ TEXT(RNG K-E_1D PLANE COUETTE FLOW :T100
+ MESG(RNG k-e model
+ TURMOD(KERNG);KELIN=1;TLSC=EP
+ STORE(ETA,ALF,GEN1)
+ OUTPUT(ALF,Y,N,P,Y,Y,Y);OUTPUT(ETA,Y,N,P,Y,Y,Y)
+ DTF=2.*ZWLAST/WAV
WHEN KW,2
+ TEXT(Wilcox k-w_1D PLANE COUETTE FLOW :T100
+ MESG(Wilcox 1988 k-w model
+ TURMOD(KWMODL);TLSC=OMEG
+ EPSIN=EPSIN/(0.09*TKEIN)
WHEN KWR,3
+ TEXT(WilcoxR k-w_1D PLANE COUETTE FLOW :T100
+ MESG(Wilcox 2008 k-w model
+ TURMOD(KWMODLR);TLSC=OMEG
WHEN KWM,3
TEXT(Menter k-w_1D PLANE COUETTE FLOW :T100
+ MESG(Menter 1992 k-w model
+ TURMOD(KWMENTER);TLSC=OMEG
+ STORE(EP);EPSIN=EPSIN/(0.09*TKEIN)
+ STORE(BF1);FIINIT(BF1)=1.0
+ STORE(DKDY,DFDY,CDWS,SIGK,SIGW)
WHEN KWS,3
TEXT(SST k-w_1D PLANE COUETTE FLOW :T100
+ MESG(Menter k-w SST model
+ TURMOD(KWSST);TLSC=OMEG
+ STORE(EP);EPSIN=EPSIN/(0.09*TKEIN)
+ STORE(BF1,BF2);;FIINIT(BF1)=1.0;;FIINIT(BF2)=1.0
+ STORE(DKDY,DFDY,CDWS,BF1,BF2,SIGK,SIGW)
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 11. Initialization of variable or porosity fields
FIINIT(:TLSC:)=EPSIN;FIINIT(KE)=TKEIN
PATCH(ICOUF,LINVLY,1,1,1,NY,1,NZ,1,1)
INIT(ICOUF,W1,WTOP/HEIGHT,0.0)
    GROUP 13. Boundary conditions and special sources
  ** moving upper wall
WALL(WALLN,NORTH,1,1,NY,NY,1,NZ,1,1);COVAL(WALLN,W1,LOGLAW,WTOP)
  ** stationary bottom wall
WALL(WALLS,SOUTH,1,1,1,1,1,NZ,1,1)
    GROUP 15. Termination of sweeps
LSWEEP=10;LITHYD=20
    GROUP 16. Termination of iterations
MASIN=RHO1*WAV*HEIGHT; RESREF(W1)=1.E-12*MASIN*WAV
RESREF(KE)=RESREF(W1)*TKEIN; RESREF(:TLSC:)=RESREF(W1)*EPSIN
    GROUP 17. Under-relaxation devices
VARMIN(W1)=1.E-10;WALPRN=T
    GROUP 22. Spot-value print-out
IYMON=NY-2;NPLT=1;NZPRIN=1;NYPRIN=1;IYPRF=1;TSTSWP=-1
SPEDAT(SET,GXMONI,PLOTALL,L,T)
    GROUP 24. Dumps for restarts
RELAX(W1,FALSDT,DTF)
  ** heavier relaxation for k-w SST model
IF(IENUTA.EQ.19) THEN
+RELAX(KE,FALSDT,DTF/20.);RELAX(OMEG,FALSDT,DTF/20.)
ELSE
+RELAX(KE,FALSDT,DTF);RELAX(:TLSC:,FALSDT,DTF)
ENDIF
  ** heavier relaxation for Wilcox 1988 k-w model
IF(IENUTA.EQ.10) THEN
+RELAX(KE,FALSDT,DTF/20.);RELAX(OMEG,FALSDT,DTF/20.)
ENDIF 
DISTIL=T
CASE :CTURB: OF
WHEN KE,2
+EX(W1  )=5.000E-01;EX(KE  )=1.426E-03
+EX(EP  )=2.179E-03;EX(LEN1)=8.333E-03
+EX(ENUT)=1.720E-04;EX(CF  )=2.252E-04
+EX(STRS)=2.815E-05
WHEN CK,2
+EX(W1  )=5.000E-01;EX(KE  )=1.393E-03
+EX(EP  )=2.148E-03;EX(CF  )=2.200E-04
+EX(STRS)=2.750E-05;EX(LEN1)=8.014E-03
+EX(ENUT)=1.635E-04
WHEN RNG,3
+EX(W1  )=5.000E-01;EX(KE  )=1.412E-03
+EX(EP  )=2.129E-03;EX(GEN1)=2.024E+02
+EX(ALF )=6.019E-01;EX(ETA )=3.614E+00
+EX(CF  )=2.162E-04;EX(STRS)=2.703E-05
+EX(LEN1)=7.704E-03;EX(ENUT)=1.558E-04
WHEN KW,2
+EX(W1  )=5.000E-01;EX(KE  )=1.383E-03
+EX(EP  )=2.149E-03;EX(OMEG)=1.719E+01
+EX(LEN1)=7.889E-03;EX(ENUT)=1.603E-04
+EX(CF  )=2.186E-04;EX(STRS)=2.733E-05
WHEN KWR,3
+EX(W1  )=5.000E-01;EX(KE  )=1.375E-03 
+EX(EP  )=2.153E-03;EX(DWDY)=4.340E+00 
+EX(FBP )=1.000E+00;EX(OMEG)=1.716E+01 
+EX(CF  )=2.151E-04;EX(STRS)=2.689E-05 
+EX(LEN1)=7.783E-03;EX(ENUT)=1.566E-04 
+EX(GEN1)=4.630E+02
WHEN KWM,3
+EX(W1  )=5.000E-01;EX(KE  )=1.386E-03
+EX(EP  )=2.152E-03;EX(CDWS)=2.355E-11
+EX(DFDY)=2.469E+03;EX(DKDY)=3.722E-03
+EX(OMEG)=1.717E+01;EX(CF  )=2.191E-04
+EX(STRS)=2.739E-05;EX(LEN1)=7.919E-03
+EX(ENUT)=1.610E-04;EX(SIGW)=2.000E+00
+EX(SIGK)=2.000E+00;EX(LTLS)=8.352E-04
+EX(WDIS)=2.503E-02;EX(BF1 )=1.000E+00
WHEN KWS,3
+EX(SIGW)=2.000E+00;EX(SIGK)=2.000E+00 
+EX(CDWS)=3.011E-11;EX(DFDY)=2.451E+03 
+EX(BF1 )=1.000E+00;EX(OMEG)=1.704E+01 
+EX(LEN1)=7.963E-03;EX(ENUT)=1.605E-04 
+EX(W1  )=5.000E-01;EX(BF2 )=1.000E+00
+EX(LTLS)=8.352E-04;EX(WDIS)=2.503E-02
+EX(KE  )=1.359E-03;EX(EP  )=2.115E-03 
+EX(DKDY)=1.173E-02;EX(GEN1)=5.892E+02 
+EX(CF  )=2.111E-04 ;EX(STRS)=2.639E-05 
ENDCASE
 LIBREF = 100
STOP