fire.f90

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

MODULE FIRE
 
! Compute combustion 
 
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
USE COMP_FUNCTIONS, ONLY: SECOND
 
IMPLICIT NONE
PRIVATE
CHARACTER(255), PARAMETER :: fireid='$Id: fire.f90 712 2007-09-28 20:20:09Z drjfloyd $'
CHARACTER(255), PARAMETER :: firerev='$Revision: 712 $'
CHARACTER(255), PARAMETER :: firedate='$Date: 2007-09-28 16:20:09 -0400 (Fri, 28 Sep 2007) $'

TYPE(REACTION_TYPE), POINTER :: RN

PUBLIC COMBUSTION, GET_REV_fire
 
CONTAINS
 
SUBROUTINE COMBUSTION(NM)

INTEGER, INTENT(IN) :: NM
REAL(EB) :: TNOW

IF( EVACUATION_ONLY(NM)) RETURN

TNOW=SECOND()

CALL POINT_TO_MESH(NM)

IF (MIXTURE_FRACTION) THEN
   CALL COMBUSTION_MF
ELSE
   CALL COMBUSTION_FR
ENDIF

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

CONTAINS

SUBROUTINE COMBUSTION_MF
USE PHYSICAL_FUNCTIONS, ONLY : GET_MASS_FRACTION2,GET_MOLECULAR_WEIGHT2
REAL(EB) :: YFU0,A,ETRM,YCOMIN,YO2MIN,YFUMIN,YO2W,YO20,YCO0,Y_SUM_W,&
            DYF,DX_FDT,HFAC_F,DTT,& 
            Y_O2_MAX,TMP_MIN,Y_O2_CORR,Q_NEW,Q_OLD,&
            F_TO_CO,DELTAH_CO,DYCO,HFAC_CO,Z_2,RHOX, &
            X_FU,X_O,X_FU_0,X_O_0,X_FU_S,X_O_S,X_FU_N,X_O_N,X_O_MIN,X_FU_MIN,COTOO2, &
            Y_F_MAX,TMP_F_MIN,Y_F_CORR,Z2_MIN
INTEGER :: NODETS,N,I,J,K,II,IC,IW,IWA(-3:3)
!LOGICAL :: BURN
REAL(EB), POINTER, DIMENSION(:,:,:) :: YO2Z,QT,R_SUM_DILUENTS
!LOGICAL, POINTER, DIMENSION(:,:,:) :: IGNITE

