兩組分顆粒系統中凝并混合程度問題的多重蒙特卡羅模擬

姜 志,沈 杰,盧志明

(上海大學上海市應用數學和力學研究所,上海200072)

氣溶膠系統廣泛存在于自然界和工程領域[1-3],通常可以采用顆粒群平衡模擬(population balance modeling,PBM)的方法描述系統中離散相顆粒群的動力學演變過程.顆粒群平衡模擬是通過構建以顆粒尺度分布函數為基本量的顆粒群平衡方程(population balance equation,PBE)來描述顆粒“狀態”,當研究對象是均勻流動或靜止流體的系統時,稱為零維空間系統,顯然零維顆粒群平衡方程模擬是多維顆粒群平衡方程模擬的基礎.

在氣溶膠顆粒系統中,僅僅采用顆粒尺度這一個內部變量難以準確量化顆粒“狀態”,顆粒表面積、組分、表面活化能、不規則度等顆粒群內部變量往往也隨時間和空間不斷演化,而這些顆粒狀態參數在某些領域往往更受關注,如:燃燒生成的煙灰顆粒團聚體的不規則度和表面積與其毒性、吸附性等的關系[4-6];制藥工業中不同組分藥物的凝并混合[7-8]以及結晶過程[9].所以,考慮多組分或多變量動力學演變過程有廣泛的工程和科學背景[2-3,10].Lushnikov[11]首先研究了兩組分離散顆粒系統的凝并問題,并且給出了初始單分散條件下常數核模型的顆粒群平衡方程的解析解,得到在此初值條件下顆粒的組分分布為二項分布.隨后,Gelbard等[12]給出了連續系統的多組分凝并問題PBE的解析解,以4種不同初始條件來闡明凝結和生長對顆粒演化組成分布的影響.Vigil等[13]研究了離散顆粒系統一個復雜加法核模型,得到其顆粒群平衡方程的漸近解,并指出顆粒的尺度分布函數是其中一個組分的顆粒尺度分布函數和高斯函數的乘積.近幾年,Fern′andez-D′?az等[14-15]先后給出了兩組分顆粒系統加法核和乘法核模型PBE的解析解,并對比了常數核、加法核、乘法核顆粒尺度分布函數解析解的表達式,得出它們和相應單組分情況具有相似的形式,以及這3種核對于小顆粒而言,只有乘法核的質量分布函數不隨時間演化.Piskunov[16]給出了兩個動力學事件(凝并和冷凝)同時作用的系統顆粒群平衡方程的解析解.Yu等[17]將文獻[13]的問題推廣到連續顆粒系統,并給出了與組分相關核模型(引入一個組分參數)的顆粒群平衡方程的解析解,分析了參數對顆粒尺度分布函數的影響.以上都是一些關于特殊核模型顆粒群平衡方程解析解的研究,由于PBE是典型的積分微分方程,對于多組分顆粒群而言,在高維空間進行積分微分十分困難,所以對方程的數值求解成為研究者關注的熱點.Matsoukas等[18]和Lee等[19]運用常數目蒙特卡羅(constant-number Monte Carlo)法模擬了兩組分顆粒系統的凝并過程,所考慮的核函數和組分無關,僅與總質量或總體積有關,并且還研究了系統組分之間混合程度參數(分散指數,segregation index,S.I.)的規律,得到分散指數將以一定形式減小.Zhao等[20-21]采用多重蒙特卡羅(multi-Monte Carlo,MMC)方法研究納米顆粒布朗凝并問題,摒棄了Monte Carlo算法中廣泛采用的“子系統”概念,引入了加權“虛擬顆粒”,模擬過程中保持虛擬顆粒數目和計算區域體積不變,從而在計算精度和效率上有很大改進[22-23].對于更復雜的核模型,Matsoukas等[24]提出了一個與組分相關的核模型,通過在核函數上引入組分“凝并效率”函數,使得核模型依賴于組分,但僅僅研究了一些簡單核模型如常數核和加法核的情況,對于更實際且有物理意義的布朗核模型沒有給出相應結論.氣溶膠布朗凝并是氣溶膠顆粒系統中由布朗運動導致顆粒間相對運動而發生相互碰撞黏結在一起的過程,該過程會改變氣溶膠顆粒尺寸分布.而氣溶膠系統的很多物理化學性質,如對光的散射性、毒性、擴散性等,都和顆粒尺寸分布有關.因此,對氣溶膠布朗凝并理論的深入研究有重要的理論和現實意義.蒙特卡羅方法的一個重要特色是能夠以一種簡單而直接擴展的方式來處理高維問題,包括多變量系統.故本工作采用多重蒙特卡羅方法模擬了顆粒布朗凝并問題,主要關注顆粒凝并混合程度等問題,以及不同“凝并效率”參數及初值條件下,兩組分顆粒系統中幾個低階矩統計量和描述不同組分混合程度的參數(即某個組分的總偏差和分散指數)的演變規律.

