Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

FDS Source: adjust mass transfer formulation for PYRO3D #6431

Merged
merged 2 commits into from May 18, 2018
Merged

FDS Source: adjust mass transfer formulation for PYRO3D #6431

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ce9ef55
FDS Source: adjust mass transfer formulation PYRO3D
rmcdermo May 18, 2018
4745b34
FDS Source: revert accidental change
rmcdermo May 18, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 0 additions & 1 deletion Source/dump.f90
View file
Open in desktop
Expand Up @@ -2573,7 +2573,6 @@ SUBROUTINE INITIALIZE_DIAGNOSTIC_FILE(DT)
DO N=1,N_TRACKED_SPECIES
SM=>SPECIES_MIXTURE(N)
WRITE(LU_OUTPUT,'(/3X,A)') TRIM(SM%ID)
IF (N==0) WRITE(LU_OUTPUT,'( 3X,A)') 'Background Species'
WRITE(LU_OUTPUT,'(A,F11.5)') ' Molecular Weight (g/mol) ',SM%MW
WRITE(LU_OUTPUT,'(A,F8.3)') ' Ambient Density (kg/m^3) ',SM%MW*P_INF/(TMPA*R0)
WRITE(LU_OUTPUT,'(A,F8.3)') ' Initial Mass Fraction ',SM%ZZ0
Expand Down
24 changes: 21 additions & 3 deletions Source/init.f90
View file
Open in desktop
Expand Up @@ -30,7 +30,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
USE RADCONS, ONLY: UIIDIM
USE CONTROL_VARIABLES
USE MATH_FUNCTIONS, ONLY:EVALUATE_RAMP
INTEGER :: N,NNN,I,J,K,IW,IC,SURF_INDEX,IOR,IERR,IIG,JJG,KKG,N_OVERLAP,IZERO
INTEGER :: N,NN,I,J,K,IW,IC,SURF_INDEX,IOR,IERR,IIG,JJG,KKG,N_OVERLAP,IZERO
REAL(EB), INTENT(IN) :: DT
INTEGER, INTENT(IN) :: NM
REAL(EB) :: MU_N,ZZ_GET(1:N_TRACKED_SPECIES),CS,DELTA,DUMMY,INTEGRAL,TEMP
Expand Down Expand Up @@ -942,18 +942,36 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
OB=>M%OBSTRUCTION(N)
IF (OB%MATL_SURF_INDEX<=0) CYCLE OBST_LOOP_4
SF=>SURFACE(OB%MATL_SURF_INDEX)
MATL_LOOP: DO NNN = 1,SF%N_MATL
MATL_LOOP: DO NN = 1,SF%N_MATL
DO K=OB%K1+1,OB%K2
DO J=OB%J1+1,OB%J2
DO I=OB%I1+1,OB%I2
OB%RHO(I,J,K,NNN) = SF%RHO_0(1,NNN)
OB%RHO(I,J,K,NN) = SF%RHO_0(1,NN)
ENDDO
ENDDO
ENDDO
ENDDO MATL_LOOP
ENDDO OBST_LOOP_4
ENDIF

! Initialize gas densities for 3D pyrolysis

IF (SOLID_MT3D) THEN
OBST_LOOP_5: DO N=1,M%N_OBST
OB=>M%OBSTRUCTION(N)
IF (.NOT.OB%MT3D) CYCLE OBST_LOOP_5
DO K=OB%K1+1,OB%K2
DO J=OB%J1+1,OB%J2
DO I=OB%I1+1,OB%I2
M%RHO(I,J,K) = 0._EB
DO NN=1,N_TRACKED_SPECIES
M%ZZ(I,J,K,NN) = 0._EB
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO OBST_LOOP_5
ENDIF

! Allocate local auto-ignition temperature

