胡善偉，劉新華

（中國科學(xué)院過程工程研究所，多相復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100190）

引 言

氣固流化床內(nèi)的氣體、顆粒與壁面之間的相互作用導(dǎo)致該類系統(tǒng)通常呈現(xiàn)復(fù)雜的時(shí)空多尺度結(jié)構(gòu)。隨著流化氣速從低到高的變化，氣固系統(tǒng)可能會(huì)依次出現(xiàn)鼓泡、湍動(dòng)、快速流態(tài)化以及稀相輸送等跨流域流動(dòng)特征[1]。在系統(tǒng)內(nèi)部，局部非均勻性表現(xiàn)為氣泡和團(tuán)聚物的動(dòng)態(tài)生成和湮滅，具有時(shí)間相依性；而整體非均勻性則表現(xiàn)為反應(yīng)器在軸向和徑向呈現(xiàn)有序分布的特征，具有空間相依性[1]。介尺度結(jié)構(gòu)被認(rèn)為是流化床反應(yīng)器內(nèi)多尺度非均勻性產(chǎn)生的根源，但目前人們對其認(rèn)識(shí)程度仍不能滿足實(shí)際工藝設(shè)計(jì)的需要。深入理解氣固兩相流介尺度動(dòng)態(tài)結(jié)構(gòu)特性是流化床反應(yīng)器定量放大和優(yōu)化設(shè)計(jì)的基礎(chǔ)和前提。

針對氣固流態(tài)化系統(tǒng)，無論是顆粒解析的直接數(shù)值模擬（DNS）或是流體解析的離散顆粒模型（DPM），均由于計(jì)算量巨大而難以應(yīng)用于工業(yè)尺度反應(yīng)器的模擬和優(yōu)化。因此，基于粗?；╟oarsegrained）的離散顆粒模擬和基于粗網(wǎng)格（coarsegrid）的連續(xù)介質(zhì)模擬仍然是現(xiàn)階段主要的工業(yè)數(shù)值模擬策略。粗網(wǎng)格處理帶來的顯著問題是平均方法的失效，此時(shí)必須考慮非解析的亞網(wǎng)格結(jié)構(gòu)對相間作用力和固相應(yīng)力等本構(gòu)關(guān)系的貢獻(xiàn)[2-3]，且相間曳力關(guān)系的合理構(gòu)建是模擬取得成功的關(guān)鍵[2-6]。目前主要有三種相間曳力建模方法。（1）經(jīng)驗(yàn)參數(shù)修正法[5-8]。通過對比實(shí)驗(yàn)和模擬數(shù)據(jù)，逆推出曳力系數(shù)的修正關(guān)系。例如，根據(jù)實(shí)驗(yàn)獲得的最小流化速度及空隙率對傳統(tǒng)均勻曳力模型進(jìn)行校正[7]；或采用恒定的團(tuán)聚物直徑（dcl）替代均勻曳力模型中的單顆粒直徑（dp）來計(jì)算實(shí)際曳力系數(shù)[8]等。這類方法的主要缺點(diǎn)是缺乏對物理機(jī)制的分析，因而模型的可拓展性較差。（2）濾波函數(shù)法[3,9-15]。通過濾波函數(shù)對守恒方程各項(xiàng)進(jìn)行卷積計(jì)算，以獲得空間平均化的控制方程。濾波方程產(chǎn)生的余項(xiàng)來源于局部非均勻結(jié)構(gòu)，需要采用亞網(wǎng)格模型進(jìn)行封閉。以曳力系數(shù)βe(x1,x2,…)為例，一般可通過高分辨率的連續(xù)介質(zhì)模擬、離散顆粒模擬或直接數(shù)值模擬獲得的數(shù)據(jù)進(jìn)行關(guān)聯(lián)[9-15]。根據(jù)研究體系或選擇參量集xi的不同，眾多的βe(x1,x2,…)關(guān)聯(lián)函數(shù)被相繼提出。雖然它們最終形式差別很大，但一般均是濾波特征尺寸、固相分率、滑移速度等變量的函數(shù)[9-17]。近期還進(jìn)一步在模型中考慮了材料性質(zhì)、各向異性、壓力梯度等因素的影響[18-20]。如何選擇合適的標(biāo)記參量（marker）組合是目前研究的熱點(diǎn)之一[20-21]。與傳統(tǒng)濾波方法不同，基于湍流的建模方法則是在質(zhì)量和動(dòng)量守恒方程的基礎(chǔ)上嚴(yán)格推導(dǎo)本構(gòu)關(guān)系[22-23]，并通過對濾波余項(xiàng)的權(quán)重分析簡化封閉關(guān)系的表達(dá)。（3）介尺度結(jié)構(gòu)建模法[24-31]。通過在亞網(wǎng)格尺度上對局部非均勻流動(dòng)結(jié)構(gòu)進(jìn)行分解，以氣泡或團(tuán)聚物為研究對象建立一系列守恒關(guān)系，并在附加約束條件下通過解析求解或數(shù)值求解方法獲得曳力系數(shù)的顯式表達(dá)。典型的代表是能量最小多尺度（EMMS）曳力模型[25-26]。它將系統(tǒng)分解為稀密兩相并分別建立質(zhì)量和動(dòng)量守恒方程，進(jìn)一步引入團(tuán)聚物方程和穩(wěn)定性約束條件，實(shí)現(xiàn)模型的封閉和數(shù)值求解[1,27]。

