import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
def K12(Dp1, Dp2):
T = 298.15
ρ = 1000
kb = 1.38064852E-23
λ = 0.0686E-6
μ = 1.83E-5
m1 = np.pi/6 * Dp1**3 * ρ
m2 = np.pi/6 * Dp2**3 * ρ
c1 = np.sqrt(8*kb*T/np.pi/m1)
c2 = np.sqrt(8*kb*T/np.pi/m2)
Kn1 = 2*λ/Dp1
Kn2 = 2*λ/Dp2
Cc1 = 1 + Kn1*(1.257 + 0.4*np.exp(-1.1/Kn1))
Cc2 = 1 + Kn2*(1.257 + 0.4*np.exp(-1.1/Kn2))
D1 = kb*T*Cc1/(3*np.pi*μ*Dp1)
D2 = kb*T*Cc2/(3*np.pi*μ*Dp2)
l1 = 8*D1/np.pi/c1
l2 = 8*D2/np.pi/c2
g1 = np.sqrt(2)/3/Dp1/l1*((Dp1+l1)**3 - (Dp1**2+l1**2)**1.5) - np.sqrt(2)*Dp1 # note, Seinfeld Table 13.1 is missing the second root 2, needed for Fig. 13.5
g2 = np.sqrt(2)/3/Dp2/l2*((Dp2+l2)**3 - (Dp2**2+l2**2)**1.5) - np.sqrt(2)*Dp2 # with 1.0 instead of root 2 in both places, the max error in the D=D curve is 5%
β = ((Dp1+Dp2)/(Dp1+Dp2+2*np.sqrt(g1**2+g2**2)) + 8*(D1+D2)/(Dp1+Dp2)/np.sqrt(c1**2+c2**2))**-1
K12 = 2*np.pi*(D1+D2)*(Dp1+Dp2) * β
return K12*1E6
############################
Dp = np.logspace(-3, 1, 100) * 1E-6
yy = np.empty(len(Dp))
for i in range(len(yy)):
yy[i] = K12(Dp[i], Dp[i])
Dp10 = Dp[Dp<10E-6]
Dp1 = Dp[Dp<1E-6]
Dp01 = Dp[Dp<0.1E-6]
Dp001 = Dp[Dp<0.01E-6]
y10 = K12(Dp10, 10E-6)
y1 = K12(Dp1, 1E-6)
y01 = K12(Dp01, 0.1E-6)
y001 = K12(Dp001, 0.01E-6)
fig, ax = plt.subplots(figsize=(5, 7))
plt.plot(Dp10 *1E6, y10, 'k-')
plt.plot(Dp1 *1E6, y1, 'b-')
plt.plot(Dp01 *1E6, y01, 'g-')
plt.plot(Dp001*1E6, y001, 'r-')
plt.plot(Dp*1E6, yy, '-', color='orange')
plt.yscale('log')
plt.xscale('log')
plt.ylim([1E-10, 1E-3])
plt.xlim([1E-3, 10]);
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend([r'D$_{p2}$=10 $\mu$m', r'D$_{p2}$=1 $\mu$m', r'D$_{p2}$=0.1 $\mu$m', r'D$_{p2}$=0.01 $\mu$m', r'D$_{p1}$=D$_{p2}$'],frameon=False);
Note, this is a particle Knudsen number based on the particle mean free path, not a gas Knudsen number based on the gas mean free path.
Note that the calculation of $K_{12}$ involves two particles and the system is symmetric.
Below, we plot the minimum $Kn$ and show vertical lines at $Kn=0.1$ and $Kn=1$, where we appear to depart from free molecular, and continuum, respectively.
In the figures, three values of $D_{p1}$ are shown, with plots of $K_{12}$ for a wide range of $D_{p2}$, including values of $D_{p2}>D_{p1}$, unlike the more standard representations. Also shown are the curves for the free molecular and continuum limits.
Note that a gas-phase Knudsen number does not make sense in the context of the diffusion equations used to describe these processes. See Seinfeld Chapter 13.3, and especially section 13.3.1.2.
When the mean free path λp of the diffusing aerosol particle is comparable to the radius of the absorbing particle, the boundary condition at the absorbing particle surface must be corrected to account for the nature of the diffusion process in the vicinity of the surface. The physical meaning of this can be explained as follows. As we have seen, the diffusion equations can be applied to Brownian motion only for time intervals that are large compared to the relaxation time τ of the particles or for distances that are large compared to the aerosol mean free path λp. Diffusion equations cannot describe the motion of particles inside a layer of thickness λp adjacent to an absorbing wall. If the size of the absorbing sphere is comparable to λp, this layer has a substantial effect on the kinetics of coagulation.
Frenklach suggests using the harmonic mean in the transition regime. $$HM = \frac{C\cdot FM}{C + FM}.$$
def K12(Dp1, Dp2):
T = 298.15
ρ = 1000
kb = 1.38064852E-23
λ = 0.0686E-6
μ = 1.83E-5
m1 = np.pi/6 * Dp1**3 * ρ
m2 = np.pi/6 * Dp2**3 * ρ
c1 = np.sqrt(8*kb*T/np.pi/m1)
c2 = np.sqrt(8*kb*T/np.pi/m2)
Kn1 = 2*λ/Dp1
Kn2 = 2*λ/Dp2
Cc1 = 1 + Kn1*(1.257 + 0.4*np.exp(-1.1/Kn1))
Cc2 = 1 + Kn2*(1.257 + 0.4*np.exp(-1.1/Kn2))
D1 = kb*T*Cc1/(3*np.pi*μ*Dp1)
D2 = kb*T*Cc2/(3*np.pi*μ*Dp2)
l1 = 8*D1/np.pi/c1
l2 = 8*D2/np.pi/c2
g1 = np.sqrt(2)/3/Dp1/l1*((Dp1+l1)**3 - (Dp1**2+l1**2)**1.5) - np.sqrt(2)*Dp1 # note, Seinfeld Table 13.1 is missing the second root 2, needed for Fig. 13.5
g2 = np.sqrt(2)/3/Dp2/l2*((Dp2+l2)**3 - (Dp2**2+l2**2)**1.5) - np.sqrt(2)*Dp2 # with 1.0 instead of root 2 in both places, the max error in the D=D curve is 5%
β = ((Dp1+Dp2)/(Dp1+Dp2+2*np.sqrt(g1**2+g2**2)) + 8*(D1+D2)/(Dp1+Dp2)/np.sqrt(c1**2+c2**2))**-1
K12 = 2*np.pi*(D1+D2)*(Dp1+Dp2) * β
K12c = 2*np.pi*(D1+D2)*(Dp1+Dp2)
K12fm = np.pi/4*(Dp1 + Dp2)**2*np.sqrt(c1**2+c2**2)
HM = K12c*K12fm/(K12c + K12fm)
return K12*1E6, K12c*1E6, K12fm*1E6, HM*1E6, l1/Dp2, l2/Dp1
#################
Dp = np.logspace(-8, 3, 5000) * 1E-6
yy = np.empty(len(Dp))
y1 = np.empty(len(Dp))
y1c = np.empty(len(Dp))
y1fm = np.empty(len(Dp))
HM1 = np.empty(len(Dp))
y04 = np.empty(len(Dp))
y04c = np.empty(len(Dp))
y04fm = np.empty(len(Dp))
HM04 = np.empty(len(Dp))
y001 = np.empty(len(Dp))
y001c = np.empty(len(Dp))
y001fm = np.empty(len(Dp))
HM001 = np.empty(len(Dp))
r1 = np.empty(len(Dp))
r2 = np.empty(len(Dp))
rr1 = np.empty(len(Dp))
rr04 = np.empty(len(Dp))
rr001 = np.empty(len(Dp))
for i in range(len(yy)):
yy[i],dmb1,dmb2,dmb3,dmb4,dmb5 = K12(Dp[i], Dp[i])
y1[i] , y1c[i] , y1fm[i], HM1[i], r1[i], r2[i] = K12(Dp[i], 1.001E-6)
rr1[i] = np.min((r1[i], r2[i]))
y04[i] , y04c[i] , y04fm[i], HM04[i], r1[i], r2[i] = K12(Dp[i], 0.040E-6)
rr04[i] = np.min((r1[i], r2[i]))
y001[i] , y001c[i] , y001fm[i], HM001[i], r1[i], r2[i] = K12(Dp[i], 0.001E-6)
rr001[i] = np.min((r1[i], r2[i]))
ir1_1 = np.abs(rr1-1.0).argmin()
ir1_01 = np.abs(rr1-0.1).argmin()
ir04_1 = np.abs(rr04-1.0).argmin()
ir04_01 = np.abs(rr04-0.1).argmin()
ir001_1 = np.abs(rr001-1.0).argmin()
ir001_01 = np.abs(rr001-0.1).argmin()
#------------------
fig, ax = plt.subplots(figsize=(5, 5))
plt.plot(Dp*1E6, y1, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y1c, 'k-.')
plt.plot(Dp*1E6, y1fm, 'k--')
plt.plot(Dp*1E6, HM1, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr1, 'r:')
plt.plot(np.array([Dp[ir1_1], Dp[ir1_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir1_01], Dp[ir1_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
plt.ylim([1E-11, 1E3])
plt.xlim([1E-8, 1E4])
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'HM', 'Dp1=Dp2', r'($\lambda_{p1}$/D$_{p2}$)$_{min}$'], frameon=False)
plt.title(r'D$_{p,1}$ = 1 $\mu$m');
#------------------
fig, ax = plt.subplots(figsize=(5, 5))
plt.plot(Dp*1E6, y04, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y04c, 'k-.')
plt.plot(Dp*1E6, y04fm, 'k--')
plt.plot(Dp*1E6, HM04, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr04, 'r:')
plt.plot(np.array([Dp[ir04_1], Dp[ir04_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir04_01], Dp[ir04_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
plt.ylim([1E-11, 1E3])
plt.xlim([1E-8, 1E4])
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'HM', 'Dp1=Dp2', r'($\lambda_{p1}$/D$_{p2}$)$_{min}$'], frameon=False)
plt.title(r'D$_{p,1}$ = 0.04 $\mu$m');
#------------------
fig, ax = plt.subplots(figsize=(5, 5))
plt.plot(Dp*1E6, y001, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y001c, 'k-.')
plt.plot(Dp*1E6, y001fm, 'k--')
plt.plot(Dp*1E6, HM001, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr001, 'r:')
plt.plot(np.array([Dp[ir001_1], Dp[ir001_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir001_01], Dp[ir001_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
plt.ylim([1E-11, 1E3])
plt.xlim([1E-8, 1E4]);
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'HM', 'Dp1=Dp2', r'($\lambda_{p1}$/D$_{p2}$)$_{min}$'], frameon=False)
plt.title(r'D$_{p,1}$ = 0.001 $\mu$m');
#------------------
fig, ax = plt.subplots(figsize=(4, 4))
plt.plot(Dp*1E6, y1, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y1c, 'k-.')
plt.plot(Dp*1E6, y1fm, 'k--')
plt.plot(Dp*1E6, HM1, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr1, 'r:')
plt.plot(np.array([Dp[ir1_1], Dp[ir1_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir1_01], Dp[ir1_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
plt.ylim([1E-5, 1E1])
plt.xlim([1E-6, 1E-2])
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'HM', 'Dp1=Dp2', r'($\lambda_{p1}$/D$_{p2}$)$_{min}$'], frameon=False)
plt.title(r'D$_{p,1}$ = 1 $\mu$m');
#------------------
fig, ax = plt.subplots(figsize=(4, 4))
plt.plot(Dp*1E6, y04, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y04c, 'k-.')
plt.plot(Dp*1E6, y04fm, 'k--')
plt.plot(Dp*1E6, HM04, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr04, 'r:')
plt.plot(np.array([Dp[ir04_1], Dp[ir04_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir04_01], Dp[ir04_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
plt.ylim([1E-9, 1E-7])
plt.xlim([1E-3, 1E0])
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'HM', 'Dp1=Dp2', r'($\lambda_{p1}$/D$_{p2}$)$_{min}$'], frameon=False)
plt.title(r'D$_{p,1}$ = 0.04 $\mu$m');
#------------------
fig, ax = plt.subplots(figsize=(4, 4))
plt.plot(Dp*1E6, y001, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y001c, 'k-.')
plt.plot(Dp*1E6, y001fm, 'k--')
plt.plot(Dp*1E6, HM001, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr001, 'r:')
plt.plot(np.array([Dp[ir001_1], Dp[ir001_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir001_01], Dp[ir001_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
plt.ylim([1E-7, 1E-3])
plt.xlim([1E-2, 1E1]);
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'HM', 'Dp1=Dp2', r'($\lambda_{p1}$/D$_{p2}$)$_{min}$'], frameon=False)
plt.title(r'D$_{p,1}$ = 0.001 $\mu$m');
def K12(Dp1, Dp2):
T = 298.15
ρ = 1000
kb = 1.38064852E-23
λ = 0.0686E-6
μ = 1.83E-5
m1 = np.pi/6 * Dp1**3 * ρ
m2 = np.pi/6 * Dp2**3 * ρ
c1 = np.sqrt(8*kb*T/np.pi/m1)
c2 = np.sqrt(8*kb*T/np.pi/m2)
Kn1 = 2*λ/Dp1
Kn2 = 2*λ/Dp2
Cc1 = 1 + Kn1*(1.257 + 0.4*np.exp(-1.1/Kn1))
Cc2 = 1 + Kn2*(1.257 + 0.4*np.exp(-1.1/Kn2))
D1 = kb*T*Cc1/(3*np.pi*μ*Dp1)
D2 = kb*T*Cc2/(3*np.pi*μ*Dp2)
l1 = 8*D1/np.pi/c1
l2 = 8*D2/np.pi/c2
g1 = np.sqrt(2)/3/Dp1/l1*((Dp1+l1)**3 - (Dp1**2+l1**2)**1.5) - np.sqrt(2)*Dp1 # note, Seinfeld Table 13.1 is missing the second root 2, needed for Fig. 13.5
g2 = np.sqrt(2)/3/Dp2/l2*((Dp2+l2)**3 - (Dp2**2+l2**2)**1.5) - np.sqrt(2)*Dp2 # with 1.0 instead of root 2 in both places, the max error in the D=D curve is 5%
Knp = np.min([l1/Dp2, l2/Dp1])
β = ((Dp1+Dp2)/(Dp1+Dp2+2*np.sqrt(g1**2+g2**2)) + 8*(D1+D2)/(Dp1+Dp2)/np.sqrt(c1**2+c2**2))**-1
K12 = 2*np.pi*(D1+D2)*(Dp1+Dp2) * β
K12c = 2*np.pi*(D1+D2)*(Dp1+Dp2)
K12fm = np.pi/4*(Dp1 + Dp2)**2*np.sqrt(c1**2+c2**2)
TR = K12c/(1+Knp) + K12fm/(1+1/Knp)
return K12*1E6, K12c*1E6, K12fm*1E6, TR*1E6, Knp
########################
Dp = np.logspace(-8, 3, 5000) * 1E-6
yy = np.empty(len(Dp))
y1 = np.empty(len(Dp))
y1c = np.empty(len(Dp))
y1fm = np.empty(len(Dp))
TR1 = np.empty(len(Dp))
y04 = np.empty(len(Dp))
y04c = np.empty(len(Dp))
y04fm = np.empty(len(Dp))
TR04 = np.empty(len(Dp))
y001 = np.empty(len(Dp))
y001c = np.empty(len(Dp))
y001fm = np.empty(len(Dp))
TR001 = np.empty(len(Dp))
rr1 = np.empty(len(Dp))
rr04 = np.empty(len(Dp))
rr001 = np.empty(len(Dp))
for i in range(len(yy)):
yy[i],dmb1,dmb2,dmb3,dmb4 = K12(Dp[i], Dp[i])
y1[i], y1c[i], y1fm[i], TR1[i], rr1[i] = K12(Dp[i], 1.001E-6)
y04[i], y04c[i], y04fm[i], TR04[i], rr04[i] = K12(Dp[i], 0.040E-6)
y001[i], y001c[i], y001fm[i], TR001[i], rr001[i] = K12(Dp[i], 0.001E-6)
#------------------
fig, ax = plt.subplots(figsize=(4, 4))
plt.plot(Dp*1E6, y1, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y1c, 'k-.')
plt.plot(Dp*1E6, y1fm, 'k--')
plt.plot(Dp*1E6, TR1, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr1, 'r:')
plt.plot(np.array([Dp[ir1_1], Dp[ir1_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir1_01], Dp[ir1_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
#plt.ylim([1E-5, 1E1])
#plt.xlim([1E-6, 1E-2])
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'TR', 'Kn', 'Dp1=Dp2'], frameon=False)
plt.title(r'D$_{p,1}$ = 1 $\mu$m');
#------------------
fig, ax = plt.subplots(figsize=(4, 4))
plt.plot(Dp*1E6, y04, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y04c, 'k-.')
plt.plot(Dp*1E6, y04fm, 'k--')
plt.plot(Dp*1E6, TR04, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr04, 'r:')
plt.plot(np.array([Dp[ir04_1], Dp[ir04_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir04_01], Dp[ir04_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
#plt.ylim([1E-9, 1E-7])
#plt.xlim([1E-3, 1E0])
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'TR', 'Kn', 'Dp1=Dp2'], frameon=False)
plt.title(r'D$_{p,1}$ = 0.04 $\mu$m');
#------------------
fig, ax = plt.subplots(figsize=(4, 4))
plt.plot(Dp*1E6, y001, '-', linewidth=8, color='lightgray')
plt.plot(Dp*1E6, y001c, 'k-.')
plt.plot(Dp*1E6, y001fm, 'k--')
plt.plot(Dp*1E6, TR001, 'k:')
plt.plot(Dp*1E6, yy)
plt.plot(Dp*1E6, rr001, 'r:')
plt.plot(np.array([Dp[ir001_1], Dp[ir001_1]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.plot(np.array([Dp[ir001_01], Dp[ir001_01]])*1E6, [1E-12, 1E4], ':', color='gray')
plt.yscale('log')
plt.xscale('log')
#plt.ylim([1E-7, 1E-3])
#plt.xlim([1E-2, 1E1]);
plt.xlabel(r'D $\mu$m')
plt.ylabel(r'K$_{12}$ cm$^3$/s')
plt.legend(['K11', 'Cont.', 'FM', 'TR', 'Kn', 'Dp1=Dp2'], frameon=False)
plt.title(r'D$_{p,1}$ = 0.001 $\mu$m');