吳 曉，王兆婷，史 瑤

（1.中國石油大學(xué)勝利學(xué)院，山東東營257061；2.中國石油大學(xué)（華東）儲運(yùn)與建筑工程學(xué)院，山東青島266580；3.中國石油新疆油田分公司工程技術(shù)研究院，新疆克拉瑪依834000）

水平管氣液兩相段塞流廣泛存在于石油、化工和核能等領(lǐng)域。段塞流流動帶來的壓力和流量波動對一些生產(chǎn)設(shè)備的設(shè)計具有重要影響，所以深入研究水平管段塞流的流動特性具有重要的學(xué)術(shù)意義和應(yīng)用價值。由于段塞流流動的復(fù)雜性，對其流動特性的研究多采用實驗平均測量的手段，其只能給出宏觀規(guī)律特征，而CFD更能細(xì)致地捕捉到瞬態(tài)流動特征、起塞機(jī)理及氣液段塞流相界面結(jié)構(gòu)特征，選取準(zhǔn)確的模型能夠保障多相流管道的安全設(shè)計和操作。

目前已有的氣液段塞流的計算模型可以分為穩(wěn)態(tài)模型和瞬態(tài)模型兩類。段塞流的穩(wěn)態(tài)模型使用較為廣泛但存在一些缺陷，使用段塞流穩(wěn)態(tài)模型前需要進(jìn)行流型判斷，段塞流穩(wěn)態(tài)模型只能給出段塞流的平均特征參數(shù)；同時Y.Taitel等[1]研究表明，穩(wěn)態(tài)模型在計算下傾管段塞流時常出現(xiàn)非物理解。

段塞流瞬態(tài)模型大致可以分為一維雙流體模型、液塞跟蹤模型和界面捕捉模型。一維雙流體模型一般都基于M.Ishii[2]推導(dǎo)的一維雙流體模型框架基礎(chǔ)上建立的。但針對不同流型建立的守恒方程和方程的本構(gòu)關(guān)系卻有很大的差異。M.Lu[3]通過實驗和模擬相結(jié)合方法對水平管段塞流的起塞、成長及合并的過程進(jìn)行研究，并對六種CFD軟件進(jìn)行對比分析，發(fā)現(xiàn)TRIOMPH下的一維雙流體模型對段塞流瞬態(tài)特征的模擬存在以下問題：只能捕捉由K‐H模型不穩(wěn)定產(chǎn)生的長波，不能考慮到界面擾動帶來的短波影響。一維雙流體模型具有較好的適用性，但存在很多缺點(diǎn)：計算范圍具有局限性，超出范圍出現(xiàn)數(shù)值發(fā)散、病態(tài)無解情況；界面阻力的本構(gòu)方程不夠成熟；能夠預(yù)測管道內(nèi)壓力波動特征，對液塞演化過程的預(yù)測不理想；一維模型不能模擬入口幾何條件，無法模擬入口擾動帶來的影響及三維下的液塞演化過程。

G.Zheng等[4]較早提出了液塞跟蹤模型，該模型是通過對每個液塞頭尾位置的Lagrangian追蹤，模擬液塞的運(yùn)動、成長和消失的演化過程，能夠較好地模擬段塞流的瞬態(tài)特征。但液塞跟蹤模型存在如下缺陷：計算時必須在入口處給定液塞長度分布，不能從根本上模擬追蹤液塞的自然演化過程；模型中氣彈移動速度與前方液塞長度的關(guān)系式是基于低氣液速度的實驗得到的，對于高氣液速的工況計算有待驗證。

最近發(fā)展起來的界面捕捉模型具有如下優(yōu)點(diǎn)：不同流型可采用相同的控制方程和本構(gòu)關(guān)系式；通過數(shù)值求解自動捕捉液塞的自然演化過程，不需引入額外擾動；能夠真實地反映出起塞機(jī)理、段塞流的瞬態(tài)流動特性及相界面結(jié)構(gòu)特征。