與經(jīng)驗(yàn)方法及濾波模型不同，EMMS 模型通過量化局部非均勻結(jié)構(gòu)來更直接地表征操作條件或流域變化時(shí)局部介尺度結(jié)構(gòu)的演變及其對曳力的影響。EMMS模型的提出最初是為了計(jì)算快速流化床的整體非均勻結(jié)構(gòu)。為了與計(jì)算流體動(dòng)力學(xué)（CFD）相結(jié)合，Yang 等[26]在EMMS 模型中引入了顆粒加速度并將其應(yīng)用于計(jì)算網(wǎng)格，從而給出了一種基于結(jié)構(gòu)的EMMS 曳力系數(shù)表達(dá)式。此后，該思想被研究者廣泛采用，許多改進(jìn)模型被相繼提出以適應(yīng)更廣泛的操作條件[28-37]。由于原模型中團(tuán)聚物方程在低顆粒通量、強(qiáng)慣性力和負(fù)顆粒速度情況下并不適用（可能導(dǎo)致dcl

盡管EMMS模型有效考慮了局部非均勻結(jié)構(gòu)對曳力本構(gòu)關(guān)系的影響，但既往研究大多圍繞靜態(tài)的氣泡和團(tuán)聚物建立模型，網(wǎng)格內(nèi)介尺度結(jié)構(gòu)特征在特定操作參數(shù)和局部空隙率情況下為一定值。這樣的處理方式忽略了氣泡和團(tuán)聚物動(dòng)態(tài)演化的歷史信息，很大可能導(dǎo)致計(jì)算網(wǎng)格內(nèi)實(shí)際的介尺度結(jié)構(gòu)信息和EMMS 預(yù)測值不匹配，因此有必要在模型中進(jìn)一步引入時(shí)間變量并考慮氣泡和團(tuán)聚物的時(shí)空動(dòng)態(tài)演化。群平衡模型（PBM）考慮了離散相的聚并破碎動(dòng)力學(xué)對數(shù)密度函數(shù)的影響，已經(jīng)被廣泛應(yīng)用于多相體系中離散相尺寸分布的時(shí)空演化特征表征[59-61]。Hu 等[62-63]通過量化稀密相質(zhì)量交換建立了團(tuán)聚物的連續(xù)生長模型以及低通量床中的團(tuán)聚物群平衡方程，并在此基礎(chǔ)上完善了連續(xù)介質(zhì)模型、EMMS 曳力和群平衡模型的耦合框架，隨后還拓展了PBM 和EMMS 耦合模型在鼓泡流化床中的應(yīng)用。由于考慮了介尺度結(jié)構(gòu)時(shí)空動(dòng)態(tài)演化行為對局部氣泡和團(tuán)聚物尺寸分布的影響，所以PBM 模型的引入使得EMMS 曳力的修正結(jié)果更加可靠。

計(jì)算流體動(dòng)力學(xué)模擬方法在被應(yīng)用于工業(yè)反應(yīng)器時(shí)可能會(huì)難以兼顧效率、精度和成本，但通過優(yōu)化模擬初場、從特定的整體分布開始迭代可以改善上述問題[64-65]。一方面，相對顆粒均勻分布和緊密堆積的初始狀態(tài)，通過宏尺度模型對擬穩(wěn)態(tài)流場進(jìn)行預(yù)測無疑更符合實(shí)際工況[30]。另一方面，人們有時(shí)更關(guān)注反應(yīng)器在不同操作條件下的整體表現(xiàn)，采用大批量的數(shù)值模擬成本高昂且效率低下，此時(shí)穩(wěn)態(tài)建模也可以為反應(yīng)器的設(shè)計(jì)優(yōu)化提供直接依據(jù)[66-70]。然而，氣固流化床操作方式多變，且伴隨著狀態(tài)多值性和流域轉(zhuǎn)變等復(fù)雜的非線性跨流域特征。隨操作氣速的變化，反應(yīng)器內(nèi)的氣固流動(dòng)還可能經(jīng)歷從移動(dòng)床到快速床等連續(xù)演化過程。這些特征給流化床反應(yīng)器的宏尺度建模帶來了挑戰(zhàn)[1]。經(jīng)驗(yàn)的穩(wěn)態(tài)模型往往很難捕捉流域轉(zhuǎn)變過程中出現(xiàn)的“噎塞”和“液泛”等跨流域特征，但EMMS 思想為氣固系統(tǒng)的跨流域統(tǒng)一建模提供了理論基礎(chǔ)[65]。Hu等[30,67]通過考慮提升管反應(yīng)器內(nèi)顆粒加速效應(yīng)和壁面的影響，在局部和全局穩(wěn)定性條件的約束下將EMMS模型進(jìn)行了軸徑向二維擴(kuò)展。通過分析不同流域條件下的控制機(jī)制及穩(wěn)定性條件表達(dá)，相繼提出了鼓泡EMMS穩(wěn)態(tài)模型[45]和并/逆流EMMS穩(wěn)態(tài)模型[68-69]。在此基礎(chǔ)上，Liu 等[69-70]建立了基于EMMS的全流域統(tǒng)一建模方法并繪制了操作相圖。Tu等[71]進(jìn)一步指出，當(dāng)對多回路循環(huán)流化床系統(tǒng)進(jìn)行全回路模擬時(shí)，必須恰當(dāng)?shù)乜紤]各子反應(yīng)器模塊的邊界條件及單元操作之間的強(qiáng)耦合作用。EMMS操作相圖為循環(huán)流化床反應(yīng)器內(nèi)的流域識(shí)別和量化表征提供了可行途徑。通過考慮各模塊之間的壓降平衡和全床物料守恒關(guān)系，所謂的EMMS 全回路建模方法可以實(shí)現(xiàn)對該類反應(yīng)器整體穩(wěn)態(tài)動(dòng)力學(xué)的準(zhǔn)實(shí)時(shí)預(yù)測[70,72]。