Expand Down
198 changes: 102 additions & 96 deletions Source/wall.f90
View file
Open in desktop
Expand Up @@ -1239,16 +1239,17 @@ SUBROUTINE SOLID_PYROLYSIS_3D(DT_SUB,T_LOC)
ZZ(I,J,K,NS) = MAX( 0._EB, RHO(I,J,K)*ZZ(I,J,K,NS) + DT_SUB*M_DOT_G_PPP_ADJUST(NS) )
ENDDO
RHO(I,J,K) = SUM(ZZ(I,J,K,1:N_TRACKED_SPECIES))
ZZ(I,J,K,1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)/RHO(I,J,K)
IF (RHO(I,J,K)>TWO_EPSILON_EB) ZZ(I,J,K,1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)/RHO(I,J,K)
ELSE
! simple model (no transport): pyrolyzed mass is ejected via wall cell index WALL_INDEX_HT3D(IC,OB%PYRO3D_IOR)
DO NS=1,N_TRACKED_SPECIES
WC%ONE_D%MASSFLUX(NS) = WC%ONE_D%MASSFLUX(NS) + M_DOT_G_PPP_ADJUST(NS)*GEOM_FACTOR*TIME_FACTOR
WC%ONE_D%MASSFLUX_SPEC(NS) = WC%ONE_D%MASSFLUX_SPEC(NS) + M_DOT_G_PPP_ACTUAL(NS)*GEOM_FACTOR*TIME_FACTOR
ENDDO
DO NN=1,SF%N_MATL
WC%ONE_D%MASSFLUX_MATL(NN) = WC%ONE_D%MASSFLUX_MATL(NN) + M_DOT_S_PPP(NN)*GEOM_FACTOR*TIME_FACTOR
ENDDO
!! MASSFLUX_MATL should not be needed in 3D
! DO NN=1,SF%N_MATL
! WC%ONE_D%MASSFLUX_MATL(NN) = WC%ONE_D%MASSFLUX_MATL(NN) + M_DOT_S_PPP(NN)*GEOM_FACTOR*TIME_FACTOR
! ENDDO
! If the fuel or water massflux is non-zero, set the ignition time
IF (WC%ONE_D%T_IGN > T) THEN
IF (SUM(WC%ONE_D%MASSFLUX(1:N_TRACKED_SPECIES)) > 0._EB) WC%ONE_D%T_IGN = T
Expand Down Expand Up @@ -1364,57 +1365,38 @@ END SUBROUTINE SOLID_PYROLYSIS_3D

SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)

! Based on C. Lautenberger, C. Fernandez-Pello / Fire Safety Journal 44 (2009) 819-839.

USE MATH_FUNCTIONS, ONLY: INTERPOLATE1D_UNIFORM
REAL(EB), INTENT(IN) :: DT_SUB
INTEGER :: I,J,K,N,IC,NN,II,JJ,KK,IOR,ICM,ICP
REAL(EB) :: D_Z_TEMP,D_Z_N(0:5000),VOLSUM,PSISUM,VOLRAT,RHO_D_F,R_RHO_D,VN_MT3D,ZZ_I,RHO_D_MAX,RHO_D_M,RHO_D_P,RHO_D_BAR
REAL(EB), POINTER, DIMENSION(:,:,:) :: PSIBAR=>NULL(),RHO_D_DZDX=>NULL(),RHO_D_DZDY=>NULL(),RHO_D_DZDZ=>NULL(),RHO_D=>NULL()
INTEGER :: I,J,K,N,IC,II,JJ,KK,IOR,ICM,ICP,IIG,JJG,KKG
REAL(EB) :: D_Z_TEMP,D_Z_N(0:5000),D_F,R_D,VN_MT3D,RHO_ZZ_I,D_MAX,D_M,D_P,D_BAR,RDN
REAL(EB), POINTER, DIMENSION(:,:,:) :: D_DRHOZDX=>NULL(),D_DRHOZDY=>NULL(),D_DRHOZDZ=>NULL(),D_Z_P=>NULL(),RHO_ZZ_P=>NULL()
LOGICAL :: CONT_MATL_PROP
TYPE(OBSTRUCTION_TYPE), POINTER :: OB=>NULL(),OBM=>NULL(),OBP=>NULL()
TYPE(SURFACE_TYPE), POINTER :: MS=>NULL()
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()

! compute cell porosity, LFP Eq. (5)

