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 Press  and 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