华中科技大学煤燃烧国家重点实验室, 武汉 430074
2016年04月18日 收稿; 2016年06月30日 收修改稿
基金项目: 国家自然科学基金（51390494，51276077）资助
通信作者: 赵海波, E-mail:klinsmannzhb@163.com
摘要: 多组分颗粒凝并是颗粒长大过程的主要物理机制之一。对于双组分凝并，混合程度χ为一重要衡量标准，是预测组分分布的关键参数。针对双组分颗粒团聚的非稳态过程研究混合程度随时间的变化关系，采用颗粒群平衡模拟异权值快速Monte Carlo方法，进行模拟，最终预测出χ与初始条件的关系，得到自由分子区布朗凝并与连续区布朗凝并的指数形式预测公式，并进行验证，从而可以通过合理选择初始参数，优化控制组分分布和实现颗粒定向功能制备。
关键词: 颗粒群平衡模拟双组分凝并Monte Carlo方法混合程度
Predictions of compositional mixing degree in two-component aggregation
ZHAO Hanqing, XU Zuwei, ZHAO Haibo
State Key Laboratory of Coal Combustion, Huazhong University of Science and Technology, Wuhan 430074, China
Abstract: Multi-component particle aggregation is one of the main physical mechanisms in the process of particle growth. For two-component aggregation, the compositional mixing degree χ is an important criterion and the key to determination of compositional distribution. Now the dynamic evolution of χ before the attainment of a steady-state value is beyond numerical modeling and theoretical research. In this work, the fast differentially-weighted Monte Carlo method for population balance modeling was used to predict the dependence of time-varied χon initial feeding conditions through hundreds of systematically varied simulations. It is found that χis subject to exponential decay in the free molecular regime and the continuum regime. By using the explored exponential formulas for the dynamic mixing degree, it is possible to achieve an optimum control over the compositional distributions during two-component aggregation processes through selecting the initial feeding parameters.
Key words: population balance modelingtwo-component aggregationMonte Carlo methodcompositional mixing degree
在高分子聚集、湿法造粒、纳米颗粒合成等许多过程中，颗粒凝并/团聚是颗粒体积增大的最主要物理机制之一[1-5]，并且通常发生于不同化学组分的颗粒之间，颗粒内部组分分布对于颗粒物理化学性能均有重要甚至决定性的影响[6-8]，因此对于多组分颗粒凝并及组分分布的研究十分重要。对于初始双分散的颗粒群的凝并过程，组分分布可以通过高斯函数来描述：组分A的组分分布以概率密度函数g(mA|m) dmA表示，代表在质量为m的颗粒中组分A含量在质量区间 (mA, mA+dmA) 内的颗粒所占比例。基于随机混合理论和中心限制定理，已证明概率密度函数g(mA|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) |
| $\phi = \frac{{{M_{\text{A}}}}}{{{M_{\text{A}}} + {M_{\text{B}}}}},$ | (2) |
从公式 (1) 可以看到，组分分布函数g(mA|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_{\text{A}}} - \phi m = {m_{\text{A}}} - \phi \left( {{m_{\text{A}}} + {m_{\text{B}}}} \right).$ | (4) |
1 数值方法与计算工况这里采用颗粒群平衡随机模拟的方法 (PBMC)，它能够直接模拟大量颗粒。但因为传统的MC方法对凝并问题的计算代价与模拟颗粒数目Ns的平方成正比，而计算精度与Ns的平方根成正比，异权值MC方法 (DWMC) 能够协调每个尺度区间内模拟颗粒数目，使其保持在一个适当的范围内，在保证足够计算精度的同时使计算更加高效[15-19]。并且为了避免双重遍历，在模拟颗粒凝并时，采用快速DWMC方法加速PBMC模拟[20-21]。具体数值方法详见参考文献[14]。
如表 1所示，我们考虑3种凝并核[21]，且本文仅考虑组分无关工况即双组分密度相同，具体计算工况如表 2所示。3种工况均为经典的模拟工况，线性凝并为一种理想工况，而自由分子区布朗凝并。连续区布朗凝并为实际工况，适用于日常生活中常见的气溶胶系统研究如尘、烟、雾、霾等现象，制药工业、纳米粉体制造业、矿物燃烧过程中超细颗粒物生成。布朗凝并表示凝并过程是由由布朗运动引起的颗粒间的相对运动导致的。我们通过克努德森数来区分 (Kn=2λ/dp，λ为气体分子平均自由程，即相邻两次碰撞之间气体分子移动的平均距离，dp为颗粒直径)。亚微米颗粒 (Kn>10) 处于自由分子区，这时气体分子与颗粒表面的碰撞占据主导地位，气溶胶颗粒的行为接近于自由运动的气体分子，可以使用分子运动论描述气溶胶颗粒间的碰撞凝并现象。微米尺度颗粒或者超微米颗粒 (Kn < 0.1) 处于连续区，此时假设气体为连续介质，与颗粒表面接触处气体速度为0，适用于纳维-斯托克斯方程。
Table 1
 表 1 凝并核与加权强核[21]Table 1 Normal kernels and weighted majorant kernels[21]
