基于介尺度穩定性條件的多相流曳力與群體平衡模型

管小平，楊寧

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

引 言

多相流廣泛存在于化工的單元操作，比如輸送、精餾、結晶、干燥、反應等，是傳遞-反應的重要媒介?；ざ嘞嗔鞒尸F非線性、多尺度、多流型的特征，直接影響著傳質、傳熱和反應速率，是工藝過程放大困難的根源。比如，對于氣固兩相流，隨著操作條件的變化，出現固定床、鼓泡床、湍動床、快速床以及輸送床流型，床內涉及小尺度單顆粒、介尺度氣泡或者顆粒團、大尺度全床的流動行為。因此，復雜多尺度結構的解析及形成機理是科學認識化工多相流的難點。

近些年，計算流體力學(computational fluid dynamics,CFD)成為認識多相流的重要手段，包括直接數值模擬(DNS)、歐拉-拉格朗日方法以及歐拉-歐拉方法[1]。歐拉-歐拉方法的模型是對單相流模型采用平均化方法得到的連續介質多流體模型[2]。模型中含有額外的非封閉項，用于描述相內或者相間微觀相互作用對宏觀流動的影響。常見的封閉子模型包括相間作用力模型、聚并破碎核函數以及湍流模型。相間作用力模型用于封閉由于相界面上黏性應力和壓力分布不均產生的平均作用力，包括曳力、升力、虛擬質量力、湍流擴散力以及壁面潤滑力等[3]。聚并破碎核函數針對的是分散相內部的相互作用，描述分散相聚并破碎過程的動力學，包括聚并速率模型、破碎速率模型以及子氣泡分布模型[4,5]。湍流模型用于封閉速度脈動產生的額外應力，包括大渦模擬和雷諾平均湍流模型。盡管多流體模型框架經過嚴格的數學推導，解決了宏觀輸運尺度的動量傳遞問題，但將多相流的物理復雜性轉嫁到封閉子模型中，關鍵科學問題并沒有解決。

從介科學角度發展介尺度物理模型，解析多相流非均勻結構演化的控制機制，是認識多相流復雜性的新途徑。介尺度模型提出于20 世紀80 年代末，針對氣固流態化體系，認為在物理上流態化系統存在稀、密相兩種不同的運動機制。兩種機制的競爭協調構成了系統的介尺度穩定性條件，是系統介尺度結構形成及演化的推動力[6]。21 世紀初，該方法拓展至氣液鼓泡體系，認為體系存在有利于形成大氣泡和小氣泡的兩種形態機制，兩種機制的競爭協調決定了氣泡形態的演化[7]。進一步將介尺度模型與CFD 相結合，為發展多相流模型提供了新思路。

通過介尺度理論模型，導出兩類介尺度封閉模型：介尺度曳力模型和介尺度群體平衡模型，本質上是通過介尺度穩定性條件決定系統的介尺度結構和封閉模型。將該封閉模型與多流體模型耦合，可用于多相流裝備的模擬計算，實現了介尺度模型從理論到實用的跨越。本文將從能量最小多尺度方法、介尺度曳力模型、介尺度群體平衡模型以及模型應用四個方面進行綜述，討論模型的發展歷程以及模型在實際體系的應用，最后對該方向的發展進行展望。

1 能量最小多尺度方法

圖1 EMMS方法的總體思路[9]Fig.1 General framework of EMMS approach[9]

圖2 氣液體系的介尺度理論模型[7]Fig.2 A mesoscale model for gas-liquid bubbly flow[7]

介尺度理論模型結果表明[12-14]，對于空氣-水體系鼓泡塔，隨著表觀氣速增加，氣含率增大，但在表觀氣速0.101 m/s 時，氣含率發生跳躍，由0.275 降低至0.248，出現動力學系統中的“分岔”現象，反映了系統本征穩定性的變化規律，流域由此過渡到湍流鼓泡區。數學上，在過渡點附近Nst存在兩個極小值，而Nst=min 在兩個極小值間發生了轉化，導致了系統結構參數的跳躍，物理上表現為系統在兩個穩定狀態之間發生轉變，宏觀流域發生了過渡?？梢?，理論模型能很好地闡明宏觀流域過渡的介尺度物理機制[7,12,15]。

