FDS Source: add LS SPREAD RATE output for storing local level set spread rates

Author: ericvmueller

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -11320,6 +11320,8 @@ \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.
 
+For level set simulations another output is availble, called \ct{LS SPREAD RATE}. This output stores the magnitude of the local spread rate calculated as the fire reaches a given point along the surface. The dependence on slope and local wind is related to the mode of level set which is activated (Sec.~\ref{info:level_set}). 
+
 
 \newpage
 
@@ -11552,6 +11554,7 @@ \section{Solid Phase Output Quantities}
 \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{LS SPREAD RATE}                            & Section \ref{info:fire_spread_output}         & \si{m/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

Source/cons.f90

@@ -274,6 +274,7 @@ MODULE GLOBAL_CONSTANTS
 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 :: STORE_LS_SPREAD_RATE=.FALSE.             !< Flag for outputting local level set spread rate magnitude
 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

Source/data.f90

@@ -1642,6 +1642,11 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
 OUTPUT_QUANTITY(-91)%SHORT_NAME = 't_r'
 OUTPUT_QUANTITY(-91)%PART_APPROPRIATE = .FALSE.
 
+OUTPUT_QUANTITY(-92)%NAME = 'LS SPREAD RATE'
+OUTPUT_QUANTITY(-92)%UNITS = 'm/s'
+OUTPUT_QUANTITY(-92)%SHORT_NAME = 'r'
+OUTPUT_QUANTITY(-92)%PART_APPROPRIATE = .FALSE.
+
 ! Condensation
 OUTPUT_QUANTITY(-100)%NAME = 'CONDENSATION HEAT FLUX'
 OUTPUT_QUANTITY(-100)%UNITS= 'kW/m2'

Source/dump.f90

@@ -3645,6 +3645,7 @@ SUBROUTINE DUMP_RESTART(T,DT,NM)
 
 IF (STORE_FIRE_ARRIVAL) WRITE(LU_CORE(NM)) FIRE_ARRIVAL_TIME
 IF (STORE_FIRE_RESIDENCE) WRITE(LU_CORE(NM)) FIRE_RESIDENCE_TIME
+IF (STORE_LS_SPREAD_RATE) WRITE(LU_CORE(NM)) LS_SPREAD_RATE
 
 CLOSE(LU_CORE(NM))
 
@@ -3871,6 +3872,7 @@ SUBROUTINE READ_RESTART(T,DT,NM)
 
 IF (STORE_FIRE_ARRIVAL) READ(LU_RESTART(NM)) FIRE_ARRIVAL_TIME
 IF (STORE_FIRE_RESIDENCE) READ(LU_RESTART(NM)) FIRE_RESIDENCE_TIME
+IF (STORE_LS_SPREAD_RATE) READ(LU_RESTART(NM)) LS_SPREAD_RATE
 
 CLOSE(LU_RESTART(NM))
 
@@ -9387,20 +9389,17 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_IN
    CASE(82) ! BLOWING CORRECTION
       SOLID_PHASE_OUTPUT = 0._EB
       IF (SF%INCLUDE_BOUNDARY_PROP2_TYPE) SOLID_PHASE_OUTPUT = B2%BLOWING_CORRECTION
-   CASE(90) ! FIRE ARRIVAL TIME
+   CASE(90:92) ! FIRE ARRIVAL TIME, FIRE RESIDENCE TIME, LS SPREAD RATE
       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)
+      SELECT CASE(INDX)
+         CASE(90); SOLID_PHASE_OUTPUT = FIRE_ARRIVAL_TIME(OUTPUT_INDEX)
+         CASE(91); SOLID_PHASE_OUTPUT = FIRE_RESIDENCE_TIME(OUTPUT_INDEX)
+         CASE(92); SOLID_PHASE_OUTPUT = LS_SPREAD_RATE(OUTPUT_INDEX)
+      END SELECT
 
    CASE(100) ! CONDENSATION HEAT FLUX
       SOLID_PHASE_OUTPUT = B1%Q_CONDENSE * 0.001_EB

