有限系統(tǒng)中兩組分顆粒凝并問題的異權(quán)值蒙特卡羅模擬

左 昊， 沈 杰， 盧志明

(上海大學(xué)上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所，上海 200072)

自然界與工程領(lǐng)域中都廣泛存在著氣溶膠系統(tǒng)[1-3]。研究氣溶膠系統(tǒng)，一般采用顆粒群平衡模擬(population balance modeling, PBM)的方法來(lái)描述該系統(tǒng)中離散相顆粒群的動(dòng)力學(xué)演變過(guò)程，最早是通過(guò)構(gòu)建一個(gè)以顆粒尺度分布函數(shù)為基礎(chǔ)變量的單變量顆粒群平衡方程(popolation balance equation, PBE)來(lái)描述顆粒系統(tǒng)的狀態(tài)。

然而，僅以顆粒尺度為單變量的PBE 方程描述顆粒分布狀態(tài)和凝并過(guò)程并不能考慮凝并過(guò)程中的許多重要因素，如藥物凝并混合過(guò)程中組分之間的相互影響，顆粒不規(guī)則度、表面積、表面活化能、核電特性等對(duì)顆粒群演化過(guò)程的影響。因此，多組分多變量顆粒系統(tǒng)的動(dòng)力學(xué)演變過(guò)程的研究有著重要的科學(xué)意義和廣泛的工業(yè)用途[2-9]。Lushnikov[10]最早進(jìn)行了兩組分離散顆粒系統(tǒng)的研究，并得到了常數(shù)核模型的解析解。之后，Gelbard等[11]研究了多組分連續(xù)系統(tǒng)，并得到了顆粒平衡方程的解析解。Vigil等[12]發(fā)展了一個(gè)復(fù)雜加法核模型，并得到了該模型的漸進(jìn)解。近年來(lái)，F(xiàn)ern′andez-D′?a等[13-14]進(jìn)行了兩組分顆粒系統(tǒng)的研究，得到了加法核以及乘法核模型PBE 的解析解。Psikunov[15]研究了顆粒系統(tǒng)中同時(shí)存在凝并和冷凝兩種動(dòng)力學(xué)行為時(shí)顆粒凝并的情況。Yu等[16]則將Vigil 的復(fù)雜加法核模型推廣到了連續(xù)顆粒系統(tǒng)。同時(shí)，Barrasso等[17]研究了多組分系統(tǒng)的連續(xù)凝并過(guò)程，并與相關(guān)實(shí)驗(yàn)進(jìn)行了對(duì)比。Hashemian等[18]利用Laguerre多項(xiàng)式發(fā)展了一個(gè)兩組分顆粒系統(tǒng)的降階模型。而Burgalat等[19]發(fā)展了球體和不規(guī)則體兩種不同的顆粒凝并模型。以上理論研究或建模都有很好的成效，但對(duì)于兩組分顆粒系統(tǒng)的研究而言，數(shù)值模擬方法的發(fā)展也不可忽視。

Matsoukas等[20]和Lee等[21]發(fā)展了常數(shù)目蒙特卡羅(constant-number Monte Carlo)方法，并提出了一個(gè)與組分相關(guān)的核模型，即在核函數(shù)的基礎(chǔ)上加上一個(gè)“凝并效率”函數(shù)，以此研究凝并過(guò)程依賴組分時(shí)的情況，并提出用總偏差和分散指數(shù)研究雙組分顆粒系統(tǒng)的混合凝并狀態(tài)。而Zhao等[22-23]發(fā)展了異權(quán)值蒙特卡羅(multi-Monte Carlo, MMC)方法，該方法摒棄了之前的蒙特卡羅方法中的“子系統(tǒng)”的概念，引入了“虛擬顆?！钡母拍?，通過(guò)附加的權(quán)值的變化保證模擬過(guò)程中模擬顆粒數(shù)目和計(jì)算域體積都保持不變，從而在計(jì)算精度和計(jì)算效率上大大提高[24-25]。此后，Zhao等[26]又進(jìn)一步研究發(fā)現(xiàn)了系統(tǒng)凝并混合達(dá)到穩(wěn)定狀態(tài)前的演化規(guī)律。針對(duì)以前研究都假定系統(tǒng)顆粒數(shù)目無(wú)限，系統(tǒng)能夠無(wú)限演化的問題，Matsoukas等[27-28]提出了有限顆粒系統(tǒng)的概念，并研究了有限顆粒系統(tǒng)中顆粒凝并的演化規(guī)律，指出影響顆粒系統(tǒng)凝并和混合狀態(tài)的是核函數(shù)的齊次性指數(shù)。雖然Matsoukas 等提出了有限系統(tǒng)的概念，并指出核函數(shù)的齊次性指數(shù)的重要性，但對(duì)于有限系統(tǒng)以及核函數(shù)齊次性的研究還不全面。因此，本工作采用異權(quán)值蒙特卡羅方法對(duì)有限的兩組分系統(tǒng)開展研究，同時(shí)考慮組分間相互作用，對(duì)兩組分顆粒系統(tǒng)的演化規(guī)律進(jìn)行深入研究。

