read.f90

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

   IF (CO_PRODUCTION) THEN
      RN => REACTION(2)
   ELSE
      RN => REACTION(1)   
   ENDIF
   DO IPC=1,N_PART
      PC=>PARTICLE_CLASS(IPC)
      IF (PC%HEAT_OF_COMBUSTION > 0._EB) PC%ADJUST_EVAPORATION = PC%HEAT_OF_COMBUSTION/RN%HEAT_OF_COMBUSTION
   ENDDO
ENDIF

 
! Get gas properties
 
DO N=0,N_SPECIES
   ABSORBING = .FALSE.
   CALL GAS_PROPS(SPECIES(N)%NAME,SIG(N),EPSK(N),SPECIES(N)%MW,ABSORBING)
   IF (ABSORBING) SPECIES(N)%ABSORBING = .TRUE.
ENDDO
 
! Compute the initial value of R0/MW_AVG
 
SS0 => SPECIES(0)
 
SS0%YY0 = 1._EB
DO N=1,N_SPECIES
   SS => SPECIES(N)
   IF (SS%MODE==GAS_SPECIES) SS0%YY0 = SS0%YY0 - SS%YY0
ENDDO
 
RSUM0  = SS0%YY0*R0/SS0%MW
SS0%RCON = R0/SS0%MW
IF (MIXTURE_FRACTION) THEN
   RCON_MF(FUEL_INDEX) = R0/REACTION(1)%MW_FUEL
   RCON_MF(O2_INDEX) = R0/MW_O2
   RCON_MF(N2_INDEX) = R0/MW_N2
   RCON_MF(H2O_INDEX) = R0/MW_H2O
   RCON_MF(CO2_INDEX) = R0/MW_CO2
   RCON_MF(CO_INDEX) = R0/MW_CO
   RCON_MF(H2_INDEX) = R0/MW_H2
   RCON_MF(SOOT_INDEX) = R0/MW_SOOT
   RSUM0 = SPECIES(I_FUEL)%RSUM_MF(0)
   SS0%MW = SPECIES(I_FUEL)%MW_MF(0)
   SS0%RCON = SPECIES(I_FUEL)%RSUM_MF(0)   
ENDIF

MW_MIN = SS0%MW
MW_MAX = SS0%MW
 
N_SPEC_DILUENTS = 0
RSUM_DILUENTS = 0._EB
YSUM_DILUENTS = 0._EB
SPECIES_LOOP_0: DO N=1,N_SPECIES
   SS  => SPECIES(N)
   SS%RCON = R0/SS%MW
   IF (SS%MODE==GAS_SPECIES) THEN
      N_SPEC_DILUENTS = N_SPEC_DILUENTS + 1
      MW_MIN = MIN(MW_MIN,SS%MW)
      MW_MAX = MAX(MW_MAX,SS%MW)
      YSUM_DILUENTS = YSUM_DILUENTS + SS%YY0
      RSUM_DILUENTS = RSUM_DILUENTS + SS%YY0*R0/SS%MW
   ELSE
      DO I=0,10000
         MW_MIN = MIN(MW_MIN,SS%MW_MF(I))
         MW_MAX = MAX(MW_MAX,SS%MW_MF(I))
      ENDDO
   ENDIF
ENDDO SPECIES_LOOP_0

IF (MIXTURE_FRACTION) THEN
   RSUM0 = RSUM0 * (1._EB-YSUM_DILUENTS) +RSUM_DILUENTS
ELSE
   RSUM0 = RSUM0 + RSUM_DILUENTS
ENDIF
 
! Compute background density from other background quantities
 
RHOA = P_INF/(TMPA*RSUM0)
 
! Compute constant-temperature specific heats
 
CP_GAMMA = SS0%RCON*GAMMA/(GAMMA-1._EB)
CPOPR = CP_GAMMA/PR
 
! Compute variable-temperature specific heats for specific species
 
