左 昊， 沈 杰， 卢志明

(上海大学上海市应用数学和力学研究所，上海 200072)

自然界与工程领域中都广泛存在着气溶胶系统[1-3]。研究气溶胶系统，一般采用颗粒群平衡模拟(population balance modeling, PBM)的方法来描述该系统中离散相颗粒群的动力学演变过程，最早是通过构建一个以颗粒尺度分布函数为基础变量的单变量颗粒群平衡方程(popolation balance equation, PBE)来描述颗粒系统的状态。

然而，仅以颗粒尺度为单变量的PBE 方程描述颗粒分布状态和凝并过程并不能考虑凝并过程中的许多重要因素，如药物凝并混合过程中组分之间的相互影响，颗粒不规则度、表面积、表面活化能、核电特性等对颗粒群演化过程的影响。因此，多组分多变量颗粒系统的动力学演变过程的研究有着重要的科学意义和广泛的工业用途[2-9]。Lushnikov[10]最早进行了两组分离散颗粒系统的研究，并得到了常数核模型的解析解。之后，Gelbard等[11]研究了多组分连续系统，并得到了颗粒平衡方程的解析解。Vigil等[12]发展了一个复杂加法核模型，并得到了该模型的渐进解。近年来，Fern´andez-D´ıa等[13-14]进行了两组分颗粒系统的研究，得到了加法核以及乘法核模型PBE 的解析解。Psikunov[15]研究了颗粒系统中同时存在凝并和冷凝两种动力学行为时颗粒凝并的情况。Yu等[16]则将Vigil 的复杂加法核模型推广到了连续颗粒系统。同时，Barrasso等[17]研究了多组分系统的连续凝并过程，并与相关实验进行了对比。Hashemian等[18]利用Laguerre多项式发展了一个两组分颗粒系统的降阶模型。而Burgalat等[19]发展了球体和不规则体两种不同的颗粒凝并模型。以上理论研究或建模都有很好的成效，但对于两组分颗粒系统的研究而言，数值模拟方法的发展也不可忽视。

Matsoukas等[20]和Lee等[21]发展了常数目蒙特卡罗(constant-number Monte Carlo)方法，并提出了一个与组分相关的核模型，即在核函数的基础上加上一个“凝并效率”函数，以此研究凝并过程依赖组分时的情况，并提出用总偏差和分散指数研究双组分颗粒系统的混合凝并状态。而Zhao等[22-23]发展了异权值蒙特卡罗(multi-Monte Carlo, MMC)方法，该方法摒弃了之前的蒙特卡罗方法中的“子系统”的概念，引入了“虚拟颗粒”的概念，通过附加的权值的变化保证模拟过程中模拟颗粒数目和计算域体积都保持不变，从而在计算精度和计算效率上大大提高[24-25]。此后，Zhao等[26]又进一步研究发现了系统凝并混合达到稳定状态前的演化规律。针对以前研究都假定系统颗粒数目无限，系统能够无限演化的问题，Matsoukas等[27-28]提出了有限颗粒系统的概念，并研究了有限颗粒系统中颗粒凝并的演化规律，指出影响颗粒系统凝并和混合状态的是核函数的齐次性指数。虽然Matsoukas 等提出了有限系统的概念，并指出核函数的齐次性指数的重要性，但对于有限系统以及核函数齐次性的研究还不全面。因此，本工作采用异权值蒙特卡罗方法对有限的两组分系统开展研究，同时考虑组分间相互作用，对两组分颗粒系统的演化规律进行深入研究。

1 理论和方法

本工作的研究对象是一个颗粒数目为N, 总质量为M的两组分颗粒系统。为了简便，用A组分和B组分表示颗粒系统中的两个组分。两变量的PBE 方程可以由Smoluchowski 方程推广得到，

式中:r= (v,m)表示颗粒群内部变量集，其中m代表颗粒中组分A的质量(或体积)，v代表颗粒的总质量(或总体积)；F(r,t)表示颗粒群内部变量集r的分布函数，F(r,t)dr表示t时刻、变量范围在r和r+dr之间的颗粒在单位体积内的数目浓度；K(r,r′)表示核函数，即颗粒系统凝并过程中动力学事件发生服从的某种概率。因此，组分A在颗粒系统中的总质量分数可以定义为

组分A在单个颗粒内的质量分数c=m/v，在系统中的平均质量分数为φ。当颗粒系统达到完全混合均匀状态时，c=φ；而当系统还没有到达完全混合均匀状态时，二者会存在偏差。定义组分A的偏差x=v(c-φ)，则其在系统中的总偏差为

因为颗粒的平均质量v=M/N，所以

式中：〈·〉表示对所有颗粒取平均。

研究颗粒凝并过程中系统χ的变化，有利于了解系统的混合均匀程度。Matsoukas等[27]提出了χ和系统平均体积v的函数关系：

对于加法性质的核函数或者初始部分混合均匀的组分独立核函数，χ和v的关系式可以进一步简化为

对于一个近似有限系统(N有限)而言，Matsoukas等[27]给出了常数核的精确解：

