吳志豪， 凡鳳仙， 胡曉紅

(上海理工大學能源與動力工程學院， 上海200093)

氣固兩相流廣泛存在于能源、環(huán)境、化工等領域,如煤粉的氣力輸送[1]、 固體燃料顆粒在燃燒室內(nèi)的流動與燃燒[2]、 細顆粒物在除塵器內(nèi)的流動與分離[3]、 大氣顆粒物的擴散和演化[4]、 流化床內(nèi)顆粒物料的流化[5-6]等。近年來，隨著計算機硬件的更新?lián)Q代和數(shù)值模擬技術的迅速發(fā)展，計算機數(shù)值模擬已成為研究氣固兩相流過程和機理的重要手段。建立可靠的描述顆粒運動、 碰撞、 凝并及反彈等物理過程的顆粒動力學模型，對氣固兩相流的高精度模擬顯得非常重要和迫切。離散元方法(discrete element method，DEM)是顆粒動力學行為建模和模擬的最直接的方法，該方法基于牛頓第二定律，根據(jù)顆粒與顆粒之間的相互作用力和顆粒與流體之間的耦合作用來求解每個顆粒的速度與運動軌跡，通過跟蹤所有顆粒的運動軌跡來判斷顆粒間碰撞的發(fā)生。 盡管DEM對顆粒動力學過程的描述最為接近實際， 但其在處理涉及大量顆粒的氣固兩相流時， 尤其是當固相為粉體顆粒時，計算成本極高。直接模擬蒙特卡洛(direct simulation Monte Carlo, DSMC)方法是對氣固兩相流中顆粒動力學過程進行建模和模擬的另一種方法。DSMC方法在硬球模型的基礎上結合了判斷顆粒碰撞的蒙特卡洛方法，在稀相氣固兩相流數(shù)值模擬中得到廣泛應用。該方法立足于氣體分子運動論，采用取樣顆粒代替真實顆粒，為每個取樣顆粒賦予一定的權重，使其代表一定數(shù)目的真實顆粒，跟蹤取樣顆粒在流場中的受力與運動，通過概率判斷顆粒碰撞的發(fā)生，從而有效地節(jié)省了計算時間與計算機內(nèi)存，適宜于對大量顆粒構成的稀疏氣固兩相流中顆粒動力學模擬。

DSMC方法在氣固兩相流模擬的實際應用中往往將顆粒碰撞的結果處理為單一事件，即所有碰撞都導致顆粒反彈或都導致顆粒凝并。例如，在顆粒碰撞導致顆粒反彈的假設下，Tsuji等[7]對循環(huán)流化床上升管內(nèi)的團簇分布進行了數(shù)值預測，F(xiàn)an等[8]研究了駐波聲場中顆粒的碰撞頻率，F(xiàn)ede等[9]探究了湍流中的顆粒碰撞。在相關顆粒碰撞-反彈的研究中，往往對所有取樣顆粒采用相同的數(shù)目權重，這種等權重算法便于應用，但當研究對象為寬分布多分散顆粒體系或為含有少量特定成分的混合顆粒體系時，等權重算法將可能導致代表少量顆粒的取樣顆粒數(shù)目不足，甚至缺失，從而造成統(tǒng)計噪聲增加，計算精度降低。為了解決這一問題，He等[10]研究出一種處理顆粒碰撞-反彈行為的異權重DSMC方法。DSMC方法在黏性顆粒凝并過程模擬上得到廣泛應用，目前已發(fā)展出常體積蒙特卡洛(constant-volume Monte Carlo，CVMC)方法[11]、 常數(shù)目蒙特卡洛(constant-number Monte Carlo，CNMC)方法[12]、 多重蒙特卡洛(multi-Monte Carlo，MMC)方法[13]。其中，CVMC方法和CNMC方法均采用等權重算法。CVMC方法保持計算區(qū)域不變，隨著顆粒凝并的進行，取樣顆粒數(shù)目不斷減少，使得計算精度不斷降低，因此，為了保證計算精度，往往采用大量取樣顆粒，導致計算成本很高。CNMC方法通過增大計算區(qū)域體積以保持取樣顆粒數(shù)目不變，但該方法難以考慮顆粒的空間擴散和計算區(qū)域的邊界條件，使其應用范圍受到限制。針對這2種方法的不足，Zhao等[13]提出了一種異權重算法，即MMC方法，該方法在不改變計算區(qū)域體積的同時，通過調(diào)整取樣顆粒的數(shù)目權重，以保持計算過程中取樣顆粒數(shù)目不變。MMC方法應用于亞微米顆粒的布朗凝并和聲凝并建模和模擬，展現(xiàn)出對亞微米顆粒凝并行為優(yōu)良的數(shù)值預測能力[13-16]。