目前段塞流模擬多數(shù)依靠一維、二維模型，不能很好地反映實際情況下相界面結(jié)構(gòu)和瞬態(tài)流動特征，M.Lu[3]、S.C.K.De Schepper[5]、M.Abdulkadir[6]等通過 Ansys、Fluent、STAR‐CD 等商業(yè)軟件建立三維模型對段塞流進(jìn)行模擬，結(jié)果并不理想。作為開源軟件，OpenFoam擁有強(qiáng)大的底層類庫，充分利用了C++封裝、繼承、多態(tài)等特點(diǎn)，減小二次開發(fā)難度，用戶可以自由的對算法、邊界條件、求解器等進(jìn)行研究開發(fā)，根據(jù)模擬情況選擇最適合自己的模型，為CFD領(lǐng)域的研究者們提供了一個強(qiáng)大的開源平臺。

綜上所述，本文基于OpenFoam[7]開源軟件，利用界面捕捉方法建立統(tǒng)一的兩相流流動三維數(shù)值模型，無需外加擾動，對段塞流的自然起塞機(jī)理和液塞的成長、演化過程進(jìn)行深入細(xì)致的研究，并與理論、實驗相比較驗證模型的有效性。

1 兩相流數(shù)值模型

基于三維兩相流數(shù)值模擬的復(fù)雜性，對其進(jìn)行假設(shè)：氣、液相為不可壓縮流體；兩相之間無傳熱和傳質(zhì)現(xiàn)象；氣、液相與管壁之間為無滑移邊界條件。

1.1 VOF模型

VOF模型最早由C.V.Hirt等[8]提出，之后由E.Berberovi'c等[9]針對兩相流進(jìn)行改進(jìn)，通過求解單獨(dú)的動量方程和關(guān)于流體體積函數(shù)的微分方程，間接實現(xiàn)單元體的界面追蹤，基本原理如圖1所示。

圖1 VOF界面追蹤原理Fig.1 Schematic diagram of VOF interface tracking

函數(shù)α守恒方程為：

體積分?jǐn)?shù)守恒對于流動來說非常重要，微小的體積分?jǐn)?shù)誤差都會導(dǎo)致界面位置的較大變化，為了減少氣液相界面之間的模糊邊界，對α的傳輸方程加入人工壓縮項，改進(jìn)后方程為：

式中，u為瞬時速度，ur為壓縮界面的速度場，即氣、液相的相對速度，限制了α的求解區(qū)間為0～1，即人工壓縮項僅對相界面有作用，對氣液相內(nèi)部不產(chǎn)生影響。局部流體物理性質(zhì)，如ρ、μ通過體積分?jǐn)?shù)對兩種流體加權(quán)得到。

1.2 控制方程

模型采用單流體法，氣、液共用一套質(zhì)量、動量方程，與體積分?jǐn)?shù)傳輸方程組成封閉方程，控制方程為：

式中：P為壓強(qiáng)；g為重力加速度；ρ為密度；δ為表面張力系數(shù)；κ為交界面曲率；τ為黏性力，本文采用牛頓流體黏度模型。動量方程最后一項為界面張力源項，采用J.U.Brackbill[10]連續(xù)界面力模型(CFS)。

1.3 幾何模型、數(shù)值方法和邊界初始條件

幾何模型采用笛卡爾坐標(biāo)系對內(nèi)徑50 mm，長16 m的三維水平圓管進(jìn)行模擬，其中壁厚極薄忽略不計，與顧漢洋[11]實驗條件相一致，幾何模型如圖2所示。

圖2 幾何模型Fig.2 Geometric model

圖2 中，入口分層高度比為0.5，空氣、水分別通過圖中的藍(lán)、紅色入口截面恒流速進(jìn)入水平管道；出口壓力控制恒為大氣壓；對可編譯的算例文件進(jìn)行α、u、初場壓強(qiáng)以及空氣、水物性的設(shè)置，參考王鑫[12]的液塞追蹤模型，定義單位長度的相界面張力為0.07 N/m，湍流模型的k、ε設(shè)置根據(jù)不同工況通過下式計算得到。湍流模型的選擇將在后面進(jìn)行詳細(xì)討論。