PSIBAR=>WORK1; PSIBAR=0._EB
DO K=0,KBP1
DO J=0,JBP1
DO I=0,IBP1
IC = CELL_INDEX(I,J,K); IF (.NOT.SOLID(IC)) CYCLE
OB => OBSTRUCTION(OBST_INDEX_C(IC)); IF (.NOT.OB%MT3D) CYCLE
MS => SURFACE(OB%MATL_SURF_INDEX)
VOLSUM = 0._EB
DO NN=1,MS%N_MATL
ML => MATERIAL(MS%MATL_INDEX(NN))
VOLRAT = OB%RHO(I,J,K,NN)/ML%RHO_S
VOLSUM = VOLSUM + VOLRAT
PSISUM = PSISUM + VOLRAT*(1._EB - VOLRAT)
ENDDO
IF (VOLSUM>TWO_EPSILON_EB) PSIBAR(I,J,K) = PSISUM/VOLSUM
ENDDO
ENDDO
ENDDO

SPECIES_LOOP: DO N=1,N_TRACKED_SPECIES
! initialize work arrays

D_Z_P =>WORK1; D_Z_P=0._EB
D_DRHOZDX=>WORK2; D_DRHOZDX=0._EB
D_DRHOZDY=>WORK3; D_DRHOZDY=0._EB
D_DRHOZDZ=>WORK4; D_DRHOZDZ=0._EB
RHO_ZZ_P =>WORK5; RHO_ZZ_P=0._EB

RHO_D_DZDX=>WORK2
RHO_D_DZDY=>WORK3
RHO_D_DZDZ=>WORK4
! loop over all tracked gas species

! get gas phase diffusivity
SPECIES_LOOP: DO N=1,N_TRACKED_SPECIES

! get gas phase diffusivity and density

RHO_D => WORK5
RHO_D = 0._EB
D_Z_N = D_Z(:,N)
DO K=0,KBP1
DO J=0,JBP1
DO I=0,IBP1
IC = CELL_INDEX(I,J,K); IF (.NOT.SOLID(IC)) CYCLE
OB => OBSTRUCTION(OBST_INDEX_C(IC)); IF (.NOT.OB%MT3D) CYCLE
CALL INTERPOLATE1D_UNIFORM(LBOUND(D_Z_N,1),D_Z_N,TMP(I,J,K),D_Z_TEMP)
RHO_D(I,J,K) = PSIBAR(I,J,K)*RHO(I,J,K)*D_Z_TEMP
D_Z_P(I,J,K) = D_Z_TEMP
RHO_ZZ_P(I,J,K) = RHO(I,J,K)*ZZ(I,J,K,N)
ENDDO
ENDDO
ENDDO
Expand All @@ -1436,11 +1418,11 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)
! 2. we assume that if either OBM%MT3D .OR. OBP%MT3D we should process the boundary
IF (.NOT.(OBM%MT3D.OR.OBP%MT3D)) CYCLE

RHO_D_M = RHO_D(I,J,K)
RHO_D_P = RHO_D(I+1,J,K)
D_M = D_Z_P(I,J,K)
D_P = D_Z_P(I+1,J,K)

IF (RHO_D_M<TWO_EPSILON_EB .OR. RHO_D_P<TWO_EPSILON_EB) THEN
RHO_D_DZDX(I,J,K) = 0._EB
IF (D_M<TWO_EPSILON_EB .OR. D_P<TWO_EPSILON_EB) THEN
D_DRHOZDX(I,J,K) = 0._EB
CYCLE
ENDIF

Expand All @@ -1458,17 +1440,17 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)