Ruzicka 等[16-17]的實驗表明，黏度或者表面張力對氣液流動穩定性呈雙重影響：在低黏度或者表面張力時，增加黏度或者表面張力有利于流動穩定，流域過渡延遲；而在高黏度或者表面張力時恰恰相反。Yang 等[12]發現介尺度理論模型預測的跳躍點能反映此雙重效應：當黏度由1.0 mPa·s 增加至3.0 mPa·s 或者表面張力由20 mN/m 增加至40 mN/m 時氣含率跳躍點向高表觀氣速移動，流域過渡延遲；進一步增加黏度或者表面張力，跳躍點向低表觀氣速移動，流域過渡提前。

需要指出的是，基于EMMS 方法的介尺度理論模型是零維概念模型，反映的是體系多相流的內在本征穩定性，并不包含裝置結構(如塔徑、氣體分布器以及內構件等)外在因素的影響。對于有外在因素影響的實際問題，需要求解完整的介尺度穩定性條件約束的多流體模型(stability-constrained multifluid model, SCMF)。但是，該完整模型的求解需要在每個時間步內進行優化計算，在技術上和計算量上目前都還比較困難。可以采取一種簡化的穩定性條件約束的多流體模型，如圖3 所示。對于計算單元建立介尺度穩定性條件約束的穩態守恒方程，獲得單元結構信息參數，然后重構介尺度曳力模型或者介尺度群體平衡模型，再耦合到平均化的多流體模型控制方程中去。該方法通過介尺度本構模型這一橋梁間接實現了穩定性條件對多流體模型的約束。下面將詳細介紹介尺度曳力和介尺度群體平衡模型及模型的應用。

圖3 簡化的穩定性條件約束的多流體模型Fig.3 Simplified multi-fluid model constrained by stability condition

2 介尺度曳力模型

曳力是多相流中連續相和離散相之間重要的相間作用力，分散相受到的單位體積平均曳力大小為：

其中，β為平均相間曳力系數；uc為連續相速度；ud為分散相速度；αc為連續相含率；ρc為連續相密度；CD為有效曳力系數；dd為分散相粒徑。Yang 等[18]首次從介尺度理論模型的結構參數重構了考慮顆粒團聚物影響的氣固曳力模型，并用于循環流化床的CFD 模擬。模擬結果表明，與Wen&Yu/Ergun 曳力模型相比，該介尺度曳力模型能很好地捕捉顆粒聚團現象，床層顆粒夾帶量與實驗值更加符合，還可以預測快速流態化的噎塞現象。Wang 等[19]采用兩步法建立了曳力系數與計算網格內相含率及氣固速度的關系，模型能捕捉S 型軸向固含率分布。近些年氣固流態化體系的介尺度曳力模型在介尺度結構參數模型化[20-21]、拓寬應用流域范圍[22]、與人工智能結合[23-24]以及實際應用[25-26]等方面取得重要進展，大幅推進了該研究領域的發展及影響力。

在氣液體系，有效曳力系數和氣泡粒徑都是體系的未知參數，隨著操作參數而改變。但是，如式(1)所示，有效曳力系數CD和氣泡粒徑db比值的大小決定了多流體模型中的平均相間曳力系數β。因此，從介尺度理論模型結構參數重構這一比值足以實現與多流體模型的耦合，并不需要具體指定曳力系數的修正因子或者氣泡粒徑的大小。Chen 等[27]、Yang等[13]以及Xiao等[14]計算了不同表觀氣速下鼓泡塔的CD/db的大小，并建立了CD/db和表觀氣速Ug之間的數學關聯式，如表1所示。

Yang等[13]以及Xiao等[14]將該介尺度曳力模型應用到鼓泡塔的CFD 模擬。如圖4 所示,將介尺度曳力模型和CFD 模型耦合后，總體氣含率隨表觀氣速增加而增大，但進入湍動流域后總體氣含率隨表觀氣速變化不顯著，出現平臺現象。傳統曳力模型只能預測氣含率單調遞增，而介尺度曳力模型與實驗結果更加相符。僅用基于EMMS的介尺度理論模型總體氣含率有跳躍現象，表明流型發生了過渡。但理論模型是零維概念模型，僅能反映體系內在的本征穩定性，不能考慮設備結構外在因素的影響。將理論模型導出的介尺度曳力模型與CFD 結合后，可以考慮外在因素的影響，提高了總體氣含率的預測。Xiao 等[32]還考察了將鼓泡塔系統分成氣相及水相的兩相模型、稀相(大氣泡)及密相(小氣泡和水)的兩相模型和大氣泡相、小氣泡相及水相的三相模型，分別結合相應的介尺度曳力模型，均能很好地預測總體氣含率、氣含率及液速徑向分布。

