TALK=T;RUN( 1, 1)
** LOAD(254) from the PHOENICS Input Library
------------------------------------------------------------------
PHOTON USE
p
up z;vec x 1 sh;pau;clear
con p1 x 1 fi;.01;pau;clear
con tem1 x 1 fi;.01;pau;cl
ENDUSE
------------------------------------------------------------------
GROUP 1. Run title and other preliminaries
TEXT(Free Convection Between Vertical Plates
#cls
TITLE
DISPLAY
The problem considered is laminar free convection of air
between a pair of parallel, vertically-mounted flat plates of
length H. These plates are maintained either at uniform
heat fluxes (UHF) qw or uniform wall temperatures (UWT) tw.
The solution domain has one plane of symmetry at the centre
of the gap B between the two plates. The external pressure
is the same at both the top and bottom of the plate, so the
flow is driven entirely by buoyancy forces.
The Boussinesq assumption is made in the simulations, and the
flow is characterised by the Grashof number GrB and the
gap-to-length ratio B/H.
#pause
The solution domain comprises the gap between the
two channel plates, as well as the free environment upstream
and downstream of the channel.
For economy, the solution domain may be restricted to the gap
alone, although this approach is less accurate for developing
flow because of the empirical representation of total pressure
losses at the channel inlet and outlet.
The practical application of free convection in vertical
channels may be found in the electronics and telecommunications
industry.
#pause
The flow is fully developed if GrB*(B/H) < 1 for a UWT
condition, and when GrB*(B/H) < 100 for a UHF condition (see
Aung, Int.J. Heat Mass Transfer, Vol.15,p1577, 1972).
Developing flow has been considered by Aung et al, Int.J.Heat
Mass Transfer, Vol.15, p2293, 1972 & Morrone et al, Int.J.Heat
Mass Transfer, Vol.40, No.5, p993, 1997).
#pause
The dimensionless equations solved are:
dW/dZ + dV/dY = 0 ;
W*dZ/dZ + V*dV/dY = d/dY(dW/dY) + d/dZ(dW/dZ)
- dP/DZ + GrB*T
W*dT/dZ + V*dT/dY = (1./Pr)*d/dY(dT/dY)
+ (1./Pr)*d/dZ(dT/dZ)
where
Z = z/B ; Y = y/B ; W = w*B/enul ; V = v*B/enul
P = (p-p,ref)*(B**2)/(rho*enul**2)
T = k*(t-tref)/(qw*B) ; Pr = rho*enul*cp/k
GrB = g*beta*qw*B**4/(k*enul**2)
#pause
This Q1 may be used to run 6 cases which are documented in
the literature. The results are summarised below in terms of
the dimensionless volumetric flow rate PSI(=wbulk*B/enul):
Developing Flow (UHF)
---------------------
phoenics phoenics
reduced domain full domain Morrone
GrH=1.E5 B/H=0.5 19.0 27.0 30.0
GrH=1.E4 B/H=0.8 12.5 19.0 21.5
#pause
Fully-Developed Flow
--------------------
phoenics Aung
reduced & full domain
GrB=1.E2 B/H=0.04 (UWT) 8.16 8.33
GrB=1.E2 B/H=0.04 (UHF) 14.80 17.13
where GrH=GrB*(H/B)**4a n The dimensionless volumetric flow
rate produced by PHOENICS can be obtained from the RESULT file
by taking TWICE the printed NETT source of R1 at INLET.
#pause
ENDDIS
BOOLEAN(ISOWAL,FDEV,FDOM)
MESG( Enter y for isothermal wall and n for a uniform heat flux
MESG( Default = uniform heat flux
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
+ ISOWAL=T
+ MESG( Isothermal wall condition
ELSE
+ ISOWAL=F
+ MESG( Uniform heat flux condition
ENDIF
MESG( Do you want to simulate fully-developed flow? (y/n)
MESG( Default = developing flow
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
+ FDEV=T
+ MESG( Fully-developed flow
ELSE
+ FDEV=F
+ MESG( Developing flow
ENDIF
MESG( Do you want to solve in the gap only? (y/n)
MESG( Default = solution domain is gap + free environment
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
+ FDOM=F
+ MESG( Solution domain restricted to gap
ELSE
+ FDOM=T
+ MESG( Solution domain includes free environment
ENDIF
REAL(PLATEB,PLATEL,PBDL,GRASH,GPRL,TAMB,QWAL,KOND,VBUOY,PIN,KLOSS)
REAL(PSI,PSI2,TMAX,PMAX,GRASL)
GPRL=0.71;TAMB=0.
IF(FDEV) THEN
+ MESG(Enter Grashof number - default 1.E2
+ READVDU(GRASH,REAL,1.E2)
+ MESG(Enter Plate Width/Length ratio - default 0.04
+ READVDU(PBDL,REAL,0.04)
ELSE
+ MESG(Enter Grashof number - default 1.E4
+ READVDU(GRASL,REAL,1.E4)
+ MESG(Enter Plate Width/Length ratio - default 0.8
+ READVDU(PBDL,REAL,0.8)
+ GRASH = GRASL*(PBDL)**4
ENDIF
PLATEB=1.0 ; PLATEL=PLATEB/PBDL
GROUP 4. Y-direction grid specification
IF(FDOM) THEN
+ NREGY=2
+ IREGY=1;GRDPWR(Y,25,PLATEB, 0.8)
+ IREGY=2;GRDPWR(Y,30,0.5*PLATEB,1.2)
ELSE
+ GRDPWR(Y,30,0.5*PLATEB,1.0)
ENDIF
GROUP 5. Z-direction grid specification
IF(FDOM) THEN
+ NREGZ=3
+ IREGZ=1;GRDPWR(Z, 20,PLATEL,0.8)
+ IREGZ=2;GRDPWR(Z,-30,PLATEL,1.2)
+ IREGZ=3;GRDPWR(Z, 15,PLATEL,1.8)
ELSE
+ GRDPWR(Z,30,PLATEL,1.0)
ENDIF
GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1,TEM1);SOLUTN(TEM1,Y,Y,Y,P,P,N)
SOLUTN(V1,P,P,P,P,P,N);SOLUTN(W1,P,P,P,P,P,N)
SOLUTN(P1,Y,Y,Y,P,P,P)
GROUP 8. Terms (in differential equations) & devices
TERMS(TEM1,N,Y,Y,N,Y,N)
IF(FDOM) THEN
+ SCHEME(LUS,V1,W1)
ENDIF
GROUP 9. Properties of the medium (or media)
CP1=1.0 ; RHO1=1.0 ; ENUL=1.0
KOND=RHO1*CP1*ENUL/GPRL
PRNDTL(TEM1)=-KOND
GROUP 11. Initialization of variable or porosity fields
VBUOY=(GRASH*PLATEL)**0.5
FIINIT(TEM1)=TAMB
IF(FDOM) THEN
+ FIINIT(W1)=0.0
ELSE
+ FIINIT(W1)=0.25*VBUOY
+ PIN=-0.5*VBUOY*VBUOY
**Linear pressure variation set as initial field.
+ PATCH(PINIT,LINVLZ,1,1,1,NY,1,NZ,1,1)
+ COVAL(PINIT,P1,-PIN/ZWLAST,PIN)
ENDIF
CONPOR(PLATE,0.0,CELL,1,1,-#1,-#1,-#2,-#2)
GROUP 13. Boundary conditions and special sources
** Inlet boundary condition
PATCH(INLET,LOW,1,1,1,NY,1,1,1,1)
IF(FDEV) THEN
+ COVAL(INLET,P1,1.E2,0.)
ELSE
** Inflow proportional to square root of pressure difference.
+ KLOSS=1.0 ; COVAL(INLET,P1,-2./(1.+KLOSS),0.)
ENDIF
IF(FDOM) THEN
+ COVAL(INLET,P1,1.E2,0.)
ENDIF
** momentum convected into cell assumed equal to that leaving
COVAL(INLET,TEM1,ONLYMS,TAMB)
IF(.NOT.FDOM) THEN
** Wall friction boundary condition.
+ WALL (HOTWALL,SOUTH,1,1,1,1,1,NZ,1,1)
+ COVAL(HOTWALL,W1,1.0,0.)
ENDIF
** Wall thermal boundary condition.
QWAL=1.0
IF(ISOWAL) THEN
** isothermal wall
+ PATCH(HEATFLUX,SWALL,1,1,$2,$2,#2,#2,1,1)
+ COVAL(HEATFLUX,TEM1,1.0,1.0)
ELSE
** uniform heat flux wall
+ PATCH(HEATFLUX,SWALL,1,1,$2,$2,#2,#2,1,1)
+ COVAL(HEATFLUX,TEM1,FIXFLU,QWAL)
ENDIF
** Exit boundary condition.
PATCH(EXIT,HIGH,1,1,1,NY,NZ,NZ,1,1)
** Total pressure in external environment = 0
IF(FDEV) THEN
+ COVAL(EXIT,P1,1.E2,0.)
ELSE
** Inflow proportional to square root of pressure difference.
+ KLOSS=1.0 ; COVAL(EXIT,P1,-2./(1.+KLOSS),0.)
ENDIF
IF(FDOM) THEN
+ COVAL(EXIT,P1,1.E2,0.)
ENDIF
COVAL(EXIT,TEM1,ONLYMS,TAMB)
IF(FDOM) THEN
+ PATCH(FREEBH,SOUTH,1,1,1,1,#3,#3,1,1)
+ COVAL(FREEBH,P1,1.E2,0.)
+ COVAL(FREEBH,TEM1,ONLYMS,TAMB)
ENDIF
**Boussinesq approximation for buoyancy force.
BUOYE=TAMB; DVO1DT=GRASH; BUOYC=-1.0
PATCH(BUOY,PHASEM,1,1,1,NY,1,NZ,1,1)
COVAL(BUOY,W1,FIXFLU,BOUSS)
GROUP 15. Termination of sweeps
LSWEEP=800 ; LITER(TEM1)=50 ; LITER(P1)=50
ENDIT(P1)=GRND1
GROUP 17. Under-relaxation devices
REAL(DTF);DTF=50.*PLATEL/(0.25*VBUOY)
RELAX(V1,FALSDT,DTF); RELAX(W1,FALSDT,DTF)
RELAX(TEM1,FALSDT,1.E9)
GROUP 22. Spot-value print-out
IYMON=NY-1; IZMON=NZ-1
GROUP 23. Field print-out and plot control
ITABL=3; NPLT=10; NZPRIN=2
TSTSWP=-1
** fully-developed analytical solutions of Aung (1972)
psi = dimensionless volume flow rate
pmax = dimensionless maximum pressure
tmax = dimensionless peak temperature
IF(FDEV) THEN
IF(ISOWAL) THEN
** uniform temperature condition
+ PSI=GRASH/12.
ELSE
** uniform heat flux condition
+ PSI = (GRASH/(12.*GPRL*PBDL))**0.5
+ PMAX = -17.525*PSI**3/(GPRL*GRASH)
+ TMAX = PSI*24./GRASH
+ PMAX
+ TMAX
ENDIF
ELSE
** numerical correlation of Morrone et al (1997) for
developing flow valid for 1E2 < GrL < 1E5 and
0.2 < B/H < 3.0
+ PSI=0.81*(GRASL)**0.384;PSI=PSI*(PBDL)**1.17
ENDIF
MESG(Dimensionless flow rate for 1/2 channel
PSI=0.5*PSI
PSI
DISTIL=T
EX(P1 )=1.254E+02;EX(V1 )=1.595E+00;EX(W1 )=9.028E+00
STORE(PRPS);EX(PRPS)=7.902E-01;EX(VPOR)=7.902E-01;EX(TEM1)=4.486E-02
LIBREF = 254
STOP