diff --git a/Source/dump.f90 b/Source/dump.f90 index e8feaab6462..1f68921e584 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -8523,7 +8523,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_IN INTEGER, INTENT(IN), OPTIONAL :: OPT_WALL_INDEX,OPT_LP_INDEX,OPT_CFACE_INDEX,OPT_BNDF_INDEX,OPT_DEVC_INDEX,OPT_CUT_FACE_INDEX,& OPT_NODE_INDEX,OPT_PROF_INDEX INTEGER, INTENT(IN) :: INDX,Y_INDEX,Z_INDEX,PART_INDEX -REAL(EB) :: Q_CON,RHOSUM,VOLSUM,ZZ_GET(1:N_TRACKED_SPECIES),Y_SPECIES,DEPTH,CHAR_FRONT,UN,H_S,RHO_D_DYDN,U_CELL,V_CELL,W_CELL,& +REAL(EB) :: Q_CON,RHOSUM,VOLSUM,ZZ_GET(1:N_TRACKED_SPECIES),Y_SPECIES,DEPTH,ASH_DEPTH,UN,H_S,RHO_D_DYDN,U_CELL,V_CELL,W_CELL,& LTMP,ATMP,CTMP,H_W_EFF,X0,VOL,DN,PRESS,& NVEC(3),PVEC(3),TAU_IJ(3,3),VEL_CELL(3),VEL_WALL(3),MU_WALL,RHO_WALL,FVEC(3),SVEC(3),TVEC1(3),TVEC2(3),& PR1,PR2,Z1,Z2,RADIUS,CUT_FACE_AREA,SOLID_PHASE_OUTPUT_CTF,AAA,BBB,CCC,ALP,BET,GAM,MMM,DTMP @@ -9278,9 +9278,9 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_IN ! for the moment this assumes there is only one char reaction IF (ML%N_O2(1)>0._EB) THEN DEPTH = 0.5_EB*(ONE_D%X(I_DEPTH-1)+ONE_D%X(I_DEPTH)) - CHAR_FRONT = 0._EB - IF (TEST_NEW_CHAR_MODEL) CHAR_FRONT = ONE_D%X(B2%I_CHAR_FRONT-1) - SOLID_PHASE_OUTPUT = B2%Y_O2_F*EXP(-MAX(0._EB,DEPTH-CHAR_FRONT)/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(1))) + ASH_DEPTH = 0._EB + IF (TEST_NEW_CHAR_MODEL) ASH_DEPTH = ONE_D%X(B2%I_ASH_DEPTH-1) + SOLID_PHASE_OUTPUT = B2%Y_O2_F*EXP(-MAX(0._EB,DEPTH-ASH_DEPTH)/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(1))) ENDIF ENDIF CASE(82) ! BLOWING CORRECTION diff --git a/Source/type.f90 b/Source/type.f90 index a5361d84e2e..d5144a86162 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -363,7 +363,7 @@ MODULE TYPES INTEGER :: SURF_INDEX=-1 !< Surface index INTEGER :: HEAT_TRANSFER_REGIME=0 !< 1=Forced convection, 2=Natural convection, 3=Impact convection, 4=Resolved INTEGER :: Y_O2_ITER=0 !< Number of iterations for surface O2 solve - INTEGER :: I_CHAR_FRONT=1 !< Index of solid cell where char formation starts + INTEGER :: I_ASH_DEPTH=1 !< Index of solid cell where char formation starts END TYPE BOUNDARY_PROP2_TYPE diff --git a/Source/wall.f90 b/Source/wall.f90 index bf10653df35..03cdebe6df2 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -2883,7 +2883,7 @@ SUBROUTINE PERFORM_PYROLYSIS USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION,GET_SPECIFIC_HEAT REAL(EB), DIMENSION(N_TRACKED_SPECIES) :: M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,ZZ_GET REAL(EB), DIMENSION(MAX_MATERIALS) :: RHO_TEMP,M_DOT_S_PPP -REAL(EB) :: Q_DOT_G_PPP,Q_DOT_O2_PPP,CP_FILM,TMP_FILM,H_MASS,CHAR_FRONT,& +REAL(EB) :: Q_DOT_G_PPP,Q_DOT_O2_PPP,CP_FILM,TMP_FILM,H_MASS,ASH_DEPTH,& Y_O2_G,Y_O2_F,M_DOT_O2_PP,Y_LOWER,Y_UPPER,M_DOT_ERROR,M_DOT_ERROR_OLD,Y_O2_F_OLD,DY,DE REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART,M_DOT_PART INTEGER :: ITER,MAX_ITER @@ -2901,17 +2901,17 @@ SUBROUTINE PERFORM_PYROLYSIS ! Determine char front position -CHAR_FRONT=0._EB +ASH_DEPTH=0._EB IF (TEST_NEW_CHAR_MODEL .AND. Y_O2_F>TWO_EPSILON_EB .AND. CHAR_INDEX>0) THEN - ! The new char model starts the exp profile of O2 at the CHAR_FRONT + ! The new char model starts the exp profile of O2 at the ASH_DEPTH ! Find first cell where char has not been consumed - CHAR_FRONT_POINTS_LOOP: DO I=B2%I_CHAR_FRONT,NWP + ASH_DEPTH_POINTS_LOOP: DO I=B2%I_ASH_DEPTH,NWP IF (ONE_D%MATL_COMP(CHAR_INDEX)%RHO(I)>CHAR_DENSITY_THRESHOLD) THEN - CHAR_FRONT=ONE_D%X(I-1) - B2%I_CHAR_FRONT=I ! store last position to save time on next time step - EXIT CHAR_FRONT_POINTS_LOOP + ASH_DEPTH=ONE_D%X(I-1) + B2%I_ASH_DEPTH=I ! store last position to save time on next time step + EXIT ASH_DEPTH_POINTS_LOOP ENDIF - ENDDO CHAR_FRONT_POINTS_LOOP + ENDDO ASH_DEPTH_POINTS_LOOP ENDIF ! Initialize iterations for OXPYRO_MODEL @@ -2964,13 +2964,13 @@ SUBROUTINE PERFORM_PYROLYSIS IF (PRESENT(PARTICLE_INDEX)) THEN CALL PYROLYSIS(ONE_D%N_MATL,ONE_D%MATL_INDEX,SURF_INDEX,BC%IIG,BC%JJG,BC%KKG,ONE_D%TMP(I),B1%TMP_F,Y_O2_F,BC%IOR,& - RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),ONE_D%X(I-1),CHAR_FRONT,DX_S,DT_BC_SUB,& + RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),ONE_D%X(I-1),ASH_DEPTH,DX_S,DT_BC_SUB,& M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_S(I),Q_DOT_G_PPP,Q_DOT_O2_PPP,& Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B1%B_NUMBER,LAYER_INDEX(I),REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX=I,& R_DROP=R_SURF,LPU=U_SURF,LPV=V_SURF,LPW=W_SURF) ELSE CALL PYROLYSIS(ONE_D%N_MATL,ONE_D%MATL_INDEX,SURF_INDEX,BC%IIG,BC%JJG,BC%KKG,ONE_D%TMP(I),B1%TMP_F,Y_O2_F,BC%IOR,& - RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),0.5*(ONE_D%X(I-1)+ONE_D%X(I)),CHAR_FRONT,DX_S,DT_BC_SUB,& + RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),0.5*(ONE_D%X(I-1)+ONE_D%X(I)),ASH_DEPTH,DX_S,DT_BC_SUB,& M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_S(I),Q_DOT_G_PPP,Q_DOT_O2_PPP,& Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B1%B_NUMBER,LAYER_INDEX(I),REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX=I) ENDIF @@ -3057,7 +3057,7 @@ END SUBROUTINE SOLID_HEAT_TRANSFER !> \param RHO_DOT_OUT (1:N_MATS) Array of component reaction rates (kg/m3/s) !> \param RHO_S (1:N_MATS) Array of component densities (kg/m3) !> \param DEPTH Distance from surface (m) -!> \param CHAR_FRONT Distance from surface to start of char front (m) +!> \param ASH_DEPTH Distance from surface to start of char front (m) !> \param DX_S Array of node sizes (m) !> \param DT_BC Time step used by the solid phase solver (s) !> \param M_DOT_G_PPP_ADJUST (1:N_TRACKED_SPECIES) Adjusted mass generation rate per unit volume of the gas species @@ -3081,7 +3081,7 @@ END SUBROUTINE SOLID_HEAT_TRANSFER !> \param LPW (OPTIONAL) z component of droplet velocity (m/s) SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F,IOR,& - RHO_DOT_OUT,RHO_S,DEPTH,CHAR_FRONT,DX_S,DT_BC,& + RHO_DOT_OUT,RHO_S,DEPTH,ASH_DEPTH,DX_S,DT_BC,& M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_DOT_S_PPP,Q_DOT_G_PPP,Q_DOT_O2_PPP,& Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B_NUMBER,LAYER_INDEX,REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX,& R_DROP,LPU,LPV,LPW) @@ -3094,7 +3094,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F INTEGER, INTENT(IN), OPTIONAL :: SOLID_CELL_INDEX LOGICAL, INTENT(IN) :: REMOVE_LAYER REAL(EB), INTENT(OUT), DIMENSION(:,:) :: RHO_DOT_OUT(N_MATS) -REAL(EB), INTENT(IN) :: TMP_S,TMP_F,DT_BC,DEPTH,RHO_S(N_MATS),Y_O2_F,CHAR_FRONT +REAL(EB), INTENT(IN) :: TMP_S,TMP_F,DT_BC,DEPTH,RHO_S(N_MATS),Y_O2_F,ASH_DEPTH REAL(EB), INTENT(IN), OPTIONAL :: R_DROP,LPU,LPV,LPW REAL(EB), INTENT(IN), DIMENSION(NWP_MAX) :: DX_S REAL(EB), DIMENSION(:) :: ZZ_GET(1:N_TRACKED_SPECIES),Y_ALL(1:N_SPECIES) @@ -3377,7 +3377,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F ! Calculate oxygen volume fraction at the surface X_O2 = SPECIES(O2_INDEX)%RCON*Y_O2_F/RSUM(IIG,JJG,KKG) ! Calculate oxygen concentration inside the material, assuming decay function - X_O2 = X_O2 * EXP(-MAX(0._EB,DEPTH-CHAR_FRONT)/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(J))) + X_O2 = X_O2 * EXP(-MAX(0._EB,DEPTH-ASH_DEPTH)/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(J))) REACTION_RATE = REACTION_RATE * X_O2**ML%N_O2(J) ENDIF REACTION_RATE = MIN(REACTION_RATE,ML%MAX_REACTION_RATE(J)) ! User-specified limit