1 理論和方法

本工作的研究對(duì)象是一個(gè)顆粒數(shù)目為N, 總質(zhì)量為M的兩組分顆粒系統(tǒng)。為了簡(jiǎn)便，用A組分和B組分表示顆粒系統(tǒng)中的兩個(gè)組分。兩變量的PBE 方程可以由Smoluchowski 方程推廣得到，

式中:r= (v,m)表示顆粒群內(nèi)部變量集，其中m代表顆粒中組分A的質(zhì)量(或體積)，v代表顆粒的總質(zhì)量(或總體積)；F(r,t)表示顆粒群內(nèi)部變量集r的分布函數(shù)，F(xiàn)(r,t)dr表示t時(shí)刻、變量范圍在r和r+dr之間的顆粒在單位體積內(nèi)的數(shù)目濃度；K(r,r′)表示核函數(shù)，即顆粒系統(tǒng)凝并過(guò)程中動(dòng)力學(xué)事件發(fā)生服從的某種概率。因此，組分A在顆粒系統(tǒng)中的總質(zhì)量分?jǐn)?shù)可以定義為

組分A在單個(gè)顆粒內(nèi)的質(zhì)量分?jǐn)?shù)c=m/v，在系統(tǒng)中的平均質(zhì)量分?jǐn)?shù)為φ。當(dāng)顆粒系統(tǒng)達(dá)到完全混合均勻狀態(tài)時(shí)，c=φ；而當(dāng)系統(tǒng)還沒有到達(dá)完全混合均勻狀態(tài)時(shí)，二者會(huì)存在偏差。定義組分A的偏差x=v(c-φ)，則其在系統(tǒng)中的總偏差為

因?yàn)轭w粒的平均質(zhì)量v=M/N，所以

式中：〈·〉表示對(duì)所有顆粒取平均。

研究顆粒凝并過(guò)程中系統(tǒng)χ的變化，有利于了解系統(tǒng)的混合均勻程度。Matsoukas等[27]提出了χ和系統(tǒng)平均體積v的函數(shù)關(guān)系：

對(duì)于加法性質(zhì)的核函數(shù)或者初始部分混合均勻的組分獨(dú)立核函數(shù)，χ和v的關(guān)系式可以進(jìn)一步簡(jiǎn)化為

對(duì)于一個(gè)近似有限系統(tǒng)(N有限)而言，Matsoukas等[27]給出了常數(shù)核的精確解：

顯然，當(dāng)N趨于無(wú)窮大時(shí)，χ是一個(gè)常數(shù)。對(duì)于一個(gè)齊次性指數(shù)為λ的凝并核，Matsoukas等[28]推導(dǎo)得到了

式(8)表明，齊次性指數(shù)λ是影響顆粒凝并過(guò)程的重要物理量。為了進(jìn)一步研究顆粒系統(tǒng)的不均勻性，通常引入一個(gè)無(wú)量綱的參數(shù)，即分散指數(shù)(segregation index, S.I.)，

顯然，對(duì)于常數(shù)核，根據(jù)式(7)可以得到

對(duì)于兩組分顆粒系統(tǒng)凝并的混合程度問題，核函數(shù)中的齊次指數(shù)對(duì)顆粒凝并起著關(guān)鍵作用。因此，本工作將采用異權(quán)值蒙特卡羅[22-23]方法深入研究齊次指數(shù)的影響。異權(quán)值蒙特卡羅方法脫胎于蒙特卡羅方法，發(fā)展了一個(gè)異數(shù)目權(quán)值虛擬顆粒的策略，采用數(shù)目權(quán)值不等的虛擬顆粒群來(lái)代表實(shí)際顆粒群，盡可能多地繼承了顆粒群的尺度分布信息，具有高精度、高效率等優(yōu)點(diǎn)[24-25]。

2 模擬結(jié)果和討論

2.1 凝并核模型

為了深入研究凝并核齊次指數(shù)對(duì)凝并的影響，本工作針對(duì)4 種常見的核函數(shù)，構(gòu)造了一個(gè)簡(jiǎn)單的核函數(shù)與之對(duì)比。

