radi.f90

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

               MW_RADCAL = 32
            ENDIF
         END IF
         BAND_LOOP_FR1: DO IBND = 1,NSB
            IF (NSB>1) THEN
               OMMIN = REAL(NINT(1.E4_EB/WL_HIGH(IBND)),EB)
               OMMAX = REAL(NINT(1.E4_EB/WL_LOW(IBND)),EB)
            ELSE
               OMMIN = 50._EB
               OMMAX = 10000._EB
            ENDIF
            CALL INIT_RADCAL
            T_LOOP_MF2: DO K = 0,SS%NKAP_TEMP
               RCT(1) = 300._EB + K*(2400._EB-300._EB)/SS%NKAP_TEMP
               Z_LOOP_MF2: DO J=0,SS%NKAP_MASS
                  YY = REAL(J)*SS%MAXMASS/REAL(SS%NKAP_MASS)
                  MTOT = YY/MW_RADCAL+(1._EB-YY)/MW_N2
                  SPECIE = 0._EB
                  SPECIE(I_RADCAL) = YY
                  P(:,1) = 0._EB
                  SVF = 0._EB
                  P(I_RADCAL,1) = YY/MW_RADCAL/MTOT
                  P(6,1) = (1-YY)/MW_N2/MTOT
                  CALL RADCAL(AMEAN,AP0)
                  IF (NSB==1 .AND. PATH_LENGTH > 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(1),WL_HIGH(1),RCT(1))
                     ENDIF
                     SS%KAPPA(J,K,IBND) = AP0/BBF
                  ENDIF
               ENDDO Z_LOOP_MF2
            ENDDO T_LOOP_MF2
         ENDDO BAND_LOOP_FR1

      CASE (AEROSOL_SPECIES) GAS_TYPE

         SS%NKAP_TEMP = 21
         SS%NKAP_MASS = 200  
         SS%MAXMASS=0.2_EB   
         ALLOCATE(SS%KAPPA(0:SS%NKAP_MASS,0:SS%NKAP_TEMP,1:NSB),STAT=IZERO)
         CALL ChkMemErr('RADI','KAPPA',IZERO)
         BAND_LOOP_FR2: DO IBND = 1,NSB
            IF (NSB>1) THEN
               OMMIN = REAL(NINT(1.E4_EB/WL_HIGH(IBND)),EB)
               OMMAX = REAL(NINT(1.E4_EB/WL_LOW(IBND)),EB)
            ELSE
               OMMIN = 50._EB
               OMMAX = 10000._EB
            ENDIF
            CALL INIT_RADCAL
            I_RADCAL=5
            T_LOOP_MF3: DO K = 0,SS%NKAP_TEMP
               RCT(1) = 300._EB + K*(2400._EB-300._EB)/SS%NKAP_TEMP
               Z_LOOP_MF3: DO J=0,SS%NKAP_MASS
                  YY = REAL(J)*SS%MAXMASS/REAL(SS%NKAP_MASS)
                  RCRHO = 29._EB * P_INF/(R0*RCT(1)) !Good enough
                  SPECIE = 0._EB               
                  P(:,1) = 0._EB               
                  P(6,1) = 1._EB
                  SPECIE(I_RADCAL) = YY*RCRHO/RHO_SOOT                  
                  SVF(1) = YY*RCRHO/RHO_SOOT                                    
                  CALL RADCAL(AMEAN,AP0)
                  IF (NSB==1 .AND. PATH_LENGTH > 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(1),WL_HIGH(1),RCT(1))
                     ENDIF
                     SS%KAPPA(J,K,IBND) = AP0/BBF
                  ENDIF
               ENDDO Z_LOOP_MF3
            ENDDO T_LOOP_MF3
         ENDDO BAND_LOOP_FR2

   END SELECT GAS_TYPE

ENDDO SPECIES_LOOP
 
CALL RCDEALLOC  ! Deallocate RadCal arrays
 
 
!-----------------------------------------------------
!
!     Tables for droplet absorption coefficients
! 
!-----------------------------------------------------

