Merge pull request #13937 from rmcdermo/master

FDS Source: add TEST_NEW_CHAR_MODEL

Source/cons.f90

@@ -165,6 +165,7 @@ MODULE GLOBAL_CONSTANTS
 INTEGER :: NO2_INDEX=0                     !< Index for NO2
 INTEGER :: ZETA_INDEX=0                    !< Index for unmixed fuel fraction, ZETA
 INTEGER :: MOISTURE_INDEX=0                !< Index for MATL MOISTURE
+INTEGER :: CHAR_INDEX=0                    !< Index for MATL CHAR
 
 INTEGER :: STOP_STATUS=NO_STOP             !< Indicator of whether and why to stop the job
 INTEGER :: INPUT_FILE_LINE_NUMBER=0        !< Indicator of what line in the input file is being read
@@ -270,9 +271,10 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: TENSOR_DIFFUSIVITY=.FALSE.               !< If true, use experimental tensor diffusivity model for spec and tmp
 LOGICAL :: OXPYRO_MODEL=.FALSE.                     !< Flag to use oxidative pyrolysis mass transfer model
 LOGICAL :: OUTPUT_WALL_QUANTITIES=.FALSE.           !< Flag to force call to WALL_MODEL
-LOGICAL :: FLUX_LIMITER_MW_CORRECTION=.FALSE.
+LOGICAL :: FLUX_LIMITER_MW_CORRECTION=.FALSE.       !< Flag for MW correction ensure consistent equation of state at face
 LOGICAL :: STORE_FIRE_ARRIVAL=.FALSE.               !< Flag for tracking arrival of spreading fire front over a surface
 LOGICAL :: STORE_FIRE_RESIDENCE=.FALSE.             !< Flag for tracking residence time of spreading fire front over a surface
+LOGICAL :: TEST_NEW_CHAR_MODEL=.FALSE.              !< Flag to envoke new char model
 
 INTEGER, ALLOCATABLE, DIMENSION(:) :: CHANGE_TIME_STEP_INDEX      !< Flag to indicate if a mesh needs to change time step
 INTEGER, ALLOCATABLE, DIMENSION(:) :: SETUP_PRESSURE_ZONES_INDEX  !< Flag to indicate if a mesh needs to keep searching for ZONEs

Source/dump.f90

@@ -8524,7 +8524,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_IN
 INTEGER, INTENT(IN), OPTIONAL :: OPT_WALL_INDEX,OPT_LP_INDEX,OPT_CFACE_INDEX,OPT_BNDF_INDEX,OPT_DEVC_INDEX,OPT_CUT_FACE_INDEX,&
                                  OPT_NODE_INDEX,OPT_PROF_INDEX
 INTEGER, INTENT(IN) :: INDX,Y_INDEX,Z_INDEX,PART_INDEX
-REAL(EB) :: Q_CON,RHOSUM,VOLSUM,ZZ_GET(1:N_TRACKED_SPECIES),Y_SPECIES,DEPTH,UN,H_S,RHO_D_DYDN,U_CELL,V_CELL,W_CELL,&
+REAL(EB) :: Q_CON,RHOSUM,VOLSUM,ZZ_GET(1:N_TRACKED_SPECIES),Y_SPECIES,DEPTH,CHAR_FRONT,UN,H_S,RHO_D_DYDN,U_CELL,V_CELL,W_CELL,&
             LTMP,ATMP,CTMP,H_W_EFF,X0,VOL,DN,PRESS,&
             NVEC(3),PVEC(3),TAU_IJ(3,3),VEL_CELL(3),VEL_WALL(3),MU_WALL,RHO_WALL,FVEC(3),SVEC(3),TVEC1(3),TVEC2(3),&
             PR1,PR2,Z1,Z2,RADIUS,CUT_FACE_AREA,SOLID_PHASE_OUTPUT_CTF,AAA,BBB,CCC,ALP,BET,GAM,MMM,DTMP
@@ -9277,7 +9277,12 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_IN
       IF (SF%INCLUDE_BOUNDARY_PROP2_TYPE .AND. MATL_INDEX>0) THEN
          ML => MATERIAL(MATL_INDEX)
          ! for the moment this assumes there is only one char reaction
