FDS Source: Fix Kn number calculation and associated verification Issue #13926

Manuals/FDS_Technical_Reference_Guide/Aerosol_Chapter.tex

@@ -111,11 +111,11 @@ \section{Aerosol Agglomeration}
 \be
 \frac{\d N_i}{\d t} = \sum_{k=1}^{M} \sum_{j=k}^{M} \left( 1- \frac{1}{2} \delta_{j,k} \right) \eta \Phi(j,k) N_j N_k - N_i \sum_{k=1}^M \Phi(j,k) N_k - R_i + S_i
 \ee
-where $\eta$, shown below, is a function for apportioning mass between adjacent particle bins.
+where $\eta$, shown below, is a function for apportioning the new particle mass, $m_{j+k}=x_j+x_k$, between adjacent particle bins.
 \be
 \eta = \begin{cases}
-   \frac{x_{i+1}-m_i}{x_{i+1}-x_i} & x_i< m_i \leq x_{i+1} \\
-   \frac{m_i-x_{i-1}}{x_i-x_{i-1}} & x_{i-1}< m_i \leq x_i
+   \frac{x_{i+1}-m_{j+k}}{x_{i+1}-x_i} & x_i< m_{j+k} \leq x_{i+1} \\
+   \frac{m_{j+k}-x_{i-1}}{x_i-x_{i-1}} & x_{i-1}< m_{j+k} \leq x_i
 \end{cases}
 \ee
 The agglomeration kernel is given by the sum of the Brownian (Br), gravitational (Gr), shear (Sh), and inertial (In) agglomeration  terms.

Source/func.f90

@@ -6172,4 +6172,4 @@ SUBROUTINE ACCUMULATE_STRING(STRING_SIZE,MYSTR,ACCSTR,ACCSTR_T_LEN,ACCSTR_USE_LE
 END SUBROUTINE ACCUMULATE_STRING
 
 END MODULE MISC_FUNCTIONS
-
+   

Source/soot.f90

@@ -67,9 +67,10 @@ SUBROUTINE SETTLING_VELOCITY(NM)
             PRES_G=0.5_EB*(PBAR(K,PRESSURE_ZONE(I,J,K))+PBAR(K,PRESSURE_ZONE(I+1,J,K)))
             DTDN = -(TMP(I+1,J,K)-TMP(I,J,K))*RDXN(I)
             CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_G)
-            KN_FAC = SQRT(0.5_EB*PI/(PRES_G*RHO_G))*MU_G
+            ! Kn=2 lambda/d, lambda has sqrt(1/2). kn_fac has 2*sqrt(1/2)=sqrt(4/2)=sqrt(2)
+            KN_FAC = SQRT(2._EB*PI/(PRES_G*RHO_G))*MU_G 
             IF (THERMOPHORETIC_SETTLING) THEN
-               KN = KN_FAC/SPECIES_MIXTURE(N)%THERMOPHORETIC_DIAMETER
+               KN = KN_FAC/SPECIES_MIXTURE(N)%THERMOPHORETIC_DIAMETER 
                CN = CUNNINGHAM(KN)
                CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP_G)
                ALPHA = K_G/SPECIES_MIXTURE(N)%CONDUCTIVITY_SOLID
@@ -89,7 +90,7 @@ SUBROUTINE SETTLING_VELOCITY(NM)
             PRES_G=0.5_EB*(PBAR(K,PRESSURE_ZONE(I,J,K))+PBAR(K,PRESSURE_ZONE(I,J+1,K)))
             DTDN = -(TMP(I,J+1,K)-TMP(I,J,K))*RDYN(J)
             CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_G)
-            KN_FAC = SQRT(0.5_EB*PI/(PRES_G*RHO_G))*MU_G
+            KN_FAC = SQRT(2._EB*PI/(PRES_G*RHO_G))*MU_G
             IF (THERMOPHORETIC_SETTLING) THEN
                KN = KN_FAC/SPECIES_MIXTURE(N)%THERMOPHORETIC_DIAMETER
                CN = CUNNINGHAM(KN)
@@ -111,7 +112,7 @@ SUBROUTINE SETTLING_VELOCITY(NM)
             PRES_G=0.5_EB*(PBAR(K,PRESSURE_ZONE(I,J,K))+PBAR(K,PRESSURE_ZONE(I,J,K+1)))
             DTDN = -(TMP(I,J,K+1)-TMP(I,J,K))*RDZN(K)
             CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_G)