本文將重點(diǎn)介紹和總結(jié)本課題組在氣固系統(tǒng)介尺度建模和EMMS模型擴(kuò)展方面的工作。首先介紹基于團(tuán)聚物或氣泡的EMMS曳力模型的改進(jìn)和發(fā)展；第二部分論述介尺度結(jié)構(gòu)時(shí)空動(dòng)態(tài)演化行為對本構(gòu)關(guān)系的影響，以及群平衡模型在氣固流態(tài)化系統(tǒng)中的應(yīng)用；第三部分介紹基于EMMS 的反應(yīng)器宏尺度建模方法，包括模型的擴(kuò)展、相圖繪制和全回路建模等；最后是對未來工作的展望。

1 EMMS曳力模型

如圖1 所示，EMMS 模型將非均勻氣固系統(tǒng)分解為顆粒富集的密相和氣體富集的稀相，兩者之間通過相界面相互作用[1]。假設(shè)兩相內(nèi)的氣固相互作用可以用均勻模型描述，則特定操作條件系統(tǒng)狀態(tài)可以用8個(gè)參數(shù)X={Upf,Uf,εf,Upc,Uc,εc,f,dcl}來定義。通過對稀密兩相分別建立質(zhì)量和動(dòng)量守恒方程，再加上兩者之間的壓降平衡方程及團(tuán)聚物直徑關(guān)聯(lián)式，總共6 個(gè)守恒方程得以建立。Li 等[1,27]認(rèn)為在氣固流態(tài)化系統(tǒng)存在兩種控制機(jī)制：當(dāng)流體控制時(shí)，Wst=min；顆?？刂茣r(shí)，ε=min。兩者的競爭協(xié)調(diào)關(guān)系構(gòu)成了系統(tǒng)的穩(wěn)定性條件，即

圖1 EMMS模型的物理圖譜和表達(dá)式[65]Fig.1 Physical concept and formulation of the EMMS model[65]

以上穩(wěn)定性條件和守恒方程共同構(gòu)成EMMS模型。在給定的宏尺度操作條件下時(shí)，X可以通過數(shù)值方法求解。然而，由于EMMS 模型最初用于描述流化床系統(tǒng)的整體流動(dòng)結(jié)構(gòu)，將其拓展到計(jì)算網(wǎng)格求解有效曳力時(shí)至少需要解決和考慮以下問題：（1）局部結(jié)構(gòu)遠(yuǎn)離平衡態(tài)，模型方程如何考慮顆粒加速度的影響；（2）原穩(wěn)定性條件在拓展到其他流域時(shí)應(yīng)該如何改進(jìn)；（3）曳力關(guān)聯(lián)式的形式及其計(jì)算方法。

1.1 改進(jìn)的基于團(tuán)聚物的EMMS曳力

EMMS模型基于能量平衡假設(shè)建立了團(tuán)聚物方程，但該方程難以適用于低通量、高加速度和順重力情形，這限制了模型在曳力計(jì)算中的應(yīng)用。為了改進(jìn)對顆粒團(tuán)聚動(dòng)力學(xué)的描述以提高EMMS曳力模型的普適性，Hu 等[30]通過量化團(tuán)聚物和稀相顆粒之間的質(zhì)量交換關(guān)系，提出了團(tuán)聚物動(dòng)態(tài)平衡方程，并以此替換原團(tuán)聚物方程對EMMS 模型進(jìn)行了改進(jìn)。如圖2(a)所示，通常情況下，流體計(jì)算網(wǎng)格中呈現(xiàn)稀密相共存的非均勻兩相結(jié)構(gòu)。在團(tuán)聚物存續(xù)時(shí)間內(nèi)，稀相和密相顆粒之間的質(zhì)量交換可以認(rèn)為由兩個(gè)可逆的過程組成。如圖2（b）所示，對于球形顆粒物，在迎風(fēng)面，稀相顆粒和團(tuán)聚物發(fā)生非彈性碰撞導(dǎo)致顆粒被聚團(tuán)捕捉；而在背風(fēng)面，團(tuán)聚物表面顆粒在隨機(jī)脈動(dòng)作用下進(jìn)入稀相。兩者滿足動(dòng)態(tài)平衡時(shí)可得

圖2 EMMS曳力模型中亞網(wǎng)格非均勻性（a）和稀密相質(zhì)量傳遞過程（b）Fig.2 Sub-grid heterogeneity(a)and solids mass transfer between the dilute and dense phases(b)in the EMMS drag model

其中，w是常數(shù)，本研究取1/12。密相顆粒溫度Θc的計(jì)算可參考文獻(xiàn)[30]。根據(jù)圖2(a)可以分別對稀相和密相建立守恒方程。連續(xù)性方程同原EMMS模型（圖1）一致，但由于局部結(jié)構(gòu)往往遠(yuǎn)離平衡態(tài)，顆粒的加速效應(yīng)常不能忽略，稀密兩相的力平衡方程應(yīng)改寫為

非均勻性因子定義為EMMS曳力系數(shù)和均勻曳力系數(shù)的比值，即

