高一博 耿琳琳 王 振 張振楠 張德勝 )

*(江蘇大學(xué)流體機械工程技術(shù)研究中心，江蘇鎮(zhèn)江 212013)

?(大連理工大學(xué)數(shù)學(xué)科學(xué)學(xué)院，遼寧大連 116024)

**(中國石油大學(xué)（華東）儲運與建筑工程學(xué)院，山東青島 266580)

海洋中蘊藏著豐富的石油、天然氣等能源資源，約占油氣資源總量的34%[1]，而其中70%處于深海區(qū)域，開發(fā)潛力以及開采難度十分巨大。為了解決我國油氣對外依存度過高的問題，需要盡快推進我國深海油氣勘探開發(fā)進程，深海油氣開采已成為保障我國能源安全的重要舉措。

現(xiàn)如今，在深海油氣輸運中，氣液混輸技術(shù)大幅簡化了開采工藝流程，節(jié)省了管道鋪設(shè)成本、提高了油井采收率，具有高效、節(jié)能、環(huán)保的優(yōu)點，已逐漸成為深海油氣輸運的首選方案。深海油氣的主要成分為天然氣和石油的氣液混合物，因此油氣輸送本質(zhì)上是一種氣液兩相流混輸?shù)目茖W(xué)問題。

對于長距離輸送管道，內(nèi)部多相流受氣液流量、管道尺寸及空間布置等多因素影響，呈現(xiàn)出多種流型。如圖1所示，在豎直管道中存在五種常見流型：分散氣泡流、段塞流、混狀流、環(huán)形流、環(huán)形液滴流。而在水平管道中由于重力的作用，又表現(xiàn)為與垂直管道不同的流型，具體又可細分為分散氣泡流、環(huán)形氣泡流、拉伸氣泡流、段塞流、分層流、分層波浪流等[2]。而且，在深海油氣開采進程中，由于復(fù)雜的工況環(huán)境、瞬變的來流條件，以及含氣率的日益升高，導(dǎo)致混輸管內(nèi)氣液兩相流動狀態(tài)不斷發(fā)生變化，具體表現(xiàn)為多種氣液流型并存，相間結(jié)構(gòu)及其分布狀態(tài)時空各異，各種動力學(xué)參數(shù)、流動特性隨流型轉(zhuǎn)換頻繁變化等等。因此，如何準確地預(yù)測這種復(fù)雜環(huán)境下的氣液兩相流動是油氣輸運數(shù)值仿真領(lǐng)域的關(guān)鍵科學(xué)問題。

圖1 管內(nèi)氣液兩相流型示意圖[2]Fig.1 Gas-liquid two-phase flow pattern in pipe[2]

目前，在數(shù)值上處理氣液兩相流的方法主要分為歐拉-拉格朗日和歐拉-歐拉兩種。歐拉-拉格朗日方法將液相視為連續(xù)的流體，將氣相視為拉格朗日質(zhì)點，通過求解質(zhì)點的運動微分方程得到離散相在流場中的運動軌跡與運動狀態(tài)。該方法的最大優(yōu)點是可以得到離散氣泡在流場中完整的運動信息，以及由于受力引起的尺寸演變；缺點是其僅適用于離散相體積分數(shù)比較低的多相流動，并且往往需要消耗巨量的計算資源。而歐拉-歐拉方法則引入體積分數(shù)的概念，將每一相都假定為可互相穿透的連續(xù)介質(zhì)，每一相都擁有各自獨立的壓力場、速度場以及溫度場等物理屬性，在保證一定精度的前提下，相較于歐拉-拉格朗日方法大幅降低了計算成本，目前在多相流領(lǐng)域中被廣泛應(yīng)用。

此外，隨著氣液兩相中含氣率的升高，氣泡與氣泡之間會頻繁發(fā)生相互作用，具體表現(xiàn)為由于氣泡的聚并和破碎引起的氣泡產(chǎn)生與消亡，這導(dǎo)致了氣泡尺寸在流場中呈現(xiàn)多尺度分布特性。而在常規(guī)的歐拉-歐拉方法中，往往采用的是均一氣泡直徑，忽略了氣泡的多尺度分布及泡間相互作用，從而導(dǎo)致常規(guī)方法對于高含氣率氣液流場的解析度有限。目前，能夠較好預(yù)測氣泡多尺度分布的方法是通過引入群體平衡模型（population balance model，PBM）來模擬氣泡之間的相互作用。該模型也將在下文中進行詳細介紹。