1 理論和方法

本工作的研究對象是一個顆粒總質量為M、顆粒總數目為N的兩組分顆粒系統,為方便起見,把兩個組分分別記為組分A和組分B.兩變量PBE方程可由Smoluchowski方程推廣而來[18,25]:

式中:r為顆粒群內部變量集,定義r=(v,m),其中m代表顆粒中組分A的質量(或體積),而v代表顆粒的總質量(或體積);F(r,t)表示顆粒群內部變量r的分布函數,所以F(r,t)dr表示t時刻、變量范圍在r～r+dr內的顆粒在單位體積內的數目濃度;K(r,r0)為核(kernel)函數,表示動力學事件發生服從的某種概率.故組分A在所有顆粒中的總質量分數(φ)可以定義為

而組分A在某個顆粒中所占的質量分數為c=m/v.當顆粒系統達到一個完全混合均勻的狀態時,兩種不同質量分數雖然在表達形式上不同,但在數值上是相等的.但是,當系統還沒有到達一個混合均勻的狀態時,二者會有偏差.組分A的偏差定義為x=v(c?φ)=m?φv,則顆粒系統的組分A總偏差定義為[18]

事實上,顆粒系統的組分A的總偏差χ和分散指數S.I.都是描述顆粒系統組分A的混合程度或均勻程度的參考指標.很多領域往往只需要研究顆粒群的幾個重要平均量或整體量,即幾個低階矩來描述顆粒群的動力學演變.兩組分顆粒的(k,l)階的矩定義如下[26]:

顯然M00表示顆粒總數目濃度,M10表示顆粒的總質量濃度,M01表示組分A的總質量濃度,M11,M02,M20與顆粒系統的分散性有關.在引入矩Mkl后,很容易將組分A的總質量分數φ、分散指數S.I.改寫為

由于系統質量守恒,即M10,M01為常數,因此對兩組分問題的凝并演化描述需要計算以下幾個混合矩:M00,M11,M02,M20.

為了求解方程(1),本工作采用Zhao等[20-21]的多重蒙特卡羅方法,該方法發展了異數目權值虛擬顆粒策略,采用數目權值不等的虛擬顆粒群來代表實際物理顆粒群,盡可能多地繼承了多分散顆粒群的尺度分布信息,具有高精度、高效率等優點[22-23].

2 結果與討論

2.1 模擬工況

本工作研究的是顆粒凝并問題,考慮的核模型是由Matsoukas等[24]首次提出的,即

式中,ci表示顆粒i中組分A的質量分數,所以ψ(c1,c2)可以看成在一些常見核函數的基礎上添加一個“凝并效率”函數.圖1給出了“凝并效率”函數中的參數aAB對顆粒系統凝并的影響.由圖1(a)可以發現,當aAB為正值(排斥)時,只含組分A或只含組分B的顆粒發生凝并的效率最高,而既含組分A又含組分B的顆粒發生凝并的效率最低,即顆粒組分不同不利于混合;在圖1(b)中,當aAB為負值(吸引)時,結論恰恰相反,成分不同的顆粒凝并效率較低,有利于均勻混合.值得注意的是,當aAB=0時,核函數K(r1,r2)為原核函數k(v1,v2).Matsoukas等[24]研究了兩個簡單的核函數k(v1,v2)即常數核與加法核,但沒有研究布朗核模型的情形.

圖1 效率函數隨aAB的變化Fig.1 Collision efficiency as a function of aAB

布朗凝并核有兩種形式,顆粒直接在納米級的自由分子區布朗凝并核為(無量綱形式)

顆粒直徑在微米級以上的連續區布朗凝并核為(無量綱形式)

圖2 不同兩組分比例示意圖(φ=0.5,0.4,0.3)Fig.2 Schematic of different proportions of two components(φ=0.5,0.4,0.3)

表1 初始分布(單分散)Table 1 Initial distribution(monodisperse)

2.2 低階矩統計量的演化規律

