在高分子聚集、湿法造粒、纳米颗粒合成等许多过程中，颗粒凝并/团聚是颗粒体积增大的最主要物理机制之一^[1-5]，并且通常发生于不同化学组分的颗粒之间，颗粒内部组分分布对于颗粒物理化学性能均有重要甚至决定性的影响^[6-8]，因此对于多组分颗粒凝并及组分分布的研究十分重要。对于初始双分散的颗粒群的凝并过程，组分分布可以通过高斯函数来描述：组分A的组分分布以概率密度函数g(m_A|m) dm_A表示，代表在质量为m的颗粒中组分A含量在质量区间 (m_A, m_A+dm_A) 内的颗粒所占比例。基于随机混合理论和中心限制定理，已证明概率密度函数g(m_A|m) 满足公式 (1) 所示的高斯函数的形式^[9-11]，并且该结论也被颗粒群平衡模拟中的常数目法和异权值Monte Carlo (MC) 方法验证^[11-12]

$ g\left( {{m_{\text{A}}}\left| {m,t} \right.} \right) = \frac{1}{{\sqrt {2{\text{π}}m\chi } }}\exp \left[ { - \frac{{{{\left( {{m_{\text{A}}} - \phi m} \right)}^2}}}{{2m\chi }}} \right], $

(1)

其中φ是组分A的质量分数，随着凝并过程的进行保持不变

$ \phi = \frac{{{M_{\text{A}}}}}{{{M_{\text{A}}} + {M_{\text{B}}}}}, $

(2)

组分A的初始质量M_A=N_A0m_A0。

从公式 (1) 可以看到，组分分布函数g(m_A|m) 取决于混合程度χ(量化描述组分A分散程度的无量纲参数，未知变量) 和组分A的总质量分数φ(已知变量)。因此一旦能预测混合程度，即可依据已有理论预测颗粒内部组分分布。

而颗粒群的混合程度由χ表示^[7], 如下定义:

$ \begin{array}{*{20}{c}} {\chi = \frac{{{X^2}}}{M} = \frac{{\int\limits_0^\infty {{\text{d}}{m_{\text{A}}}} \int\limits_0^\infty {{\text{d}}{m_{\text{B}}}{x^2}n\left( {{m_{\text{A}}},{m_{\text{B}}},t} \right)} }}{{\int\limits_0^\infty {{\text{d}}{m_{\text{A}}}} \int\limits_0^\infty {{\text{d}}{m_{\text{B}}}\left( {{m_{\text{A}}} + {m_{\text{B}}}} \right)n\left( {{m_{\text{A}}},{m_{\text{B}}},t} \right)} }}} \\ { = \frac{{\int\limits_0^\infty {{\text{d}}m} \int\limits_0^m {{\text{d}}{m_{\text{A}}}{x^2}f\left( {m,t} \right)g\left( {{m_{\text{A}}}\left| {m,t} \right.} \right)} }}{{\int\limits_0^\infty {{\text{d}}m} \int\limits_0^m {{\text{d}}{m_{\text{A}}}mf\left( {m,t} \right)g\left( {{m_{\text{A}}}\left| {m,t} \right.} \right)} }}.} \end{array} $

(3)

其中：X²量化描述双组分混合程度；M为颗粒总质量 (M=M_A+M_B)；n(m_A, m_B, t) 是时刻t的数密度函数；n(m_A, m_B, t) dm_Adm_B代表组分A质量在m_A到m_A+dm_A之间，组分B质量在m_B到m_B+dm_B之间的颗粒数目浓度；f(m，t) 是颗粒尺度分布函数；x为组分A对于平均含量φm的过余量：

$ x = {m_{\text{A}}} - \phi m = {m_{\text{A}}} - \phi \left( {{m_{\text{A}}} + {m_{\text{B}}}} \right). $

(4)

χ值越小，即双组分混合得越好，有学者发现随凝并过程进行，χ总会达到一个稳定值^{[11, 13]}，而我们发现在χ稳定时，达到自保持分布状态^[12]，并且可进一步根据初始工况，确定出稳态χ的值^[14]，而本文关注如何通过初始状态预测出χ在到达稳态之前随时间的演变过程。

