#$r002 PHOTON USE AUTOPLOT file PHI 5 cl;d 1 h1;d 1 ha;col3 1;blb4 2;redr msg temperature profile; pressto continue ENDUSE TEXT(1D RADIATIVE EQUILIBRIUM IN SLAB : 121 #cls TITLE DISPLAY The problem considered is radiative heat transfer in a stationary emitting and absorbing gray medium bounded by two infinite, plane parallel walls. The boundary surfaces at y=0 and y=L are kept at fixed temperatures Ts and Tn. The energy transfer is by pure radiation, so that the energy equation is simply: -d/dy (Qr) = 0 where Qr is the radiative heat transfer per unit area. Thermal radiation is modelled by solving an the equation for the radiosity R, as follows: d/dy ( 1/(a+s) dR/dy ) + 4*a*(E - R) = 0 where a is the absorption coefficient, s is the scattering coefficient, and E is the black-body emissive power. The problem is to solve for the temperature distribution, and hence Qr, given Ts, Tn, the optical thickness Kr*L; and the wall emissivities emws and emwn. Kr is the Rosseland mean absorption coefficient which is given by: Kr=(a+s). #pause This problem has been solved by Deissler [ASME J.Heat Transfer P241-246, (1964)] who used the Diffusion approximation with jump boundary conditions to obtain the following solutions: Qr = (Ews - Ewn)/(0.75*Kr*L + 1./emws + 1./emwn - 1.0) (Ews - Egs)/Qrad = 1./emws - 0.5 (Egn - Ewn)/Qrad = 1./emwn - 0.5 E = Egs - 0.75*Kr*Qr*y where Ews and Ewn are the emissive powers at the wall, and Egs and Egn are the emissive powers in the gas at the wall. The results of the radiosity model show excellent agreement with Deissler's results for the whole range of optical thickness considered by Deissler, i.e. 0 < Kr*L < 10.0. ENDDIS #pause REAL(EMWS,EMWN,TWS,TWN,EWS,EWN,ACON,GY) REAL(QRAD,KRAD,OTHICK,LENGTH,QRADA,TGNA,TGSA,EGS,EG,TA) CHAR(CH1);INTEGER(JJM1) MESG( Enter optical thickness 0 < Kr*L << 10.0 (default 5) READVDU(OTHICK,REAL,5.0) LENGTH=1.0;KRAD=OTHICK/LENGTH SCATT=0.0;EMISS=KRAD-SCATT EMWS=1.0;EMWN=1.0;TWS=1500.0;TWN=1000.0 EWS=SIGMA*TWS**4;EWN=SIGMA*TWN**4;QRAD=EWS-EWN ** analytical solution QRADA=1.0/(0.75*OTHICK+1./EMWS+1./EMWN-1.0) QRADA=QRADA*(EWS-EWN) TGSA =((EWS-QRADA*(1./EMWS-0.5))/SIGMA)**0.25 TGNA =((EWN+QRADA*(1./EMWN-0.5))/SIGMA)**0.25 GROUP 3,4,5. X,Y,Z-direction grid specification GRDPWR(Y,50,LENGTH,1.0) GROUP 6. Body-fitted coordinates or grid distortion GROUP 7. Variables stored, solved & named CP1=1.0 MESG( Enter required energy variable ? (TEM1 or H1) READVDU(CH1,CHAR,H1) IF(:CH1:.EQ.TEM1) THEN + MESG( TEM1 solution selected ELSE + MESG( H1 solution selected + TMP1=LINH;TMP1B=1.0/CP1 ENDIF RADIAT(RADI,EMISS,SCATT,:CH1:);STORE(EMPO) GROUP 8. Terms (in differential equations) & devices ** Deactive conduction & any built-in sources TERMS(:CH1:,N,N,N,N,P,P) GROUP 9. Properties of the medium (or media) GROUP 10. Inter-phase-transfer processes and properties GROUP 11. Initialization of variable or porosity fields FIINIT(:CH1:)=TWS;FIINIT(SRAD)=EWS ** analytical solution STORE(HA);ACON=0.75*KRAD*QRADA;EGS=SIGMA*TGSA**4 DO JJ=1,NY +PATCH(IN:JJ:,INIVAL,1,NX,JJ,JJ,1,NZ,1,1) +GY=0.5*YFRAC(JJ) IF(JJ.NE.1) THEN +JJM1=JJ-1;GY=YFRAC(JJM1)+0.5*(YFRAC(JJ)-YFRAC(JJM1)) ENDIF +GY=GY*YVLAST;EG=EGS-ACON*GY +TA=(EG/SIGMA)**0.25;INIT(IN:JJ:,HA,ZERO,TA) ENDDO GROUP 13. Boundary conditions and special sources ** Net radiation flux from wall PATCH(WALLRA,SOUTH,1,NX,1,1,1,NZ,1,1) COVAL(WALLRA,SRAD,2.*EMWS/(2.0-EMWS),EWS) PATCH(WALLRB,NORTH,1,NX,NY,NY,1,NZ,1,1) COVAL(WALLRB,SRAD,2.*EMWN/(2.0-EMWN),EWN) GROUP 15. Termination of sweeps LSWEEP=500;LITC=5 GROUP 16. Termination of iterations SELREF=F;RESREF(:CH1:)=1.E-4*QRAD;RESREF(SRAD)=RESREF(:CH1:) GROUP 17. Under-relaxation devices RELAX(:CH1:,LINRLX,0.3);RELAX(SRAD,LINRLX,0.6) GROUP 18. Limits on variables or increments to them VARMIN(:CH1:)=0.5*TWS GROUPS 19 to 21 GROUP 22. Spot-value print-out IYMON=1;NPLT=20;NYPRIN=1;ITABL=3 GROUP 23. Field print-out and plot control OUTPUT(:CH1:,Y,N,N,Y,Y,Y);PATCH(YWISE,PROFIL,1,1,1,NY,1,1,1,1) PLOT(YWISE,:CH1:,0.0,0.0);PLOT(YWISE,SRAD,0.0,0.0) GROUP 24. Dumps for restarts TSTSWP=-1