在本文中，將對歐拉-歐拉框架下氣液多相流研究中所使用的主要相間力模型及其發(fā)展情況進行介紹。并針對該方法在深海油氣混輸過程中的高壓高含氣率特殊環(huán)境中的局限性進行闡述，展望其未來的發(fā)展趨勢，以期為我國深海油氣輸運數(shù)值模擬研究的發(fā)展提供有益的參考。

1 相間力模型

正確描述液體中氣泡的受力情況，對于研究氣液兩相流的復(fù)雜系統(tǒng)具有重要意義，而歐拉-歐拉模型模擬氣液兩相流動狀況需要解析相間作用力來封閉控制方程。在油氣混輸過程中，一般忽略氣液兩相的相間傳質(zhì)作用，控制方程為

其中αi，ρi，ui分別表示為第i相的宏觀體積分數(shù)、密度和速度，F(xiàn)ij表示氣液兩相的相間作用力，主要用于表示相間動量的傳遞，它實際上是多種力共同作用的結(jié)果[3]

其中，等號右邊分別代表五種不同的相間作用力，其產(chǎn)生機理如圖2所示。FD代表曳力，是由氣相與液相之間的相對運動引起的；FL代表升力，是由垂直于氣泡運動方向的壓差所引起的側(cè)向力；FWL代表壁面潤滑力，是由壁面附近的滑移速度產(chǎn)生的；FTD代表湍流擴散力，代表由液體湍動對氣泡引起的牽引力；FVM代表虛擬質(zhì)量力，是由相對加速度的變化而產(chǎn)生的力。相間力的選用會直接影響到數(shù)值計算的收斂性，所以進一步優(yōu)化相間作用力模型是該方向的突破口[4]。

圖2 氣液相間作用力示意圖[5]Fig.2 Illustration of the interphase forces in bubbly flow[5]

1.1 曳力

目前普遍認為曳力是氣液兩相間最重要的作用力，它來源于氣泡與液體的相對運動，是兩相間動量傳遞的主要作用力。曳力的大小取決于氣泡的形狀、流體的特性以及氣液相對速度[6]。單個氣泡所受到的曳力表達式為

式中，CD，ρl，ur和Ad分別表示氣泡的曳力系數(shù)、液體密度、兩相間的相對速度和氣相投影面積。由式（4）可以看出，若要計算出曳力的數(shù)值，需要對曳力系數(shù)CD進行明確定義。但是對于不同流動狀態(tài)的流場，曳力系數(shù)存在不確定性，因此研究者對于曳力模型的研究也就主要集中于曳力系數(shù)上。對氣泡的曳力系數(shù)模型的描述通常引入一些無量綱參數(shù)，如氣泡Reb（Reynolds數(shù)）、Eo（E?tv?s數(shù)）、Mo（Morton數(shù)）、We（Weber數(shù)），并通常將曳力模型定義為它們的函數(shù)。表1中總結(jié)了文獻中較為常用的氣泡曳力系數(shù)模型。

表1 氣泡曳力模型（續(xù)）Table 1 Bubble drag force models (continued)

表1 氣泡曳力模型Table 1 Bubble drag force models

式中，db和ρg分別表示氣泡的直徑和密度，σ表示氣液界面張力，g表示重力加速度。

現(xiàn)如今，例如Schiller-Naumann模型[7]、Grace模型[8]、Ishii-Zuber模型[9]、Tomiyama 模型[10]等曳力模型已被廣泛應(yīng)用于氣液兩相流數(shù)值模擬[11]，并已植入商用CFD軟件中供用戶選擇。Schiller等[7]在1935年提出的模型是最為經(jīng)典的曳力模型，其根據(jù)氣泡雷諾數(shù)的兩個不同范圍來分別定義曳力系數(shù)，對于Reb≤ 1 000的情況，CD定義為氣泡雷諾數(shù)的函數(shù)，而當Reb> 1 000時，CD則定義為常數(shù)值0.44。Schiller-Naumann模型通常適用于尺寸小于2.5 mm的球形氣泡[12]，而對于較大的氣泡尺寸則適合使用Ishii-Zuber模型。1979年Ishii等[9]利用混合物黏度的概念對液滴和顆粒的曳力關(guān)系進行了深入的研究，發(fā)現(xiàn)曳力系數(shù)可以用不同形狀的氣泡和不同的氣泡流型加以區(qū)分，如扭曲氣泡區(qū)域、帽狀氣泡區(qū)域和橢球氣泡區(qū)域，此外，Ishii 等[9]、Miao等[23]和Xu等[24]還引入了氣泡Eo，從而更準確地描述氣泡形狀與曳力的關(guān)系。同樣是考慮了氣泡形狀，Grace曳力模型在考慮Eo的同時，還考慮了氣泡Mo的影響[8]。

