TALK=T;RUN( 1, 1)
** PHOENICS VALIDATION CASE
GROUP 1. Run title and other preliminaries
TEXT(S-A Low Re_1D DEVELOPED PIPE FLOW
TITLE
DISPLAY
The case considered is fully-developed turbulent flow and heat
transfer in a circular pipe at Re=1.E5 and Pr=3.0. The tube wall
is held at a constant temperature, and the calculation integrates
down to the wall through the viscous sublayer. A non-uniform grid
is employed so as to concentrate cells very close to the wall.
For this purpose a grid is generated which is a geometric
progression with the property that the ratio of any two adjacent
cell lengths is a constant. The turbulence is simulated using
the low Reynolds version of the Spalart-Allmaras model.
This model predicts the following results compared to experimental
data:
f vdp Nu
SA 0.018 4.07 361
Data 0.018 3.75 392
where f =8.*(wstar/wbulk)**2 is the friction factor,
vdp=(wcl-wbulk)/u* is the velocity defect parameter and
Nu=h*d/k is the Nusselt number.
In-Form is used to print the f, vdp and Nu to the PHIDA, RESULT
and inforout files.
ENDDIS
REAL(DIAM,WIN,REY,MIXL,FRIC,DPDZ,MASIN,DTF,TKEIN,EPSIN)
REAL(DELT1,DELT2,DELT3,US,EPSPLS,KFAC,DELY,AA,GR,GYPLUS)
REAL(GWPLUS,GWR,QIN,DTDZ,COND,CP,TW,AIN,AWAL,NUSS,XR,RIN)
DIAM=0.1;WIN=1.0; REY=1.E5;FRIC=1./(1.82*LOG10(REY)-1.64)**2
FRIC
DPDZ=FRIC*RHO1*WIN*WIN/(2.*DIAM);US=WIN*(FRIC/8.)**0.5
** estimate initial values from K+=2 & EP+=1./(ak*Y+)
EPSPLS=1./(0.41*100.); RIN=0.5*DIAM
TKEIN=2.*US*US;ENULA=WIN*DIAM/REY;EPSIN=EPSPLS*US**4/ENULA
** the grid-expansion factor Kfac defines a constant ratio of
lengths of two adjacent cells.
KFAC=1.1
GROUP 3. X-direction grid specification
CARTES=F;XULAST=0.1;AIN=0.5*RIN*RIN*XULAST
GROUP 4. Y-direction grid specification
ENULA=WIN*DIAM/REY
** define first dely from wall
DELT1=1.*ENULA/US;DELY=DELT1/(0.5*DIAM)
** calculate NY from dely & Kfac
AA=(YVLAST/DELY)*(KFAC-1.0)+1.0;AA=LOG(AA)/LOG(KFAC)+1.0001
NY=AA
** define uniform grid initially
IREGY=1;GRDPWR(Y,NY,YVLAST,1.0)
** compute expanding grid from north boundary
YFRAC(NY)=1.0;INTEGER(JJM,JJM1)
DO JJ=NY,2,-1
+ JJM=JJ-1
+ YFRAC(JJM)=YFRAC(JJ)-DELY
+ DELY=KFAC*DELY
ENDDO
YVLAST=0.5*DIAM
GROUP 5. Z-direction grid specification
ZWLAST=0.1*DIAM
GROUP 7. Variables stored, solved & named
SOLVE(W1,ENTI,TEM1);STORE(V1,ENUT,LEN1)
SOLUTN(W1,P,P,P,P,P,N)
STORE(STRS,YPLS,SKIN)
TURMOD(SPALART-ALLMARAS-LOWRE)
KELIN=3
GROUP 8. Terms (in differential equations) & devices
TERMS(W1,N,N,P,P,P,P)
TERMS(TEM1,N,N,P,P,P,P)
GROUP 9. Properties of the medium (or media)
ENUL=ENULA;PRNDTL(H1)=3.0;MASIN=RHO1*WIN*AIN
** 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
NUSS=0.023*REY**0.8*PRNDTL(H1)**0.4;CP=1.0;AWAL=RIN*XULAST
COND=RHO1*CP*ENUL/PRNDTL(H1);QIN=NUSS*5.0*COND/DIAM
NUSS
** compute d(tbulk)/dz for input to single-slab
thermal solver and prescribe wall temperature
DTDZ=QIN*AWAL/MASIN;TW=10.
STORE(KOND,AREH)
GROUP 11. Initialization of variable or porosity fields
FIINIT(TEM1)=0.5*TW
FIINIT(ENTI)=1.4025E-4
** use log-law for initial W profile
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
+GYPLUS=YVLAST*(1.-GR)*US/ENULA
+GWPLUS=LOG(GYPLUS)/0.41+5.25
IF(GYPLUS.LE.11.5) THEN
+ GWPLUS=GYPLUS
ENDIF
+INIT(IN:JJ:,W1,ZERO,GWPLUS*US)
ENDDO
GROUP 13. Boundary conditions and special sources
WALL(WALLN,NORTH,1,1,NY,NY,1,NZ,1,1);COVAL(WALLN,W1,LOGLAW,0.0)
COVAL(WALLN,TEM1,LOGLAW,TW)
** 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,MASIN,GRND1)
** temperature source/sink term for fully-developed flow
PATCH(FDFCWT,PHASEM,1,NX,1,NY,1,NZ,1,1)
COVAL(FDFCWT,TEM1,DTDZ,TW)
GROUP 15. Termination of sweeps
TSTSWP=-1;LITHYD=6
GROUP 16. Termination of iterations
GROUP 17. Under-relaxation devices
DTF=5.0*ZWLAST/WIN; RELAX(W1,FALSDT,DTF)
LSWEEP=70
RELAX(ENTI,LINRLX,0.3)
SPEDAT(SET,GXMONI,CLASSIC,L,T)
GROUP 18. Limits on variables or increments to them
VARMIN(W1)=1.E-10
GROUP 22. Spot-value print-out
IYMON=NY-2;NZPRIN=1;NYPRIN=2;IYPRF=1;NUMCLS=5
GROUP 24. Dumps for restarts
** store(stan) returns incorrect stanton number in RESULT file
STORE(YPLS,SKIN,STAN,STRS,HTCO)
** compute expected Nusselt number from Petukhov
XR=1.07+12.7*(PRNDTL(H1)**.666-1.)*(FRIC/8.)**0.5
NUSS=REY*PRNDTL(H1)*FRIC/(8.*XR)
NUSS
STORE(PRPS);EX(PRPS)=33;FIINIT(PRPS)=33;PRNDTL(TEM1)=CONDFILE
** mat no. rho enul cp kond expan
** 1 air
CSG10='q1'
MATFLG=T;NMAT=1
33 1. 1.E-6 1.0 3.333E-7 0
(stored of FRIC is 8.*STRS/(:WIN:*:WIN:))
(stored of UTAU is STRS^0.5)
(STORED of USTR is UTAU[,NY,])
(stored of VDP is (W1[&1&]-WIN)/UTAU[&NY&])
(STORED of UPL is W1/USTR)
(STORED of YPL is ((YVLAST-YG)*USTR)/ENULA)
** compute Tbulk & Nusselt
(stored AREH is AHIGH)
(stored TSUM is 0.0)
(stored TSUM at FDFCWT is SUM(W1*AREH*TEM1))
(stored ASUM is 0.0)
(stored ASUM at FDFCWT is SUM(W1*AREH))
(stored TB at FDFCWT is TSUM/ASUM)
(stored of NUSS is HTCO*:DIAM:*(:TW:-TEM1[&NY&1])/(:COND:*(:TW:-TB)))
(make ffac is 0.0)
(make ustar is 0.0)
(make vdus is 0.0)
(make nus is 0.0)
(store1 of ffac at walln is fric)
(store1 of ustar at walln is utau)
(store1 of vdus at walln is vdp)
(store1 of nus at walln is nuss)
** print to inforout file
(print of f is ffac)
(print of u* is ustar)
(print of (ucl-ub)/u* is vdus)
(print of Nu is nus)
DISTIL=T
EX(V1 )= 1.000E-10
EX(W1 )= 7.849E-01
EX(NUSS)= 6.222E+00
EX(TB )= 5.055E+00
EX(ASUM)= 1.250E-04
EX(TSUM)= 6.318E-04
EX(YPL )= 4.395E+02
EX(UPL )= 1.668E+01
EX(VDP )= 4.074E+00
EX(USTR)= 4.706E-02
EX(UTAU)= 8.212E-04
EX(FRIC)= 3.055E-04
EX(HTCO)= 5.447E-04
EX(STAN)= 2.331E-02
EX(AREH)= 2.155E-06
EX(KOND)= 3.333E-07
EX(EPKE)= 1.000E-10
EX(DWDY)= 3.540E+02
EX(LTLS)= 1.714E-04
EX(WDIS)= 7.667E-03
EX(VOR1)= 3.540E+02
EX(SKIN)= 6.994E-02
EX(YPLS)= 8.560E-03
EX(STRS)= 3.819E-05
EX(ENUT)= 7.414E-05
EX(TEM1)= 5.893E+00
EX(ENTI)= 7.494E-05
LIBREF = 0
STOP