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