圖3 改進(jìn)EMMS曳力[（a）,(c),(e),(g)]和均勻曳力[（b）,(d),(f),(h)]對快速床提升管[(a)]、湍動(dòng)床[（c）]、變徑提升管[(e)）及下行床[(g)]的預(yù)測結(jié)果[30,51,73]Fig.3 Snapshot of solids concentration in the fast fluidized riser[(a)],turbulent bed[(c)],tapered riser[(e)]and co-current downer[(g)]by using the improved EMMS drag model[(a),(c),(e),(g)]and the homogeneous drag model[(b),(d),（f）,(h)][30,51,73]

1.2 具有伽利略不變性的EMMS曳力

傳統(tǒng)的EMMS曳力計(jì)算依賴系統(tǒng)宏觀操作參數(shù)（如Ug、Gs等），構(gòu)建滿足伽利略不變性的EMMS曳力模型可以避免當(dāng)操作條件改變時(shí)必須重構(gòu)Hd系數(shù)關(guān)聯(lián)式的問題，并可以適應(yīng)具有復(fù)雜幾何結(jié)構(gòu)的流化床反應(yīng)器的曳力計(jì)算。建立滿足伽利略不變性EMMS曳力模型的重要步驟是將方程中的絕對速度項(xiàng)轉(zhuǎn)化為相對速度項(xiàng)。根據(jù)滑移速度的定義，式（2）可以重組為

借助相間壓降平衡方程消去式（3）、式（4）中的Fi可得稀密相力平衡方程（k=c 代表密相；k=f 代表稀相）

為了獲得滿足伽利略不變性的穩(wěn)定性條件的等價(jià)表達(dá)，Hu 等[54]分析了不同操作條件下EMMS 模型的解集，發(fā)現(xiàn)當(dāng)氣固物性不變時(shí)，密相空隙率可以近似表達(dá)為平均空隙率的函數(shù)。對于典型A類顆粒（dp=60 μm,ρp=1490 kg/m3）有

質(zhì)量守恒方程（氣相、固相和固含率）、模型基本假設(shè)（εf→εmax,af→0）及式（7）～式（9）共同構(gòu)成滿足伽利略不變性的EMMS 曳力模型?？梢姡鄬υ瓉淼姆蔷€性優(yōu)化求解，上述模型的求解復(fù)雜度大大下降，可以容易地構(gòu)建Hd關(guān)于（ε，Re）的關(guān)聯(lián)式或在模擬過程中直接對模型進(jìn)行在線求解。圖4(a)展示了采用該模型計(jì)算得到的非均勻因子，其隨空隙率和Reynolds 數(shù)的變化符合物理規(guī)律[54]。圖4(b)進(jìn)一步比較了本模型、改進(jìn)EMMS 模型和實(shí)驗(yàn)數(shù)據(jù)的對比。可見式（9）的簡化并未降低模型的準(zhǔn)確性，但模型的普適性有了很大提高。

圖4 滿足伽利略不變性的EMMS曳力模型[54]Fig.4 EMMS drag model following Galilean invariance[54]

1.3 改進(jìn)的基于氣泡的EMMS曳力

EMMS 曳力模型也被拓展到氣固鼓泡流態(tài)化系統(tǒng)[43-44,46,49-50]。但由于經(jīng)驗(yàn)性的氣泡關(guān)聯(lián)式往往基于特定的實(shí)驗(yàn)條件獲得而難以被推廣到其他情況，因而Liu 等[45]提出了采用改進(jìn)穩(wěn)定性條件來約束氣泡尺寸的方法。

如圖5 所示，氣固鼓泡流態(tài)化系統(tǒng)可以分解為氣泡相和乳化相，共采用9 個(gè)結(jié)構(gòu)參數(shù)X={Ub,fb,εb,db,ab,Uge,Upe,εe,ape}來定義。但為了簡化計(jì)算，一般假設(shè)氣泡中不含顆粒，即εb=1.0。分別以顆粒和氣泡為研究對象可以建立一系列連續(xù)性方程和力平衡方程。與上述基于團(tuán)聚物的EMMS 模型不同，EMMS鼓泡模型一開始就考慮了顆粒和氣泡的加速效應(yīng)，且在穩(wěn)定性條件的表達(dá)上也有很大的變化。這是因?yàn)闅馀菰陂L大過程中會(huì)因?yàn)橄嘟缑娴恼龎毫Χa(chǎn)生能量耗散，因此在計(jì)算單位質(zhì)量床層的懸浮輸送能耗（Nst）時(shí)，除了曳力做功外，還應(yīng)額外考慮氣泡克服乳化相壓力膨脹所做的體積功。系統(tǒng)的穩(wěn)定性條件定義為相對懸浮輸送能耗趨于最小，其具體形式見圖5。改進(jìn)后的穩(wěn)定性條件中包含了氣泡直徑的影響，這為氣泡直徑的取值提供了自然約束，避免了經(jīng)驗(yàn)關(guān)聯(lián)式的引入。

圖5 EMMS鼓泡曳力模型的物理圖解和公式Fig.5 Physical diagram and formulation of the bubbling EMMS drag model

在鼓泡床中，非均勻性因子可以表達(dá)為