-         IF (ML%N_O2(1)>0._EB) SOLID_PHASE_OUTPUT = B2%Y_O2_F*EXP(-ONE_D%X(I_DEPTH-1)/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(1)))
+         IF (ML%N_O2(1)>0._EB) THEN
+            DEPTH = 0.5_EB*(ONE_D%X(I_DEPTH-1)+ONE_D%X(I_DEPTH))
+            CHAR_FRONT = 0._EB
+            IF (TEST_NEW_CHAR_MODEL) CHAR_FRONT = ONE_D%X(B2%I_CHAR_FRONT-1)
+            SOLID_PHASE_OUTPUT = B2%Y_O2_F*EXP(-MAX(0._EB,DEPTH-CHAR_FRONT)/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(1)))
+         ENDIF
       ENDIF
    CASE(90) ! FIRE ARRIVAL TIME
       IF (PRESENT(OPT_WALL_INDEX)) THEN

Source/read.f90

@@ -1732,7 +1732,8 @@ SUBROUTINE READ_MISC
                 CFL_MAX,CFL_MIN,CFL_VELOCITY_NORM,CHECK_FO,CHECK_HT,CHECK_VN, &
                 CNF_CUTOFF,CONSTANT_SPECIFIC_HEAT_RATIO,&
                 C_SMAGORINSKY,C_VREMAN,C_WALE,DEPOSITION,EXTERNAL_FILENAME,&
-                FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,FREEZE_VELOCITY,FYI,GAMMA,GRAVITATIONAL_DEPOSITION,&
+                FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,FLUX_LIMITER_MW_CORRECTION,&
+                FREEZE_VELOCITY,FYI,GAMMA,GRAVITATIONAL_DEPOSITION,&
                 GRAVITATIONAL_SETTLING,GVEC,H_F_REFERENCE_TEMPERATURE,&
                 HUMIDITY,HVAC_LOCAL_PRESSURE,HVAC_MASS_TRANSPORT_CELL_L,HVAC_PRES_RELAX,HVAC_QFAN,IBLANK_SMV,I_MAX_TEMP,&
                 LES_FILTER_TYPE,LEVEL_SET_ELLIPSE,LEVEL_SET_MODE,&
@@ -1745,7 +1746,7 @@ SUBROUTINE READ_MISC
                 RAMP_UX,RAMP_UY,RAMP_UZ,RAMP_VX,RAMP_VY,RAMP_VZ,RAMP_WX,RAMP_WY,RAMP_WZ,&
                 RESTART,RESTART_CHID,SC,&
                 RND_SEED,SIMULATION_MODE,SMOKE3D_16,SMOKE_ALBEDO,SOLID_PHASE_ONLY,SOOT_DENSITY,SOOT_OXIDATION,&
-                TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,FLUX_LIMITER_MW_CORRECTION,TEXTURE_ORIGIN,&
+                TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,TEST_NEW_CHAR_MODEL,TEXTURE_ORIGIN,&
                 THERMOPHORETIC_DEPOSITION,THERMOPHORETIC_SETTLING,THICKEN_OBSTRUCTIONS,&
                 TMPA,TURBULENCE_MODEL,TURBULENT_DEPOSITION,UVW_FILE,&
                 VERBOSE,VISIBILITY_FACTOR,VN_MAX,VN_MIN,Y_CO2_INFTY,Y_O2_INFTY,&
@@ -7431,6 +7432,9 @@ SUBROUTINE PROC_MATL
    H_ADJUST = ML%REFERENCE_ENTHALPY - ANS
    ML%H = ML%H + H_ADJUST
 
+   ! Store index of char material
+   IF (TRIM(ML%ID)=='CHAR') CHAR_INDEX=N
+
 ENDDO PROC_MATL_LOOP
 
 ! Construct and solve linear system to adjust MATL enthalpies to conserve energy. For reaction i for MATL m the equation is
@@ -7674,7 +7678,6 @@ SUBROUTINE PROC_MATL
    ENDIF
 ENDDO
 
-
 END SUBROUTINE PROC_MATL
 
 