盡管DSMC方法在氣固兩相流數(shù)值模擬中發(fā)揮著重要作用，并已得到廣泛應用，然而該方法在對碰撞結果的處理上存在一定的局限性。分析顆粒碰撞的物理過程可以發(fā)現(xiàn)，大粒徑和高碰撞速度將會導致顆粒反彈，而小粒徑和低碰撞速度條件下顆粒碰撞的結果則傾向于凝并[17]。當對粒徑分布寬廣、 氣相場復雜的氣固兩相流進行建模與模擬時，如果仍對顆粒碰撞的后果進行單一化處理(全部為反彈或凝并)，將導致數(shù)值預測結果與實際有較大的偏差。為解決這一問題，本文中對已有研究進行繼承和發(fā)展，構建了一種能夠求解顆粒碰撞后凝并與反彈共存的問題的DSMC方法，并基于該方法對寬粒徑雙峰分布的燃煤飛灰聲凝并效果進行預測，獲得了與實驗吻合良好的結果，驗證了模型的可靠性。

1 顆粒凝并與反彈共存的DSMC方法

1.1 顆粒碰撞概率

DSMC方法認為顆粒碰撞為隨機的二元碰撞，并利用概率的方法來判斷顆粒碰撞的發(fā)生。通常對計算區(qū)域劃分網(wǎng)格，認為只有處于同一網(wǎng)格內(nèi)的顆粒才可能發(fā)生碰撞，則取樣顆粒i與取樣顆粒j所代表的真實顆粒發(fā)生碰撞的概率[16]可寫為

Pij=βijwjΔtcol/Vci，

(1)

式中：Pij為碰撞概率；βij為碰撞核函數(shù)， m3/s；wj為取樣顆粒j的數(shù)目權重； Δtcol為求解顆粒間碰撞的時間間隔， s；Vci為取樣顆粒i所在網(wǎng)格的體積， m3。

取樣顆粒i和所有顆粒的碰撞概率Pi可寫為

(2)

式中，N為取樣顆粒i所在網(wǎng)格內(nèi)取樣顆?？倲?shù)目。

1.2 顆粒碰撞的判斷

采用修正的南布方法[7]判斷顆粒碰撞事件的發(fā)生。在Pij<1/N的前提下，利用計算機產(chǎn)生的[0, 1)范圍均勻分布的隨機數(shù)R，選擇同一網(wǎng)格內(nèi)的取樣顆粒j作為候選被碰顆粒：

j=int[RN]+1，

(3)

式中，int[RN]表示RN的整數(shù)部分。

如果滿足關系式

R>j/N-Pij，

(4)

則認為取樣顆粒i與j之間發(fā)生碰撞。

1.3 顆粒碰撞后果的確定

顆粒碰撞后果由顆粒碰撞后的動能和顆粒間的黏附力共同決定。對于不帶電的固體干顆粒，范德華力為顆粒間唯一的黏附力。基于能量平衡關系式，可得臨界碰撞速度vc的表達式[17]為

(5)

式中：A為Hamaker常數(shù)， J；e為恢復系數(shù)；z0為最小接觸距離， m；pc為極限接觸壓力， Pa；ρ為顆粒密度， kg/m3；ds和dl分別為較小和較大顆粒的直徑，m。

當碰撞顆粒之間的法向相對速度不大于臨界速度，即滿足式(6)時，碰撞后果為凝并；反之為反彈。

(6)

式中：v為顆粒速度， m/s；n為由顆粒i指向顆粒j的法向單位向量； 上標0表示碰前值。

1.4 顆粒凝并的處理

若取樣顆粒i和j發(fā)生凝并，則意味著它們代表的真實顆粒發(fā)生凝并。分析顆粒凝并的物理過程可知，真實顆粒發(fā)生碰撞-凝并的次數(shù)為min(wi,wj)，形成了min(wi,wj)個顆粒團聚體，同時數(shù)目權重較大的取樣顆粒所代表的真實顆粒中有|wi-wj|個未參與碰撞-凝并過程。根據(jù)質(zhì)量和動量守恒定律，可得：

(7)

(8)

式中：m為顆粒質(zhì)量， kg；Vagg為團聚體體積， m3；vagg為團聚體的速度。

