Merge pull request #13167 from rmcdermo/master

FDS Source: add optional CHECK_FO flag to enforce tighter time step c…

Source/cons.f90

@@ -210,6 +210,7 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: BNDF_DEFAULT=.TRUE.              !< Output boundary output files
 LOGICAL :: SPATIAL_GRAVITY_VARIATION=.FALSE.!< Assume gravity varies as a function of the \f$ x \f$ coordinate
 LOGICAL :: CHECK_VN=.TRUE.                  !< Check the Von Neumann number
+LOGICAL :: CHECK_FO=.FALSE.                 !< Check the solid phase Fourier number
 LOGICAL :: SOLID_PARTICLES=.FALSE.          !< Indicates the existence of solid particles
 LOGICAL :: HVAC=.FALSE.                     !< Perform an HVAC calculation
 LOGICAL :: BAROCLINIC=.TRUE.                !< Include the baroclinic terms in the momentum equation

Source/read.f90

@@ -1717,7 +1717,7 @@ SUBROUTINE READ_MISC
 
 NAMELIST /MISC/ AEROSOL_AL2O3,AEROSOL_SCRUBBING,AGGLOMERATION,ALIGNMENT_TOLERANCE, &
                 BNDF_DEFAULT,CC_IBM,CCVOL_LINK,C_DEARDORFF,&
-                CFL_MAX,CFL_MIN,CFL_VELOCITY_NORM,CHECK_HT,CHECK_VN, &
+                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,&

Source/velo.f90

@@ -3051,7 +3051,7 @@ SUBROUTINE CHECK_STABILITY(DT,DT_NEW,NM)
       IIG = BC%IIG
       JJG = BC%JJG
       KKG = BC%KKG
-      UVW = (ABS(B1%Q_CON_F)/B1%RHO_F)**ONTH * 2._EB*B1%RDN
+      UVW = (ABS(B1%Q_RAD_IN-B1%Q_RAD_OUT+B1%Q_CON_F)/B1%RHO_F)**ONTH * 2._EB*B1%RDN
       IF (UVW>=UVWMAX) THEN
          UVWMAX = UVW
          ICFL=IIG

Source/wall.f90

@@ -1686,7 +1686,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
             Q_WATER_F,Q_WATER_B,LAYER_DIVIDE,TMP_GAS_BACK,GEOM_FACTOR,DT_BC_SUB_OLD,&
             DEL_DOT_Q_SC,Q_DOT_G_PP,Q_DOT_G_PP_NET,Q_DOT_O2_PP,Q_DOT_O2_PP_NET,R_SURF,U_SURF,V_SURF,W_SURF,T_BC_SUB,DT_BC_SUB,&
             Q_NET_F,Q_NET_B,TMP_RATIO,KODXF,KODXB,H_S,T_NODE,C_S,H_NODE,VOL,T_BOIL_EFF,&
-            RADIUS,HTC_LIMIT,CP1,CP2,DENOM,SF_HTC_F,SF_HTC_B,DDSUM,THICKNESS
+            RADIUS,HTC_LIMIT,CP1,CP2,DENOM,SF_HTC_F,SF_HTC_B,DDSUM,THICKNESS,DT_FO
 REAL(EB), DIMENSION(N_TRACKED_SPECIES) :: M_DOT_G_PP_ADJUST,M_DOT_G_PP_ADJUST_NET,M_DOT_G_PP_ACTUAL,M_DOT_G_PP_ACTUAL_NET
 REAL(EB), DIMENSION(MAX_MATERIALS) :: M_DOT_S_PP,M_DOT_S_PP_NET
 REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART_S,M_DOT_PART_S
@@ -1927,6 +1927,34 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
 B1%N_SUBSTEPS = 1
 ONE_D%DELTA_TMP = 0._EB
 