DROPLETS: IF (N_EVAP_INDICIES>0) THEN
   GET_PC_RADI: DO J=1,N_PART
      PC => PARTICLE_CLASS(J)
      IF ((.NOT. PC%FUEL) .AND. (.NOT. PC%WATER)) CYCLE GET_PC_RADI
      IF (PC%FUEL) THEN  ! Tuntomo, Tien and Park, Comb. Sci. and Tech., 84, 133-140, 1992
         PC%WQABS(:,1) = (/0,10,16,28,52,98,191,386,792,1629,3272,6163, &
            10389,15588,20807,23011,22123,22342,22200,22241,21856, &
            22795,23633,24427,25285,26207,27006,27728,28364,28866, &
            29260/)!,29552,29748,30000/)
         PC%KWR(0) = 0._EB
         DO I=1,NDG
            PC%KWR(I)=EXP(I/2.5_EB-4._EB)
         ENDDO
         PC%WQABS=PC%WQABS/10000._EB
         PC%KWR=PC%KWR/1000000._EB
      ENDIF
      IF (PC%WATER) CALL MEAN_CROSS_SECTIONS(J)
   ENDDO GET_PC_RADI
ENDIF DROPLETS
 
! A few miscellaneous constants

FOUR_SIGMA = 4._EB*SIGMA
RPI_SIGMA  = RPI*SIGMA
 
! In axially symmetric case, each angle represents two symmetric angles. So weight the intensities by two.

W_AXI = 1._EB
IF (CYLINDRICAL) W_AXI = 2._EB
 
END SUBROUTINE INIT_RADIATION


 
SUBROUTINE COMPUTE_RADIATION(NM)

! Call radiation routine or simply specify the radiative loss

USE MESH_POINTERS
USE COMP_FUNCTIONS, ONLY : SECOND  
REAL(EB) :: TNOW
INTEGER, INTENT(IN) :: NM 

IF (EVACUATION_ONLY(NM)) RETURN

TNOW=SECOND()
 
CALL POINT_TO_MESH(NM)

IF (RADIATION) THEN
   CALL RADIATION_FVM(NM)
ELSE
   IF (N_REACTIONS>0) QR = -RADIATIVE_FRACTION*Q
ENDIF

TUSED(9,NM)=TUSED(9,NM)+SECOND()-TNOW

CONTAINS 
 
SUBROUTINE RADIATION_FVM(NM)
USE MIEV
USE MATH_FUNCTIONS, ONLY : INTERPOLATE1D

REAL(EB) :: ZZ, RAP, AX, AXU, AXD, AY, AYU, AYD, AZ, VC, RU, RD, RP, &
            ILXU, ILYU, ILZU, QVAL, BBF, BBFA, NCSDROP, RSA_RAT, WAXIDLN , KAPPA_1, Z_2
INTEGER  :: N,IIG,JJG,KKG,I,J,K,IW,II,JJ,KK,IOR,IC,IWUP,IWDOWN, &
            ISTART, IEND, ISTEP, JSTART, JEND, JSTEP, &
            KSTART, KEND, KSTEP, NSTART, NEND, NSTEP, &
            I_UIID, N_UPDATES, IBND, TYY, NOM, NS, IBC,EVAP_INDEX
LOGICAL :: UPDATE_INTENSITY
REAL(EB), POINTER, DIMENSION(:,:,:) :: KFST4, IL, UIIOLD, KAPPAW, KFST4W, EXTCOE, SCAEFF
REAL(EB), POINTER, DIMENSION(:)     :: OUTRAD_W, INRAD_W
INTEGER, INTENT(IN) :: NM
TYPE (OMESH_TYPE), POINTER :: M2
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(PARTICLE_CLASS_TYPE), POINTER :: PC

KFST4    => WORK1
IL       => WORK2
UIIOLD   => WORK3
EXTCOE   => WORK4
KAPPAW   => WORK5
SCAEFF   => WORK6
KFST4W   => WORK7
OUTRAD_W => WALL_WORK1
INRAD_W  => WALL_WORK2
 

! Ratio of solid angle, used in scattering
 
RSA_RAT = 1._EB/(1._EB-1._EB/NRA)
 
! Check if it time to update radiation intensity field
 
RAD_CALL_COUNTER  = RAD_CALL_COUNTER+1
IF ( MOD(RAD_CALL_COUNTER,TIME_STEP_INCREMENT)==1 .OR. TIME_STEP_INCREMENT==1) THEN
   UPDATE_INTENSITY = .TRUE.
ELSE
   UPDATE_INTENSITY = .FALSE.
ENDIF

IF (WIDE_BAND_MODEL) QR = 0._EB
IF (UPDATE_INTENSITY) QRADIN = 0._EB
 
