read.f

一个开源的火灾动力模拟的系统

C Adjust units of input parameters
C
      IF (DELTAH.GT.0.) DELTAH_FUEL = DELTAH*1000.
      EPUMO2  = EPUMO2*1000.
      TMP_LOWER = TMP_LOWER + TMPM
      TMP_CRIT = CRITICAL_FLAME_TEMPERATURE + TMPM
      DTHRR   = DTSAM
      DTMINT  = DTSAM
      Y_O2_LL = X_O2_LL*32./(X_O2_LL*32. + (1.-X_O2_LL)*28.)
      EC_LL   = VISIBILITY_FACTOR/MAXIMUM_VISIBILITY
      E_GAS   = E*1000.
C
C Set up Mixture Fraction variables
C
      IF (MIXTURE_FRACTION) CALL STATE_RELATIONSHIPS
C
C Set chi_r
C
      CHI_R = RADIATIVE_FRACTION
C
C Get molecular weights
C
      DO N=0,NSPEC
      CALL GAS_PROPS(SPECIES_ID(N),SIG(N),EPSK(N),MWN(N))
      ENDDO
C
C Compute average molecular weight
C
      YY0(0)  = 1.
      DO N=0,NSPEC
      RCON(N) = R0/MWN(N)
      ENDDO
C
      IF (MIXTURE_FRACTION) THEN
C
      RCON_STATE(1) = R0/MW_FUEL
      RCON_STATE(2) = R0/MW_O2
      RCON_STATE(3) = R0/MW_N2
      RCON_STATE(4) = R0/MW_H2O
      RCON_STATE(5) = R0/MW_CO2
      RCON_STATE(6) = R0/MW_CO
      RCON_STATE(7) = R0/MW_H2
      RCON_STATE(8) = R0/MW_C
      RSUM0         = R0/MW_AVG(0)
C
      MW_MIN = MW_AVG(0)
      MW_MAX = MW_AVG(0)
      DO I=1,10000
      MW_MIN = MIN(MW_MIN,MW_AVG(I))
      MW_MAX = MAX(MW_MAX,MW_AVG(I))
      ENDDO
C
      IF (NSPEC.GT.1) THEN
         RCON_MF =  RSUM0
         EXTRA_SPECIES_LOOP: DO N=1,NSPEC
         IF (N.EQ.IFUEL) CYCLE EXTRA_SPECIES_LOOP
         RSUM0 = RSUM0 + YY0(N)*(RCON(N)-RCON_MF)
         MW_MIN = MIN(MW_MIN,MWN(N))
         MW_MAX = MAX(MW_MAX,MWN(N))
         ENDDO EXTRA_SPECIES_LOOP
      ENDIF
C
      RHOA = PINF/(TMPA*RSUM0)
C
      ENDIF
C
      IF (.NOT.MIXTURE_FRACTION) THEN
C
      DO N=1,NSPEC
      YY0(0) = YY0(0) - YY0(N)
      ENDDO
C
      MW_MIN = MWN(0)
      MW_MAX = MWN(0)
      RSUM0  = YY0(0)*R0/MWN(0)
      DO N=1,NSPEC
      MW_MIN = MIN(MW_MIN,MWN(N))
      MW_MAX = MAX(MW_MAX,MWN(N))
      RSUM0 = RSUM0 + YY0(N)*R0/MWN(N)
      ENDDO
C
      RHOA = PINF/(TMPA*RSUM0)
C
      ENDIF
C
C Compute constant-temperature specific heats
C
      CP_GAMMA = RCON(0)*GAMMA/(GAMMA-1.)
      IF (DT0DZ.EQ.-1000000.) DT0DZ= -GRAV/CP_GAMMA
      CPOPR = CP_GAMMA/PR
C
C Compute variable-temperature specific heats for specific species
C
      IF (.NOT.MIXTURE_FRACTION) THEN
C
      SMOKE3D = .FALSE.
C
      SPECIES_LOOP: DO N=0,NSPEC
C
      HH(N,0) = 0.
C
      TEMPERATURE_LOOP: DO J=1,500
      TE = MAX(J*0.01_EB,0.301_EB)
C
      SELECT CASE(SPECIES_ID(N))
      CASE DEFAULT   !  Use Molecular Weight
      CP(N,J) = 0.001*RCON(N)*GAMMA/(GAMMA-1.)
