氣固流化系統多尺度跨流域EMMS建模

胡善偉，劉新華

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

引 言

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

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

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

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

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

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

1 EMMS曳力模型

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

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

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

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

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

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

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

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

圖3 改進EMMS曳力[（a）,(c),(e),(g)]和均勻曳力[（b）,(d),(f),(h)]對快速床提升管[(a)]、湍動床[（c）]、變徑提升管[(e)）及下行床[(g)]的預測結果[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曳力

傳統的EMMS曳力計算依賴系統宏觀操作參數（如Ug、Gs等），構建滿足伽利略不變性的EMMS曳力模型可以避免當操作條件改變時必須重構Hd系數關聯式的問題，并可以適應具有復雜幾何結構的流化床反應器的曳力計算。建立滿足伽利略不變性EMMS曳力模型的重要步驟是將方程中的絕對速度項轉化為相對速度項。根據滑移速度的定義，式（2）可以重組為

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

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

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

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

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

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

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

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

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

鼓泡EMMS 曳力的計算及與CFD 耦合的方式同1.1 節類似，此處不再贅述。如圖6 所示，相對均勻曳力模型，EMMS 鼓泡曳力模型可以更好地預測反應器的非均勻結構，與實驗數據呈現更好的一致性，這也說明了模型的合理性。關于模型方程的推導和實驗驗證的詳細信息可以見文獻[45]。

圖6 采用EMMS鼓泡曳力和均勻曳力對A類顆粒鼓泡床的模擬結果對比[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 介尺度結構時空動態演化的群平衡建模

2.1 基于團聚物動態演化的群平衡建模

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

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

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

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

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

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

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

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

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

2.2 基于氣泡動態演化的群平衡建模

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

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

式中，λr和λw分別表示由于氣泡追趕和尾流效應導致的聚并效率[63]；θw表示尾流的平均張角；urel,w表示由于尾流誘導產生的碰撞速度，實驗發現該附加速度和碰撞發生時的氣泡上升速度大小相當。氣泡的破碎機制則更加復雜，一般認為是頂部顆粒快速下落導致的氣泡切割和分裂。通過對實驗數據回歸發現氣泡的破碎頻率和其直徑的2.9 次冪成正比。因此，直徑為di的氣泡破碎頻率可以通過式(16)關聯獲得

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

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

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

3 基于EMMS的全流域操作相圖和全回路穩態建模

3.1 穩態EMMS模型的擴展和相圖繪制

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

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

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

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

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

3.2 循環流化床全回路穩態建模及軟件開發

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

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

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

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

圖13 采用Virtual Fluidization軟件對循環流化床系統的穩態動力學進行模擬[72]Fig.13 Simulation of steady-state hydrodynamics in a CFB system by using the software Virtual Fluidization[72]

4 總結與展望

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

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

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

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

符 號 說 明

a——加速度，m/s2

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

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

CD——單顆粒曳力系數

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

dcl——團聚物直徑，m

dp——顆粒直徑，m

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

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

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

f——密相體積分數

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

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

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

Hd——非均勻性因子

L——內基變量，此處為單體直徑，m

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

NT——氣體流過單位質量顆粒床層的總能耗，J/(kg·s)

n——基于長度的數密度函數,m-4

p——破碎元

Re——Reynolds數

U——表觀速度，m/s

Umf——起始流態化速度，m/s

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

u,u——速度，m/s

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

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

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

ε——體積分數

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

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

ρ——密度，kg/m3

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

θw——尾流的平均張角

ω——概率因子

下角標

b——氣泡

c——密相

e——乳化相

f——稀相

g——氣體

p——顆粒

s——滑移