-            KN_FAC = SQRT(0.5_EB*PI/(PRES_G*RHO_G))*MU_G
+            KN_FAC = SQRT(2._EB*PI/(PRES_G*RHO_G))*MU_G
             IF (THERMOPHORETIC_SETTLING) THEN
                KN = KN_FAC/SPECIES_MIXTURE(N)%THERMOPHORETIC_DIAMETER
                CN = CUNNINGHAM(KN)
@@ -177,9 +178,9 @@ END SUBROUTINE SETTLING_VELOCITY
 
 
 SUBROUTINE INITIALIZE_AGGLOMERATION
+USE GLOBAL_CONSTANTS, ONLY: GRAVITATIONAL_SETTLING
 INTEGER :: N,I,II,III
 REAL(EB) :: E_PK,MIN_PARTICLE_MASS,MAX_PARTICLE_MASS
-
 ALLOCATE(BIN_S(N_AGGLOMERATION_SPECIES))
 ALLOCATE(BIN_M(N_AGGLOMERATION_SPECIES,0:MAXVAL(N_PARTICLE_BINS)))
 ALLOCATE(BIN_X(N_AGGLOMERATION_SPECIES,1:MAXVAL(N_PARTICLE_BINS)))
@@ -211,8 +212,6 @@ SUBROUTINE INITIALIZE_AGGLOMERATION
       BIN_X(N,I) = 2._EB*BIN_M(N,I)/(1._EB+BIN_S(N))
    ENDDO
 
-
-
    DO I=1,N_PARTICLE_BINS(N)
       PARTICLE_RADIUS(N,I) = (BIN_X(N,I) / FOTHPI / SPECIES(AGGLOMERATION_SPEC_INDEX(N))%DENSITY_SOLID)**ONTH
       MOBILITY_FAC(N,I) = 1._EB/(6._EB*PI*PARTICLE_RADIUS(N,I))
@@ -222,7 +221,11 @@ SUBROUTINE INITIALIZE_AGGLOMERATION
       DO II=1,N_PARTICLE_BINS(N)
          PARTICLE_MASS(N,I,II) = BIN_X(N,I) + BIN_X(N,II)
          E_PK = MIN(PARTICLE_RADIUS(N,I),PARTICLE_RADIUS(N,II))**2/(2._EB*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))**2)
-         PHI_G_FAC(N,I,II) = E_PK*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))**2
+         IF (GRAVITATIONAL_SETTLING) THEN
+            PHI_G_FAC(N,I,II) = E_PK*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))**2
+         ELSE
+            PHI_G_FAC(N,I,II) = 0._EB
+         ENDIF
          PHI_B_FAC(N,I,II)= 4._EB*PI*K_BOLTZMANN*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))
          !Check what Re and r are for PHI_I and _S
          PHI_S_FAC(N,I,II) = E_PK*(PARTICLE_RADIUS(N,I)+PARTICLE_RADIUS(N,II))**3*SQRT(8._EB*PI/15._EB)
@@ -234,21 +237,21 @@ SUBROUTINE INITIALIZE_AGGLOMERATION
          BINDO:DO III=2,N_PARTICLE_BINS(N)
             IF (PARTICLE_MASS(N,I,II) > BIN_X(N,N_PARTICLE_BINS(N))) THEN
                BIN_ETA_INDEX(N,I,II,:) = N_PARTICLE_BINS(N)
-               BIN_ETA(N,I,II,:) = 0.5_EB*BIN_X(N,N_PARTICLE_BINS(N))/PARTICLE_MASS(N,I,II)
+               BIN_ETA(N,I,II,:) = 0.5_EB*PARTICLE_MASS(N,I,II)/BIN_X(N,N_PARTICLE_BINS(N))
                EXIT BINDO
             ELSE
                IF (PARTICLE_MASS(N,I,II) > BIN_X(N,III-1) .AND. PARTICLE_MASS(N,I,II) < BIN_X(N,III)) THEN
                   BIN_ETA_INDEX(N,I,II,1) = III-1
                   BIN_ETA(N,I,II,1) = (BIN_X(N,III)-PARTICLE_MASS(N,I,II))/(BIN_X(N,III)-BIN_X(N,III-1))
                   BIN_ETA_INDEX(N,I,II,2) = III
                   BIN_ETA(N,I,II,2) = (PARTICLE_MASS(N,I,II)-BIN_X(N,III-1))/(BIN_X(N,III)-BIN_X(N,III-1))
-                  IF (I==II) BIN_ETA(N,I,II,:) = BIN_ETA(N,I,II,:) *0.5_EB
                   EXIT BINDO
                ENDIF
             ENDIF
          ENDDO BINDO
       ENDDO
    ENDDO
