Merge pull request #14723 from ericvmueller/master

FDS Source: improve HRR curves for LEVEL_SET_MODE=4

Author: ericvmueller

Source/func.f90

@@ -4368,6 +4368,7 @@ SUBROUTINE PACK_BOUNDARY_PROP2(NM,IC,RC,OS,B2_INDEX,UNPACK_IT,COUNT_ONLY)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B2%Y_PLUS,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B2%Z_STAR,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B2%PHI_LS,UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B2%TAU_LS,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B2%WORK1,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B2%WORK2,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B2%WORK3,UNPACK_IT)

Source/type.f90

@@ -351,6 +351,7 @@ MODULE TYPES
    REAL(EB) :: Y_PLUS=1._EB          !< Dimensionless boundary layer thickness unit
    REAL(EB) :: Z_STAR=1._EB          !< Dimensionless boundary layer unit
    REAL(EB) :: PHI_LS=-1._EB         !< Level Set value for output only
+   REAL(EB) :: TAU_LS=-1._EB         !< Level Set HRR ramp characteristic time
    REAL(EB) :: WORK1=0._EB           !< Work array
    REAL(EB) :: WORK2=0._EB           !< Work array
    REAL(EB) :: WORK3=0._EB           !< Work array

Source/vege.f90

@@ -555,18 +555,18 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
                   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)
+                  IF (PHI_LS(IIG,JJG)>=0._EB .AND. B1%T_IGN>9.E5_EB) CALL IGNITE_GRID_CELL                  
                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)
+                     IF (PHI_LS(IIG,JJG)>=0._EB .AND. B1%T_IGN>9.E5_EB) CALL IGNITE_GRID_CELL                     
                   ENDDO
                ENDIF
             ENDDO
@@ -583,9 +583,9 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
             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)
+            IF (PHI_LS(IIG,JJG)>=0._EB .AND. B1%T_IGN>9.E5_EB) CALL IGNITE_GRID_CELL            
          ENDDO
       ENDDO
    ENDIF
@@ -600,10 +600,10 @@ SUBROUTINE LEVEL_SET_FIRESPREAD(T,DT,NM)
 
 SUBROUTINE IGNITE_GRID_CELL
 
-REAL(EB) :: CROSSING_DISTANCE
+REAL(EB) :: CROSSING_DISTANCE,FIRE_DEPTH,TAU_2
 
 B1%T_IGN = T
-ROS_MAG = MAX(0.01_EB,SQRT(SR_X_LS(IIG,JJG)**2 + SR_Y_LS(IIG,JJG)**2))  ! Rate Of Spread magnitude
+ROS_MAG = MAX(0.0001_EB,SQRT(SR_X_LS(IIG,JJG)**2 + SR_Y_LS(IIG,JJG)**2))  ! Rate Of Spread magnitude
 IF (ABS(SR_X_LS(IIG,JJG))<TWO_EPSILON_EB) THEN
    CROSSING_DISTANCE = DY(JJG)
 ELSEIF (ABS(SR_Y_LS(IIG,JJG))<TWO_EPSILON_EB) THEN
@@ -612,11 +612,23 @@ SUBROUTINE IGNITE_GRID_CELL
    CROSSING_DISTANCE = SQRT( (MIN(DX(IIG),ABS(SR_X_LS(IIG,JJG)/SR_Y_LS(IIG,JJG))*DY(JJG)))**2 + &
                              (MIN(DY(JJG),ABS(SR_Y_LS(IIG,JJG)/SR_X_LS(IIG,JJG))*DX(IIG)))**2 )
 ENDIF
-B1%BURN_DURATION = SF%BURN_DURATION + CROSSING_DISTANCE/ROS_MAG
-B1%AREA_ADJUST = SF%BURN_DURATION/B1%BURN_DURATION
-! Adjust mass flux for TAU ramp to conserve total energy output
+FIRE_DEPTH = SF%BURN_DURATION*ROS_MAG
+IF (FIRE_DEPTH >= CROSSING_DISTANCE) THEN
+   ! tau_1
+   B2%TAU_LS = CROSSING_DISTANCE/ROS_MAG
+   ! tau_2
+   TAU_2 = (FIRE_DEPTH - CROSSING_DISTANCE)/ROS_MAG
+ELSE
+   ! tau_1
+   B2%TAU_LS = FIRE_DEPTH/ROS_MAG
+   ! tau_2
+   TAU_2 = (CROSSING_DISTANCE - FIRE_DEPTH)/ROS_MAG
+ENDIF
+
+B1%BURN_DURATION = 2._EB*B2%TAU_LS + TAU_2
+! Adjust mass flux to conserve total energy output
 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)
+      SF%BURN_DURATION/(B2%TAU_LS+TAU_2)
 
 IF (STORE_LS_SPREAD_RATE) LS_SPREAD_RATE(OUTPUT_INDEX) = SQRT(SR_X_LS(IIG,JJG)**2 + SR_Y_LS(IIG,JJG)**2)
 

Source/wall.f90

@@ -1322,6 +1322,15 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)
                      B1%M_DOT_G_PP_ACTUAL(N) = 0._EB
                   ENDIF
                ENDIF
+            ! Trapezoidal HRR curve for level set spread
+            ELSEIF (SF%VEG_LSET_SPREAD) THEN
+               IF (TSI <= B2%TAU_LS) THEN
+                  B1%M_DOT_G_PP_ACTUAL(N) = TSI/B2%TAU_LS*SF%MASS_FLUX(N)
+               ELSEIF (TSI <= B1%BURN_DURATION-B2%TAU_LS) THEN
+                  B1%M_DOT_G_PP_ACTUAL(N) = SF%MASS_FLUX(N)
+               ELSE
+                  B1%M_DOT_G_PP_ACTUAL(N) = (1-(TSI-B1%BURN_DURATION+B2%TAU_LS)/B2%TAU_LS)*SF%MASS_FLUX(N)
+               ENDIF
             ELSE
                B1%M_DOT_G_PP_ACTUAL(N) = EVALUATE_RAMP(TSI,SF%RAMP(N)%INDEX,TAU=SF%RAMP(N)%TAU)*SF%MASS_FLUX(N)
             ENDIF