Source/type.f90

@@ -362,6 +362,7 @@ MODULE TYPES
    INTEGER :: SURF_INDEX=-1          !< Surface index
    INTEGER :: HEAT_TRANSFER_REGIME=0 !< 1=Forced convection, 2=Natural convection, 3=Impact convection, 4=Resolved
    INTEGER :: Y_O2_ITER=0            !< Number of iterations for surface O2 solve
+   INTEGER :: I_CHAR_FRONT=1         !< Index of solid cell where char formation starts
 
 END TYPE BOUNDARY_PROP2_TYPE
 

Source/wall.f90

@@ -2883,12 +2883,12 @@ SUBROUTINE PERFORM_PYROLYSIS
 USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION,GET_SPECIFIC_HEAT
 REAL(EB), DIMENSION(N_TRACKED_SPECIES) :: M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,ZZ_GET
 REAL(EB), DIMENSION(MAX_MATERIALS) :: RHO_TEMP,M_DOT_S_PPP
-REAL(EB) :: Q_DOT_G_PPP,Q_DOT_O2_PPP,CP_FILM,TMP_FILM,H_MASS,&
+REAL(EB) :: Q_DOT_G_PPP,Q_DOT_O2_PPP,CP_FILM,TMP_FILM,H_MASS,CHAR_FRONT,&
             Y_O2_G,Y_O2_F,M_DOT_O2_PP,Y_LOWER,Y_UPPER,M_DOT_ERROR,M_DOT_ERROR_OLD,Y_O2_F_OLD,DY,DE
 REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART,M_DOT_PART
 INTEGER :: ITER,MAX_ITER
 LOGICAL :: REMOVE_LAYER
-REAL(EB), PARAMETER :: M_DOT_ERROR_TOL=1.E-6_EB
+REAL(EB), PARAMETER :: M_DOT_ERROR_TOL=1.E-6_EB, CHAR_DENSITY_THRESHOLD=5.7_EB ! approx. ASH_DENSITY
 
 ! Get surface oxygen mass fraction
 
@@ -2899,6 +2899,23 @@ SUBROUTINE PERFORM_PYROLYSIS
    Y_O2_F = 0._EB
 ENDIF
 
+! Determine char front position
+
+CHAR_FRONT=0._EB
+IF (TEST_NEW_CHAR_MODEL .AND. Y_O2_F>TWO_EPSILON_EB .AND. CHAR_INDEX>0) THEN
+   ! The new char model starts the exp profile of O2 at the CHAR_FRONT
+   ! Find first cell where char has not been consumed
+   CHAR_FRONT_POINTS_LOOP: DO I=B2%I_CHAR_FRONT,NWP
+      IF (ONE_D%MATL_COMP(CHAR_INDEX)%RHO(I)>CHAR_DENSITY_THRESHOLD) THEN
+         CHAR_FRONT=ONE_D%X(I-1)
+         B2%I_CHAR_FRONT=I ! store last position to save time on next time step
+         EXIT CHAR_FRONT_POINTS_LOOP
+      ENDIF
+   ENDDO CHAR_FRONT_POINTS_LOOP
+ENDIF
+
+! Initialize iterations for OXPYRO_MODEL
+
 MAX_ITER=1
 
 ! Initialize parameters for oxidative pyrolysis mass transfer model
