read.f90

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

   IF (REACTION(N)%MODE==MIXTURE_FRACTION_REACTION) THEN
      N_MIX  = N_MIX  + 1
      N_SPECIES = N_SPECIES + 1
      IF (N_MIX==1) THEN
         I_FUEL = N_SPECIES
      ELSEIF (N_MIX== 2) THEN
         IF(CO_PRODUCTION) THEN 
            I_PROG_CO = N_SPECIES
         ELSE
            I_PROG_F = N_SPECIES
         ENDIF
      ELSEIF (N_MIX==3) THEN
         I_PROG_F = N_SPECIES
      ENDIF
      YYMIN(N_MIX)  = 0.0_EB
      YYMIN(I_FUEL) = 0.0_EB
   ENDIF
ENDDO
 
IF (WATER_EVAPORATION .AND. I_WATER==0) THEN
   N_SPECIES  = N_SPECIES + 1
   I_WATER = N_SPECIES
ENDIF
 
! Allocate species-related arrays
 
ALLOCATE(SPECIES(0:N_SPECIES),STAT=IZERO)
CALL ChkMemErr('READ','SPECIES',IZERO)
SPECIES(:)%YY0 = 0._EB
SPECIES(:)%MW  = 0._EB
SPECIES(0)%MW  = MW_BACKGROUND
SPECIES(0)%NAME= BACKGROUND_SPECIES
EPSK           = 0._EB
SIG            = 0._EB
SPECIES(:)%ABSORBING=.FALSE.
SPECIES(:)%MODE=GAS_SPECIES
 
N_MIX=0
NN=N_SPEC_EXTRA
DO N=1,N_REACTIONS
   IF (REACTION(N)%MODE==MIXTURE_FRACTION_REACTION) THEN
      NN    = NN + 1
      N_MIX = N_MIX + 1
      SPECIES(NN)%MODE = MIXTURE_FRACTION_SPECIES
      SPECIES(NN)%REAC_INDEX    = N
      SPECIES(NN)%ABSORBING=.TRUE.
      IF (N_REACTIONS > 1) WRITE(SPECIES(NN)%NAME,'(A,I1.1)') 'MIXTURE_FRACTION_',N_MIX
      IF (N_REACTIONS ==1) WRITE(SPECIES(NN)%NAME,'(A)'     ) 'MIXTURE_FRACTION'
   ENDIF
ENDDO
 
! Initialize initial WATER VAPOR mass fraction, if necessary
 
IF (I_WATER>0) THEN
   XVAP  = MIN(1._EB,EXP(2259.E3_EB*MW_H2O/R0*(1._EB/373.15_EB-1._EB/ MIN(TMPA,373.15_EB))))
   SPECIES(I_WATER)%MODE = GAS_SPECIES
   SPECIES(I_WATER)%NAME = 'WATER VAPOR'
   SPECIES(I_WATER)%YY0 =HUMIDITY*0.01_EB*XVAP/(MW_AIR/MW_H2O+(1._EB-MW_AIR/MW_H2O)*XVAP)
ENDIF
IF (I_SOOT>0) SPECIES(I_SOOT)%MODE=AEROSOL_SPECIES
 
! Read SPEC lines from input file
 
REWIND(LU_INPUT)
SPEC_LOOP: DO N=1,N_SPEC_READ
 
   CONDUCTIVITY    = -1._EB
   DIFFUSIVITY     = -1._EB
   EPSILONKLJ      = 0._EB
   MASS_FRACTION_0 = -1._EB
   MW              = 0._EB
   SIGMALJ         = 0._EB
   VISCOSITY       = -1._EB
   ABSORBING       = .FALSE.
   ID              = 'null'
 
   CALL CHECKREAD('SPEC',LU_INPUT,IOS)
   READ(LU_INPUT,NML=SPEC)
   IF (MIXTURE_FRACTION) THEN
      SELECT CASE(ID)
         CASE('MIXTURE_FRACTION_1')
            NN = I_FUEL
            ABSORBING  = .TRUE.
         CASE('MIXTURE_FRACTION_2')
            ABSORBING  = .TRUE.
            IF (CO_PRODUCTION) THEN
               NN = I_PROG_CO
            ELSE
               NN = I_PROG_F
            ENDIF
         CASE('MIXTURE_FRACTION_3')
            IF (CO_PRODUCTION) THEN
               NN = I_PROG_CO
               ABSORBING  = .TRUE.
            ELSE
               NN = N
            ENDIF
         CASE DEFAULT
            NN = N
      END SELECT
   ELSE
      NN = N
   ENDIF
   SS => SPECIES(NN) 
   MU_USER(NN) = VISCOSITY
   K_USER(NN)  = CONDUCTIVITY
   D_USER(NN)  = DIFFUSIVITY
   EPSK(NN)    = EPSILONKLJ
   SIG(NN)     = SIGMALJ
   SS%MW        = MW
   SS%ABSORBING = ABSORBING
   IF (MASS_FRACTION_0 < 0._EB) THEN
      IF (N/=I_WATER)  SS%YY0 = 0._EB
   ELSE
      SS%YY0 = MASS_FRACTION_0
   ENDIF
   IF (ID/='null') SS%NAME = ID
 
