液-液流動出現(xiàn)在各種工業(yè)過程中，特別是在石油工業(yè)中，氣體與伴生液體(油、凝析油和/或水)的混合物一起生產和運輸。當可變形界面在管道中共流時，它們獲得了flow regimes/flow patterns的各種特征分布。理解從分層流到段塞流的轉變在油氣輸送系統(tǒng)中是很重要的，因此其不斷地推動著這一方向的研究。段塞流是水平和近水平氣液流動中常見的一種流型，其處于大相干擾動狀態(tài)，從而引起大的壓力波動和流量變化，將影響工藝設備。其另一個特點是間歇出現(xiàn)充氣液體團塊，這些液體團塊完全填滿了管道的橫截面，并以很大的速度向下游流動。

本文通過兩類CMFD方法：CFX軟件的雙流體模型與VirtualFlow的Level Set方法進行數(shù)值模擬以解決段塞流形成的問題。首先比較兩種預測二維管道中分層流動的方法，隨后在只使用Level Set方法的計算結果討論了在湍流條件下圓形三維管道中分層流向段塞流過渡的結果。湍流通過由Liovic和Lakehal(2007)擴展到界面兩相流的LES方法來建模，即自由界面大渦模擬(LEIS)。

1.模型描述

1.1多物理場計算軟件VirtualFlow

VirtualFlow是基于結構化多塊網(wǎng)格求解多相流體N-S方程的多物理場、有限體積法流體力學計算軟件。針對多塊問題，采用基于MPI并行算法，網(wǎng)格排列并置。采用壓力基求解器，使用Karki-Patankar技術對弱可壓縮流進行校正。采用三階龍格-庫塔顯式時間積分格式求解N-S方程和Level Set方程。對流通量采用Quick格式離散，擴散通量采用二階中心差分格式離散。對多相流的建模采用Level Set和VOF層流/湍流模型實現(xiàn)。VirtualFlow的數(shù)值計算過程幾乎是完全自動化實現(xiàn)，前處理使用浸沒邊界方法(IST)將復雜系統(tǒng)組件映射到簡單矩形笛卡爾網(wǎng)格中。近墻區(qū)域采用BMR(基于塊的網(wǎng)格細化)處理，其中細化網(wǎng)格塊放置在適當?shù)奈恢?。塊之間的連接可以并行實現(xiàn)(使用MPI)，最多可實現(xiàn)8對1的單元映射。IST/BMR組合可以在三維模型中節(jié)省高達70%的網(wǎng)格單元。該技術也允許直接進行共軛傳熱建模。而本文用于測試的雙流體模型是在CFX-5 CFD軟件中構建。

?1.2界面跟蹤vs.雙流體概念

界面流動是指兩種或兩種以上的不混相流體被隨時間變化的尖銳界面所分離的多相流動問題。通常，當界面一側的流體是在界面上施加切向應力的氣體時，則該界面則被稱為自由表面。ITM方法最適合對這類流動進行處理，因為其相當準確地表示了接口拓撲。當界面的確切形狀無法確定或相關性較弱時，可借助雙流體方法，其中每個相需要單獨的守恒方程，并考慮兩相之間的作用力。這兩類方法都可處理界面流，但精度不同。ITM是本文介紹的單流體方法的基礎算法。

單流體方法需要求解一組具有可變材料性質和表面力的守恒方程(Lakehal等 2002)。因此，該方法比雙流體模型更準確，因為它最大限度地減少了建模假設。在單流體形式中表示的不可壓縮多相流動方程采用以下形式：

其中ρ是密度，p是壓力，μ是粘度，

。方程(2)中最后的一項表示表面張力，n表示到界面的法向量，κ為表面曲率，γ為表面張力系數(shù)，δ表示以界面為中心的光滑狄拉克函數(shù)。在這里采用的Level Set方法(Sussman et al 1994)中，互不摻混流體之間的界面用連續(xù)函數(shù)φ表示，φ表示到界面的距離，在界面上設置為零，在一側為正，在另一側為負。材料屬性、體力和表面力作為此函數(shù)的函數(shù)在局部更新，并使用平滑的Heaviside函數(shù)在界面上進行平滑：

在實踐中，Level Set函數(shù)在方程(4)的單步迭代后不再是距離界面的帶符號距離。為了恢復其在界面附近的正確分布，使用三階或五階WENO格式將re-distancing方程迭代到穩(wěn)態(tài)。

不可壓縮流動條件下不存在相變時，k相層流運動方程用平均雙流體形式表示為:

式中εk=ρkαk是表觀密度，αk為孔隙率(αL+αG=1)，實際上為當?shù)氐臍庖后w積比。假定系統(tǒng)壓力唯一，通過界面交換力對每個相交替進行耦合求解。

需要類似于(4)的拓撲方程來更新空隙率隨時間的變化，但不能與ITM的(例如Level Set和VOF)進行相似。方程(6)的RHS項反映了動量方程中相間相互作用力，即FM,G+FM,L=0。例如在散體流動中，RHS項將包括增加的質量、曳力、升力(以及固體顆粒的薩弗曼力)。對于自由表面流動，該項包括基于滑移速度定義的界面剪切力