圖4 高表觀氣速下總體氣含率的平臺現象[14]Fig.4 The shoulder of total holdup at high superficial gas velocity[14]

Guan 等[33]建立了不同操作模式下的介尺度曳力模型，應用于不同操作模式的鼓泡塔的模擬，在總體氣含率及氣含率分布方面獲得很好的預測效果。Zhou 等[34-35]將介尺度曳力模型拓展應用至含固體顆粒的漿態床模擬，表明低固含率下可以直接采用氣液體系的介尺度曳力模型，而高固含率情況下需要針對性的發展氣液介尺度曳力模型。Yan等[36-37]在常壓空氣-水體系的介尺度曳力模型的基礎上，考慮氣體密度、黏度或者表面張力的修正，應用于壓力或者物性不同的鼓泡塔體系模擬。

上述工作將CD/db和整體宏觀操作參數直接關聯，然后應用至整個氣液體系，稱此類模型為“整體式介尺度曳力模型”。然而，在氣液鼓泡體系內，尤其是處于湍流鼓泡區時，流動的時空非均勻性強烈，強烈的局部流動效應使得介尺度結構也呈現時空動態變化，致使有效曳力系數也將隨之動態變化。因此，發展考慮局部流動效應的“局部式介尺度曳力模型”更為合適。蔣雪冬[28]首先嘗試關聯了CD/db和局部表觀氣速及表觀液速的關系，然后將每個計算網格的表觀速度代入得到計算網格的CD/db。Guan 等[29]計算了介尺度理論模型在不同表觀氣速下的氣含率和CD/db，然后將氣含率和CD/db一一對應，進而實現CD/db和計算網格內氣含率的關聯，如表1 所示。Guan 等[29]比較了不同曳力模型對鼓泡塔總體氣含率、氣含率和液速分布的預測，表明局部式介尺度曳力模型能進一步提高對液速的預測。

表1 氣液體系介尺度曳力模型數學關聯式Table 1 Mathematical correlations for the mesoscale drag model in the gas-liquid system

介尺度曳力模型的發展經歷了從氣固體系到氣液及氣液固體系、從整體式模型到局部式模型的發展歷程。該曳力模型不含有可調參數，彌補了傳統多流體模型丟失介尺度結構信息的缺陷，突破了傳統方法的薄弱環節，從而極大地提高了多流體模型模擬的準確性。

3 介尺度群體平衡模型

聚并破碎核函數用于描述氣泡聚并破碎過程動力學，是群體平衡模型的核心。介尺度群體平衡模型的思路如圖5 所示，求解介尺度理論模型獲得系統的結構參數以及基于結構分解的介尺度能耗。將結構參數代入簡化群體平衡模型，可獲得分散相粒徑分布，進而得到基于該粒徑分布的介尺度能耗。該介尺度能耗與介尺度理論模型得到的介尺度能耗應該相等，采用此物理限制條件確定聚并破碎速率。換而言之，將介尺度能耗作為系統聚并破碎過程的物理限制，用于改進傳統聚并破碎核函數。體系的介尺度能耗是由穩定性條件決定的，因此改進的聚并破碎核函數也受穩定性條件的約束。

圖5 介尺度群體平衡模型[38]Fig.5 Mesoscale population balance model[38]

Yang 等[38]以及An 等[39]針對鼓泡塔，確定了采用Luo-Svendsen 破碎模型分別結合Prince-Blanch 聚并模型以及Luo 聚并模型的聚并修正因子，并建立了修正因子與表觀氣速的關聯。如圖6 所示，無論是Luo 聚并模型還是Prince-Blanch 聚并模型，采用基于介尺度理論模型的修正因子均能給出較好的氣泡粒徑分布，而原始模型或者采用Bhole 修正因子時的準確性對聚并模型有依賴性，適用性有限。王玨等[40]發現介尺度群平衡模型能更好地預測鼓泡塔內不同高度的氣泡粒徑分布和軸向液速。An等[39]將介尺度群體平衡模型應用至漿態床的氣泡粒徑分布模擬，并分析了固體顆粒對物性、液膜排干時間和湍流的影響，表明固體顆粒對湍流的衰減作用是改變氣泡粒徑分布的主要原因。