(1) 常數(shù)核(無(wú)量綱形式)：

(2) 加法核(無(wú)量綱形式)：

(3) 顆粒直徑在納米級(jí)的自由分子區(qū)布朗核(無(wú)量綱形式)：

(4) 顆粒直徑在微米級(jí)的連續(xù)區(qū)布朗核(無(wú)量綱形式)：

(5) 構(gòu)造的核函數(shù)(無(wú)量綱形式)：

各個(gè)凝并核的齊次性指數(shù)見表1。

表1 不同核函數(shù)的齊次性Table 1 Degree of homogeneity for different kernels

另外，本工作還考慮了組分間吸引排斥的顆粒凝并問題（見表2）。采用Matsoukas等[27]首次提出的模型：

式中：ci表示顆粒i中組分A所占的質(zhì)量分?jǐn)?shù)，ψ(c1,c2)可以看成在核函數(shù)的基礎(chǔ)上添加一個(gè)“凝并效率”函數(shù)；α是吸引排斥指數(shù)。當(dāng)α為正時(shí)，兩顆粒同時(shí)含有組分A或都沒有組分A的幾率較高，一個(gè)顆粒全部是組分A而一個(gè)顆粒沒有組分A的幾率較低，不利于混合，稱這種情況為“排斥”，α值越大，“排斥”作用越強(qiáng)；當(dāng)α為負(fù)時(shí)，結(jié)論相反，兩個(gè)顆粒中一顆全部是組分A，另一個(gè)沒有組分A的幾率較高，利于混合，稱這種情況為“吸引”，α值越小，“吸引”作用越強(qiáng)。

Jiang等[29]研究發(fā)現(xiàn)兩組分系統(tǒng)顆粒凝并結(jié)果受初始時(shí)A組分體積占比影響很小。因此，本工作選取模擬的初值條件是顆粒單分散分布, 即所有顆粒初始尺度大小v相同,A組分顆粒和B組分顆粒各占一半。

表2 初始工況Table 2 Initial distribution

2.2 模擬驗(yàn)證

圖1 給出了常數(shù)核情況下本工作的模擬結(jié)果與常數(shù)目法以及解析式(7)的對(duì)比情況，其中橫坐標(biāo)是平均顆粒體積與初始平均顆粒體積的比值，以此作為顆粒系統(tǒng)凝并的時(shí)間尺度；縱坐標(biāo)是系統(tǒng)的χ值。為說(shuō)明本工作模擬結(jié)果的合理性(N代表異權(quán)值蒙特卡羅方法中的虛擬顆粒數(shù)目，雖然N值有限，但系統(tǒng)會(huì)無(wú)限凝并)，圖(1)僅給出了N= 10 和100 兩種情況。由圖1 可以發(fā)現(xiàn)：對(duì)比常數(shù)目蒙特卡羅方法的模擬結(jié)果，本工作的模擬結(jié)果精度更好，與理論解吻合得更好，但仍有細(xì)小誤差；隨著N的增大，如解析式(7)所示，χ值下降越來(lái)越慢；當(dāng)N值無(wú)限大時(shí)，χ應(yīng)該無(wú)限趨近于水平線，這意味著系統(tǒng)的混合狀態(tài)不發(fā)生改變，即是一個(gè)無(wú)限系統(tǒng)。

圖1 常數(shù)核模擬結(jié)果Fig.1 Simulation results for the constant kernel

2.3 有限系統(tǒng)中不同核函數(shù)的顆粒凝并

當(dāng)N值較大時(shí)，核函數(shù)對(duì)系統(tǒng)凝并的影響被掩蓋。因此，本工作在N=100 的系統(tǒng)(近似為一個(gè)有限系統(tǒng))中研究核函數(shù)齊次性指數(shù)λ對(duì)系統(tǒng)凝并的影響, 即雙組分系統(tǒng)中χ和S.I. 的影響。圖2 給出了N=100 不同核函數(shù)的模擬結(jié)果。由圖2 可以發(fā)現(xiàn)：常數(shù)核和布朗核連續(xù)區(qū)的χ值演化結(jié)果基本一致；布朗核自由分子區(qū)χ值下降略快于前者，但與λ=1/6 的構(gòu)造核演化結(jié)果接近一致；λ= 1/2 構(gòu)造核與加法核的χ值下降更快。模擬結(jié)果表明，凝并核函數(shù)齊次性指數(shù)λ值決定了χ值的演化，且λ值越大，χ值下降越快，系統(tǒng)混合程度越好。