c     CASE('AIR')  ! Nitrogen
c     CP(N,J) = 0.931857+0.293529*TE-0.0705765*TE**2+0.00568836*TE**3+
c    .          0.001589693/TE**2
      CASE('NITROGEN')
      CP(N,J) = 0.931857+0.293529*TE-0.0705765*TE**2+0.00568836*TE**3+
     .          0.001589693/TE**2
      CASE('OXYGEN')
      CP(N,J) = 0.926844+0.191789*TE-0.0370788*TE**2+0.00299313*TE**3-
     .          0.00686447/TE**2
      CASE('METHANE')
      IF (J.LE.130) THEN
      CP(N,J) = -0.0439393+6.77983*TE-2.6576*TE**2+0.366424*TE**3+
     .          0.0424103/TE**2
      ELSE
      CP(N,J) = 5.36326+0.704042*TE-0.132134*TE**2+0.00863688*TE**3-
     .          1.65139/TE**2
      ENDIF
      CASE('CARBON DIOXIDE')
      IF (J.LE.120) THEN
      CP(N,J) = 0.568122+1.25425*TE-0.765713*TE**2+0.180645*TE**3-
     .          0.00310541/TE**2
      ELSE
      CP(N,J) = 1.32196+0.0618199*TE-0.0111884*TE**2+0.00082818*TE**3-
     .          0.146529/TE**2
      ENDIF
      CASE('WATER VAPOR')
      IF (J.LE.170) THEN
      CP(N,J) = 1.67178+0.379583*TE+0.377413*TE**2-0.140804*TE**3+
     .          0.00456328/TE**2
      ELSE
      CP(N,J) = 2.33135+0.479003*TE-0.00833212*TE**2+0.00545106*TE**3-
     .          0.619869/TE**2
      ENDIF
      CASE('HELIUM')
      CP(N,J) = 0.001*RCON(N)*GAMMA/(GAMMA-1.)
      END SELECT
      CP(N,J) = CP(N,J)*1000.            ! J/kg/K
      HH(N,J) = HH(N,J-1) + CP(N,J)*10.  ! J/kg
      RCP(N,J) = 1./CP(N,J)
      ENDDO TEMPERATURE_LOOP
C
      ENDDO SPECIES_LOOP
C
      ENDIF
C
C Compute temperature-dependent specific heat based on mixture fraction
C
      VARIABLE_CP: IF (MIXTURE_FRACTION) THEN
C
      HH(:,0) = 0.
      T_LOOP: DO J=1,500
C
      TE  = MAX(0.01_EB*J,0.301_EB)
C
      CP_O2 = 0.926844+0.191789*TE-0.0370788*TE**2+0.00299313*TE**3-
     .        0.00686447/TE**2
C
      CP_N2 = 0.931857+0.293529*TE-0.0705765*TE**2+0.00568836*TE**3+
     .        0.001589693/TE**2
C
      IF (J.LE.130) THEN   ! Ethylene
      CP_F  = -0.0439393 + 6.77983*TE - 2.6576*TE**2 + 0.366424*TE**3 +
     .         0.0424103/TE**2
      ELSE
      CP_F  = 5.36326+0.704042*TE-0.132134*TE**2+0.00863688*TE**3-
     .        1.65139/TE**2
      ENDIF
C
      IF (J.LE.120) THEN
      CP_CO2 = 0.568122+1.25425*TE-0.765713*TE**2+0.180645*TE**3-
     .         0.00310541/TE**2
      ELSE
      CP_CO2 = 1.32196+0.0618199*TE-0.0111884*TE**2+0.00082818*TE**3-
     .         0.146529/TE**2
      ENDIF
      IF (J.LE.170) THEN
      CP_H2O = 1.67178+0.379583*TE+0.377413*TE**2-0.140804*TE**3+
     .         0.00456328/TE**2
      ELSE
      CP_H2O = 2.33135+0.479003*TE-0.00833212*TE**2+0.00545106*TE**3-
     .         0.619869/TE**2
      ENDIF
