Skip to content

Commit

Permalink
Merge pull request #13947 from drjfloyd/master
Browse files Browse the repository at this point in the history
FDS Source: Fix Kn number calculation and associated verification Issue #13926
  • Loading branch information
drjfloyd authored Dec 25, 2024
2 parents 287cc5e + 303b752 commit 61efeed
Show file tree
Hide file tree
Showing 11 changed files with 48 additions and 44 deletions.
6 changes: 3 additions & 3 deletions Manuals/FDS_Technical_Reference_Guide/Aerosol_Chapter.tex
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
2 changes: 1 addition & 1 deletion Source/func.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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

45 changes: 25 additions & 20 deletions Source/soot.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)))
Expand Down Expand Up @@ -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))
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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))
Expand Down
9 changes: 4 additions & 5 deletions Source/wall.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -3894,4 +3894,3 @@ SUBROUTINE TGA_ANALYSIS(NM)
END SUBROUTINE TGA_ANALYSIS

END MODULE WALL_ROUTINES

2 changes: 1 addition & 1 deletion Verification/Aerosols/aerosol_agglomeration.csv
Original file line number Diff line number Diff line change
@@ -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
4 changes: 2 additions & 2 deletions Verification/Aerosols/aerosol_agglomeration.fds
Original file line number Diff line number Diff line change
Expand Up @@ -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 /

Expand All @@ -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./
Expand Down
2 changes: 1 addition & 1 deletion Verification/Aerosols/aerosol_agglomeration_2.csv
Original file line number Diff line number Diff line change
@@ -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
4 changes: 2 additions & 2 deletions Verification/Aerosols/aerosol_agglomeration_2.fds
Original file line number Diff line number Diff line change
Expand Up @@ -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 /

Expand All @@ -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./
Expand Down
8 changes: 4 additions & 4 deletions Verification/Aerosols/aerosol_gravitational_deposition.csv
Original file line number Diff line number Diff line change
@@ -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
8 changes: 4 additions & 4 deletions Verification/Aerosols/aerosol_gravitational_deposition_2.csv
Original file line number Diff line number Diff line change
@@ -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
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 61efeed

Please sign in to comment.