IF (CONT_MATL_PROP) THEN
! use linear average from inverse lever rule
RHO_D_BAR = ( RHO_D_M*DX(I+1) + RHO_D_P*DX(I) )/( DX(I) + DX(I+1) )
RHO_D_MAX = MAX(RHO_D_MAX,RHO_D_BAR)
RHO_D_DZDX(I,J,K) = RHO_D_BAR*(ZZ(I+1,J,K,N)-ZZ(I,J,K,N))*2._EB/(DX(I+1)+DX(I))
D_BAR = ( D_M*DX(I+1) + D_P*DX(I) )/( DX(I) + DX(I+1) )
D_MAX = MAX(D_MAX,D_BAR)
D_DRHOZDX(I,J,K) = D_BAR*(RHO_ZZ_P(I+1,J,K)-RHO_ZZ_P(I,J,K))*2._EB/(DX(I+1)+DX(I))
ELSE
! for discontinuous material properties maintain continuity of flux, C0 continuity of composition
! (allow C1 discontinuity of composition due to jump in properties across interface)
R_RHO_D = RHO_D_P/RHO_D_M * DX(I)/DX(I+1)
ZZ_I = (ZZ(I,J,K,N) + R_RHO_D*ZZ(I+1,J,K,N))/(1._EB + R_RHO_D) ! interface concentration
!! RHO_D_DZDX(I,J,K) = RHO_D_P * (ZZ(I+1,J,K,N)-ZZ_I) * 2._EB/DX(I+1) !! should be identical
RHO_D_DZDX(I,J,K) = RHO_D_M * (ZZ_I-ZZ(I,J,K,N)) * 2._EB/DX(I)
RHO_D_MAX = MAX(RHO_D_MAX,MAX(RHO_D_M,RHO_D_P))
R_D = D_P/D_M * DX(I)/DX(I+1)
RHO_ZZ_I = (RHO_ZZ_P(I,J,K) + R_D*RHO_ZZ_P(I+1,J,K))/(1._EB + R_D) ! interface concentration
!! D_DRHOZDX(I,J,K) = D_P * (RHO_ZZ_P(I+1,J,K)-RHO_ZZ_I) * 2._EB/DX(I+1) !! should be identical
D_DRHOZDX(I,J,K) = D_M * (RHO_ZZ_I-RHO_ZZ_P(I,J,K)) * 2._EB/DX(I)
D_MAX = MAX(D_MAX,MAX(D_M,D_P))
ENDIF
ENDDO
ENDDO
Expand All @@ -1484,11 +1466,11 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)
OBP => OBSTRUCTION(OBST_INDEX_C(ICP))
IF (.NOT.(OBM%MT3D.OR.OBP%MT3D)) CYCLE

RHO_D_M = RHO_D(I,J,K)
RHO_D_P = RHO_D(I,J+1,K)
D_M = D_Z_P(I,J,K)
D_P = D_Z_P(I,J+1,K)

IF (RHO_D_M<TWO_EPSILON_EB .OR. RHO_D_P<TWO_EPSILON_EB) THEN
RHO_D_DZDY(I,J,K) = 0._EB
IF (D_M<TWO_EPSILON_EB .OR. D_P<TWO_EPSILON_EB) THEN
D_DRHOZDY(I,J,K) = 0._EB
CYCLE
ENDIF

Expand All @@ -1504,20 +1486,20 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)
ENDIF

IF (CONT_MATL_PROP) THEN
RHO_D_BAR = ( RHO_D_M*DY(J+1) + RHO_D_P*DY(J) )/( DY(J) + DY(J+1) )
RHO_D_MAX = MAX(RHO_D_MAX,RHO_D_BAR)
RHO_D_DZDY(I,J,K) = RHO_D_BAR*(ZZ(I,J+1,K,N)-ZZ(I,J,K,N))*2._EB/(DY(J+1)+DY(J))
D_BAR = ( D_M*DY(J+1) + D_P*DY(J) )/( DY(J) + DY(J+1) )
D_MAX = MAX(D_MAX,D_BAR)
D_DRHOZDY(I,J,K) = D_BAR*(RHO_ZZ_P(I,J+1,K)-RHO_ZZ_P(I,J,K))*2._EB/(DY(J+1)+DY(J))
ELSE
R_RHO_D = RHO_D_P/RHO_D_M * DY(J)/DY(J+1)
ZZ_I = (ZZ(I,J,K,N) + R_RHO_D*ZZ(I,J+1,K,N))/(1._EB + R_RHO_D)
RHO_D_DZDY(I,J,K) = RHO_D_M * (ZZ_I-ZZ(I,J,K,N)) * 2._EB/DY(J)
RHO_D_MAX = MAX(RHO_D_MAX,MAX(RHO_D_M,RHO_D_P))
R_D = D_P/D_M * DY(J)/DY(J+1)
RHO_ZZ_I = (RHO_ZZ_P(I,J,K) + R_D*RHO_ZZ_P(I,J+1,K))/(1._EB + R_D)
D_DRHOZDY(I,J,K) = D_M * (RHO_ZZ_I-RHO_ZZ_P(I,J,K)) * 2._EB/DY(J)
D_MAX = MAX(D_MAX,MAX(D_M,D_P))
ENDIF
ENDDO
ENDDO
ENDDO
ELSE TWO_D_IF
RHO_D_DZDY(I,J,K) = 0._EB
D_DRHOZDY(I,J,K) = 0._EB
ENDIF TWO_D_IF
DO K=0,KBAR
DO J=1,JBAR
Expand All @@ -1529,11 +1511,11 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)
OBP => OBSTRUCTION(OBST_INDEX_C(ICP))
IF (.NOT.(OBM%MT3D.OR.OBP%MT3D)) CYCLE