SPECIES_LOOP: DO N=0,N_SPECIES
 
   SS => SPECIES(N)
   SS%H_G(0) = 0.

   T_LOOP: DO J=1,500
      TE  = MAX(0.01_EB*J,0.301_EB)
      CP_O2 = 0.926844_EB+0.191789_EB*TE-0.0370788_EB*TE**2+0.00299313_EB*TE**3- 0.00686447_EB/TE**2
      CP_N2 = 0.931857_EB+0.293529_EB*TE-0.0705765_EB*TE**2+0.00568836_EB*TE**3+ 0.001589693_EB/TE**2
      IF (J<=120) THEN   ! Ethylene
         CP_F  = -0.2218014_EB + 6.402844_EB*TE - 3.922632_EB*TE**2 + 0.9894420_EB*TE**3 + 0.01095625_EB/TE**2
      ELSE
         CP_F  = 3.698278_EB+0.4768264_EB*TE-0.0912667_EB*TE**2+0.006062326_EB*TE**3- 0.9078017_EB/TE**2
      ENDIF
      IF (J<=120) THEN
         CP_CO2 = 0.568122_EB+1.25425_EB*TE-0.765713_EB*TE**2+0.180645_EB*TE**3- 0.00310541_EB/TE**2
      ELSE
         CP_CO2 = 1.32196_EB+0.0618199_EB*TE-0.0111884_EB*TE**2+0.00082818_EB*TE**3- 0.146529_EB/TE**2
      ENDIF
      IF (J<=50) THEN
         CP_H2O = 1.95657565_EB
      ELSEIF (50 < J .AND. J <=170) THEN
         CP_H2O = 1.67178_EB+0.379583_EB*TE+0.377413_EB*TE**2-0.140804_EB*TE**3+ 0.00456328_EB/TE**2
      ELSE
         CP_H2O = 2.33135_EB+0.479003_EB*TE-0.00833212_EB*TE**2+0.00545106_EB*TE**3- 0.619869_EB/TE**2
      ENDIF
      IF (J<=130) THEN
         CP_CO = 0.913128_EB+0.217719_EB*TE+0.144809_EB*TE**2-0.0954036_EB*TE**3+ 0.00467932_EB/TE**2
      ELSE
         CP_CO = 1.255382_EB+0.046432_EB*TE-0.00735432_EB*TE**2+0.00483928_EB*TE**3- 0.1172421_EB/TE**2
      ENDIF
      IF (J<=100) THEN
         CP_H2 = 16.53309_EB-5.681708_EB*TE+5.716408_EB*TE**2-1.386437_EB*TE**3- 0.079279_EB/TE**2
      ELSEIF (J > 100 .AND. J <= 250) THEN
         CP_H2 = 9.281542_EB+6.1286785_EB*TE-1.429893_EB*TE**2+0.134119_EB*TE**3+ 0.988995_EB/TE**2
      ELSE
         CP_H2 = 21.70678_EB-2.14654_EB*TE+0.636214_EB*TE**2-0.048438_EB*TE**3-10.266931_EB/TE**2            
      ENDIF
      SELECT CASE (SPECIES(N)%MODE)
         CASE (MIXTURE_FRACTION_SPECIES)
            SS%CP_MF2(FUEL_INDEX,J)  = CP_F*1000._EB
            SS%CP_MF2(O2_INDEX,J)    = CP_O2*1000._EB
            SS%CP_MF2(N2_INDEX,J)    = CP_N2*1000._EB
            SS%CP_MF2(OTHER_INDEX,J) = CP_N2*1000._EB
            SS%CP_MF2(H2O_INDEX,J)   = CP_H2O*1000._EB
            SS%CP_MF2(CO2_INDEX,J)   = CP_CO2*1000._EB
            SS%CP_MF2(CO_INDEX,J)    = CP_CO*1000._EB
            SS%CP_MF2(SOOT_INDEX,J)  = 904._EB !cp carbon
            SS%CP_MF2(H2_INDEX,J)    = CP_H2*1000._EB
            SS%RCP_MF2(:,J)          = 1._EB/SS%CP_MF2(:,J)               

         CASE (GAS_SPECIES)
            SELECT CASE(SPECIES(N)%NAME)
               CASE DEFAULT
                  SS%CP(J) = SS%RCON*GAMMA/(GAMMA-1._EB)
               CASE('AIR')
                  SS%CP(J) = (0.77_EB*CP_N2+0.23_EB*CP_O2)*1000._EB
               CASE('NITROGEN')
                  SS%CP(J) = CP_N2 *1000._EB
               CASE('OXYGEN')
                  SS%CP(J) = CP_O2 *1000._EB
               CASE('METHANE')
                  SS%CP(J) = CP_F  *1000._EB
               CASE('CARBON DIOXIDE')
                  SS%CP(J) = CP_CO2*1000._EB
               CASE('WATER VAPOR')
                  SS%CP(J) = CP_H2O*1000._EB
               CASE('CARBON MONOXIDE')
                  SS%CP(J) = CP_CO*1000._EB
               CASE('HYDROGEN')
                  SS%CP(J) = CP_H2*1000._EB
            END SELECT
            SS%H_G(J) = SS%H_G(J-1) + SS%CP(J)*10._EB  ! J/kg
            SS%RCP(J) = 1._EB/SS%CP(J)
      END SELECT
   ENDDO T_LOOP
