c
c file name GXRADIAT.HTM 160504
C**** SUBROUTINE GXRADI provides coding for the composite-flux
C and radiosity models for describing thermal radiation in a
C participating media. The temperature may be solved either directly
C (via TEM1) and indirectly (via H1) .
C
C Subroutine GXRADI is called from the following Groups of
C subroutine GREX3:-
C
C * Group 1, Section 1: in order to define various constants
C in the radiation models and to allocate any storage.
C
C * Group 9, Section 7: in order to define the diffusivity for
C the diffusion fluxes in the equations for the composite fluxes
C and radiosity.
C
C * Group 13, Section 13: in order to introduce the source terms
C in the energy equation (TEM1 or H1) and the equations for the
C composite fluxes and the radiosity.
C
C How the subroutine is activated:
C
C In the Q1 file, the user sets the following PIL commands:
C RADIAT(MODEL,ABSORB,SCAT,ENERGY)
C where MODEL = FLUX for the composite-flux model or
C = RADI for the radiosity model;
C ABSORB = absorption coefficient in m**-1 (=RADIA)
C SCAT = scattering coefficient in m**-1 (=RADIB)
C ENERGY = H1 for energy eqn solved via enthalpy or
C = TEM1 for energy eqn solved via temperature.
C
C Limitations:
C
C The radiation models are coded for ONEPHS=T only, and a gray
C medium is presumed. The composite-flux model is not valid for
C BFC=T, although the radiosity model can be used in such
C applications.
C
C For further information see PHOENICS Encyclopedia entry:
C 'Radiation Heat Transfer'
C
SUBROUTINE GXRADI(CARTES,TEMP1,DEN1,RHO1)
C----------------------------------------------------------------
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
COMMON/IDATA/NX,NY,NZ,LUPR1,IGFILL(116)
COMMON/RDATA/TINY,RDFIL1(38),HUNIT,RDFIL2(45)
COMMON/GENI/NXNY,IGFIL(48),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,
1 IPRPS,IRADX,IRADY,IRADZ,IVFOL
LOGICAL CARTES,NEZ
INTEGER TEMP1,DEN1,TEMVAR
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
common/dbe/dbear
logical dbear
SAVE CORAD,COTEM,FACRAD,VALTEM,IEMPO,ISRAD,NFLX,SIGMA,TEMVAR
C
NAMSUB = 'GXRADI'
IF(IGR.EQ.1) THEN ! Group 1
C.... SIGMA is the Stefan-Boltzmann constant in watts m**-2 K**-4
SIGMA = 5.6697E-8
TEMVAR= H1
IF(SOLVE(ITEM1)) TEMVAR=ITEM1
CALL SUB2(ISRAD,LBNAME('SRAD'),IEMPO,LBNAME('EMPO'))
IF(ISRAD.GT.0) THEN ! radiosity model
CALL SUB3R(COTEM,4.*RADIA*SIGMA,VALTEM,4.*RADIA,
1 CORAD,4.*RADIA)
ELSE ! composite-flux (i.e. six-flux) model
NFLX=0 ! count number of flux-pairs, nflx
IF(SOLVE(IRADX)) NFLX = NFLX + 1
IF(SOLVE(IRADY)) NFLX = NFLX + 1
IF(SOLVE(IRADZ)) NFLX = NFLX + 1
IF(NFLX.EQ.0) THEN
WRITE(LUPR1,*)'GXRADI has been called because'
WRITE(LUPR1,*)'a patch named RADI.. exists, but'
WRITE(LUPR1,*)'no radiation flux is being solved'
CALL SET_ERR(227,'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
CORAD = RADIA + FLOAT(NFLX-1)*RADIB/FLOAT(NFLX)
COTEM = FLOAT(2*NFLX)*RADIA*SIGMA
IF(CORAD.GT.0.0) THEN
FACRAD = RADIB/(FLOAT(NFLX)*CORAD)
ELSE
FACRAD=0.0
ENDIF
ENDIF
ELSEIF(IGR.EQ.9) THEN ! Group 9 for exchange coefficient
IF(INDVAR.EQ.IRADX.OR.INDVAR.EQ.IRADY.OR.INDVAR.EQ.IRADZ
1 .OR.INDVAR.EQ.ISRAD) THEN
C.... PIL variables RADIA and RADIB are the absorption coefficient and
c scattering coefficients respectively
CONST=RADIA+RADIB
IF(INDVAR.EQ.ISRAD.OR.CARTES.OR.INDVAR.NE.IRADY) THEN
RCONST=1.0/CONST
IF(ISRAD.GT.0) RCONST=4.0*RCONST/3.0
CALL FN1(LAMPR,RCONST) ! fn1(y,a) sets y = a
ELSE
CALL FN7(LAMPR,RV2D,0.0,1.0,CONST) ! fn7(y,x,a,b,c) sets
! y = (a + b*x)/(1.+c*x)
ENDIF
CALL FN27(LAMPR,DEN1)
ENDIF
ELSEIF(IGR.EQ.13) THEN ! Group 13: Sources for composite-flux variables
! and enthalpy or temperature
IF(INDVAR.EQ.IRADX.OR.INDVAR.EQ.IRADY.OR.INDVAR.EQ.IRADZ
1 .OR.INDVAR.EQ.ISRAD) THEN
IF(ISC.EQ.2) THEN ! CO
CALL FN1(CO,CORAD)
CALL FNZERO_IN_SLD_OR_POR(CO)
ELSEIF(ISC.EQ.13) THEN ! VAL
ITEMP=ITWO(TEMP1,TEMVAR,TEMVAR.EQ.H1)
IF(NEZ(TEMP0)) THEN ! put absolute temperature in EASP1
CALL FN2(EASP1,ITEMP,TEMP0,1.0)
ITEMP=EASP1
ENDIF
CALL FN63(VAL,ITEMP,SIGMA) ! fn63(y,x,a) sets y = a*x**4
IF(IEMPO.GT.0) THEN ! store emissive power 3D
CALL FN0(IEMPO,VAL)
ENDIF
IF(ISRAD.EQ.0.AND.NFLX.GT.1) THEN
CALL FN25(VAL,RADIA/(CORAD+TINY)) ! fn25(y,a) sets y = a*y
IF(INDVAR.EQ.IRADX) THEN
CALL SUB2(MRJ,IRADY,MRK,IRADZ)
ELSEIF(INDVAR.EQ.IRADY) THEN
CALL SUB2(MRJ,IRADX,MRK,IRADZ)
ELSE
CALL SUB2(MRJ,IRADX,MRK,IRADY)
ENDIF
IF(MRJ.GT.0) CALL FN34(VAL,MRJ,FACRAD) ! fn34(y,x,a) sets
! y = y + a * x
IF(MRK.GT.0) CALL FN34(VAL,MRK,FACRAD)
CALL L0F1(VAL,I,IADD,'GXRADI')
ENDIF
CALL FNZERO_IN_SLD_OR_POR(VAL)
ENDIF
ELSEIF(INDVAR.EQ.TEMVAR) THEN ! Thermal energy source from radiation
IF(TEMVAR.EQ.H1) THEN ! Enthalpy
IF(NEZ(TEMP0)) THEN ! Put absolute temperature in EASP1
CALL FN2(EASP1,TEMP1,TEMP0,1.0)
ITEMP=EASP1
ELSE ! temp0 = 0
ITEMP=TEMP1
ENDIF
IF(ISC.EQ.2) THEN ! CO
CALL FN22(ISPH1,1.E-10)
CALL FN28(CO,ISPH1,4.*COTEM*HUNIT)
CALL FN37(CO,ITEMP,3.0)
ELSEIF(ISC.EQ.13) THEN ! VAL
IF(ISRAD.GT.0) THEN
CALL FN2(VAL,ISRAD,0.0,VALTEM) ! fn2(y,x,a,b) sets y = a + b*x
ELSE
CALL FN1(VAL,0.0)
IF(IRADX.GT.0) CALL FN34(VAL,IRADX,2.*RADIA)
IF(IRADY.GT.0) CALL FN34(VAL,IRADY,2.*RADIA)
IF(IRADZ.GT.0) CALL FN34(VAL,IRADZ,2.*RADIA)
ENDIF
CALL L0F4(VAL,H1,ITEMP,CO,I,IH1,ITM,ICO,IADD,'GXRADI')
DO IIX = IXF,IXL
I = I + IADD
DO IIY = IYF,IYL
I = I + 1
TEM4= F(I+ITM)**4
F(I)= F(I+IH1) + HUNIT*(F(I)-COTEM*TEM4)/
1 (F(I+ICO) + TINY)
ENDDO
ENDDO
CALL FNZERO_IN_SLD_OR_POR(VAL)
ENDIF
ELSE ! true temperature
IF(NEZ(TEMP0)) THEN
CALL FN2(EASP1,TEMVAR,TEMP0,1.0)
ITEMP=EASP1
ELSE
ITEMP=TEMVAR
ENDIF
IF(ISC.EQ.2) THEN ! CO
CALL FN61(CO,ITEMP,COTEM)
CALL FNZERO_IN_SLD_OR_POR(CO)
ELSEIF(ISC.EQ.13) THEN ! VAL
IF(ISRAD.GT.0) THEN
CALL FN2(VAL,ISRAD,0.0,VALTEM)
ELSE
CALL FN1(VAL,0.0)
IF(IRADX.GT.0) CALL FN34(VAL,IRADX,2.*RADIA)
IF(IRADY.GT.0) CALL FN34(VAL,IRADY,2.*RADIA)
IF(IRADZ.GT.0) CALL FN34(VAL,IRADZ,2.*RADIA)
ENDIF
CALL FN62(VAL,ITEMP,1.0/(COTEM+TINY))! FN62(Y,X,A) sets
! Y = A * Y/X**3
IF(NEZ(TEMP0)) CALL FN33(VAL,-TEMP0) ! FN33(Y,A) sets Y = Y + A
ENDIF
CALL FNZERO_IN_SLD_OR_POR(VAL)
ENDIF
ENDIF
ENDIF
NAMSUB = 'gxradi'
END
c