+
    BIN_ETA(N,N_PARTICLE_BINS(N),N_PARTICLE_BINS(N),1) = 1._EB
    BIN_ETA_INDEX(N,N_PARTICLE_BINS(N),N_PARTICLE_BINS(N),1) = N_PARTICLE_BINS(N)
    BIN_ETA(N,N_PARTICLE_BINS(N),N_PARTICLE_BINS(N),2) = 0._EB
@@ -265,8 +268,8 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
 INTEGER, INTENT(IN) :: NM
 REAL(EB), INTENT(IN) :: DT
 REAL(EB) :: DUDX,DVDY,DWDZ,ONTHDIV,S11,S22,S33,DUDY,DUDZ,DVDX,DVDZ,DWDX,DWDY,S12,S23,S13,STRAIN_RATE,DISSIPATION_RATE
-REAL(EB) :: KN,MFP,N_I(MAXVAL(N_PARTICLE_BINS)),N0(MAXVAL(N_PARTICLE_BINS)),N1(MAXVAL(N_PARTICLE_BINS)),&
-            N2(MAXVAL(N_PARTICLE_BINS)),N3(MAXVAL(N_PARTICLE_BINS)),&
+REAL(EB) :: KN,KN_FAC,N_I(MAXVAL(N_PARTICLE_BINS)),N0(MAXVAL(N_PARTICLE_BINS)),N1(MAXVAL(N_PARTICLE_BINS)),&
+            N2(MAXVAL(N_PARTICLE_BINS)),N3(MAXVAL(N_PARTICLE_BINS)),DN,&
             RHO_G,TMP_G,MUG,TERMINAL(MAXVAL(N_PARTICLE_BINS)),&
             FU,MOBILITY(MAXVAL(N_PARTICLE_BINS)),ZZ_GET(1:N_TRACKED_SPECIES),AM,AMT(MAXVAL(N_PARTICLE_BINS)),&
             PHI_B,PHI_S,PHI_G,PHI_I,PHI(MAXVAL(N_PARTICLE_BINS),MAXVAL(N_PARTICLE_BINS)),VREL,FU1,FU2,DT_SUBSTEP,DT_SUM,TOL
@@ -287,7 +290,7 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
             ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)
             TMP_G = TMP(I,J,K)
             CALL GET_VISCOSITY(ZZ_GET,MUG,TMP_G)
-            MFP = MUG*SQRT(0.5_EB*PI/(PBAR(K,PRESSURE_ZONE(I,J,K))*RHO_G))
+            KN_FAC = MUG*SQRT(PI/(2._EB*PBAR(K,PRESSURE_ZONE(I,J,K))*RHO_G))
             IM1 = MAX(0,I-1)
             JM1 = MAX(0,J-1)
             KM1 = MAX(0,K-1)
@@ -316,8 +319,9 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
             STRAIN_RATE = 2._EB*(S11**2 + S22**2 + S33**2 + 2._EB*(S12**2 + S13**2 + S23**2))
             DISSIPATION_RATE = MU(I,J,K)/RHO_G*STRAIN_RATE
             N_I = N0
+
             DO N=1,N_PARTICLE_BINS(NS)
-               KN=0.5_EB*MFP/PARTICLE_RADIUS(NS,N)
+               KN=KN_FAC/PARTICLE_RADIUS(NS,N)
                !Verify CN
                MOBILITY(N) = CUNNINGHAM(KN)*MOBILITY_FAC(NS,N)/MUG
                TERMINAL(N) = MOBILITY(N)*GRAV*BIN_X(NS,N)
@@ -326,6 +330,7 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
                         (3._EB*PARTICLE_RADIUS(NS,N)*AM)-PARTICLE_RADIUS(NS,N)
             ENDDO
             PHI = 0._EB
+
             DO N=1,N_PARTICLE_BINS(NS)
                DO NN=1,N_PARTICLE_BINS(NS)
                   IF (NN<N) CYCLE
@@ -356,15 +361,15 @@ SUBROUTINE CALC_AGGLOMERATION(DT,NM)
                   DO NN=N,N_PARTICLE_BINS(NS)
                      IF (N0(N)<MIN_AGGLOMERATION .OR. N0(NN)<MIN_AGGLOMERATION) CYCLE
                      !Remove particles that agglomerate
-                     N1(N)=N1(N)-PHI(NN,N)*N0(N)*N0(NN)*DT_SUBSTEP
-                     IF (NN/=N) N1(NN)=N1(NN)-PHI(NN,N)*N0(N)*N0(NN)*DT_SUBSTEP
+                     DN = PHI(NN,N)*N0(N)*N0(NN)*DT_SUBSTEP
+                     N1(N) = N1(N) - DN
+                     N1(NN) = N1(NN) - DN
                      ! Create new particles from agglomeration
