diff --git a/Source/cons.f90 b/Source/cons.f90 index 62254cecee..79a0343519 100644 --- a/Source/cons.f90 +++ b/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 diff --git a/Source/dump.f90 b/Source/dump.f90 index 0fbf9d55a2..b9e88e2db8 100644 --- a/Source/dump.f90 +++ b/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 diff --git a/Source/read.f90 b/Source/read.f90 index 486c1d3374..d853e80aca 100644 --- a/Source/read.f90 +++ b/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 diff --git a/Source/type.f90 b/Source/type.f90 index 0f30e0d185..cd71c47d48 100644 --- a/Source/type.f90 +++ b/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 diff --git a/Source/wall.f90 b/Source/wall.f90 index fb5493d240..3e50b23181 100644 --- a/Source/wall.f90 +++ b/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