Tomiyama等[10]針對單個氣泡的曳力系數(shù)進行了實驗研究，考慮了氣泡的形狀，提出了一個具有一般性的曳力系數(shù)關(guān)聯(lián)式，但是在將Tomiyama模型應(yīng)用于多氣泡系統(tǒng)時，應(yīng)當考慮含氣率對曳力的影響。另外還考慮了液體的污染程度，并根據(jù)液體的純度提出了不同的相關(guān)性。Tomiyama等[10]提出的曳力模型已經(jīng)在10-2

Simonnet等[13]通過實驗研究了氣泡與氣泡之間的相互影響，通過分析較大含氣率范圍內(nèi)的氣泡群的運動狀態(tài)，提出了適用于0%～30%含氣率范圍和純水-空氣體系的曳力關(guān)系。如圖3所示，由于相鄰氣泡導(dǎo)致的氣泡曳力減小，因此將局部含氣率αg,loc和氣泡平均直徑db相關(guān)聯(lián)。氣泡狀態(tài)變化的臨界含氣率為15%，在此處由均勻狀態(tài)向非均勻狀態(tài)過渡，曳力系數(shù)的比值也由線性變化轉(zhuǎn)變?yōu)榉蔷€性變化。為了更好地描述變化的趨勢，Simonnet等[13]在曳力關(guān)系式中加入了修正因子。出真實的流場。為此，研究者將壓力的影響引入到了曳力模型中[25-26]，Krishna等[25]認為壓力對流場的影響主要在兩個方面：一是壓力的增加增大了氣相密度；二是壓力的增加增大了大氣泡的破碎概率，減小了大氣泡上升速度。因此在大氣泡的曳力模型中添加了氣相密度修正項 (ρg/ρg,0) ，從而考慮了操作壓力的影響，然而，這種修正因子的引入高估了在較高表觀氣速和較高壓力下的總體含氣率。

圖3 氣泡曳力系數(shù)與局部含氣率的關(guān)系[13]Fig.3 Bubble drag coefficient vs.local void fraction[13]

Chen等[26]對壓力范圍0.1～1.0 MPa的氣泡流進行了數(shù)值模擬，考慮了氣泡聚并與破碎對氣泡尺寸分布的影響，在分析Krishna等[25]引入的氣相密度修正項的基礎(chǔ)上，在曳力模型里添加了新的氣相密度修正項 (ρg/ρg,0)0.25，得到的結(jié)果中不同壓力下的局部氣含率及軸向液速的徑向分布都與實驗結(jié)果吻合較好。Yang等[17]提出一種根據(jù)氣泡尺寸分布的平均曳力系數(shù)模型，它同時考慮了小氣泡的阻礙效應(yīng)和大氣泡的尾流加速效應(yīng)，對于小氣泡的考慮實質(zhì)上是一種對高壓環(huán)境下的曳力修正，因此可適用的范圍更廣。

Zhang等[18]在Roghair等[16]的氣泡群曳力模型的基礎(chǔ)上，基于冷態(tài)實驗測量的整體含氣率，在壓力0.5～2.0 MPa范圍內(nèi)對氣泡群曳力模型進行修正，采用了含有壓力修正因子的曳力系數(shù)計算公式，適用于加壓條件下改進的氣泡群曳力計算。Mach 等[20]實驗測定了四種不同壓力（0.14 MPa，0.46 MPa，2.00 MPa，4.00 MPa）下的氣泡尺寸分布和單個氣泡曳力系數(shù)，量化了氣泡尺寸分布與單個氣泡曳力系數(shù)之間的關(guān)系。在氣含率不變的情況下，壓力的增加導(dǎo)致單個氣泡的曳力系數(shù)減小。并且這些較小的氣泡會引起邊緣液體擾動，起到抑制湍流的效果，從而減少曳力。

Van Tran等[19]采用歐拉-歐拉模型模擬了壓力在 0.1 MPa，1.5 MPa和3.5 MPa條件下的空氣-水體系。采用了帶有壓力參數(shù)的修正因子修正曳力，并根據(jù)不同氣體流量和壓力下的實驗數(shù)據(jù)，對壓力修正因子進行了優(yōu)化。修正后的PBM模型在合理的精度范圍內(nèi)較好地反映了空氣-水加壓鼓泡塔的流體力學(xué)特性。Van Tran等[21]測試了最高16 MPa下減壓渣油與氫氣的氣液兩相流特性，測量了密度、黏度、表面張力等物理參

