PHOTON USE p view z gr ou z 1 msg secondary-flow velocity vectors set vec comp;u1 v1 0. set vec ref;1.E-2 vec z 1 sh pause msg Pressand then to END pause ENDUSE GROUP 1. Run title and other preliminaries TEXT(RSTM_3D DVLPD DUCT FLOW/HEAT TRANS :T608 TITLE DISPLAY The case considered is three-dimensional fully-developed, turbulent flow and heat transfer of an incompressible, constant-property fluid in a rectangular duct of aspect ratio alfa=1.0, where alfa=height/width. The calculations are performed for a Reynolds number of 5.E4, a Prandtl number of 3.0 and a constant temperature boundary condition is applied on all four sides of the duct. The solution exploits symmetry by performing the calculation over one quadrant of the duct cross section only. The flow Reynolds number is based on the hydraulic diameter of the duct. The RSTM calculations are made with one of two variants, namely: the IPM pressure-strain model (IRSMHM=0) or the QIM pressure-strain model (IRSMHM=2). A transport model (IRSMSM=2) is employed for closure of the turbulent heat-flux components. ENDDIS Fully-developed turbulent flows in straight ducts of non-circular cross-section are accompanied by secondary motions in the plane perpendicular to the streamwise direction. These motions, which are not present in fully-developed laminar flow, are caused by the turbulence field, i.e by gradients of the Reynolds stresses. These turbulence-driven secondary velocities are typically the order of 2-3% of the streamwise bulk velocity, and they transport high-momentum fluid towards the duct corners resulting in a bulging of the velocity contours towards the corners. Turbulence models employing an isotropic eddy viscosity, such as the k-e model, fail to reproduce the secondary motions. However, the RSTM is successful because it is capable of simulating the anisotropy of the in-plane Reynolds normal stresses. For data on the friction factor and Nusselt number the literature recommends that for most engineering calculations, it is sufficiently accurate to use the circular-tube correlations with the hydraulic diameter replacing the circular-tube diameter in the Reynolds number and Nusselt number ( see for example, 'Handbook of Heat Transfer Fundamentals', Ed. W.M.Rohsenow, J.P.Hartnett & E.N.Ganic', Chapter 7, McGraw Hill, 2nd Edition, ( 1985). Experimental data and numerical results on the secondary motions and turbulence field may be found in Demuren and Rodi ( Journal of Fluid Mechanics, P189, Vol.140, [1984] ). For the values of Re and alpha considered here, the data indicates that the friction factor f = 0.021 and the Nusselt number Nu = 221.6. The Petukhov correlation ( see below ) is used to estimate the Nusselt number. Both RSTM predictions yield essentially the same values of f and Nu, namely: f = 0.019 and Nu = 218, and the ratio of the centre-line velocity to bulk velocity is 1.19 which agrees well with the data. Experiments indicate that along the corner bisector, the ratio of the peak secondary velocity to the friction velocity is 0.2. The QIM model predicts a ratio of 0.16 whereas the IPM model yields 0.07, demonstrating that the simpler IPM produces is not adequate for an accurate prediction of the secondary motions.It should be mentioned that no grid-sensitivity studies have been performed. IRSMHM=2;BOOLEAN(HEAT);HEAT=T; REAL(AWAL,TW,NUSS,COND) REAL(HEIGHT,WIDTH,ALF,HD2,WD2,WIN,DPDZ,REY,FRIC,QIN,DTDZ) REAL(US,DELT,AIN,DHYD,TKEIN,EPSIN,MIXL,FLOWIN,DTF,CP) ** ALFA = HEIGHT/WIDTH WIDTH=2.0;ALF=1.0;HEIGHT=ALF*WIDTH HD2=0.5*HEIGHT;WD2=0.5*WIDTH;WIN=1.0 REY=5.E4;DHYD=4.*WIDTH*HEIGHT/(2.*HEIGHT+2.*WIDTH) ** compute expected pressure-drop for SATELLITE printout FRIC=1./(1.82*LOG10(REY)-1.64)**2 DPDZ=FRIC*RHO1*WIN*WIN/(2.*DHYD);US=WIN*(FRIC/8.)**0.5 FRIC DPDZ DHYD GROUP 3. X-direction grid specification AIN=HD2*WD2;ENUL=WIN*DHYD/REY;DELT=2.*40.*ENUL/US NREGX=2; REGEXT(X,WD2);IREGX=1;GRDPWR(X,24,WD2-DELT,0.8) IREGX=2;GRDPWR(X,1,DELT,1.0) GROUP 4. Y-direction grid specification NREGY=2; REGEXT(Y,HD2);IREGY=1;GRDPWR(Y,24,HD2-DELT,0.8) IREGY=2;GRDPWR(Y,1,DELT,1.0) GROUP 7. Variables stored, solved & named SOLVE(P1,U1,V1,W1);STORE(ENUT,LEN1) GROUP 8. Terms (in differential equations) & devices TERMS(W1,N,P,P,P,P,P) IF(HEAT) THEN + SOLVE(H1);TERMS(H1,N,P,P,P,P,P) + PRNDTL(H1)=3.0;IRSMSM=2 ENDIF PATCH(WALLT,NWALL,1,NX,NY,NY,1,NZ,1,1) COVAL(WALLT,W1,1.0,0.0);COVAL(WALLT,U1,1.0,0.0) PATCH(WALLS,EWALL,NX,NX,1,NY,1,NZ,1,1) COVAL(WALLS,W1,1.0,0.0);COVAL(WALLS,V1,1.0,0.0) DTF=0.04;TURMOD(REYSTRS,DTF,WALLT,WALLS) GROUP 9. Properties of the medium (or media) RHO1=1.0;FLOWIN=RHO1*WIN*AIN;TKEIN=0.25*WIN*WIN*FRIC MIXL=0.09*DHYD;EPSIN=TKEIN**1.5/MIXL*0.1643 IF(HEAT) THEN ** prescribe energy flow from slab and fluid specific heat estimated from Dittus-Boelter Nu=0.023*Re**0.8*Pr**0.4 with (Tw-Tb)=5.0 + AWAL=(WD2+HD2)*ZWLAST;NUSS=0.023*REY**0.8*PRNDTL(H1)**0.4 + CP=1.0;COND=RHO1*CP*ENUL/PRNDTL(H1) + QIN=NUSS*5.0*COND/DHYD NUSS ** compute d(tbulk)/dz for input to single-slab thermal solver + DTDZ=QIN*AWAL/(CP*FLOWIN);TW=10. ENDIF DTDZ GROUP 11. Initialization of variable or porosity fields FIINIT(W1)=WIN;FIINIT(KE)=TKEIN;FIINIT(EP)=EPSIN FIINIT(W2RS)=2.*FIINIT(KE)/3.;FIINIT(V2RS)=FIINIT(W2RS) FIINIT(U2RS)=FIINIT(W2RS);FIINIT(VWRS)=0.3*FIINIT(KE) FIINIT(UWRS)=0.3*FIINIT(KE) IF(HEAT) THEN + FIINIT(H1)=0.5*TW ENDIF GROUP 12. Convection and diffusion adjustments PATCH(GP12CONH,CELL,1,NX,1,NY,1,NZ,1,1) COVAL(GP12CONH,U1,0.0,0.0);COVAL(GP12CONH,V1,0.0,0.0) COVAL(GP12CONH,W1,0.0,0.0);COVAL(GP12CONH,EP,0.0,0.0) COVAL(GP12CONH,U2RS,0.0,0.0);COVAL(GP12CONH,V2RS,0.0,0.0) COVAL(GP12CONH,W2RS,0.0,0.0);COVAL(GP12CONH,UWRS,0.0,0.0) COVAL(GP12CONH,VWRS,0.0,0.0);COVAL(GP12CONH,UVRS,0.0,0.0) GROUP 13. Boundary conditions and special sources PATCH(RELIEF,CELL,NX/2,NX/2,NY/2,NY/2,1,NZ,1,1) COVAL(RELIEF,P1,FIXP,0.0);COVAL(RELIEF,H1,ONLYMS,SAME) FDFSOL=T;USOURC=T;PATCH(FDFW1DP,VOLUME,1,NX,1,NY,1,NZ,1,1) COVAL(FDFW1DP,W1,FLOWIN,GRND1) PATCH(SMPLS,SOUTH,1,NX,1,1,1,NZ,1,1);COVAL(SMPLS,VWRS,GRND1,0.0) COVAL(SMPLS,VTRS,GRND1,0.0) PATCH(SMPLW,WEST,1,1,1,NY,1,NZ,1,1);COVAL(SMPLW,UWRS,GRND1,0.0) COVAL(SMPLW,UTRS,GRND1,0.0) IF(HEAT) THEN ** constant wall-temperature boundary condition + PATCH(FDFCWT,PHASEM,1,NX,1,NY,1,NZ,1,1) + COVAL(FDFCWT,H1,DTDZ,TW) + COVAL(WALLS,H1,LOGLAW,TW);COVAL(WALLT,H1,LOGLAW,TW) + COVAL(GP12CONH,H1,0.0,0.0) + COVAL(GP12CONH,VTRS,0.0,0.0);COVAL(GP12CONH,UTRS,0.0,0.0) ENDIF GROUP 15. Termination of sweeps LSWEEP=120;LITHYD=6;LITER(W1)=2;LITER(H1)=5 GROUP 16. Termination of iterations RESREF(P1)=1.E-12*WIN*AIN RESREF(W1)=1.E-12*DPDZ*ZWLAST*AIN RESREF(U1)=RESREF(W1); RESREF(V1)=RESREF(W1) RESREF(KE)=RESREF(P1)*RHO1*WIN*TKEIN RESREF(EP)=RESREF(P1)*RHO1*WIN*EPSIN IF(HEAT) THEN + RELAX(H1,FALSDT,5.0); RELAX(UTRS,FALSDT,0.05) + RELAX(VTRS,FALSDT,0.05); RESREF(H1)=1.E-12*QIN*ZWLAST*AWAL + QIN + COND=RHO1*CP*ENUL/PRNDTL(H1) + COND ENDIF GROUP 17. Under-relaxation devices DTF=4.*ZWLAST/WIN; RELAX(U1,FALSDT,DTF) RELAX(V1,FALSDT,DTF); RELAX(W1,FALSDT,DTF) GROUP 22. Spot-value print-out IXMON=NX-2;IYMON=NY-2;TSTSWP=-1 GROUP 23. Field print-out and plot control NPLT=5;NYPRIN=3;NXPRIN=3;WALPRN=T ** compute expected Nusselt number from Petukhov correlation and printout from satellite REAL(XR) XR=1.07+12.7*(PRNDTL(H1)**.666-1.)*(FRIC/8.)**0.5 NUSS=REY*PRNDTL(H1)*FRIC/(8.*XR) NUSS