有限系統中兩組分顆粒凝并問題的異權值蒙特卡羅模擬

左 昊， 沈 杰， 盧志明

(上海大學上海市應用數學和力學研究所，上海 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.的變化規律值得進一步深入研究。