PHOTON USE
p;;;;
up x; view z; Gr ou z 1
msg temperature contours
con tem1 z 1 fi;0.001; Pause;con off;red
msg velocity vectors
vec z 1 sh; pause
ENDUSE
GROUP 1. RUN TITLE AND OTHER PRELIMINARIES
TEXT(2D Rayleigh-Benard Convection
TITLE
DISPLAY
The problem considered is the 2d thermal convection of air
between isothermal horizontal walls in an enclosure of 5:1 aspect
ratio with adiabatic side walls. The lower wall is at a higher
temperature than the upper wall. This configuration frequently
occurs in technology, for example, in energy storage systems and
in many industrial heating and cooling systems. For Rayleigh
numbers Ra less than the critical value of 1708, no motion
occurs. For Ra in the range 1708 < Ra < 4E5, the fluid motion
consists of regularly spaced roll cells, while for larger Ra,
the cells break down and the motion is turbulent.
PHOENICS simulations may be performed for laminar flow (the
default) at a Rayleigh number Ra of 6.468E4 or turbulent flow at
Ra=2.35E6. The Boussinesq approximation is used and the energy
equation is solved with the temperature TEM1 as dependent
variable. For turbulent flow, the low-Re k-e model of Lam and
Bremhorst is used with buoyancy extensions and a turbulence
buoyancy coefficient of C3e=1.44.
ENDDIS
The main interest is the heat transfer rate which for large
aspect ratios is given by the correlation Nu=c*Gr**n where,
according to Jakob (Heat Transfer, Vol.1, Chap.25, John Wiley,
1949): c=0.195 and n=0.25 for laminar flow (1E44E5). The Nusselt
number and Grashof number are based on the distance between the
two plates. The comparison between the correlations and PHOENICS
is given below:
Correlations PHOENICS
turbulent flow Nu = 9.951 Nu = 6.345 with C3e=1.44
Nu = 11.852 with C3e=1.35
laminar flow Nu = 3.319 Nu = 3.515
The laminar flow field displays a 4-cell structure, whereas the
turbulent solution with C3e=1.44 exhibits a 2-cell structure with
small counter-rotating cells at the bottom corners of the
enclosure. For C3e=1.35, the 2-cell structure breaks down into a
more complex flow pattern which results in closer agreement with
the heat transfer data. No grid-dependency studies have been
carried out.
BOOLEAN(LTURB)
MESG( Do you want a turbulent simulation? (y/n)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
+ LTURB=T
+ MESG( Turbulent simulation
ELSE
+ LTURB=F
+ MESG( Laminar simulation
ENDIF
REAL(TCOLD,THOT,TDIF,CP,YLEN,HEIGHT,MU,BETA,TREF,AGRAV,GRAS,GPRL)
REAL(DTF,VBUOY,RAY,RHO,TBRUNT,KOND,GNUSS,GCE3); CP=1.008E3
IF(LTURB) THEN
+ TDIF=25.0;HEIGHT=0.1
ELSE
+ TDIF=5.0;HEIGHT=0.05
ENDIF
YLEN=5.*HEIGHT;TCOLD=30.0+273.;THOT=TCOLD+TDIF
TREF=0.5*(THOT+TCOLD)
AGRAV=-9.81;BETA=1./TREF; RHO=353.505/TREF;MU=2.0383E-5
GROUP 2. Transience; time-step specification
GROUP 3. X-direction grid specification
** for a fine mesh use NX=NY=56 and ALF=3.5
REAL(ALF,GAM,ALFM)
IF(LTURB) THEN
+ NX=28;NY=48;ALF=2.5;ALFM=-ALF
ELSE
+ NX=28;NY=42;ALF=1.5;ALFM=-ALF
ENDIF
DO JJ=1,NX
+ GAM=ALF*(2.*(JJ/NX)-1.)
+ XFRAC(JJ)=-0.5*TANH(GAM)/TANH(ALFM)+0.5
+ XFRAC(JJ)=XFRAC(JJ)*HEIGHT
ENDDO
group 4. y-direction grid specification
DO JJ=1,NY
+ GAM=ALF*(2.*(JJ/NY)-1.)
+ YFRAC(JJ)=-0.5*TANH(GAM)/TANH(ALFM)+0.5
+ YFRAC(JJ)=YFRAC(JJ)*YLEN
ENDDO
GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1,TEM1);STORE(PRPS)
IF(LTURB) THEN
+ TURMOD(KEMODL-LOWRE);STORE(ENUT);PRT(TEM1)=0.9
ENDIF
group 8. terms (in differential equations) & devices
TERMS(TEM1,N,P,P,P,P,P)
group 9. properties of the medium (or media)
RHO1=RHO;KOND=3.5E-3+7.58E-5*TREF;ENUL=MU/RHO
PRT(TEM1)=0.9
RAY=-AGRAV*BETA*HEIGHT**3*TDIF*PRNDTL(H1)/(MU/RHO)**2
PRNDTL(TEM1)=-KOND
CP1=CP
GPRL=CP1*MU/KOND
GRAS=RAY/GPRL
GROUP 11. Initialization of variable or porosity fields
VBUOY=(-2.*BETA*0.5*TDIF*AGRAV*0.1*HEIGHT)**0.5
IF(LTURB) THEN
+ FIINIT(ENUT)=40.*MU/RHO
+ FIINIT(KE)=(0.1*VBUOY)**2
+ FIINIT(EP)=0.09*FIINIT(KE)**2/FIINIT(ENUT)
ENDIF
FIINIT(TEM1)=0.5*(THOT+TCOLD)
group 13. boundary conditions and special sources
1. Hot wall boundary: Constant Temperature of 65.8 deg C.
WALL (HOT,WEST,1,1,1,NY,1,NZ,1,1)
COVAL(HOT,TEM1,LOGLAW,THOT)
2. Cold wall boundary: Constant Temperature of 20 deg C.
WALL (COLD,EAST,NX,NX,1,NY,1,NZ,1,1)
COVAL(COLD,TEM1,LOGLAW,TCOLD)
3. Side wall boundary: adiabatic but with friction
WALL(TOPWAL,NORTH,1,NX,NY,NY,1,NZ,1,1)
4. Side wall boundary: adiabatic but with friction
WALL(BOTWAL,SOUTH,1,NX,1,1,1,NZ,1,1)
5. Buoyancy force
BUOYA=AGRAV;BUOYB=0.;BUOYC=0.
PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,1)
COVAL(BUOY,U1,FIXFLU,BOUSS)
BUOYD=BETA; BUOYE=TREF
6. Reference pressure at the centre of the cavity
PATCH(REFP,CELL,NX/2,NX/2,NY/2,NY/2,1,NZ,1,1)
COVAL(REFP,P1,FIXP,0.0); COVAL(REFP,V1,ONLYMS,0.0)
COVAL(REFP,TEM1,ONLYMS,SAME)
7. Buoyancy terms in the k-e model
IF(LTURB) THEN
+ PATCH(KEBUOY,PHASEM,1,NX,1,NY,1,NZ,1,1)
+ COVAL(KEBUOY,KE,KESOURCE,KESOURCE)
+ COVAL(KEBUOY,EP,KESOURCE,KESOURCE)
MESG( Enter value of turbulence buoyancy constant C3e
MESG( Default 1.44
READVDU(GCE3,REAL,1.44)
The GETSDR is in gxgenk.htm
+ SPEDAT(SET,KEBUOY,GCEB,R,GCE3)
ENDIF
GROUP 16. Termination of iterations
LITER(V1)=2;LITER(U1)=2;LITER(TEM1)=50
IF(LTURB) THEN
+ LITER(KE)=2;LITER(EP)=2
ENDIF
GROUP 17. Under-relaxation devices
VBUOY=(-BETA*TDIF*AGRAV*HEIGHT)**0.5; TBRUNT=HEIGHT/VBUOY
DTF=TBRUNT/NX
RELAX(V1,FALSDT,5.*DTF); RELAX(U1,FALSDT,5.*DTF)
IF(LTURB) THEN
+ RELAX(KE,FALSDT,2.*DTF); RELAX(EP,FALSDT,2.*DTF)
+ RELAX(TEM1,FALSDT,50.*DTF)
** for ce3=1.35 only 2000 sweeps are needed
+ KELIN = 3 ; LSWEEP = 3500 ; NPLT = 100
** Jakob (1949) correlation for air with Gr > 4.E5
+ GNUSS = 0.068*GRAS**(1./3.)
ELSE
+ LSWEEP = 800 ; NPLT = 50
** Jakob (1949) correlation for air with 1E4 < Gr < 4.E5
+ GNUSS = 0.195*GRAS**(1./4.)
ENDIF
GRAS
RAY
GNUSS
GROUP 18. Limits on variables or increments to them
VARMIN(TEM1)=TCOLD; VARMAX(TEM1)=THOT
GROUP 22. Spot-value print-out
IYMON=NY/2; IXMON=5; NYPRIN=NY/5; NXPRIN=NX/5
TSTSWP=-1;ITABL=3
GROUP 23. Field print-out and plot control
DISTIL=T
IF(LTURB) THEN
+ EX(P1)=2.532e-03;EX(U1)=1.009E-02;EX(V1)=2.360E-02;EX(KE)=5.328E-4
+ EX(EP)=4.214E-04;EX(EPKE)=3.548E+05;EX(TEM1)=3.145E+02
ELSE
+ EX(P1)=3.703E-04;EX(U1)=6.746E-03
+ EX(V1)=1.186E-02;EX(TEM1)=3.053E+02
ENDIF
libref=297