Merge pull request #14823 from shostikk/master

Surf layer delamination (fall-off), based on density or temperature

Author: shostikk

Source/read.f90

@@ -7809,7 +7809,8 @@ SUBROUTINE READ_SURF(QUICK_READ)
             DRAG_COEFFICIENT,MINIMUM_BURNOUT_TIME,DELTA_TMP_MAX,BURN_DURATION,&
             REFERENCE_HEAT_FLUX(MAX_QDOTPP_REF),REFERENCE_HEAT_FLUX_TIME_INTERVAL,MINIMUM_SCALING_HEAT_FLUX,&
             MAXIMUM_SCALING_HEAT_FLUX,REFERENCE_THICKNESS(MAX_QDOTPP_REF),&
-            AREA_MULTIPLIER,Z_0,PARTICLE_EXTRACTION_VELOCITY,NEAR_WALL_EDDY_VISCOSITY
+            AREA_MULTIPLIER,Z_0,PARTICLE_EXTRACTION_VELOCITY,NEAR_WALL_EDDY_VISCOSITY,&
+            FALLOFF_TMP(MAX_LAYERS),FALLOFF_DENSITY(MAX_LAYERS)
 INTEGER :: NPPC,N,IOS,NL,NN,NNN,NNNN,N_LIST,LEAK_PATH(2),DUCT_PATH(2),RGB(3),NR,IL
 INTEGER ::  N_LAYER_CELLS_MAX(MAX_LAYERS),VEG_LSET_FUEL_INDEX,INDEX_LIST(MAX_MATERIALS**2)
 INTEGER :: CHILD_LAYER(MAX_LAYERS,N_MATL),CHILD_SURF(N_MATL)
@@ -7828,6 +7829,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
                 DSC_CONVERSION_FACTOR,DT_INSERT,E_COEFFICIENT,&
                 EMBER_GENERATION_HEIGHT,EMBER_IGNITION_POWER_MEAN,EMBER_IGNITION_POWER_SIGMA,EMBER_TRACKING_RATIO,EMBER_YIELD,&
                 EMISSIVITY,EMISSIVITY_BACK,EXTERNAL_FLUX,EXTINCTION_TEMPERATURE,&
+                FALLOFF_DENSITY,FALLOFF_TMP,&
                 FREE_SLIP,INERT_Q_REF,FYI,GEOMETRY,HEAT_OF_VAPORIZATION,&
                 HEAT_TRANSFER_COEFFICIENT,HEAT_TRANSFER_COEFFICIENT_BACK,HEAT_TRANSFER_COEFFICIENT_SIGMA,&
                 HEAT_TRANSFER_MODEL,HORIZONTAL,HRRPUA,VARIABLE_THICKNESS,HT3D,ID,IGNITION_TEMPERATURE,&
@@ -8761,6 +8763,16 @@ SUBROUTINE READ_SURF(QUICK_READ)
       ENDIF
    ENDDO
 
+   ! Set delamination temperature and/or density
+
+   SF%FALLOFF_TMP = HUGE_EB
+   SF%FALLOFF_DENSITY = -1._EB
+   DO NL = 1,MAX_LAYERS
+      IF (FALLOFF_TMP(NL) > 0._EB) SF%FALLOFF_TMP(NL) = FALLOFF_TMP(NL) + TMPM
+      IF (FALLOFF_DENSITY(NL) > 0._EB) SF%FALLOFF_DENSITY(NL) = FALLOFF_DENSITY(NL)
+   ENDDO
+   
+      
    ! Adjust min and max temperature and convert C to K
 
    IF (TMP_FRONT >= -TMPM) TMPMIN = MIN(TMPMIN,TMP_FRONT+TMPM)
@@ -8961,6 +8973,8 @@ SUBROUTINE SET_SURF_DEFAULTS
 EXTERNAL_FLUX           = 0._EB
 EXTINCTION_TEMPERATURE  = -273._EB
 INERT_Q_REF             = .FALSE.
+FALLOFF_TMP             = -1._EB
+FALLOFF_DENSITY         = -1._EB
 FREE_SLIP               = .FALSE.
 FYI                     = 'null'
 GEOMETRY                = 'CARTESIAN'

Source/type.f90

@@ -964,7 +964,8 @@ MODULE TYPES
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: CELL_SIZE           !< Specified constant cell size (m)
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: STRETCH_FACTOR
    REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_DENSITY,&
-                                      MOISTURE_CONTENT,SURFACE_VOLUME_RATIO,PACKING_RATIO,KAPPA_S=-1._EB,RENODE_DELTA_T
+        MOISTURE_CONTENT,SURFACE_VOLUME_RATIO,PACKING_RATIO,KAPPA_S=-1._EB,RENODE_DELTA_T
+   REAL(EB), DIMENSION(MAX_LAYERS) :: FALLOFF_TMP,FALLOFF_DENSITY !< Temp. and density criteria for layer fall-off (delam.)
    REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: DENSITY_ADJUST_FACTOR=1._EB,RHO_S
    CHARACTER(LABEL_LENGTH), ALLOCATABLE, DIMENSION(:) :: MATL_NAME
    CHARACTER(LABEL_LENGTH), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_ID

