#$r002
   PHOTON USE
   AUTOPLOT
   file
   PHI 5
 
   cl;d 1 h1;d 1 ha;col3 1;blb4 2;redr
   msg    temperature profile; press  to continue
   ENDUSE
TEXT(1D RADIATION+HEAT SOURCE IN A TUBE:123
TITLE
  DISPLAY
   The problem considered is similar to case 122 above, except
   that the radiative heat transfer takes place in a tube with
   a symmetry plane at r=0 and an isothermal wall at r=D/2.
 
   The problem is to solve for the temperature distribution given
   the wall temperature Tw, the optical thickness Kr*D; and the wall
   emissivity emw. Kr is the Rosseland mean absorption coefficient
   which is given by: Kr=(a+s).
  ENDDIS
   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.25*D*Qv
 
     (Egc - Ew)/Qw = (3./16)*Kr*D+1./emw-0.5+9./(8.*Kr*D)
 
     (Egw - Ew)/Qw = (1./emw-0.5+9./(8.*Kr*D))
 
     Eg  = Egc - 0.75*Kr*Qw*r**2/D
 
   where Qw is the wall heat flux, and Egc and Egw are the gas
   emissive powers at the centre line and wall.
REAL(EMWN,TWN,EWN,TA,GY,EG,EGC,EGW,ACON)
REAL(QRAD,QVOL,KRAD,OTHICK,DIAM,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)
DIAM=1.0;KRAD=OTHICK/DIAM
SCATT=0.0;EMISS=KRAD-SCAT
EMWN=1.0;TWN=1000.0;EWN=SIGMA*TWN**4;QRAD=EWN;QVOL=EWN/DIAM
  ** analytical solution
QWNA=QVOL*0.25*DIAM
EGC=EWN+QWNA*(3.*OTHICK/16.+1./EMWN-0.5+9./(8.*OTHICK))
EGW=EWN+QWNA*(1./EMWN-0.5+9./(8.*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
CARTES=F;GRDPWR(Y,50,0.5*DIAM,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:)=TWN;FIINIT(SRAD)=EWN
  ** analytical solution
STORE(HA);ACON=0.75*KRAD*QWNA/DIAM
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. Data communicated by satellite to GROUND
    GROUP 20. Preliminary print-out
    GROUP 21. Print-out of variables
    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