圖2 不同核函數(shù)的顆粒凝并(N =100)Fig.2 Simulation results for different kernels (N =100)

為了更好地揭示核函數(shù)的齊次性對(duì)顆粒凝并、組分混合過(guò)程的影響，考察顆粒系統(tǒng)分散指數(shù)的變化。圖3 給出了核函數(shù)S.I.的演化規(guī)律?？梢钥闯觯弘S著齊次性指數(shù)的增大，S.I.下降更快；當(dāng)齊次性指數(shù)相同時(shí)，S.I.曲線非常接近；S.I.呈直線下降，其斜率與齊次性指數(shù)有關(guān)，齊次性越大，斜率越大。

圖3 不同核函數(shù)的分散指數(shù)Fig.3 Segregation index of different kernels

圖4 k 與λ 的關(guān)系Fig.4 Relationship between k and λ

從圖4 可以看出：當(dāng)λ≤0.5 時(shí)，模擬結(jié)果較好；當(dāng)λ≥0.5 時(shí)，結(jié)果偏差較大。這是由于S.I.的演化速度加劇，在雙對(duì)數(shù)坐標(biāo)下不完全呈線性下降，此階段需要單獨(dú)研究，而不適合λ≤0.5 時(shí)的演化規(guī)律。由式(17)可得，經(jīng)過(guò)長(zhǎng)時(shí)間演化后，不同核函數(shù)的S.I.下降滿足

顯然，與常數(shù)核(λ=0)相比, 其他的核函數(shù)多出了一項(xiàng)

2.4 有限系統(tǒng)中考慮吸引或排斥作用的顆粒凝并

在有限系統(tǒng)中，考慮兩組分之間存在吸引或排斥作用時(shí)，核函數(shù)采用式(16)，選取不同的α值表示顆粒組分對(duì)凝并效率的影響。

以常數(shù)核(式(16)中k取常數(shù))為例，圖5 給出了有限系統(tǒng)中(N= 10)χ值演化的模擬結(jié)果?？梢园l(fā)現(xiàn)：當(dāng)α值為正，即兩組分間存在排斥作用時(shí)，描述系統(tǒng)混合狀態(tài)的參數(shù)χ值變化可以分作兩個(gè)階段。第一個(gè)階段，χ值隨時(shí)間迅速增加，再緩慢減小，這個(gè)階段認(rèn)為核函數(shù)和組分間的排斥作用在共同影響著χ的演化；第二個(gè)階段，χ是線性減小的，且無(wú)論α值是多少，χ值演化曲線的標(biāo)度率和α= 0 時(shí)基本一致。因此認(rèn)為，第二階段主要是核函數(shù)在影響著系統(tǒng)χ值的演化，而組分間的排斥作用對(duì)χ值的影響基本可以忽略不計(jì)。由于當(dāng)N值較大時(shí)，顆粒系統(tǒng)可以認(rèn)為是一個(gè)無(wú)限系統(tǒng)。Matsoukas等[27]發(fā)現(xiàn)，在無(wú)限系統(tǒng)中考慮組分間吸引排斥作用時(shí)的顆粒凝并情況，χ會(huì)在足夠演化時(shí)間后保持不變，系統(tǒng)達(dá)到“自保持分布狀態(tài)”。

當(dāng)α值為負(fù)時(shí)，即兩組分間存在吸引作用時(shí)，描述系統(tǒng)混合狀態(tài)的參數(shù)χ值變化同樣可以分為兩個(gè)階段：第一個(gè)階段，χ值隨時(shí)間迅速減小，這個(gè)階段認(rèn)為是核函數(shù)以及組分間的吸引作用在共同影響著χ值的演化；第二個(gè)階段，χ值線性減小，且無(wú)論α值是多少，χ值演化曲線的標(biāo)度率和α= 0 時(shí)基本一致。所以認(rèn)為，第二階段主要是核函數(shù)在影響著系統(tǒng)χ值的演化，而組分間的吸引作用對(duì)χ值的影響基本可以忽略不計(jì)。同理，如果N值較大，系統(tǒng)會(huì)達(dá)到“自保持分布狀態(tài)”。

圖5 常數(shù)核χ 值演化結(jié)果(N =10)Fig.5 Evolution results for constant kernel (N =10)

