氣固流化系統(tǒng)多尺度跨流域EMMS建模

胡善偉，劉新華

（中國科學院過程工程研究所，多相復雜系統(tǒng)國家重點實驗室，北京 100190）

引 言

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

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

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

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

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

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

1 EMMS曳力模型

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

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

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

1.1 改進的基于團聚物的EMMS曳力

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

圖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的計算可參考文獻[30]。根據(jù)圖2(a)可以分別對稀相和密相建立守恒方程。連續(xù)性方程同原EMMS模型（圖1）一致，但由于局部結(jié)構(gòu)往往遠離平衡態(tài)，顆粒的加速效應常不能忽略，稀密兩相的力平衡方程應改寫為

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

圖3 改進EMMS曳力[（a）,(c),(e),(g)]和均勻曳力[（b）,(d),(f),(h)]對快速床提升管[(a)]、湍動床[（c）]、變徑提升管[(e)）及下行床[(g)]的預測結(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曳力計算依賴系統(tǒng)宏觀操作參數(shù)（如Ug、Gs等），構(gòu)建滿足伽利略不變性的EMMS曳力模型可以避免當操作條件改變時必須重構(gòu)Hd系數(shù)關聯(lián)式的問題，并可以適應具有復雜幾何結(jié)構(gòu)的流化床反應器的曳力計算。建立滿足伽利略不變性EMMS曳力模型的重要步驟是將方程中的絕對速度項轉(zhuǎn)化為相對速度項。根據(jù)滑移速度的定義，式（2）可以重組為

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

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

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

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

1.3 改進的基于氣泡的EMMS曳力

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

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

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

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

鼓泡EMMS 曳力的計算及與CFD 耦合的方式同1.1 節(jié)類似，此處不再贅述。如圖6 所示，相對均勻曳力模型，EMMS 鼓泡曳力模型可以更好地預測反應器的非均勻結(jié)構(gòu)，與實驗數(shù)據(jù)呈現(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)時空動態(tài)演化的群平衡建模

2.1 基于團聚物動態(tài)演化的群平衡建模

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 CFD-PBM-EMMS耦合模型在鼓泡流態(tài)化中的應用[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模型的擴展和相圖繪制

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

圖10 展示了徑向EMMS 模型的物理框架。自上而下，在截面尺度，各物理量必須滿足截面穩(wěn)定性條件和邊界約束的守恒方程；在單元介尺度，各徑向微元體內(nèi)參數(shù)服從于穩(wěn)定性條件和質(zhì)量動量守恒方程；而在顆粒微尺度，氣固相互作用在稀密兩相內(nèi)呈現(xiàn)均勻結(jié)構(gòu)，可以用均勻曳力計算。與提升管反應器類似，氣固下行床中也存在顆粒團聚現(xiàn)象。因此，上述軸徑向建模思路也被應用到氣固并/逆流下行床的計算。將氣固并/逆流下行床進行結(jié)構(gòu)分解，分別建立各相的質(zhì)量和動量守恒方程以及不同發(fā)展階段的穩(wěn)定性條件，最終可建立EMMS 并/逆流下行床數(shù)學模型[68-69]。至于鼓泡流態(tài)化系統(tǒng)，氣泡的形成和演化代替了提升管中顆粒團聚物的聚并和破碎而成為系統(tǒng)的主要非均勻結(jié)構(gòu)特征。假設乳化相的顆粒加速度為0，則圖5所示方程可以直接用于鼓泡床穩(wěn)態(tài)動力學的求解。模型無須經(jīng)驗氣泡直徑的輔助就可對鼓泡床動力學進行預測[45]。 更多研究細節(jié)和實驗驗證可參考文獻[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模型在二維空間以及不同流域的擴展使廣義流態(tài)化相圖的繪制成為可能。根據(jù)Kwauk 廣義流態(tài)化相圖的繪制思路[79]，遍歷不同操作范圍的氣體速度（=Ug/Ut）和顆粒速度(=Up/Ut)，判斷其所處流域，進而采用相應的EMMS 模型或擴展模型計算床層的平均空隙率，最終可得如圖11所示的操作相圖。值得說明的是，該相圖給出了典型流域轉(zhuǎn)變時的特征速度曲線，如噎塞曲線、液泛曲線和湍動流態(tài)化速度的參考值。噎塞的預測作為EMMS模型的重要特征已經(jīng)被廣泛討論，這里不再贅述。而湍動流域的確定一直存在廣泛爭議，不同研究者對湍動氣速的確定標準持有不同的觀點。本研究根據(jù)基于團聚物的原始EMMS模型和基于氣泡的鼓泡EMMS 模型聯(lián)立求解，所獲得的二者的交叉點被認為是連續(xù)相從密相轉(zhuǎn)變?yōu)橄∠嗟南噢D(zhuǎn)變點，并將其所對應的流體速度定義為湍動流態(tài)化速度，如圖中黑色虛線Uc所示。而液泛作為另外一個重要的流域轉(zhuǎn)折特征，一般發(fā)生在當逆流床的氣速或顆粒通量過大時，顆粒被氣流帶出從而中斷系統(tǒng)操作的狀態(tài)。通過考察EMMS 逆流下行床模型解的特點，本研究將空隙率對氣速二階導的最大值點設為液泛發(fā)生速度，而EMMS 無解點設為液泛終止速度[65]，并通過和實驗對比驗證了這種定義的合理性。

圖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ā)

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

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

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

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

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

4 總結(jié)與展望

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

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

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

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

符 號 說 明

a——加速度，m/s2

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

B——生成源項，m-4/s

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

D——消亡源項，m-4/s

dcl——團聚物直徑，m

dp——顆粒直徑，m

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

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

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

f——密相體積分數(shù)

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

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

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

Hd——非均勻性因子

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

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

NT——氣體流過單位質(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——氣體流過單位體積床層的懸浮輸送能耗，J/(m3·s)

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

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

ε——體積分數(shù)

εmax——團聚物存在的最大空隙率

λr,λw——氣泡追趕和尾流導致的碰撞效率

ρ——密度，kg/m3

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

θw——尾流的平均張角

ω——概率因子

下角標

b——氣泡

c——密相

e——乳化相

f——稀相

g——氣體

p——顆粒

s——滑移