PRODUCE_CO: IF (.NOT. CO_PRODUCTION) THEN !Combustion without CO formation and destruction
   YO2Z => WORK1
   Q =  0._EB
   Z_2 = 0._EB
   DO K=1,KBAR
      DO J=1,JBAR
         ILOOPA: DO I=1,IBAR
            IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ILOOPA
            CALL GET_MASS_FRACTION2(YY(I,J,K,I_FUEL),Z_2,YY(I,J,K,I_PROG_F),O2_INDEX,Y_SUM(I,J,K),YO2Z(I,J,K))   
         ENDDO ILOOPA
      ENDDO
   ENDDO
   YO2MIN = 1.E-10_EB
   YFUMIN = 1.E-10_EB
   !Loop and do F -> CO
   RN => REACTION(1)
   HFAC_F   = RN%HEAT_OF_COMBUSTION/DT
   SUPPRESSIONIF: IF (SUPPRESSION) THEN !TMP and Y_O2 dependent combustion
      DO K=1,KBAR
         DO J=1,JBAR
            ILOOPB: DO I=1,IBAR
               IC = CELL_INDEX(I,J,K)
               IF (SOLID(IC)) CYCLE ILOOPB
               IWA = WALL_INDEX(IC,:)
               YO20  = YO2Z(I,J,K)
               YFU0  = MAX(0._EB,MIN(1._EB,YY(I,J,K,I_FUEL)))*RN%Y_F_INLET
               IF (YFU0<=YFUMIN .OR. YO20<=YO2MIN) CYCLE ILOOPB

               !Get min O2
               Y_O2_MAX = YO20
               Y_F_MAX = YFU0/RN%Y_F_INLET
               TMP_MIN = TMP(I,J,K)
               TMP_F_MIN = TMP(I,J,K)
               !Check neighboring cells for fuel and oxygen
               !X direction
               IF (IWA(-1)==0) THEN
                  IF (YO2Z(I-1,J,K)>Y_O2_MAX) THEN
                     Y_O2_MAX = YO2Z(I-1,J,K)
                     TMP_MIN = TMP(I-1,J,K)
                  ENDIF
                  IF (YY(I-1,J,K,I_FUEL)>Y_F_MAX) THEN
                     Y_F_MAX = YY(I-1,J,K,I_FUEL)
                     TMP_F_MIN = TMP(I-1,J,K)
                  ENDIF
               ELSE
                  IW = IWA(-1)
                  IF (SURFACE(IJKW(5,IW))%SPECIES_BC_INDEX/=NO_MASS_FLUX) THEN
                     Y_SUM_W = 0._EB           
                     DO N=1,N_SPECIES
                        IF (SPECIES(N)%MODE==GAS_SPECIES) THEN
                           Y_SUM_W = Y_SUM_W + YY_W(IW,N)
                        ENDIF
                     ENDDO
                     CALL GET_MASS_FRACTION2(YY_W(IW,I_FUEL),Z_2,YY_W(IW,I_PROG_F),O2_INDEX,Y_SUM_W,YO2W)   
                     IF (YO2W>Y_O2_MAX) THEN
                        Y_O2_MAX = YO2W
                        TMP_MIN = TMP_F(IW)
                     ENDIF
                     IF (MAX(0._EB,MIN(1._EB,YY_W(IW,I_FUEL)))>Y_F_MAX) THEN
                        Y_F_MAX = YY_W(IW,I_FUEL)
                        TMP_F_MIN = TMP_F(IW)
                     ENDIF
                  ENDIF
               ENDIF
               IF (IWA(1)==0) THEN
                  IF (YO2Z(I+1,J,K)>Y_O2_MAX) THEN
                     Y_O2_MAX = YO2Z(I+1,J,K)
                     TMP_MIN = TMP(I+1,J,K)
                  ENDIF
                  IF (YY(I+1,J,K,I_FUEL)>Y_F_MAX) THEN
                     Y_F_MAX = YY(I+1,J,K,I_FUEL)
                     TMP_F_MIN = TMP(I+1,J,K)
                  ENDIF
               ELSE
                  IW = IWA(1)
                  IF (SURFACE(IJKW(5,IW))%SPECIES_BC_INDEX/=NO_MASS_FLUX) THEN    
                     Y_SUM_W = 0._EB           
                     DO N=1,N_SPECIES
                        IF (SPECIES(N)%MODE==GAS_SPECIES) THEN
                           Y_SUM_W = Y_SUM_W + YY_W(IW,N)
                        ENDIF
                     ENDDO
                     CALL GET_MASS_FRACTION2(YY_W(IW,I_FUEL),Z_2,YY_W(IW,I_PROG_F),O2_INDEX,Y_SUM_W,YO2W)   
                     IF (YO2W>Y_O2_MAX) THEN
                        Y_O2_MAX = YO2W
                        TMP_MIN = TMP_F(IW)
                     ENDIF
                     IF (MAX(0._EB,MIN(1._EB,YY_W(IW,I_FUEL)))>Y_F_MAX) THEN
                        Y_F_MAX = YY_W(IW,I_FUEL)
                        TMP_F_MIN = TMP_F(IW)
                     ENDIF
                  ENDIF
               ENDIF
               !Y direction            
               IF (IWA(-2)==0) THEN
                  IF (YO2Z(I,J-1,K)>Y_O2_MAX) THEN
                     Y_O2_MAX = YO2Z(I,J-1,K)
                     TMP_MIN = TMP(I,J-1,K)
                  ENDIF
                  IF (YY(I,J-1,K,I_FUEL)>Y_F_MAX) THEN
                     Y_F_MAX = YY(I,J-1,K,I_FUEL)
                     TMP_F_MIN = TMP(I,J-1,K)
                  ENDIF
               ELSE
                  IW = IWA(-2)
                  IF (SURFACE(IJKW(5,IW))%SPECIES_BC_INDEX/=NO_MASS_FLUX) THEN               
                     Y_SUM_W = 0._EB           
                     DO N=1,N_SPECIES
                        IF (SPECIES(N)%MODE==GAS_SPECIES) THEN
                           Y_SUM_W = Y_SUM_W + YY_W(IW,N)
                        ENDIF
                     ENDDO
                     CALL GET_MASS_FRACTION2(YY_W(IW,I_FUEL),Z_2,YY_W(IW,I_PROG_F),O2_INDEX,Y_SUM_W,YO2W)   
                     IF (YO2W>Y_O2_MAX) THEN
                        Y_O2_MAX = YO2W
                        TMP_MIN = TMP_F(IW)
                     ENDIF
                     IF (MAX(0._EB,MIN(1._EB,YY_W(IW,I_FUEL)))>Y_F_MAX) THEN
                        Y_F_MAX = YY_W(IW,I_FUEL)
                        TMP_F_MIN = TMP_F(IW)
                     ENDIF
                  ENDIF
               ENDIF
               IF (IWA(2)==0) THEN
                  IF (YO2Z(I,J+1,K)>Y_O2_MAX) THEN
                     Y_O2_MAX = YO2Z(I,J+1,K)
                     TMP_MIN = TMP(I,J+1,K)
                  ENDIF
                  IF (YY(I,J+1,K,I_FUEL)>Y_F_MAX) THEN
                     Y_F_MAX = YY(I,J+1,K,I_FUEL)
                     TMP_F_MIN = TMP(I,J+1,K)
                  ENDIF
               ELSE
                  IW = IWA(2)
                  IF (SURFACE(IJKW(5,IW))%SPECIES_BC_INDEX/=NO_MASS_FLUX) THEN               
                     Y_SUM_W = 0._EB           
                     DO N=1,N_SPECIES
                        IF (SPECIES(N)%MODE==GAS_SPECIES) THEN
                           Y_SUM_W = Y_SUM_W + YY_W(IW,N)
                        ENDIF
                     ENDDO
                     CALL GET_MASS_FRACTION2(YY_W(IW,I_FUEL),Z_2,YY_W(IW,I_PROG_F),O2_INDEX,Y_SUM_W,YO2W)   
                     IF (YO2W>Y_O2_MAX) THEN
                        Y_O2_MAX = YO2W
                        TMP_MIN = TMP_F(IW)
                     ENDIF
                     IF (MAX(0._EB,MIN(1._EB,YY_W(IW,I_FUEL)))>Y_F_MAX) THEN
                        Y_F_MAX = YY_W(IW,I_FUEL)
                        TMP_F_MIN = TMP_F(IW)
                     ENDIF
                  ENDIF
               ENDIF                        
               !Z direction
               IF (IWA(-3)==0) THEN
                  IF (YO2Z(I,J,K-1)>Y_O2_MAX) THEN
                     Y_O2_MAX = YO2Z(I,J,K-1)
                     TMP_MIN = TMP(I,J,K-1)
                  ENDIF
                  IF (YY(I,J,K-1,I_FUEL)>Y_F_MAX) THEN
                     Y_F_MAX = YY(I,J,K-1,I_FUEL)
                     TMP_F_MIN = TMP(I,J,K-1)
                  ENDIF
               ELSE
                  IW = IWA(-3)
                  IF (SURFACE(IJKW(5,IW))%SPECIES_BC_INDEX/=NO_MASS_FLUX) THEN               
                     Y_SUM_W = 0._EB           
                     DO N=1,N_SPECIES
                        IF (SPECIES(N)%MODE==GAS_SPECIES) THEN
                           Y_SUM_W = Y_SUM_W + YY_W(IW,N)
                        ENDIF
                     ENDDO
                     CALL GET_MASS_FRACTION2(YY_W(IW,I_FUEL),Z_2,YY_W(IW,I_PROG_F),O2_INDEX,Y_SUM_W,YO2W)   
                     IF (YO2W>Y_O2_MAX) THEN
                        Y_O2_MAX = YO2W
                        TMP_MIN = TMP_F(IW)
                     ENDIF
                     IF (MAX(0._EB,MIN(1._EB,YY_W(IW,I_FUEL)))>Y_F_MAX) THEN
                        Y_F_MAX = YY_W(IW,I_FUEL)
                        TMP_F_MIN = TMP_F(IW)
                     ENDIF
                  ENDIF
               ENDIF
               IF (IWA(3)==0) THEN
                  IF (YO2Z(I,J,K+1)>Y_O2_MAX) THEN
                     Y_O2_MAX = YO2Z(I,J,K+1)
                     TMP_MIN = TMP(I,J,K+1)
                  ENDIF
                  IF (YY(I,J,K+1,I_FUEL)>Y_F_MAX) THEN
                     Y_F_MAX = YY(I,J,K+1,I_FUEL)
                     TMP_F_MIN = TMP(I,J,K+1)
                  ENDIF
               ELSE
                  IW = IWA(3)
                  IF (SURFACE(IJKW(5,IW))%SPECIES_BC_INDEX/=NO_MASS_FLUX) THEN
                     Y_SUM_W = 0._EB           
                     DO N=1,N_SPECIES
                        IF (SPECIES(N)%MODE==GAS_SPECIES) THEN
                           Y_SUM_W = Y_SUM_W + YY_W(IW,N)
                        ENDIF
                     ENDDO
                     CALL GET_MASS_FRACTION2(YY_W(IW,I_FUEL),Z_2,YY_W(IW,I_PROG_F),O2_INDEX,Y_SUM_W,YO2W)   
                     IF (YO2W>Y_O2_MAX) THEN
                        Y_O2_MAX = YO2W
                        TMP_MIN = TMP_F(IW)
                     ENDIF
                     IF (MAX(0._EB,MIN(1._EB,YY_W(IW,I_FUEL)))>Y_F_MAX) THEN
                        Y_F_MAX = YY_W(IW,I_FUEL)
                        TMP_F_MIN = TMP_F(IW)
                     ENDIF
                  ENDIF
               ENDIF
               Y_O2_CORR = RN%Y_O2_LL*(RN%CRIT_FLAME_TMP-TMP_MIN)/(RN%CRIT_FLAME_TMP-TMPA)
               Y_F_CORR  = RN%Y_F_LFL*(RN%CRIT_FLAME_TMP-TMP_F_MIN)/(RN%CRIT_FLAME_TMP-TMPA)
               IF (Y_O2_MAX < Y_O2_CORR .OR. Y_F_MAX*RN%Y_F_INLET < Y_F_CORR) CYCLE ILOOPB
               DYF = MIN(YFU0,YO20/RN%O2_F_RATIO)
               Q_NEW = MIN(Q_UPPER,DYF*RHO(I,J,K)*HFAC_F)
               DYF = Q_NEW /(RHO(I,J,K)*HFAC_F*RN%Y_F_INLET)
               Q(I,J,K) = Q_NEW
               YY(I,J,K,I_FUEL) = YY(I,J,K,I_FUEL) - DYF
               YY(I,J,K,I_PROG_F) = YY(I,J,K,I_PROG_F) + DYF
            ENDDO ILOOPB
         ENDDO
      ENDDO
      
   ELSE SUPPRESSIONIF !No suppression
      DO K=1,KBAR
         DO J=1,JBAR
            ILOOPC: DO I=1,IBAR
               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ILOOPC
               YO20  = YO2Z(I,J,K)
               YFU0  = MAX(0._EB,MIN(1._EB,YY(I,J,K,I_FUEL)))*RN%Y_F_INLET
               IF (YFU0<=YFUMIN .OR. YO20<=YO2MIN) CYCLE ILOOPC
               DYF = MIN(YFU0,YO20/RN%O2_F_RATIO)
               Q_NEW = MIN(Q_UPPER,DYF*RHO(I,J,K)*HFAC_F)
               DYF = Q_NEW /(RHO(I,J,K)*HFAC_F*RN%Y_F_INLET)
               Q(I,J,K) = Q_NEW
               YY(I,J,K,I_FUEL) = YY(I,J,K,I_FUEL) - DYF
               YY(I,J,K,I_PROG_F) = YY(I,J,K,I_PROG_F) + DYF
            ENDDO ILOOPC
         ENDDO
      ENDDO
   ENDIF SUPPRESSIONIF
       