Roghair等[16]認為通過實驗同時測量多個相關(guān)參數(shù)存在大量的困難，因此采用直接數(shù)值模擬方法和歐拉-拉格朗日法，研究了密實氣泡對氣泡曳力的影響，并通過考慮局部含氣率和單個上升氣泡的相互作用提出了改進的曳力模型，適用的氣含率最高可達45%。Van Tran等[19]認為在均相流區(qū)，氣泡以小氣泡為主，小氣泡合并形成大氣泡。因此為了綜合考慮氣泡群和尾跡對曳力系數(shù)的影響，將局部氣含率作為影響因素加入到曳力系數(shù)中。

Gillissen等[14]考慮了大范圍復(fù)雜泡狀流和氣泡變形的共同作用，采用格子Boltzmann法和浸沒邊界法，對三維周期域內(nèi)浮力驅(qū)動的無質(zhì)量球形氣泡的懸浮進行了直接數(shù)值模擬。研究了在液體中上升的分散氣泡群所受到的曳力，并得到了氣泡曳力系數(shù)與氣泡韋伯數(shù)和雷諾數(shù)的關(guān)系。

Yang等[15]通過氣固能量最小多尺度模型考慮非均勻結(jié)構(gòu)對相間曳力的影響，創(chuàng)建了雙氣泡尺寸（dual-bubble-size, DBS）模型，將氣泡分為大小兩種，通過大小氣泡的氣含率，直徑和表觀速度來描述氣相的狀態(tài)，通過不同尺度的能量耗散來封閉控制方程，并將能量最小原則作為穩(wěn)定性判據(jù)。最終提出DBS曳力模型，該模型無須調(diào)整模型參數(shù)即可合理地預(yù)測徑向含氣量分布和兩相流場。Guan等[22]在Yang 等[15]的基礎(chǔ)上，對不同操作模式下的DBS曳力模型進行研究，計算了不同表觀氣速下的氣含率和曳力系數(shù)，并將其一一對應(yīng)，進一步提高了模型的預(yù)測精度。

上文中所介紹曳力模型通常只適用于常壓條件下的氣液兩相流，對于高壓環(huán)境下則無法反映數(shù)，通過擬合經(jīng)驗方程得到經(jīng)驗修正因子，并將修正因子加入到了Ishii-Zuber曳力模型中。

1.2 升力

升力受流體剪切和氣泡大小的影響，是一種沿著運動氣泡的橫向方向作用于氣泡的力。由于其考慮了流體剪切對氣相運動的影響，可用于穩(wěn)定流動中球形氣泡的橫向遷移。有學(xué)者研究表明，為了減少計算量，在數(shù)值模擬中可以忽略升力的影響[26-28]。氣泡的剪切誘導(dǎo)升力FL為

其中，CL表示升力系數(shù)，αg表示氣相的體積分數(shù)，ul表示液相速度。與曳力一樣，升力對氣泡在流場中的橫向分布也起著重要作用。數(shù)值和實驗證明，升力方向的變化與變形氣泡的尾流不穩(wěn)定性有關(guān)，氣泡的變形可以誘發(fā)升力符號的改變[29-30]。Hibiki等[31]對在各種流動條件下開發(fā)的升力現(xiàn)有模型進行了廣泛的文獻調(diào)查。與曳力的計算類似，針對不同升力模型的研究主要針對的是升力系數(shù)CL，現(xiàn)有的主要升力模型如表2所示。

表2 氣泡升力模型Table 2 Bubble lift force models

Saffman[32]首先對均勻剪切流場中單個球體的升力進行研究，并將升力系數(shù)定義為渦度雷諾數(shù)Reω的函數(shù)，適用于低雷諾數(shù)流的固體顆?；驓馀?。Mei等[33]研究了均勻剪切流動中的球形氣泡。通過在氣泡表面施加不同的邊界條件計算慣性遷移速度。利用Saffman給出的升力與遷移速度之間的關(guān)系，提出了在一定雷諾數(shù)范圍下升力的表達式，適用于湍流自由剪切流中的小球形氣泡。缺點是不能精確表示壁面上的非球形氣泡的剪切升力。Moraga等[34]在此基礎(chǔ)上提出了一種主要應(yīng)用于球形氣泡的升力模型，將氣泡渦度誘導(dǎo)升力和剪切升力結(jié)合起來，從而綜合考慮了尾流產(chǎn)生的渦流效應(yīng)。