RHO_D_M = RHO_D(I,J,K)
RHO_D_P = RHO_D(I,J,K+1)
D_M = D_Z_P(I,J,K)
D_P = D_Z_P(I,J,K+1)

IF (RHO_D_M<TWO_EPSILON_EB .OR. RHO_D_P<TWO_EPSILON_EB) THEN
RHO_D_DZDZ(I,J,K) = 0._EB
IF (D_M<TWO_EPSILON_EB .OR. D_P<TWO_EPSILON_EB) THEN
D_DRHOZDZ(I,J,K) = 0._EB
CYCLE
ENDIF

Expand All @@ -1549,14 +1531,14 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)
ENDIF

IF (CONT_MATL_PROP) THEN
RHO_D_BAR = ( RHO_D_M*DZ(K+1) + RHO_D_P*DZ(K) )/( DZ(K) + DZ(K+1) )
RHO_D_MAX = MAX(RHO_D_MAX,RHO_D_BAR)
RHO_D_DZDZ(I,J,K) = RHO_D_BAR*(ZZ(I,J,K+1,N)-ZZ(I,J,K,N))*2._EB/(DZ(K+1)+DZ(K))
D_BAR = ( D_M*DZ(K+1) + D_P*DZ(K) )/( DZ(K) + DZ(K+1) )
D_MAX = MAX(D_MAX,D_BAR)
D_DRHOZDZ(I,J,K) = D_BAR*(RHO_ZZ_P(I,J,K+1)-RHO_ZZ_P(I,J,K))*2._EB/(DZ(K+1)+DZ(K))
ELSE
R_RHO_D = RHO_D_P/RHO_D_M * DZ(K)/DZ(K+1)
ZZ_I = (ZZ(I,J,K,N) + R_RHO_D*ZZ(I,J,K+1,N))/(1._EB + R_RHO_D)
RHO_D_DZDZ(I,J,K) = RHO_D_M * (ZZ_I-ZZ(I,J,K,N)) * 2._EB/DZ(K)
RHO_D_MAX = MAX(RHO_D_MAX,MAX(RHO_D_M,RHO_D_P))
R_D = D_P/D_M * DZ(K)/DZ(K+1)
RHO_ZZ_I = (RHO_ZZ_P(I,J,K) + R_D*RHO_ZZ_P(I,J,K+1))/(1._EB + R_D)
D_DRHOZDZ(I,J,K) = D_M * (RHO_ZZ_I-RHO_ZZ_P(I,J,K)) * 2._EB/DZ(K)
D_MAX = MAX(D_MAX,MAX(D_M,D_P))
ENDIF
ENDDO
ENDDO
Expand All @@ -1573,30 +1555,55 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)
II = WC%ONE_D%II
JJ = WC%ONE_D%JJ
KK = WC%ONE_D%KK
IIG = WC%ONE_D%IIG
JJG = WC%ONE_D%JJG
KKG = WC%ONE_D%KKG
IOR = WC%ONE_D%IOR

