Merge pull request #15046 from cxp484/master

FDS Source: Adiabatic Flame Temperature calculation in Mixing zone

Author: cxp484

Source/chem.f90

@@ -783,7 +783,8 @@ END FUNCTION DDTMP_TROE
 !> \brief cvode interface for ODE integrator. Call sundials cvode in serial mode.
 !> \param CC species concentration (kmol/m3) array
 !> \param ZZ_0 initial species mass fraction array (of the unbuned zone),needed for mixing+chem
-!> \param TMP_IN is the temperature of the cell (unburned zone)
+!> \param TMP_IN is the temperature of the mixed zone (For DNS it is cell temp)
+!> \param TMP_UNMIX is the unmix zone temperature
 !> \param PR_IN is the pressure
 !> \param ZETA0 is the initial unmixed fraction 
 !> \param TAU_MIX is Mixing timescale
@@ -792,13 +793,13 @@ END FUNCTION DDTMP_TROE
 !> \param TEND is the end time in seconds
 !> \param RTOL is the relative error for all the species (REAL_EB)
 !> \param ATOL is the absolute error tolerance array for the species (REAL_EB)
-!> \param TMP_OUT reactor calculated temperature at the end
+!> \param TMP_OUT reactor calculated mixing zone temperature at the end
 !> \param CHEM_TIME Chemical time scale
 !> \param WRITE_SUBSTEPS Whether to write cvode substeps. Only write for first cfd step.
 !> \param CVODE_CALL_OPTION 1:CV_NORMAL, 2=CV_ONE_STEP
 !> \details This is the interface subroutine to the other modules.
 
-SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,TEND, RTOL, ATOL, &
+SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,TEND, RTOL, ATOL, &
                         TMP_OUT, CHEM_TIME, WRITE_SUBSTEPS, CVODE_CALL_OPTION)
 USE PHYSICAL_FUNCTIONS, ONLY : MOLAR_CONC_TO_MASS_FRAC, CALC_EQUIV_RATIO, GET_ENTHALPY, GET_MOLECULAR_WEIGHT
 USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER
@@ -813,7 +814,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
 USE FSUNDIALS_NVECTOR_MOD      ! FORTRAN INTERFACE TO GENERIC N_VECTOR
 
 REAL(EB), INTENT(INOUT) :: CC(N_TRACKED_SPECIES)
-REAL(EB), INTENT(IN)    :: ZZ_0(N_TRACKED_SPECIES),TMP_IN,PR_IN,ZETA0,TAU_MIX,CELL_MASS,TCUR,TEND
+REAL(EB), INTENT(IN)    :: ZZ_0(N_TRACKED_SPECIES),TMP_IN,TMP_UNMIX,PR_IN,ZETA0,TAU_MIX,CELL_MASS,TCUR,TEND
 REAL(EB), INTENT(IN)    :: ATOL(N_TRACKED_SPECIES)
 REAL(EB), INTENT(IN)    :: RTOL
 REAL(EB), INTENT(OUT)   :: TMP_OUT,CHEM_TIME
@@ -841,7 +842,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
 TYPE(C_PTR)                    :: USERDATAPTR  ! USER DATA CONTAINS MIXING INFORMATION
 
 REAL(EB) :: ZZ(N_TRACKED_SPECIES), EQUIV, H_IN
-INTEGER :: CVODE_TASK, NS, NTRY, MAXTRY, SUBSTEP_COUNT
+INTEGER :: CVODE_TASK, NS, NTRY, MAXTRY, SUBSTEP_COUNT, MAXTRY_FAC
 REAL(EB) :: H_G
 TYPE(USERDATA), TARGET :: USER_DATA
 LOGICAL :: ONLY_FIRST_STEP=.TRUE. ! Needed in CV_ONE_STEP
@@ -948,7 +949,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
 USER_DATA%CELL_MASS = CELL_MASS
 ALLOCATE(USER_DATA%ZZ_0(N_TRACKED_SPECIES))
 USER_DATA%ZZ_0 = ZZ_0
-CALL GET_ENTHALPY(ZZ_0,H_IN,TMP_IN)
+CALL GET_ENTHALPY(ZZ_0,H_IN,TMP_UNMIX)
 USER_DATA%H_IN = H_IN 
 USERDATAPTR = C_LOC(USER_DATA)
 IERR_C = FCVODESETUSERDATA(CVODE_MEM, USERDATAPTR)
@@ -987,12 +988,14 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
       ENDIF  
    END DO
 
-
 ENDIF
 
 IF (IERR_C /= 0) THEN
-   MAXTRY = CVODE_MAX_TRY
    ! If all internal substeps are taken try two more times. This will allow larger CFD timestep.