! Loop over spectral bands
 
BAND_LOOP: DO IBND = 1,NUMBER_SPECTRAL_BANDS
 
   KAPPAW = 0._EB
   KFST4W = 0._EB
   SCAEFF = 0._EB
 
! Calculate fraction on ambient black body radiation
 
   IF (NUMBER_SPECTRAL_BANDS==1) THEN
      BBFA = 1._EB
   ELSE
      BBFA = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TMPA)
   ENDIF      

! Generate water absorption and scattering coefficients
 
   IF_DROPLETS_INCLUDED: IF (NLP>0 .AND. N_EVAP_INDICIES>0) THEN
      IF (NUMBER_SPECTRAL_BANDS==1) THEN
         BBF = 1._EB
      ELSE
         BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),RADTMP)
      ENDIF      
      DO K=1,KBAR
         DO J=1,JBAR
            ZLOOPM: DO I=1,IBAR
               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ZLOOPM
                  PC_LOOP: DO N = 1,N_PART
                     PC => PARTICLE_CLASS(N)
                     EVAP_INDEX = PC%EVAP_INDEX
                     IF (EVAP_INDEX==0) CYCLE PC_LOOP
                     IF (AVG_DROP_DEN(I,J,K,EVAP_INDEX)==0._EB) CYCLE PC_LOOP
                     NCSDROP = THFO*AVG_DROP_DEN(I,J,K,EVAP_INDEX)/ (PC%DENSITY*AVG_DROP_RAD(I,J,K,EVAP_INDEX))
                     ! Absorption and scattering efficiency
                     CALL INTERPOLATE1D(PC%KWR,PC%WQABS(:,IBND), 0.95_EB*AVG_DROP_RAD(I,J,K,EVAP_INDEX),QVAL) 
                     KAPPAW(I,J,K) = KAPPAW(I,J,K) + NCSDROP*QVAL
                     IF (PC%WATER) THEN
                        CALL INTERPOLATE1D(PC%KWR,PC%WQSCA(:,IBND),0.95_EB*AVG_DROP_RAD(I,J,K,EVAP_INDEX),QVAL)
                        SCAEFF(I,J,K) = SCAEFF(I,J,K) + NCSDROP*QVAL
                     ENDIF
                     KFST4W(I,J,K) = KFST4W(I,J,K)+ BBF*KAPPAW(I,J,K)*FOUR_SIGMA*AVG_DROP_TMP(I,J,K,EVAP_INDEX)**4
                  END DO PC_LOOP
            ENDDO ZLOOPM
         ENDDO
      ENDDO
      QR_W = 0._EB
   ENDIF IF_DROPLETS_INCLUDED
 
! Compute absorption coefficient KAPPA and source term KAPPA*4*SIGMA*TMP**4
 
   BBF = 1._EB
   KAPPA=KAPPA0
   DO K=1,KBAR
      DO J=1,JBAR
         ZLOOP: DO I=1,IBAR
            IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ZLOOP
            SUM_SPECIES: DO NS=1,N_SPECIES
               IF (SPECIES(NS)%MODE == MIXTURE_FRACTION_SPECIES) CYCLE SUM_SPECIES
               IF (.NOT. SPECIES(NS)%ABSORBING)                  CYCLE SUM_SPECIES
               TYY = NINT(SPECIES(NS)%NKAP_TEMP*(TMP(I,J,K)-300._EB)/2100._EB)
               TYY = MAX(0,MIN(SPECIES(NS)%NKAP_TEMP,TYY))
               ZZ = MIN(1._EB,MAX(0._EB,YY(I,J,K,NS)))
               KAPPA(I,J,K) = KAPPA(I,J,K) + ZZ2KAPPA(ZZ,SPECIES(NS)%MAXMASS,2*SPECIES(NS)%NKAP_MASS,TYY,IBND,NS)
            ENDDO SUM_SPECIES
            IF (MIXTURE_FRACTION) THEN
               TYY = NINT(SPECIES(I_FUEL)%NKAP_TEMP*(TMP(I,J,K)-300._EB)/2100._EB)
               TYY = MAX(0,MIN(SPECIES(I_FUEL)%NKAP_TEMP,TYY))
               IF (CO_PRODUCTION) THEN
                  Z_2 = YY(I,J,K,I_PROG_CO)
               ELSE
                  Z_2 = 0._EB
               ENDIF
               CALL GET_KAPPA(YY(I,J,K,I_FUEL),Z_2,YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),KAPPA_1,TYY,IBND)              
               KAPPA(I,J,K) = KAPPA(I,J,K) + KAPPA_1
            ENDIF
            IF (WIDE_BAND_MODEL)  BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TMP(I,J,K))
            KFST4(I,J,K) = BBF*KAPPA(I,J,K)*FOUR_SIGMA*TMP(I,J,K)**4
            IF (LES) THEN
               IF (BBF*RADIATIVE_FRACTION*Q(I,J,K) > KFST4(I,J,K)) THEN
                  KFST4(I,J,K) = BBF*RADIATIVE_FRACTION*Q(I,J,K)
                  KAPPA(I,J,K) = 0._EB
               ENDIF
            ENDIF
         ENDDO ZLOOP
      ENDDO
   ENDDO