1 数值方法与计算工况

这里采用颗粒群平衡随机模拟的方法 (PBMC)，它能够直接模拟大量颗粒。但因为传统的MC方法对凝并问题的计算代价与模拟颗粒数目N_s的平方成正比，而计算精度与N_s的平方根成正比，异权值MC方法 (DWMC) 能够协调每个尺度区间内模拟颗粒数目，使其保持在一个适当的范围内，在保证足够计算精度的同时使计算更加高效^[15-19]。并且为了避免双重遍历，在模拟颗粒凝并时，采用快速DWMC方法加速PBMC模拟^[20-21]。具体数值方法详见参考文献[14]。

如表 1所示，我们考虑3种凝并核^[21]，且本文仅考虑组分无关工况即双组分密度相同，具体计算工况如表 2所示。3种工况均为经典的模拟工况，线性凝并为一种理想工况，而自由分子区布朗凝并。连续区布朗凝并为实际工况，适用于日常生活中常见的气溶胶系统研究如尘、烟、雾、霾等现象，制药工业、纳米粉体制造业、矿物燃烧过程中超细颗粒物生成。布朗凝并表示凝并过程是由由布朗运动引起的颗粒间的相对运动导致的。我们通过克努德森数来区分 (Kn=2λ/d_p，λ为气体分子平均自由程，即相邻两次碰撞之间气体分子移动的平均距离，d_p为颗粒直径)。亚微米颗粒 (Kn>10) 处于自由分子区，这时气体分子与颗粒表面的碰撞占据主导地位，气溶胶颗粒的行为接近于自由运动的气体分子，可以使用分子运动论描述气溶胶颗粒间的碰撞凝并现象。微米尺度颗粒或者超微米颗粒 (Kn < 0.1) 处于连续区，此时假设气体为连续介质，与颗粒表面接触处气体速度为0，适用于纳维-斯托克斯方程。

2 结果与讨论 2.1 混合程度的相关参数

在整个演变过程中，χ(t)/χ₀一直保持相同的数量级，考虑初始工况凝并核函数中的相关参数，我们通过控制变量法已证温度T、初始数目N_A0、粒径d_A0、和密度ρ_A对混合程度没有影响^[14]。在此基础上，为探讨实时的χ(t)/χ₀与初始工况的关系，引入如下参数

$ \alpha = \frac{{{N_{{\text{A}}0}}}}{{{N_{{\text{B}}0}}}},\beta = \frac{{{d_{{\text{A}}0}}}}{{{d_{{\text{B}}0}}}}, $

(5)

$ \varphi = \frac{{{N_{{\text{A}}0}}}}{{{N_{{\text{A}}0}} + {N_{{\text{B}}0}}}} = \frac{\alpha }{{1 + \alpha }}, $

(6)

$ \phi = \frac{{{M_{{\text{A}}0}}}}{{{M_{{\text{A}}0}} + {M_{{\text{B}}0}}}} = \frac{{\alpha {\beta ^3}}}{{1 + \alpha {\beta ^3}}}. $

(7)

φ、φ别为初始时刻组分A的数目浓度和总质量分数。

初始混合程度χ₀^[9-11]如下表示：

$ \begin{array}{*{20}{c}} {{\chi _0} = \left( {1 - \phi } \right)\phi \left( {{m_{{\text{A}}0}}\left( {1 - \phi } \right) + {m_{{\text{B}}0}}\phi } \right) = } \\ {\frac{{\alpha {\beta ^3}\gamma \left( {1 + \alpha } \right)}}{{{{\left( {1 + \alpha {\beta ^3}\gamma } \right)}^3}}}{m_{{\text{A}}0}} = \frac{{\phi {{\left( {1 - \phi } \right)}^2}}}{{\left( {1 - \phi } \right)}}{m_{{\text{A}}0}}.} \end{array} $

(8)

基于大量MC模拟数据，在初始双分散工况下，发现χ随无量纲时间τ的变化规律，即χ(t)/χ₀满足指数增长函数^[14]