在單組分系統中,氣溶膠凝并過程存在自保持分布,各階矩(3個)隨時間的推移存在漸近行為[17,27-28],如隨著時間以固定的標度律上升或下降. 在兩組分顆粒系統中,研究的矩統計量有6個,其中兩個矩統計量(M10,M01)保持不變,所以關注的矩統計量有M00,M11,M20,M02.圖3給出了自由分子區和連續區初始分布φ=0.5系統中的無量綱顆粒總數濃度M00/M00ini隨時間的演化,圖中M00ini為初始M00,黑色表示只含組分A的顆粒.由圖3(a)可以看出,對于自由分子區氣溶膠布朗凝并過程,不論aAB取什么值,當凝并時間足夠長時,無量綱顆粒總數目濃度與無量綱時間有如下關系:M00/M00ini～τ?6/5;而由圖3(b)可以看出,對于連續區氣溶膠布朗凝并過程,不論aAB取什么值,當凝并時間足夠長時,M00/M00ini～τ.但是由圖1可以發現,雖然aAB不會影響標度率的大小,但是會影響達到穩定標度率的時間,對于高排斥的“凝并效率”參數aAB=4或5時,需要更久的時間才能觀察到隨時間的標度率關系.同時,本工作也研究了其他幾個無量綱矩統計量M11,M20,M02隨無量綱時間的演化規律,發現其他幾個無量綱矩也具有相同的漸近行為,但是在標度率上有所區別.對于自由分子區氣溶膠布朗凝并過程來說,M11,M20,M02隨時間的標度率為6/5;對于連續區氣溶膠布朗凝并過程來說,M11,M20,M02隨時間的標度率為1.需要指出的是,對于其他初始條件如φ=0.4和0.3下所得的結論和φ=0.5相同,在這里不一一指出.

2.3 組分A總偏差χ和分散指數S.I.的演化規律

下面研究兩組分顆粒系統中,不同組分之間的凝并混合程度問題,其中組分A的總偏差χ和分散指數S.I.都是用來描述系統混合狀態的參數,這里以初始φ=0.5為例,其他初始條件下結論類似.圖4給出了自由分子區和連續區初始分布(φ=0.5)系統中組分A的總偏差χ隨時間尺度的演化過程.可以發現,不論在自由分子區還是連續區布朗凝并過程,當aAB為負值時,χ迅速下降,然后隨著時間的推移趨于一個穩定的值.這與預想情況相一致.由圖1(a)可知:顆粒成分為組分A的(純A)更容易與成分為組分B的(純B)顆粒發生凝并,再加上給定的初始條件為純A顆粒與純B顆粒以1∶1混合,所以在初始階段,混合較快.而當aAB為正值時,χ迅速上升,然后隨著時間的推移也趨于一個穩定的值.這是由在初始狀態下顆粒間的排斥作用導致的.值得注意的是,在排斥效應很強如aAB=4或5時,需要更久的演化才能發現χ趨于穩定值.圖5給出了分散指數S.I.隨時間尺度的演化過程,可以發現在不同aAB條件下,經過足夠長時間的演化后,在自由分子區和連續區,分散指數S.I.和時間尺度存在如下關系:S.I.～

圖3 aAB不同時自由分子區和連續區氣溶膠布朗凝并過程中無量綱顆粒數目M00/M00ini隨無量綱時間τ的演變(φ=0.5)Fig.3 Evolution of normalized M00for Brownian kernel in the continuum regime and free molecular regime with different aAB(φ=0.5)

圖4 aAB不同時自由分子區和連續區氣溶膠布朗凝并過程中組分A總偏差χ隨時間尺度的演變Fig.4 Evolution of power density of excess component A for Brownian kernel in the continuum regime and free molecular regime with di ff erent aABas a function of the time scale

圖5 aAB不同時自由分子區和連續區氣溶膠布朗凝并過程中分散指數S.I.隨時間尺度的演變(φ =0.5)Fig.5 Evolution of S.I.for Brownian kernel in the continuum regime and free molecular regime with di ff erent aABas a function of the time scale(φ=0.5)

由圖4和5可以看出,系統最后都趨于一個穩定狀態.為了進一步研究顆粒從初始階段到穩定階段的過程,以φ=0.5,aAB=4連續區的布朗凝并為例,將顆粒布朗凝并過程分為三個階段(見圖6):第一階段(初始階段),由于初始條件的影響S.I.隨時間尺度迅速上升或下降;第二階段(混合階段)S.I.隨時間尺度呈緩慢下降趨勢,顆粒尺度變寬,排斥作用不占據主導地位,系統向混合均勻狀態推進;第三階段(穩定階段)S.I.隨時間尺度以下降,系統趨于自保持狀態.