C
      Z_LOOP: DO I=0,100
      IYY = 100*I
      CP(I,J) = ( Y_STATE(IYY,1)*CP_F   + 
     .            Y_STATE(IYY,2)*CP_O2  +
     .            Y_STATE(IYY,3)*CP_N2  + 
     .            Y_STATE(IYY,4)*CP_H2O +
     .            Y_STATE(IYY,5)*CP_CO2 )*1000.
      RCP(I,J) = 1./CP(I,J)
C
      HH(I,J) = HH(I,J-1) + CP(I,J)*10.
C
      ENDDO Z_LOOP
C
      ENDDO T_LOOP
C
      ENDIF VARIABLE_CP
C
C For mixture fraction model, get the fuel heat of combustion
C
      IF (MIXTURE_FRACTION) THEN
         COMBUSTION_MODEL = 2
         DELTAH_FUEL = EPUMO2*MW_O2*NU_O2/(MW_FUEL*NU_FUEL)
         DO IPC=1,NPC
         LP=>LAGRANGIAN(IPC)
         IF (LP%DELTAH.GT.0.) 
     .   LP%ADJUST_EVAPORATION = LP%DELTAH/DELTAH_FUEL
         ENDDO
         ENDIF
C
C For finite rate reaction, set parameters
C
      DO N=1,NSPEC
      IF (FUEL.EQ.SPECIES_ID(N)) IFUEL = N
      ENDDO
C
      IF (IFUEL.GT.0 .AND. IOXYGEN.GT.0) THEN
         COMBUSTION_MODEL = 3
         MW_FUEL          = MWN(IFUEL)
         IF (NUN(IFUEL)  .EQ.0.) NUN(IFUEL)   = -1.
         IF (NUN(IFUEL)  .GT.0.) NUN(IFUEL)   = -NUN(IFUEL)
         IF (NUN(IOXYGEN).GT.0.) NUN(IOXYGEN) = -NUN(IOXYGEN)
         IF (NUN(IOXYGEN).EQ.0.) THEN
            WRITE(MESSAGE,'(A)') 
     .     'ERROR: Specify a stoichiometric coefficient for oxygen'
            CALL SHUTDOWN(MESSAGE)
            ENDIF
         EPUMO2 = DELTAH_FUEL*MWN(IFUEL)*NUN(IFUEL)/
     .                       (MWN(IOXYGEN)*NUN(IOXYGEN))
      ENDIF
C
C Compute viscosities and thermal conductivities for species 0-N
C
      IF (.NOT.MIXTURE_FRACTION) THEN
C
      ALLOCATE(MU_SPEC(0:NSPEC,1:500),STAT=IZERO)
      CALL ChkMemErr('READ','MU_SPEC',IZERO)
      ALLOCATE(K_SPEC(0:NSPEC,1:500),STAT=IZERO)
      CALL ChkMemErr('READ','K_SPEC',IZERO)
C
      DO N=0,NSPEC
      SIGMA2 = SIG(N)**2
      DO ITMP=1,500
      TSTAR = ITMP*10/EPSK(N)
      OMEGA = 1.16145*TSTAR**(-0.14874) + 
     .        0.52487*EXP(-0.77320*TSTAR) +
     .        2.16178*EXP(-2.43787*TSTAR)
      MU_N = 26.69E-7*(MWN(N)*ITMP*10)**0.5/(SIGMA2*OMEGA)
      IF (MU_USER(N).GE.0.) MU_N = MU_USER(N)*(ITMP*10/TMPA)**0.75
      K_N = MU_N*CP(N,ITMP)/PR
      IF (K_USER(N).GE.0.)  K_N = K_USER(N)*(ITMP*10/TMPA)**0.75
      MU_SPEC(N,ITMP) = MU_N
      K_SPEC(N,ITMP)  = K_N
      ENDDO
      ENDDO
C
      ALLOCATE(D_SPEC(0:NSPEC,1:500),STAT=IZERO)
      CALL ChkMemErr('READ','D_SPEC',IZERO)