圖6 Luo-Svendsen破碎模型結合不同聚并模型時的氣泡粒徑分布預測[38]Fig.6 Bubble size distribution predicted by different coalescence model combined with Luo-Svendsen breakage model[38]

聚并破碎過程復雜，影響因素很多，造成現有的聚并破碎核函數需要針對不同操作條件、裝備結構以及物化性質調整模型參數。介尺度群體平衡模型采用介尺度能耗對群體平衡模型進行了額外的物理限制，避免了隨意調整模型參數，提高了對粒徑分布的預測。

4 模型的應用

介尺度曳力和群體平衡模型已經拓展應用至氣升式環流反應器、氣液攪拌槽、液液乳化設備等不同工藝過程核心設備的模擬、優化以及放大。氣升式環流反應器包括內環流反應器和外環流反應器兩類。根據下降管內氣體存在形態劃分為三種流動形態[41]：無氣體攜入流域、氣體攜入而不循環流域以及氣體完全循環流域。能模擬下降管內氣體的形態是準確預測環流反應器流動規律的前提。Xu等[42]以及張佳寶等[43]將局部式介尺度曳力模型應用至內環流反應器，表明模型能很好地模擬下降管內氣體形態以及氣含率，而其他曳力模型預測結果表明下降管內未觀察到氣體，處于無氣體攜入流域，與實驗現象不一致。Jiang 等[44]將該介尺度曳力模型應用至外環流反應器，表明模型能很好地預測上升管內的總氣含率、局部氣含率、液速和氣泡速度。

氣液攪拌槽是一類帶機械攪拌的氣液反應器，根據氣體的運動形態劃分為氣泛區、載氣區和完全循環區三種流型[45]。肖頎等[46]比較了整體式介尺度曳力模型和Schiller-Naumann 曳力模型在三種流型下攪拌槳附近排出區的流動預測效果。表明在完全循環區流域，Schiller-Naumann 曳力模型明顯低估槳葉下方的氣含率，而介尺度曳力模型能更準確描述此流域的槽內流動情況。李新菊等[47]將局部式介尺度曳力模型應用至氣液攪拌槽，并與Tomiyama曳力模型預測結果相比較。模擬結果表明Tomiyama 曳力模型只能可靠地預測低轉速(140 r/min，280 r/min)條件下循環區局部氣含率，而介尺度曳力模型能適用于整個實驗測量的轉速范圍。Guan 等[48]比較了整體式和局部式介尺度曳力模型、Brucato 曳力模型以及Schiller-Naumann 曳力模型在預測流型方面的可靠性，表明Brucato曳力模型預測的流域過渡提前，而Schiller-Naumann 曳力模型流域過渡延遲，而介尺度曳力模型能很好地預測不同操作條件下的操作流域。如圖7 所示，在400 r/min的轉速下，實驗處于載氣區，即氣體在攪拌槳上方能有效分散，但氣體不能被攜帶至攪拌槳下方，介尺度模型能給出準確的流型，而Schiller-Naumann模型預測是氣泛區，氣體在攪拌槳上方未得到有效分散。在600 r/min 的轉速時，實驗處于完全循環區，即氣體能帶到攪拌槳下方形成循環，介尺度模型的預測確實有氣體能攜帶入攪拌槳下方，處于完全循環區，而Schiller-Naumann 模型預測的是載氣區流型，攪拌槳下方未觀察到氣體。

圖7 氣液攪拌槽操作流域的預測[48-49](模擬圖左邊：氣含率云圖；右邊：氣體速度矢量圖)Fig.7 Predicted flow regime in gas-liquid stirred tanks[48-49](left:gas holdup contour;right:gas velocity vector)

