Merge pull request #14119 from drjfloyd/master

FDS Source: Fix enthalpy adjust for N_O2 plus some TGA enhancements

Author: drjfloyd

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -3373,6 +3373,8 @@ \subsection{Simulating Bench-scale Measurements like the TGA, DSC, and MCC}
 
 The result of the \ct{TGA_ANALYSIS} is a single comma-delimited file called \ct{CHID_tga.csv}. The first and second columns of the file consist of the time and sample temperature. The third column is the normalized sample mass; that is, the sample mass divided by its initial mass. The following columns list the mass fractions of the individual material components. The next column is the total mass loss rate, in units of s$^{-1}$, followed by the mass loss rates of the individual material components. The next column is the heat release rate per unit mass of the sample in units of W/g, typical of an MCC measurement. The final column is the rate of heat absorbed by the sample normalized by its {\it original} mass, also in units of W/g, typical of a DSC measurement. 
 
+The timestep used in the TGA analysis can be controlled with \ct{TGA_DT} on the \ct{SURF} line. The output spacing for temperature can be set with \ct{TGA_DUMP} on the \ct{SURF} line.
+
 Details of the output quantities are discussed in Sec.~\ref{info:material_components}. Further details on these measurement techniques and how to interpret them are found in the FDS Verification Guide~\cite{FDS_Verification_Guide}.
 
 \subsubsection{Example}
@@ -3428,6 +3430,8 @@ \section{Pyrolysis and Energy Conservation}
 
 If the sum of the yields is less than 1, then for the purpose of solving for the $H_{\rm adj}$ values, FDS will assume that the missing mass is a material with the same specific heat as the original material.
 
+If a material has a reaction where \ct{N_O2} is set and the adjustment process needs to determine the \ct{REFERENCE_TEMPERATURE}, then a value for the oxygen concentration is needed. This value will default to the ambient oxygen concentration for the simulation; however, it can be overridden on a material basis by specifying \ct {X_O2_PYRO} for that material. This may be needed if, for example, the kinetics and heat of reaction(s) for a material were developed using a specific oxygen concentration.
+
 \chapter{Ventilation}
 
 This chapter explains how to model a ventilation system. There are two ways to do this. First, is to explicitly specify air flow rates into and out of
@@ -12511,6 +12515,7 @@ \section{\texorpdfstring{{\tt MATL}}{MATL} (Material Properties)}
 \ct{SPECIFIC_HEAT_RAMP}          & Character       & Section~\ref{info:thermal_properties} &                   &               \\ \hline
 \ct{SPEC_ID(:,:)}                 & Char.~Array     & Section~\ref{info:solid_pyrolysis}    &                   &               \\ \hline
 \ct{SURFACE_OXIDATION_MODEL}     & Logical         & Section~\ref{vegetation_model}        &                   & \ct{F}       \\ \hline
+\ct{X_O2_PYRO}                    & Logical         & Section~\ref{solid_phase_energy_conservation}       &                   &        \\ \hline
 \end{longtable}
 
 
@@ -13491,7 +13496,9 @@ \section{\texorpdfstring{{\tt SURF}}{SURF} (Surface Properties)}
 \ct{TEXTURE_HEIGHT}                                     & Real            & Section~\ref{info:texture_map}            & m                   & 1.                      \\ \hline
 \ct{TEXTURE_MAP}                                        & Character       & Section~\ref{info:texture_map}            &                     &                         \\ \hline
 \ct{TEXTURE_WIDTH}                                      & Real            & Section~\ref{info:texture_map}            & m                   & 1.                      \\ \hline
-\ct{TGA_ANALYSIS}                                       & Logical         & Section~\ref{info:TGA_DSC_MCC}            &                     & \ct{F}                 \\ \hline
+\ct{TGA_ANALYSIS}                                       & Logical         & Section~\ref{info:TGA_DSC_MCC}            &                     & \ct{F}                  \\ \hline
+\ct{TGA_DT}                                             & Real            & Section~\ref{info:TGA_DSC_MCC}            & s                   & 0.1                     \\ \hline
+\ct{TGA_DUMP}                                           & Real            & Section~\ref{info:TGA_DSC_MCC}            & K                   & 1                       \\ \hline
 \ct{TGA_FINAL_TEMPERATURE}                             & Real            & Section~\ref{info:TGA_DSC_MCC}            & $^\circ$C           & 800.                    \\ \hline
 \ct{TGA_HEATING_RATE}                                  & Real            & Section~\ref{info:TGA_DSC_MCC}            & $^\circ$C/min       & 5.                      \\ \hline
 \ct{THICKNESS(:)}                                        & Real Array      & Section~\ref{info:SURF_MATL_Basics}       & m                   &                         \\ \hline

Source/cons.f90