$ \chi \left( t \right)/{\chi _0} = {\chi _\infty }/{\chi _0} + {C_1}\exp \left( { - t/{C_2}} \right). $

(9)

其中，C₂为时间常数，和变化速率相关：

$ \frac{{{\text{d}}\frac{{\chi \left( t \right)}}{{{\chi _0}}}}}{{{\text{d}}t}} = - \frac{1}{{{C_2}}}\frac{{\chi \left( t \right)}}{{{\chi _0}}}. $

(10)

而如图 1所示，当稳定状态的χ相同时，C₂并不相同，所以还与其他因素有关。

首先，根据初始状态t=0，得

$ {C_1} + {\chi _\infty }/{\chi _0} = 1. $

(11)

2.2 自由分子区布朗凝并

已证自由分子区稳定状态的χ^[11]满足

$ {\chi _\infty }/{\chi _0} = \exp \left[ { - 2{{\left( {\phi - {\phi _{\beta = 1}}} \right)}^2}} \right]. $

(12)

当γ=1，α=0.1，β变化时, 发现C₂=f(χ_∞/χ₀)，且经快速DWMC方法模拟结果拟合，该函数满足如下所示的线性关系：

$ {C_2} = {C_3}\ln \left( {{\chi _\infty }/{\chi _0}} \right) + {C_4}. $

(13)

通过7 × 19种工况验证，其中α=0.03, 0.067, 0.1, 0.2, 0.3, 0.43, 0.67；而φ分别取{0.05, 0.1, 0.15, 0.2, …, 0.90, 0.95}，发现对于不同α，斜率C₃基本不变，如图 2所示。经具体拟合每个α所对应的线性关系，得到C₃=2.5。

而关于每一条直线截距满足C₄=f(α)， 具体计算结果如图 3，C₄符合如下关系：

$ {C_4} = 2 \times {\left( {\alpha + \frac{1}{\alpha }} \right)^{0.44}}. $

(14)

综上所述，χ_∞/χ₀作为φ的函数可以表示成如下形式

$ \frac{{\chi \left( t \right)}}{{{\chi _0}}} = \left( {1 - \frac{{{\chi _\infty }}}{{{\chi _0}}}} \right)\exp \left\{ { - t/\left[ \begin{gathered} 2.5\ln \left( {\frac{{{\chi _\infty }}}{{{\chi _0}}}} \right) + \hfill \\ 2{\left( {\alpha + 1/\alpha } \right)^{0.44}} \hfill \\ \end{gathered} \right]} \right\} + \frac{{{\chi _\infty }}}{{{\chi _0}}}. $

(15)

根据上述稳态参数与初始参数的函数关系，可得

$ \begin{array}{*{20}{c}} {\frac{{\chi \left( t \right)}}{{{\chi _0}}} = \left\{ {1 - \exp \left[ { - 2{{\left( {\phi - {\phi _{\beta = 1}}} \right)}^2}} \right]} \right\}} \\ {\exp \left\{ { - t/\left[ { - 5{{\left( {\phi - {\phi _{\beta = 1}}} \right)}^2} + 2{{\left( {\alpha + 1/\alpha } \right)}^{0.44}}} \right]} \right\} + } \\ {\exp \left[ { - 2{{\left( {\phi - {\phi _{\beta = 1}}} \right)}^2}} \right].} \end{array} $

(16)

对于该形式的预测公式，进行进一步验证，以α为定值10，β为变化为例，φ分别为0.077, 0.2, 0.345, 0.46, 0.58和0.69。在初始双分散工况下，发现公式所预测的结果同计算结果的拟合曲线吻合良好，如图 4所示。

混合程度的实时演变χ(t)/χ₀为一单调递减过程，并且下降的速率逐渐变慢，随时间推进越来越接近稳态值。在多种初等函数中，我们发现χ可以满足公式 (9)，χ(t)/χ₀以一个与当前值成正比的速度下降$\text{d}\frac{\chi \left( t \right)}{{{\chi }_{0}}}/\text{d}t\propto-\frac{1}{{{C}_{2}}}\frac{\chi \left( t \right)}{{{\chi }_{0}}} $， 这符合模拟方法中的马尔可夫过程，即下一步的状态只和现在状态有关而同之前状态无关。C₂则为这个指数形式的决定性参数，因此我们进一步讨论C₂的准确性，如图 5所示。

