diff --git a/Manuals/FDS_Technical_Reference_Guide/Aerosol_Chapter.tex b/Manuals/FDS_Technical_Reference_Guide/Aerosol_Chapter.tex index 814585c921a..4f835f175d2 100644 --- a/Manuals/FDS_Technical_Reference_Guide/Aerosol_Chapter.tex +++ b/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. diff --git a/Source/func.f90 b/Source/func.f90 index bfe47c4e2e9..f0e51205646 100644 --- a/Source/func.f90 +++ b/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 - + \ No newline at end of file diff --git a/Source/soot.f90 b/Source/soot.f90 index 735aba1a278..a836f76aa9a 100644 --- a/Source/soot.f90 +++ b/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,7 +237,7 @@ 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 @@ -242,13 +245,13 @@ SUBROUTINE INITIALIZE_AGGLOMERATION 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 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 - diff --git a/Verification/Aerosols/aerosol_agglomeration.csv b/Verification/Aerosols/aerosol_agglomeration.csv index 47a1d9957e2..eaab85ba68a 100644 --- a/Verification/Aerosols/aerosol_agglomeration.csv +++ b/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 diff --git a/Verification/Aerosols/aerosol_agglomeration.fds b/Verification/Aerosols/aerosol_agglomeration.fds index b4a91e9008f..1bd11f2a418 100644 --- a/Verification/Aerosols/aerosol_agglomeration.fds +++ b/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./ diff --git a/Verification/Aerosols/aerosol_agglomeration_2.csv b/Verification/Aerosols/aerosol_agglomeration_2.csv index 49f1ed811fb..3d515d1066e 100644 --- a/Verification/Aerosols/aerosol_agglomeration_2.csv +++ b/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 diff --git a/Verification/Aerosols/aerosol_agglomeration_2.fds b/Verification/Aerosols/aerosol_agglomeration_2.fds index 58c5dbe9c3b..14a89b1a19f 100644 --- a/Verification/Aerosols/aerosol_agglomeration_2.fds +++ b/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./ diff --git a/Verification/Aerosols/aerosol_gravitational_deposition.csv b/Verification/Aerosols/aerosol_gravitational_deposition.csv index 29fefacbd4d..e5328a5554c 100644 --- a/Verification/Aerosols/aerosol_gravitational_deposition.csv +++ b/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 diff --git a/Verification/Aerosols/aerosol_gravitational_deposition_2.csv b/Verification/Aerosols/aerosol_gravitational_deposition_2.csv index 424d6af13c3..e5328a5554c 100644 --- a/Verification/Aerosols/aerosol_gravitational_deposition_2.csv +++ b/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 diff --git a/Verification/Aerosols/aerosol_thermophoretic_deposition.csv b/Verification/Aerosols/aerosol_thermophoretic_deposition.csv index fa48d0bf2e3..6b6eb54f549 100644 --- a/Verification/Aerosols/aerosol_thermophoretic_deposition.csv +++ b/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