显然，当N趋于无穷大时，χ是一个常数。对于一个齐次性指数为λ的凝并核，Matsoukas等[28]推导得到了

式(8)表明，齐次性指数λ是影响颗粒凝并过程的重要物理量。为了进一步研究颗粒系统的不均匀性，通常引入一个无量纲的参数，即分散指数(segregation index, S.I.)，

显然，对于常数核，根据式(7)可以得到

对于两组分颗粒系统凝并的混合程度问题，核函数中的齐次指数对颗粒凝并起着关键作用。因此，本工作将采用异权值蒙特卡罗[22-23]方法深入研究齐次指数的影响。异权值蒙特卡罗方法脱胎于蒙特卡罗方法，发展了一个异数目权值虚拟颗粒的策略，采用数目权值不等的虚拟颗粒群来代表实际颗粒群，尽可能多地继承了颗粒群的尺度分布信息，具有高精度、高效率等优点[24-25]。

2 模拟结果和讨论

2.1 凝并核模型

为了深入研究凝并核齐次指数对凝并的影响，本工作针对4 种常见的核函数，构造了一个简单的核函数与之对比。

(1) 常数核(无量纲形式)：

(2) 加法核(无量纲形式)：

(3) 颗粒直径在纳米级的自由分子区布朗核(无量纲形式)：

(4) 颗粒直径在微米级的连续区布朗核(无量纲形式)：

(5) 构造的核函数(无量纲形式)：

各个凝并核的齐次性指数见表1。

表1 不同核函数的齐次性Table 1 Degree of homogeneity for different kernels

另外，本工作还考虑了组分间吸引排斥的颗粒凝并问题（见表2）。采用Matsoukas等[27]首次提出的模型：

式中：ci表示颗粒i中组分A所占的质量分数，ψ(c1,c2)可以看成在核函数的基础上添加一个“凝并效率”函数；α是吸引排斥指数。当α为正时，两颗粒同时含有组分A或都没有组分A的几率较高，一个颗粒全部是组分A而一个颗粒没有组分A的几率较低，不利于混合，称这种情况为“排斥”，α值越大，“排斥”作用越强；当α为负时，结论相反，两个颗粒中一颗全部是组分A，另一个没有组分A的几率较高，利于混合，称这种情况为“吸引”，α值越小，“吸引”作用越强。

Jiang等[29]研究发现两组分系统颗粒凝并结果受初始时A组分体积占比影响很小。因此，本工作选取模拟的初值条件是颗粒单分散分布, 即所有颗粒初始尺度大小v相同,A组分颗粒和B组分颗粒各占一半。

表2 初始工况Table 2 Initial distribution

2.2 模拟验证

图1 给出了常数核情况下本工作的模拟结果与常数目法以及解析式(7)的对比情况，其中横坐标是平均颗粒体积与初始平均颗粒体积的比值，以此作为颗粒系统凝并的时间尺度；纵坐标是系统的χ值。为说明本工作模拟结果的合理性(N代表异权值蒙特卡罗方法中的虚拟颗粒数目，虽然N值有限，但系统会无限凝并)，图(1)仅给出了N= 10 和100 两种情况。由图1 可以发现：对比常数目蒙特卡罗方法的模拟结果，本工作的模拟结果精度更好，与理论解吻合得更好，但仍有细小误差；随着N的增大，如解析式(7)所示，χ值下降越来越慢；当N值无限大时，χ应该无限趋近于水平线，这意味着系统的混合状态不发生改变，即是一个无限系统。

图1 常数核模拟结果Fig.1 Simulation results for the constant kernel

2.3 有限系统中不同核函数的颗粒凝并

当N值较大时，核函数对系统凝并的影响被掩盖。因此，本工作在N=100 的系统(近似为一个有限系统)中研究核函数齐次性指数λ对系统凝并的影响, 即双组分系统中χ和S.I. 的影响。图2 给出了N=100 不同核函数的模拟结果。由图2 可以发现：常数核和布朗核连续区的χ值演化结果基本一致；布朗核自由分子区χ值下降略快于前者，但与λ=1/6 的构造核演化结果接近一致；λ= 1/2 构造核与加法核的χ值下降更快。模拟结果表明，凝并核函数齐次性指数λ值决定了χ值的演化，且λ值越大，χ值下降越快，系统混合程度越好。

图2 不同核函数的颗粒凝并(N =100)Fig.2 Simulation results for different kernels (N =100)

为了更好地揭示核函数的齐次性对颗粒凝并、组分混合过程的影响，考察颗粒系统分散指数的变化。图3 给出了核函数S.I.的演化规律。可以看出：随着齐次性指数的增大，S.I.下降更快；当齐次性指数相同时，S.I.曲线非常接近；S.I.呈直线下降，其斜率与齐次性指数有关，齐次性越大，斜率越大。

图3 不同核函数的分散指数Fig.3 Segregation index of different kernels

图4 k 与λ 的关系Fig.4 Relationship between k and λ