@@ -743,6 +743,8 @@ MODULE GLOBAL_CONSTANTS
 INTEGER :: TGA_SURF_INDEX=-100               !< Surface properties to use for special TGA calculation
 INTEGER :: TGA_WALL_INDEX=-100               !< Wall index to use for special TGA calculation
 INTEGER :: TGA_PARTICLE_INDEX=-100           !< Particle index to use for special TGA calculation
+REAL(EB) :: TGA_DT=0.1_EB                    !< Time step (s) to use for special TGA calculation
+REAL(EB) :: TGA_DUMP=1._EB                   !< Temperature output interval (K), starting at TMPA, to use for special TGA calculation
 REAL(EB) :: TGA_HEATING_RATE=5._EB           !< Heat rate (K/min) to use for special TGA calculation
 REAL(EB) :: TGA_FINAL_TEMPERATURE=800._EB    !< Final Temperature (C) to use for special TGA calculation
 

Source/func.f90

@@ -2830,15 +2830,16 @@ SUBROUTINE GET_TMP_REF(N_MATL,NR)
 ELSE
    RR_MAX = 0._EB
 ENDIF
+
 DO WHILE (INT(TMP)<I_MAX_TEMP)
    TMP = TMP + DTDT * DT
-   REACTION_RATE = ML%A(NR)*RHO_S**ML%N_S(NR)*EXP(-ML%E(NR)/(R0*TMP))*TMP**ML%N_T(NR)
+   REACTION_RATE = ML%A(NR)*RHO_S**ML%N_S(NR)*EXP(-ML%E(NR)/(R0*TMP))*TMP**ML%N_T(NR)*ML%X_O2_PYRO**ML%N_O2(NR)
    IF (REACTION_RATE > RR_MAX) THEN
       ML%TMP_REF(NR) = TMP
       RR_MAX = REACTION_RATE
    ENDIF
    RHO_S = RHO_S - REACTION_RATE * DT
-   IF (RHO_S<0._EB) EXIT
+   IF (RHO_S<TWO_EPSILON_EB) EXIT
 ENDDO
 
 END SUBROUTINE GET_TMP_REF

Source/init.f90

@@ -2894,15 +2894,15 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)
 USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_STORAGE
 USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
 USE COMP_FUNCTIONS, ONLY: SHUTDOWN
-USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT
+USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT,GET_MASS_FRACTION
 USE CONTROL_VARIABLES, ONLY : CONTROL
 USE DEVICE_VARIABLES, ONLY : DEVICE
 INTEGER, INTENT(IN) :: I,J,K,NM,OBST_INDEX,IW,IOR,SURF_INDEX
 INTEGER  :: NOM_FOUND,NOM=0,ITER,IIO_MIN,IIO_MAX,JJO_MIN,JJO_MAX,KKO_MIN,KKO_MAX,VENT_INDEX
 INTEGER, INTENT(OUT) :: IERR
 REAL(EB), INTENT(IN) :: TT
 REAL(EB) :: PX,PY,PZ,T_ACTIVATE,XIN,YIN,ZIN,DIST,XW,YW,ZW,RDN,AW,TSI,&
-            ZZ_GET(1:N_TRACKED_SPECIES),RSUM_F,R1,RR,DELTA
+            ZZ_GET(1:N_TRACKED_SPECIES),RSUM_F,R1,RR,DELTA,Y_O2_F
 INTEGER  :: N,SURF_INDEX_NEW,IIG,JJG,KKG,IIO,JJO,KKO,IC,ICG,ICO,NOM_CHECK(0:1),BOUNDARY_TYPE,FI,VENT_INDEX_FOUND
 LOGICAL :: ALIGNED
 TYPE (MESH_TYPE), POINTER :: M,MM
@@ -3312,6 +3312,11 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)
    M%RSUM(I,J,K) = M%RSUM(IIG,JJG,KKG)
    B1%ZZ_F(1:N_TRACKED_SPECIES)  = M%ZZ(IIG,JJG,KKG,1:N_TRACKED_SPECIES)
    M%ZZ(I,J,K,1:N_TRACKED_SPECIES) = M%ZZ(IIG,JJG,KKG,1:N_TRACKED_SPECIES)
+   IF (O2_INDEX > 0) THEN
+      ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) 
+      CALL GET_MASS_FRACTION(ZZ_GET,O2_INDEX,Y_O2_F)
+      B2%Y_O2_F = Y_O2_F
+   ENDIF
 ENDIF
 
 ! Compute the mass of the grid cell corresponding to the wall cell

Source/read.f90

@@ -6871,11 +6871,12 @@ END SUBROUTINE PROC_PROP
 SUBROUTINE READ_MATL
 
 USE COMP_FUNCTIONS, ONLY : SEARCH_INPUT_FILE
+USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION
 CHARACTER(LABEL_LENGTH) :: CONDUCTIVITY_RAMP,SPECIFIC_HEAT_RAMP
 CHARACTER(LABEL_LENGTH) :: SPEC_ID(MAX_SPECIES,MAX_REACTIONS),PART_ID(MAX_LPC,MAX_REACTIONS)
 REAL(EB) :: EMISSIVITY,CONDUCTIVITY,SPECIFIC_HEAT,DENSITY,ABSORPTION_COEFFICIENT,BOILING_TEMPERATURE, &
-            PEAK_REACTION_RATE,MW,&
-            REFERENCE_ENTHALPY,REFERENCE_ENTHALPY_TEMPERATURE,REAC_RATE_DELTA
+            PEAK_REACTION_RATE,MW,ZZ_GET(1:N_TRACKED_SPECIES),&
+            REFERENCE_ENTHALPY,REFERENCE_ENTHALPY_TEMPERATURE,REAC_RATE_DELTA,X_O2_PYRO
 REAL(EB), DIMENSION(MAX_MATERIALS,MAX_REACTIONS) :: NU_MATL
 REAL(EB), DIMENSION(MAX_REACTIONS) :: A,E,HEATING_RATE,PYROLYSIS_RANGE,HEAT_OF_REACTION, &
           N_S,N_T,N_O2,REFERENCE_RATE,REFERENCE_TEMPERATURE, &
@@ -6894,7 +6895,7 @@ SUBROUTINE READ_MATL
                 MAX_REACTION_RATE,MW,N_O2,N_REACTIONS,N_S,N_T,NU_MATL,NU_PART,NU_SPEC,PART_ID,&
                 PYROLYSIS_RANGE,REAC_RATE_DELTA,REFERENCE_ENTHALPY,REFERENCE_ENTHALPY_TEMPERATURE,&
                 REFERENCE_RATE,REFERENCE_TEMPERATURE,&
-                SPECIFIC_HEAT,SPECIFIC_HEAT_RAMP,SPEC_ID,SURFACE_OXIDATION_MODEL
+                SPECIFIC_HEAT,SPECIFIC_HEAT_RAMP,SPEC_ID,SURFACE_OXIDATION_MODEL,X_O2_PYRO
 
 ! Count the MATL lines in the input file
 
@@ -7118,6 +7119,7 @@ SUBROUTINE READ_MATL
    ALLOCATE(ML%N_LPC(N_REACTIONS),STAT=IZERO)
    CALL ChkMemErr('READ','ML%N_LPC',IZERO)
    ML%N_LPC=0
+   ML%X_O2_PYRO = X_O2_PYRO
 
    ! Decide which pyrolysis model to use
 
@@ -7145,6 +7147,12 @@ SUBROUTINE READ_MATL
       IF (O2_INDEX <= 0) THEN
          WRITE(MESSAGE,'(A,A,A)') 'ERROR(256): MATL ',TRIM(ID),', oxidation reaction set but OXYGEN not a defined species.'
          CALL SHUTDOWN(MESSAGE) ; RETURN
+      ELSE
+         IF (ML%X_O2_PYRO < 0._EB) THEN
+            ZZ_GET(1:N_TRACKED_SPECIES) = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0
+            CALL GET_MASS_FRACTION(ZZ_GET,O2_INDEX,X_O2_PYRO)
+            ML%X_O2_PYRO = X_O2_PYRO
+         ENDIF
       ENDIF
       OXIDATION_REACTION = .TRUE.
    ENDIF
@@ -7316,6 +7324,7 @@ SUBROUTINE SET_MATL_DEFAULTS
 SURFACE_OXIDATION_MODEL= .FALSE.
 HEATING_RATE           = 5._EB       ! K/min
 PYROLYSIS_RANGE        = 80._EB      ! K or C
+X_O2_PYRO              = -1._EB
 
 END SUBROUTINE SET_MATL_DEFAULTS
 
@@ -7748,7 +7757,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
                 RGB,ROUGHNESS,SHAPE_FACTOR,SPEC_ID,&
                 SPREAD_RATE,STRETCH_FACTOR,SURFACE_VOLUME_RATIO,&
                 TAU_EF,TAU_MF,TAU_PART,TAU_Q,TAU_T,TAU_V,TEXTURE_HEIGHT,TEXTURE_MAP,TEXTURE_WIDTH,&
-                TGA_ANALYSIS,TGA_FINAL_TEMPERATURE,TGA_HEATING_RATE,THICKNESS,TIME_STEP_FACTOR,&
+                TGA_ANALYSIS,TGA_DT,TGA_DUMP,TGA_FINAL_TEMPERATURE,TGA_HEATING_RATE,THICKNESS,TIME_STEP_FACTOR,&
                 TMP_BACK,TMP_FRONT,TMP_FRONT_INITIAL,TMP_GAS_BACK,TMP_GAS_FRONT,TMP_INNER,TRANSPARENCY,&
                 VEG_LSET_BETA,VEG_LSET_CHAR_FRACTION,VEG_LSET_FIREBASE_TIME,VEG_LSET_FUEL_INDEX,VEG_LSET_HT,VEG_LSET_IGNITE_TIME,&
                 VEG_LSET_M1,VEG_LSET_M10,VEG_LSET_M100,VEG_LSET_MLW,VEG_LSET_MLH,VEG_LSET_QCON,&