對顆粒凝并的處理方案為保持取樣顆粒數(shù)目不變，調(diào)整取樣顆粒的數(shù)目權重、體積和速度以體現(xiàn)碰撞-凝并的發(fā)生[13-15]，具體如下：

當wi=wj=2l(l=1, 2, 3, …)時，

(9)

當wi=wj=2l-1(l=1, 2, 3, …)時，

(10)

當wi>wj時，

(11)

當wi

(12)

式中,上標new表示碰撞后。

1.5 顆粒反彈的處理

對于每次真實顆粒碰撞-反彈過程，顆粒反彈速度可根據(jù)動量方程計算。若不考慮顆粒轉(zhuǎn)動，并認為顆?；七^程所受的摩擦力遵守庫侖摩擦定律，且顆粒一旦停止滑移，不會再次發(fā)生滑移，顆粒反彈速度[7]為：

(13)

(14)

(15)

(16)

(17)

(18)

分析碰撞-反彈過程可知，對于等權重的取樣顆粒的碰撞-反彈，僅需按式(13)—(16)更新顆粒速度，以表征顆粒反彈的結果。對于不等權重的兩取樣顆粒的碰撞-反彈，采用“分離-合并”方法處理碰撞-反彈過程[10]。 以wi>wj時為例，真實顆粒碰撞-反彈次數(shù)為wj，即取樣顆粒i所代表的wj個真實顆粒與取樣顆粒j所代表的所有真實顆粒發(fā)生碰撞并產(chǎn)生反彈；而取樣顆粒i所代表的剩余的(wi-wj)個真實顆粒未參與到碰撞-反彈過程中?；诖?，首先將取樣顆粒i分離為子取樣顆粒i*和i**，其數(shù)目權重分別為wj和(wi-wj)；接著，按照等權重取樣顆粒碰撞-反彈處理方法，獲得i*與j的反彈速度；最后，將碰撞-反彈后的子取樣顆粒i*和不參與碰撞-反彈過程的子取樣顆粒i**進行合并，并將合并后的顆粒重新存儲為取樣顆粒i，其數(shù)目權重恢復為wi，其速度可由動量守恒定律計算：

(19)

式中： 上標mer表示合并后；下標i*和i**分別表示顆粒i*和i**。

在合并過程中，不可避免地產(chǎn)生了能量損失，即

(20)

式中, ΔEi為能量損失。

為保持合并前后系統(tǒng)的能量不變，將能量損失補償?shù)絽⑴c后續(xù)碰撞的某一取樣顆粒s上，根據(jù)能量平衡方程，補償后顆粒s的速度為

(21)

式中, 上標com表示補償后。

為在能量補償對系統(tǒng)擾動盡量小的前提下能夠快速搜尋到顆粒s，采用的搜索規(guī)則為

(22)

2 顆粒凝并與反彈共存的DSMC方法在聲凝并中的應用

聲凝并是指利用高強度外加聲場作用促使顆粒發(fā)生相互作用，進而碰撞接觸，凝并為粒徑較大顆粒的現(xiàn)象。聲凝并作為富有工業(yè)應用潛力的顆粒預處理技術，為燃燒源、工業(yè)源細顆粒物的高效脫除提供了一種重要技術途徑。聲凝并過程中涉及的顆粒既有亞微米尺度黏附顆粒，又有微米尺度慣性顆粒，適用于本文中建立的顆粒凝并與反彈共存的DSMC方法的驗證。

2.1 取樣顆粒運動方程

在聲場中，顆粒除發(fā)生碰撞外，還將受到曳力、重力、浮力作用發(fā)生運動，顆粒i的運動方程[18-19]可寫為：

(23)

式中：μg為氣體的動力黏度， Pa·s；u為氣體速度， m/s；ρg為氣體密度， kg/m3；g為重力加速度， m/s2；Cine為慣性修正系數(shù)；CCun為肯寧漢滑移修正系數(shù)。

(24)

CCun=1+Kn[1.257+0.400exp(-1.100/Kn)]，

(25)

式中:Kn為克努森數(shù);Re為雷諾數(shù)。 其表達式分別為

(26)

(27)

式中,λg為氣體分子平均自由程，m。

水平平面行波聲場的波動方程為

ux=uampsin[ω(t-x/c)]，

(28)

式中：ux為聲波引起的氣體介質(zhì)振動速度，m/s；uamp為聲波速度振幅，m/s；ω為聲波的角頻率， rad/s；t為時間， s；x為距聲源的距離， m；c為聲速，m/s。