鼓泡EMMS 曳力的計(jì)算及與CFD 耦合的方式同1.1 節(jié)類似，此處不再贅述。如圖6 所示，相對均勻曳力模型，EMMS 鼓泡曳力模型可以更好地預(yù)測反應(yīng)器的非均勻結(jié)構(gòu)，與實(shí)驗(yàn)數(shù)據(jù)呈現(xiàn)更好的一致性，這也說明了模型的合理性。關(guān)于模型方程的推導(dǎo)和實(shí)驗(yàn)驗(yàn)證的詳細(xì)信息可以見文獻(xiàn)[45]。

圖6 采用EMMS鼓泡曳力和均勻曳力對A類顆粒鼓泡床的模擬結(jié)果對比[45]Fig.6 Comparison between the bubbling EMMS model and homogeneous model in predicting the gas-solid bubbling beds of Geldart A particles[45]

2 介尺度結(jié)構(gòu)時(shí)空動(dòng)態(tài)演化的群平衡建模

2.1 基于團(tuán)聚物動(dòng)態(tài)演化的群平衡建模

在傳統(tǒng)的EMMS 曳力模型中，所取團(tuán)聚物直徑可以認(rèn)為是當(dāng)前操作條件（Ug,Up,ε）下的最概然值或統(tǒng)計(jì)均值，而團(tuán)聚物的動(dòng)態(tài)特征及尺寸分布（CSD）信息未被考慮。為了定量描述團(tuán)聚物的時(shí)空動(dòng)態(tài)分布及其對網(wǎng)格內(nèi)曳力的影響，Hu 等[62]提出了基于團(tuán)聚物的群平衡方程，并進(jìn)一步建立了CFD、PBM及EMMS模型的耦合計(jì)算框架。

以團(tuán)聚物直徑（L）為內(nèi)坐標(biāo)，團(tuán)聚物尺寸演化的群平衡方程可以表達(dá)為

式中，n(L,x,t)表示尺寸為L的團(tuán)聚物在t時(shí)刻和x位置的數(shù)密度；B(L,x,t)和D(L,x,t)分別表示由于聚并和破碎效應(yīng)導(dǎo)致的團(tuán)聚物生成項(xiàng)和消亡項(xiàng)，一般表達(dá)為

由于流化床中顆粒湍流及團(tuán)聚物聚并、破碎的機(jī)制尚不明確，因此初步研究主要聚焦于團(tuán)聚物和稀相顆粒的碰撞和離析等連續(xù)演化過程。通過分析單位時(shí)間顆粒和團(tuán)聚物的非彈性碰撞及稀相流體對聚團(tuán)表面顆粒的侵蝕效應(yīng)，Hu 等[62]建立了團(tuán)聚物的連續(xù)生長方程，代入PBM可得

式(14)理論上更適用于床層較稀的工況。通過離散法或矩方法可以對上述方程進(jìn)行求解，進(jìn)而可以獲得團(tuán)聚物的數(shù)密度函數(shù)及網(wǎng)格內(nèi)的平均直徑分布。

前述信息可以用來封閉EMMS曳力模型中的團(tuán)聚物直徑項(xiàng)。圖7 展示了CFD-PBM-EMMS 的耦合計(jì)算框架：首先利用CFD 模型計(jì)算得到的t時(shí)刻所有網(wǎng)格內(nèi)的相速度、空隙率、湍流特征量并傳遞給PBM；求解PBM 可以得到t時(shí)刻網(wǎng)格內(nèi)團(tuán)聚物的尺寸分布（CSD）；將以上CFD 和PBM 的求解信息傳遞給EMMS 曳力模型，進(jìn)而獲得粗網(wǎng)格的有效曳力系數(shù)，并提高CFD的計(jì)算精度。

圖7 CFD、PBM和EMMS模型的耦合Fig.7 Integration of the CFD,PBM and EMMS models

耦合模型可以有效避免實(shí)際團(tuán)聚物直徑和傳統(tǒng)EMMS預(yù)測值不一致的現(xiàn)象，具有更好的自洽性。更重要的是，該算法可以直接計(jì)算反應(yīng)器內(nèi)不同位置處團(tuán)聚物的尺寸分布隨時(shí)間的演化規(guī)律，這對反應(yīng)器的設(shè)計(jì)和優(yōu)化無疑具有重要的意義。圖8(a)展示了基于PBM 修正的EMMS 曳力校正系數(shù)，可見該模型合理地預(yù)測了團(tuán)聚物尺寸這一關(guān)鍵參數(shù)對Hd的影響。通過圖8(b)可見，相對均勻曳力和傳統(tǒng)EMMS 曳力，本模型可以更好地預(yù)測流化床內(nèi)的非均勻分布。另外，群平衡模型還可以對團(tuán)聚物的尺寸分布進(jìn)行直接預(yù)測，如圖8(c)所示。

圖8 CFD-PBM-EMMS耦合模型在快速流態(tài)化中的應(yīng)用[62]Fig.8 Implement of the CFD-PBM-EMMS model in the fast fluidized beds[62]

2.2 基于氣泡動(dòng)態(tài)演化的群平衡建模

在氣固鼓泡流態(tài)化中，氣泡自分布板上部生成之后，往往經(jīng)歷復(fù)雜的聚并、破碎等動(dòng)態(tài)演化過程[74-75]。很多研究采用實(shí)驗(yàn)和經(jīng)驗(yàn)建模的方式對氣泡的尺寸分布進(jìn)行研究和建模，不僅難以適應(yīng)各種復(fù)雜操作條件，也無法對氣泡的動(dòng)態(tài)時(shí)空分布進(jìn)行準(zhǔn)確計(jì)算。采用這些經(jīng)驗(yàn)方程進(jìn)行封閉也在很大程度上影響了曳力模型的準(zhǔn)確性和合理性。