| 表 1 凝并核与加权强核[21]Table 1 Normal kernels and weighted majorant kernels[21] |
Table 2
表 2 快速DWMC方法模拟工况Table 2 Conditions used in simulations
| 表 2 快速DWMC方法模拟工况Table 2 Conditions used in simulations |
2 结果与讨论2.1 混合程度的相关参数在整个演变过程中，χ(t)/χ0一直保持相同的数量级，考虑初始工况凝并核函数中的相关参数，我们通过控制变量法已证温度T、初始数目NA0、粒径dA0、和密度ρA对混合程度没有影响[14]。在此基础上，为探讨实时的χ(t)/χ0与初始工况的关系，引入如下参数
| $\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) |
初始混合程度χ0[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) |
| $\chi \left( t \right)/{\chi _0} = {\chi _\infty }/{\chi _0} + {C_1}\exp \left( { - t/{C_2}} \right).$ | (9) |
| $\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) |
Fig. 1
 | Download: JPG larger image |
图 1 时间常数C2与χ∞/χ0的关系 Fig. 1 Time constant C2 vs. χ∞/χ0 | |
首先，根据初始状态t=0，得
| ${C_1} + {\chi _\infty }/{\chi _0} = 1.$ | (11) |
| ${\chi _\infty }/{\chi _0} = \exp \left[ { - 2{{\left( {\phi - {\phi _{\beta = 1}}} \right)}^2}} \right].$ | (12) |
| ${C_2} = {C_3}\ln \left( {{\chi _\infty }/{\chi _0}} \right) + {C_4}.$ | (13) |
Fig. 2
 | Download: JPG larger image |
图 2 C3与χ∞/χ0的关系 Fig. 2 C3 vs. χ∞/χ0 | |
而关于每一条直线截距满足C4=f(α)， 具体计算结果如图 3，C4符合如下关系：
Fig. 3
 | Download: JPG larger image |
图 3 C4与α的关系 Fig. 3 C4 vs. α | |
| ${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) |
Fig. 4
 | Download: JPG larger image |
图 4 MC模拟值与公式预测值 Fig. 4 Comparison between MC simulations and predictions | |
混合程度的实时演变χ(t)/χ0为一单调递减过程，并且下降的速率逐渐变慢，随时间推进越来越接近稳态值。在多种初等函数中，我们发现χ可以满足公式 (9)，χ(t)/χ0以一个与当前值成正比的速度下降
Fig. 5
 | Download: JPG larger image |
图 5 C2—φ公式预测值与模拟值对照 (自由分子区) Fig. 5 Comparison between MC simulations andC2-φ predictions (in free molecule regime) | |
关于该公式的误差，在同时改变φ和φ(分别取{0.1, 0.2, …, 0.90}) 时共9×9个工况，对于指数函数速率时间常数C2公式预测值与MC计算结果的相对误差在0.2以内，如图 5所示。误差产生的原因主要有3点，首先作为经验公式，这种指数形式对于实时χ的预测公式是一种规律的总结，此处会引入一定程度误差，并且常数目的MC方法单次模拟也会与平均值有±5%的波动。其次自由分子区与连续区工况，对于稳定状态的χ的公式预测值与模拟值有6.5%的误差[14]。同时在图 4中，当β=0.441时，C2的公式预测值3.94与MC模拟值4.53之比为1.15，此时图 4混合程度随时间演变的方差为2.96×10-4，因此我们认为图 5的分散度可以接受。
2.3 连续区布朗凝并为了验证χ(t)/χ0公式的适用范围，我们还研究了连续区布朗凝并，由于该核函数考虑颗粒的体积而非质量，所以本质上满足组分无关性。采用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) |
Fig. 6
 | Download: JPG larger image |