关于该公式的误差，在同时改变φ和φ(分别取{0.1, 0.2, …, 0.90}) 时共9×9个工况，对于指数函数速率时间常数C₂公式预测值与MC计算结果的相对误差在0.2以内，如图 5所示。误差产生的原因主要有3点，首先作为经验公式，这种指数形式对于实时χ的预测公式是一种规律的总结，此处会引入一定程度误差，并且常数目的MC方法单次模拟也会与平均值有±5%的波动。其次自由分子区与连续区工况，对于稳定状态的χ的公式预测值与模拟值有6.5%的误差^[14]。同时在图 4中，当β=0.441时，C₂的公式预测值3.94与MC模拟值4.53之比为1.15，此时图 4混合程度随时间演变的方差为2.96×10^-4，因此我们认为图 5的分散度可以接受。

2.3 连续区布朗凝并

为了验证χ(t)/χ₀公式的适用范围，我们还研究了连续区布朗凝并，由于该核函数考虑颗粒的体积而非质量，所以本质上满足组分无关性。采用DWMC模拟25个工况 (表 2工况2) 发现满足公式 (17) 的形式，误差范围同图 5所示。

$ \begin{array}{*{20}{c}} {\frac{{\chi \left( t \right)}}{{{\chi _0}}} = \left\{ {1 - \exp \left[ { - {{\left( {\phi - {\phi _{\beta = 1}}} \right)}^2}/\sqrt 2 } \right]} \right\}} \\ {\exp \left\{ { - t/\left[ { - 10.6{{\left( {\phi - {\phi _{\beta = 1}}} \right)}^2} + {{\left( {\alpha + 1/\alpha } \right)}^{0.75}} + 6} \right]} \right\} + } \\ {\exp \left[ { - {{\left( {\phi - {\phi _{\beta = 1}}} \right)}^2}/\sqrt 2 } \right].} \end{array} $

(17)

2.4 线性凝并

线性凝并核工况表示两颗粒的凝并概率正比于两颗粒的体积总和，通常可用于近似描述湍流重力机制作用下的颗粒凝并行为。但是当采用线性凝并核 (表 2工况3) 模拟时并不满足指数函数关系。从理论上考虑，因为该核函数满足，$k\left( {{m}_{\text{A}}}, {{m}_{\text{B}}};m{{'}_{\text{A}}}, m{{'}_{\text{B}}} \right)-k\left( {{m}_{\text{A}}}, {{m}_{\text{B}}} \right)+k\left( m{{'}_{\text{A}}}, m{{'}_{\text{B}}} \right)$, 所以公式 (5) 的结果为 dχ/dt=0 ^[9]，如图 6所示χ为一常数一直保持初始值。

3 结论与展望

发现χ随着时间的变化满足指数增长函数，该函数取决于稳定状态的混合程度和其初始值之比 (χ_∞/χ₀)、以及初始参数如数目比 (α)、粒径比 (β) 等，根据模拟结果拟合经验模型的过程中，得到相关时间常数 (C₂) 与初始参数的关系。获得自由分子区布朗凝并χ的指数形式预测公式，同时发现在连续区布朗凝并也满足类似的指数关系，并验证模型预测与PBMC模拟结果的误差在0.2以内，该混合程度预测公式对于选择合理的初始参数、优化控制组分分布和实现颗粒定向制备有重要作用。但因为本文只考虑组分无关，即γ=1时颗粒粒径的影响。密度不同的工况 (如γ=ρ_Fe/ρ_Pt) 进行了一定研究，证明稳定状态的混合程度可以由公式 (12) 进行预测，但是准确性有所下降^[14]。我们验算了γ≠1的其他工况，发现只有部分工况满足对于混合程度的稳态值、实时值的具体预测公式，但是混合程度随时间演变过程全部符合公式 (12) 的指数函数形式。因此需要进一步去讨论不同范围的γ值，探究是否存在其他规律以进行预测公式的修正。