从图4 可以看出：当λ≤0.5 时，模拟结果较好；当λ≥0.5 时，结果偏差较大。这是由于S.I.的演化速度加剧，在双对数坐标下不完全呈线性下降，此阶段需要单独研究，而不适合λ≤0.5 时的演化规律。由式(17)可得，经过长时间演化后，不同核函数的S.I.下降满足

显然，与常数核(λ=0)相比, 其他的核函数多出了一项

2.4 有限系统中考虑吸引或排斥作用的颗粒凝并

在有限系统中，考虑两组分之间存在吸引或排斥作用时，核函数采用式(16)，选取不同的α值表示颗粒组分对凝并效率的影响。

以常数核(式(16)中k取常数)为例，图5 给出了有限系统中(N= 10)χ值演化的模拟结果。可以发现：当α值为正，即两组分间存在排斥作用时，描述系统混合状态的参数χ值变化可以分作两个阶段。第一个阶段，χ值随时间迅速增加，再缓慢减小，这个阶段认为核函数和组分间的排斥作用在共同影响着χ的演化；第二个阶段，χ是线性减小的，且无论α值是多少，χ值演化曲线的标度率和α= 0 时基本一致。因此认为，第二阶段主要是核函数在影响着系统χ值的演化，而组分间的排斥作用对χ值的影响基本可以忽略不计。由于当N值较大时，颗粒系统可以认为是一个无限系统。Matsoukas等[27]发现，在无限系统中考虑组分间吸引排斥作用时的颗粒凝并情况，χ会在足够演化时间后保持不变，系统达到“自保持分布状态”。

当α值为负时，即两组分间存在吸引作用时，描述系统混合状态的参数χ值变化同样可以分为两个阶段：第一个阶段，χ值随时间迅速减小，这个阶段认为是核函数以及组分间的吸引作用在共同影响着χ值的演化；第二个阶段，χ值线性减小，且无论α值是多少，χ值演化曲线的标度率和α= 0 时基本一致。所以认为，第二阶段主要是核函数在影响着系统χ值的演化，而组分间的吸引作用对χ值的影响基本可以忽略不计。同理，如果N值较大，系统会达到“自保持分布状态”。

图5 常数核χ 值演化结果(N =10)Fig.5 Evolution results for constant kernel (N =10)

其他形式的核函数（如布朗核、构造核和加法核等）的模拟结果与图5 结果相似，即χ值都可以分成两个阶段。系统何时达到均匀混合状态，对药物混合过程是一个重要的物理量。因此，系统达到第二个阶段的临界时间对于研究颗粒系统有着重要意义。这里的临界时间指系统凝并仅由核函数主导而与吸引排斥作用无关的时间。模拟结果表明临界时间与模拟颗粒数目N基本无关，但与吸引排斥指数α以及齐次性λ密切相关。

图6 给出了N= 10 时不同核函数临界时间随α的变化情况。由图6 可以看出：临界时间随着α的增大而增加；临界时间与α在双对数坐标下呈直线，且斜率与λ无关。这说明：当α为正时，临界时间随α呈幂函数增长；当α为负时，临界时间随-α呈幂函数增长，但增长幅度不如α为正时明显。拟合得到的标度指数分别为2.4 和0.8，即有

图6 6 种核函数在不同的α 下的临界时间Fig.6 Critical time with different α for 6 kernels

从图6 还可以看出，当α固定时，临界时间随着λ的增加而减小。图7 给出了临界时间与λ的定量关系，其中纵坐标为确定的一个α值下不同λ的临界时间除以λ= 0 的临界时间。从图7 可以看出：对不同的α值，曲线几乎完全重合；无量纲临界时间与λ呈指数函数减小,且当α为正时λ对临界时间的影响比较明显,而当α为负时λ对临界时间的影响稍小。通过拟合可得

图7 临界时间与λ 的关系Fig.7 Relationship between the critical time and λ

3 结束语

本工作采用异权值蒙特卡罗方法模拟了数目有限的两组分颗粒系统中的凝并混合问题。首先，进一步验证了凝并核函数齐次性指数λ对系统混合凝并的重要影响，并得到了有限系统中分散指数S.I. 的演化与λ的关系式。其次，发现并总结了有限系统中考虑吸引或排斥作用时颗粒凝并混合情况，即χ值的演化规律：临界时间之前系统在吸引排斥因素以及核函数的共同作用下进行颗粒凝并，且排斥作用的影响远远大于吸引作用的影响；临界时间之后系统在核函数的单独作用下进行颗粒凝并，排斥或吸引作用的影响可以忽略不计。最后，通过数据拟合得到临界时间与α的幂函数关系式，以及与齐次性λ的指数函数关系。本工作定量分析了齐次性指数对双组分颗粒凝并和组分混合过程的影响，所得结果对药物混合工业过程的设计具有一定的指导意义。

本工作中采用异权值蒙特卡罗方法进行颗粒系统凝并研究，由于其波动性较大，结果均为10～10 000 次运算的平均结果。当N值较大(近似无限系统)时，不同核函数对系统混合凝并的影响难以体现，所以在N= 10(或100)的系统(近似有限系统)中研究颗粒的凝并情况。当λ≥0.5 时，S.I.的变化规律值得进一步深入研究。