圖6 連續區氣溶膠布朗凝并過程中分散指數S.I.隨時間尺度的演變(aAB=4,φ=0.5)Fig.6 Evolution of S.I.for Brownian kernel in the continuum regime as a function of the time scale(aAB=4,φ=0.5)

綜合以上規律,可以得到χ經過長時間演化后會保持一個穩定的值,記為χ∞.事實上,將方程(3)兩邊同時對時間求導,再結合方程(1)可得[24]:

式中,x(ri)=(mi?),由Matsoukas等[18]得到的對于初始單組分的顆粒系統,無論方程的核函數是什么形式,方程(11)的右側都將趨于零,即χ將趨于一個穩定值,這也驗證了本模擬的可靠性.由方程(7)以及圖1可知,隨著aAB的增大,顆粒間的排斥作用越明顯,凝并效率越低,并且結合方程(7)及χ∞的值,猜想χ∞隨aAB存在指數函數關系.由于aAB的正負對系統凝并混合起完全相反的作用,故分別對這兩種情況下χ∞與aAB進行擬合.圖7和8分別給出了組分A總偏差的穩定值χ∞與“凝并效率”參數aAB的函數關系,圖中FR為自由分子區,CR為連續區.由圖可知,在不同的初值條件下,不論自由分子區還是連續區,χ∞與aAB可以用指數函數來確定,具體函數式如下:

具體參數如表2所示(表中Ci為擬合參數),可以發現隨著aAB的增大,χ趨于穩定的值將呈指數型增長.從增長趨勢來看,不論在自由分子區還是連續區,當aAB為負值即組分相同時不利于混合,且隨著aAB的繼續增大,初始組分A的質量分數越大,χ∞值越大,顆粒系統質量分布越不均勻;當aAB為正值即組分相同時有利于混合,且隨著aAB的繼續增大,初始組分A的質量分數越大,χ∞值越小,顆粒系統質量分布越均勻.此外,由圖4可知,當aAB=4或5時,需要消耗大量的計算機成本才能發現χ趨于穩定值,方程(12)可以預測在排斥效應比較強時組分A總偏差的穩定值,相應地可以預測顆粒系統的分布均勻程度.

圖7 不同初值條件下自由分子區氣溶膠布朗凝并過程中組分A總偏差χ∞與aAB的關系Fig.7 Dependence of power density of excess component A,χ∞ on aABfor Brownian kernel in free molecular regime

圖8 不同初值條件下連續區氣溶膠布朗凝并過程中組分A總偏差χ∞與aAB的關系Fig.8 Dependence of power density of excess component A,χ∞ on aABfor Brownian kernel in continuum regime

表2 不同初始單分散工況下的擬合參數Table 2 Coefficients of the fitted curves of different initial conditions

3 結束語

本工作利用多重蒙特卡羅方法模擬了兩組分顆粒系統中與組分有關的核的布朗凝并問題.核模型中引入了“凝并效率”參數aAB,aAB值的正負會加快或者阻礙顆粒凝并過程.在某種機制或某幾種機制作用下,顆粒群尺度譜將達到自保持分布[17,27-28],即顆粒尺度分布函數的所有特征變量在進行無量綱后具有穩定的表達式.本工作發現在自由分子區和連續區布朗核下,經充分長時間的演化后兩組分顆粒系統的矩統計量也存在漸近行為,并且在自由分子區,M00～ τ?6/5,M11,M20,M02～ τ6/5;在連續區,M00～ τ?1,M11,M20,M02～ τ.當系統達到自保持狀態時,系統混合度達到一個穩定值并且分散指數以下降,根據分散指數的發展過程可將系統演變分為三個階段,即初始階段、混合階段、穩定階段.在初始階段,由于凝并效率參數的影響,S.I.隨時間尺度迅速上升或下降;隨著時間的推移,S.I.隨時間尺度呈緩慢下降趨勢,顆粒尺度變寬,排斥作用不占據主導地位,系統向混合均勻狀態推進;經長時間演化后,S.I.隨時間尺度以下降,系統趨于自保持狀態.反之,當S.I.隨時間尺度以1/v下降時,認為此時系統達到穩定階段或者處于自保持狀態.為了研究當系統達到自保持狀態下組分A的不均勻性,通過擬合得到不同初值條件下組分A的總偏差的穩定值與“凝并效率”參數aAB的函數關系,發現二者之間為指數函數關系,此關系式可以用來預測排斥作用比較強的情況下組分A總偏差的穩定值.