Source/type.f90

@@ -776,6 +776,7 @@ MODULE TYPES
    REAL(EB) :: REFERENCE_ENTHALPY_TEMPERATURE           !< Temperature for the reference enthalpy (K)
    REAL(EB) :: REAC_RATE_DELTA                          !< Allowable jump condition for pyrolysis rate. Determines RENODE_DELTA_T
    REAL(EB) :: RENODE_DELTA_T=1.E9_EB                   !< Temperature change (K) above which no resmeshing
+   REAL(EB) :: X_O2_PYRO                                !< Oxygen mass fraction to use for enthalpy adjustment
    INTEGER :: PYROLYSIS_MODEL                           !< Type of pyrolysis model (SOLID, LIQUID, VEGETATION)
    CHARACTER(LABEL_LENGTH) :: ID                        !< Identifier
    CHARACTER(LABEL_LENGTH) :: RAMP_K_S                  !< Name of RAMP for thermal conductivity of solid

Source/wall.f90

@@ -3802,7 +3802,7 @@ SUBROUTINE TGA_ANALYSIS(NM)
 
 USE PHYSICAL_FUNCTIONS, ONLY: SURFACE_DENSITY
 USE COMP_FUNCTIONS, ONLY: SHUTDOWN
-REAL(EB) :: DT_TGA=0.01_EB,T_TGA,SURF_DEN_0,HRR
+REAL(EB) :: T_TGA,SURF_DEN_0,HRR,TMP_DUMP
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: SURF_DEN
 INTEGER, INTENT(IN) :: NM
 INTEGER :: N_TGA,I,IW,IP,N,NR
@@ -3820,7 +3820,7 @@ SUBROUTINE TGA_ANALYSIS(NM)
 RADIATION = .FALSE.
 TGA_HEATING_RATE = TGA_HEATING_RATE/60._EB  ! K/min --> K/s
 TGA_FINAL_TEMPERATURE = TGA_FINAL_TEMPERATURE + TMPM  ! C --> K
-N_TGA = NINT((TGA_FINAL_TEMPERATURE-TMPA)/(TGA_HEATING_RATE*DT_TGA))
+N_TGA = NINT((TGA_FINAL_TEMPERATURE-TMPA)/(TGA_HEATING_RATE*TGA_DT))
 T_TGA = 0._EB
 
 IF (TGA_WALL_INDEX>0) THEN
@@ -3846,17 +3846,18 @@ SUBROUTINE TGA_ANALYSIS(NM)
 
 SURF_DEN_0 = SF%SURFACE_DENSITY
 WRITE(TCFORM,'(A,I3.3,5A)') "(",2*SF%N_MATL+5,"(",TRIM(FMT_R),",','),",TRIM(FMT_R),")"
-
+TMP_DUMP = TMPA + TGA_DUMP
 DO I=1,N_TGA
    IF (ONE_D%LAYER_THICKNESS(1)<TWO_EPSILON_EB) EXIT
-   T_TGA = I*DT_TGA
+   T_TGA = REAL(I,EB)*TGA_DT
    B1%TMP_G = TMPA + TGA_HEATING_RATE*T_TGA
    IF (TGA_WALL_INDEX>0) THEN
-      CALL SOLID_HEAT_TRANSFER(NM,T_TGA,DT_TGA,WALL_INDEX=IW)
+      CALL SOLID_HEAT_TRANSFER(NM,T_TGA,TGA_DT,WALL_INDEX=IW)
    ELSE
-      CALL SOLID_HEAT_TRANSFER(NM,T_TGA,DT_TGA,PARTICLE_INDEX=IP)
+      CALL SOLID_HEAT_TRANSFER(NM,T_TGA,TGA_DT,PARTICLE_INDEX=IP)
    ENDIF
-   IF (I==1 .OR. MOD(I,NINT(1._EB/(TGA_HEATING_RATE*DT_TGA)))==0) THEN
+   IF (I==1 .OR. B1%TMP_G >= TMP_DUMP) THEN
+      IF (B1%TMP_G >= TMP_DUMP) TMP_DUMP = TMP_DUMP + TGA_DUMP
       IF (TGA_WALL_INDEX>0) THEN
          SURF_DEN(0) = SURFACE_DENSITY(0,SF,ONE_D)
          DO N=1,SF%N_MATL