Tomiyama等[10]對甘油溶液中的單個氣泡進行了實驗研究，提出了CL關(guān)于氣泡雷諾數(shù)Reb和修正的Eo′數(shù)的關(guān)系式，Tomiyama升力模型適用于橢球狀到球帽狀范圍內(nèi)的可變形氣泡，并驗證了當Eo′>10時升力方向會發(fā)生轉(zhuǎn)變。

Legendre等[35]對CL進行了數(shù)值研究，并提出了適用于球形氣泡的經(jīng)驗關(guān)聯(lián)式，作為無量綱剪切速率（Sr）和氣泡雷諾數(shù)Reb的函數(shù)。并提出了低雷諾數(shù)時氣泡升力的變化不同于高雷諾數(shù)的情況，Legendre-Magnaudet升力模型的適用范圍為 0.1≤Reb≤500 。

1.3 壁面潤滑力

壁面潤滑力是由靠近壁面的氣泡周圍速度分布的變化所產(chǎn)生的力[36]，受氣液表面張力的影響，防止氣泡附著在壁面上[37]。其數(shù)學(xué)表達式為

其中，CWL表示壁面潤滑力系數(shù)，urw表示相對速度在壁面處的切向分量，nw表示壁面上的單位外法向量。表3總結(jié)了研究者現(xiàn)已提出的壁面潤滑力系數(shù)計算方法。

表3 氣泡壁面潤滑力模型Table 3 Bubble wall lubrication force models

Antal等[36]首先將壁面潤滑力的概念引入到球形氣泡，假設(shè)了球形氣泡和層流條件，將壁面潤滑力與氣泡直徑和壁面距離相關(guān)聯(lián)。Antal模型適用于在氣泡到壁面的距離ywall≤5db的范圍。Tomiyama等[38]對甘油流動中的氣泡進行研究，提出了基于管道直徑Dpipe和氣泡Eo的壁面潤滑力模型，并將氣泡Eo分為四個區(qū)間區(qū)分了不同大小的壁面潤滑力，具有很好的計算精度，是一種經(jīng)典的壁面潤滑力模型。在2002年Tomiyama等[30]又對氣泡的橫向運動機理進行研究，在1 ≤Eo≤ 33的范圍內(nèi)重新定義了壁面潤滑力，并使其擺脫了管道直徑的限制。Hosokawa等[39-40]對不同濃度甘油溶液中的氣泡進行了實驗研究，發(fā)現(xiàn)壁面潤滑力系數(shù)不僅取決于氣泡Reb和Eo[39]，還與氣泡Mo數(shù)有關(guān)[40]。Hosokawa模型的適用范圍僅限于氣液兩相泡狀流，并且已經(jīng)在2.0 ≤Eo≤10.0，Reb< 10，以及 -6.0 ≤ lgMo≤ -2.5，0.4 ≤Reb≤ 65的范圍內(nèi)得到了驗證。

Frank等[41]修正了Tomiyama的壁面潤滑力模型，使該模型脫離了流型的限制，此外還引入了兩個系數(shù)，一個是阻尼系數(shù)Cdamp， 另一個是截斷系數(shù)Ccut-off。Cdamp決定了壁面潤滑力的相對大小，Ccut-off決定了壁面潤滑力起作用的相對距離，這兩個系數(shù)的引入使得Frank模型可以更準確模擬壁面潤滑力。

1.4 湍流擴散力

當連續(xù)相速度的連續(xù)隨機波動作用于氣泡時，就產(chǎn)生了湍流擴散力。它受液體中湍流渦對分散氣泡的影響，由連續(xù)相湍流渦流對分散相的綜合作用組成，在氣泡體積分數(shù)的徑向分布中起著重要作用[42]?，F(xiàn)行的湍流擴散力模型分為兩類：Lopez de Bertodano模型[43]和 Favre-averageddrag（FAD）模型[44]。

Lopez de Bertodano[43]通過假設(shè)氣泡運動近似于空氣分子在大氣中的熱擴散而對簡單湍流彌散力建立模型

其中，kl表示液相的湍動能，湍流擴散力系數(shù)CTD的取值在0.1～0.5之間[44]。

Burns等[45]通過將Favre平均值應(yīng)用于流體曳力來建立的湍流分散力模型為

其中，Prα表示體積分數(shù)離散的湍流普朗特數(shù)，γl表示液相的運動渦流黏度。

其中，斯托克斯數(shù)St是由氣泡的松弛時間tb和渦松弛時間te共同定義的

1.5 虛擬質(zhì)量力

虛擬質(zhì)量力是由于加速的氣泡對液體所做的功而產(chǎn)生的。此外，虛擬質(zhì)量力還可以表現(xiàn)為由于氣泡在連續(xù)相中加速時所產(chǎn)生的額外曳力而產(chǎn)生的力。對于加速流動的兩相系統(tǒng)，當兩相之間密度相差很大時，虛擬質(zhì)量力起到十分顯著的作用[46]。虛擬質(zhì)量力可表示為