根據(jù)Hu等[63]的研究，氣泡的聚并主要由大小氣泡的上升速度差異和氣泡的尾流加速引起的，因此根據(jù)兩者的物理過程可分別建立其聚并核函數(shù)。大小為di和dj的兩個(gè)氣泡的聚并頻率可以表達(dá)為

式中，λr和λw分別表示由于氣泡追趕和尾流效應(yīng)導(dǎo)致的聚并效率[63]；θw表示尾流的平均張角；urel,w表示由于尾流誘導(dǎo)產(chǎn)生的碰撞速度，實(shí)驗(yàn)發(fā)現(xiàn)該附加速度和碰撞發(fā)生時(shí)的氣泡上升速度大小相當(dāng)。氣泡的破碎機(jī)制則更加復(fù)雜，一般認(rèn)為是頂部顆?？焖傧侣鋵?dǎo)致的氣泡切割和分裂。通過對實(shí)驗(yàn)數(shù)據(jù)回歸發(fā)現(xiàn)氣泡的破碎頻率和其直徑的2.9 次冪成正比。因此，直徑為di的氣泡破碎頻率可以通過式(16)關(guān)聯(lián)獲得

而破裂導(dǎo)致的子氣泡概率分布函數(shù)則采用實(shí)驗(yàn)統(tǒng)計(jì)的擬合函數(shù)。體積為V?的氣泡破碎生成的子氣泡近似滿足以下高斯分布

將式（15）～式（17）代入式（11）～式（13）可以得到氣泡的群平衡方程。采用和圖7一致的建模思路建立鼓泡床的CFD-PBM-EMMS 模型。如圖9(a)所示，借助PBM 封閉的鼓泡EMMS 模型合理地考慮了氣泡直徑對非均勻曳力的影響。在氣泡直徑趨于0時(shí)，曳力自然趨于均勻曳力。圖9(b)是模型預(yù)測結(jié)果和實(shí)驗(yàn)的對比，可見該模型可以較好地預(yù)測氣泡尺寸沿床層高度的分布規(guī)律。

圖9 CFD-PBM-EMMS耦合模型在鼓泡流態(tài)化中的應(yīng)用[63]Fig.9 Implementation of the CFD-PBM-EMMS model in the bubbling fluidized beds[63]

3 基于EMMS的全流域操作相圖和全回路穩(wěn)態(tài)建模

3.1 穩(wěn)態(tài)EMMS模型的擴(kuò)展和相圖繪制

研究表明氣固流態(tài)化和物質(zhì)相變具有相似的多狀態(tài)特征，隨操作條件改變呈現(xiàn)不同的非均勻分布和連續(xù)的流域轉(zhuǎn)變[76-78]。量化流化床的宏尺度穩(wěn)態(tài)動(dòng)力學(xué)不僅可以快速了解反應(yīng)器的操作狀況，亦可以作為初始分布加速工業(yè)級(jí)反應(yīng)的數(shù)值模擬過程。Hu 等[67]發(fā)現(xiàn)在提升管底部團(tuán)聚物具有強(qiáng)烈的碰撞、聚并作用，而在頂部團(tuán)聚物組織松散，稀密相差異很小。在此基礎(chǔ)上，通過引入顆粒加速度建立了EMMS 軸向模型，首次實(shí)現(xiàn)了對提升管反應(yīng)器中軸向S 形非均勻分布的EMMS 預(yù)測[67]。根據(jù)軸向模型獲得各截面處的平均空隙率和氣固速度信息，Hu等[30]進(jìn)一步改進(jìn)和完善了徑向EMMS 模型，并通過函數(shù)逼近的方法實(shí)現(xiàn)了對提升管中徑向環(huán)核結(jié)構(gòu)的數(shù)學(xué)預(yù)測。

圖10 展示了徑向EMMS 模型的物理框架。自上而下，在截面尺度，各物理量必須滿足截面穩(wěn)定性條件和邊界約束的守恒方程；在單元介尺度，各徑向微元體內(nèi)參數(shù)服從于穩(wěn)定性條件和質(zhì)量動(dòng)量守恒方程；而在顆粒微尺度，氣固相互作用在稀密兩相內(nèi)呈現(xiàn)均勻結(jié)構(gòu)，可以用均勻曳力計(jì)算。與提升管反應(yīng)器類似，氣固下行床中也存在顆粒團(tuán)聚現(xiàn)象。因此，上述軸徑向建模思路也被應(yīng)用到氣固并/逆流下行床的計(jì)算。將氣固并/逆流下行床進(jìn)行結(jié)構(gòu)分解，分別建立各相的質(zhì)量和動(dòng)量守恒方程以及不同發(fā)展階段的穩(wěn)定性條件，最終可建立EMMS 并/逆流下行床數(shù)學(xué)模型[68-69]。至于鼓泡流態(tài)化系統(tǒng)，氣泡的形成和演化代替了提升管中顆粒團(tuán)聚物的聚并和破碎而成為系統(tǒng)的主要非均勻結(jié)構(gòu)特征。假設(shè)乳化相的顆粒加速度為0，則圖5所示方程可以直接用于鼓泡床穩(wěn)態(tài)動(dòng)力學(xué)的求解。模型無須經(jīng)驗(yàn)氣泡直徑的輔助就可對鼓泡床動(dòng)力學(xué)進(jìn)行預(yù)測[45]。 更多研究細(xì)節(jié)和實(shí)驗(yàn)驗(yàn)證可參考文獻(xiàn)[30,45,67-70]。

