PHOTON USE
p;;;;
view z; up x
msg velocity vectors
vec z m sh; pause
msg temperature contours
con tem1 z m fi;0.001; Pause
ENDUSE
GROUP 1. RUN TITLE AND OTHER PRELIMINARIES
TEXT(3D Laminar Free Convection In A Cavity
TITLE
DISPLAY
This case, created by jzw and mrm at CHAM in August 2003,
concerns laminar buoyancy-driven air flow in a cubical cavity
which has the side-length of 0.127m.
The Rayleigh number, RAY is set to 4E4 for this calculation.
The user may change RAY to 1E5, 1E6 or 1E7.
The tilt angle, ANG is set to 90 as the default setting. The user
may change ANG to 45 or 0 so that the results can be compared
with experimental data.
The uniform grid of 15x15x7 is employed. It can be increased to
30x30x14 or 60x60x28 in order to produce more accurate results.
Specific Heat and conductivity are specified by a polynomial
formula via InForm.
Viscosity is calculated by Sutherland's law.
Density is calculated from Rayleigh number.
THOT(= 307 K) and TCOLD(= 300 K) are the vertical-wall
temperatures.
A linear profile of temperature between THOT and TCOLD is set
for the side-walls. A symmetry condition is applied to the high
boundary.
The Boussinesq approximation is employed with variable physical
properties.
Nusselt number based on the cold wall is calculated using
In-Form.
This problem has been studied experimentally and numerically
by Leong et al (Int. J. of Heat & Mass Transfer 42 1999).
The experimental data are as follows:
Nusselt number
Tilt angle
Ray 0 45 90
4E4 2.018 2.650 2.337
1E5 3.509 3.492 3.097
1E7 15.38 17.50 12.98
ENDDIS
REAL(TCOLD,THOT,TDIF,CP,HEIGHT,WIDTH,MU,BETA,TREF,AGRAV)
REAL(HCOLD,HHOT,HREF,DTF,VBUOY,RAY,RHO,TBRUNT,KOND,DEPTH)
REAL(LREF3,ANG,TWALL,PI)
REAL(AK1,AK2,AK3,AK4,AC1,AC2,AC3,AC4,AC5,TREF2,TREF3,TREF4)
HEIGHT=0.1272;WIDTH=HEIGHT;DEPTH=HEIGHT;LREF3=HEIGHT**3
TCOLD=27.+273.;TDIF=7.;THOT=TCOLD+TDIF
THOT
TCOLD
TREF=0.5*(THOT+TCOLD)
AGRAV=-9.81;BETA=1./TREF
** thermal conductivity for Rayleigh number
TREF2=TREF*TREF;TREF3=TREF2*TREF;TREF4=TREF3*TREF
** the following constants are supplied by CGCRI India
AK1=1.5207E-11;AK2=-4.8574E-8;AK3=1.0184E-4;AK4=-3.9333E-4
KOND=AK1*TREF3 + AK2*TREF2 + AK3*TREF + AK4
** set RG(1) for Nusselt number computed by In-Form
RG(1)=KOND
KOND
** specific heat for Rayleigh number
** the following constants are supplied by CGCRI India
AC1=1.9327E-10;AC2=-7.9999E-7;AC3=1.1407E-3;AC4=-4.489E-1
AC5=1.0575e3
CP=AC1*TREF4 + AC2*TREF3 + AC3*TREF2 + AC4*TREF + AC5
** dynamic viscosity from Sutherland's law
ENULA=1.458E-6;ENULB=110.4
MU=ENULA*TREF**1.5/(TREF+ENULB)
MU
CP
** Define Rayleigh number
RAY=4.E4
** Define tilt angle 90 deg for vertical heating
0 deg for bottom heating
ANG=90.0
PI=3.1415926;ANG=ANG*PI/180.
ANG
** Compute density for given Rayleigh number
RHO=(RAY*MU*KOND/(-AGRAV*BETA*TDIF*CP*LREF3))**0.5
RHO
** Implied press0 (but not used in simulations)
PRESS0=RHO*286.7*TREF
PRESS0
** check on Rayleigh number
RAY=-AGRAV*BETA*LREF3*TDIF*CP*RHO**2/(MU*KOND)
RAY
GROUP 2. Transience; time-step specification
GROUP 3. X-direction grid specification
group 4. y-direction grid specification
** As in paper of Leong et al (1999) solve for one-half
of cavity in the z-direction
NX=15;NY=15;NZ=7
GRDPWR(X,NX,WIDTH,1.0)
GRDPWR(Y,NY,HEIGHT,1.0)
GRDPWR(Z,NZ,0.5*DEPTH,1.0)
group 7. variables stored, solved & named
SOLVE(P1,U1,V1,W1,TEM1);SOLUTN(P1,Y,Y,Y,P,P,P)
SOLUTN(TEM1,Y,Y,Y,P,P,P)
STORE(KOND,SPH1,VISL,QWAL,NUSS,HTCB,PART)
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
** variable thermal conductivity at present
INFORM9BEGIN
(PROPERTY PRNDTL(TEM1) IS POL3(TEM1,-3.9333E-4,1.0184E-4,-4.8574E-8$
,1.5207E-11))
** variable specific heat at present
(PROPERTY CP1 IS POL4(TEM1, AC5, AC4, AC3, AC2, AC1) )
INFORM9END
** variable kinematic viscosity from Sutherlands Law
ENUL=SUTHLAND;ENULA=1.458E-6;ENULB=110.4
GROUP 11. Initialization of variable or porosity fields
VBUOY=(-2.*BETA*0.5*TDIF*AGRAV*0.1*HEIGHT)**0.5
VBUOY
FIINIT(TEM1)=0.5*(THOT+TCOLD)
group 13. boundary conditions and special sources
** Hot wall boundary
WALL (HOT,WEST,1,1,1,NY,1,NZ,1,LSTEP)
COVAL(HOT,TEM1,1.0,THOT)
** Cold wall boundary
WALL (COLD,EAST,NX,NX,1,NY,1,NZ,1,LSTEP)
COVAL(COLD,TEM1,1.0,TCOLD)
** N/S/L wall boundaries - linear temperature variation
REAL(XP,XM)
XM=0.0
DO JJ=1,NX ! exemplifies use of PIL DO loop to set wall temps
+ XP=0.5*(XM+XFRAC(JJ)) ! however InForm could provide the same
+ TWALL=THOT-TDIF*XP ! input more neatly
+ XM=XFRAC(JJ)
+ WALL(NWAL:JJ:,NORTH,:JJ:,:JJ:,NY,NY,1,NZ,1,LSTEP)
+ COVAL(NWAL:JJ:,TEM1,1.0,TWALL)
+ WALL(SWAL:JJ:,SOUTH,:JJ:,:JJ:,1,1,1,NZ,1,LSTEP)
+ COVAL(SWAL:JJ:,TEM1,1.0,TWALL)
+ WALL(LWAL:JJ:,LOW,:JJ:,:JJ:,1,NY,1,1,1,LSTEP)
+ COVAL(LWAL:JJ:,TEM1,1.0,TWALL)
ENDDO
** Buoyancy force
** assume uniform thermal expansion coefficient
BUOYA=AGRAV*COS(ANG)
BUOYB=AGRAV*SIN(ANG)
BUOYC=0.
PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
COVAL(BUOY,U1,FIXFLU,BOUSS)
COVAL(BUOY,V1,FIXFLU,BOUSS)
DVO1DT=BETA; BUOYE=-BETA*TREF
** Reference pressure at the centre of the cavity
PATCH(REFP,CELL,NX/2,NX/2,NY/2,NY/2,NZ/2,NZ/2,1,LSTEP)
COVAL(REFP,P1,FIXP,0.0); COVAL(REFP,V1,ONLYMS,0.0)
COVAL(REFP,TEM1,ONLYMS,SAME)
GROUP 16. Termination of iterations
LITER(TEM1)=50
GROUP 17. Under-relaxation devices
VBUOY=(-BETA*TDIF*AGRAV*HEIGHT)**0.5; TBRUNT=HEIGHT/VBUOY
DTF=TBRUNT/NY
RELAX(V1,FALSDT,DTF); RELAX(U1,FALSDT,DTF)
RELAX(W1,FALSDT,DTF)
RELAX(KE,FALSDT,DTF); RELAX(EP,FALSDT,DTF)
RELAX(TEM1,FALSDT,50.*DTF)
LSWEEP=300
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=1
NXPRIN=1; IXPRF=NX-3;IXPRL=NX
NPLT=20
TSTSWP=-1;YPLS=T;WALPRN=T; ITABL=3
GROUP 23. Field print-out and plot control
OUTPUT(PART,N,N,N,N,N,N)
OUTPUT(HTCB,N,N,N,N,N,N)
OUTPUT(SPH1,N,N,N,N,N,N)
** Compute Nusselt number distribution along the cold wall
on the last sweep
INFORM19BEGIN
(STORED of PART at COLD is (KOND/(0.5*(HEIGHT/NY))) with IF(ISWEEP.$
EQ.LSWEEP))
** compute surface heat flux, using function SUM
(STORED OF QWAL at COLD is SUM(PART*(TEM1-TCOLD))/(NZ*NY) with IF(I$
SWEEP.EQ.LSWEEP))
** compute bulk heat-transfer coefficient
(STORED OF HTCB at COLD is ABS(QWAL/TDIF) with IF(ISWEEP.EQ.LSWEEP))
** compute average nusselt number;
** RG(1) is the KOND based on 0.5*(THOT+TCOLD)
(STORED OF NUSS at COLD is HTCB*:XULAST:/:RG(1): with IF(ISWEEP.EQ.$
LSWEEP))
INFORM19END