CVM是虛擬質(zhì)量力系數(shù)，它取決于流型、離散相的幾何形狀以及體積分數(shù)。在一般的數(shù)值研究中，采用的是固定的虛擬質(zhì)量力系數(shù)。如球形氣泡的虛擬質(zhì)量力系數(shù)通常為0.5[47]。然而，對于特定的條件，例如管直徑大于0.15 m的情況，虛擬質(zhì)量力對模擬結(jié)果的影響可以忽略[47-48]。

綜上所述，歐拉-歐拉方法，通過建立不同的相間力模型研究氣液兩相的流動形態(tài)，數(shù)值模擬的效率和準確性取決于合適的界面力模型，一些研究者在求解歐拉-歐拉控制方程時同時使用了不同的相間力[19-21]。由于曳力對氣泡運動的影響明顯大于其他相間力的作用[49]，有研究者僅使用曳力模型來預(yù)測氣液兩相流體動力學(xué)性質(zhì)[26,50-51]。目前盡管這些相間力的研究已經(jīng)取得了很大進展，但在數(shù)值模擬中相間力模型的普適性仍有待進一步探究。

目前，在氣液兩相流的模擬中，較常使用的歐拉-歐拉方法往往采用的是平均氣泡直徑的處理方式，然而在現(xiàn)實的氣液兩相流動過程中，氣泡的產(chǎn)生和消失、氣泡的合并和破碎現(xiàn)象非常普遍，這導(dǎo)致了氣泡尺寸呈現(xiàn)多尺度分布，并且氣泡的大小分布還會隨著多相體系的相間傳遞現(xiàn)象而發(fā)生變化。因此，在模擬氣液兩相流時，需要對氣液兩相流中氣泡尺寸分布進行更合理的預(yù)測，進而能夠獲取界面區(qū)域濃度和動量曳力項來更新計算界面?zhèn)鬟f項，從而對氣液兩相流的動態(tài)行為進行更準確的計算[52]。針對以上流場特征和仿真目標，群體平衡模型通過求解氣泡數(shù)密度輸運方程來考慮氣泡破碎聚并的影響，且通過耦合歐拉-歐拉模型對氣液兩相流進行模擬已成為氣液兩相流領(lǐng)域研究的熱點。以下將對群體平衡模型基本理論進行介紹與討論。

2 群體平衡模型

群體平衡模型最早由Hulburt等[53]提出，近年來該模型在氣液多相流的模擬中應(yīng)用越來越廣泛。對于氣液兩相流體系，群體平衡模型的有效利用能夠描述不同尺寸氣泡在液相中的分布，對于更深入了解氣液多相流運動機理有很好的幫助?？紤]氣泡破碎聚并的群平衡模型的輸運方程為

其中，Gv表示氣泡的生長速率，BA和DA是與氣泡合并有關(guān)的產(chǎn)生項和消失項，BB和DB分別是與氣泡破碎有關(guān)的產(chǎn)生項和消失項，V表示破碎后子氣泡體積，n(V)表示體積為V的氣泡數(shù)密度函數(shù)。

其中，V′表示破碎前氣泡體積，gb(V′) 表示破碎頻率，βb(V|V′) 表示氣泡破碎的概率密度函數(shù)，Γc(V,V′)表示氣泡聚并速率，現(xiàn)有的聚并破碎模型一般需要碰撞頻率hc和聚并效率λc的模型來描述聚并頻率。表示為

破碎模型一般表示為氣泡破碎頻率gb和子氣泡尺寸分布βb的模型來模擬破碎，這些模型的建立與流體和流動屬性有關(guān)[54]。表4和表5中列出了一些適用于氣泡流的聚并和破碎模型。

表4 氣泡聚并模型Table 4 Bubble coalescence models

表5 氣泡破碎模型Table 5 Bubble breakage models

2.1 氣泡聚并模型

氣泡碰撞主要有四種機制：湍流誘導(dǎo)碰撞、浮力誘導(dǎo)碰撞、尾跡卷吸碰撞和黏性剪切碰撞。Coulaloglou等[55]首先提出了氣泡的湍流誘導(dǎo)碰撞機制，假設(shè)聚并主要在平緩的相界面上進行，還必須滿足碰撞的條件，并且充分接觸以便在這些過程中發(fā)生膜排出、膜破裂和合并過程。Coulaloglou-Tavlarides聚并模型同樣適用于圓球狀的液滴和小于湍流微尺度的液滴，該建模思想被后來的研究者所采用。Prince 等[56]詳細定義了氣液體系中發(fā)生氣泡聚并的必要因素，考慮了由湍流、浮力和黏性剪切引起的碰撞，并假設(shè)不同機制的碰撞是累積的，用氣體分子運動論的方法在Coulaloglou等[55]的基礎(chǔ)上建立了適用于氣泡-水體系的碰撞頻率模型。