-                     N2(BIN_ETA_INDEX(NS,N,NN,1)) = N2(BIN_ETA_INDEX(NS,N,NN,1)) + &
-                        BIN_ETA(NS,N,NN,1)*PHI(N,NN)*N0(N)*N0(NN)*DT_SUBSTEP
-                     N2(BIN_ETA_INDEX(NS,N,NN,2)) = N2(BIN_ETA_INDEX(NS,N,NN,2)) + &
-                        BIN_ETA(NS,N,NN,2)*PHI(N,NN)*N0(N)*N0(NN)*DT_SUBSTEP
+                     N2(BIN_ETA_INDEX(NS,N,NN,1)) = N2(BIN_ETA_INDEX(NS,N,NN,1)) + BIN_ETA(NS,N,NN,1) * DN
+                     N2(BIN_ETA_INDEX(NS,N,NN,2)) = N2(BIN_ETA_INDEX(NS,N,NN,2)) + BIN_ETA(NS,N,NN,2) * DN
                   ENDDO
                ENDDO AGGLOMERATE_LOOP
+
                N3 = N1 + N2
                N3 = N3 + N0
                TOL = MAXVAL((N0-N3)/(N0+TINY_EB))

Source/wall.f90

@@ -1160,8 +1160,8 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)
          DT_SPYRO(1:SF%N_THICK_REF) = DT * SF%SPYRO_TH_FACTOR(1:SF%N_THICK_REF)
          TSI = MIN(T-B1%T_IGN+DT, SF%REFERENCE_HEAT_FLUX_TIME_INTERVAL+DT)
          IF (SOLID_PHASE_ONLY .AND. .NOT. SF%INERT_Q_REF) THEN
-            B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + &
-               DT*Q_REF_FIT(SUM(B1%M_DOT_G_PP_ACTUAL)*SF%HOC_EFF,SF%HOC_EFF,SF%Y_S_EFF,B1%Q_RAD_IN/B1%EMISSIVITY))/TSI
+            B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + DT * &
+               Q_REF_FIT(SUM(B1%M_DOT_G_PP_ACTUAL)*SF%HOC_EFF,SF%HOC_EFF,SF%Y_S_EFF,B1%Q_RAD_IN/B1%EMISSIVITY))/TSI
          ELSE
             IF (B1%EMISSIVITY > 0._EB) THEN
                B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + DT*(B1%Q_CON_F+B1%Q_RAD_IN/B1%EMISSIVITY))/TSI
@@ -1623,7 +1623,8 @@ SUBROUTINE CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX,CFACE_INDEX)
    TMP_FILM = 0.5_EB*(B1%TMP_G+B1%TMP_F)
    CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_FILM)
    CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP_FILM)
-   KN_FAC = MU_G*SQRT(0.5_EB*PI/(PBAR(BC%KKG,PRESSURE_ZONE(BC%IIG,BC%JJG,BC%KKG))*B1%RHO_G))
+   ! Kn=2 lambda/d, lambda has sqrt(1/2). kn_fac has 2*sqrt(1/2)=sqrt(4/2)=sqrt(2)
+   KN_FAC = MU_G*SQRT(2._EB*PI/(PBAR(BC%KKG,PRESSURE_ZONE(BC%IIG,BC%JJG,BC%KKG))*B1%RHO_G))
    ALPHA = K_G/SM%CONDUCTIVITY_SOLID
    DTMPDX = B1%HEAT_TRANS_COEF*(B1%TMP_G-B1%TMP_F)/K_G
    U_THERM = 0._EB
@@ -1634,7 +1635,6 @@ SUBROUTINE CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX,CFACE_INDEX)
       KN = KN_FAC/SM%THERMOPHORETIC_DIAMETER
       U_THERM = CS2*(ALPHA+CT*KN)*CUNNINGHAM(KN)/((1._EB+CM3*KN)*(1+2*ALPHA+CT2*KN)) * MU_G/(B1%TMP_G*B1%RHO_G)*DTMPDX
    ENDIF
-
    IF (GRAVITATIONAL_DEPOSITION) THEN
       KN = KN_FAC/SM%MEAN_DIAMETER
       U_GRAV = - DOT_PRODUCT(GVEC,BC%NVEC)*CUNNINGHAM(KN)*SM%MEAN_DIAMETER**2*SM%DENSITY_SOLID/(18._EB*MU_G)
@@ -3894,4 +3894,3 @@ SUBROUTINE TGA_ANALYSIS(NM)
 END SUBROUTINE TGA_ANALYSIS
 
 END MODULE WALL_ROUTINES