其他形式的核函數(shù)（如布朗核、構(gòu)造核和加法核等）的模擬結(jié)果與圖5 結(jié)果相似，即χ值都可以分成兩個(gè)階段。系統(tǒng)何時(shí)達(dá)到均勻混合狀態(tài)，對(duì)藥物混合過(guò)程是一個(gè)重要的物理量。因此，系統(tǒng)達(dá)到第二個(gè)階段的臨界時(shí)間對(duì)于研究顆粒系統(tǒng)有著重要意義。這里的臨界時(shí)間指系統(tǒng)凝并僅由核函數(shù)主導(dǎo)而與吸引排斥作用無(wú)關(guān)的時(shí)間。模擬結(jié)果表明臨界時(shí)間與模擬顆粒數(shù)目N基本無(wú)關(guān)，但與吸引排斥指數(shù)α以及齊次性λ密切相關(guān)。

圖6 給出了N= 10 時(shí)不同核函數(shù)臨界時(shí)間隨α的變化情況。由圖6 可以看出：臨界時(shí)間隨著α的增大而增加；臨界時(shí)間與α在雙對(duì)數(shù)坐標(biāo)下呈直線，且斜率與λ無(wú)關(guān)。這說(shuō)明：當(dāng)α為正時(shí)，臨界時(shí)間隨α呈冪函數(shù)增長(zhǎng)；當(dāng)α為負(fù)時(shí)，臨界時(shí)間隨-α呈冪函數(shù)增長(zhǎng)，但增長(zhǎng)幅度不如α為正時(shí)明顯。擬合得到的標(biāo)度指數(shù)分別為2.4 和0.8，即有

圖6 6 種核函數(shù)在不同的α 下的臨界時(shí)間Fig.6 Critical time with different α for 6 kernels

從圖6 還可以看出，當(dāng)α固定時(shí)，臨界時(shí)間隨著λ的增加而減小。圖7 給出了臨界時(shí)間與λ的定量關(guān)系，其中縱坐標(biāo)為確定的一個(gè)α值下不同λ的臨界時(shí)間除以λ= 0 的臨界時(shí)間。從圖7 可以看出：對(duì)不同的α值，曲線幾乎完全重合；無(wú)量綱臨界時(shí)間與λ呈指數(shù)函數(shù)減小,且當(dāng)α為正時(shí)λ對(duì)臨界時(shí)間的影響比較明顯,而當(dāng)α為負(fù)時(shí)λ對(duì)臨界時(shí)間的影響稍小。通過(guò)擬合可得

圖7 臨界時(shí)間與λ 的關(guān)系Fig.7 Relationship between the critical time and λ

3 結(jié)束語(yǔ)

本工作采用異權(quán)值蒙特卡羅方法模擬了數(shù)目有限的兩組分顆粒系統(tǒng)中的凝并混合問題。首先，進(jìn)一步驗(yàn)證了凝并核函數(shù)齊次性指數(shù)λ對(duì)系統(tǒng)混合凝并的重要影響，并得到了有限系統(tǒng)中分散指數(shù)S.I. 的演化與λ的關(guān)系式。其次，發(fā)現(xiàn)并總結(jié)了有限系統(tǒng)中考慮吸引或排斥作用時(shí)顆粒凝并混合情況，即χ值的演化規(guī)律：臨界時(shí)間之前系統(tǒng)在吸引排斥因素以及核函數(shù)的共同作用下進(jìn)行顆粒凝并，且排斥作用的影響遠(yuǎn)遠(yuǎn)大于吸引作用的影響；臨界時(shí)間之后系統(tǒng)在核函數(shù)的單獨(dú)作用下進(jìn)行顆粒凝并，排斥或吸引作用的影響可以忽略不計(jì)。最后，通過(guò)數(shù)據(jù)擬合得到臨界時(shí)間與α的冪函數(shù)關(guān)系式，以及與齊次性λ的指數(shù)函數(shù)關(guān)系。本工作定量分析了齊次性指數(shù)對(duì)雙組分顆粒凝并和組分混合過(guò)程的影響，所得結(jié)果對(duì)藥物混合工業(yè)過(guò)程的設(shè)計(jì)具有一定的指導(dǎo)意義。

本工作中采用異權(quán)值蒙特卡羅方法進(jìn)行顆粒系統(tǒng)凝并研究，由于其波動(dòng)性較大，結(jié)果均為10～10 000 次運(yùn)算的平均結(jié)果。當(dāng)N值較大(近似無(wú)限系統(tǒng))時(shí)，不同核函數(shù)對(duì)系統(tǒng)混合凝并的影響難以體現(xiàn)，所以在N= 10(或100)的系統(tǒng)(近似有限系統(tǒng))中研究顆粒的凝并情況。當(dāng)λ≥0.5 時(shí)，S.I.的變化規(guī)律值得進(jìn)一步深入研究。