圖10 提升管徑向非均勻分布的EMMS穩(wěn)態(tài)建模[64-65]Fig.10 Modelling of radial hydrodynamics of risers by using the EMMS model[64-65]

EMMS模型在二維空間以及不同流域的擴(kuò)展使廣義流態(tài)化相圖的繪制成為可能。根據(jù)Kwauk 廣義流態(tài)化相圖的繪制思路[79]，遍歷不同操作范圍的氣體速度（=Ug/Ut）和顆粒速度(=Up/Ut)，判斷其所處流域，進(jìn)而采用相應(yīng)的EMMS 模型或擴(kuò)展模型計(jì)算床層的平均空隙率，最終可得如圖11所示的操作相圖。值得說明的是，該相圖給出了典型流域轉(zhuǎn)變時(shí)的特征速度曲線，如噎塞曲線、液泛曲線和湍動(dòng)流態(tài)化速度的參考值。噎塞的預(yù)測作為EMMS模型的重要特征已經(jīng)被廣泛討論，這里不再贅述。而湍動(dòng)流域的確定一直存在廣泛爭議，不同研究者對湍動(dòng)氣速的確定標(biāo)準(zhǔn)持有不同的觀點(diǎn)。本研究根據(jù)基于團(tuán)聚物的原始EMMS模型和基于氣泡的鼓泡EMMS 模型聯(lián)立求解，所獲得的二者的交叉點(diǎn)被認(rèn)為是連續(xù)相從密相轉(zhuǎn)變?yōu)橄∠嗟南噢D(zhuǎn)變點(diǎn)，并將其所對應(yīng)的流體速度定義為湍動(dòng)流態(tài)化速度，如圖中黑色虛線Uc所示。而液泛作為另外一個(gè)重要的流域轉(zhuǎn)折特征，一般發(fā)生在當(dāng)逆流床的氣速或顆粒通量過大時(shí)，顆粒被氣流帶出從而中斷系統(tǒng)操作的狀態(tài)。通過考察EMMS 逆流下行床模型解的特點(diǎn)，本研究將空隙率對氣速二階導(dǎo)的最大值點(diǎn)設(shè)為液泛發(fā)生速度，而EMMS 無解點(diǎn)設(shè)為液泛終止速度[65]，并通過和實(shí)驗(yàn)對比驗(yàn)證了這種定義的合理性。

圖11 Geldart B類顆粒的EMMS廣義流態(tài)化相圖[69]Fig.11 EMMS-based generalized fluidization phase diagrams for Geldart B particles[69]

3.2 循環(huán)流化床全回路穩(wěn)態(tài)建模及軟件開發(fā)

基于以上各擴(kuò)展模型和操作相圖，本課題組進(jìn)一步完善了基于EMMS方法的循環(huán)流化床全回路穩(wěn)態(tài)計(jì)算模型[70,72]。

圖12(a)展示了循環(huán)流化床全回路建模的基本流程：（1）首先根據(jù)系統(tǒng)的邊界和接口情況將其劃分為不同的循環(huán)通路和相應(yīng)的子反應(yīng)器組合；（2）遍歷各循環(huán)通路內(nèi)的顆粒循環(huán)量等參數(shù)，并根據(jù)所處的氣固流速選擇合適的計(jì)算模型，如果操作處于廣義流態(tài)化狀態(tài)，則選取EMMS 擴(kuò)展模型計(jì)算或通過相圖查閱該單元內(nèi)的動(dòng)力學(xué)參數(shù)，其他單元模塊如旋風(fēng)、出口管等則使用經(jīng)驗(yàn)關(guān)聯(lián)式求解；（3）根據(jù)系統(tǒng)的壓降平衡和物料平衡迭代計(jì)算直到收斂。圖12(b)是一套多產(chǎn)異構(gòu)烷烴的催化裂化過程冷態(tài)實(shí)驗(yàn)裝置。圖12(c)是對該反應(yīng)器進(jìn)行全回路計(jì)算的結(jié)果，可見模型可以合理預(yù)測反應(yīng)器內(nèi)流動(dòng)參數(shù)的宏觀分布。

圖12 基于EMMS的循環(huán)流化床統(tǒng)一建模[72]Fig.12 EMMS-based general method for global hydrodynamic calculation of typical CFB systems[72]

為實(shí)現(xiàn)流化床反應(yīng)器整體穩(wěn)態(tài)動(dòng)力學(xué)的快速計(jì)算，EMMS 全回路建模方法還被封裝成軟件Virtual Fluidization（2015SRBJ0279）。軟件采用模塊化設(shè)計(jì)，可以根據(jù)實(shí)際的工業(yè)裝置自定義搭建、更改和優(yōu)化反應(yīng)器設(shè)計(jì)。軟件系統(tǒng)通過調(diào)用全循環(huán)計(jì)算代碼可以在幾秒到幾分鐘之內(nèi)完成穩(wěn)態(tài)動(dòng)力學(xué)的計(jì)算和顯示。如圖13所示，系統(tǒng)合理預(yù)測了不同操作條件下反應(yīng)器內(nèi)的各模塊的顆粒濃度非均勻分布。通過導(dǎo)出整體分布數(shù)據(jù)作為后續(xù)CFD 計(jì)算的初場，可以提高CFD模擬的效率和速度。