Luo等[57]基于能量守恒分析提出了估算氣泡接觸時間的模型。模型只能對小尺寸的液滴給出很好的預(yù)測。Lehr 等[58]根據(jù)前人的理論并通過實驗發(fā)現(xiàn)，當兩個氣泡之間的相對速度大到一定程度時，發(fā)生碰撞后并不會產(chǎn)生聚并行為而是相互彈開。根據(jù)這一現(xiàn)象假設(shè)氣泡聚并發(fā)生的條件是碰撞速度低于一個臨界值ucrit，當氣泡碰撞速度高于臨界值時氣泡發(fā)生反彈而非聚并，Lehr模型同時考慮了兩個氣泡相對速度和流體湍流渦兩種碰撞機理。

Zhang等[18]針對壓力高達2 MPa的空氣-水體系，考察了壓力對氣液兩相流動的影響，認為壓力對氣泡的聚并也有著一定的影響，根據(jù)冷態(tài)實驗數(shù)據(jù)，在Luo聚并模型的基礎(chǔ)上新增了修正系數(shù)Ce=0.319ln(ρ/ρ0)+0.665 ，由于是通過數(shù)據(jù)擬合得到，該模型屬于一種半理論半經(jīng)驗的模型。

2.2 氣泡破碎模型

氣泡破碎主要有四種機制：湍流渦碰撞，黏性剪切力，尾渦剪切脫落，大氣泡表面不穩(wěn)定。這四種機制中，以湍流渦碰撞為主，Coulaloglou等[55]根據(jù)各向同性湍流理論，通過假定子液滴的運動與湍流的運動具有類似的性質(zhì)，從而確定了破碎時間，認為液滴破碎的決定條件是液滴湍動能大于表面能，并且預(yù)測破碎產(chǎn)生的子液滴尺寸呈現(xiàn)正態(tài)分布。然而對于氣泡-水體系，Coulaloglou等[55]的模型并不能很好地預(yù)測。Prince等[56]對水中氣泡的破碎進行實驗研究，忽略了小于0.2db的湍流渦體對氣泡破碎的影響，將氣泡和湍流渦的碰撞過程類比為氣體分子運動，認為氣泡破碎速率等于氣泡和湍流渦體碰撞頻率與破碎效率之積。

Luo等[57]基于Prince等[56]的理論，認為在氣液兩相中湍流渦的湍動能量大于氣泡破碎后表面能的增量時才會引起破碎，并直接由氣泡破碎速率函數(shù)導(dǎo)出子氣泡尺寸分布的表達式，因此當子氣泡的尺寸趨近0時，氣泡的破碎概率會無限增加。Laakkonen等[59]為了避免這一現(xiàn)象的發(fā)生，提出一種在數(shù)學(xué)上更簡潔的尺寸分布模型，使無限小的子氣泡和等于原始氣泡的子氣泡的破碎概率為0，更有利于計算。Lehr等[58,61]和Laakkonen等[59]將壓力的影響加入到破碎模型中，認為只有在湍流渦產(chǎn)生的動壓大于破碎后最小氣泡的內(nèi)部壓力時，母氣泡才會發(fā)生破碎。氣泡破碎概率密度為湍流渦動壓與氣泡表面力的比值。

自從群體平衡模型被引入到多相流模擬中來，模型的建立和求解等方面都得到了較大的發(fā)展，并根據(jù)對氣泡破碎和聚并機理的研究，建立了許多針對不同機理的氣泡破碎或聚并模型，但對于適用于加壓條件的模型尚處于理論研究階段，文獻中報道壓力的變化主要影響氣泡的破碎行為，對氣泡聚并行為影響較小[62]。

Wang等[63]認為湍流渦碰撞和大氣泡表面不穩(wěn)定性導(dǎo)致的氣泡破裂是主要的。通過將這兩種機制產(chǎn)生的速率相加，計算出總的氣泡破裂速率。Yang等[17]在此基礎(chǔ)上考慮操作壓力對氣泡破裂的影響，通過實驗測量氣體密度對不同大小氣泡破碎率的影響，將擬合得到的經(jīng)驗壓力修正因子加入到氣泡破碎模型（圖4）。

