radi.f90

FDS为火灾动力学模拟软件源代码,该软件为开源项目,代码语言主要为FORTRAN,可在WINDOWS和LINUX下编译运行,详细说明可参考http://fire.nist.gov/fds/官方网址

MODULE RAD
 
! Radiation heat transfer
 
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_VARIABLES
USE RADCONS

IMPLICIT NONE
PRIVATE

CHARACTER(255), PARAMETER :: radiid='$Id: radi.f90 712 2007-09-28 20:20:09Z drjfloyd $'
CHARACTER(255), PARAMETER :: radirev='$Revision: 712 $'
CHARACTER(255), PARAMETER :: radidate='$Date: 2007-09-28 16:20:09 -0400 (Fri, 28 Sep 2007) $'

PUBLIC INIT_RADIATION,COMPUTE_RADIATION,NSB,NRA,UIIDIM,NRT,RSA,NRP,GET_REV_radi
 
CONTAINS
 
 
SUBROUTINE INIT_RADIATION
USE MEMORY_FUNCTIONS, ONLY : CHKMEMERR
USE MIEV
USE RADCALV
REAL(EB) :: THETAUP,THETALOW,PHIUP,PHILOW,F_THETA,MW_RADCAL,PLANCK_C2,KSI,LT, &
            ZZ,RCRHO,YY,BBF,AP0,AMEAN,MTOT,XLENG,YLENG,ZLENG
INTEGER  :: N,I,J,K,IZERO,NN,NI,II,JJ,IIM,JJM,IBND,NS,I_RADCAL,IZ
TYPE(SPECIES_TYPE), POINTER :: SS
TYPE(REACTION_TYPE), POINTER :: RN
TYPE(PARTICLE_CLASS_TYPE), POINTER :: PC
 
! Determine the number of polar angles (theta)
 
NRA = NUMBER_RADIATION_ANGLES
IF (CYLINDRICAL) THEN
   NRT = NINT(SQRT(REAL(NRA)))
ELSEIF (TWO_D) THEN
   NRT = 1
ELSE
   NRT = 2*NINT(0.5_EB*1.17*REAL(NRA)**(1._EB/2.26))
ENDIF      
 
ALLOCATE(NRP(1:NRT),STAT=IZERO)
CALL ChkMemErr('INIT','NRP',IZERO)
 
! Determine number of azimuthal angles (phi)
 
N = 0
DO I=1,NRT
   IF (CYLINDRICAL) THEN
      NRP(I) = NINT(REAL(NRA)/(REAL(NRT)))
      ELSEIF (TWO_D) THEN
      NRP(I) = 4*NINT(0.25_EB*REAL(NRA))
      ELSE
      THETALOW = PI*REAL(I-1)/REAL(NRT)
      THETAUP  = PI*REAL(I)/REAL(NRT)
      NRP(I) = NINT(0.5_EB*REAL(NRA)*(COS(THETALOW)-COS(THETAUP)))
      NRP(I) = MAX(4,NRP(I))
      NRP(I) = 4*NINT(0.25_EB*REAL(NRP(I)))
      ENDIF
   N = N + NRP(I)
ENDDO
NRA = N
NUMBER_RADIATION_ANGLES = NRA
 
! Set the opening angle of the cylindrical geometry equal to the azimuthal angle
 
IF (CYLINDRICAL) DPHI0 = PI/REAL(NRP(1))
 
NSB = NUMBER_SPECTRAL_BANDS
 
ALLOCATE(RSA(1:NRA),STAT=IZERO)
CALL ChkMemErr('INIT','RSA',IZERO)
ALLOCATE(DLX(1:NRA),STAT=IZERO)
CALL ChkMemErr('INIT','DLX',IZERO)
ALLOCATE(DLY(1:NRA),STAT=IZERO)
CALL ChkMemErr('INIT','DLY',IZERO)
ALLOCATE(DLZ(1:NRA),STAT=IZERO)
CALL ChkMemErr('INIT','DLZ',IZERO)
IF (CYLINDRICAL) THEN
   ALLOCATE(DLB(1:NRA),STAT=IZERO)
   CALL ChkMemErr('INIT','DLB',IZERO)