ENDDO SPEC_LOOP

END SUBROUTINE READ_SPEC
 
 
SUBROUTINE READ_REAC
 
CHARACTER(30) :: FUEL,OXIDIZER
LOGICAL :: IDEAL
INTEGER :: NN,N_REAC_READ
REAL(EB) :: Y_O2_INFTY,Y_F_INLET, &
            H2_YIELD,SOOT_YIELD,CO_YIELD, Y_F_LFL,X_O2_LL,EPUMO2,BOF, SOOT_H_FRACTION,&
            CRITICAL_FLAME_TEMPERATURE,HEAT_OF_COMBUSTION,NU(MAX_SPECIES),E,N_S(MAX_SPECIES),C,H,N,O,OTHER,MW_OTHER, &
            FUEL_HEAT_OF_FORMATION,MAXIMUM_VISIBILITY 
NAMELIST /REAC/ E,BOF,HEAT_OF_COMBUSTION,FYI,FUEL,OXIDIZER,EPUMO2,ID, N_S,&
                Y_O2_INFTY,Y_F_INLET,HRRPUA_SHEET, &
                H2_YIELD,SOOT_YIELD,CO_YIELD,Y_F_LFL,X_O2_LL,CRITICAL_FLAME_TEMPERATURE,NU,SOOT_H_FRACTION, &
                C,H,N,O,OTHER,MW_OTHER,IDEAL,MASS_EXTINCTION_COEFFICIENT,VISIBILITY_FACTOR,MAXIMUM_VISIBILITY
 
N_REACTIONS = 0
REWIND(LU_INPUT)
 
COUNT_REAC_LOOP: DO
   ID   = 'null'
   FUEL = 'null'      
   CALL CHECKREAD('REAC',LU_INPUT,IOS)
   IF (IOS==1) EXIT COUNT_REAC_LOOP
   READ(LU_INPUT,REAC,ERR=434,IOSTAT=IOS)
   N_REACTIONS = N_REACTIONS + 1
   IF (FUEL=='null') THEN
      MIXTURE_FRACTION = .TRUE.
   ENDIF
   IF (FUEL/='null' .AND. MIXTURE_FRACTION) THEN
      WRITE(MESSAGE,'(A)') 'ERROR: cannot use both finite rate REAC and mixture fraction REAC'
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   IF (MIXTURE_FRACTION .AND. N_REACTIONS > 1) THEN
      WRITE(MESSAGE,'(A)') 'ERROR: can not have more than one reaction when using mixture fraction'
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   434 IF (IOS>0) THEN
      WRITE(MESSAGE,'(A,I3)') 'ERROR: Problem with REAC ',N_REACTIONS+1
      CALL SHUTDOWN(MESSAGE)
   ENDIF
ENDDO COUNT_REAC_LOOP
 
IF (FUEL_EVAPORATION) MIXTURE_FRACTION = .TRUE.
IF (MIXTURE_FRACTION .AND. N_REACTIONS==0) N_REACTIONS = 1

IF (N_REACTIONS==0) RETURN

IF (MIXTURE_FRACTION) THEN
   N_REAC_READ = 1
ELSE
   N_REAC_READ = N_REACTIONS
ENDIF

IF (MIXTURE_FRACTION .AND.      CO_PRODUCTION) N_REACTIONS = 3
IF (MIXTURE_FRACTION .AND. .NOT.CO_PRODUCTION) N_REACTIONS = 2
ALLOCATE(REACTION(N_REACTIONS),STAT=IZERO)
CALL ChkMemErr('READ','REACTION',IZERO)

! Read the input file looking for REAC lines
 