ELSE PRODUCE_CO !Combustion with suppression and CO production
   YO2Z   => WORK1
!   IGNITE => LOGICAL_WORK
!   IGNITE = .FALSE.
   QT     => WORK3
   DO K=1,KBAR
      DO J=1,JBAR
         ILOOPX: DO I=1,IBAR
            IC = CELL_INDEX(I,J,K)
            IF (SOLID(IC)) CYCLE ILOOPX
!            IF (Q(I,J,K) > 0._EB) IGNITE(I,J,K) = .TRUE.
!            IF (MAXVAL(WALL_INDEX(IC,:)) > 0) THEN
!               DO N = -3,3
!                  IF (WALL_INDEX(IC,N) > 0) THEN
!                     IF (SURFACE(IJKW(5,WALL_INDEX(IC,N)))%HRRPUA>0._EB)  IGNITE(I,J,K) = .TRUE.
!                  ENDIF
!               END DO
!            ENDIF
            CALL GET_MASS_FRACTION2(YY(I,J,K,I_FUEL),YY(I,J,K,I_PROG_CO),YY(I,J,K,I_PROG_F),O2_INDEX,Y_SUM(I,J,K),YO2Z(I,J,K))
         ENDDO ILOOPX
      ENDDO
   ENDDO
   QT = 0._EB
   Q =  0._EB
   RN => REACTION(1)
   HFAC_F   = RN%HEAT_OF_COMBUSTION/DT
   F_TO_CO = RN%MW_FUEL/(RN%NU_CO*MW_CO)  
   YO2MIN = 1.E-10_EB
   YFUMIN = 1.E-10_EB
   !Loop and do F -> CO
   DO K=1,KBAR
      DO J=1,JBAR
         ILOOPY: DO I=1,IBAR
            IC = CELL_INDEX(I,J,K)
            IF (SOLID(IC)) CYCLE ILOOPY
            IWA = WALL_INDEX(IC,:)
            YO20  = YO2Z(I,J,K)
            YFU0  = MIN(1._EB,MAX(0._EB,YY(I,J,K,I_FUEL)))*RN%Y_F_INLET
            IF (YFU0<=YFUMIN .OR. YO20<=YO2MIN) CYCLE ILOOPY
            !Get min O2
            Y_O2_MAX = YO20
            Y_F_MAX = YFU0/RN%Y_F_INLET
            TMP_MIN = TMP(I,J,K)
            TMP_F_MIN = TMP(I,J,K)
            !Check neighboring cells for fuel and oxygen
            !X direction
            IF (IWA(-1)==0) THEN
               IF (YO2Z(I-1,J,K)>Y_O2_MAX) THEN
                  Y_O2_MAX = YO2Z(I-1,J,K)
                  TMP_MIN = TMP(I-1,J,K)
               ENDIF
               IF (YY(I-1,J,K,I_FUEL)>Y_F_MAX) THEN
                  Y_F_MAX = YY(I-1,J,K,I_FUEL)
                  TMP_F_MIN = TMP(I-1,J,K)
               ENDIF
            ELSE
               IW = IWA(-1)
               IF (SURFACE(IJKW(5,IW))%SPECIES_BC_INDEX/=NO_MASS_FLUX) THEN
                  Y_SUM_W = 0._EB           
                  DO N=1,N_SPECIES
                     IF (SPECIES(N)%MODE==GAS_SPECIES) THEN
                        Y_SUM_W = Y_SUM_W + YY_W(IW,N)
                     ENDIF
                  ENDDO
                  CALL GET_MASS_FRACTION2(YY_W(IW,I_FUEL),YY_W(IW,I_PROG_CO),YY_W(IW,I_PROG_F),O2_INDEX,Y_SUM_W,YO2W)   
                  IF (YO2W>Y_O2_MAX) THEN
                     Y_O2_MAX = YO2W
                     TMP_MIN = TMP_F(IW)
                  ENDIF
                  IF (MAX(0._EB,MIN(1._EB,YY_W(IW,I_FUEL)))>Y_F_MAX) THEN
                     Y_F_MAX = YY_W(IW,I_FUEL)
                     TMP_F_MIN = TMP_F(IW)
                  ENDIF
               ENDIF
            ENDIF
            IF (IWA(1)==0) THEN
               IF (YO2Z(I+1,J,K)>Y_O2_MAX) THEN
                  Y_O2_MAX = YO2Z(I+1,J,K)
                  TMP_MIN = TMP(I+1,J,K)
               ENDIF
               IF (YY(I+1,J,K,I_FUEL)>Y_F_MAX) THEN
                  Y_F_MAX = YY(I+1,J,K,I_FUEL)
                  TMP_F_MIN = TMP(I+1,J,K)
               ENDIF
            ELSE
               IW = IWA(1)
               IF (SURFACE(IJKW(5,IW))%SPECIES_BC_INDEX/=NO_MASS_FLUX) THEN               
                  Y_SUM_W = 0._EB