ENDIF
ALLOCATE(DLN(-3:3,1:NRA),STAT=IZERO)
CALL ChkMemErr('INIT','DLN',IZERO)
ALLOCATE(DLM(1:NRA,3),STAT=IZERO)
CALL ChkMemErr('INIT','DLM',IZERO)
 
! Determine mean direction normals and sweeping orders
 
N = 0
DO I=1,NRT
   DO J=1,NRP(I)
      N = N + 1
      THETALOW  = PI*REAL(I-1)/REAL(NRT)
      THETAUP   = PI*REAL(I)/REAL(NRT)
      F_THETA   = 0.5_EB*(THETAUP-THETALOW  - COS(THETAUP)*SIN(THETAUP) + COS(THETALOW)*SIN(THETALOW))
      IF (CYLINDRICAL) THEN
         PHILOW = PI*REAL(J-1)/REAL(NRP(I))
         PHIUP  = PI*REAL(J)/REAL(NRP(I))
      ELSEIF (TWO_D) THEN
         PHILOW = TWOPI*REAL(J-1)/REAL(NRP(I)) + PIO2
         PHIUP  = TWOPI*REAL(J)/REAL(NRP(I))   + PIO2
      ELSE
         PHILOW = TWOPI*REAL(J-1)/REAL(NRP(I))
         PHIUP  = TWOPI*REAL(J)/REAL(NRP(I))
      ENDIF
      RSA(N) = (PHIUP-PHILOW)*(COS(THETALOW)-COS(THETAUP))
      IF (CYLINDRICAL) THEN
         DLX(N) = 2._EB*SIN(DPHI0/2.)*(SIN(PHIUP)-SIN(PHILOW)) *F_THETA
         DLY(N) =  (-SIN(DPHI0/2.)*(SIN(PHIUP)-SIN(PHILOW))  +COS(DPHI0/2.)*(COS(PHILOW)-COS(PHIUP)))*F_THETA
         DLB(N) =  (-SIN(DPHI0/2.)*(SIN(PHIUP)-SIN(PHILOW))  -COS(DPHI0/2.)*(COS(PHILOW)-COS(PHIUP)))*F_THETA
         DLZ(N)    = 0.5_EB*(PHIUP-PHILOW)   * ((SIN(THETAUP))**2-(SIN(THETALOW))**2)
      ELSEIF (TWO_D) THEN
         DLX(N) = (SIN(PHIUP)-SIN(PHILOW))*F_THETA
         DLY(N) = 0._EB
         DLZ(N) = (COS(PHILOW)-COS(PHIUP))*F_THETA
      ELSE
         DLX(N) = (SIN(PHIUP)-SIN(PHILOW))*F_THETA
         DLY(N) = (COS(PHILOW)-COS(PHIUP))*F_THETA
         DLZ(N)    = 0.5_EB*(PHIUP-PHILOW)      * ((SIN(THETAUP))**2-(SIN(THETALOW))**2)
      ENDIF
   ENDDO
ENDDO
 
! Set (wall normal)*(angle vector) value
 
DO N = 1,NRA
   DLN(-1,N) = -DLX(N)
   DLN( 1,N) =  DLX(N)
   DLN(-2,N) = -DLY(N)
   DLN( 2,N) =  DLY(N)
   DLN(-3,N) = -DLZ(N)
   DLN( 3,N) =  DLZ(N)
ENDDO
 
! Calculate mirroring matrix
 
