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