REWIND(LU_INPUT)
READ_REACTION_LOOP: DO NN=1,N_REAC_READ
   RN => REACTION(NN)
   CALL CHECKREAD('REAC',LU_INPUT,IOS) 
   CALL SET_REAC_DEFAULTS
   IF (IOS==0) READ(LU_INPUT,REAC)
   RN%BOF                = BOF
   RN%E                  = E*1000._EB
   RN%EPUMO2             = EPUMO2*1000._EB
   RN%FUEL               = FUEL        
   RN%HEAT_OF_COMBUSTION = HEAT_OF_COMBUSTION*1000._EB
   IF (FUEL=='null') THEN
      RN%MODE = MIXTURE_FRACTION_REACTION
   ELSE
      RN%MODE = FINITE_RATE_REACTION
   ENDIF
   RN%N(:)               = N_S(:)
   RN%NAME               = ID
   RN%NU(:)              = NU(:)
   RN%OXIDIZER           = OXIDIZER
ENDDO READ_REACTION_LOOP

SET_MIXTURE_FRACTION: IF (MIXTURE_FRACTION) THEN
   !Set reaction variable constants
   REACTION%MW_OTHER        = MW_OTHER
   REACTION%MW_FUEL         = C * MW_C + H * MW_H + O * MW_O + N * MW_N + OTHER * MW_OTHER
   IF (H==0._EB) THEN
      REACTION%SOOT_H_FRACTION = 0._EB
   ELSE
      REACTION%SOOT_H_FRACTION = SOOT_H_FRACTION
   ENDIF
   REACTION%Y_O2_INFTY      = Y_O2_INFTY
   REACTION%Y_O2_LL         = X_O2_LL*MW_O2/ (X_O2_LL*MW_O2+(1._EB-X_O2_LL)*MW_N2)
   REACTION%Y_F_LFL         = Y_F_LFL
   REACTION%CRIT_FLAME_TMP  = CRITICAL_FLAME_TEMPERATURE + TMPM
   REACTION%Y_F_INLET       = Y_F_INLET
   REACTION%MODE            = MIXTURE_FRACTION_REACTION
   REACTION%NAME            = ID
   MW_SOOT                  = MW_C * (1._EB - SOOT_H_FRACTION) + MW_H * SOOT_H_FRACTION   
   
   SET_THREE_PARAMETER: IF (CO_PRODUCTION) THEN !Set reaction variables for three parameter mixture fraction 

      !Set reaction variables for complete reaction
      RN => REACTION(2) 
      RN%IDEAL      = IDEAL
      IF (RN%IDEAL) THEN
         !Compute fuel heat of formation
         RN%NU_O2 = C + 0.5_EB * H - 0.5_EB * O
         IF (HEAT_OF_COMBUSTION < 0._EB) HEAT_OF_COMBUSTION = EPUMO2*(MW_O2*RN%NU_O2)/RN%MW_FUEL
         FUEL_HEAT_OF_FORMATION =  HEAT_OF_COMBUSTION * RN%MW_FUEL - &
                                   (C * CO2_HEAT_OF_FORMATION + 0.5_EB * H * H2O_HEAT_OF_FORMATION)
      ENDIF
      RN%CO_YIELD   = CO_YIELD
      RN%SOOT_YIELD = SOOT_YIELD
      RN%H2_YIELD   = H2_YIELD
      RN%NU_H2      = H2_YIELD * RN%MW_FUEL / MW_H2
      RN%NU_SOOT    = SOOT_YIELD * RN%MW_FUEL / MW_SOOT
      RN%NU_CO      = CO_YIELD * RN%MW_FUEL / MW_CO
      RN%NU_CO2     = C - RN%NU_CO - RN%NU_SOOT * (1._EB - SOOT_H_FRACTION)
      IF (RN%NU_CO2 < 0._EB) THEN
         WRITE(MESSAGE,'(A)') 'Values for SOOT_YIELD, CO_YIELD, and SOOT_H_FRACTION result in negative CO2 yield'
         CALL SHUTDOWN(MESSAGE)
      ENDIF
      RN%NU_H2O     = 0.5_EB * H - RN%NU_H2 - 0.5_EB * RN%NU_SOOT * SOOT_H_FRACTION
      IF (RN%NU_H2O < 0._EB) THEN
         WRITE(MESSAGE,'(A)') 'Values for SOOT_YIELD, H2_YIELD, and SOOT_H_FRACTION result in negative H2O yield'
         CALL SHUTDOWN(MESSAGE)
      ENDIF
      RN%NU_O2      = RN%NU_CO2 + 0.5_EB * (RN%NU_CO + RN%NU_H2O - O)
      RN%NU_N2      = 0.5_EB * N
      RN%NU_OTHER   = OTHER
      RN%BOF        = 2.53E12_EB
      RN%E          = 199547._EB*1000._EB
      RN%O2_F_RATIO    = (MW_O2*RN%NU_O2)/RN%MW_FUEL
      IF (.NOT. RN%IDEAL) THEN
         IF(HEAT_OF_COMBUSTION < 0._EB) HEAT_OF_COMBUSTION = EPUMO2*RN%O2_F_RATIO
         RN%HEAT_OF_COMBUSTION = HEAT_OF_COMBUSTION*1000._EB
      ELSE
         ! Correct heat of combustion for minor products of combustion
         RN%HEAT_OF_COMBUSTION = (FUEL_HEAT_OF_FORMATION + RN%NU_CO2 * CO2_HEAT_OF_FORMATION + &
                                                         RN%NU_CO  * CO_HEAT_OF_FORMATION + &
                                                         RN%NU_H2O * H2O_HEAT_OF_FORMATION) * 1000._EB /RN%MW_FUEL
      ENDIF
      
      ! Set reaction variables for incomplete reaction

      RN => REACTION(1)
      RN%IDEAL              = .TRUE.
      RN%CO_YIELD           = 0._EB
      RN%H2_YIELD           = H2_YIELD
      RN%SOOT_YIELD         = SOOT_YIELD
      RN%NU_H2              = H2_YIELD * RN%MW_FUEL / MW_H2
      RN%NU_SOOT            = SOOT_YIELD * RN%MW_FUEL / MW_SOOT
      RN%NU_CO              = C - RN%NU_SOOT * (1._EB - SOOT_H_FRACTION)
      RN%NU_CO2             = 0._EB
      IF (RN%NU_CO2 < 0._EB) THEN
         WRITE(MESSAGE,'(A)') 'Values for SOOT_YIELD, CO_YIELD, and SOOT_H_FRACTION result in negative CO2 yield'
         CALL SHUTDOWN(MESSAGE)
      ENDIF
      RN%NU_H2O             = 0.5_EB * H - RN%NU_H2 - 0.5_EB * RN%NU_SOOT * SOOT_H_FRACTION
      IF (RN%NU_H2O < 0._EB) THEN
         WRITE(MESSAGE,'(A)') 'Values for SOOT_YIELD, H2_YIELD, and SOOT_H_FRACTION result in negative H2O yield'
         CALL SHUTDOWN(MESSAGE)
      ENDIF
      RN%NU_O2              = RN%NU_CO2 + 0.5_EB * (RN%NU_CO + RN%NU_H2O - O)
      RN%NU_N2              = 0.5_EB * N
      RN%NU_OTHER           = OTHER
      RN%O2_F_RATIO    = (MW_O2*RN%NU_O2)/RN%MW_FUEL      
      RN%HEAT_OF_COMBUSTION = REACTION(2)%HEAT_OF_COMBUSTION - 1.E3_EB * REACTION(2)%NU_CO2 /RN%MW_FUEL * &
                             (CO2_HEAT_OF_FORMATION - CO_HEAT_OF_FORMATION)
      NN = 3
      
   ELSE SET_THREE_PARAMETER !Set reaction variables for two parameter mixture fraction 
   
      !Set reaction variables for complete reaction
      RN => REACTION(1) 
      RN%IDEAL      = IDEAL
      IF (RN%IDEAL) THEN
         !Compute fuel heat of formation
         RN%NU_O2 = C + 0.5_EB * H - 0.5_EB * O
         IF (HEAT_OF_COMBUSTION < 0._EB) HEAT_OF_COMBUSTION = EPUMO2*(MW_O2*RN%NU_O2)/RN%MW_FUEL
         FUEL_HEAT_OF_FORMATION =  HEAT_OF_COMBUSTION * RN%MW_FUEL - &
                                   (C * CO2_HEAT_OF_FORMATION + 0.5_EB * H * H2O_HEAT_OF_FORMATION)
      ENDIF
      RN%CO_YIELD   = CO_YIELD
      RN%SOOT_YIELD = SOOT_YIELD
      RN%H2_YIELD   = H2_YIELD
      RN%NU_H2      = H2_YIELD * RN%MW_FUEL / MW_H2
      RN%NU_SOOT    = SOOT_YIELD * RN%MW_FUEL / MW_SOOT
      RN%NU_CO      = CO_YIELD * RN%MW_FUEL / MW_CO
      RN%NU_CO2     = C - RN%NU_CO - RN%NU_SOOT * (1._EB - SOOT_H_FRACTION)
      RN%NU_H2O     = 0.5_EB * H - RN%NU_H2 - 0.5_EB * RN%NU_SOOT * SOOT_H_FRACTION
      RN%NU_O2      = RN%NU_CO2 + 0.5_EB * (RN%NU_CO + RN%NU_H2O - O)
      RN%NU_N2      = 0.5_EB * N
      RN%NU_OTHER   = OTHER
      RN%O2_F_RATIO    = (MW_O2*RN%NU_O2)/RN%MW_FUEL      
      IF (.NOT. RN%IDEAL) THEN
         IF(HEAT_OF_COMBUSTION < 0._EB) HEAT_OF_COMBUSTION = EPUMO2*RN%O2_F_RATIO
         RN%HEAT_OF_COMBUSTION = HEAT_OF_COMBUSTION*1000._EB
      ELSE
         !Correct heat of combustion for minor products of combustion
         RN%HEAT_OF_COMBUSTION = (FUEL_HEAT_OF_FORMATION + RN%NU_CO2 * CO2_HEAT_OF_FORMATION + &
                                                         RN%NU_CO  * CO_HEAT_OF_FORMATION + &
                                                         RN%NU_H2O * H2O_HEAT_OF_FORMATION) * 1000._EB /RN%MW_FUEL
      ENDIF
      NN = 2
   ENDIF SET_THREE_PARAMETER
   
   !Set reaction variables for extinction reaction   

   RN => REACTION(NN)
   RN%CO_YIELD           = 0._EB
   RN%NU_CO              = 0._EB
   RN%NU_CO2             = 0._EB
   RN%NU_H2              = 0._EB
   RN%NU_H2O             = 0._EB
   RN%NU_O2              = 0._EB
   RN%NU_N2              = 0._EB   
   RN%NU_OTHER           = 0._EB      
   RN%NU_SOOT            = 0._EB
   RN%SOOT_YIELD         = 0._EB
   RN%H2_YIELD           = 0._EB   
   RN%HEAT_OF_COMBUSTION = 0._EB   
   RN%EPUMO2             = 0._EB      
   RN%IDEAL              = .TRUE.
