diff --git a/Manuals/FDS_User_Guide/FDS_User_Guide.tex b/Manuals/FDS_User_Guide/FDS_User_Guide.tex index 4470056f531..8d2de0cc10a 100644 --- a/Manuals/FDS_User_Guide/FDS_User_Guide.tex +++ b/Manuals/FDS_User_Guide/FDS_User_Guide.tex @@ -11222,6 +11222,10 @@ \subsection{Extinction} \label{info:extinct_out} In combustion, knowing if, when, or where chemical reactions have been extinguished is important. The output quantity \ct{EXTINCTION} tells the user whether or not combustion has been prevented by the extinction routine. By default, \ct{EXTINCTION} = 0, which means that the FDS extinction routine has not prevented combustion. An \ct{EXTINCTION} value of 1 means that the routine has prevented combustion. The criteria for an \ct{EXTINCTION} value of 1 is the presence of fuel and oxidizer without any energy release. An \ct{EXTINCTION} value of -1 means that there is either no fuel or no oxidizer present. +\subsection{Fire spread over a surface} +\label{info:fire_spread_output} +In some cases it can be useful to output timings related to the spread of fire over a surface. For example, when modeling wildland fires the shape and spread of the fire front can be very important. For this reason, two specialized output quantities are available as boundary files (Sec.~\ref{info:BNDF}) or as measurements from devices places on a solid boundary (Sec.~\ref{info:DEVC2}). These are \ct{FIRE ARRIVAL TIME} and \ct{FIRE RESIDENCE TIME}. The \ct{FIRE ARRIVAL TIME} quantity outputs the time at which the gas-phase cell adjacent to the solid exceeds a threshold for heat release rate per volume (\ct{HRRPUV}) and the \ct{FIRE RESIDENCE TIME} gives the cummulative time over which the threshold is exceeded during the simulation. The chosen threshold is the same as used by smokeview for rendering \ct{HRRPUV}, as described in Sec.~\ref{info:SMOKE3D}. In the case of a level set simulation (Sec.~\ref{info:level_set}), the \ct{FIRE ARRIVAL TIME} can be computed directly from the level set value and the \ct{FIRE RESIDENCE TIME} comes from the spread-rate adjusted burning duration of the fuel, as described in Sec.~\ref{level_set_fuel_model_1}. These two quantities are cummulatively populated over time such that a full picture of the fire spread can be obtained from relatively infrequent outputs - theoretically only one snapshot at the end of the simulation is required. + \newpage @@ -11449,6 +11453,8 @@ \section{Solid Phase Output Quantities} \ct{ENTHALPY FLUX WALL} & Section~\ref{info:enthalpy_flux} & \unit{kW/m^2} & B,D \\ \hline \ct{TOTAL HEAT FLUX} & Section~\ref{info:heat_flux} & \unit{kW/m^2} & B,D \\ \hline \ct{EMISSIVITY} & Surface emissivity (usually constant) & & B,D \\ \hline +\ct{FIRE ARRIVAL TIME} & Section \ref{info:fire_spread_output} & \si{s} & B,D \\ \hline +\ct{FIRE RESIDENCE TIME} & Section \ref{info:fire_spread_output} & \si{s} & B,D \\ \hline \ct{GAS DENSITY} & Gas Density near wall & \si{kg/m^3} & B,D \\ \hline \ct{GAS TEMPERATURE} & Gas Temperature near wall & $^\circ$C & B,D \\ \hline \ct{HEAT TRANSFER COEFFICIENT} & Section \ref{info:convection} & \si{W/(m^2.K)} & B,D \\ \hline diff --git a/Source/cons.f90 b/Source/cons.f90 index d84d038d1ce..53b0a48dd25 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -271,6 +271,8 @@ MODULE GLOBAL_CONSTANTS 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 :: 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 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/data.f90 b/Source/data.f90 index 2b0eeead39c..21a9325680a 100644 --- a/Source/data.f90 +++ b/Source/data.f90 @@ -1605,14 +1605,14 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES OUTPUT_QUANTITY(-80)%SHORT_NAME = 'hrrpua_O2' ! Fire spread -OUTPUT_QUANTITY(-90)%NAME = 'IGNITION TIME' +OUTPUT_QUANTITY(-90)%NAME = 'FIRE ARRIVAL TIME' OUTPUT_QUANTITY(-90)%UNITS = 's' -OUTPUT_QUANTITY(-90)%SHORT_NAME = 't_ig' +OUTPUT_QUANTITY(-90)%SHORT_NAME = 't_a' OUTPUT_QUANTITY(-90)%PART_APPROPRIATE = .FALSE. -OUTPUT_QUANTITY(-91)%NAME = 'BURN DURATION' +OUTPUT_QUANTITY(-91)%NAME = 'FIRE RESIDENCE TIME' OUTPUT_QUANTITY(-91)%UNITS = 's' -OUTPUT_QUANTITY(-91)%SHORT_NAME = 't_b' +OUTPUT_QUANTITY(-91)%SHORT_NAME = 't_r' OUTPUT_QUANTITY(-91)%PART_APPROPRIATE = .FALSE. ! Condensation diff --git a/Source/dump.f90 b/Source/dump.f90 index 0a080cce34f..d668a782686 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -62,6 +62,7 @@ MODULE DUMP SUBROUTINE UPDATE_GLOBAL_OUTPUTS(T,DT,NM) USE COMP_FUNCTIONS, ONLY : CURRENT_TIME +USE VEGE, ONLY : UPDATE_FIRE_SPREAD_OUTPUTS REAL(EB) :: TNOW INTEGER, INTENT(IN) :: NM REAL(EB),INTENT(IN) :: T,DT @@ -72,6 +73,8 @@ SUBROUTINE UPDATE_GLOBAL_OUTPUTS(T,DT,NM) CALL UPDATE_HRR(DT,NM) CALL UPDATE_MASS(DT,NM) +! Update fire spread outputs +IF (STORE_FIRE_ARRIVAL .OR. STORE_FIRE_RESIDENCE) CALL UPDATE_FIRE_SPREAD_OUTPUTS(T,DT,NM) CALL UPDATE_DEVICES_1(T,DT,NM) T_USED(7) = T_USED(7) + CURRENT_TIME() - TNOW @@ -3540,6 +3543,9 @@ SUBROUTINE DUMP_RESTART(T,DT,NM) WRITE(LU_CORE(NM)) PATCH ENDIF +IF (STORE_FIRE_ARRIVAL) WRITE(LU_CORE(NM)) FIRE_ARRIVAL_TIME +IF (STORE_FIRE_RESIDENCE) WRITE(LU_CORE(NM)) FIRE_RESIDENCE_TIME + CLOSE(LU_CORE(NM)) END SUBROUTINE DUMP_RESTART @@ -3763,6 +3769,9 @@ SUBROUTINE READ_RESTART(T,DT,NM) READ(LU_RESTART(NM)) PATCH ENDIF +IF (STORE_FIRE_ARRIVAL) READ(LU_RESTART(NM)) FIRE_ARRIVAL_TIME +IF (STORE_FIRE_RESIDENCE) READ(LU_RESTART(NM)) FIRE_RESIDENCE_TIME + CLOSE(LU_RESTART(NM)) ! Keep track of whether the output timing intervals are specified by the user or not @@ -8519,7 +8528,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_IN 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 -INTEGER :: I_DEPTH,II2,IIG,JJG,KKG,NN,IWX,SURF_INDEX,I,J,NWP,M_INDEX,ICC,IND1,IND2,IC2,ITMP,ICF,JCF,NFACE,NR,MATL_INDEX +INTEGER :: I_DEPTH,II2,IIG,JJG,KKG,NN,IWX,SURF_INDEX,I,J,NWP,M_INDEX,ICC,IND1,IND2,IC2,ITMP,ICF,JCF,NFACE,NR,MATL_INDEX,OUTPUT_INDEX TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1=>NULL() TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2=>NULL() TYPE(BOUNDARY_RADIA_TYPE), POINTER :: BR=>NULL() @@ -9263,10 +9272,20 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_IN IF (SF%INCLUDE_BOUNDARY_PROP2_TYPE) SOLID_PHASE_OUTPUT = B2%Y_O2_ITER CASE(80) ! OXIDATIVE HRRPUA SOLID_PHASE_OUTPUT = B1%Q_DOT_O2_PP*0.001_EB - CASE(90) ! IGNITION TIME - SOLID_PHASE_OUTPUT = B1%T_IGN - CASE(91) ! BURN DURATION - SOLID_PHASE_OUTPUT = B1%BURN_DURATION + CASE(90) ! FIRE ARRIVAL TIME + IF (PRESENT(OPT_WALL_INDEX)) THEN + OUTPUT_INDEX = OPT_WALL_INDEX + ELSEIF (PRESENT(OPT_CFACE_INDEX)) THEN + OUTPUT_INDEX = OPT_CFACE_INDEX-INTERNAL_CFACE_CELLS_LB+N_INTERNAL_WALL_CELLS+N_EXTERNAL_WALL_CELLS + ENDIF + SOLID_PHASE_OUTPUT = FIRE_ARRIVAL_TIME(OUTPUT_INDEX) + CASE(91) ! FIRE RESIDENCE TIME + IF (PRESENT(OPT_WALL_INDEX)) THEN + OUTPUT_INDEX = OPT_WALL_INDEX + ELSEIF (PRESENT(OPT_CFACE_INDEX)) THEN + OUTPUT_INDEX = OPT_CFACE_INDEX-INTERNAL_CFACE_CELLS_LB+N_INTERNAL_WALL_CELLS+N_EXTERNAL_WALL_CELLS + ENDIF + SOLID_PHASE_OUTPUT = FIRE_RESIDENCE_TIME(OUTPUT_INDEX) CASE(100) ! CONDENSATION HEAT FLUX SOLID_PHASE_OUTPUT = B1%Q_CONDENSE * 0.001_EB diff --git a/Source/init.f90 b/Source/init.f90 index bf416922314..bc26bd0ebd6 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -1246,6 +1246,16 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM) CALL INITIALIZE_INTERPOLATION +! If fire spread outputs requested allocate arrays +IF (STORE_FIRE_ARRIVAL) THEN + ALLOCATE(M%FIRE_ARRIVAL_TIME(1:M%N_INTERNAL_WALL_CELLS+M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_CFACE_CELLS)) + CALL ChkMemErr('INIT','FIRE_ARRIVAL_TIME',IZERO) ; M%FIRE_ARRIVAL_TIME = 1.E6_EB +ENDIF +IF (STORE_FIRE_RESIDENCE) THEN + ALLOCATE(M%FIRE_RESIDENCE_TIME(1:M%N_INTERNAL_WALL_CELLS+M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_CFACE_CELLS)) + CALL ChkMemErr('INIT','FIRE_RESIDENCE_TIME',IZERO) ; M%FIRE_RESIDENCE_TIME = 0._EB +ENDIF + CONTAINS diff --git a/Source/mesh.f90 b/Source/mesh.f90 index ad9388f7008..6a5f122cc14 100644 --- a/Source/mesh.f90 +++ b/Source/mesh.f90 @@ -120,6 +120,7 @@ MODULE MESH_VARIABLES REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: TURB_WORK9,TURB_WORK10 REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: CCVELDIV,CARTVELDIV REAL(EB), ALLOCATABLE, DIMENSION(:) :: WALL_WORK1,WALL_WORK2,FACE_WORK1,FACE_WORK2,FACE_WORK3 + REAL(EB), ALLOCATABLE, DIMENSION(:) :: FIRE_ARRIVAL_TIME,FIRE_RESIDENCE_TIME REAL(FB), ALLOCATABLE, DIMENSION(:,:,:,:) :: QQ, QQ2 REAL(FB), ALLOCATABLE, DIMENSION(:,:) :: PP,PPN,BNDF_TIME_INTEGRAL INTEGER, ALLOCATABLE, DIMENSION(:,:) :: IBK @@ -377,6 +378,8 @@ MODULE MESH_POINTERS REAL(EB), POINTER, DIMENSION(:,:,:) :: CCVELDIV,CARTVELDIV REAL(EB), POINTER, DIMENSION(:) :: WALL_WORK1,WALL_WORK2,FACE_WORK1,FACE_WORK2,FACE_WORK3 +REAL(EB), POINTER, DIMENSION(:) :: FIRE_ARRIVAL_TIME,FIRE_RESIDENCE_TIME + REAL(FB), POINTER, DIMENSION(:,:,:,:) :: QQ, QQ2 REAL(FB), POINTER, DIMENSION(:,:) :: PP,PPN,BNDF_TIME_INTEGRAL INTEGER, POINTER, DIMENSION(:,:) :: IBK @@ -631,6 +634,8 @@ SUBROUTINE POINT_TO_MESH(NM) FACE_WORK1=>M%FACE_WORK1 FACE_WORK2=>M%FACE_WORK2 FACE_WORK3=>M%FACE_WORK3 +FIRE_ARRIVAL_TIME=>M%FIRE_ARRIVAL_TIME +FIRE_RESIDENCE_TIME=>M%FIRE_RESIDENCE_TIME QQ=>M%QQ QQ2=>M%QQ2 PP=>M%PP diff --git a/Source/read.f90 b/Source/read.f90 index 48b5a07de20..729972deb4b 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -15039,6 +15039,12 @@ SUBROUTINE PROC_DEVC CALL SHUTDOWN(MESSAGE) ; RETURN ENDIF + CASE ('FIRE ARRIVAL TIME') + STORE_FIRE_ARRIVAL = .TRUE. + + CASE ('FIRE RESIDENCE TIME') + STORE_FIRE_RESIDENCE = .TRUE. + END SELECT SPECIAL_QUANTITIES IF ((DV%SPATIAL_STATISTIC/='null' .OR. DV%TEMPORAL_STATISTIC/='null') .AND. .NOT.DV%UNITS_SPECIFIED) THEN @@ -15782,6 +15788,10 @@ SUBROUTINE READ_BNDF CALL CHANGE_UNITS(QUANTITY,BF%UNITS,DUMMY,TEMPORAL_STATISTIC,LU_ERR) ENDIF + ! Set flags for fire spread outputs + IF (TRIM(QUANTITY)=='FIRE ARRIVAL TIME') STORE_FIRE_ARRIVAL = .TRUE. + IF (TRIM(QUANTITY)=='FIRE RESIDENCE TIME') STORE_FIRE_RESIDENCE = .TRUE. + ENDDO READ_BNDF_LOOP REWIND(LU_INPUT) ; INPUT_FILE_LINE_NUMBER = 0 diff --git a/Source/vege.f90 b/Source/vege.f90 index 1f162339261..4f586fbc9d6 100644 --- a/Source/vege.f90 +++ b/Source/vege.f90 @@ -9,7 +9,7 @@ MODULE VEGE USE MEMORY_FUNCTIONS, ONLY: CHKMEMERR IMPLICIT NONE (TYPE,EXTERNAL) PRIVATE -PUBLIC INITIALIZE_LEVEL_SET_FIRESPREAD_1,INITIALIZE_LEVEL_SET_FIRESPREAD_2,LEVEL_SET_FIRESPREAD +PUBLIC INITIALIZE_LEVEL_SET_FIRESPREAD_1,INITIALIZE_LEVEL_SET_FIRESPREAD_2,LEVEL_SET_FIRESPREAD,UPDATE_FIRE_SPREAD_OUTPUTS INTEGER :: IZERO INTEGER :: LIMITER_LS REAL(EB) :: B_ROTH,BETA_OP_ROTH,C_ROTH,E_ROTH,T_NOW @@ -1279,4 +1279,87 @@ REAL(EB) FUNCTION ROS_NO_WIND_NO_SLOPE(ROTHERMEL_FUEL_INDEX,SURF_INDEX) END FUNCTION ROS_NO_WIND_NO_SLOPE +!> \brief Update quantities which are used to output metrics for quantifying fire spread +!> +!> \param T Current time (s) +!> \param DT Time step (s) +!> \param NM Mesh number +!> \details: Cycles through each solid wall cell and CFACE and populates output arrays as needed + +SUBROUTINE UPDATE_FIRE_SPREAD_OUTPUTS(T,DT,NM) + +INTEGER, INTENT(IN) :: NM +REAL(EB), INTENT(IN) :: T,DT +INTEGER :: ICF,IW,OUTPUT_INDEX +TYPE(WALL_TYPE), POINTER :: WC +TYPE(CFACE_TYPE), POINTER :: CFA +TYPE(SURFACE_TYPE), POINTER :: SF +TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC + +CALL POINT_TO_MESH(NM) + +DO IW = 1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS + WC => WALL(IW) + IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE + BC => BOUNDARY_COORD(WC%BC_INDEX) + B1 => BOUNDARY_PROP1(WC%B1_INDEX) + SF => SURFACE(WC%SURF_INDEX) + OUTPUT_INDEX = IW + CALL COMPUTE_FIRE_SPREAD_OUTPUTS +ENDDO + +DO ICF = INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS + CFA => CFACE(ICF) + IF (CFA%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE + BC => BOUNDARY_COORD(CFA%BC_INDEX) + B1 => BOUNDARY_PROP1(CFA%B1_INDEX) + SF => SURFACE(CFA%SURF_INDEX) + OUTPUT_INDEX = ICF-INTERNAL_CFACE_CELLS_LB+N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS + CALL COMPUTE_FIRE_SPREAD_OUTPUTS +ENDDO + +CONTAINS + +!> \brief Compute output metrics for quantifying fire spread +!> +!> For level set, fire arrival and residence time is already stored in B1%T_IGN and B1%BURN_DURATION. +!> In other cases, fire presence is detected in the cell adjacent to the solid face by evaluating +!> HRRPUV against the default smokeview threshold. + +SUBROUTINE COMPUTE_FIRE_SPREAD_OUTPUTS + +REAL(EB) :: HRRPUVCUT +LOGICAl :: FIRE_PRESENT + +FIRE_PRESENT = .FALSE. +HRRPUVCUT = MIN(200._EB,20._EB/CHARACTERISTIC_CELL_SIZE) +IF (Q(BC%IIG,BC%JJG,BC%KKG)>HRRPUVCUT) FIRE_PRESENT = .TRUE. + +IF (STORE_FIRE_ARRIVAL) THEN + IF (SF%VEG_LSET_SPREAD .AND. B1%T_IGN<9.E5_EB) THEN + IF (FIRE_ARRIVAL_TIME(OUTPUT_INDEX)>9.E5_EB) & + FIRE_ARRIVAL_TIME(OUTPUT_INDEX) = B1%T_IGN + ELSE + IF (FIRE_PRESENT .AND. FIRE_ARRIVAL_TIME(OUTPUT_INDEX)>9.E5_EB) FIRE_ARRIVAL_TIME(OUTPUT_INDEX) = T + ! reset if threshold is not met for 1/100th of a second + IF (.NOT. FIRE_PRESENT .AND. ABS(T-FIRE_ARRIVAL_TIME(OUTPUT_INDEX)) < 0.01_EB) FIRE_ARRIVAL_TIME(OUTPUT_INDEX) = 1.E6_EB + ENDIF +ENDIF + +IF (STORE_FIRE_RESIDENCE) THEN + IF (SF%VEG_LSET_SPREAD .AND. B1%T_IGN<9.E5_EB) THEN + IF (FIRE_RESIDENCE_TIME(OUTPUT_INDEX)