+! Compute initial thermal properties for Fo number calculation
+! CHECK_FO=F by default
+! CHECK_FO=T enforces near explicit time step accuracy (Forward-Euler would require 1/2 DT_FO)
+
+CHECK_FO_IF: IF (CHECK_FO) THEN
+   POINT_LOOP0: DO I=1,NWP
+      VOLSUM = 0._EB
+      ITMP = MIN(I_MAX_TEMP-1,INT(ONE_D%TMP(I)))
+      MATERIAL_LOOP0: DO N=1,ONE_D%N_MATL
+         IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE MATERIAL_LOOP0
+         ML  => MATERIAL(ONE_D%MATL_INDEX(N))
+         VOLSUM = VOLSUM + ONE_D%MATL_COMP(N)%RHO(I)/RHO_ADJUSTED(LAYER_INDEX(I),N)
+         ONE_D%K_S(I) = ONE_D%K_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%K_S(ITMP)/RHO_ADJUSTED(LAYER_INDEX(I),N)
+         ONE_D%RHO_C_S(I) = ONE_D%RHO_C_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%C_S(ITMP)
+
+         IF (.NOT.E_FOUND) B1%EMISSIVITY = B1%EMISSIVITY + ONE_D%MATL_COMP(N)%RHO(I)*ML%EMISSIVITY/RHO_ADJUSTED(LAYER_INDEX(I),N)
+         RHO_S(I) = RHO_S(I) + ONE_D%MATL_COMP(N)%RHO(I)
+      ENDDO MATERIAL_LOOP0
+      IF (SF%PACKING_RATIO(LAYER_INDEX(I))>0._EB) ONE_D%K_S(I) = ONE_D%K_S(I)*SF%PACKING_RATIO(LAYER_INDEX(I))
+
+      IF (VOLSUM > 0._EB) THEN
+         ONE_D%K_S(I) = ONE_D%K_S(I)/VOLSUM
+      ENDIF
+      IF (ONE_D%K_S(I)<=TWO_EPSILON_EB)      ONE_D%K_S(I)      = 10000._EB
+      IF (ONE_D%RHO_C_S(I)<=TWO_EPSILON_EB)  ONE_D%RHO_C_S(I)  = 0.001_EB
+   ENDDO POINT_LOOP0
+ENDIF CHECK_FO_IF
+
 SUB_TIMESTEP_LOOP: DO
 
    ! Compute grid for reacting nodes
@@ -1950,6 +1978,14 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
       MF_FRAC(1:NWP)       = SF%MF_FRAC(1:NWP)
    ENDIF COMPUTE_GRID
 
+   ! Calculate minimum DT based on Fourier number, Fo
+
+   IF (CHECK_FO) THEN
+      DT_FO = MINVAL( DX_S(1:NWP)**2 * ONE_D%RHO_C_S(1:NWP) / ONE_D%K_S(1:NWP) )
+   ELSE
+      DT_FO = HUGE(1._EB)
+   ENDIF
+
    ! Compute convective heat flux at the surface
 
    DTMP = B1%TMP_G - B1%TMP_F
@@ -2173,7 +2209,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
       TMP_RATIO = MAX(TWO_EPSILON_EB,MAXVAL(ABS(DELTA_TMP(1:NWP)))/SF%DELTA_TMP_MAX)
       DT_BC_SUB_OLD = DT_BC_SUB
       DT_BC_SUB = DT_BC/REAL(MIN(NINT(SF%TIME_STEP_FACTOR*WALL_INCREMENT),MAX(1,NINT(TMP_RATIO))),EB)
-      DT_BC_SUB = MIN( DT_BC-T_BC_SUB , DT_BC_SUB )
+      DT_BC_SUB = MIN( DT_BC-T_BC_SUB , DT_BC_SUB , DT_FO )
       IF (SF%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED .AND. DT_BC_SUB_OLD/=DT_BC_SUB) CALL PERFORM_PYROLYSIS
    ENDIF
 
@@ -2615,12 +2651,11 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
       IF (B1%EMISSIVITY>=0._EB) E_FOUND = .TRUE.
       IF (ONE_D%K_S(I)<=TWO_EPSILON_EB)      ONE_D%K_S(I)      = 10000._EB
       IF (ONE_D%RHO_C_S(I)<=TWO_EPSILON_EB)  ONE_D%RHO_C_S(I)  = 0.001_EB
-
    ENDDO POINT_LOOP3
 
    IF (SF%EMISSIVITY_SPECIFIED) B1%EMISSIVITY = SF%EMISSIVITY
 
-   ! Calculate average K_S between at grid cell boundaries. Store result in K_S
+   ! Calculate average K_S at grid cell boundaries. Store result in K_S
 
    ONE_D%K_S(0)     = ONE_D%K_S(1)
    ONE_D%K_S(NWP+1) = ONE_D%K_S(NWP)