式中，I為湍流強(qiáng)度，l為湍流尺度；L為水力直徑；Cμ值為0.09；U為平均速度。

在模型中，時間的離散格式采用一階隱式歐拉格式；梯度格式采用有限元高斯線性插值；對流項格式采取高斯迎風(fēng)和高斯線性格式，這些是二階插值格式，可以保證計算的精度和穩(wěn)定性；擴(kuò)散項格式采用二階高斯守恒格式；表面插值采用中心差分格式；壓力的計算采用PCG預(yù)處理共軛梯度對稱矩陣求解器；速度、湍動能、耗散率的計算采用PBiCG預(yù)處理共軛梯度非對稱矩陣求解器；體積分?jǐn)?shù)采用光順求解器下的高斯賽德爾迭代計算。

本文模型采用了OpenFoam獨(dú)有的PIMPLE算法，即在每個時間步內(nèi)用SIMPLE穩(wěn)態(tài)算法求解，將每個時間步內(nèi)看成穩(wěn)態(tài)流動，而時間步長的步進(jìn)用PISO算法來完成。外循環(huán)的壓力修正次數(shù)設(shè)為2，內(nèi)循環(huán)體積分?jǐn)?shù)求解次數(shù)設(shè)為4，相界面人工壓縮項設(shè)為保守壓縮。

1.4 湍流模型