ENDIF SET_MIXTURE_FRACTION

! Set the lower limit of the extinction coefficient

EC_LL = VISIBILITY_FACTOR/MAXIMUM_VISIBILITY

CONTAINS
 

SUBROUTINE SET_REAC_DEFAULTS
 
BOF                         = 0._EB       ! cm**3/mol-s
CO_YIELD                    = 0._EB
CRITICAL_FLAME_TEMPERATURE  = 1427._EB    ! C
E                           = 0._EB       ! kJ/kmol
EPUMO2                      = 13100._EB   ! kJ/kg
FUEL                        = 'null'
H2_YIELD                    = 0._EB
HEAT_OF_COMBUSTION          = -1._EB
HRRPUA_SHEET                = 200._EB     ! kW/m2
ID                          = 'null'
Y_F_INLET                   = 1._EB
N_S                         = -999._EB
NU                          = 0._EB
OXIDIZER                    = 'null'
IF (LES) SOOT_YIELD         = 0.01_EB
IF (DNS) SOOT_YIELD         = 0.0_EB
SOOT_H_FRACTION             = 0.1_EB
X_O2_LL                     = 0.15_EB    ! %
Y_F_LFL                     = 0.0_EB 
Y_F_INLET                   = 1._EB
Y_O2_INFTY                  = 0.23_EB
IDEAL                       = .FALSE.
C                           = 3._EB
H                           = 8._EB
O                           = 0._EB
N                           = 0._EB
OTHER                       = 0._EB
MW_OTHER                    = 28_EB        ! Nitrogen MW
MASS_EXTINCTION_COEFFICIENT = 8700._EB     ! m2/kg
MAXIMUM_VISIBILITY          = 30._EB       ! m
VISIBILITY_FACTOR           = 3._EB
 
END SUBROUTINE SET_REAC_DEFAULTS
 
END SUBROUTINE READ_REAC
 
 
SUBROUTINE PROC_SPEC
 
REAL(EB) :: EPSIJ,SIGMA2,AMW,OMEGA,TSTAR,MU_N,K_N,D_N0,TE,EPSK_N,SIG_N,MW_N, &
            CP_N2,CP_F,CP_CO2,CP_O2,CP_H2O,CP_CO,CP_H2,RSUM_DILUENTS,YSUM_DILUENTS
INTEGER :: IPC,NN,ITMP,N
CHARACTER(30) :: STATE_SPECIES(10)
TYPE(PARTICLE_CLASS_TYPE), POINTER :: PC
LOGICAL :: ABSORBING

! Compute state relations for mixture fraction model
 
IF (MIXTURE_FRACTION) THEN
    CALL STATE_RELATIONSHIPS
   ! For mixture fraction model, get the fuel heat of combustion