C
      DO N=1,NSPEC
      SIGMA2 = (0.5*(SIG(N)+SIG(0)))**2
      EPSIJ  = SQRT(EPSK(N)*EPSK(0))
      AMW    = SQRT( (MWN(N)+MWN(0))/(2.*MWN(N)*MWN(0)) )
      DO ITMP=1,500
      TSTAR = ITMP*10/EPSIJ
      OMEGA = 1.06036/TSTAR**(0.15610) +
     .        0.19300/EXP(0.47635*TSTAR) +
     .        1.03587/EXP(1.52996*TSTAR) +
     .        1.76474/EXP(3.89411*TSTAR)
      D_N0 = 0.00266E-4*AMW*(10*ITMP)**1.5/(SIGMA2*OMEGA)
      IF (D_USER(N).GE.0.) D_N0 = D_USER(N)*(ITMP*10/TMPA)**1.75
      D_SPEC(N,ITMP) = D_N0
      ENDDO
      ENDDO
C
      ENDIF
C
      IF (MIXTURE_FRACTION) THEN
C
      ALLOCATE(MU_SPEC(0:100,1:500),STAT=IZERO)
      CALL ChkMemErr('READ','MU_SPEC',IZERO) ; MU_SPEC = 0.
      ALLOCATE(K_SPEC(0:100,1:500),STAT=IZERO)
      CALL ChkMemErr('READ','K_SPEC',IZERO)  ; K_SPEC = 0.
      ALLOCATE(D_SPEC(0:100,1:500),STAT=IZERO)
      CALL ChkMemErr('READ','D_SPEC',IZERO)    ; D_SPEC = 0.
C
      STATE_SPECIES(1) = 'ETHYLENE'
      STATE_SPECIES(2) = 'OXYGEN'
      STATE_SPECIES(3) = 'NITROGEN'
      STATE_SPECIES(4) = 'WATER VAPOR'
      STATE_SPECIES(5) = 'CARBON DIOXIDE'
C
      DO N=1,5
      SIG_N = -1. ; EPSK_N = -1. ; MW_N = -.1
      CALL GAS_PROPS(STATE_SPECIES(N),SIG_N,EPSK_N,MW_N)
      SIGMA2 = SIG_N**2
      DO ITMP=1,500
      TSTAR = ITMP*10/EPSK_N
      OMEGA = 1.16145*TSTAR**(-0.14874) +
     .        0.52487*EXP(-0.77320*TSTAR) +
     .        2.16178*EXP(-2.43787*TSTAR)
      MU_N = 26.69E-7*(MW_N*ITMP*10)**0.5/(SIGMA2*OMEGA)
      DO IYY=0,100
      K_N  = MU_N*CP(IYY,ITMP)/PR
      WGT  = Y_STATE(100*IYY,N)*MW_AVG(100*IYY)/MW_N
      MU_SPEC(IYY,ITMP) = MU_SPEC(IYY,ITMP) + WGT*MU_N
      K_SPEC(IYY,ITMP)  = K_SPEC(IYY,ITMP)  + WGT*K_N
      ENDDO
      ENDDO
C
      SIGMA2 = (0.5*(SIG_N+SIG(0)))**2
      EPSIJ  = SQRT(EPSK_N*EPSK(0))
      AMW    = SQRT( (MW_N+MWN(0))/(2.*MW_N*MWN(0)) )
      DO ITMP=1,500
      TSTAR = ITMP*10/EPSIJ
      OMEGA = 1.06036/TSTAR**(0.15610) +
     .        0.19300/EXP(0.47635*TSTAR) +
     .        1.03587/EXP(1.52996*TSTAR) +
     .        1.76474/EXP(3.89411*TSTAR)
      D_N0 = 0.00266E-4*AMW*(10*ITMP)**1.5/(SIGMA2*OMEGA)
      DO IYY=0,100
      D_SPEC(IYY,ITMP) = D_SPEC(IYY,ITMP) + Y_STATE(100*IYY,N)*D_N0
      ENDDO
      ENDDO
      ENDDO
C
      ENDIF
C
C Assign Mass Per Unit Energy
C
      IF (IOXYGEN.GT.0) THEN
      ALLOCATE(MPUE(NSPEC))
      ALLOCATE(EPUM(NSPEC))
      CALL ChkMemErr('READ','MPUE',IZERO)
      MPUE(1:NSPEC) = NUN(1:NSPEC)*MWN(1:NSPEC)/
     .          ABS(EPUMO2*MWN(IOXYGEN)*NUN(IOXYGEN))
      EPUM = 1./MPUE
      ENDIF