+   ! Make MAXTRY at least CVODE_MAX_TRY or for larger timestep (>1E-3) scale it proportionaly.
+   MAXTRY_FAC = CEILING((TEND-TCUR)/1.0E-3_EB)
+   MAXTRY_FAC = MIN(MAX(MAXTRY_FAC,1),50)
+   MAXTRY = CVODE_MAX_TRY*MAXTRY_FAC
    IF (IERR_C == CV_TOO_MUCH_WORK) THEN !CV_TOO_MUCH_WORK == all internal substeps are taken
       NTRY = 0
       DO WHILE (NTRY < MAXTRY)
@@ -1008,8 +1011,8 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,
 
    IF (IERR_C .NE. CV_SUCCESS) THEN
       IF (IERR_C == CV_TOO_MUCH_WORK) THEN
-         WRITE(LU_ERR,'(A, 2E18.8, A)')" WARN: CVODE took all internal substeps. CUR_CFD_TIME, DT=", CUR_CFD_TIME, (TEND-TCUR), &
-                     ". If the warning persists, reduce the timestep."
+         WRITE(LU_ERR,'(A, 2E18.8, I8, A)')" WARN: CVODE took all internal substeps. CUR_CFD_TIME, DT, MAXTRY=", CUR_CFD_TIME, &
+                     (TEND-TCUR), MAXTRY, ". If the warning persists, reduce the timestep."
       ELSE
          WRITE(LU_ERR,'(A, I4, A, 2E18.8, A)')" WARN: CVODE didn't finish ODE solution with message code:", IERR_C, &
                         " and CUR_CFD_TIME, DT=", CUR_CFD_TIME, (TEND-TCUR), ". If the warning persists, reduce the timestep."

Source/cons.f90

@@ -886,14 +886,25 @@ END MODULE RADCONS
 MODULE CHEMCONS
 USE PRECISION_PARAMETERS
 
+!> \brief Parameters associated with IGNITION_ZONES
+TYPE IGNITION_ZONE_TYPE
+   REAL(EB) :: X1            !< Lower x bound of Ignition Zone
+   REAL(EB) :: X2            !< Upper x bound of Ignition Zone
+   REAL(EB) :: Y1            !< Lower y bound of Ignition Zone
+   REAL(EB) :: Y2            !< Upper y bound of Ignition Zone
+   REAL(EB) :: Z1            !< Lower z bound of Ignition Zone
+   REAL(EB) :: Z2            !< Upper z bound of Ignition Zone
+   INTEGER :: DEVC_INDEX=0   !< Index of device controlling the status of the zone
+   CHARACTER(LABEL_LENGTH) :: DEVC_ID='null'  !< Name of device controlling the status of the zone
+END TYPE IGNITION_ZONE_TYPE
+
 INTEGER, ALLOCATABLE, DIMENSION(:) :: YP2ZZ
 REAL(EB) :: ODE_MIN_ATOL= -1._EB
 LOGICAL  :: EQUIV_RATIO_CHECK = .FALSE.
-REAL(EB) :: MIN_EQUIV_RATIO=0.0_EB
+REAL(EB) :: MIN_EQUIV_RATIO=0.1_EB
 REAL(EB) :: MAX_EQUIV_RATIO=20.0_EB
 LOGICAL  :: DO_CHEM_LOAD_BALANCE = .FALSE.
 INTEGER  :: MAX_CVODE_SUBSTEPS=100000
-REAL(EB) :: MAX_CHEM_TIME=1.E-6_EB
 INTEGER  :: CVODE_MAX_TRY=4
 LOGICAL  :: IS_EXPONENT_LT_1 = .FALSE.
 
@@ -902,5 +913,18 @@ MODULE CHEMCONS
 REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: CVODE_SUBSTEP_DATA
 INTEGER :: TOTAL_SUBSTEPS_TAKEN
 
-END MODULE CHEMCONS
+! Adiabatic flame temperature calculation
+CHARACTER(LABEL_LENGTH) :: FUEL_ID_FOR_AFT='null'
+INTEGER :: I_FUEL,I_CO2,I_H2O,I_O2,I_N2 ! Store the index of the species in the ZZ array.
+LOGICAL  :: USE_MIXED_ZN_AFT_TMP = .TRUE.
+
+! Mixing
+REAL(EB) :: ZETA_ARTIFICAL_MIN_LIMIT=0.99_EB
+REAL(EB) :: ZETA_ARTIFICAL_MAX_LIMIT=0.9999_EB
+REAL(EB) :: ZETA_FIRST_STEP_DIV=10._EB
 
+! IGNITION ZONES (mainly for premixed flame)
+INTEGER :: N_IGNITION_ZONES = 0
+TYPE(IGNITION_ZONE_TYPE), DIMENSION(MAX_IGNITION_ZONES) :: IGNITION_ZONES      !< Coordinates of ignition zones used for detailed chemistry
+
+END MODULE CHEMCONS

Source/data.f90

