diff --git a/Source/init.f90 b/Source/init.f90 index 45c7295465d..3445d53a29f 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -1587,9 +1587,8 @@ SUBROUTINE REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL,THIN_WALL_CELL) INTEGER, ALLOCATABLE, DIMENSION(:) :: LAYER_INDEX INTEGER, ALLOCATABLE, DIMENSION(:) :: N_LAYER_CELLS_OLD REAL(EB) :: MIN_DENSITY -REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_DENSITY TYPE(MATERIAL_TYPE), POINTER :: ML -REAL(EB), ALLOCATABLE, DIMENSION(:) :: X_S_OLD +REAL(EB), ALLOCATABLE, DIMENSION(:) :: X_S_OLD,LAYER_DENSITY LOGICAL, ALLOCATABLE, DIMENSION(:) :: REMESH_LAYER TYPE(WALL_TYPE), POINTER :: WC TYPE(THIN_WALL_TYPE), POINTER :: TW @@ -1648,7 +1647,7 @@ SUBROUTINE REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL,THIN_WALL_CELL) ONE_D%N_CELLS_INI = 0 ONE_D%N_CELLS_MAX = 0 -LAYER_DENSITY = 0._EB +ALLOCATE(LAYER_DENSITY(ONE_D%N_LAYERS)) ; LAYER_DENSITY = 0._EB LAYER_LOOP: DO NL=1,ONE_D%N_LAYERS @@ -1754,6 +1753,7 @@ SUBROUTINE REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL,THIN_WALL_CELL) ENDDO MATERIAL_LOOP3 ENDDO POINT_LOOP3 +DEALLOCATE(LAYER_DENSITY) DEALLOCATE(LAYER_INDEX) ! Reset emissivity at the surface @@ -3828,13 +3828,13 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW) INTEGER :: II,JJ,KK,IC,ICG,IOR,NOM,ITER,OBST_INDEX,OBST_INDEX_PREVIOUS,NN,NNN,NL,N_LAYERS_OBST,& N_MATL_OBST,N_LAYERS,N_MATLS,IIF,JJF,KKF INTEGER, DIMENSION(MAX_MATERIALS) :: MATL_INDEX_OBST,MATL_INDEX -REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION_OBST,MATL_MASS_FRACTION -REAL(EB), DIMENSION(0:MAX_LAYERS) :: LAYER_THICKNESS,MIN_LAYER_THICKNESS -REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_THICKNESS_OBST,HEAT_SOURCE,HEAT_SOURCE_OBST,& +REAL(EB), DIMENSION(MAX_LAYERS_HT3D,MAX_MATERIALS) :: MATL_MASS_FRACTION_OBST,MATL_MASS_FRACTION +REAL(EB), DIMENSION(MAX_LAYERS_HT3D) :: LAYER_THICKNESS,MIN_LAYER_THICKNESS,& + LAYER_THICKNESS_OBST,HEAT_SOURCE,HEAT_SOURCE_OBST,& STRETCH_FACTOR,STRETCH_FACTOR_OBST,CELL_SIZE,CELL_SIZE_OBST,& CELL_SIZE_FACTOR,CELL_SIZE_FACTOR_OBST -INTEGER, DIMENSION(MAX_LAYERS) :: N_LAYER_CELLS_MAX,N_LAYER_CELLS_MAX_OBST,RAMP_IHS_INDEX,RAMP_IHS_INDEX_OBST -LOGICAL, DIMENSION(MAX_LAYERS) :: HT3D_LAYER +INTEGER, DIMENSION(MAX_LAYERS_HT3D) :: N_LAYER_CELLS_MAX,N_LAYER_CELLS_MAX_OBST,RAMP_IHS_INDEX,RAMP_IHS_INDEX_OBST +LOGICAL, DIMENSION(MAX_LAYERS_HT3D) :: HT3D_LAYER REAL(EB) :: XXC,YYC,ZZC,THICKNESS,OLD_THICKNESS,FRONT_LINING_THICKNESS,BACK_LINING_THICKNESS,LAYER_THICKNESS_OBST_TOTAL CHARACTER(MESSAGE_LENGTH) :: MESSAGE LOGICAL :: THIN_OBSTRUCTION @@ -3994,10 +3994,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW) IF (OBST_INDEX/=OBST_INDEX_PREVIOUS .AND. OBST_INDEX_PREVIOUS>0 .AND. OBST_INDEX>0) THEN OB_PREV => OM_PREV%OBSTRUCTION(OBST_INDEX_PREVIOUS) - IF ( (ANY(OB%MATL_MASS_FRACTION(:)/=OB_PREV%MATL_MASS_FRACTION(:),DIM=1)) .OR. & - (ANY(OB%MATL_INDEX(:) /=OB_PREV%MATL_INDEX(:) ,DIM=1)) .OR. & - (OB%HEAT_SOURCE/=OB_PREV%HEAT_SOURCE) .OR. & - (OB%RAMP_IHS_INDEX/=OB_PREV%RAMP_IHS_INDEX) ) THEN + IF (OB%ORDINAL/=OB_PREV%ORDINAL) THEN N_LAYERS_OBST = N_LAYERS_OBST + 1 LAYER_THICKNESS_OBST(N_LAYERS_OBST) = 0._EB HEAT_SOURCE_OBST(N_LAYERS_OBST) = OB%HEAT_SOURCE @@ -4238,7 +4235,7 @@ SUBROUTINE FIND_THIN_WALL_BACK_INDEX(NM,ITW) INTEGER, INTENT(IN) :: NM,ITW INTEGER :: II,JJ,KK,IC,IOR,IEC,ITW2,NOM,NN,NNN,NL,N_MATLS,IIGM,IIGP,JJGM,JJGP,KKGM,KKGP,ICM,ICP INTEGER, DIMENSION(MAX_MATERIALS) :: MATL_INDEX -REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION +REAL(EB), DIMENSION(MAX_LAYERS_HT3D,MAX_MATERIALS) :: MATL_MASS_FRACTION REAL(EB) :: XXC,YYC,ZZC TYPE (MESH_TYPE), POINTER :: M TYPE (THIN_WALL_TYPE), POINTER :: TW diff --git a/Source/prec.f90 b/Source/prec.f90 index c9e7512f9ef..babd4e30025 100644 --- a/Source/prec.f90 +++ b/Source/prec.f90 @@ -10,6 +10,7 @@ MODULE PRECISION_PARAMETERS INTEGER, PARAMETER :: MAX_LPC=20 !< Maximum number of declared particle classes INTEGER, PARAMETER :: MAX_SPECIES=20 !< Maximum number of declared species INTEGER, PARAMETER :: MAX_LAYERS=20 !< Maximum number of solid material layers +INTEGER, PARAMETER :: MAX_LAYERS_HT3D=500 !< Maximum number of solid material layers for an HT3D solid INTEGER, PARAMETER :: MAX_MATERIALS=20 !< Maximum number of solid material components INTEGER, PARAMETER :: MAX_MATERIALS_TOTAL=400 !< Dimension of material work array INTEGER, PARAMETER :: MAX_CONE_CURVES=10 !< Maximum number of cone calorimeter curves diff --git a/Source/wall.f90 b/Source/wall.f90 index 92c779b915f..b852ffff069 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -1593,7 +1593,7 @@ SUBROUTINE DEPOSIT_PARTICLE_MASS(LP,LPC) IF (SF%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN ONE_D => BOUNDARY_ONE_D(LP%OD_INDEX) - LP%RADIUS = SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS)) + LP%RADIUS = SUM(ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS)) ENDIF END SUBROUTINE DEPOSIT_PARTICLE_MASS @@ -1838,18 +1838,18 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART_S,M_DOT_PART_S REAL(EB), DIMENSION(NWP_MAX) :: TMP_S,RHO_H_S REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: RHO_DOT,INT_WGT -REAL(EB), DIMENSION(MAX_LAYERS) :: DX_MIN -REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: RHO_ADJUSTED +REAL(EB), DIMENSION(MAX_LAYERS_HT3D) :: DX_MIN +REAL(EB), DIMENSION(MAX_LAYERS_HT3D,MAX_MATERIALS) :: RHO_ADJUSTED REAL(EB), DIMENSION(NWP_MAX) :: AAS,BBS,CCS,DDS,DDT,Q_S,Q_IR,Q_ADD,TWO_DX_KAPPA_S,DX_S,MF_FRAC,REGRID_FACTOR REAL(EB), DIMENSION(0:NWP_MAX+1) :: RHO_S,DELTA_TMP,RDX_S REAL(EB), DIMENSION(0:NWP_MAX) :: X_S_NEW,RDXN_S,R_S,R_S_NEW,DX_WGT_S INTEGER, DIMENSION(0:NWP_MAX+1) :: LAYER_INDEX -INTEGER, DIMENSION(MAX_LAYERS) :: N_LAYER_CELLS_NEW +INTEGER, DIMENSION(MAX_LAYERS_HT3D) :: N_LAYER_CELLS_NEW INTEGER :: NWP_NEW,I_GRAD,IZERO,SURF_INDEX,BACKING,NWP,I,NL,N,OBST_INDEX,& N_CELLS,ITMP,ITER,BACK_MESH,BACK_INDEX,BACK_WALL_INDEX CHARACTER(MESSAGE_LENGTH) :: MESSAGE -LOGICAL :: ISOLATED_THIN_WALL,ISOLATED_THIN_WALL_BACK,REMESH_LAYER(MAX_LAYERS),REMESH_CHECK,DIRICHLET_BACK,& - CELL_ZERO(MAX_LAYERS),TMP_CHECK(MAX_LAYERS) +LOGICAL :: ISOLATED_THIN_WALL,ISOLATED_THIN_WALL_BACK,REMESH_LAYER(MAX_LAYERS_HT3D),REMESH_CHECK,DIRICHLET_BACK,& + CELL_ZERO(MAX_LAYERS_HT3D),TMP_CHECK(MAX_LAYERS_HT3D) TYPE(WALL_TYPE), POINTER :: WC,WC_BACK TYPE(THIN_WALL_TYPE), POINTER :: TW,TW_BACK TYPE(CFACE_TYPE), POINTER :: CFA,CFA_BACK @@ -2014,7 +2014,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, DO NL=1,ONE_D%N_LAYERS DO N=1,ONE_D%N_MATL - RHO_ADJUSTED(NL,N) = MATERIAL(ONE_D%MATL_INDEX(N))%RHO_S*SF%DENSITY_ADJUST_FACTOR(NL,N) + RHO_ADJUSTED(NL,N) = MATERIAL(ONE_D%MATL_INDEX(N))%RHO_S*SF%DENSITY_ADJUST_FACTOR(MIN(MAX_LAYERS,NL),N) ENDDO ENDDO @@ -2071,7 +2071,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, 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) 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 (SF%BOUNDARY_FUEL_MODEL) 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 @@ -2121,7 +2121,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, ELSEIF (PRESENT(CFACE_INDEX)) THEN HTCF = B1%HEAT_TRANS_COEF ELSEIF (PRESENT(PARTICLE_INDEX)) THEN - RADIUS = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS)) + RADIUS = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS)) SELECT CASE(SF%GEOMETRY) CASE (SURF_CARTESIAN) ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/( DT_BC_SUB) CASE (SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/(2._EB*DT_BC_SUB) @@ -2158,7 +2158,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, Q_CON_B = HTCB*DTMP Q_RAD_IN_B = EMISSIVITY_BACK*SIGMA*TMP_GAS_BACK**4 Q_LIQUID_B = 0._EB - LAYER_DIVIDE = REAL(SF%N_LAYERS+1) + LAYER_DIVIDE = REAL(ONE_D%N_LAYERS+1) MF_FRAC = 1._EB CASE(INSULATED) ! No heat transfer out the back @@ -2774,7 +2774,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, 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_LOOP3 - 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 (SF%BOUNDARY_FUEL_MODEL) ONE_D%K_S(I) = ONE_D%K_S(I)*SF%PACKING_RATIO(LAYER_INDEX(I)) IF (VOLSUM > 0._EB) ONE_D%K_S(I) = ONE_D%K_S(I)/VOLSUM 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 @@ -3387,7 +3387,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F MFLUX = MIN(MFLUX_MAX,MFLUX + RHO_DOT_EXTRA*DX_S(SOLID_CELL_INDEX)) RHO_DOT_REAC(J) =MFLUX/DX_S(SOLID_CELL_INDEX) ! kg/m3/s CASE(SURF_SPHERICAL) - NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS)) + NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS)) R_S_0 = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(0) R_S_1 = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(1) DR = (R_S_0**3-R_S_1**3)/(3._EB*R_S_0**2) @@ -3445,7 +3445,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F IF (SF%BOUNDARY_FUEL_MODEL) THEN LENGTH_SCALE = 1._EB/(SF%SURFACE_VOLUME_RATIO(LAYER_INDEX)*SF%PACKING_RATIO(LAYER_INDEX)) ELSE - LENGTH_SCALE = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS)) + LENGTH_SCALE = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS)) SELECT CASE(SF%GEOMETRY) CASE(SURF_SPHERICAL) LENGTH_SCALE = LENGTH_SCALE/3._EB @@ -3808,7 +3808,7 @@ SUBROUTINE HT3D_TEMPERATURE_EXCHANGE(NM) BC => M%BOUNDARY_COORD(TW%BC_INDEX) IF (ABS(BC%IOR)==M%HT_3D_SWEEP_DIRECTION) CYCLE THIN_WALL_LOOP ONE_D => M%BOUNDARY_ONE_D(TW%OD_INDEX) - NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS)) + NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS)) THR_D => M%BOUNDARY_THR_D(TW%TD_INDEX) IF (.NOT.ALLOCATED(THR_D%NODE)) CYCLE THIN_WALL_LOOP