IC = CELL_INDEX(II,JJ,KK); IF (.NOT.SOLID(IC)) CYCLE MT3D_WALL_LOOP
OB => OBSTRUCTION(OBST_INDEX_C(IC)); IF (.NOT.OB%MT3D ) CYCLE MT3D_WALL_LOOP
MS => SURFACE(OB%MATL_SURF_INDEX)

CALL INTERPOLATE1D_UNIFORM(LBOUND(D_Z_N,1),D_Z_N,WC%ONE_D%TMP_F,D_Z_TEMP)
RHO_D_F = PSIBAR(II,JJ,KK)*RHO(II,JJ,KK)*D_Z_TEMP
RHO_D_MAX = MAX(RHO_D_MAX,RHO_D_F)
D_F = D_Z_TEMP
D_MAX = MAX(D_MAX,D_F)

METHOD_OF_MASS_TRANSFER: SELECT CASE(SF%THERMAL_BC_INDEX)
SELECT CASE(IOR)
CASE( 1); RDN=RDX(II)
CASE(-1); RDN=RDX(II)
CASE( 2); RDN=RDY(JJ)
CASE(-2); RDN=RDY(JJ)
CASE( 3); RDN=RDZ(KK)
CASE(-3); RDN=RDZ(KK)
END SELECT

METHOD_OF_MASS_TRANSFER: SELECT CASE(SF%SPECIES_BC_INDEX)

CASE DEFAULT METHOD_OF_MASS_TRANSFER

! for simplicity, assume surface mass fraction is equivalent to gas phase adjacent to wall cell

WC%ONE_D%ZZ_F(N) = ZZ(IIG,JJG,KKG,N)

! compute mass flux at the surface

WC%ONE_D%MASSFLUX_SPEC(N) = - D_F * 2._EB*(WC%ONE_D%RHO_F*WC%ONE_D%ZZ_F(N)-RHO_ZZ_P(II,JJ,KK))*RDN
WC%ONE_D%MASSFLUX(N) = WC%ONE_D%MASSFLUX_SPEC(N) * WC%ONE_D%AREA_ADJUST

SELECT CASE(IOR)
CASE( 1); RHO_D_DZDX(II,JJ,KK) = RHO_D_F * 2._EB*(WC%ONE_D%ZZ_F(N)-ZZ(II,JJ,KK,N))*RDX(II)
CASE(-1); RHO_D_DZDX(II-1,JJ,KK) = RHO_D_F * 2._EB*(ZZ(II,JJ,KK,N)-WC%ONE_D%ZZ_F(N))*RDX(II)
CASE( 2); RHO_D_DZDY(II,JJ,KK) = RHO_D_F * 2._EB*(WC%ONE_D%ZZ_F(N)-ZZ(II,JJ,KK,N))*RDY(JJ)
CASE(-2); RHO_D_DZDY(II,JJ-1,KK) = RHO_D_F * 2._EB*(ZZ(II,JJ,KK,N)-WC%ONE_D%ZZ_F(N))*RDY(JJ)
CASE( 3); RHO_D_DZDZ(II,JJ,KK) = RHO_D_F * 2._EB*(WC%ONE_D%ZZ_F(N)-ZZ(II,JJ,KK,N))*RDZ(KK)
CASE(-3); RHO_D_DZDZ(II,JJ,KK-1) = RHO_D_F * 2._EB*(ZZ(II,JJ,KK,N)-WC%ONE_D%ZZ_F(N))*RDZ(KK)
CASE( 1); D_DRHOZDX(II,JJ,KK) = WC%ONE_D%MASSFLUX(N)
CASE(-1); D_DRHOZDX(II-1,JJ,KK) = -WC%ONE_D%MASSFLUX(N)
CASE( 2); D_DRHOZDY(II,JJ,KK) = WC%ONE_D%MASSFLUX(N)
CASE(-2); D_DRHOZDY(II,JJ-1,KK) = -WC%ONE_D%MASSFLUX(N)
CASE( 3); D_DRHOZDZ(II,JJ,KK) = WC%ONE_D%MASSFLUX(N)
CASE(-3); D_DRHOZDZ(II,JJ,KK-1) = -WC%ONE_D%MASSFLUX(N)
END SELECT