。與前面提出的單流體方法相比，這個方程的問題在于需要確定特定系數(shù)(例如Cd)來確定動量交換項。

?

2.二維段塞流的形成

在二維平面管道中模擬二維段塞流動的形成過程，兩相流體是液態(tài)和氣態(tài)氨，與參考對比的實驗工況相同(Martin等, 2005)。計算域為一根長6.3m，內徑146.3mm，壁厚11mm的碳鋼管。管道兩端關閉，但在管道頂部設置有兩個直徑0.3m的孔是開放的，分別設置為入口和出口，流體從這兩個孔進入和離開管道。純氨最初被放置在隔熱的水平管道中，在指定深度達到熱平衡，初始深度從38.1mm到127mm不等。在頂部入口處，流速在1-10m/s之間變化。使用的計算網(wǎng)格由400×20個節(jié)點組成。兩側壁面采用無滑移邊界條件，在管道的右端出口處采用outflow邊界條件。

2.1Level Sets(LS)結果

圖1展示了管道流動中波浪的發(fā)展情況。入口流速為5m/s，氣液交界面選擇空隙率為50%。該圖顯示了在經(jīng)過不同時間步后，數(shù)值計算得到的自由表面變化情況。段塞的形成過程清楚地表明是一階模態(tài)生長的結果。結果顯示其是由壓力造成的，而非由施加的剪切造成。在第一模態(tài)(高振幅長波)的尾跡中形成的第二模態(tài)(短波)，其在第4幅圖中與一階模態(tài)合并導致管道密封。圖1展示的結果很好地說明了這種模態(tài)間的相互作用或波-波間的相互作用過程。

段塞流在距離入口6個通道高度的下游形成(第5幅圖)。隨著時間的推移，液體形成的堵塞會導致更多的液體積聚(最后3幅圖)，這似乎進一步隔離了氣體段塞。該阻塞機制不斷重復，最終沿著管道形成連續(xù)的段塞。第5幅圖中顯示出兩種新模式的生長，這隨后引起波的合并，產生第二個阻塞，段塞最終在管道的右端消散。

2.2雙流體模型結果

圖2是在入流速度為5m/s工況下，采用雙流體模型得到的自由液面演化過程。波生長機制的開始方式與Level Set模擬的結果非常相似。但與Level Set方法的模擬結果相反，在第一幅圖上已經(jīng)出現(xiàn)了各種其他的模式，在較遠的下游處的自由表面開始出現(xiàn)表面的不穩(wěn)定現(xiàn)象。這是因為在Level Set模型中，表面張力阻止了短波的形成，而在雙流體模型中，表面張力無法平衡界面剪切的作用。

雙流體模型本質上并不能預測尖銳界面。段塞在接近管道的右側時出現(xiàn)消散，如圖2的最后3幅圖所示。在這種情況下，并沒有觀察到LevelSet模擬得到的流動結構。相反，這兩種液體傾向于逐漸混合均勻。圖2中的綠色等高線反映了雙流體模型計算結果中混合帶的存在。這是雙流體方法與ITM方法相比的另一個缺點，ITM方法現(xiàn)在被證明是最適合這類流動的。當雙流體模型應用于自由表面流動時，另一個誤差來源是缺少一個嚴格的界面阻力模型

?

3.三維圓管段塞流的形成

3.1問題建立

建立一個完整的三維計算域。部分網(wǎng)格如圖3所示。管道長6.3m，管徑0.14m。采用多塊網(wǎng)格策略，用相鄰子域覆蓋域，處理器12個，MPI并行執(zhí)行，網(wǎng)格數(shù)量近36萬。采用引言中提到的One-Fluid-LES方法(LEIS方法)。具體來說，其利用Level Set技術進行界面跟蹤。湍流采用由Liovic和Lakehal(2007)擴展到界面兩相流的LES方法建模。

流入邊界條件包括固定的表面氣液速度和空隙率，條件與Martin等(2005)的實驗保持相同，管道右端保持敞開。入口處設置氣體速度為14m/s，液體速度為0.5m/s，氣液以入口中線為界。入口流入的液體溫度為-33℃，氣體溫度為20℃?；诒诿婕羟欣字Z數(shù)，在整個流域內施加初始流動擾動。

3.2整體水動力學特征

下圖為進口氣體速度為14m/s，空隙率為50%的氣液管道流動中波浪的發(fā)展情況。圖4顯示了在不同時間步后自由液面的位置變化情況。在壓力和界面剪切作用下，在距離入口約1.5倍管徑處形成一階波模態(tài)，波面具有三維特征。較小的斜坡波在下游進一步發(fā)展。二階波模態(tài)在一階波模態(tài)之后形成。在段塞形成前，二階波試圖進一步向下游的一階波合并。段塞流在距離入口6倍管徑的下游形成。與上文中的二維模擬的不同之處在于，密封前形成的短、低振幅波本質上是三維的，自由液面沿展向分布存在不均勻的特征。

