diff --git a/Source/func.f90 b/Source/func.f90 index 03336f54f1f..437a3300cf3 100644 --- a/Source/func.f90 +++ b/Source/func.f90 @@ -1840,7 +1840,7 @@ REAL(EB) FUNCTION GET_PARTICLE_ENTHALPY(I_LPC,TMP_S) USE MATH_FUNCTIONS, ONLY: INTERPOLATE1D_UNIFORM REAL(EB), INTENT(IN) :: TMP_S REAL(EB) :: RHO,RHO_H,VOL,DTMP,H_S,THICKNESS -INTEGER :: I,N,ITMP,I_GRAD +INTEGER :: I,N,ITMP INTEGER, INTENT(IN) :: I_LPC TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC TYPE(SURFACE_TYPE), POINTER :: SF @@ -1852,11 +1852,6 @@ REAL(EB) FUNCTION GET_PARTICLE_ENTHALPY(I_LPC,TMP_S) CALL INTERPOLATE1D_UNIFORM(LBOUND(SPECIES(LPC%Y_INDEX)%H_L,1),SPECIES(LPC%Y_INDEX)%H_L,TMP_S,GET_PARTICLE_ENTHALPY) ELSE SF=>SURFACE(LPC%SURF_INDEX) - SELECT CASE(SF%GEOMETRY) - CASE(SURF_CARTESIAN) ; I_GRAD = 1 - CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; I_GRAD = 2 - CASE(SURF_SPHERICAL) ; I_GRAD = 3 - END SELECT RHO_H = 0._EB RHO = 0._EB ITMP = MIN(I_MAX_TEMP-1,INT(TMP_S)) @@ -1864,9 +1859,9 @@ REAL(EB) FUNCTION GET_PARTICLE_ENTHALPY(I_LPC,TMP_S) THICKNESS = SUM(SF%LAYER_THICKNESS) DO I=1,SUM(SF%N_LAYER_CELLS) IF (SF%GEOMETRY==SURF_INNER_CYLINDRICAL) THEN - VOL = (SF%INNER_RADIUS+SF%X_S(I))**I_GRAD - (SF%INNER_RADIUS+SF%X_S(I-1))**I_GRAD + VOL = (SF%INNER_RADIUS+SF%X_S(I))**SF%I_GRAD - (SF%INNER_RADIUS+SF%X_S(I-1))**SF%I_GRAD ELSE - VOL = (THICKNESS+SF%INNER_RADIUS-SF%X_S(I-1))**I_GRAD - (THICKNESS+SF%INNER_RADIUS-SF%X_S(I))**I_GRAD + VOL = (THICKNESS+SF%INNER_RADIUS-SF%X_S(I-1))**SF%I_GRAD - (THICKNESS+SF%INNER_RADIUS-SF%X_S(I))**SF%I_GRAD ENDIF MATL_REMESH: DO N=1,SF%N_MATL IF (SF%RHO_0(I,N)<=TWO_EPSILON_EB) CYCLE MATL_REMESH @@ -2112,7 +2107,7 @@ REAL(EB) FUNCTION SURFACE_DENSITY(MODE,SF,ONE_D,MATL_INDEX) INTEGER, INTENT(IN) :: MODE INTEGER, INTENT(IN), OPTIONAL :: MATL_INDEX -INTEGER :: I_GRAD,NWP,II2,N,ITMP +INTEGER :: NWP,II2,N,ITMP REAL(EB) :: WGT,R_S(0:NWP_MAX),EPUM,DTMP TYPE(BOUNDARY_ONE_D_TYPE), INTENT(IN), POINTER :: ONE_D TYPE(SURFACE_TYPE), INTENT(IN), POINTER :: SF @@ -2124,12 +2119,6 @@ REAL(EB) FUNCTION SURFACE_DENSITY(MODE,SF,ONE_D,MATL_INDEX) ELSE THERMALLY_THICK_IF - SELECT CASE(SF%GEOMETRY) - CASE(SURF_CARTESIAN) ; I_GRAD = 1 - CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; I_GRAD = 2 - CASE(SURF_SPHERICAL) ; I_GRAD = 3 - END SELECT - NWP = SUM(ONE_D%N_LAYER_CELLS) IF (SF%GEOMETRY==SURF_INNER_CYLINDRICAL) THEN R_S(0:NWP) = SF%INNER_RADIUS + ONE_D%X(0:NWP) @@ -2140,8 +2129,9 @@ REAL(EB) FUNCTION SURFACE_DENSITY(MODE,SF,ONE_D,MATL_INDEX) SURFACE_DENSITY = 0._EB NUMBER_WALL_POINTS_LOOP: DO II2=1,NWP AREA_VOLUME_SELECT: SELECT CASE(MODE) - CASE(0,2); WGT = ABS(R_S(II2-1)**I_GRAD-R_S(II2)**I_GRAD)/(REAL(I_GRAD,EB)*(SF%INNER_RADIUS+SF%THICKNESS)**(I_GRAD-1)) - CASE(1,3); WGT = ABS(R_S(II2-1)**I_GRAD-R_S(II2)**I_GRAD)/(SF%INNER_RADIUS+SF%THICKNESS)**I_GRAD + CASE(0,2); WGT = ABS(R_S(II2-1)**SF%I_GRAD-R_S(II2)**SF%I_GRAD)/ & + (REAL(SF%I_GRAD,EB)*(SF%INNER_RADIUS+SF%THICKNESS)**(SF%I_GRAD-1)) + CASE(1,3); WGT = ABS(R_S(II2-1)**SF%I_GRAD-R_S(II2)**SF%I_GRAD)/(SF%INNER_RADIUS+SF%THICKNESS)**SF%I_GRAD END SELECT AREA_VOLUME_SELECT EPUM = 1._EB ! energy per unit mass @@ -4063,6 +4053,7 @@ SUBROUTINE PACK_BOUNDARY_ONE_D(NM,IC,RC,LC,OS,OD_INDEX,UNPACK_IT,COUNT_ONLY,CHEC I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%LAYER_THICKNESS(1:RC-I1+1) , UNPACK_IT) I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%LAYER_THICKNESS_OLD(1:RC-I1+1), UNPACK_IT) I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%MIN_LAYER_THICKNESS(1:RC-I1+1), UNPACK_IT) +I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%MIN_LAYER_MASS(1:RC-I1+1), UNPACK_IT) I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%MIN_DIFFUSIVITY(1:RC-I1+1) , UNPACK_IT) I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%DDSUM(1:RC-I1+1) , UNPACK_IT) I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%SMALLEST_CELL_SIZE(1:RC-I1+1) , UNPACK_IT) @@ -4106,6 +4097,7 @@ SUBROUTINE REALLOCATE_BOUNDARY_ONE_D(ONE_D) ALLOCATE(ONE_D%LAYER_THICKNESS_OLD(ONE_D%N_LAYERS)) IF (ALLOCATED(ONE_D%MIN_LAYER_THICKNESS)) DEALLOCATE(ONE_D%MIN_LAYER_THICKNESS) ALLOCATE(ONE_D%MIN_LAYER_THICKNESS(ONE_D%N_LAYERS)) +IF (ALLOCATED(ONE_D%MIN_LAYER_MASS)) DEALLOCATE(ONE_D%MIN_LAYER_MASS) ; ALLOCATE(ONE_D%MIN_LAYER_MASS(ONE_D%N_LAYERS)) IF (ALLOCATED(ONE_D%HT3D_LAYER)) DEALLOCATE(ONE_D%HT3D_LAYER) ; ALLOCATE(ONE_D%HT3D_LAYER(ONE_D%N_LAYERS)) IF (ALLOCATED(ONE_D%MIN_DIFFUSIVITY)) DEALLOCATE(ONE_D%MIN_DIFFUSIVITY) ; ALLOCATE(ONE_D%MIN_DIFFUSIVITY(ONE_D%N_LAYERS)) IF (ALLOCATED(ONE_D%RHO_C_S)) DEALLOCATE(ONE_D%RHO_C_S) ; ALLOCATE(ONE_D%RHO_C_S(ONE_D%N_CELLS_MAX)) @@ -4132,12 +4124,12 @@ END SUBROUTINE REALLOCATE_BOUNDARY_ONE_D SUBROUTINE INITIALIZE_BOUNDARY_ONE_D(NM,OD_INDEX,SURF_INDEX) -USE GLOBAL_CONSTANTS, ONLY: RADIATION +USE GLOBAL_CONSTANTS, ONLY: RADIATION,NWP_MAX INTEGER, INTENT(IN) :: NM,OD_INDEX,SURF_INDEX TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D TYPE(SURFACE_TYPE), POINTER :: SF INTEGER :: NN,I,II -REAL(EB) :: RAMP_POSITION,TF,TB +REAL(EB) :: RAMP_POSITION,TF,TB,R(0:NWP_MAX) ONE_D => MESHES(NM)%BOUNDARY_ONE_D(OD_INDEX) SF => SURFACE(SURF_INDEX) @@ -4153,6 +4145,12 @@ SUBROUTINE INITIALIZE_BOUNDARY_ONE_D(NM,OD_INDEX,SURF_INDEX) DO I=1,MIN(ONE_D%N_CELLS_MAX,ONE_D%N_CELLS_INI) ONE_D%DX_OLD(I) = ONE_D%X(I)-ONE_D%X(I-1) ENDDO +R(0:ONE_D%N_CELLS_INI) = SF%THICKNESS + SF%INNER_RADIUS - ONE_D%X(0:ONE_D%N_CELLS_INI) +I = 0 +DO NN = 1,ONE_D%N_LAYERS + ONE_D%MIN_LAYER_MASS(NN) = SF%MIN_LAYER_MASS(NN)*(R(I)**SF%I_GRAD - R(I+ONE_D%N_LAYER_CELLS(NN))**SF%I_GRAD) + I = I + ONE_D%N_LAYER_CELLS(NN) +ENDDO IF (ALLOCATED(ONE_D%LAYER_THICKNESS_OLD)) ONE_D%LAYER_THICKNESS_OLD(1:ONE_D%N_LAYERS) = SF%LAYER_THICKNESS(1:SF%N_LAYERS) IF (SF%RAMP_T_I_INDEX > 0) THEN !NOTE: Replicating EVALUATE_RAMP since MODULE MATH_FUNCTIONS uses the MODULE which contains this routine @@ -4310,6 +4308,7 @@ SUBROUTINE PACK_BOUNDARY_PROP1(NM,IC,RC,LC,OS,B1_INDEX,UNPACK_IT,COUNT_ONLY,SURF RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RHO_G,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RDN,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%K_G,UNPACK_IT) +RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%M_DOT_LAYER_PP,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%Q_DOT_G_PP,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%Q_DOT_O2_PP,UNPACK_IT) RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%Q_CONDENSE,UNPACK_IT) @@ -5254,24 +5253,14 @@ END SUBROUTINE GET_WALL_NODE_WEIGHTS !> \param X_S_NEW Array of interior cell edge positions after shrinkage or swelling (m) !> \param INT_WGT Array of weighting factors for new arrangement of interior cells -SUBROUTINE GET_INTERPOLATION_WEIGHTS(GEOMETRY,NWP,NWP_NEW,INNER_RADIUS,X_S,X_S_NEW,INT_WGT) +SUBROUTINE GET_INTERPOLATION_WEIGHTS(GEOMETRY,I_GRAD,NWP,NWP_NEW,INNER_RADIUS,X_S,X_S_NEW,INT_WGT) -INTEGER, INTENT(IN) :: GEOMETRY,NWP,NWP_NEW +INTEGER, INTENT(IN) :: GEOMETRY,I_GRAD,NWP,NWP_NEW REAL(EB), INTENT(IN) :: X_S(0:NWP), X_S_NEW(0:NWP_NEW), INNER_RADIUS REAL(EB), INTENT(OUT) :: INT_WGT(NWP_NEW,NWP) REAL(EB) XUP,XLOW,XUP_NEW,XLOW_NEW,VOL_NEW,VOL,THICKNESS -INTEGER I_NEW, I_OLD, I_GRAD - - -SELECT CASE(GEOMETRY) - CASE(SURF_CARTESIAN) - I_GRAD = 1 - CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) - I_GRAD = 2 - CASE(SURF_SPHERICAL) - I_GRAD = 3 -END SELECT +INTEGER I_NEW, I_OLD I_OLD = 1 INT_WGT = 0._EB diff --git a/Source/init.f90 b/Source/init.f90 index 49b16a91ef4..2c0737921e9 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -4188,6 +4188,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW) DEALLOCATE(ONE_D%MATL_INDEX) ; ALLOCATE(ONE_D%MATL_INDEX(ONE_D%N_MATL)) DEALLOCATE(ONE_D%LAYER_THICKNESS) ; ALLOCATE(ONE_D%LAYER_THICKNESS(ONE_D%N_LAYERS)) DEALLOCATE(ONE_D%MIN_LAYER_THICKNESS) ; ALLOCATE(ONE_D%MIN_LAYER_THICKNESS(ONE_D%N_LAYERS)) + DEALLOCATE(ONE_D%MIN_LAYER_MASS) ; ALLOCATE(ONE_D%MIN_LAYER_MASS(ONE_D%N_LAYERS)) DEALLOCATE(ONE_D%HT3D_LAYER) ; ALLOCATE(ONE_D%HT3D_LAYER(ONE_D%N_LAYERS)) ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS) = LAYER_THICKNESS(1:ONE_D%N_LAYERS) ONE_D%MIN_LAYER_THICKNESS(1:ONE_D%N_LAYERS) = MIN_LAYER_THICKNESS(1:ONE_D%N_LAYERS) diff --git a/Source/read.f90 b/Source/read.f90 index 8503b306167..8f2fc10d1fc 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -7430,7 +7430,7 @@ SUBROUTINE PROC_MATL USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM,LINEAR_SYSTEM_SOLVE USE PHYSICAL_FUNCTIONS, ONLY: GET_TMP_REF -INTEGER :: I,N,NN,NL,NLPC,NS,NS2,NR,NR2,Z_INDEX(N_TRACKED_SPECIES,MAX_REACTIONS),IERR,ITMP,I_GRAD,& +INTEGER :: I,N,NN,NL,NLPC,NS,NS2,NR,NR2,Z_INDEX(N_TRACKED_SPECIES,MAX_REACTIONS),IERR,ITMP,& MATL_MATRIX_SIZE,REAC_COUNTER,TEMP_COUNTER,TEMP_MATL(N_MATL,MAX_REACTIONS) REAL(EB) :: ANS,H_ADJUST,NU_INERT,H_R_CALC(0:I_MAX_TEMP),SUM_NU(N_MATL,MAX_REACTIONS),DTMP,THICKNESS,VOL,X1 INTEGER, ALLOCATABLE, DIMENSION(:) :: MATL_MATRIX_POINTER @@ -7648,11 +7648,6 @@ SUBROUTINE PROC_MATL MATL_COEF_VECTOR(REAC_COUNTER) = MATL_COEF_VECTOR(REAC_COUNTER) + ML%NU_LPC(NLPC,NR)*ANS*ML%TMP_REF(NR) ELSE LIQUID_IF SF=>SURFACE(LPC%SURF_INDEX) - SELECT CASE(SF%GEOMETRY) - CASE(SURF_CARTESIAN) ; I_GRAD = 1 - CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; I_GRAD = 2 - CASE(SURF_SPHERICAL) ; I_GRAD = 3 - END SELECT ALLOCATE(RHO_H(SF%N_MATL)) RHO_H = 0._EB ALLOCATE(RHO(SF%N_MATL)) @@ -7662,7 +7657,7 @@ SUBROUTINE PROC_MATL THICKNESS = SUM(SF%LAYER_THICKNESS) X1 = THICKNESS+SF%INNER_RADIUS DO I=1,SF%N_LAYERS - VOL = X1**I_GRAD-(X1 - SF%LAYER_THICKNESS(I))**I_GRAD + VOL = X1**SF%I_GRAD-(X1 - SF%LAYER_THICKNESS(I))**SF%I_GRAD X1 = X1 - SF%LAYER_THICKNESS(I) MATL_RHO: DO NN=1,SF%N_MATL IF (SF%MATL_MASS_FRACTION(I,NN)<=TWO_EPSILON_EB) CYCLE MATL_RHO @@ -7815,7 +7810,7 @@ SUBROUTINE READ_SURF(QUICK_READ) 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,& - DELAMINATION_TMP(MAX_LAYERS),DELAMINATION_DENSITY(MAX_LAYERS) + DELAMINATION_TMP(MAX_LAYERS),DELAMINATION_DENSITY(MAX_LAYERS),MINIMUM_LAYER_MASS_FRACTION(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) @@ -7843,7 +7838,7 @@ SUBROUTINE READ_SURF(QUICK_READ) LEAK_PATH,LEAK_PATH_ID,LENGTH,MASS_FLUX,MASS_FLUX_TOTAL,MASS_FLUX_VAR,MASS_FRACTION,& MASS_TRANSFER_COEFFICIENT, & MATL_ID,MATL_MASS_FRACTION,MASS_PER_VOLUME,MCC_CONVERSION_FACTOR,MINIMUM_BURNOUT_TIME,& - MINIMUM_LAYER_THICKNESS,MLRPUA,MOISTURE_CONTENT,MOISTURE_FRACTION,& + MINIMUM_LAYER_MASS_FRACTION,MINIMUM_LAYER_THICKNESS,MLRPUA,MOISTURE_CONTENT,MOISTURE_FRACTION,& N_LAYER_CELLS_MAX,NEAR_WALL_EDDY_VISCOSITY,NEAR_WALL_TURBULENCE_MODEL,NET_HEAT_FLUX,& NO_SLIP,NODE_ID,NPPC,NUSSELT_C0,NUSSELT_C1,NUSSELT_C2,NUSSELT_M,& PARTICLE_EXTRACTION_VELOCITY,PARTICLE_MASS_FLUX,PARTICLE_SURFACE_DENSITY,PART_ID,& @@ -8387,19 +8382,23 @@ SUBROUTINE READ_SURF(QUICK_READ) SELECT CASE(GEOMETRY) CASE('CARTESIAN') SF%GEOMETRY = SURF_CARTESIAN + SF%I_GRAD = 1 IF (SF%WIDTH>0._EB) SF%BACKING = INSULATED CASE('CYLINDRICAL') SF%GEOMETRY = SURF_CYLINDRICAL IF (SF%INNER_RADIUS 0._EB) SF%LAYER_DENSITY(NL) = 1./SF%LAYER_DENSITY(NL) SF%THICKNESS = SF%THICKNESS + SF%LAYER_THICKNESS(NL) ENDDO COUNT_LAYERS - + SF%MIN_LAYER_MASS = SF%LAYER_DENSITY*MINIMUM_LAYER_MASS_FRACTION ! Define mass flux division point if the user does not specify. For all but ! surfaces with exposed backing, all pyrolyzed mass migrates to the front surface. @@ -9014,6 +9014,7 @@ SUBROUTINE SET_SURF_DEFAULTS MINIMUM_SCALING_HEAT_FLUX = 0._EB MAXIMUM_SCALING_HEAT_FLUX = 1.E9_EB MINIMUM_BURNOUT_TIME = 1.E6_EB +MINIMUM_LAYER_MASS_FRACTION = 1.E-12_EB MINIMUM_LAYER_THICKNESS = -1.E-6_EB ! The absolute value is the default, m MLRPUA = 0._EB MOISTURE_CONTENT = 0._EB @@ -9210,6 +9211,7 @@ SUBROUTINE PROC_SURF_1 CALL SHUTDOWN(MESSAGE) ; RETURN ENDIF ENDIF + CASE(SURF_SPHERICAL) END SELECT ENDIF ENDDO @@ -9251,7 +9253,7 @@ SUBROUTINE PROC_SURF_2 LOGICAL :: INDEX_LIST(MAX_LPC) REAL(EB) :: R_L(0:MAX_LAYERS), FUEL_MF,HRRPUA,E,MF_SUM,HR_SUM REAL(EB), ALLOCATABLE, DIMENSION(:) :: Q_REF_INT -INTEGER :: I_GRAD,I_CONE_RAMP,E_PTR +INTEGER :: I_CONE_RAMP,E_PTR INTEGER :: PROCESSED(N_RAMP) LOGICAL :: BURNING,BLOWING,SUCKING TYPE(RAMPS_TYPE),POINTER :: RP,RP_E2T,RP_INT,RP_QREF @@ -9265,12 +9267,6 @@ SUBROUTINE PROC_SURF_2 SF => SURFACE(N) IF (SF%THERMAL_BC_INDEX==THERMALLY_THICK) SF%INCLUDE_BOUNDARY_ONE_D_TYPE = .TRUE. - SELECT CASE(SF%GEOMETRY) - CASE(SURF_CARTESIAN) ; I_GRAD = 1 - CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; I_GRAD = 2 - CASE(SURF_SPHERICAL) ; I_GRAD = 3 - END SELECT - ! Particle Information SF%PART_INDEX = 0 @@ -9560,7 +9556,7 @@ SUBROUTINE PROC_SURF_2 DO NL=1,SF%N_LAYERS R_L(NL) = R_L(NL-1)-SF%LAYER_THICKNESS(NL) SF%SURFACE_DENSITY = SF%SURFACE_DENSITY + SF%LAYER_DENSITY(NL) * & - (R_L(NL-1)**I_GRAD-R_L(NL)**I_GRAD)/(REAL(I_GRAD,EB)*SF%THICKNESS**(I_GRAD-1)) + (R_L(NL-1)**SF%I_GRAD-R_L(NL)**SF%I_GRAD)/(REAL(SF%I_GRAD,EB)*SF%THICKNESS**(SF%I_GRAD-1)) ENDDO IF ((ABS(SF%SURFACE_DENSITY) <= TWO_EPSILON_EB) .AND. SF%BURN_AWAY) THEN diff --git a/Source/type.f90 b/Source/type.f90 index 0df02ff6556..4972317e29c 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -241,6 +241,7 @@ MODULE TYPES REAL(EB), ALLOCATABLE, DIMENSION(:) :: LAYER_THICKNESS !< (1:N_LAYERS) Thickness of layer (m) REAL(EB), ALLOCATABLE, DIMENSION(:) :: LAYER_THICKNESS_OLD !< (1:N_LAYERS) Thickness of layer (m) at last remesh REAL(EB), ALLOCATABLE, DIMENSION(:) :: MIN_LAYER_THICKNESS !< (1:N_LAYERS) Minimum thickness of layer (m) + REAL(EB), ALLOCATABLE, DIMENSION(:) :: MIN_LAYER_MASS !< (1:N_LAYERS) Layer mass for removal (m) REAL(EB), ALLOCATABLE, DIMENSION(:) :: MIN_DIFFUSIVITY !< (1:N_LAYERS) Min diffusivity of all matls in layer (m2/s) REAL(EB), ALLOCATABLE, DIMENSION(:) :: DDSUM !< Scaling factor to get minimum cell size REAL(EB), ALLOCATABLE, DIMENSION(:) :: SMALLEST_CELL_SIZE !< Minimum cell size (m) @@ -321,6 +322,7 @@ MODULE TYPES REAL(EB) :: RHO_G !< Gas density in near wall cell (kg/m3) REAL(EB) :: RDN=1._EB !< \f$ 1/ \delta n \f$ at the surface (1/m) REAL(EB) :: K_G=0.025_EB !< Thermal conductivity of gas in adjacent gas phase cell near wall + REAL(EB) :: M_DOT_LAYER_PP=0._EB !< Mass of a layer removed in kg used to update OB%MASS REAL(EB) :: Q_DOT_G_PP=0._EB !< Heat release rate per unit area (W/m2) REAL(EB) :: Q_DOT_O2_PP=0._EB !< Heat release rate per unit area (W/m2) due to oxygen consumption REAL(EB) :: Q_CONDENSE=0._EB !< Heat release rate per unit area (W/m2) due to gas condensation @@ -922,6 +924,7 @@ MODULE TYPES REAL(EB), ALLOCATABLE, DIMENSION(:) :: MASS_FRACTION,MASS_FLUX,DDSUM,SMALLEST_CELL_SIZE,SWELL_RATIO REAL(EB), ALLOCATABLE, DIMENSION(:) :: REFERENCE_HEAT_FLUX !< Reference flux for the flux scaling pyrolysis model (kW/m2) REAL(EB), ALLOCATABLE, DIMENSION(:) :: SPYRO_TH_FACTOR !< Thickness scaling factor for the flux scaling pyrolysis model + REAL(EB), ALLOCATABLE, DIMENSION(:) :: MIN_LAYER_MASS !< Initial mass in each layer TYPE(RAMP_ID_TYPE), ALLOCATABLE, DIMENSION(:) :: RAMP INTEGER, DIMENSION(3) :: RGB @@ -938,6 +941,7 @@ MODULE TYPES INTEGER :: N_CONE_CURVES=0 !< Total number of time vs. heat release rate curves specified for the S-pyro model INTEGER :: N_MATL !< Total number of unique materials over all layers INTEGER :: NODE_INDEX=0 !< Ductnode index + INTEGER :: I_GRAD = 1 !< Geometry factor for volume INTEGER, DIMENSION(30) :: ONE_D_REALS_ARRAY_SIZE=0,ONE_D_INTEGERS_ARRAY_SIZE=0,ONE_D_LOGICALS_ARRAY_SIZE=0 INTEGER, ALLOCATABLE, DIMENSION(:) :: N_LAYER_CELLS !< Number of wall cells per material layer INTEGER, ALLOCATABLE, DIMENSION(:) :: LAYER_INDEX !< The layer each wall cell belongs to diff --git a/Source/wall.f90 b/Source/wall.f90 index f87bd022447..5ff406bf216 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -1403,13 +1403,13 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) LL = 2*OMESH(EWC%NOM)%N_EXTERNAL_OBST OMESH(EWC%NOM)%REAL_SEND_PKG8(LL-1) = REAL(OTHER_MESH_OBST_INDEX,EB) OMESH(EWC%NOM)%REAL_SEND_PKG8(LL) = & - (B1%M_DOT_PART_ACTUAL+SUM(B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES)))*DT*B1%AREA + (B1%M_DOT_PART_ACTUAL+B1%M_DOT_LAYER_PP+SUM(B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES)))*DT*B1%AREA ENDIF !$OMP END CRITICAL ELSE !$OMP CRITICAL IF (OBST_INDEX>0) OBSTRUCTION(OBST_INDEX)%MASS = OBSTRUCTION(OBST_INDEX)%MASS - & - (B1%M_DOT_PART_ACTUAL+SUM(B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES)))*DT*B1%AREA + (B1%M_DOT_PART_ACTUAL+B1%M_DOT_LAYER_PP+SUM(B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES)))*DT*B1%AREA !$OMP END CRITICAL ENDIF ENDIF CONSUME_MASS @@ -1832,7 +1832,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, Q_LIQUID_F,Q_LIQUID_B,LAYER_DIVIDE,TMP_BACK,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,& - RADIUS,HTC_LIMIT,CP1,CP2,DENOM,THICKNESS,DT_FO,DDSUM,NODE_RDT(NWP_MAX) + RADIUS,HTC_LIMIT,CP1,CP2,DENOM,THICKNESS,DT_FO,DDSUM,NODE_RDT(NWP_MAX),RMR 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,T_BOIL_EFF REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART_S,M_DOT_PART_S @@ -1840,16 +1840,17 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: RHO_DOT,INT_WGT 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(NWP_MAX) :: AAS,BBS,CCS,DDS,DDT,Q_S,Q_IR,Q_ADD,TWO_DX_KAPPA_S,DX_S,MF_FRAC,REGRID_FACTOR,DX_S_NEW +REAL(EB), DIMENSION(0:NWP_MAX+1) :: RHO_S,DELTA_TMP,RDX_S,TMP_NEW REAL(EB), DIMENSION(0:NWP_MAX) :: X_S_NEW,RDXN_S,R_S,R_S_NEW,DX_WGT_S +REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_MASS INTEGER, DIMENSION(0:NWP_MAX+1) :: LAYER_INDEX 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 +INTEGER :: NWP_NEW,IZERO,SURF_INDEX,BACKING,NWP,I,NL,N,NN,OBST_INDEX,& + N_CELLS,ITMP,ITER,BACK_MESH,BACK_INDEX,BACK_WALL_INDEX,FIRST_CELL,LAST_CELL CHARACTER(MESSAGE_LENGTH) :: MESSAGE 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) + CELL_ZERO(MAX_LAYERS_HT3D),TMP_CHECK(MAX_LAYERS_HT3D),CELL_ZERO_CELL(NWP_MAX) TYPE(WALL_TYPE), POINTER :: WC,WC_BACK TYPE(THIN_WALL_TYPE), POINTER :: TW,TW_BACK TYPE(CFACE_TYPE), POINTER :: CFA,CFA_BACK @@ -2002,14 +2003,6 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, RETURN ENDIF -! Exponents for cylindrical or spherical coordinates - -SELECT CASE(SF%GEOMETRY) - CASE(SURF_CARTESIAN) ; I_GRAD = 1 - CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; I_GRAD = 2 - CASE(SURF_SPHERICAL) ; I_GRAD = 3 -END SELECT - ! Create array of material densities that are adjusted in special cases DO NL=1,ONE_D%N_LAYERS @@ -2035,6 +2028,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, Q_DOT_G_PP_NET = 0._EB Q_DOT_O2_PP_NET = 0._EB B1%M_DOT_PART_ACTUAL = 0._EB + B1%M_DOT_LAYER_PP = 0._EB ENDIF ! Start time iterations here @@ -2053,34 +2047,6 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, 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 - ONE_D%RHO_C_S = 0._EB - 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) - ENDDO MATERIAL_LOOP0 - 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 - 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 @@ -2103,6 +2069,34 @@ 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 + ! 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 + ONE_D%RHO_C_S = 0._EB + 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) + ENDDO MATERIAL_LOOP0 + 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 + 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 + ! Calculate minimum DT based on Fourier number, Fo IF (CHECK_FO) THEN @@ -2317,9 +2311,14 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, Q_NET_F = B1%Q_RAD_IN - B1%EMISSIVITY *SIGMA*B1%TMP_F**4 + Q_CON_F Q_NET_B = Q_RAD_IN_B - EMISSIVITY_BACK*SIGMA*B1%TMP_B**4 + Q_CON_B ENDIF - DO I=1,NWP - DELTA_TMP(I) = DT_BC*Q_S(I)/ONE_D%RHO_C_S(I) + DO I=2,NWP-1 + DELTA_TMP(I) = (DT_BC/ONE_D%RHO_C_S(I))*(RDX_S(I)*(ONE_D%K_S(I)* RDXN_S(I)* (ONE_D%TMP(I+1)-ONE_D%TMP(I))-& + ONE_D%K_S(I-1)*RDXN_S(I-1)*(ONE_D%TMP(I)- ONE_D%TMP(I-1))) + Q_S(I)) ENDDO + DELTA_TMP(1) = (DT_BC/ONE_D%RHO_C_S(1))*& + (RDX_S(1)*(ONE_D%K_S(1)*RDXN_S(1)*(ONE_D%TMP(2)-ONE_D%TMP(1))+Q_NET_F) + Q_S(1)) + DELTA_TMP(NWP) = (DT_BC/ONE_D%RHO_C_S(NWP))*& + (RDX_S(NWP)*(-Q_NET_B-ONE_D%K_S(NWP-1)*RDXN_S(NWP-1)*(ONE_D%TMP(NWP)-ONE_D%TMP(NWP-1))) + Q_S(NWP)) 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) @@ -2360,12 +2359,14 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, REMESH_LAYER = .FALSE. TMP_CHECK = .FALSE. + LAYER_MASS = HUGE(1._EB) PYROLYSIS_PREDICTED_IF_2: IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN + LAYER_MASS = 0._EB ! Convert Q_S to kW DO I=1,NWP - Q_S(I) = Q_S(I)*ABS(R_S(I-1)**I_GRAD-R_S(I)**I_GRAD) + Q_S(I) = Q_S(I)*ABS(R_S(I-1)**SF%I_GRAD-R_S(I)**SF%I_GRAD) ENDDO POINT_LOOP2: DO I=1,NWP @@ -2376,6 +2377,10 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, IF (ALL(ABS(RHO_DOT(:,I)) TWO_EPSILON_EB) NODE_RDT(I) = MIN(NODE_RDT(I),MATERIAL(ONE_D%MATL_INDEX(N))%RENODE_DELTA_T) REGRID_FACTOR(I) = REGRID_FACTOR(I) + ONE_D%MATL_COMP(N)%RHO(I)/RHO_ADJUSTED(LAYER_INDEX(I),N) ENDDO MATERIAL_LOOP1a @@ -2421,96 +2428,69 @@ 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%DELAMINATION_TMP(NL)) .OR.& - ( RHO_S(I+ONE_D%N_LAYER_CELLS(NL)) < SF%DELAMINATION_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 R_S_NEW(NWP) = 0._EB DO I=NWP-1,0,-1 - R_S_NEW(I) = ( R_S_NEW(I+1)**I_GRAD + (R_S(I)**I_GRAD-R_S(I+1)**I_GRAD)*REGRID_FACTOR(I+1) )**(1./REAL(I_GRAD,EB)) + R_S_NEW(I) = ( R_S_NEW(I+1)**SF%I_GRAD + (R_S(I)**SF%I_GRAD-R_S(I+1)**SF%I_GRAD)*REGRID_FACTOR(I+1) )**& + (1./REAL(SF%I_GRAD,EB)) ENDDO X_S_NEW(0) = 0._EB CELL_ZERO = .FALSE. + CELL_ZERO_CELL = .FALSE. - DO I=1,NWP + RDT_CHECK: DO I=1,NWP X_S_NEW(I) = R_S_NEW(0) - R_S_NEW(I) ! If Cell disappears we must remesh IF ((X_S_NEW(I)-X_S_NEW(I-1)) < TWO_EPSILON_EB) THEN - REMESH_LAYER(LAYER_INDEX(I)) = .TRUE. + IF (ONE_D%N_LAYER_CELLS(LAYER_INDEX(I))>1) REMESH_LAYER(LAYER_INDEX(I)) = .TRUE. CELL_ZERO(LAYER_INDEX(I)) = .TRUE. + CELL_ZERO_CELL(I) = .TRUE. ELSE + IF (ONE_D%N_LAYER_CELLS(LAYER_INDEX(I))==1) CYCLE RDT_CHECK ! If cell size changes enough compared to prior remseh size, remesh - IF (ABS((X_S_NEW(I)-X_S_NEW(I-1))/ONE_D%DX_OLD(I)-1._EB) > SF%REMESH_RATIO) REMESH_LAYER(LAYER_INDEX(I)) = .TRUE. - IF (I > 1 .AND. REMESH_LAYER(LAYER_INDEX(I))) THEN - IF (LAYER_INDEX(I-1)==LAYER_INDEX(I) .AND. ABS(ONE_D%TMP(I)-ONE_D%TMP(I-1)) > NODE_RDT(I)) THEN - IF (.NOT. CELL_ZERO(LAYER_INDEX(I))) REMESH_LAYER(LAYER_INDEX(I)) = .FALSE. - TMP_CHECK(LAYER_INDEX(I)) = .TRUE. - ENDIF - ENDIF - IF (I < NWP .AND. REMESH_LAYER(LAYER_INDEX(I))) THEN - IF (LAYER_INDEX(I+1)==LAYER_INDEX(I) .AND. ABS(ONE_D%TMP(I+1)-ONE_D%TMP(I)) > NODE_RDT(I)) THEN - IF (.NOT. CELL_ZERO(LAYER_INDEX(I))) REMESH_LAYER(LAYER_INDEX(I)) = .FALSE. - TMP_CHECK(LAYER_INDEX(I)) = .TRUE. - ENDIF - ENDIF - ENDIF - ENDDO - - ! Check for layers that are too small, layer thickness dropping too much, or couldn't drop enough nodes last remesh - I = 0 - DO NL=1,ONE_D%N_LAYERS - IF (ONE_D%N_LAYER_CELLS(NL) == 0) CYCLE - ONE_D%LAYER_THICKNESS(NL) = X_S_NEW(I+ONE_D%N_LAYER_CELLS(NL)) - X_S_NEW(I) - IF (ONE_D%LAYER_THICKNESS(NL) < 0.1_EB*ONE_D%MIN_LAYER_THICKNESS(NL)) THEN - REMESH_LAYER(NL) = .TRUE. - ELSE - IF (.NOT. TMP_CHECK(NL) .AND. ONE_D%LAYER_THICKNESS_OLD(NL)-ONE_D%LAYER_THICKNESS(NL) > & - 1.5_EB*MINVAL(ONE_D%DX_OLD(I+1:I+ONE_D%N_LAYER_CELLS(NL)))) THEN - REMESH_LAYER(NL) = .TRUE. - ELSE - IF (.NOT. TMP_CHECK(NL) .AND. ONE_D%REMESH_NWP(NL) > ONE_D%N_LAYER_CELLS(NL)) REMESH_LAYER(NL) = .TRUE. + RMR = (X_S_NEW(I)-X_S_NEW(I-1))/ONE_D%DX_OLD(I) + IF (ABS(RMR - 1._EB) > TWO_EPSILON_EB) THEN + IF (ABS(RMR-1._EB) > SF%REMESH_RATIO) THEN + REMESH_LAYER(LAYER_INDEX(I)) = .TRUE. + ELSEIF (RMR < 1._EB) THEN + IF (I>1 .AND. LAYER_INDEX(I-1)==LAYER_INDEX(I)) THEN + IF (ABS(ONE_D%TMP(I-1)*(1._EB-RMR)+ONE_D%TMP(I)*RMR - ONE_D%TMP(I)) > NODE_RDT(I)) THEN + REMESH_LAYER(LAYER_INDEX(I)) = .TRUE. + TMP_CHECK(LAYER_INDEX(I)) = .TRUE. + ENDIF + ENDIF + IF (I NODE_RDT(I)) THEN + REMESH_LAYER(LAYER_INDEX(I)) = .TRUE. + TMP_CHECK(LAYER_INDEX(I)) = .TRUE. + ENDIF + ENDIF + ENDIF ENDIF ENDIF - I = I + ONE_D%N_LAYER_CELLS(NL) - ENDDO + ENDDO RDT_CHECK ! 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 - IF (X_S_NEW(1)-X_S_NEW(0) < TWO_EPSILON_EB) Q_S(2) = Q_S(2) + Q_S(1) - IF (X_S_NEW(NWP)-X_S_NEW(NWP-1) < TWO_EPSILON_EB) Q_S(NWP-1) = Q_S(NWP-1) + Q_S(NWP) + IF (CELL_ZERO_CELL(1)) Q_S(2) = Q_S(2) + Q_S(1) + IF (CELL_ZERO_CELL(NWP)) Q_S(NWP-1) = Q_S(NWP-1) + Q_S(NWP) DO I=2,NWP-1 - IF (X_S_NEW(I) - X_S_NEW(I-1) < TWO_EPSILON_EB) THEN + IF (CELL_ZERO_CELL(I)) THEN N = 0 - IF (X_S_NEW(I-1) - X_S_NEW(I-2) > TWO_EPSILON_EB) N=N+1 - IF (X_S_NEW(I+1) - X_S_NEW(I) > TWO_EPSILON_EB) N=N+2 + IF (CELL_ZERO_CELL(I-1)) N=N+1 + IF (CELL_ZERO(I+1)) N=N+2 SELECT CASE (N) CASE(1) Q_S(I-1) = Q_S(I-1) + Q_S(I) CASE(2) Q_S(I+1) = Q_S(I+1) + Q_S(I) CASE(3) - VOL = (R_S_NEW(I-1)**I_GRAD-R_S_NEW(I)**I_GRAD) / & - ((R_S_NEW(I-1)**I_GRAD-R_S_NEW(I)**I_GRAD)+(R_S_NEW(I)**I_GRAD-R_S_NEW(I+1)**I_GRAD)) + VOL = (R_S_NEW(I-1)**SF%I_GRAD-R_S_NEW(I)**SF%I_GRAD) / & + ((R_S_NEW(I-1)**SF%I_GRAD-R_S_NEW(I)**SF%I_GRAD)+(R_S_NEW(I)**SF%I_GRAD-R_S_NEW(I+1)**SF%I_GRAD)) Q_S(I-1) = Q_S(I-1) + Q_S(I) * VOL Q_S(I+1) = Q_S(I+1) + Q_S(I) * (1._EB-VOL) END SELECT @@ -2518,7 +2498,88 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, ENDDO ENDIF - REMESH_CHECK = ANY(ABS(REGRID_FACTOR(1:NWP)-1._EB)>TWO_EPSILON_EB) + ! Check for layers that are too small, layer thickness dropping too much, or couldn't drop enough nodes last remesh + I = 0 + CELL_ZERO_CELL = .FALSE. + LAYER_THICKNESS_CHECK: DO NL=1,ONE_D%N_LAYERS + ONE_D%LAYER_THICKNESS(NL) = X_S_NEW(I+ONE_D%N_LAYER_CELLS(NL)) - X_S_NEW(I) + IF (ONE_D%N_LAYER_CELLS(NL) == 0) CYCLE + IF (ONE_D%LAYER_THICKNESS(NL) <= ONE_D%MIN_LAYER_THICKNESS(NL) .OR. LAYER_MASS(NL) <= ONE_D%MIN_LAYER_MASS(NL)) THEN + CELL_ZERO_CELL(I+1:I+ONE_D%N_LAYER_CELLS(NL)) = .TRUE. + ONE_D%N_LAYER_CELLS(NL)=0 + ONE_D%LAYER_THICKNESS(NL) = 0._EB + REMESH_LAYER(NL) = .FALSE. + REMESH_CHECK = .TRUE. + B1%M_DOT_LAYER_PP = B1%M_DOT_LAYER_PP + LAYER_MASS(NL) + CYCLE LAYER_THICKNESS_CHECK + ENDIF + + ! Delamination layer fall-off + + IF (ONE_D%TMP(I+ONE_D%N_LAYER_CELLS(NL)) > SF%DELAMINATION_TMP(NL) .OR. & + RHO_S(I+ONE_D%N_LAYER_CELLS(NL)) < SF%DELAMINATION_DENSITY(NL)) THEN + CELL_ZERO_CELL(I+1:I+ONE_D%N_LAYER_CELLS(NL)) = .TRUE. + ONE_D%N_LAYER_CELLS(NL)=0 + ONE_D%LAYER_THICKNESS(NL) = 0._EB + REMESH_LAYER(NL) = .FALSE. + REMESH_CHECK = .TRUE. + B1%M_DOT_LAYER_PP = B1%M_DOT_LAYER_PP + LAYER_MASS(NL) + CYCLE LAYER_THICKNESS_CHECK + ENDIF + + ! Overall layer thickness drops significantly compared to smallest cell or prior renoding stopped at one + ! wall cell removal for the layer + IF (.NOT. TMP_CHECK(NL) .AND. ONE_D%LAYER_THICKNESS_OLD(NL)-ONE_D%LAYER_THICKNESS(NL) > & + 1.5_EB*MINVAL(ONE_D%DX_OLD(I+1:I+ONE_D%N_LAYER_CELLS(NL)))) THEN + REMESH_LAYER(NL) = .TRUE. + ELSE + IF (.NOT. TMP_CHECK(NL) .AND. ONE_D%REMESH_NWP(NL) > ONE_D%N_LAYER_CELLS(NL)) REMESH_LAYER(NL) = .TRUE. + ENDIF + I = I + ONE_D%N_LAYER_CELLS(NL) + ENDDO LAYER_THICKNESS_CHECK + + LAYER_REMOVE: IF (ANY(CELL_ZERO_CELL)) THEN + N = 0 + FIRST_CELL = -1 + X_S_NEW(0) = 0._EB + DO I=1,NWP + IF (CELL_ZERO_CELL(I)) CYCLE + IF (FIRST_CELL < 0) FIRST_CELL = I + LAST_CELL = I + N = N + 1 + DX_S_NEW(N) = DX_S(I) + X_S_NEW(N) = X_S_NEW(N-1) + DX_S(I) + TMP_NEW(N) = ONE_D%TMP(I) + Q_S(N) = Q_S(I) + DO NN=1,ONE_D%N_MATL + ONE_D%MATL_COMP(NN)%RHO(N) = ONE_D%MATL_COMP(NN)%RHO(I) + ENDDO + ENDDO + IF (N==0) EXIT LAYER_REMOVE + IF (FIRST_CELL > 1) THEN + B1%TMP_F = ONE_D%TMP(FIRST_CELL-1) + DX_S(FIRST_CELL-1)/(DX_S(FIRST_CELL-1)+DX_S(FIRST_CELL)) * & + (ONE_D%TMP(FIRST_CELL)-ONE_D%TMP(FIRST_CELL-1)) + ONE_D%TMP(0) = 2*B1%TMP_F - ONE_D%TMP(FIRST_CELL) + ENDIF + IF (LAST_CELL < NWP) THEN + B1%TMP_F = ONE_D%TMP(LAST_CELL) + DX_S(LAST_CELL)/(DX_S(LAST_CELL)+DX_S(LAST_CELL+1)) * & + (ONE_D%TMP(LAST_CELL+1)-ONE_D%TMP(LAST_CELL)) + ONE_D%TMP(N+1) = 2*B1%TMP_B - ONE_D%TMP(LAST_CELL) + ENDIF + ONE_D%X(1:N) = X_S_NEW(1:N) + ONE_D%X(N+1:NWP) = 0._EB + DX_S(1:N) = DX_S_NEW(1:N) + DX_S(N+1:NWP) = 0._EB + ONE_D%TMP(1:N) = TMP_NEW(1:N) + ONE_D%TMP(N+2:NWP) = 0._EB + Q_S(N+1:NWP) = 0._EB + DO I=NN,ONE_D%N_MATL + ONE_D%MATL_COMP(NN)%RHO(N+1:NWP) = 0._EB + ENDDO + NWP = N + ENDIF LAYER_REMOVE + + REMESH_CHECK = REMESH_CHECK .OR. ANY(ABS(REGRID_FACTOR(1:NWP)-1._EB)>TWO_EPSILON_EB) ! Some node changes size but no layer trips remesh check. Just redo node weight with X_S_NEW @@ -2528,7 +2589,6 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, 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) ENDIF - N_LAYER_CELLS_NEW(1:ONE_D%N_LAYERS) = ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS) ! Some layer needs to be checked for a remesh. @@ -2539,20 +2599,10 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, I = 0 DX_MIN = 0._EB ONE_D%X(0:NWP) = X_S_NEW(0:NWP) - - LAYER_LOOP: DO NL=1,ONE_D%N_LAYERS - ! Layer is too small. Delete it and shift any following nodes up in x-distance - - IF (ONE_D%LAYER_THICKNESS(NL) < 0.1_EB*ONE_D%MIN_LAYER_THICKNESS(NL)) THEN - N_LAYER_CELLS_NEW(NL) = 0 - ONE_D%X(I+ONE_D%N_LAYER_CELLS(NL):NWP) = ONE_D%X(I+ONE_D%N_LAYER_CELLS(NL):NWP)-ONE_D%LAYER_THICKNESS(NL) - ONE_D%DX_OLD(I+1:I+ONE_D%N_LAYER_CELLS(NL)) = 0._EB - ONE_D%LAYER_THICKNESS(NL) = 0._EB - ONE_D%LAYER_THICKNESS_OLD(NL) = 0._EB - ONE_D%REMESH_NWP(NL) = 0 - CYCLE LAYER_LOOP - ENDIF + LAYER_LOOP: DO NL=1,ONE_D%N_LAYERS + + IF (ONE_D%N_LAYER_CELLS(NL) == 0) CYCLE LAYER_LOOP ! Layer is one cell. Only invoke remeshing if the layer grows in size which could add a wall node @@ -2573,6 +2623,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, IF (.NOT. REMESH_LAYER(NL)) THEN N_LAYER_CELLS_NEW(NL) = ONE_D%N_LAYER_CELLS(NL) + ONE_D%REMESH_NWP(NL) = N_LAYER_CELLS_NEW(NL) NWP_NEW = NWP_NEW + N_LAYER_CELLS_NEW(NL) DO N = I,I+ONE_D%N_LAYER_CELLS(NL)-1 IF (ONE_D%X(I+1)-ONE_D%X(I) < ONE_D%SMALLEST_CELL_SIZE(NL)) ONE_D%SMALLEST_CELL_SIZE(NL) = ONE_D%X(I+1)-ONE_D%X(I) @@ -2622,7 +2673,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, THICKNESS = SUM(ONE_D%LAYER_THICKNESS) ONE_D%X(0:NWP) = X_S_NEW(0:NWP) - + ENDIF REMESH_LAYER_1 ! Shrinking wall has gone to zero thickness. @@ -2644,9 +2695,9 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, ENDIF EXIT SUB_TIMESTEP_LOOP ENDIF - + ! Re-generate grid for a wall changing thickness - + REMESH_LAYER_2: IF (ANY(REMESH_LAYER(1:ONE_D%N_LAYERS))) THEN RHO_H_S = 0._EB TMP_S = 0._EB @@ -2654,7 +2705,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, ! Store wall enthalpy for later temperature extraction. DO I=1,NWP - VOL = (THICKNESS+SF%INNER_RADIUS-ONE_D%X(I-1))**I_GRAD-(THICKNESS+SF%INNER_RADIUS-ONE_D%X(I))**I_GRAD + VOL = (THICKNESS+SF%INNER_RADIUS-ONE_D%X(I-1))**SF%I_GRAD-(THICKNESS+SF%INNER_RADIUS-ONE_D%X(I))**SF%I_GRAD MATL_REMESH: DO N=1,ONE_D%N_MATL IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE MATL_REMESH ML => MATERIAL(ONE_D%MATL_INDEX(N)) @@ -2683,13 +2734,14 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, ! Interpolate densities and temperature from old grid to new grid ALLOCATE(INT_WGT(NWP_NEW,NWP),STAT=IZERO) - CALL GET_INTERPOLATION_WEIGHTS(SF%GEOMETRY,NWP,NWP_NEW,SF%INNER_RADIUS,ONE_D%X(0:NWP),X_S_NEW(0:NWP_NEW),INT_WGT) + CALL GET_INTERPOLATION_WEIGHTS(SF%GEOMETRY,SF%I_GRAD,NWP,NWP_NEW,SF%INNER_RADIUS,ONE_D%X(0:NWP),& + X_S_NEW(0:NWP_NEW),INT_WGT) N_CELLS = MAX(NWP,NWP_NEW) 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 + VOL = (THICKNESS+SF%INNER_RADIUS-X_S_NEW(I-1))**SF%I_GRAD-(THICKNESS+SF%INNER_RADIUS-X_S_NEW(I))**SF%I_GRAD TMP_S(I) = TMP_S(I) / VOL ENDDO @@ -2736,8 +2788,9 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, T_NODE = ONE_D%TMP(I) ENDDO T_SEARCH DO N=1,ONE_D%N_MATL - ONE_D%MATL_COMP(N)%RHO(I) = ONE_D%MATL_COMP(N)%RHO(I) /& - ((SF%INNER_RADIUS+X_S_NEW(NWP_NEW)-X_S_NEW(I-1))**I_GRAD-(SF%INNER_RADIUS+X_S_NEW(NWP_NEW)-X_S_NEW(I))**I_GRAD) + ONE_D%MATL_COMP(N)%RHO(I) = ONE_D%MATL_COMP(N)%RHO(I) / & + ((SF%INNER_RADIUS+X_S_NEW(NWP_NEW)-X_S_NEW(I-1))**SF%I_GRAD - & + (SF%INNER_RADIUS+X_S_NEW(NWP_NEW)-X_S_NEW(I))**SF%I_GRAD) ENDDO ENDDO @@ -2752,7 +2805,8 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, ! Convert Q_S back to kW/m^3 DO I=1,NWP - Q_S(I) = Q_S(I)/((SF%INNER_RADIUS+ONE_D%X(NWP)-ONE_D%X(I-1))**I_GRAD-(SF%INNER_RADIUS+ONE_D%X(NWP)-ONE_D%X(I))**I_GRAD) + Q_S(I) = Q_S(I)/((SF%INNER_RADIUS+ONE_D%X(NWP)-ONE_D%X(I-1))**SF%I_GRAD - & + (SF%INNER_RADIUS+ONE_D%X(NWP)-ONE_D%X(I))**SF%I_GRAD) ENDDO ENDIF PYROLYSIS_PREDICTED_IF_2 @@ -2875,6 +2929,11 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, B1%TMP_F = MIN(TMPMAX,MAX(TMPMIN,B1%TMP_F)) B1%TMP_B = MIN(TMPMAX,MAX(TMPMIN,B1%TMP_B)) + ! If clipped make ghost consistent with new value + + IF (B1%TMP_F==TMPMIN .OR. B1%TMP_F==TMPMAX) ONE_D%TMP(0) = 2._EB*B1%TMP_F - ONE_D%TMP(1) + IF (B1%TMP_B==TMPMIN .OR. B1%TMP_B==TMPMAX) ONE_D%TMP(NWP+1) = 2._EB*B1%TMP_B - ONE_D%TMP(NWP) + ! Updated particle production IF (ONE_D%N_LPC > 0) THEN @@ -2891,6 +2950,8 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, B1%N_SUBSTEPS = B1%N_SUBSTEPS + 1 ENDDO SUB_TIMESTEP_LOOP +IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) & + B1%M_DOT_LAYER_PP = B1%M_DOT_LAYER_PP / (SF%I_GRAD*(SF%THICKNESS+SF%INNER_RADIUS)**(SF%I_GRAD-1)*DT_BC) IF (ALLOCATED(B1%M_DOT_G_PP_ACTUAL) .AND. ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN B1%M_DOT_G_PP_ADJUST = M_DOT_G_PP_ADJUST_NET @@ -2999,7 +3060,7 @@ SUBROUTINE PERFORM_PYROLYSIS RHO_TEMP(N) = ONE_D%MATL_COMP(N)%RHO(I) ENDDO - IF (ONE_D%LAYER_THICKNESS(LAYER_INDEX(I))