通常采用聲壓級L(單位dB)描述聲場，聲波速度振幅可根據(jù)聲壓級計算，

(29)

式中,Pr為參考聲壓，Pr=2×10-5Pa。

2.2 碰撞核函數(shù)

碰撞核函數(shù)是碰撞概率求解中的一個重要參數(shù)(見式(1))， 用于描述顆粒相互作用機理對顆粒碰撞的作用效果， 定義為單位時間內(nèi)顆粒因某種機理而能夠與其它顆粒發(fā)生碰撞的體積。 聲凝并中的碰撞核函數(shù)包括同向相互作用、 重力沉降作用、 聲尾流效應、 互輻射壓力效應以及布朗擴散效應核函數(shù)。

顆粒i和j之間的同向相互作用和重力沉降作用碰撞核函數(shù)βij，Orth-Grav可基于顆粒i和j運動方程的解vi和vj[16]給出，

(30)

聲尾流效應碰撞核函數(shù)βij, AWE[20]可表示為:

(31)

式中：r為顆粒i與j的中心距， m；li和lj分別為顆粒i和j的滑移系數(shù)。

li的表達式為

(32)

其中，

(33)

(34)

(35)

(36)

互輻射壓力效應碰撞核函數(shù)βij,MRPE的表達式[21]為

(37)

式中：gijr和gjir為顆粒i與顆粒j之間的流體力學相互作用函數(shù)，具體見文獻[21]。

布朗擴散效應碰撞核函數(shù)βij,Brow可表示[22]為

(38)

式中：T為溫度， K；kB為玻爾茲曼常數(shù)，kB= 1.38×10-23J/K。

總的碰撞核函數(shù)βij可寫為

βij=βij,Orth-Grav+βij,AWE+βij,MRPE+βij,Brow

。

(39)

2.3 聲凝并數(shù)值模擬條件和流程

2.3.1 時間步長

求解顆粒運動方程的時間步長Δt[18]應滿足

(40)

(41)

在求解顆粒碰撞的時間步長Δtcol的選擇上，考慮到修正的南布方法的要求Pij<1/N，則有

(42)

基于上述選擇原則，為加速計算過程，本文中采用Δtcol=100Δt，即顆粒運動每求解100步，顆粒碰撞進行搜索和求解一遍。

2.3.2 初始條件和邊界條件

沿聲波傳播方向取一個波長λ范圍即λ×0.01 m×0.01 m的三維區(qū)域為計算區(qū)域，并將其沿x向劃分成160個網(wǎng)格。依據(jù)每個粒徑區(qū)間內(nèi)顆粒的數(shù)目濃度，對取樣顆粒的初始數(shù)目權重進行分配。將取樣顆粒隨機均勻地添加到計算區(qū)域，并為每個顆粒分配一個隨機的微小初始速度。采用循環(huán)邊界條件處理流出流進計算區(qū)域的顆粒。

2.3.3 計算流程

數(shù)值計算流程如下：

1)初始化取樣顆粒場。基于初始粒徑分布產(chǎn)生取樣顆粒，為每個顆粒賦予一定的數(shù)目權重、位置和速度。

2)求解取樣顆粒的運動。采用四階龍格-庫塔算法，基于2.1節(jié)給出的顆粒運動方程，計算每一個取樣顆粒經(jīng)過一個時間步長Δt后的速度和位置，并基于循環(huán)邊界條件對流出計算區(qū)域的顆粒進行處理。該步驟重復100次。

3)判斷顆粒的碰撞并對碰撞后果進行處理。按照第1部分結合第2.2節(jié)的碰撞核函數(shù)，對顆粒的碰撞及碰撞后的凝并與反彈進行處理。

4)重復第2)、3)步的計算，直到達到設定的聲凝并時間。

2.4 數(shù)值模擬與實驗結果的對比驗證

為驗證模型的準確性，采用和Zhang等[23]的實驗報道一致的燃煤飛灰顆粒粒徑分布和物性參數(shù)、氣相參數(shù)、聲場條件和聲波作用時間，開展數(shù)值模擬。數(shù)值模擬參數(shù)見表1。

表1 數(shù)值模擬參數(shù)