N = 0
DO I=1,NRT
   DO J=1,NRP(I)
      N = N + 1
      DO K=1,3
         IF (TWO_D .AND. .NOT.CYLINDRICAL) THEN
            SELECT CASE(K)
               CASE(1)             ! X-surfaces
                  IIM = 1
                  JJM = NRP(I) - J + 1
               CASE(2)             ! Y-surfaces
                  IIM = 1
                  JJM = J
               CASE(3)             ! Z-surfaces
                  IIM = 1
                  JJM = NRP(I)/2 - J + 1
            END SELECT
            JJM = MODULO(JJM,NRP(I))
            IF (JJM==0) JJM = NRP(I)
         ELSE
            SELECT CASE(K)
               CASE(1)             ! X-surfaces
                  IIM = I
                  JJM = NRP(I)/2 - J + 1
               CASE(2)             ! Y-surfaces
                  IIM = I
                  JJM = NRP(I) - J + 1
               CASE(3)             ! Z-surfaces
                  IIM = NRT - I + 1
                  JJM = J
            END SELECT
            IIM = MODULO(IIM,NRT)
            JJM = MODULO(JJM,NRP(I))
            IF (IIM==0) IIM = NRT
            IF (JJM==0) JJM = NRP(I)
         ENDIF
 
         NN = 0
         DO II = 1,IIM
            DO JJ = 1,NRP(II)
               NN = NN + 1
               IF ((II==IIM).AND.(JJ==JJM)) NI = NN
            ENDDO
         ENDDO
         DLM(N,K) = NI
      ENDDO
   ENDDO
ENDDO
 
!-----------------------------------------------------
!
!            Spectral information
!
!-----------------------------------------------------

INIT_WIDE_BAND: IF (WIDE_BAND_MODEL) THEN
 
! Fraction of blackbody emission in a wavelength interval
 
   PLANCK_C2 = 14387.69_EB            ! micron.K
   NLAMBDAT  = 4000_EB
   LTSTEP    = 25.0_EB           ! maximum LAMBDA*T = NLANBDAT*LTSTEP
   ALLOCATE(BBFRAC(0:NLAMBDAT),STAT=IZERO)
   CALL ChkMemErr('INIT','BBFRAC',IZERO)
   BBFRAC = 0._EB
   LT = 0._EB
   DO I = 1,NLAMBDAT
      LT =  LT + LTSTEP
      KSI = PLANCK_C2/LT
      DO J = 1,50
         BBFRAC(I) = BBFRAC(I) + EXP(-KSI*REAL(J))/REAL(J) * (KSI**3 + 3.*KSI**2/REAL(J) + 6.*KSI/REAL(J)**2 + 6./REAL(J)**3)
      ENDDO
   ENDDO
   BBFRAC =  BBFRAC * 15._EB/PI**4
 
! Define band limit wave lengths in micrometers
 
   ALLOCATE(WL_LOW(1:NSB),STAT=IZERO)
   CALL ChkMemErr('INIT','WL_LOW',IZERO)
   ALLOCATE(WL_HIGH(1:NSB),STAT=IZERO)
   CALL ChkMemErr('INIT','WL_HIGH',IZERO)
   IF (CH4_BANDS) THEN
      WL_LOW(1:NSB) =(/1.00_EB, 2.63_EB, 2.94_EB, 3.57_EB, 4.17_EB, 4.6_EB, 7.00_EB, 8.62_EB, 10.0_EB /)
      WL_HIGH(1:NSB)=(/2.63_EB, 2.94_EB, 3.57_EB, 4.17_EB, 4.60_EB, 7.0_EB, 8.62_EB, 10.0_EB, 200._EB /) 
   ELSE
      WL_LOW(1:NSB) =(/1.00_EB, 2.63_EB, 2.94_EB, 4.17_EB, 4.6_EB, 10.0_EB /)
      WL_HIGH(1:NSB)=(/2.63_EB, 2.94_EB, 4.17_EB, 4.6_EB, 10.0_EB, 200.0_EB /)
   ENDIF
 
ENDIF INIT_WIDE_BAND
 