Source/init.f90

@@ -1261,7 +1261,10 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM)
    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
-
+IF (STORE_LS_SPREAD_RATE) THEN
+   ALLOCATE(M%LS_SPREAD_RATE(1:M%N_INTERNAL_WALL_CELLS+M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_CFACE_CELLS))
+   CALL ChkMemErr('INIT','LS_SPREAD_RATE',IZERO) ; M%LS_SPREAD_RATE = 0._EB
+ENDIF
 
 CONTAINS
 

Source/mesh.f90

@@ -123,7 +123,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(EB), ALLOCATABLE, DIMENSION(:) :: FIRE_ARRIVAL_TIME,FIRE_RESIDENCE_TIME,LS_SPREAD_RATE
    REAL(FB), ALLOCATABLE, DIMENSION(:,:,:,:) :: QQ, QQ2
    REAL(FB), ALLOCATABLE, DIMENSION(:,:) :: PP,PPN,BNDF_TIME_INTEGRAL
    INTEGER, ALLOCATABLE, DIMENSION(:,:) :: IBK
@@ -382,7 +382,7 @@ 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(EB), POINTER, DIMENSION(:) :: FIRE_ARRIVAL_TIME,FIRE_RESIDENCE_TIME,LS_SPREAD_RATE
 
 REAL(FB), POINTER, DIMENSION(:,:,:,:) :: QQ, QQ2
 REAL(FB), POINTER, DIMENSION(:,:) :: PP,PPN,BNDF_TIME_INTEGRAL
@@ -644,6 +644,7 @@ SUBROUTINE POINT_TO_MESH(NM)
 FACE_WORK3=>M%FACE_WORK3
 FIRE_ARRIVAL_TIME=>M%FIRE_ARRIVAL_TIME
 FIRE_RESIDENCE_TIME=>M%FIRE_RESIDENCE_TIME
+LS_SPREAD_RATE=>M%LS_SPREAD_RATE
 QQ=>M%QQ
 QQ2=>M%QQ2
 PP=>M%PP

Source/read.f90

@@ -15130,6 +15130,9 @@ SUBROUTINE PROC_DEVC
       CASE ('FIRE RESIDENCE TIME')
          STORE_FIRE_RESIDENCE = .TRUE.
 
+      CASE ('LS SPREAD RATE')
+         STORE_LS_SPREAD_RATE = .TRUE.
+
    END SELECT SPECIAL_QUANTITIES
 
    IF ((DV%SPATIAL_STATISTIC/='null' .OR. DV%TEMPORAL_STATISTIC/='null') .AND. .NOT.DV%UNITS_SPECIFIED) THEN
@@ -15876,6 +15879,7 @@ SUBROUTINE READ_BNDF
    ! 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.
+   IF (TRIM(QUANTITY)=='LS SPREAD RATE') STORE_LS_SPREAD_RATE = .TRUE.
 
 ENDDO READ_BNDF_LOOP
 REWIND(LU_INPUT) ; INPUT_FILE_LINE_NUMBER = 0

Source/vege.f90

@@ -330,7 +330,7 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
 USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
 INTEGER, INTENT(IN) :: NM
 REAL(EB), INTENT(IN) :: T,DT
-INTEGER :: IIG,IW,JJG,IC
+INTEGER :: IIG,IW,JJG,IC,OUTPUT_INDEX
 INTEGER :: KDUM,KWIND,ICF,IKT
 REAL(EB) :: UMF_TMP,PHX,PHY,MAG_PHI,PHI_W_X,PHI_W_Y,UMF_X,UMF_Y,UMAG,ROS_MAG,UMF_MAG,ROTH_FACTOR,&
             SIN_THETA,COS_THETA,THETA
@@ -520,13 +520,16 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
                   IW = CELL(IC)%WALL_INDEX(-3)
                   WC => WALL(IW)
                   B1 => BOUNDARY_PROP1(WC%B1_INDEX)
