GROUP 1. Run title and other preliminaries
 
TEXT(2D Turb Duct Flow And Heat Trans     
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=0.5, where alfa=height/width. The calculations are performed
  for a Reynolds number of 1.E5, a Prandtl number of 3.0 and a
  constant temperature boundary condition is applied on all four
  sides ofthe 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.
  ENDDIS
 
  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).
 
  For the Reynolds number and duct aspect ratio considered here,
  the data indicates that the friction factor f = 0.018 and the
  Nusselt number Nu = 392.0. The Petukhov correlation ( see below )
  is used to estimate the Nusselt number. The PHOENICS prediction
  yields f = 0.017 and Nu = 373.0. These values agree fairly well
  with the data. It should be mentioned that no grid-sensitivity
  studies have been performed.
 
  It should be noted that the k-e turbulence model, which is used
  in the present calculation, is unable to predict the turbulence-
  driven secondary motions that are present in fully-developed
  turbulent flow in non-circular ducts. The present calculation
  predicts zero secondary flows, and one would have to resort to a
  suitable, but more complex turbulence model to capture the
  secondary motions, e.g. a Reynolds stress transport model.
 
  Finally, one may BFC=T with NONORT=T or F to verify that very
  similar results are obtained with the body-fitted-coordinates
  facility.
 
BOOLEAN(HEAT);HEAT=T
REAL(HEIGHT,WIDTH,ALF,HD2,WD2,WIN,DPDZ,REY,FRIC)
REAL(US,DELT,AIN,DHYD,TKEIN,EPSIN,MIXL,FLOWIN)
REAL(QIN,DTDZ,CP,COND,AWAL,TW,NUSS)
   ** ALFA = HEIGHT/WIDTH
WIDTH=2.0;ALF=0.5;HEIGHT=ALF*WIDTH
HD2=0.5*HEIGHT;WD2=0.5*WIDTH
WIN=1.0
REY=1.E5;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,14,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,14,HD2-DELT,0.8)
IREGY=2;GRDPWR(Y,1,DELT,1.0)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1,W1)
TURMOD(KEMODL);KELIN=1;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
ENDIF
    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
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=WIN;FIINIT(KE)=TKEIN;FIINIT(EP)=EPSIN
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,KE,0.0,0.0)
COVAL(GP12CONH,EP,0.0,0.0)
    GROUP 13. Boundary conditions and special sources
PATCH(WALLT,NWALL,1,NX,NY,NY,1,NZ,1,1)
COVAL(WALLT,W1,LOGLAW,0.0);COVAL(WALLT,U1,LOGLAW,0.0)
COVAL(WALLT,KE,LOGLAW,LOGLAW);COVAL(WALLT,EP,LOGLAW,LOGLAW)
 
PATCH(WALLS,EWALL,NX,NX,1,NY,1,NZ,1,1)
COVAL(WALLS,W1,LOGLAW,0.0);COVAL(WALLS,V1,LOGLAW,0.0)
COVAL(WALLS,KE,LOGLAW,LOGLAW);COVAL(WALLS,EP,LOGLAW,LOGLAW)
 
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)
 
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)
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP=50;LITHYD=1;LITER(W1)=15;LITER(H1)=10
    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
+ RESREF(H1)=1.E-12*QIN*ZWLAST*AWAL
+ QIN
+ COND=RHO1*CP*ENUL/PRNDTL(H1)
+ COND
ENDIF
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=5.*ZWLAST/WIN
RELAX(U1,FALSDT,DTF);RELAX(V1,FALSDT,DTF)
RELAX(W1,FALSDT,DTF);RELAX(KE,FALSDT,DTF)
RELAX(EP,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
    GROUP 24. Dumps for restarts
YPLS=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