!----------------------------------------------------------------------------
!
!     Tables for gas phase absorption coefficient
!
!     CONTROLLING PROGRAM FOR SUBROUTINE "RADCAL", A NARROW-BAND
!     MODEL FOR CALCULATING SPECTRAL INTENSITY (W/M-2/SR/MICRON) AND
!     SPECTRAL TRANSMITTANCE VERSUS WAVELENGTH (MICRONS) IN A NONISO-
!     THERMAL, VARIABLE COMPOSITION  MIXTURE OF CO2, H2O, CO, N2, O2,
!     CH4, AND SOOT. FOR A HOMOGENEOUS PATH, THE PROGRAM ALSO COMPUTES
!     THE PLANCK-MEAN ABSORPTION COEF, AP0, THE INCIDENT-MEAN ABSORPTION
!     COEFFICIENT, AIWALL, AND THE EFFECTIVE-MEAN ABSORPTION COEFFICIENT,
!     AMEAN, ALL IN UNITS OF INVERSE METERS.
!
!     INPUT PARAMETERS:
!          NPT=NUMBER OF HOMOGENEOUS ELEMENTS
!          DD(J)=THICKNESS OF J TH ELEMENT, M
!          RCT(J)=TEMPERATURE OF J TH ELEMENT, K.
!          P(I,J)=PARTIAL PRESSURE OF GASEOUS COMPONENTS, kPa:
!                  I   GASEOUS SPECIES
!                  1        CO2
!                  2        H2O
!                  3        CH4
!                  4        CO
!                  5        O2
!                  6        N2
!          SVF(J)=SOOT VOLUME FRACTION OF J TH ELEMENT
!          OMMIN=MINIMUM WAVE NUMBER IN SPECTRUM, CM-1.
!          OMMAX=MAXIMUM WAVE NUMBER IN SPECTRUM, CM-1.
!
!-------------------------------------------------------------------------
 
CALL RCALLOC  ! Allocate arrays for RadCal
 
! 20% of mean beam length, Eq 8-51, Holman, 7th Ed. Heat Transfer. Length = 3.6*Volume/Area
 
XLENG = MESHES(1)%XF-MESHES(1)%XS
YLENG = MESHES(1)%YF-MESHES(1)%YS
ZLENG = MESHES(1)%ZF-MESHES(1)%ZS
IF (PATH_LENGTH < 0.0_EB) THEN  ! default was -1.0
   IF (TWO_D) THEN ! calculate based on the geometry
      PATH_LENGTH = MIN( 10._EB , 0.1_EB*3.6_EB*XLENG*ZLENG/(XLENG+ZLENG) )
   ELSE
      PATH_LENGTH = MIN( 10._EB , 0.1_EB*3.6_EB*XLENG*YLENG*ZLENG/(XLENG*YLENG+XLENG*ZLENG+YLENG*ZLENG) )
   ENDIF
ENDIF
DD(1) = MAX(PATH_LENGTH,1.0E-4_EB)
 
! Using RadCal, create look-up tables for the absorption coefficients for all gas species, mixture fraction or aerosols