介尺度群體平衡模型用于預測分散相粒徑分布。Qin 等[50]以及Chen 等[51-52]將介尺度群體平衡模型拓展應用至轉子-定子乳化設備，表明該方法能改善模型對不同操作條件下液滴尺寸分布、索特平均直徑、中位直徑以及跨度的預測。Chen 等[52]將介尺度群體平衡模型與表面活性劑在體相及界面上分配的輸運模型相耦合，建立了跨層次的介尺度模型，揭示了相界面層次和設備層次介尺度結構的耦合機制，深化了對乳化過程流動-傳遞的科學認識。如圖8 所示，僅考慮相界面層次或者設備層次的介尺度機制均明顯高估液滴尺寸，而將兩個層次的介尺度耦合起來將大幅度提高液滴尺寸分布的預測準確性，與實驗結果更加相符。

圖8 跨層次介尺度模型對液滴尺寸分布的預測[52][30%（質量）分散相，乳化劑0.1%（質量），轉速9000 r/min]Fig.8 Prediction of droplet size distribution with levelcrossing meso-scale model[52][30%(mass)dispersed phase,0.1%(mass)emulsifier,rotational speed 9000 r/min]

5 結論及展望

相間作用力模型、湍流模型以及聚并破碎核函數用于描述微觀的相內或者相間的相互作用對宏觀流動的影響，這些封閉子模型直接決定了多流體模型模擬結果的準確性。本文首先概述了能量最小多尺度方法及其在氣固流態化和氣液多相流中的應用，然后總結了基于介尺度穩定條件的兩類介尺度封閉模型(介尺度曳力模型和介尺度群體平衡模型)的構建方法，進而綜述了這些封閉模型在流化床、鼓泡塔、氣升式環流反應器、攪拌槽、轉定子乳化器等多相流設備中的應用，主要結論如下。

（1）從介科學角度發展介尺度物理模型，解析多相流非均勻結構演化的控制機制，并將介尺度模型與CFD 相結合，為認識多相流復雜性提供了新思路以及新途徑；

（2）對多相體系進行尺度分解以及控制機制分解，建立不同控制機制競爭協調的介尺度穩定性條件，能深入認識體系宏觀流域過渡的物理機制，比如氣固體系的噎塞現象，氣液體系向湍流鼓泡區的過渡；

（3）基于穩定性條件約束的介尺度曳力模型和介尺度群體平衡模型由于彌補了介尺度結構信息，突破了傳統模型的薄弱環節，使得模型在不含有可調參數的情況下，能夠準確預測流場、流域過渡以及分散相粒徑分布等；

（4）介尺度封閉模型在流化床、鼓泡塔、氣升式環流反應器、攪拌槽、轉定子乳化器等多相流設備的應用取得很好的效果，大幅度提高了CFD 的預測能力。

基于以上文獻綜述和總結，建議未來的研究方向應包括以下方面。

（1）介尺度穩定性條件的基礎科學問題：穩定性條件反映了控制機制間的競爭協調，數學上表達為多目標優化，如何求解該數學問題，如何分離體系的控制機制，是否存在普世性的控制機制建模方法，氣液及氣固體系的穩定性條件是否存在相似性，相界面層次的穩定性條件如何構建，是否可能通過穩定性條件實現跨層次的耦合。

（2）氣液多相湍流包括剪切湍流產生機制和氣泡誘導湍流產生機制，發展氣液體系的介尺度湍流模型，并實現該模型的介尺度穩定性條件與CFD 的結合是未來的重要研究方向。

（3）目前通過介尺度封閉模型實現穩定性條件和CFD 的結合，進一步探索其他的動態、穩態或者半動態的結合方式。

（4）研究介尺度封閉模型和DNS、過濾模型之間的數學物理關系，能否實現三者的有機統一。

（5）進一步強化介尺度模型與人工智能的結合。目前在介尺度封閉模型的關聯方面實現了與人工智能的初步結合。在介尺度模型中的控制機制分離以及穩定性條件構建能否實現與人工智能的結合。

符 號 說 明

CD——曳力系數

db——氣泡粒徑，m

FD——平均曳力，N/m3

f——氣含率

U——表觀速度，m/s

u——速度，m/s

α——相含率

β——平均相間曳力系數，N·s/m2

ρ——密度，kg/m3

下角標

c——連續相

d——分散相