ENDDO SPECIES_LOOP

! For finite rate reaction, set parameters
 
FINITE_RATE_REACTION_LOOP: DO NN=1,N_REACTIONS
   RN => REACTION(NN)
   IF (RN%MODE/=FINITE_RATE_REACTION)  CYCLE FINITE_RATE_REACTION_LOOP
   DO N=1,N_SPECIES
      IF (RN%FUEL    ==SPECIES(N)%NAME) THEN
         RN%I_FUEL     = N
         IF (RN%N(N) /=-999._EB) THEN
            RN%N_F        = RN%N(N)
         ELSE
            RN%N_F        = 1._EB         
         ENDIF
      ENDIF
      IF (RN%OXIDIZER==SPECIES(N)%NAME) THEN
         RN%I_OXIDIZER = N
         IF (RN%N(N) /=-999._EB) THEN
            RN%N_O        = RN%N(N)
         ELSE
            RN%N_O        = 1._EB         
         ENDIF
      ENDIF
   ENDDO
   IF (NN==1) I_FUEL = RN%I_FUEL
   RN%MW_FUEL = SPECIES(RN%I_FUEL)%MW
   IF (RN%NU(RN%I_FUEL)     == 0._EB) RN%NU(RN%I_FUEL)     = -1._EB
   IF (RN%NU(RN%I_FUEL)     >  0._EB) RN%NU(RN%I_FUEL)     = -RN%NU(RN%I_FUEL)
   IF (RN%NU(RN%I_OXIDIZER) >  0._EB) RN%NU(RN%I_OXIDIZER) = -RN%NU(RN%I_OXIDIZER)
   IF (RN%NU(RN%I_OXIDIZER) == 0._EB) THEN
      WRITE(MESSAGE,'(A)')  'ERROR: Specify a stoichiometric coefficient for oxidizer'
      CALL SHUTDOWN(MESSAGE)
   ENDIF
   RN%O2_F_RATIO = SPECIES(RN%I_OXIDIZER)%MW*RN%NU(RN%I_OXIDIZER)/(SPECIES(RN%I_FUEL)%MW*RN%NU(RN%I_FUEL))
   RN%EPUMO2     = RN%HEAT_OF_COMBUSTION/RN%O2_F_RATIO
ENDDO FINITE_RATE_REACTION_LOOP
 
! Compute viscosity, thermal conductivity for species 0 to N_SPECIES. Diffusivity for species 1 to N_SPECIES. 
! These terms are used for DNS, and as the lower limits in LES calcs (via SPECIES(0)).
! Source: Poling, Prausnitz and O'Connell. Properties of Gases and Liquids, 5th ed, 2000.
 
SPECIES_LOOP_2: DO N=0,N_SPECIES
!
   SS  => SPECIES(N)
   SS0 => SPECIES(0)
!
   SPEC_MODEIF: IF (SPECIES(N)%MODE==GAS_SPECIES) THEN
!
      SIGMA2 = SIG(N)**2
      DO ITMP=1,500
         TSTAR = ITMP*10/EPSK(N)
         OMEGA = 1.16145_EB*TSTAR**(-0.14874_EB) +  0.52487_EB*EXP(-0.77320_EB*TSTAR) + 2.16178_EB*EXP(-2.43787_EB*TSTAR)
         MU_N = 26.69E-7_EB*(SS%MW*ITMP*10)**0.5_EB/(SIGMA2*OMEGA)
         IF (MU_USER(N)>=0._EB) MU_N = MU_USER(N)*(ITMP*10/TMPA)**0.75_EB
         K_N = MU_N*SS%CP(ITMP)/PR
         IF (K_USER(N)>=0._EB)  K_N = K_USER(N)*(ITMP*10/TMPA)**0.75_EB
         SS%MU(ITMP) = MU_N
         SS%K(ITMP)  = K_N
      ENDDO
      IF (N>0) THEN
         SIGMA2 = (0.5_EB*(SIG(N)+SIG(0)))**2
         EPSIJ  = SQRT(EPSK(N)*EPSK(0))
         AMW    = SQRT( (SS%MW+SS0%MW)/(2._EB*SS%MW*SS0%MW) )
         DO ITMP=1,500
            TSTAR = ITMP*10/EPSIJ
            OMEGA = 1.06036_EB/TSTAR**(0.15610_EB) + 0.19300_EB/EXP(0.47635_EB*TSTAR) + 1.03587_EB/EXP(1.52996_EB*TSTAR) +&
                  1.76474_EB/EXP(3.89411_EB*TSTAR)
            D_N0 = 0.00266E-4_EB*AMW*(10*ITMP)**1.5_EB/(SIGMA2*OMEGA)
            IF (D_USER(N)>=0._EB) D_N0 = D_USER(N)*(ITMP*10/TMPA)**1.75_EB
            SS%D(ITMP) = D_N0
         ENDDO
      ENDIF
 
   ELSE SPEC_MODEIF
      STATE_SPECIES(1) = 'METHANE'