+                  OUTPUT_INDEX = IW
                   IF (PHI_LS(IIG,JJG)>=0._EB .AND. B1%T_IGN>9.E5_EB) CALL IGNITE_GRID_CELL
                   B2 => BOUNDARY_PROP2(WC%B2_INDEX)
                   B2%PHI_LS = PHI_LS(IIG,JJG)
                ELSE  
                   DO IW=1,CUT_FACE(ICF)%NFACE ! All CC_INBOUNDARY CFACES on this cell.
                      CFA => CFACE(CUT_FACE(ICF)%CFACE_INDEX(IW))
                      B1  => BOUNDARY_PROP1(CFA%B1_INDEX)
+                     OUTPUT_INDEX = CUT_FACE(ICF)%CFACE_INDEX(IW) - &
+                        INTERNAL_CFACE_CELLS_LB+N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
                      IF (PHI_LS(IIG,JJG)>=0._EB .AND. B1%T_IGN>9.E5_EB) CALL IGNITE_GRID_CELL
                      B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
                      B2%PHI_LS = PHI_LS(IIG,JJG)
@@ -545,6 +548,7 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
             IW = CELL(IC)%WALL_INDEX(-3)
             WC => WALL(IW)
             B1 => BOUNDARY_PROP1(WC%B1_INDEX)
+            OUTPUT_INDEX = IW
             IF (PHI_LS(IIG,JJG)>=0._EB .AND. B1%T_IGN>9.E5_EB) CALL IGNITE_GRID_CELL
             B2 => BOUNDARY_PROP2(WC%B2_INDEX)
             B2%PHI_LS = PHI_LS(IIG,JJG)
@@ -580,6 +584,8 @@ SUBROUTINE IGNITE_GRID_CELL
 IF (LEVEL_SET_COUPLED_FIRE) B1%AREA_ADJUST =  B1%AREA_ADJUST * &
       B1%BURN_DURATION/(B1%BURN_DURATION - TWTH*SF%RAMP(REACTION(1)%FUEL_SMIX_INDEX)%TAU)
 
+IF (STORE_LS_SPREAD_RATE) LS_SPREAD_RATE(OUTPUT_INDEX) = SQRT(SR_X_LS(IIG,JJG)**2 + SR_Y_LS(IIG,JJG)**2)
+
 END SUBROUTINE IGNITE_GRID_CELL
 
 END SUBROUTINE LEVEL_SET_FIRESPREAD

Verification/WUI/LS_ellipse_0ms_00deg.fds

@@ -23,6 +23,7 @@
 &SLCF AGL_SLICE=0.1, QUANTITY='LEVEL SET VALUE' /
 &BNDF QUANTITY='FIRE ARRIVAL TIME'/
 &BNDF QUANTITY='FIRE RESIDENCE TIME'/
+&BNDF QUANTITY='LS SPREAD RATE'/
 
 &DEVC XYZ=0,4,0, ID='000', QUANTITY='FIRE ARRIVAL TIME', TEMPORAL_STATISTIC='INSTANT VALUE', IOR=3/
 &DEVC XYZ=2.8284,2.8284,0, ID='045', QUANTITY='FIRE ARRIVAL TIME', TEMPORAL_STATISTIC='INSTANT VALUE', IOR=3/

Verification/WUI/LS_ellipse_0ms_30deg.fds

@@ -26,6 +26,7 @@
 &SLCF AGL_SLICE=0.1, QUANTITY='LEVEL SET VALUE' /
 &BNDF QUANTITY='FIRE ARRIVAL TIME'/
 &BNDF QUANTITY='FIRE RESIDENCE TIME'/
+&BNDF QUANTITY='LS SPREAD RATE'/
 
 &DEVC XYZ=0,4,6, ID='000', QUANTITY='FIRE ARRIVAL TIME', TEMPORAL_STATISTIC='INSTANT VALUE', IOR=3/
 &DEVC XYZ=2.8284,2.8284,7, ID='045', QUANTITY='FIRE ARRIVAL TIME', TEMPORAL_STATISTIC='INSTANT VALUE', IOR=3/