基于歐拉-歐拉方法的氣液兩相流數(shù)值模型發(fā)展綜述1)

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

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

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

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

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

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

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

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

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

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

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

1 相間力模型

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

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

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

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

1.1 曳力

目前普遍認(rèn)為曳力是氣液兩相間最重要的作用力，它來(lái)源于氣泡與液體的相對(duì)運(yùn)動(dòng)，是兩相間動(dòng)量傳遞的主要作用力。曳力的大小取決于氣泡的形狀、流體的特性以及氣液相對(duì)速度[6]。單個(gè)氣泡所受到的曳力表達(dá)式為

式中，CD，ρl，ur和Ad分別表示氣泡的曳力系數(shù)、液體密度、兩相間的相對(duì)速度和氣相投影面積。由式（4）可以看出，若要計(jì)算出曳力的數(shù)值，需要對(duì)曳力系數(shù)CD進(jìn)行明確定義。但是對(duì)于不同流動(dòng)狀態(tài)的流場(chǎng)，曳力系數(shù)存在不確定性，因此研究者對(duì)于曳力模型的研究也就主要集中于曳力系數(shù)上。對(duì)氣泡的曳力系數(shù)模型的描述通常引入一些無(wú)量綱參數(shù)，如氣泡Reb（Reynolds數(shù)）、Eo（E?tv?s數(shù)）、Mo（Morton數(shù)）、We（Weber數(shù)），并通常將曳力模型定義為它們的函數(shù)。表1中總結(jié)了文獻(xiàn)中較為常用的氣泡曳力系數(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ù)的兩個(gè)不同范圍來(lái)分別定義曳力系數(shù)，對(duì)于Reb≤ 1 000的情況，CD定義為氣泡雷諾數(shù)的函數(shù)，而當(dāng)Reb> 1 000時(shí)，CD則定義為常數(shù)值0.44。Schiller-Naumann模型通常適用于尺寸小于2.5 mm的球形氣泡[12]，而對(duì)于較大的氣泡尺寸則適合使用Ishii-Zuber模型。1979年Ishii等[9]利用混合物黏度的概念對(duì)液滴和顆粒的曳力關(guān)系進(jìn)行了深入的研究，發(fā)現(xiàn)曳力系數(shù)可以用不同形狀的氣泡和不同的氣泡流型加以區(qū)分，如扭曲氣泡區(qū)域、帽狀氣泡區(qū)域和橢球氣泡區(qū)域，此外，Ishii 等[9]、Miao等[23]和Xu等[24]還引入了氣泡Eo，從而更準(zhǔn)確地描述氣泡形狀與曳力的關(guān)系。同樣是考慮了氣泡形狀，Grace曳力模型在考慮Eo的同時(shí)，還考慮了氣泡Mo的影響[8]。

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

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

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

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

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

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

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

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

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

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

1.2 升力

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

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

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

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

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

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

1.3 壁面潤(rùn)滑力

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

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

表3 氣泡壁面潤(rùn)滑力模型Table 3 Bubble wall lubrication force models

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

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

1.4 湍流擴(kuò)散力

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

Lopez de Bertodano[43]通過(guò)假設(shè)氣泡運(yùn)動(dòng)近似于空氣分子在大氣中的熱擴(kuò)散而對(duì)簡(jiǎn)單湍流彌散力建立模型

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

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

其中，Prα表示體積分?jǐn)?shù)離散的湍流普朗特?cái)?shù)，γl表示液相的運(yùn)動(dòng)渦流黏度。

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

1.5 虛擬質(zhì)量力

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

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

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

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

2 群體平衡模型

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

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

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

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

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

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

2.1 氣泡聚并模型

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

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

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

2.2 氣泡破碎模型

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

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

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

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

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

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

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

3 結(jié)語(yǔ)

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

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

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

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

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