值得注意的是，在段塞的下游，由于液面與氣體之間發(fā)生剪切，以及湍流結構的強烈相互作用，形成了各種坡度和波長的自由液面變化。此外，密封前形成的初始波型呈側向傾斜，并非完整的正弦波形。圖4中最后一幅圖顯示波峰在密封時的不穩(wěn)定性，導致強烈的表面起皺，正如Liovic和Lakehal(2007)在破碎波中發(fā)現(xiàn)的那樣。

隨著時間的推移，由液體形成的堵塞會導致更多的液體積聚，這進一步隔離了氣體段塞。該機制不斷重復，最終沿著管道形成連續(xù)的段塞。

?3.3表面湍流相互作用

圖5展示自由表面以及速度在x,y,z三個方向上的分量大小（云圖顏色），分別顯示在二維中心y-z平面上。圖5顯示在液體發(fā)生密封后，流場中立即形成的旋渦并脫落。三個方向上的瞬時速度分量證實了這一點，尤其是沿流向（z方向）的速度表現(xiàn)最為明顯。

圖6展示了段塞密封后自由表面的破裂過程，圖中結合展示自由表面和速度等值線。圖6的第3幅圖展示了破碎過程對分層液體自由液面的影響，該區(qū)域具有非常高的湍流強度。

圖7展示了正和負渦度水平的情況，結果表明，即使在段塞破裂之前，湍流在段塞的下游也特別活躍。在波舌的撞擊處產生的湍流非常高，最終會穿透到表面之下。這些流動結構以及湍流的產生和傳輸與破碎波中揭示的規(guī)律非常相似(Liovic & Lakehal, 2007)。

3.4段塞速度

圖8討論了段塞尾部和中心的位置隨時間的變化，顯示了段塞位置與時間的關系，并揭示了二者對時間的線性依賴關系，這與Martin等人(2005)的測量結果一致，他們在該工作中再現(xiàn)的相同流動條件下獲得了平均段塞流速度Us=9.4m/s。在該工況下，段塞速度保持恒定，與上文中二維模擬的結果相同。對于水動力學中的段塞流動，如果管道中含氣率已知，可根據(jù)氣液流速計算段塞流速度。在水平線上，段塞流內液體的平均速度近似等于混合速度，或者也可以用解析法估計(Collins et al, 1978):

其中Um代表混合速度。圖8結果表明，LES模擬預測的段塞速度(尾部和中心)與理論解一致，給出的Us=9.33 m/s，與Martin等人的實驗(Us=9.4 m/s)一致。

對于工程實際應用而言，計算結果表明，即便在使用粗網(wǎng)格的條件下，LES方法依舊能夠預測管道流動的重要流動特征。而更精細的網(wǎng)格將有助于獲得段塞流在流動自由表面上俯沖后的流動和兩相分布的詳細圖像。

3.5段塞的起源:與理論解比較

Wallis和Dobson(1973)、Taitel和Dukler(1976)以及其他很多前人的研究已經(jīng)形成了從分層流動到段塞流過渡的閾值標準。這里將本文CFD結果與Taitel和Dukler(1976)的理論進行對比，該理論解認為段塞流形成的閾值依賴于無量綱體積氣體通量和空隙率(α=AG/A)：

式中，為氣體所占面積，D為管道直徑，g為重力，A為管道橫流面積，M為質量通量。Taitel和Dukler(1976)的模型也與Martin等人(2005)提出的模型進行了比較，基于他們的測量，根據(jù)

將發(fā)生段塞的臨界體積通量與空隙率聯(lián)系起來。圖9中展示的CFD結果與理論解之間的對比包括二維層流與三維湍流的結果。

從二維結果來看，綠線對應的氣體名義速度為5m/s，確實可以預測段塞流的形成，盡管該值略低于Taitel & Dukler模型的閾值。雙流體模型的二維結果不包括在內，因為與ITM的模擬相比，很難從圖像上定義雙流體模型結果中的段塞流判斷。三維模擬結果表明，閾值與Martin等人的模型相當吻合，這實際上遠遠高于Taitel & Dukler的模型。

?

結論

本文研究氣液二維通道和三維管道流動中波浪的發(fā)展與入口空隙率和氣體雷諾數(shù)的關系。CFD結果基于VirtualFlow計算流體力學軟件，使用Level Set方法進行模擬。結果表明，該方法在段塞流動預報方面優(yōu)于雙流體方法。雙流體方法將會導致氣液界面模糊，在密封區(qū)下游產生非物理的混合流區(qū)。VirtualFlow對三維湍流的模擬提供了更好的流動細節(jié)圖像，特別是關于段塞形成和湍流界面相互作用的開始，同時對段塞區(qū)速度的預測進行了實驗和理論對比驗證。

?從工程應用角度來看，基于VirtualFlow的數(shù)值計算結果在管道石油運輸領域與使用一維集總參數(shù)模型相比可以更嚴格地接近真實解。數(shù)值結果證明，基于VirtualFlow可以使用相同的方法處理復雜地形上的管道多相流動包括含固體顆粒（散體）的復雜多相流動過程。實踐表明，如果對1.6km長的管道進行數(shù)值仿真，網(wǎng)格量約在500萬左右，預計將需要在32個CPU的Linux集群上進行1周的數(shù)值計算。??