!     STATE_SPECIES(1) = 'ETHYLENE'
      STATE_SPECIES(2) = 'OXYGEN'
      STATE_SPECIES(3) = 'NITROGEN'
      STATE_SPECIES(4) = 'WATER VAPOR'
      STATE_SPECIES(5) = 'CARBON DIOXIDE'
      STATE_SPECIES(6) = 'CARBON MONOXIDE'
      STATE_SPECIES(7) = 'HYDROGEN'
      STATE_SPECIES(8) = 'CARBON MONOXIDE'
      STATE_SPECIES(9) = 'NITROGEN'
      SS%D_MF2  = 0._EB
      SS%MU_MF2 = 0._EB
      SS%K_MF2  = 0._EB
      SUB_SPECIES_LOOP: DO NN=1,9
         SIG_N = -1._EB
         EPSK_N = -1._EB 
         MW_N = -.1_EB
         CALL GAS_PROPS(STATE_SPECIES(NN),SIG_N,EPSK_N,MW_N,ABSORBING)
         SIGMA2 = SIG_N**2
         DO ITMP=1,500
            TSTAR = ITMP*10/EPSK_N
            OMEGA = 1.16145_EB*TSTAR**(-0.14874_EB) + 0.52487_EB*EXP(-0.77320_EB*TSTAR) + 2.16178_EB*EXP(-2.43787_EB*TSTAR)
            MU_N = 26.69E-7_EB*(MW_N*ITMP*10)**0.5_EB/(SIGMA2*OMEGA)
            SS%MU_MF2(NN,ITMP) = MU_N
            SS%K_MF2(NN,ITMP)  = MU_N*SS%CP_MF2(NN,ITMP)/PR
         ENDDO
         SIGMA2 = (0.5_EB*(SIG_N+SIG(0)))**2
         EPSIJ  = SQRT(EPSK_N*EPSK(0))
         AMW    = SQRT( (MW_N+SS0%MW)/(2._EB*MW_N*SS0%MW) )
         DO ITMP=1,500
            TSTAR = ITMP*10/EPSIJ
            OMEGA = 1.06036_EB/TSTAR**(0.15610_EB) + 0.19300_EB/EXP(0.47635_EB*TSTAR) + 1.03587_EB/EXP(1.52996_EB*TSTAR) + &
                    1.76474_EB/EXP(3.89411_EB*TSTAR)
            D_N0 = 0.00266E-4_EB*AMW*(10*ITMP)**1.5_EB/(SIGMA2*OMEGA)
            SS%D_MF2(NN,ITMP) = D_N0
         ENDDO
      ENDDO SUB_SPECIES_LOOP
   ENDIF SPEC_MODEIF
 
ENDDO SPECIES_LOOP_2
 
! Define all possible output quantities
CALL SPECIES_OUTPUT_QUANTITIES
 
CONTAINS
  
 
SUBROUTINE GAS_PROPS(GAS_NAME,SIGMA,EPSOK,MW,ABSORBING)
 
! Molecular weight (g/mol) and Lennard-Jones properties
 
REAL(EB) SIGMA,EPSOK,MW,SIGMAIN,EPSOKIN,MWIN
CHARACTER(30) GAS_NAME
LOGICAL ABSORBING
 
SIGMAIN = SIGMA
EPSOKIN = EPSOK
MWIN    = MW
 