-

Verification/Aerosols/aerosol_agglomeration.csv

@@ -1,3 +1,3 @@
 Time,Exact Total,Exact Bin 1,Exact Bin 2
 0,0.00001,0.00001,0
-60,0.00001,9.9989E-06,1.06E-09
+60,0.00001,1.00E-05,2.26E-09

Verification/Aerosols/aerosol_agglomeration.fds

@@ -2,7 +2,7 @@
 
 &MESH IJK=10,10,10, XB=0.0,0.01,0.0,0.01,0.00,0.01 /
 
-&TIME T_END=60., DT=0.01/
+&TIME T_END=60.,DT=0.004/
 
 &DUMP FLUSH_FILE_BUFFERS=T, NFRAMES=20 /
 
@@ -13,7 +13,7 @@
       TURBULENT_DEPOSITION     =.FALSE.
       STRATIFICATION           =.FALSE.
       NOISE                    =.FALSE./
-       
+
 &RADI RADIATION=.FALSE./
 
 &SPEC ID='AIR',MW=28.8,CONDUCTIVITY=0.025,VISCOSITY=2.E-5,SPECIFIC_HEAT=1.,BACKGROUND=.TRUE./

Verification/Aerosols/aerosol_agglomeration_2.csv

@@ -1,3 +1,3 @@
 Time,Exact Total,Exact Bin 1,Exact Bin 2
 0,1.00E-05,1.00E-05,0.00E+00
-60,1.00E-05,9.9989E-06,1.08E-09
+60,1.00E-05,1.00E-05,2.30E-09

Verification/Aerosols/aerosol_agglomeration_2.fds

@@ -2,7 +2,7 @@
 
 &MESH IJK=10,10,10, XB=0.0,0.01,0.0,0.01,0.00,0.01 /
 
-&TIME T_END=60./
+&TIME T_END=60.,DT=0.004/
 
 &DUMP FLUSH_FILE_BUFFERS=T, NFRAMES=20 /
 
@@ -14,7 +14,7 @@
       STRATIFICATION           =.FALSE.
       NOISE                    =.FALSE.
       SIMULATION_MODE='DNS'/
-       
+
 &RADI RADIATION=.FALSE./
 
 &SPEC ID='AIR',MW=28.8,CONDUCTIVITY=0.025,VISCOSITY=2.E-5,SPECIFIC_HEAT=1.,BACKGROUND=.TRUE./

Verification/Aerosols/aerosol_gravitational_deposition.csv

@@ -1,6 +1,6 @@
 Time,Exact Gas 1,Exact Gas 2,Exact Wall Top,Exact Wall Bottom
 0,1.20E-05,1.20E-05,0.00E+00,0.00E+00
-0.728,1.20E-05,1.20E-05,0.00E+00,4.79E-08
-0.8,7.21E-06,1.20E-05,0.00E+00,5.26E-08
-0.91,0.00E+00,1.20E-05,0.00E+00,5.99E-08
-1,0.00E+00,6.02E-06,0.00E+00,6.58E-08
+0.72,1.20E-05,1.20E-05,0.00E+00,4.79E-08
+0.8,6.67E-06,1.20E-05,0.00E+00,5.32E-08
+0.9,0.00E+00,1.20E-05,0.00E+00,5.99E-08
+1,0.00E+00,5.51E-06,0.00E+00,6.64E-08

Verification/Aerosols/aerosol_gravitational_deposition_2.csv

@@ -1,6 +1,6 @@
 Time,Exact Gas 1,Exact Gas 2,Exact Wall Top,Exact Wall Bottom
 0,1.20E-05,1.20E-05,0.00E+00,0.00E+00
-0.728,1.20E-05,1.20E-05,4.79E-08,0.00E+00
-0.8,1.20E-05,7.21E-06,5.26E-08,0.00E+00
-0.91,1.20E-05,0.00E+00,5.99E-08,0.00E+00
-1,6.02E-06,0.00E+00,6.58E-08,0.00E+00
+0.72,1.20E-05,1.20E-05,0.00E+00,4.79E-08
+0.8,6.67E-06,1.20E-05,0.00E+00,5.32E-08
+0.9,0.00E+00,1.20E-05,0.00E+00,5.99E-08
+1,0.00E+00,5.51E-06,0.00E+00,6.64E-08

Verification/Aerosols/aerosol_thermophoretic_deposition.csv

@@ -1,3 +1,3 @@
 Time,Exact Gas,Exact Wall
 0,9.04E-06,0.00E+00
-2,6.01E-06,3.69E-09
+2,5.07E-06,5.05E-09