diff --git a/Source/cons.f90 b/Source/cons.f90 index abb4f456164..de61b6d5902 100644 --- a/Source/cons.f90 +++ b/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 diff --git a/Source/read.f90 b/Source/read.f90 index d55f9df2cab..2fd79c495d1 100644 --- a/Source/read.f90 +++ b/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,& diff --git a/Source/velo.f90 b/Source/velo.f90 index cd419dff888..76e2be735c5 100644 --- a/Source/velo.f90 +++ b/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 diff --git a/Source/wall.f90 b/Source/wall.f90 index 4f98dd0cd5c..3bdb09829cf 100644 --- a/Source/wall.f90 +++ b/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)