C
C Define all possible output quantities
C
      CALL DEFINE_OUTPUT_QUANTITIES
C
      CONTAINS
C 
C
      SUBROUTINE GAS_PROPS(SPECIES,SIGMA,EPSOK,MW)
C
C Molecular weight (g/mol) and Lennard-Jones properties
C
      REAL(EB) SIGMA,EPSOK,MW,SIGMAIN,EPSOKIN,MWIN
      CHARACTER(26) SPECIES
C
      SIGMAIN = SIGMA
      EPSOKIN = EPSOK
      MWIN    = MW
C
      SELECT CASE(SPECIES)
      CASE('AIR')             ; SIGMA=3.711 ; EPSOK= 78.6  ; MW=29.
      CASE('CARBON MONOXIDE') ; SIGMA=3.690 ; EPSOK= 91.7  ; MW=28.
      CASE('CARBON DIOXIDE')  ; SIGMA=3.941 ; EPSOK=195.2  ; MW=44.
      CASE('ETHYLENE')        ; SIGMA=4.163 ; EPSOK=224.7  ; MW=28.
      CASE('HELIUM')          ; SIGMA=2.551 ; EPSOK= 10.22 ; MW= 4.
      CASE('HYDROGEN')        ; SIGMA=2.827 ; EPSOK= 59.7  ; MW= 2.
      CASE('METHANE')         ; SIGMA=3.758 ; EPSOK=148.6  ; MW=16.
      CASE('NITROGEN')        ; SIGMA=3.798 ; EPSOK= 71.4  ; MW=28.
      CASE('OXYGEN')          ; SIGMA=3.467 ; EPSOK=106.7  ; MW=32.
      CASE('PROPANE')         ; SIGMA=5.118 ; EPSOK=237.1  ; MW=44.
      CASE('WATER VAPOR')     ; SIGMA=2.641 ; EPSOK=809.1  ; MW=18.
      CASE DEFAULT            ; SIGMA=3.711 ; EPSOK= 78.6  ; MW=29.
      END SELECT  
C
      IF (SIGMAIN.GT.0.) SIGMA = SIGMAIN
      IF (EPSOKIN.GT.0.) EPSOK = EPSOKIN
      IF (MWIN   .GT.0.) MW    = MWIN 
C
      IF (SPECIES.EQ.'MIXTURE_FRACTION') MW = MW_FUEL
C
      END SUBROUTINE GAS_PROPS
C
C
      SUBROUTINE STATE_RELATIONSHIPS
C
C Calculate state relations, Y_STATE(Z), and average molecular weight.
C Y_STATE(0:10000,I): I=1-Fuel,2-O2,3-N2,4-H2O,5-CO2,6-CO,7-H2,8-SOOT
C
      REAL(EB) Z_TERM,S_TERM,TOTAL_MASS,NU_N2_AIR,NU_N2_DILUENT,ETA,CHI,
     .         Y_N2_INLET
C
C Normalize ideal stoichiometric coefficients
C
      IF (NU_FUEL.NE.1.) THEN
         NU_O2   = NU_O2/NU_FUEL
         NU_H2O  = NU_H2O/NU_FUEL
         NU_CO2  = NU_CO2/NU_FUEL
         NU_FUEL = 1.
         ENDIF
C
      Y_N2_INFTY = 1.-Y_O2_INFTY
C
C If nitrogen is in the fuel, correct MW_FUEL and Y_F_INLET
C
      Y_F_INLET  = Y_F_INLET*(1.-MW_N2*FUEL_N2/MW_FUEL)
      MW_FUEL    = MW_FUEL - MW_N2*FUEL_N2
      Y_N2_INLET = 1.-Y_F_INLET 
C
C Apply Faeth's correlation to get CO yield from soot yield
C
      IF (CO_YIELD.LT.0.) THEN
      GEN_TO_YLD = NU_CO2*MW_C/MW_FUEL
      CO_YIELD   = GEN_TO_YLD*(0.0014 + 0.37*SOOT_YIELD/GEN_TO_YLD)
      ENDIF
C
C Adjust ideal stoichiometric coefficients to account for soot, CO, etc.
C
      NU_SOOT = SOOT_YIELD*MW_FUEL/MW_C
      NU_H2   = H2_YIELD*MW_FUEL/MW_H2