圖1給出了數(shù)值模擬得到的聲凝并后顆粒粒徑分布與實驗結果的對比情況。其中，數(shù)值模擬結果包括2種情況：1)將顆粒碰撞結果全部視為凝并；2)將顆粒碰撞結果視為凝并和反彈共存的情況。由圖可見，將顆粒碰撞結果全部視為凝并時，數(shù)值模擬得到的聲凝并后顆粒粒徑分布與實驗差異顯著；當利用臨界碰撞速度將顆粒碰撞結果考慮為凝并和反彈共存時，所得結果和實驗吻合良好。這表明，對于粒徑涉及亞微米至數(shù)微米的飛灰顆粒，聲凝并過程中顆粒碰撞后的反彈不能忽略，此時若將顆粒碰撞后果全部視為凝并，將導致數(shù)值模擬的預測結果優(yōu)于實際情況，從而很大程度上高估了聲凝并的效果。

a)粒徑分布b)模擬結果的局部放大圖1 數(shù)值模擬得到的聲凝并后顆粒粒徑分布與實驗結果的對比Fig.1 Comparison of particle size distribution after acoustic agglomeration obtained from numerical simulation with experimental data

由圖1中數(shù)值模擬結果可知，聲凝并后粒徑在10 μm以下顆粒的數(shù)量濃度降低，同時出現(xiàn)了粒徑在10 μm以上的顆粒，這些新顆粒是由聲凝并形成的團聚體。在外加聲場作用下，粒徑在2.5 μm以下的細顆粒因慣性很小而幾乎被聲波完全夾帶，使得細顆粒間相對運動很弱，碰撞概率很低，凝并效果很差；粒徑在2.5 μm以上的粗顆粒，慣性逐漸增強，而越來越難以被夾帶，使得粗顆粒間相對運動較明顯，碰撞概率大，凝并效果也較好[16,19]。同時，細顆粒和粗顆粒之間夾帶程度的不同導致的顯著的顆粒間相對運動，也使得細顆粒易于與粗顆粒發(fā)生碰撞，從而細顆粒數(shù)目減小，因此，聲凝并過程中生成的10 μm以上的大顆粒往往是由捕獲了細顆粒的粗顆粒之間發(fā)生凝并或粗顆粒間凝并后再捕獲細顆粒形成的團聚體。在今后的研究中，可進一步挖掘DSMC方法的優(yōu)勢，在數(shù)值模擬過程中追蹤顆粒團聚體的構成信息，以直觀展示顆粒團聚體的構成和演變特性。此外，需要說明的是：實驗結果未能給出10 μm以上顆粒的粒徑分布，這是由于實驗中采用靜電低壓沖擊器(ELPI)測量聲凝并前后顆粒的粒徑分布，而ELPI的測量范圍是0.007～10 μm[23]。

為了進一步說明數(shù)值模擬結果的可靠性，圖2給出了聲凝并前后顆粒的累積質(zhì)量濃度分布。由圖可見，數(shù)值模擬前后顆粒的總質(zhì)量濃度相等，這表明數(shù)值模擬中顆粒質(zhì)量是守恒的。綜合圖1和圖2可知，本文中建立的求解顆粒凝并與反彈的DSMC方法用于聲凝并過程模擬能夠準確預測出顆粒粒徑分布的演變。該模型為準確、高效模擬包含寬粒徑分布顆粒的稀疏氣固兩相流行為提供了重要基礎。

圖2 聲凝并前后顆粒的累積質(zhì)量濃度分布Fig.2 Cumulative particle mass concentration distribution before and after acoustic agglomeration

3 結論

DSMC方法作為一種高效、經(jīng)濟的數(shù)值模擬方法，在稀疏氣固兩相流數(shù)值模擬研究中得到了廣泛應用，然而，已有相關研究往往認為顆粒碰撞的結果全部是顆粒反彈，或者全部是顆粒凝并，使得DSMC方法在預測復雜氣相場或?qū)捔椒植碱w粒條件下的氣固兩相流時，具有明顯的局限性。為解決已有研究中存在的問題，本文中做了以下工作。

1)基于顆粒碰撞后發(fā)生凝并和反彈的臨界速度，以及異權重DSMC方法的基本思想，發(fā)展了顆粒凝并與反彈共存的DSMC方法。

2)基于所發(fā)展的方法，建立了顆粒聲凝并模型，提出了相應的數(shù)值模擬方法。

3)對聲凝并效果進行了數(shù)值模擬，發(fā)現(xiàn)對顆粒碰撞結果進行單一凝并處理時，數(shù)值模擬得到的聲凝并后顆粒粒徑分布與實驗結果存在顯著差異；而將顆粒凝并后果處理為凝并與反彈共存時，則數(shù)值模擬與實驗結果吻合良好。