GROUP 1. Run title and other preliminaries TEXT(LAM-BRE KE_1D DEVEL CHANNL FLOW :T206 TITLE DISPLAY The case considered is 2d fully-developed turbulent flow in a plane closed channel at a Re=1.E5. The turbulence is simulated by use of the Lam-Bremhorst low-Re k-e model. The calculation integrates down to the wall and the 1d solution is obtained by use of the single-slab solver with a specified mass flow rate. A non-uniform grid is employed so as to concentrate cells very close to the wall. For this purpose a grid is generated which is a geometric progression with the property that the ratio of any two adjacent cell lengths is a constant. The calculation may also be performed with the low-Re Chen-kim k-e model and the Wilcox low-Re k-omega model. The Lam-Bremhorst model produces a friction factor f=8*tauw/(rho*wb**2)=0.0178, whereas the Chen-Kim k-e and k-omega models produce values of f=0.0163, which are in closer agreement with the measured value of 0.0156. No grid-dependency studies have been carried out. All 3 models produce similar velocity and turbulence profiles, but the k-omega model has the advantage of very smooth and rapid convergence behaviour and not requiring the additional calculation of a wall distance. ENDDIS The following AUTOPLOT use file produces three plots; the first is the axial velocity profile; the second is the turbulence energy profile; and the third is the eddy-viscosity profile. AUTOPLOT USE file phi 5 da 1 w1;col9 1 msg Velocity (W1) profile msg Press RETURN to continue pause clear da 1 ke;col9 1 msg KE profile msg Press RETURN to continue pause clear da 1 ENUT;col9 1 msg Press e to END ENDUSE CHAR(CTURB,TLSC);BOOLEAN(VARLAM);VARLAM=T REAL(HGHT,WIN,REY,REYH,DHYDR,MIXL,FRIC,DPDZ,MASIN,DTF) REAL(TKEIN,EPSIN,DELT1,DELT2,DELT3,US,EPSPLS) REAL(KFAC,DELY,GR,GYPLUS,GWPLUS,GWR);INTEGER(JJM,JJM1) ** NB: The the hydraulic diameter is equal to 2.*duct height, so that pipe-flow correlations still apply with diameter replaced by 2.*height HGHT=0.1;WIN=1.0; REY=1.E5;DHYDR=2.*HGHT; REYH=2.*REY FRIC=1./(1.82*LOG10(REYH)-1.64)**2 US=WIN*(FRIC/8.)**0.5;DPDZ=0.5*RHO1*WIN*WIN*FRIC/DHYDR ** estimate initial values from K+=2 & EP+=1./(ak*Y+) EPSPLS=1./(0.41*100.) TKEIN=2.*US*US;ENULA=WIN*HGHT/REY;EPSIN=EPSPLS*US**4/ENULA GROUP 4. Y-direction grid specification ENULA=WIN*HGHT/REY ** define first dely from wall and the grid-expansion factor Kfac which defines a constant ratio of lengths of two adjacent cells. DELT1=1.*ENULA/US;KFAC=1.1;DELY=DELT1/(0.5*HGHT) ** calculate NY from dely & Kfac REAL(AA) AA=(YVLAST/DELY)*(KFAC-1.0)+1.0;AA=LOG(AA)/LOG(KFAC)+1.0001 NY=AA ** define uniform grid initially IREGY=1;GRDPWR(Y,NY,YVLAST,1.0) ** compute expanding grid from north boundary YFRAC(NY)=1.0 DO JJ=NY,2,-1 + JJM=JJ-1 + YFRAC(JJM)=YFRAC(JJ)-DELY + DELY=KFAC*DELY ENDDO YVLAST=0.5*HGHT GROUP 5. Z-direction grid specification ZWLAST=0.1*HGHT GROUP 7. Variables stored, solved & named SOLVE(W1);STORE(ENUT,LEN1);SOLUTN(W1,P,P,P,P,P,N) MESG( Enter the required turbulence model: MESG( CK - Chen-Kim low-Re k-e model MESG( LB - Lam-Bemhorst low-Re k-e model (default) MESG( KO - Wilcox low-re k-omega model MESG( READVDU(CTURB,CHAR,LB) CASE :CTURB: OF WHEN CK,2 + TEXT(CHEN-KIM KE_1D DEVEL CHANNL FLOW :T206 + MESG(Chen-Kim low-Re k-e model + TURMOD(KECHEN-LOWRE);STORE(REYN);TLSC=EP WHEN LB,2 + MESG(Lam-Bremhorst low-Re k-e model + TURMOD(KEMODL-LOWRE);STORE(REYN);TLSC=EP WHEN KO,2 + TEXT(K-OMEGA_1D DEVEL CHANNL FLOW :T206 + MESG(k-omega low-Re model + TURMOD(KOMODL-LOWRE);TLSC=OMEG + STORE(EP);EPSIN=EPSIN/(0.09*TKEIN) 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 9. Properties of the medium (or media) ENUL=ENULA ** test for ground-set enul IF(VARLAM) THEN + TMP1=CONST;TMP1A=0.0;ENUL=LINTEM ENDIF GROUP 11. Initialization of variable or porosity fields FIINIT(:TLSC:)=EPSIN;FIINIT(KE)=TKEIN ** use log-law for initial W profile DO JJ=1,NY +PATCH(IN:JJ:,INIVAL,1,NX,JJ,JJ,1,NZ,1,1) +GR=0.5*YFRAC(JJ) IF(JJ.NE.1) THEN + JJM1=JJ-1 + GR=YFRAC(JJM1)+0.5*(YFRAC(JJ)-YFRAC(JJM1)) ENDIF +GYPLUS=YVLAST*(1.-GR)*US/ENULA +GWPLUS=LOG(GYPLUS)/0.41+5.25 IF(GYPLUS.LE.11.5) THEN + GWPLUS=GYPLUS ENDIF +INIT(IN:JJ:,W1,ZERO,GWPLUS*US) ENDDO GROUP 13. Boundary conditions and special sources WALLCO=GRND3 WALL(WALLN,NORTH,1,1,NY,NY,1,NZ,1,1);COVAL(WALLN,W1,LOGLAW,0.0) MASIN=RHO1*WIN*0.5*HGHT ** Activate single-slab solver with specified mass flow rate FDSOLV(FLOW,MASIN) GROUP 15. Termination of sweeps LSWEEP=20;LITHYD=10 GROUP 16. Termination of iterations GROUP 17. Under-relaxation devices DTF=5.0*ZWLAST/WIN; RELAX(W1,FALSDT,DTF) IF(:TLSC:.EQ.OMEG) THEN + DTF=DTF*10.; RELAX(KE,FALSDT,DTF); RELAX(:TLSC:,FALSDT,DTF) ELSE + KELIN=3; RELAX(KE,LINRLX,0.6); RELAX(EP,LINRLX,1.0) ENDIF GROUP 18. Limits on variables or increments to them VARMIN(W1)=1.E-10 GROUP 22. Spot-value print-out ITABL=3;IYMON=NY-2;NZPRIN=1;NYPRIN=2;IYPRF=1;NUMCLS=5 TSTSWP=-1 GROUP 24. Dumps for restarts STORE(FMU,REYT,FONE,FTWO,STRS,YPLS,SKIN)