#$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 RADIATION+HEAT SOURCE IN A SLAB:122 #cls TITLE DISPLAY The problem considered is that of radiative heat transfer in a stationary emitting and absorbing gray medium containing a uniformly-distributed volumetric heat source. The medium is bounded by two infinite, plane parallel walls at y=0 and y=L, which are kept at the same uniform temperatures. Thus, symmetry is exploited in the calculations. The energy transfer is by pure radiation, so that the energy equation is simply: -d/dy (Qr) + Qv = 0 where Qr is the radiative heat flux and Qv is the uniform volumetric heat generation rate in the medium. Thermal radiation is modelled by solving 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. #pause The problem is to solve for the temperature distribution given the wall temperature Tw, the optical thickness Kr*L; and the wall emissivity emw. Kr is the Rosseland mean absorption coefficient which is given by: Kr=(a+s). 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: Qw = 0.5*L*Qv (Egc - Ew)/Qw = (3./16)*Kr*L+1./emw-0.5+3./(4.*Kr*L) (Egw - Ew)/Qw = (1./emw-0.5+3.0/(4.*Kr*L)) Eg = Egc - (3./8)*Kr*Qv*y**2 where Qw is the wall heat flux, and Egc and Egw are the gas emissive powers at the centre line and wall. #pause ENDDIS REAL(EMWN,TWN,EWN,TA,GY,EG,EGC,EGW,ACON) REAL(QRAD,QVOL,KRAD,OTHICK,LENGTH,QRADA,TWNA,QWNA,TGCA,TGWA) 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 EMWN=1.0;TWN=1000.0 EWN=SIGMA*TWN**4;QRAD=EWN;QVOL=EWN/LENGTH ** analytical solution QWNA=QVOL*0.5*LENGTH EGC=EWN+QWNA*(3.*OTHICK/16.+1./EMWN-0.5+3./(4.*OTHICK)) EGW=EWN+QWNA*(1./EMWN-0.5+3./(4.*OTHICK)) TGCA=(EGC/SIGMA)**0.25;TGWA=(EGW/SIGMA)**0.25 ACON=QWNA/(EGC-EWN) GROUP 3,4,5. X,Y,Z-direction grid specification GRDPWR(Y,50,0.5*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 ** Deactivate 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:)=TWN;FIINIT(SRAD)=EWN ** analytical solution STORE(HA);ACON=3.0*KRAD*QVOL/8.0 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=EGC-ACON*GY*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(WALLRB,NORTH,1,NX,NY,NY,1,NZ,1,1) COVAL(WALLRB,SRAD,2.*EMWN/(2.0-EMWN),EWN) ** uniformly-distributed volumetric heat source PATCH(QHEAT,VOLUME,1,NX,1,NY,1,NZ,1,1) COVAL(QHEAT,:CH1:,FIXFLU,QVOL) 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*TWN GROUP 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