SELECT CASE(GAS_NAME)
   CASE('AIR')             
      SIGMA=3.711_EB 
      EPSOK= 78.6_EB  
      MW=28.8_EB
   CASE('CARBON MONOXIDE') 
      SIGMA=3.690_EB 
      EPSOK= 91.7_EB  
      MW=28._EB
      ABSORBING = .TRUE.
   CASE('CARBON DIOXIDE')  
      SIGMA=3.941_EB 
      EPSOK=195.2_EB  
      MW=44._EB
      ABSORBING = .TRUE.      
   CASE('ETHYLENE')        
      SIGMA=4.163_EB 
      EPSOK=224.7_EB 
      MW=28._EB
      ABSORBING = .TRUE.
   CASE('HELIUM')          
      SIGMA=2.551_EB 
      EPSOK= 10.22_EB 
      MW= 4._EB
   CASE('HYDROGEN')        
      SIGMA=2.827_EB 
      EPSOK= 59.7_EB
      MW= 2._EB
   CASE('METHANE')         
      SIGMA=3.758_EB 
      EPSOK=148.6_EB  
      MW=16._EB
      ABSORBING = .TRUE.
   CASE('NITROGEN')        
      SIGMA=3.798_EB 
      EPSOK= 71.4_EB  
      MW=28._EB
   CASE('OXYGEN')          
      SIGMA=3.467_EB 
      EPSOK=106.7_EB  
      MW=32._EB
   CASE('PROPANE')         
      SIGMA=5.118_EB 
      EPSOK=237.1_EB  
      MW=44._EB
      ABSORBING = .TRUE.
   CASE('WATER VAPOR')     
      SIGMA=2.641_EB 
      EPSOK=809.1_EB  
      MW=18._EB
      ABSORBING = .TRUE.
   CASE DEFAULT            
      SIGMA=3.711_EB 
      EPSOK= 78.6_EB
      MW=28.8_EB
END SELECT  
 
IF (SIGMAIN>0._EB) SIGMA = SIGMAIN
IF (EPSOKIN>0._EB) EPSOK = EPSOKIN
IF (MWIN   >0._EB) MW    = MWIN 
 
IF (GAS_NAME=='MIXTURE_FRACTION') MW = REACTION(1)%MW_FUEL
 
END SUBROUTINE GAS_PROPS
 
 
SUBROUTINE STATE_RELATIONSHIPS
 
! Calculate state relations, SPECIES(N)%Y_MF, and average molecular weight, SPECIES(N)%MW_MF, for major species tied to MIX_FRAC
! Y_MF(0:10000,I): I=1-Fuel,2-O2,3-N2,4-H2O,5-CO2,6-CO,7-H2,8-SOOT
 
REAL(EB) :: TOTAL_MASS,FUEL_BURNED,ZZ,DZZ,FUEL_TO_CO2
INTEGER :: N,REAC_INDEX(1:N_REACTIONS)

!Calc heat of formation if ideal heat of combustion data was input on REAC
REAC_INDEX = 0._EB
DO N = 1, N_SPECIES
   IF (SPECIES(N)%MODE==MIXTURE_FRACTION_SPECIES) REAC_INDEX(SPECIES(N)%REAC_INDEX) = N
ENDDO

REACTION_LOOP: DO N=1,N_REACTIONS
   RN => REACTION(N)
   SS => SPECIES(REAC_INDEX(N))

  ! Miscellaneous constants
 
   RN%Y_N2_INFTY = 1._EB-RN%Y_O2_INFTY  ! Assumes that air is made up of only oxygen and nitrogen
   RN%Y_N2_INLET = 1._EB-RN%Y_F_INLET   ! Assumes that the fuel is only diluted by nitrogen
   RN%Z_F_CONS = RN%NU_O2 * MW_O2/RN%MW_FUEL * (1._EB + RN%Y_N2_INFTY/RN%Y_O2_INFTY)
 
! Stoichiometric value of mixture fraction
 
   RN%Z_F = RN%Y_O2_INFTY/(RN%Y_O2_INFTY+RN%Y_F_INLET*RN%O2_F_RATIO)
 
! Compute Y_MF and MW_MF
   SS%Z_MAX = 0._EB
   SS%Y_MF(10000,:) = 0._EB
   SS%Y_MF(10000,FUEL_INDEX) = RN%Y_F_INLET
   SS%Y_MF(10000,N2_INDEX) = RN%Y_N2_INLET
   SS%MW_MF(10000) = 1/(RN%Y_F_INLET/RN%MW_FUEL+RN%Y_N2_INLET/MW_N2)
   DZZ = 1.E-4_EB
 
   Z_LOOP: DO I=0,9999
 
      ZZ=REAL(I,EB)*DZZ
      IF (RN%NU_O2 == 0) THEN
         SS%Y_MF(I,FUEL_INDEX)    =  ZZ*RN%Y_F_INLET
      ELSE
         SS%Y_MF(I,FUEL_INDEX)    =  ZZ*RN%Y_F_INLET - MIN( ZZ*RN%Y_F_INLET , (1._EB-ZZ)*RN%Y_O2_INFTY/RN%O2_F_RATIO )
      ENDIF
      SS%Y_MF(I,FUEL_INDEX)  =