From ce9ef55f1790791bb5688b2a0f80333a3a8fe8b3 Mon Sep 17 00:00:00 2001 From: rmcdermo Date: Fri, 18 May 2018 17:22:03 -0400 Subject: [PATCH 1/2] FDS Source: adjust mass transfer formulation PYRO3D --- Source/dump.f90 | 1 - Source/init.f90 | 24 +++++- Source/wall.f90 | 200 +++++++++++++++++++++++++----------------------- 3 files changed, 124 insertions(+), 101 deletions(-) diff --git a/Source/dump.f90 b/Source/dump.f90 index f28f4b7ec90..a9aca4f734f 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -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 diff --git a/Source/init.f90 b/Source/init.f90 index eedc16ca806..eba65f54a2a 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -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 @@ -942,11 +942,11 @@ 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 @@ -954,6 +954,24 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) 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 diff --git a/Source/wall.f90 b/Source/wall.f90 index 19c02b5ba16..0b2dae23d93 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -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 @@ -1364,49 +1365,29 @@ 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 @@ -1414,7 +1395,8 @@ SUBROUTINE SOLID_MASS_TRANSFER_3D(DT_SUB) 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 @@ -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 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 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 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 @@ -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 @@ -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 @@ -2069,7 +2075,7 @@ SUBROUTINE CALCULATE_ZZ_F(WALL_INDEX) CONSUME_MASS: IF (CORRECTOR .AND. SF%THERMALLY_THICK .AND. .NOT.SF%THERMALLY_THICK_HT3D) THEN OB => OBSTRUCTION(WC%OBST_INDEX) DO N=1,N_TRACKED_SPECIES - OB%MASS = OB%MASS - ONE_D%MASSFLUX_SPEC(N)*DT*ONE_D%AREA + OB%MASS = OB%MASS - ONE_D%MASSFLUX(N)*DT*ONE_D%AREA ENDDO ENDIF CONSUME_MASS From 4745b34c42332713908cda4c18a0315b1ac00a39 Mon Sep 17 00:00:00 2001 From: rmcdermo Date: Fri, 18 May 2018 17:25:03 -0400 Subject: [PATCH 2/2] FDS Source: revert accidental change --- Source/wall.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/wall.f90 b/Source/wall.f90 index 0b2dae23d93..104c7b3a14b 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -2075,7 +2075,7 @@ SUBROUTINE CALCULATE_ZZ_F(WALL_INDEX) CONSUME_MASS: IF (CORRECTOR .AND. SF%THERMALLY_THICK .AND. .NOT.SF%THERMALLY_THICK_HT3D) THEN OB => OBSTRUCTION(WC%OBST_INDEX) DO N=1,N_TRACKED_SPECIES - OB%MASS = OB%MASS - ONE_D%MASSFLUX(N)*DT*ONE_D%AREA + OB%MASS = OB%MASS - ONE_D%MASSFLUX_SPEC(N)*DT*ONE_D%AREA ENDDO ENDIF CONSUME_MASS