@@ -2947,13 +2964,13 @@ SUBROUTINE PERFORM_PYROLYSIS
 
       IF (PRESENT(PARTICLE_INDEX)) THEN
          CALL PYROLYSIS(ONE_D%N_MATL,ONE_D%MATL_INDEX,SURF_INDEX,BC%IIG,BC%JJG,BC%KKG,ONE_D%TMP(I),B1%TMP_F,Y_O2_F,BC%IOR,&
-                        RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),ONE_D%X(I-1),DX_S,DT_BC_SUB,&
+                        RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),ONE_D%X(I-1),CHAR_FRONT,DX_S,DT_BC_SUB,&
                         M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_S(I),Q_DOT_G_PPP,Q_DOT_O2_PPP,&
                         Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B1%B_NUMBER,LAYER_INDEX(I),REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX=I,&
                         R_DROP=R_SURF,LPU=U_SURF,LPV=V_SURF,LPW=W_SURF)
       ELSE
          CALL PYROLYSIS(ONE_D%N_MATL,ONE_D%MATL_INDEX,SURF_INDEX,BC%IIG,BC%JJG,BC%KKG,ONE_D%TMP(I),B1%TMP_F,Y_O2_F,BC%IOR,&
-                        RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),ONE_D%X(I-1),DX_S,DT_BC_SUB,&
+                        RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),0.5*(ONE_D%X(I-1)+ONE_D%X(I)),CHAR_FRONT,DX_S,DT_BC_SUB,&
                         M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_S(I),Q_DOT_G_PPP,Q_DOT_O2_PPP,&
                         Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B1%B_NUMBER,LAYER_INDEX(I),REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX=I)
       ENDIF
@@ -3040,6 +3057,7 @@ END SUBROUTINE SOLID_HEAT_TRANSFER
 !> \param RHO_DOT_OUT (1:N_MATS) Array of component reaction rates (kg/m3/s)
 !> \param RHO_S (1:N_MATS) Array of component densities (kg/m3)
 !> \param DEPTH Distance from surface (m)
+!> \param CHAR_FRONT Distance from surface to start of char front (m)
 !> \param DX_S Array of node sizes (m)
 !> \param DT_BC Time step used by the solid phase solver (s)
 !> \param M_DOT_G_PPP_ADJUST (1:N_TRACKED_SPECIES) Adjusted mass generation rate per unit volume of the gas species
@@ -3062,7 +3080,8 @@ END SUBROUTINE SOLID_HEAT_TRANSFER
 !> \param LPV (OPTIONAL) y component of droplet velocity (m/s)
 !> \param LPW (OPTIONAL) z component of droplet velocity (m/s)
 
-SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F,IOR,RHO_DOT_OUT,RHO_S,DEPTH,DX_S,DT_BC,&
+SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F,IOR,&
+                     RHO_DOT_OUT,RHO_S,DEPTH,CHAR_FRONT,DX_S,DT_BC,&
                      M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_DOT_S_PPP,Q_DOT_G_PPP,Q_DOT_O2_PPP,&
                      Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B_NUMBER,LAYER_INDEX,REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX,&
                      R_DROP,LPU,LPV,LPW)
@@ -3075,7 +3094,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F
 INTEGER, INTENT(IN), OPTIONAL :: SOLID_CELL_INDEX
 LOGICAL, INTENT(IN) :: REMOVE_LAYER
 REAL(EB), INTENT(OUT), DIMENSION(:,:) :: RHO_DOT_OUT(N_MATS)
-REAL(EB), INTENT(IN) :: TMP_S,TMP_F,DT_BC,DEPTH,RHO_S(N_MATS),Y_O2_F
+REAL(EB), INTENT(IN) :: TMP_S,TMP_F,DT_BC,DEPTH,RHO_S(N_MATS),Y_O2_F,CHAR_FRONT
 REAL(EB), INTENT(IN), OPTIONAL :: R_DROP,LPU,LPV,LPW
 REAL(EB), INTENT(IN), DIMENSION(NWP_MAX) :: DX_S
 REAL(EB), DIMENSION(:) :: ZZ_GET(1:N_TRACKED_SPECIES),Y_ALL(1:N_SPECIES)
@@ -3358,7 +3377,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F
                ! Calculate oxygen volume fraction at the surface
                X_O2 = SPECIES(O2_INDEX)%RCON*Y_O2_F/RSUM(IIG,JJG,KKG)
                ! Calculate oxygen concentration inside the material, assuming decay function
-               X_O2 = X_O2 * EXP(-DEPTH/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(J)))
+               X_O2 = X_O2 * EXP(-MAX(0._EB,DEPTH-CHAR_FRONT)/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(J)))
                REACTION_RATE = REACTION_RATE * X_O2**ML%N_O2(J)
             ENDIF
             REACTION_RATE = MIN(REACTION_RATE,ML%MAX_REACTION_RATE(J))  ! User-specified limit