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