! NEED TO POPULATE MASSFLUX
! If the fuel or water massflux is non-zero, set the ignition time

IF (WC%ONE_D%T_IGN > T) THEN
IF (ABS(WC%ONE_D%MASSFLUX(N)) > 0._EB) WC%ONE_D%T_IGN = T
ENDIF

END SELECT METHOD_OF_MASS_TRANSFER

Expand All @@ -1612,16 +1619,17 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)
OB => OBSTRUCTION(OBST_INDEX_C(IC)); IF (.NOT.OB%MT3D) CYCLE

IF (TWO_D) THEN
VN_MT3D = MAX( VN_MT3D, 2._EB*RHO_D_MAX/RHO(I,J,K)*( RDX(I)**2 + RDZ(K)**2 ) )
VN_MT3D = MAX( VN_MT3D, 2._EB*D_MAX*( RDX(I)**2 + RDZ(K)**2 ) )
ELSE
VN_MT3D = MAX( VN_MT3D, 2._EB*RHO_D_MAX/RHO(I,J,K)*( RDX(I)**2 + RDY(J)**2 + RDZ(K)**2 ) )
VN_MT3D = MAX( VN_MT3D, 2._EB*D_MAX*( RDX(I)**2 + RDY(J)**2 + RDZ(K)**2 ) )
ENDIF

ZZ(I,J,K,N) = RHO(I,J,K)*ZZ(I,J,K,N) + DT_SUB * ( (RHO_D_DZDX(I,J,K)*R(I)-RHO_D_DZDX(I-1,J,K)*R(I-1))*RDX(I)*RRN(I) + &
(RHO_D_DZDY(I,J,K) -RHO_D_DZDY(I,J-1,K) )*RDY(J) + &
(RHO_D_DZDZ(I,J,K) -RHO_D_DZDZ(I,J,K-1) )*RDZ(K) )
ZZ(I,J,K,N) = RHO_ZZ_P(I,J,K) + DT_SUB * ( (D_DRHOZDX(I,J,K)*R(I)-D_DRHOZDX(I-1,J,K)*R(I-1))*RDX(I)*RRN(I) + &
(D_DRHOZDY(I,J,K) -D_DRHOZDY(I,J-1,K) )*RDY(J) + &
(D_DRHOZDZ(I,J,K) -D_DRHOZDZ(I,J,K-1) )*RDZ(K) )

ZZ(I,J,K,N) = MAX(0._EB,ZZ(I,J,K,N)) ! guarantee boundedness

ENDDO
ENDDO
ENDDO
Expand All @@ -1638,18 +1646,16 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB)
OB => OBSTRUCTION(OBST_INDEX_C(IC)); IF (.NOT.OB%MT3D) CYCLE

RHO(I,J,K) = SUM(ZZ(I,J,K,1:N_TRACKED_SPECIES))
ZZ(I,J,K,1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)/RHO(I,J,K)
N=MAXLOC(ZZ(I,J,K,1:N_TRACKED_SPECIES),1) ! absorb error in most abundant local species
ZZ(I,J,K,N) = 1._EB - ( SUM(ZZ(I,J,K,1:N_TRACKED_SPECIES)) - ZZ(I,J,K,N) ) ! enforce sum(ZZ)=1
IF (RHO(I,J,K)>TWO_EPSILON_EB) ZZ(I,J,K,1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)/RHO(I,J,K)
ENDDO
ENDDO
ENDDO

! stability check
! ! stability check

IF (DT_SUB*VN_MT3D > VN_MAX) THEN
print *, 'VN:', DT_SUB*VN_MT3D
ENDIF
! IF (DT_SUB*VN_MT3D > VN_MAX) THEN
! print *, 'VN:', DT_SUB*VN_MT3D
! ENDIF

END SUBROUTINE SOLID_MASS_TRANSFER_3D

Expand Down