DNS、LES湍流模型計算效率低，兩相流模擬計算中常用的湍流模型有標(biāo)準(zhǔn)k‐ε模型、RNGk‐ε模型、可實現(xiàn)k‐ε模型、k‐ωSST 模型。以Usg=5 m/s、Usl=1 m/s為例(Usg為氣相表觀速度，Usl為液相表觀速度），在相同的網(wǎng)格數(shù)和時間步長下對這四種湍流模型進(jìn)行模擬，得到起塞距離、液塞頻率與顧漢洋[12]實驗對比，如圖3、4所示。從圖3、4可以看出，可實現(xiàn)k‐ε模型下的算例與顧漢洋[12]的實驗結(jié)果最為接近，且符合程度較好，故選可實現(xiàn)k‐ε模型和壁面函數(shù)法相結(jié)合作為湍流模型。

圖3 界面起塞位置模擬與實驗值比較Fig.3 Comparision between simulated and experimental value of forming slug position

圖4 模擬與實驗的液塞頻率沿流動方向變化比較Fig.4 Comparision between simulated and experimental value of slug frequency along the fow direction

1.5 網(wǎng)格、時間步長獨(dú)立性分析

為了在計算中捕捉氣液界面的瞬態(tài)特征，保證Courant數(shù)小于0.5，空間步長Δx為一個定值，根據(jù)每個計算時層得到的最大相速度Umax自動調(diào)整時間步長Δt。

以Usg=5 m/s、Usl=1 m/s為例，在相同的空間步 長 下對 Courant數(shù) 分別 為 0.125、0.250、0.325、0.500進(jìn)行模擬，如圖5所示。從圖5中可以看出，時間步長對模擬結(jié)果影響不大，為提高計算效率，故選擇Courant數(shù)為0.500時的時間步長。

采用邊界層加密的六面體網(wǎng)格對模型進(jìn)行劃分，為了排除網(wǎng)格稀疏程度對模擬結(jié)果的影響，以Usg=5 m/s、Usl=1 m/s為例，分別對空間步長為0.002 0、0.002 5、0.003 0、0.004 0 m 的算例進(jìn)行模擬，結(jié)果如圖6所示。由圖6可知，網(wǎng)格越小，與實驗值越接近，考慮計算效率的影響，最終選擇空間步長為0.003 0 m。

圖5 界面起塞位置隨時間步長的變化Fig.5 Variation of forming slug position with time step

圖6 起塞距離、液塞頻率隨空間步長的變化Fig.6 Variation of forming slug position and slug frequency with space step

2 計算與模型驗證

通過上述選定的兩相流數(shù)值模型，對內(nèi)徑為50 mm的管道進(jìn)行不同氣、液速條件下的氣、水兩相流模擬。表1為模擬捕捉到的5種不同流型與實驗流型的對比。由表1可以看出，模擬的流型與實際較為一致，能很好地反映兩相流在不同氣液速下的相界面結(jié)構(gòu)特征。

根據(jù)表1對不同氣液速下的流型進(jìn)行判斷，以氣相表觀流速為橫坐標(biāo)、液相表觀流速為縱坐標(biāo)，取對數(shù)坐標(biāo)系作流型圖，與Mandhane（1974）流型圖進(jìn)行對比，如圖7所示。Mandhane流型圖以水‐空氣的實驗數(shù)據(jù)為基礎(chǔ)，適用范圍為：管徑12.7～165.1 mm，液相密度 705～1 009 kg/m3，氣相密度0.8～50.5 kg/m3，液相黏度 3～900 mPa?s，氣相黏度0.10～0.22 mPa?s，表面張力24～103 mN/m，本文的計算模型完全符合條件，且模擬結(jié)果基本與Mandhane流型圖相一致，驗證了該模型對兩相流流型模擬的可靠性。

表1 5種流型模擬與實驗對比Table 1 Comparison between simulation and experiment of five flow patterns

圖7 流型模擬結(jié)果與Mandhane流型圖比較Fig.7 Comparison between of flow pattern simulation and Mandhane flow pattern

圖 8為Usg=5 m/s、Usl=1 m/s時，模擬與顧漢洋[11]實驗的液塞頻率沿流動方向變化的比較。

圖8 模擬與實驗的液塞頻率沿流動方向演化比較Fig.8 Comparison between simulation and experiment of slug frequency along the flow direction

從圖8中可以看出，模型從定量上得到了模擬與實驗相一致的結(jié)果，真實地反映了起塞、液塞合并、液塞趨于穩(wěn)定的現(xiàn)象，證明了該模型對段塞流模擬的有效性。

3 模擬與結(jié)果分析

圖9為Usg=5 m/s、Usl=1 m/s時管道入口處的起塞過程，在較快流速的氣體帶動下，形成頻率大、振幅小的一系列界面波，氣、液之間的動量傳遞引發(fā)界面不穩(wěn)定，使波迅速增長并相互合并，最終導(dǎo)致界面波頻率變小，波長、振幅增大，這與顧漢洋[10]結(jié)論相一致。當(dāng)氣體上層流動空間變小，在入口恒流量的情況下，氣體加速通過界面波最高點(diǎn)處，由于伯努利效應(yīng)導(dǎo)致界面波上面空氣壓力減小，當(dāng)壓差的作用克服了界面張力和重力時，界面波被抬高，界面波的波峰阻塞管道時即產(chǎn)生起塞現(xiàn)象。

圖9 Usg=5 m/s、Usl=1 m/s的起塞過程Fig.9 Plugging process ofUsg=5 m/s andUsl=1 m/s

模擬中得到液塞頭部向前拋出并卷吸氣體的現(xiàn)象與實際高速攝像拍到的液塞頭部結(jié)構(gòu)相一致。圖10為某一液塞局部的速度矢量、水體積分?jǐn)?shù)和壓力分布云圖。分析其原因，液體橋塞管道后，液塞后面壓力陡升，在液塞前后壓差的作用下，液塞頂部液體被快速推出，形成液塞頭部向前拋出，在重力的作用下下落，液塞體的動能大于前方液膜區(qū)動能，如圖10（a）矢量箭頭的顏色所示。當(dāng)液塞吸收液膜中液體時，液膜中液體被加速到液塞速度，在液塞前端形成渦流區(qū)，由于渦流的存在，液塞會卷吸部分氣體，最終在液塞頭部形成一個具有渦流和高動能的氣液混合區(qū)。這與A.E.Dukler等[13]提出的最短液塞長度物理模型相一致，即液塞速度大于液膜速度，假設(shè)液塞靜止不動，相當(dāng)于前方液膜射流進(jìn)入液塞，在液膜分離點(diǎn)和回流點(diǎn)之間形成渦流，并伴有大量氣泡的卷吸，這個區(qū)域即為液塞混合區(qū)。

