PHOTON USE
   AUTOPLOT
   file
   phi 5
 
   cl
   msg LAMINAR PIPE FLOW WITH HEAT TRANSFER
   msg constant wall heat flux
   msg Reynolds number = 10 Prandtl number = 1.0
   msg Velocity (W1) profile
   msg Blue line --- PHOENICS solution
   msg crosses ---   analytical solution
   da 1 w1;da 1 w1a
   col3 1;blb4 2
   msg press  to continue
   pause
   cl
 
   msg constant wall heat flux
   msg Reynolds number = 10 Prandtl number = 1.0
   msg Temperature (H1) profile
   msg Blue line --- PHOENICS solution
   msg crosses ---   analytical solution
   da 1 h1;da 1 h1a
   col3 1;blb4 2
   msg press  to continue
   pause
 
   msg press  to end
   pause
   end
   END_USE
  DISPLAY
 
TEXT(1D Laminar Pipe Flow And Heat Trans  
TITLE
 
  DISPLAY
  The case considered is laminar flow and heat transfer
  of an incompressible, constant-property fluid in the
  hydrodynamically and thermally developed region of a
  of a circular tube. The heat transfer problem can be
  solved either with a constant wall heat flux ( QFCON=T )
  or a constant wall temperature ( QFCON=F ). The Prandtl
  number is set to 1.0, but it should be noted that the Nusselt
  number is a function only of the thermal boundary condition,
  as follows: Nu = 3.656 for a constant surface temperature
  and Nu = 4.364 for a constant heat flux QIN. For the latter
  case, an analytical solution exists for the radial
  temperature profile:
 
    T = Tw-GAM*(0.1875+0.0625*(R/RO)**4-0.25*(R/RO)**2)
 
  where RO is the tube radius and GAM=4.*QIN/(K*R0) wherein
  K is the thermal conductivity of the fluid.
  ENDDIS
 
  For this case the velocity profile is parabolic:
 
      W = 2.*WIN*(1.-(R/RO)**2)
 
  where WIN is the bulk velocity and RO is the tube
  radius.
 
  The H1 equation is used as a vehicle for solving the temperature,
  so that the H1 equation is divided through by the fluid specific
  heat. For hydrodynamically and thermally developed flow, the
  following source/sink term appears in the temperature equation:
  -rho1*w1*dT/dz. For a constant wall heat flux, dT/dz=dTb/dz=
  dTw/dx where Tb is the bulk temperature and Tw is the wall
  temperature. For a constant wall temperature, dT/dz is not
  constant, rather it is equal to [(Tw-T)/(Tw-Tb)]*dTb/dz. For both
  cases the user must precribe dTb/dz, which may be evaluated from
  an overall energy balance for the slab. For the case of a uniform
  wall temperature, this means that the user must presently specify
  the energy flow out of the slab.
 
    GROUP 1. Run title and other preliminaries
  ** QFCON=T sets constant heat-flux        boundary condition
          =F  "      "     wall-temperature    "       "
 
BOOLEAN(QFCON);QFCON=T
REAL(REY,RIN,DIN,WIN,AIN,DPDZ,FLOWIN)
REAL(QIN,DTDZ,CP,COND,AWAL,TC,TB,TW,BETA,GAM)
REY=10.;RIN=0.1;DIN=2.*RIN;WIN=1.0
    GROUP 2. Transience; time-step specification
CARTES=F
    GROUP 3. X-direction grid specification
AIN=0.5*RIN*RIN*XULAST
    GROUP 4. Y-direction grid specification
NY=30;GRDPWR(Y,NY,RIN,1.0)
    GROUP 5. Z-direction grid specification
  ** AWAL is wall surface area per unit length
AWAL=RIN*XULAST
    GROUP 7. Variables stored, solved & named
   ** W1A used to store analytical solution for W1
SOLVE(W1,H1);STORE(W1A)
    GROUP 8. Terms (in differential equations) & devices
   ** deactivate convection and built-in H1 sources
TERMS(W1,N,N,P,P,P,P);TERMS(H1,N,N,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
PRNDTL(H1)=1.0;RHO1=1.0;ENUT=0.;ENUL=WIN*DIN/REY
FLOWIN=RHO1*WIN*AIN
  ** compute pressure drop for printout from satellite Q1
DPDZ=8.*RHO1*ENUL*WIN/RIN/RIN
DPDZ
  ** define input parameters for thermally-developed flow
IF(QFCON) THEN
  ** specify wall flux, fluid specific heat & datum
     temperature at pipe axis
+ QIN=1.0;CP=1.0;TC=0.0
  ** compute d(tbulk)/dz for input to single-slab
     thermal solver
+ DTDZ=QIN*AWAL/(CP*FLOWIN)
   ** H1A is a store for analytical temperature solution
      with constant heat-flux boundary condition
+ STORE(H1A)
+ COND=RHO1*CP*ENUL/PRNDTL(H1)
+ BETA=QIN*2.*RIN/COND
+ TB=TC+7.*BETA/48.;TW=TB+11.*BETA/48.
+ GAM=4.*QIN*RIN/COND
ELSE
+ QIN=1.0;CP=1.0
+ DTDZ=QIN*AWAL/(CP*FLOWIN)
+ TW=10.
ENDIF
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=WIN;FIINIT(H1)=0.4*DTDZ*ZWLAST
  ** compute analytical solutions
REAL(WA,TA,GR,GR2);INTEGER(JJM1)
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
+GR=GR*YVLAST;GR2=(GR/RIN)**2
+WA=2.*WIN*(1.-GR2)
+INIT(IN:JJ:,W1A,ZERO,WA)
IF(QFCON) THEN
+TA=TW-GAM*(0.1875+0.0625*(GR2**2)-0.25*GR2)
+INIT(IN:JJ:,H1A,ZERO,TA)
ENDIF
ENDDO
    GROUP 13. Boundary conditions and special sources
 
PATCH(WALL,NWALL,1,NX,NY,NY,1,NZ,1,1);COVAL(WALL,W1,1.0,0.0)
 
  ** activate pressure-drop calculation in single-slab solver
FDFSOL=T;USOURC=T
PATCH(FDFW1DP,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(FDFW1DP,W1,FLOWIN,GRND1)
 
  ** temperature equation sources & boundary conditions
IF(QFCON) THEN
  ** constant heat-flux boundary condition
+ PATCH(FDFCHF,PHASEM,1,NX,1,NY,1,NZ,1,1)
+ COVAL(FDFCHF,H1,DTDZ,GRND1);COVAL(WALL,H1,FIXFLU,QIN)
  ** set centre-line temperature as datum
+ PATCH(TFIX,CELL,1,NX,1,1,1,NZ,1,1)
+ COVAL(TFIX,H1,FIXP,TC)
ELSE
  ** constant wall-temperature boundary condition
+ PATCH(FDFCWT,PHASEM,1,NX,1,NY,1,NZ,1,1)
+ COVAL(FDFCWT,H1,DTDZ,TW);COVAL(WALL,H1,1.0/PRNDTL(H1),TW)
ENDIF
 
    GROUP 15. Termination of sweeps
LSWEEP=8;LITHYD=2;LITC=5;ISWC1=4
    GROUP 16. Termination of iterations
RESREF(W1)=1.E-12*DPDZ*ZWLAST*AIN
RESREF(H1)=1.E-12*QIN*ZWLAST*AWAL
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=50.*(YVLAST/NY)**2/ENUL
RELAX(W1,FALSDT,DTF)
    GROUP 19. Data communicated by satellite to GROUND
COND=RHO1*CP*ENUL/PRNDTL(H1)
COND
DIN
    GROUP 21. Print-out of variables
    GROUP 22. Spot-value print-out
IYMON=NY;TSTSWP=-1
    GROUP 23. Field print-out and plot control
NPLT=1;NYPRIN=1;NZPRIN=1
    GROUP 24. Dumps for restarts