```
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

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
```