表 2分別為Usg=5 m/s、Usl=1 m/s和Usg=3 m/s、Usl=1 m/s時，一個發(fā)展中液塞隨時間變化的運(yùn)動過程。從表2中可以看出，氣液速差越小，液位越高，液塞長度越長，導(dǎo)致兩種工況下的液塞發(fā)展過程有所差異。氣液速度差較小時，液塞頭部拋出現(xiàn)象不明顯，液塞在氣體推動下速度均衡向前發(fā)展，吸收前面液膜，液膜液位波動不大；氣液速度差較大時，液塞有明顯的拋出現(xiàn)象，液塞通過拋出、下落的周期性運(yùn)動向前推進(jìn)吸收液膜。這是因為氣、液速度差較大時，液位偏低，發(fā)生橋塞的液量較小，在高速氣流的推動下液塞頭部向前拋出，液體與管子上部出現(xiàn)空隙，氣體加速通過，液體再次被抬起，液塞以此種方式循環(huán)向前運(yùn)動。而氣、液速度差較小時，液位較高，液塞體長，在氣體的推動下，液塞體較均勻向前運(yùn)動。

圖10 液塞局部的速度矢量、水體積分?jǐn)?shù)和壓力分布Fig.10 Distribution of the velocity vector,water volume

表2 不同氣、液速度差下某一發(fā)展中液塞隨時間的變化Table 2 Variation of a developing slug with time at diffie?rent gas?liquid speed difference

圖11為Usg=5 m/s、Usl=1 m/s時某一液塞不同截面處速度分布云圖和徑向速度曲線。圖11可知，液塞頂部液體速度大于底部液體速度，液塞頭平均速度大于液塞尾平均速度，很好解釋了液塞是通過頭部的液體卷吸和尾部的液體流出不斷向前運(yùn)動的。

圖11 液塞不同截面處速度分布Fig.11 Cloud chart of velocity distribution at different sections of liquid plug

氣體的推動使液塞上半部分速度接近氣速，而下半部分由于受液體和壁面剪切力控制，液速略大于液膜速度，因此可以說明，液塞運(yùn)動速度主要由液塞頭部和上半部分的液體速度決定，這與波的傳播方式極為相似，液塞的形成、發(fā)展和運(yùn)動過程可以看作管道入口的氣液速度差引起的擾動形成界面波，界面波進(jìn)行合并、增長，由于受封閉空間限制，最終橋塞形成液塞，液塞以一種變相波的方式向前運(yùn)動。

4 結(jié) 論

本文建立的模型能夠有效預(yù)測氣、液兩相流流型及相界面結(jié)構(gòu)特征，對于流動較為復(fù)雜的段塞流流型，能夠真實地反映出起塞、液塞的成長、演化過程以及段塞流的瞬態(tài)流動特性，從定性、定量兩個角度證明了該模型的有效性。

在氣體的帶動下，管道入口形成頻率大、振幅小的一系列界面波，界面不穩(wěn)定使波迅速增長并相互合并，最終導(dǎo)致界面波頻率變小，波長、振幅增大，當(dāng)界面波的波峰阻塞管道時即產(chǎn)生起塞現(xiàn)象。在氣體推動下，液塞頭部速度大于前方液膜速度，液塞頭向前拋出并下落，在液膜、液塞交界處形成渦旋，卷吸大量氣體形成液塞混合區(qū)。

氣液速度差較小時，液塞頭部拋出現(xiàn)象不明顯，液塞在氣體推動下速度均衡向前發(fā)展吸收前面液膜，液膜液位波動不大；氣液速度差較大時，液塞有明顯的周期拋出現(xiàn)象，液塞通過拋出、下落的周期性運(yùn)動向前推進(jìn)吸收液膜。

液塞頂部液體速度大于底部液體速度，液塞頭平均速度大于液塞尾平均速度，液塞以頭部液體卷吸和尾部液體流出不斷向前運(yùn)動。液塞上半部分速度接近氣速，而下半部分速度略大于液膜速度，液塞運(yùn)動速度主要由液塞頭部和上半部分的液體速度決定，并以一種變相波的方式向前運(yùn)動。