@@ -677,6 +677,10 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
 OUTPUT_QUANTITY(144)%SHORT_NAME = 'Y_ELEM'
 OUTPUT_QUANTITY(144)%ELEM_ID_REQUIRED = .TRUE.
 
+OUTPUT_QUANTITY(145)%NAME='EQUIVALENCE RATIO'
+OUTPUT_QUANTITY(145)%UNITS=''
+OUTPUT_QUANTITY(145)%SHORT_NAME = 'equivRatio'
+
 OUTPUT_QUANTITY(150)%NAME  = 'SUM LUMPED VOLUME FRACTIONS'
 OUTPUT_QUANTITY(150)%UNITS  = ''
 OUTPUT_QUANTITY(150)%SHORT_NAME  = 'sum X'

Source/devc.f90

@@ -17,7 +17,7 @@ MODULE DEVICE_VARIABLES
                SPRAY_ANGLE(2,2),P0=0._EB,PX(3)=0._EB,PXX(3,3)=0._EB,VIEW_ANGLE,PROBE_DIAMETER
    INTEGER  :: PDPA_M=0,PDPA_N=0,N_SMOKEVIEW_PARAMETERS=0,N_SMOKEVIEW_IDS=0,N_INSERT,I_VEL=0,PARTICLES_PER_SECOND
    LOGICAL  :: PDPA_INTEGRATE=.TRUE.,PDPA_NORMALIZE=.TRUE.,HISTOGRAM_NORMALIZE=.TRUE.,HISTOGRAM=.FALSE., &
-               HISTOGRAM_CUMULATIVE=.FALSE.,SPARK=.FALSE.,TC=.TRUE.
+               HISTOGRAM_CUMULATIVE=.FALSE.,SPARK=.FALSE.,TC=.TRUE.,IGNITION_ZONE=.FALSE.
    REAL(EB) :: PDPA_START=0._EB,PDPA_END=1.E6_EB,PDPA_RADIUS=0.1_EB
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: TABLE_ROW, V_FACTOR
    INTEGER  :: PART_INDEX=-1,FLOW_RAMP_INDEX,SPRAY_PATTERN_INDEX,Z_INDEX=-999,Y_INDEX=-999,PRESSURE_RAMP_INDEX

Source/dump.f90

@@ -7316,7 +7316,7 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
 USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION,FED,FIC,GET_SPECIFIC_HEAT,RELATIVE_HUMIDITY, &
                               GET_CONDUCTIVITY,GET_MOLECULAR_WEIGHT,GET_MASS_FRACTION_ALL,GET_ENTHALPY,GET_SENSIBLE_ENTHALPY, &
                               GET_VISCOSITY,GET_POTENTIAL_TEMPERATURE,GET_SPECIFIC_GAS_CONSTANT,&
-                              SURFACE_DENSITY
+                              SURFACE_DENSITY, CALC_EQUIV_RATIO
 USE COMP_FUNCTIONS, ONLY : CURRENT_TIME,SYSTEM_MEM_USAGE
 USE RADCONS, ONLY: WL_LOW, WL_HIGH, RADTMP
 USE RAD, ONLY: BLACKBODY_FRACTION
@@ -7974,6 +7974,10 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z
          GAS_PHASE_OUTPUT_RES = GAS_PHASE_OUTPUT_RES + &
              ZZ_GET(NS)*SPECIES_MIXTURE(NS)%ATOMS(ELEM_INDX)*ELEMENT(ELEM_INDX)%MASS/SPECIES_MIXTURE(NS)%MW
       ENDDO
+   CASE(145)  ! EQUIVALENCE RATIO
+      ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES)
+      GAS_PHASE_OUTPUT_RES = 0._EB
+      CALL CALC_EQUIV_RATIO(ZZ_GET(1:N_TRACKED_SPECIES), GAS_PHASE_OUTPUT_RES)
    CASE(150) ! SUM LUMPED VOLUME FRACTIONS
       ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES)
       CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW)
@@ -11609,10 +11613,10 @@ SUBROUTINE DUMP_CVODE_SUBSTEPS()
      DO COLI = 1, NCOLS
         IF (COLI == NCOLS) THEN
            ! Writing the last column without a trailing comma
-           WRITE(LU_CVODE_SUBSTEPS, '(F18.5)') CVODE_SUBSTEP_DATA(ROWI, COLI)
+           WRITE(LU_CVODE_SUBSTEPS, '(ES12.5)') CVODE_SUBSTEP_DATA(ROWI, COLI)
         ELSE
            ! Writing columns with commas
-           WRITE(LU_CVODE_SUBSTEPS, '(F18.5, A)', ADVANCE="NO") CVODE_SUBSTEP_DATA(ROWI, COLI), ','
+           WRITE(LU_CVODE_SUBSTEPS, '(ES12.5, A)', ADVANCE="NO") CVODE_SUBSTEP_DATA(ROWI, COLI), ','
         ENDIF   
      ENDDO
   ENDDO