图 6 线性凝并工况的χ/χ0vs.t Fig. 6 χ/χ0 vs.t in the linear coagulation case | |
3 结论与展望发现χ随着时间的变化满足指数增长函数，该函数取决于稳定状态的混合程度和其初始值之比 (χ∞/χ0)、以及初始参数如数目比 (α)、粒径比 (β) 等，根据模拟结果拟合经验模型的过程中，得到相关时间常数 (C2) 与初始参数的关系。获得自由分子区布朗凝并χ的指数形式预测公式，同时发现在连续区布朗凝并也满足类似的指数关系，并验证模型预测与PBMC模拟结果的误差在0.2以内，该混合程度预测公式对于选择合理的初始参数、优化控制组分分布和实现颗粒定向制备有重要作用。但因为本文只考虑组分无关，即γ=1时颗粒粒径的影响。密度不同的工况 (如γ=ρFe/ρPt) 进行了一定研究，证明稳定状态的混合程度可以由公式 (12) 进行预测，但是准确性有所下降[14]。我们验算了γ≠1的其他工况，发现只有部分工况满足对于混合程度的稳态值、实时值的具体预测公式，但是混合程度随时间演变过程全部符合公式 (12) 的指数函数形式。因此需要进一步去讨论不同范围的γ值，探究是否存在其他规律以进行预测公式的修正。
参考文献
| [1] | Friedlander S K. Smoke, dust, and haze:fundamentals of aerosol dynamics[M].New York: Oxford University Press, 2000. |
| [2] | Hofmann S, Raisch J. Solutions to inversion problems in preferential crystallization of enantiomers. Part Ⅱ:Batch crystallization in two coupled vessels[J].Chemical Engineering Science, 2013, 88:48–68.DOI:10.1016/j.ces.2012.10.030 |
| [3] | Hosseini A, Bouaswaig A E, Engell S. Novel approaches to improve the particle size distribution prediction of a classical emulsion polymerization model[J].Chemical Engineering Science, 2013, 88:108–120.DOI:10.1016/j.ces.2012.11.021 |
| [4] | Riemer N, West M, Zaveri R, et al. Estimating black carbon aging time-scales with a particle-resolved aerosol model[J].Journal of Aerosol Science, 2010, 41(1):143–158.DOI:10.1016/j.jaerosci.2009.08.009 |
| [5] | Fox R O. Higher-order quadrature-based moment methods for kinetic equations[J].Journal of Computational Physics, 2009, 228(20):7771–7791.DOI:10.1016/j.jcp.2009.07.018 |
| [6] | Barrasso D, Ramachandran R. A comparison of model order reduction techniques for a four-dimensional population balance model describing multi-component wet granulation processes[J].Chemical Engineering Science, 2012, 80:380–392.DOI:10.1016/j.ces.2012.06.039 |
| [7] | Chauhan S S, Chakraborty J, Kumar S. On the solution and applicability of bivariate population balance equations for mixing in particle phase[J].Chemical Engineering Science, 2010, 65(13):3914–3927.DOI:10.1016/j.ces.2010.03.021 |
| [8] | Lushnikov A. Evolution of coagulating systems:Ⅲ. Coagulating mixtures[J].Journal of Colloid and Interface Science, 1976, 54(1):94–101.DOI:10.1016/0021-9797(76)90288-5 |
| [9] | Matsoukas T, Kim T, Lee K. Bicomponent aggregation with composition-dependent rates and the approach to well-mixed state[J].Chemical Engineering Science, 2009, 64(4):787–799.DOI:10.1016/j.ces.2008.04.060 |
| [10] | Matsoukas T, Lee K, Kim T. Mixing of components in two-component aggregation[J].AIChE journal, 2006, 52(9):3088–3099.DOI:10.1002/aic.v52:9 |
| [11] | Lee K, Kim T, Rajniak P, et al. Compositional distributions in multicomponent aggregation[J].Chemical Engineering Science, 2008, 63(5):1293–1303.DOI:10.1016/j.ces.2007.07.060 |
| [12] | Zhao H, Kruis F E, Zheng C. Monte Carlo simulation for aggregative mixing of nanoparticles in two-component systems[J].Industrial & engineering chemistry research, 2011, 50(18):10652–10664. |
| [13] | Marshall C L, Rajniak P, Matsoukas T. Numerical simulations of two-component granulation:comparison of three methods[J].Chemical Engineering Research and Design, 2011, 89(5):545–552.DOI:10.1016/j.cherd.2010.06.003 |
| [14] | Zhao H, Kruis F E. Dependence of steady-state compositional mixing degree on feeding conditions in two-component aggregation[J].Industrial & Engineering Chemistry Research, 2014, 53(14):6047–6055. |
| [15] | Zhao H, Kruis F E, Zheng C. A differentially weighted Monte Carlo method for two-component coagulation[J].Journal of Computational Physics, 2010, 229(19):6931–6945.DOI:10.1016/j.jcp.2010.05.031 |
| [16] | Zhao H, Kruis F E, Zheng C. Reducing statistical noise and extending the size spectrum by applying weighted simulation particles in Monte Carlo simulation of coagulation[J].Aerosol Science and Technology, 2009, 43(8):781–793.DOI:10.1080/02786820902939708 |
| [17] | Zhao H, Zheng C. A new event-driven constant-volume method for solution of the time evolution of particle size distribution[J].Journal of Computational Physics, 2009, 228(5):1412–1428.DOI:10.1016/j.jcp.2008.10.033 |
| [18] | Zhao H, Zheng C. Two-component Brownian coagulation:Monte Carlo simulation and process characterization[J].Particuology, 2011, 9(4):414–423.DOI:10.1016/j.partic.2011.04.003 |
| [19] | Zhao H, Zheng C. A population balance-Monte Carlo method for particle coagulation in spatially inhomogeneous systems[J].Computers & Fluids, 2013, 71:196–207. |
| [20] | Hao X, Zhao H, Xu Z, et al. Population balance-Monte Carlo simulation for gas-to-particle synthesis of nanoparticles[J].Aerosol science and technology, 2013, 47(10):1125–1133.DOI:10.1080/02786826.2013.823642 |
| [21] | Xu Z, Zhao H, Zheng C. Fast Monte Carlo simulation for particle coagulation in population balance[J].Journal of Aerosol Science, 2014, 74:11–25.DOI:10.1016/j.jaerosci.2014.03.006 |