圖4 氣體密度對不同大小氣泡破碎率的影響[17]Fig.4 Effect of gas density on the bubble breakup rate for different bubble sizes[17]

Xing等[60]對破碎過程中氣泡/液滴頸部的內(nèi)部流動進行了機理分析，認為氣泡/液滴內(nèi)部的流動行為強烈依賴于壓力或氣體密度，基于這一機理，提出了一種氣泡和液滴的統(tǒng)一破碎模型。Zhang等[64]系統(tǒng)地研究了群體平衡模型的通用性，通過靜態(tài)應(yīng)力分析、氣泡頸部的界面應(yīng)力和黏性流動曳力充分描述變形氣泡內(nèi)的流動，結(jié)合內(nèi)外壓差對氣泡頸部的影響，重新解釋了氣泡破裂的因素，提出一種不需要針對特定的模擬情況調(diào)整模型參數(shù)的氣泡破碎模型，并已證明在0.1～1.0 MPa的壓力范圍內(nèi)的良好適用性。

Van Tran等[65]提出氣泡的破碎率隨著壓力的增加而增加，并強調(diào)在壓力較高的情況下壓力修正因子的重要性，基于Lehr等[61]提出的破碎模型進行改進，加入了修正系數(shù) (P/P0)1.5，在0.1～3.5 MPa的壓力范圍很好地預(yù)測了均勻泡狀流區(qū)的總氣含率和平均氣泡尺寸。Van Tran等[66]比較了常用的Luo，Lehr和Wang模型在總破損率、子粒尺寸分布和計算時間等方面的優(yōu)劣。在Luo破碎核的基礎(chǔ)上，考慮了壓力對氣泡破碎率的影響，提出了一種基于密度變化的修正系數(shù)0.35(ρg/ρg,0)0.3，并驗證發(fā)現(xiàn)改進后的Luo模型更適用于加壓環(huán)境。

3 結(jié)語

高含氣率下氣液多相流的數(shù)值模擬研究對于深海油氣混輸具有重要意義。采用歐拉-歐拉模型耦合群體平衡模型的數(shù)值方法已成為該領(lǐng)域研究的熱點與重難點。然而，該數(shù)值方法的準確性往往依賴于眾多子模型的封閉與選取，在過去的幾十年里，眾多研究者在該方面做了大量的工作，本文進行了如下總結(jié)。

（1）回顧了現(xiàn)有的氣液相間力（包括曳力、升力、湍流擴散力、壁面潤滑力和虛擬質(zhì)量力）模型，對各個模型的相關(guān)理論與發(fā)展進行了系統(tǒng)總結(jié)。結(jié)果表明盡管相間力模型在過去的幾十年中已取得長足進展，但仍存在準確性依賴于經(jīng)驗性的可調(diào)參數(shù)、模型普適性有限、計算量大等問題，且面對深海高壓高含氣率的特殊氣液混輸應(yīng)用需求，相間力模型的發(fā)展仍有待深入研究。

（2）對氣液流動的群體平衡模型的子模型進行了總結(jié)，包括氣泡的聚并模型和破碎模型，大量研究集中于氣泡的破碎模型，大多數(shù)的破碎模型都基于氣泡與湍流渦碰撞的假設(shè)，但是忽略壓力的因素導(dǎo)致模型與實驗結(jié)果存在一定差距，因此，建立考慮壓力對氣泡破碎的影響機理并且能夠統(tǒng)一描述氣泡破碎的模型，對于群體平衡模型的發(fā)展將有重要的意義。

（3）壓力的變化對于氣泡的動力學(xué)模型有著很大的影響，總結(jié)了一些在高壓環(huán)境下對氣液兩相數(shù)值模型的處理方法，通常是將壓力的影響引入到曳力模型與破碎模型中，對其進行修正，但多為經(jīng)驗性的，因此還需要進行針對性研究使模型具備良好的普適性。

綜上，目前針對相間力模型以及群體平衡模型的擴展是比較前沿的，在許多文獻中都可以找到很多不同的思路，并都能得到很好的應(yīng)用和驗證。但是在油氣混輸?shù)膹?fù)雜環(huán)境下，常規(guī)相間力模型和群體平衡模型都需要進行特殊處理。盡管在建立和創(chuàng)新氣液相間力和群體平衡模型方面已有學(xué)者做了大量的工作，但對于深海環(huán)境油氣混輸?shù)忍厥夤r下的氣液兩相流數(shù)值模擬，仍需進一步努力以提高模型的適用性與準確性。