圖13 采用Virtual Fluidization軟件對循環(huán)流化床系統(tǒng)的穩(wěn)態(tài)動(dòng)力學(xué)進(jìn)行模擬[72]Fig.13 Simulation of steady-state hydrodynamics in a CFB system by using the software Virtual Fluidization[72]

4 總結(jié)與展望

能量最小多尺度方法為多相非均勻系統(tǒng)的數(shù)學(xué)建模提供了一種可行的理論框架。采用該方法能夠計(jì)算局部介尺度結(jié)構(gòu)對相間曳力等本構(gòu)關(guān)系的影響。但經(jīng)典EMMS曳力模型所采用的團(tuán)聚物方程在極端操作條件下存在求解不適定的問題，而通過量化稀密相間固相質(zhì)量傳遞可以建立更通用的顆粒團(tuán)聚動(dòng)態(tài)平衡方程，基于該方程的改進(jìn)EMMS曳力可以更合理地計(jì)算快速床、湍動(dòng)床和下行床的非均勻分布。通過改進(jìn)EMMS鼓泡模型的穩(wěn)定性條件，無需氣泡經(jīng)驗(yàn)關(guān)聯(lián)式輔助即可實(shí)現(xiàn)對亞網(wǎng)格鼓泡曳力的直接預(yù)測。

通過考慮氣泡和團(tuán)聚物的聚并、破碎演化動(dòng)力學(xué)，群平衡模型實(shí)現(xiàn)了對該類介尺度結(jié)構(gòu)時(shí)空動(dòng)態(tài)演化行為的定量表征。通過群平衡方程在線封閉EMMS 模型中的氣泡/團(tuán)聚物直徑項(xiàng)可以更準(zhǔn)確地計(jì)算局部結(jié)構(gòu)曳力系數(shù)，從而提高CFD 模擬的可靠性。

通過對EMMS 模型在軸徑向二維空間、鼓泡床和下行床進(jìn)行拓展，建立了適應(yīng)不同操作模式和操作條件以及流域的反應(yīng)器宏尺度穩(wěn)態(tài)模型，并繪制基于EMMS 的氣固廣義流態(tài)化相圖，反映不同操作模式下床層空隙率、噎塞和液泛氣速等特征變量隨操作條件的變化。原EMMS及其拓展模型被應(yīng)用到循環(huán)流化床系統(tǒng)的全回路建模，實(shí)現(xiàn)了對復(fù)雜回路氣固系統(tǒng)的準(zhǔn)實(shí)時(shí)穩(wěn)態(tài)動(dòng)力學(xué)計(jì)算。

未來的工作可聚焦介尺度結(jié)構(gòu)的動(dòng)態(tài)演化機(jī)制，通過實(shí)驗(yàn)和模擬等手段探索局部非均勻結(jié)構(gòu)的時(shí)間演化特征并量化這種動(dòng)態(tài)結(jié)構(gòu)對曳力、固相應(yīng)力及湍流耗散等本構(gòu)關(guān)系的影響。另外，復(fù)雜反應(yīng)器的快速準(zhǔn)確建模對工業(yè)應(yīng)用具有重要意義，如何結(jié)合介尺度模型、高效數(shù)值模擬方法以及神經(jīng)網(wǎng)絡(luò)方法建立完善的多尺度計(jì)算模式也是值得努力的方向。

符 號(hào) 說 明

a——加速度，m/s2

a(L,λ)——直徑為L和λ的單體聚并頻率,m3/s

B——生成源項(xiàng)，m-4/s

CD——單顆粒曳力系數(shù)

D——消亡源項(xiàng)，m-4/s

dcl——團(tuán)聚物直徑，m

dp——顆粒直徑，m

Fc——單位體積密相顆粒所受曳力，kg/(m2·s2)

Ff——單位體積稀相顆粒所受曳力，kg/(m2·s2)

Fi——單位體積床層團(tuán)聚物所受稀相曳力，kg/(m2·s2)

f——密相體積分?jǐn)?shù)

Gs——顆粒循環(huán)率，kg/(m2·s)

g——重力加速度，m/s2

g(λ)——破碎頻率，s-1

Hd——非均勻性因子

L——內(nèi)基變量，此處為單體直徑，m

Nst——?dú)怏w流過單位質(zhì)量顆粒床層的懸浮輸送能耗，J/(kg·s)

NT——?dú)怏w流過單位質(zhì)量顆粒床層的總能耗，J/(kg·s)

n——基于長度的數(shù)密度函數(shù),m-4

p——破碎元

Re——Reynolds數(shù)

U——表觀速度，m/s

Umf——起始流態(tài)化速度，m/s

Ut——顆粒終端速度，m/s

u,u——速度，m/s

Wst——?dú)怏w流過單位體積床層的懸浮輸送能耗，J/(m3·s)

βE,βhomo——分別為微元體的有效曳力系數(shù)和均勻曳力系數(shù)，kg/(m3·s)

β(L|λ)——碎片分布函數(shù)

ε——體積分?jǐn)?shù)

εmax——團(tuán)聚物存在的最大空隙率

λr,λw——?dú)馀葑汾s和尾流導(dǎo)致的碰撞效率

ρ——密度，kg/m3

Θc——密相顆粒溫度，m2/s2

θw——尾流的平均張角

ω——概率因子

下角標(biāo)

b——?dú)馀?/p>

c——密相

e——乳化相

f——稀相

g——?dú)怏w

p——顆粒

s——滑移