! Calculate extinction coefficient
 
   EXTCOE = KAPPA + KAPPAW + SCAEFF*RSA_RAT

! Update intensity field
 
   INTENSITY_UPDATE: IF (UPDATE_INTENSITY) THEN
 
      IF (WIDE_BAND_MODEL) THEN
         UIIOLD = UIID(:,:,:,IBND)
      ELSE
         UIIOLD = UII
      ENDIF
      UII = 0._EB

! Compute boundary condition intensity emissivity*sigma*Tw**4/pi 
!     or emissivity*QRADOUT/pi for wall with internal radiation
 
      BBF = 1.0_EB
      DO IW = 1,NWC
         IF (WIDE_BAND_MODEL) BBF = BLACKBODY_FRACTION(WL_LOW(IBND),WL_HIGH(IBND),TMP_F(IW))
         IBC = IJKW(5,IW)
         SF  => SURFACE(IBC)
         IF (.NOT. SF%INTERNAL_RADIATION) THEN
            QRADOUT(IW) = E_WALL(IW)*SIGMA*TMP_F(IW)**4
         ENDIF
         OUTRAD_W(IW) = BBF*RPI*QRADOUT(IW)
      ENDDO
! Compute boundary condition term incoming radiation integral
 
      DO IW = 1,NWC
         IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) CYCLE
         IOR = IJKW(4,IW)
         INRAD_W(IW) = SUM(-W_AXI*DLN(IOR,:)* WALL(IW)%ILW(:,IBND),1, DLN(IOR,:)<0._EB)
      ENDDO
 
! If updating intensities first time, sweep ALL angles
 
      N_UPDATES = 1
      IF (RAD_CALL_COUNTER==1) N_UPDATES = ANGLE_INCREMENT

      UPDATE_LOOP: DO I_UIID = 1,N_UPDATES
 
! Update counters inside the radiation routine
 
         ANGLE_INC_COUNTER = MOD(ANGLE_INC_COUNTER,ANGLE_INCREMENT) + 1
         IF (WIDE_BAND_MODEL) THEN
            UIID(:,:,:,IBND) = 0._EB
         ELSE
            UIID(:,:,:,ANGLE_INC_COUNTER) = 0._EB
         ENDIF
 
! Set the bounds and increment for the angleloop. Step downdard because in cylindrical case the Nth angle 
! boundary condition comes from (N+1)th angle.
 
         NSTART    = NRA - ANGLE_INC_COUNTER + 1
         NEND      = 1
         NSTEP     = -ANGLE_INCREMENT
         IL(:,:,:) = BBFA*RPI_SIGMA*TMPA4

         ANGLE_LOOP: DO N = NSTART,NEND,NSTEP  ! Sweep through control angles
 
! Boundary conditions: Intensities leaving the boundaries.
 
            WALL_LOOP1: DO IW=1,NWC
               IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_LOOP1
               IOR = IJKW(4,IW)
               IF (DLN(IOR,N) < 0._EB) CYCLE WALL_LOOP1
               II  = IJKW(1,IW)
               JJ  = IJKW(2,IW)
               KK  = IJKW(3,IW)
               IF (.NOT.TWO_D .OR. ABS(IOR)/=2) THEN
                  SELECT CASE (BOUNDARY_TYPE(IW))
                     CASE (OPEN_BOUNDARY) 
                           IL(II,JJ,KK) = BBFA*RPI_SIGMA*TMPA4