SPECIES_LOOP: DO NS=1,N_SPECIES

   SS => SPECIES(NS)
   IF (.NOT. SS%ABSORBING) CYCLE SPECIES_LOOP
 
   GAS_TYPE: SELECT CASE (SS%MODE) 
 
      CASE (MIXTURE_FRACTION_SPECIES) GAS_TYPE
 
         RN => REACTION(SS%REAC_INDEX)
         SS%NKAP_TEMP = 21
         SS%NKAP_MASS = 40
         SS%MAXMASS   = RN%Z_F
         ALLOCATE(SS%KAPPA(0:SS%NKAP_MASS,0:SS%NKAP_TEMP,1:NSB),STAT=IZERO)
         CALL ChkMemErr('INIT','SS%KAPPA',IZERO)
         BAND_LOOP_MF: DO IBND = 1,NSB
            IF (NSB>1) THEN
               OMMIN = REAL(NINT(1.E4/WL_HIGH(IBND)),EB)
               OMMAX = REAL(NINT(1.E4/WL_LOW(IBND)),EB)
            ELSE
               OMMIN = 50.
               OMMAX = 10000.
            ENDIF
            CALL INIT_RADCAL      
            T_LOOP_MF: DO K = 0,SS%NKAP_TEMP
               RCT(1) = 300._EB + K*(2400._EB-300._EB)/SS%NKAP_TEMP
               Z_LOOP_MF: DO J=0,SS%NKAP_MASS
                  ZZ = IZ2ZZ(J,RN%Z_F,SS%NKAP_MASS)
                  IZ = MIN(10000,MAX(0,NINT(ZZ*10000._EB)))
                  RCRHO = SS%MW_MF(IZ)*P_INF/(R0*RCT(1))
                  MTOT = SS%Y_MF(IZ,FUEL_INDEX)/RN%MW_FUEL + SS%Y_MF(IZ,O2_INDEX)/MW_O2   + SS%Y_MF(IZ,N2_INDEX)/MW_N2 +  &
                         SS%Y_MF(IZ,H2O_INDEX)/MW_H2O   + SS%Y_MF(IZ,CO2_INDEX)/MW_CO2 + SS%Y_MF(IZ,CO_INDEX)/MW_CO
                  SPECIE(1) = SS%Y_MF(IZ,CO2_INDEX)
                  SPECIE(2) = SS%Y_MF(IZ,H2O_INDEX)
                  SPECIE(3) = SS%Y_MF(IZ,FUEL_INDEX)
                  SPECIE(4) = SS%Y_MF(IZ,CO_INDEX)
                  SPECIE(5) = SS%Y_MF(IZ,SOOT_INDEX)*RCRHO/RHO_SOOT
                  P(1,1) = SS%Y_MF(IZ,CO2_INDEX)/MW_CO2/MTOT
                  P(2,1) = SS%Y_MF(IZ,H2O_INDEX)/MW_H2O/MTOT
                  P(3,1) = SS%Y_MF(IZ,FUEL_INDEX)/RN%MW_FUEL/MTOT
                  P(4,1) = SS%Y_MF(IZ,CO_INDEX)/MW_CO/MTOT
                  P(5,1) = SS%Y_MF(IZ,O2_INDEX)/MW_O2/MTOT
                  P(6,1) = SS%Y_MF(IZ,N2_INDEX)/MW_N2/MTOT
                  SVF(1) = SS%Y_MF(IZ,SOOT_INDEX)*RCRHO/RHO_SOOT
                  CALL RADCAL(AMEAN,AP0)
                  IF (NSB==1 .AND. PATH_LENGTH > 0.0_EB) THEN
                     SS%KAPPA(J,K,IBND) = MIN(AMEAN,AP0)
                  ELSE
                     IF (NSB==1) THEN
                        BBF = 1._EB
                     ELSE
                        BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),RCT(1))
                     ENDIF
                     SS%KAPPA(J,K,IBND) = AP0/BBF
                  ENDIF
               ENDDO Z_LOOP_MF
            ENDDO T_LOOP_MF
         ENDDO BAND_LOOP_MF

      CASE (GAS_SPECIES) GAS_TYPE

         SS%NKAP_TEMP = 21
         SS%NKAP_MASS = 200
         SS%MAXMASS=1._EB
         ALLOCATE(SS%KAPPA(0:SS%NKAP_MASS,0:SS%NKAP_TEMP,1:NSB),STAT=IZERO)
         CALL ChkMemErr('RADI','KAPPA',IZERO)
         IF (NS==I_CO2) THEN
            I_RADCAL = 1
            MW_RADCAL = MW_CO2
         ELSEIF (NS==I_CO) THEN
            I_RADCAL = 4
            MW_RADCAL = MW_CO
         ELSEIF (NS==I_WATER) THEN
            I_RADCAL = 2
            MW_RADCAL = MW_H2O
         ELSE
            IF (SS%ABSORBING) THEN
               I_RADCAL = 3
               MW_RADCAL = 16
            ELSE
               I_RADCAL = 5