diff --git a/Source/read.f90 b/Source/read.f90 index 6f27715b7cf..f0a6a2178b3 100644 --- a/Source/read.f90 +++ b/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' diff --git a/Source/type.f90 b/Source/type.f90 index e2148bd8ab9..dde20cdf7c4 100644 --- a/Source/type.f90 +++ b/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 diff --git a/Source/wall.f90 b/Source/wall.f90 index e6ecbaff84d..daa8be850fc 100644 --- a/Source/wall.f90 +++ b/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