Source/wall.f90

@@ -2043,13 +2043,22 @@ 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 number of cells and density vector
+NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS))
+RHO_S = 0._EB
+DO I = 1,NWP
+   DO N=1,ONE_D%N_MATL
+      RHO_S(I)=RHO_S(I)+ONE_D%MATL_COMP(N)%RHO(I)
+   ENDDO
+ENDDO
+
 ! 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
    ONE_D%K_S = 0._EB
-   RHO_S   = 0._EB
    ONE_D%RHO_C_S = 0._EB
    POINT_LOOP0: DO I=1,NWP
       VOLSUM = 0._EB
@@ -2060,7 +2069,6 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
          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)
-         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))
 
@@ -2079,12 +2087,11 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
    LAYER_DIVIDE = SF%LAYER_DIVIDE
 
    COMPUTE_GRID: IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED .OR. SF%HT_DIM>1 .OR. SF%VARIABLE_THICKNESS) THEN
-      NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS))
       CALL GET_WALL_NODE_WEIGHTS(NWP,ONE_D%N_LAYERS,ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS),ONE_D%LAYER_THICKNESS,SF%GEOMETRY, &
          ONE_D%X(0:NWP),LAYER_DIVIDE,DX_S(1:NWP),RDX_S(0:NWP+1),RDXN_S(0:NWP),DX_WGT_S(0:NWP),DXF,DXB,&
          LAYER_INDEX(0:NWP+1),MF_FRAC(1:NWP),SF%INNER_RADIUS,ONE_D%LAYER_DIVIDE_DEPTH)
    ELSE COMPUTE_GRID
-      NWP                  = SF%N_CELLS_INI
+!      NWP                  = SF%N_CELLS_INI
       DXF                  = SF%DXF
       DXB                  = SF%DXB
       DX_S(1:NWP)          = SF%DX(1:NWP)
@@ -2413,6 +2420,25 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
 
       ENDDO POINT_LOOP2
 
+      ! Delamination layer fall-off
+
+      I = 0
+      DO NL=1,ONE_D%N_LAYERS
+         IF (ONE_D%LAYER_THICKNESS(NL).GT.ONE_D%MIN_LAYER_THICKNESS(NL) .AND. &
+              (ONE_D%TMP(I+ONE_D%N_LAYER_CELLS(NL)) > SF%FALLOFF_TMP(NL)) .OR.&
+              (    RHO_S(I+ONE_D%N_LAYER_CELLS(NL)) < SF%FALLOFF_DENSITY(NL))) THEN
+            REGRID_FACTOR(1:I+ONE_D%N_LAYER_CELLS(NL))=0._EB
+            ONE_D%RHO_C_S(1:I+ONE_D%N_LAYER_CELLS(NL))=0._EB
+            Q_S(1:I+ONE_D%N_LAYER_CELLS(NL))=0._EB
+            DO N=1,ONE_D%N_MATL
+               ONE_D%MATL_COMP(N)%RHO(1:I+ONE_D%N_LAYER_CELLS(NL))=0._EB
+            ENDDO
+            ONE_D%TMP(I+1:I+ONE_D%N_LAYER_CELLS(NL))=ONE_D%TMP(I+ONE_D%N_LAYER_CELLS(NL)+1)
+            B1%TMP_F = ONE_D%TMP(I+ONE_D%N_LAYER_CELLS(NL)+1)
+         ENDIF
+         I = I + ONE_D%N_LAYER_CELLS(NL)
+      ENDDO
+
       ! Compute new coordinates if the solid changes thickness. Save new coordinates in X_S_NEW.
       ! Remesh layer if any node goes to zero thickness
 
@@ -2465,7 +2491,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
          ENDIF
          I = I + ONE_D%N_LAYER_CELLS(NL)
       ENDDO
-      
+
       ! If any nodes go to zero, apportion Q_S to surrounding nodes.
 
       IF (ANY(CELL_ZERO(1:ONE_D%N_LAYERS)) .AND. NWP > 1) THEN
@@ -2661,7 +2687,6 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
          CALL INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,Q_S(1:N_CELLS))
          CALL INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,RHO_H_S(1:N_CELLS))
          CALL INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,TMP_S(1:N_CELLS))
-
          DO I=1,NWP_NEW
             VOL = (THICKNESS+SF%INNER_RADIUS-X_S_NEW(I-1))**I_GRAD-(THICKNESS+SF%INNER_RADIUS-X_S_NEW(I))**I_GRAD
             TMP_S(I) = TMP_S(I) / VOL