李巖松, 丁鼎倩, 韓 東, 劉 靜, 梁永圖

(1. 中國石油大學(北京) 城市油氣輸配北京市重點實驗室，北京 102249； 2. 美國加州大學圣地亞哥分校 機械與航空工程系，美國 圣地亞哥 92093-0411；3. 中國石化銷售有限公司華南分公司，廣州 510000；4. 清華大學 工程力學系，北京 100084)

目前，石油產品的運輸方式主要有空運、水運、公路運輸、鐵路運輸及管道運輸.其中，管道運輸具有費用低、可靠性高、自動化程度高、輸送能力大以及油品損耗率低等優(yōu)點，已成為石油產品的最主要運輸方式[1-2].在成品油管道運行過程中，管內油品所攜帶的游離水在一定溫壓和流速條件下會從油流中分離出來，聚集在管道沿線低洼位置，形成積液和固體顆粒雜質[3].當油品流經積水位置時，管內積水會嚴重影響油品的質量指標，如油品潔凈度、電導率、水分離指數等[4-7].因此，研究成品油管道內積水的高效清除方法具有十分重要的工程意義.

國內外部分學者借鑒氣藏連續(xù)攜液[8]、漿液管道輸送[9]以及鉆井液攜帶巖屑[10]的相關理論，提出了一種利用油流的剪切作用，將管道低洼處的積水攜帶出管道的“清管”方法.該方法只需通過調節(jié)管道入口壓力和流量，而不需要依靠清管器就可實現清除管內積水和固體雜質的目的，因此這種清管方式又被稱為水力清管[11-13].

臨界攜水油速是指管內積水開始進入上傾管段時的最小油速[14-16]，為了便于數值計算，Magnini等[14]將油流攜水量為積水總量5%時的管內最小油速定義為臨界攜水油速，其是制定管道水力清管方案的必需參數.國內外學者對臨界攜水油速進行了大量研究.Xu等[15-16]和Song等[17]建立了油攜水的二維仿真模型，利用ANSYS Fluent中的Volume of Fluid(VOF)模型分別研究了管道傾角、管道直徑、積水體積、油相密度對臨界攜水油速和攜水體積的影響.Magnini等[14]指出水相與管壁的潤濕性能對初始狀態(tài)下管道底部積水的存在形態(tài)影響很大，進而影響了初始狀態(tài)下管內積水的厚度.為此，Magnini等[14]基于OpenFOAM建立了三維油攜水仿真模型，通過VOF模型研究了不同潤濕性條件、油速、積水體積對積水移動特性及油水兩相流型的影響.

雖然，臨界攜水油速對研究油攜水過程具有重要的借鑒作用，但從工程應用的角度出發(fā)，若管內積水無法被完全攜帶出管道，殘留在管內的積水依然會加速管道的內腐蝕，因此工程應用中更加關注如何獲得實現管內積水全攜帶時所需的最小油速.由此可知，研究臨界完全攜積水油速的計算方法更具有工程價值.鑒于此，本文在當前油攜水相關研究成果的基礎上，深入分析了水團在上傾管道中的運動特性，建立了上傾管道臨界完全攜積水的數值模型，提出了模型的求解算法，并通過文獻數據驗證了模型及求解算法的可靠性.最后，詳細研究了上傾管道中油攜水波狀分層流動過程，給出了上傾管道臨界完全攜積水油速的數值判據.

1 完全攜水流動型態(tài)

管內臨界完全攜積水油速是水力清除管內全部積水所需的最小攜水油速，故分析完全攜水條件下油水兩相流動型態(tài)可在一定程度上反映臨界完全攜積水條件下管內的流動特性.Magnini等[14]通過數值模擬方法研究了管內不同體積積水(15、25、40 mL)在不同潤濕角(30°、90°、120°)條件下的水力攜帶過程，獲得了清除全部積水時，管內不同油相折算速度(0.18、0.20、0.25 m/s)下，水團在上傾管道中的平均速度以及平均水相高度，詳細數據可參考文獻[14].

上傾管道中油攜水流動如圖1所示.其中：D為管道直徑；vw為水團在上傾管道中的平均移動速度；lint為油水界面在管道截面上的半長度；hw為水相平均高度.根據圖1管道截面中的幾何關系可知，水相流通面積Aw與hw之間的關系為

圖1 上傾管道中油攜水流動示意圖

(1)

油水兩相的流型一般取決于油相折算速度vso和水相折算速度vsw的相對大小，水相折算速度可表示為

(2)

式中：Ap為管道的橫截面積.由式(2)即可將文獻[14]油攜水流動過程中的水相平均速度轉換為水相折算速度.

Hanafizadeh等[18]在管徑為20 mm，管長為6 m，傾角為12° 的實驗環(huán)道上進行了油水兩相流流型的觀測，并繪制了流型圖.需要說明的是，油水兩相流流型除了與油水兩相的折算速度有關外，還受到管道傾角、流體物性、管道直徑的影響，但Hanafizadeh等[18]和Magnini等[14]的研究中油水物性、管道參數基本相同，故可近似利用Hanafizadeh等[18]的流型圖來確定上傾管道內油攜水的流動型態(tài).文獻[14]中的油攜水流動工況在油水兩相流型圖中所處的位置如圖2所示.

圖2 完全攜水油速下，上傾管道中的油攜水流動型態(tài)[14，18]

由圖2可知，在Magnini等[14]的研究中，完全攜水條件下上傾管道內油攜水屬于分層流.這主要是由于管道內積水較少，水相近似以水膜的形式吸附在管壁上.在靠近管壁附近的區(qū)域，由于流動邊界層的影響，油水界面處油相的速度較小，相界面處油水之間的剪切力較小，水相在較小的剪切作用下不足以被破碎為小水滴，所以在上傾管道油攜水過程中，油水兩相呈現分層流的流型.

此外，根據Kelvin-Helmholtz油水界面穩(wěn)定性理論可知，當油相速度低于使界面失穩(wěn)的臨界油速時，水相可穩(wěn)定沉積在管道底部，而不隨油流前行[19]；而在完全攜水條件下，油相速度大于油水界面穩(wěn)定狀態(tài)下的臨界油速，故此時油水相界面必然失穩(wěn)，相界面呈現波動形態(tài).綜上所述，在完全攜水油速下，上傾管內油攜水流動屬于局部油水兩相波狀分層流.

2 數值模型及求解算法

2.1 油攜水波狀分層流的數理模型

成品油管道潔凈程度較高，管內一般僅存在少量的積水[20]，積水高度較小，在管壁邊界層的作用下，靠近管壁附近處的油流速度較小，不足以使將附著在管道底部的積水完全破碎成分散相進入油相隨之流動，而是以偏心大水滴[16]的形式在上傾管道中隨油流向管道出口緩慢遷移，如圖3所示.其中：θ為管道傾角；Lw為水膜在上傾管道中的長度.

圖3 上傾管道中油攜水波狀分層流的物理模型

油攜水分層波狀流動的控制方程包括質量守恒方程、動量方程、波狀相界面模型和積水平均高度模型.

2.1.1質量守恒方程 由于成品油管道中積水含量較少，不考慮油水間的乳化作用，即不考慮水相對油相質量流量的影響，所以在管道軸向的每個截面上，油相總質量流量保持不變.油相質量守恒方程可以表示為

(3)

2.1.2動量方程 由于管道徑向速度遠小于管道軸向速度，可忽略徑向速度的影響.此外，管內積水一般聚集在距離管道入口較遠的低洼位置，管內油流流動已處于穩(wěn)定的充分發(fā)展狀態(tài).成品油管道中流量較大，管內油相處于高雷諾數的湍流狀態(tài)，而靠近管壁的水膜由于移動速度較小，一般處于低雷諾數流動狀態(tài).大渦模擬(LES)方法適用于研究不同雷諾數下的流動模擬[21-23].鑒于此，本文采用LES方法模擬成品油管道油攜水流動過程中油相的湍流流動.處于充分發(fā)展狀態(tài)的油、水兩相軸向濾波動量方程為

(4)

LES方法中的亞格子湍流黏度為

(5)

在近管壁和相界面附近需要考慮流動邊界層對湍流的影響，亞格子湍流黏度中可用衰減函數DS來反映壁面附近湍流的衰弱.衰減函數DS為

DS=1-e-(y+/C+)

(6)

式中：y+為離開壁面或相界面的無量綱距離；C+為經驗常數，一般取C+=25.在y+較小時，DS的計算式為

(7)

離開壁面或相界面的無量綱距離y+的計算公式為

(8)

(9)

dint=|y|

(10)

式中：n為油水界面的法向量.

綜上，無量綱質量守恒方程和動量方程為

(12)

2.1.3波狀相界面模型 上傾管道油攜水流動屬于波狀分層流.在油攜水波狀分層流動過程中，界面波動會導致油水相界面處流動阻力增加，湍流強度在相界面處無法完全衰減，這意味著需要對LES方法亞格子湍流黏度模型中油水界面附近的有效距離進行修正.

Rotta[24]提出可通過在流場計算節(jié)點到油水相界面的有效距離(見式(8))中增加一個偏移距離Δy來增加波狀相界面處的湍流強度，即

dint=|y+Δy|

(13)

將波狀相界上的界面波動等效為砂粒粗糙度可有效描述界面波動對湍流流場的影響.無量綱等效Nikuradse砂粒粗糙度Ra+和無量綱偏移距離Δy+的計算公式分別如下式所示：

(14)

(15)

油水界面的等效Nikuradse砂粒粗糙度Ra的計算方法為[25]

(16)

(17)

2.1.4管內積水平均高度模型 成品油管道在臨界完全攜積水油速下，進入上傾管道中的積水以大水滴的形式存在，此時若要求解水相流通面積上的速度分布，則首先需要計算水團的平均高度.上傾管道中水團的高度取決于油水相界面的表面張力、重力、界面剪切力以及水滴與管壁間的潤濕性能等因素的復雜綜合作用.De Gennes等[26]提出當油相對水相的攜帶作用較弱時，可利用靜態(tài)條件下水團在水平管道中的鋪展高度來近似其在油流弱剪切擾動下的平均高度，計算公式為

hw=

(18)

式中：Vw為管道內積水體積；β為水團與管壁的潤濕角；Bo為邦德數，無量綱，反映了界面張力和重力的相對大?。沪覟橛退缑鎻埩?；Rw為假設水團為球形時對應的曲率半徑；Δρ為油水兩相的密度差.

在Bo?1時，當水膜進入上傾管道后，由于傾角的影響，與油水界面正交方向上的重力分量減小，所以會導致水相高度相比于水平管道有所增厚.考慮傾角影響的上傾管道中水膜高度計算模型[26]如下所示：

(19)

式中：Ch為模型系數.當管徑較大時，管壁曲率半徑較大，與平板接近，此時取Ch=2.0.

管徑對水膜高度存在影響，但對于輸油管道而言，由于管道內介質較為清潔，管內的積水體積一般較小，所形成的水膜厚度也較小，且成品油管道的管徑一般較大(可達幾百毫米)，所以水膜在管道底部鋪展開時，水膜與管壁的接觸面積相比于管道內壁面積較小，可近似忽略管徑曲率對水膜鋪展高度的影響.因此，管徑越大，水膜在管道內的鋪展狀態(tài)越接近于平面，所提出的水膜高度模型在計算大管徑內水膜高度時計算誤差越小.

表1 邊界條件

2.2 數值求解算法

2.2.1坐標變換 Charles等[27]直接在笛卡爾坐標系下利用結構化網格對管道截面進行離散，但由于管壁的非規(guī)則特性，導致需要人為重構管壁節(jié)點，為數值計算結果引入誤差.為解決多相管流存在相界面的特點，部分學者將雙極坐標引入氣液兩相流的數值計算[28-30].在雙極坐標系中，管壁以及相界面可通過直線表示以便于計算.鑒于此，本文亦采用雙極坐標系下的結構化網格離散計算區(qū)域，將非規(guī)則的管道區(qū)域(見圖4(a))轉變?yōu)橐?guī)則的有限矩形區(qū)域(見圖4(b))，并基于坐標轉換規(guī)則，離散控制方程.其中：η和ξ分別為直角坐標系x和y軸在雙極坐標系中對應的坐標軸；P1～P4分別為兩種坐標系互為對應的特征節(jié)點.

圖4 坐標變換

控制方程采用有限容積法(FVM)進行離散，其中對流項采用二階迎風格式離散，擴散項采用二階中心差分格式離散.

2.2.2動量方程壓力梯度求解方法 管道截面流場的求解屬于壓力和速度耦合問題.壓力和速度耦合問題的求解算法主要有SIMPLE系列算法、IDEAL[31]算法以及MAC[32]算法等，上述算法主要用于求解壓力場和速度場的細節(jié)信息，而本文的主要目的是計算積水進入上傾管道后，水團在油相攜帶下遷移的軸向速度.相比于軸向速度，管道中流體徑向速度很小，故為提高計算效率，假設管道同一截面上壓力相等，不考慮管道截面二次流引起的壓力變化，即不對管道截面上的壓力場進行全分辨求解.

圖5 基于正割法的壓力梯度求解算法程序框圖

(20)

式(20)為高度非線性方程，其導數計算較為復雜，該類方程可采用正割法求解.正割法迭代公式如下：

2.2.3油攜水模型的整體數值求解算法 成品油管道油攜水兩相管流是復雜的油水兩相波狀分層流過程，需要耦合求解質量守恒方程、動量方程、波狀界面模型以及水相高度模型.管道截面的速度場、壓力梯度、剪切速率等參數需要迭代求解直至全部收斂.控制方程組的求解算法如圖6所示.

圖6 油攜水模型的整體數值求解算法程序框圖

求解步驟如下：

(1) 設定管道基礎參數，如管徑、管長、管道傾角、積水和管壁的接觸角以及管道入口油相的折算速度等；

(2) 設定當前管道截面速度場和壓力梯度的初值，并基于湍流模型計算相應的流體物性以及湍流黏度；

(3) 依次耦合迭代求解速度場和湍流黏度，直至全部收斂后，執(zhí)行步驟(4)；

(4) 基于雙極坐標下的積分公式，求解當前截面的油相體積流量.若油相體積流量守恒，即油相連續(xù)性方程收斂，則計算結束；反之，調用基于正割法的壓力梯度求解算法調整壓力梯度，并返回步驟(2)，直至油相連續(xù)性方程收斂.

3 模型驗證

3.1 管內積水水膜平均高度模型驗證

本文通過與Magnini等[14]的數值計算結果進行對比，驗證管內積水水膜平均高度模型的適用性，不同積水體積以及不同油相速度下，上傾管道中的水膜平均高度模型計算值(E1)與相同條件下Magnini等[14]的模擬結果(E2)的對比分析如圖7所示.

由圖7可知，該模型的計算結果與Magnini等[14]的模擬結果基本吻合，驗證了水膜高度模型的可靠性.從圖7中還可以看出，當積水體積和油相速度一定時，水膜與管壁的接觸角越大，模型計算誤差越小.這是由于水膜與管壁的潤濕角越大，水膜與壁面的接觸面積越小，水膜高度受管道內壁曲率的影響也越小，所以模型的計算誤差也越小.上述現象說明水膜高度模型更適用于計算管內壁為疏水特性的管道內積水平均高度.

通過對比相同油速下水膜高度的計算結果(見圖7(a)～7(c))可以看出，當積水體積越大時，模型的計算結果誤差也越大.這是因為圖7中所采用的管道直徑較小(27 mm)，積水體積越大，水膜在管道底部鋪展的高度受管壁曲率影響也越大；而式(19)中假設管道壁面為平板，因此積水體積越大，管內水膜高度的計算誤差也越大.需要說明的是，工業(yè)成品油管道中沉積在管道底部的積水體積較小，且管徑較大，這為利用式(19)計算水力清管過程中水膜平均高度提供了有利因素.

對比圖7中不同油速下的模型計算誤差還可發(fā)現，在模型計算誤差相近的條件下，隨著油速的增大，水膜高度模型系數Ch越小，即在油速較高的情況下，無需調整模型系數Ch，將其取為原始公式中Ch=2.0也可獲得較好的計算結果.這是因為管內油速越大，油相對積水的剪切作用越強，水膜在油相剪切作用下被拉長，導致水膜高度越低，油水相界面越靠近管壁區(qū)域，而管壁附近的油相剪切作用反而越小，此時水膜的形成過程越接近靜態(tài)條件下水膜自由鋪展的情形，所以模型系數Ch也越接近原始值.對于成品油管道而言，管內油速一般較大(1～2 m/s)，因此，在進行成品油管道水力清管計算時可近似取Ch=2.0.

圖7 E1與E2的對比

綜上所述，利用式(19)估算成品油管道水力清管過程中水膜的平均高度是可行的.

3.2 上傾管道內油攜水波狀分層流驗證

通過對比相同參數條件下，上傾管道內油攜水流動過程中管道截面中垂線速度分布文獻[15]計算值與本文提出的油攜水波狀分層流模型(WSFM-TWLEOS)計算值，以驗證WSFM-TWLEOS的有效性.文獻[15]中所計算的管道直徑為27 mm，管道傾角為12°，管道入口油相折算速度vso=0.03 m/s，其他參數可參考文獻[15].文獻計算結果(E3)與WSFM-TWLEOS計算結果(E4)對比如圖8所示.

由圖8可以看出，所提WSFM-TWLEOS計算值與文獻計算值基本吻合，特別是管內逆流水團的移動速度與文獻計算值近似完全一致.但油相速度的計算值與文獻值存在一定誤差，這主要是因為文獻[15]采用三維非穩(wěn)態(tài)模型進行全流場模擬，而本文在求解油攜水流動過程中忽略了管道徑向的油相速度，且未考慮水團對油相流動影響，所以WSFM-TWLEOS存在一定的計算誤差.但總體而言，WSFM-TWLEOS的計算精度基本滿足工程需求.需要特別指出，所建立的模型針對水團運動速度的計算精度較高.上傾管道中臨界完全攜積水油速的確定主要是通過水相的速度分布特征來判斷的，由此可知，基于WSFM-TWLEOS來計算上傾管道內臨界完全攜積水油速是可靠的.

圖8 上傾管道內油攜水流動過程驗證

4 臨界完全攜積水油速

4.1 基礎參數

分別計算了管道入口處油相的折算速度為vso=0.1、0.5、1.0、1.5、2.0 m/s時管道橫截面上的油水兩相的速度分布.在本節(jié)的數值研究中，管徑為100 mm，管道傾角為12°，油水界面張力為18.33 N/m.如前所述，成品管道水力清管中一般利用柴油批次攜帶管內積水，故將研究柴油攜帶水團的流動過程.柴油的密度為856 kg/m3，黏度為3.43 mPa·s，水相的密度為988.88 kg/m3，水相黏度為0.597 mPa·s.假設管壁為疏水性壁面，水團與管壁間的接觸角為120°.

4.2 低油速下上傾管中水團流動特性

低油相折算速度(vso=0.1 m/s)時，上傾管道截面上的速度場分布如圖9所示.由圖9可以看出，相比于水平管道內油水兩相流動的油相速度分布，上傾管道截面上油相速度最大值位置向相界面方向偏移(見圖9(b))，這是由于油相在經過彎管進入上傾管道后，油流受到的重力影響增大，油相速度最大值在重力的作用下向管道底部方向偏移.上述油相在上傾管道中的速度分布特性有利于管道底部積水的攜帶，上傾管道中油相速度最大值向油水界面位置偏移后，油水界面上的剪切力增大，油相對水相的剪切攜帶作用增強.

圖9 vso=0.1 m/s時管道截面速度分布

從圖9(b)中還可以看出，當油相速度較小時(vso=0.1 m/s)，水相的速度與上游來油的流動方向相反，水相區(qū)域內各點的速度均為負值，這意味著水團重新流向管道低洼處，此時無法實現管道內積水的水力驅除.這是由于油相速度較小，油相在相界面上對水相的剪切力較小，不足以克服上傾管道中水團的重力.水相速度的最大絕對值出現在油水相界面與管道底部的中間位置，即水相流動區(qū)域的中心位置，管道中垂線水相速度呈現拋物線分布，這與Song等[5]的研究結果一致，進一步驗證了所提WSFM-TWLEOS的正確性.

根據上傾管道內水相的速度場可知，水膜的平均速度為

(21)

水相的雷諾數為

(22)

式中：ρw為水相密度；Dw為水相的水力直徑；νw為水相的運動黏性系數.以水力光滑區(qū)計算可知，在數值算例參數條件下，水團流動過程中由流動引起的摩阻壓降和重力引起的壓降分別為

(23)

(24)

由式(26)和(27)可知，水相的流動過程主要受重力壓降的影響，摩阻壓降相比于重力壓降可忽略，這表明水力清除積水過程中油流的剪切作用主要用于克服重力壓降.

綜上所述，較小油速下，上傾管道中的水膜在重力作用下會重新回流至管道低洼處，故要實現積水的完全水力驅除，需要不斷調整油速，直至水膜流動區(qū)域內水相速度均大于0.

4.3 上傾管內臨界完全攜積水油速數值判據

上傾管道臨界完全攜積水流速的物理意義是指積水全部進入上傾管道，并被油流全部攜帶出管道時所需的最小油速.上述水相速度分布特征即為臨界完全攜積水油速的數值確定方法.

不同油相折算速度(vso=0.25～2.0 m/s)下，上傾管道截面的速度分布如圖10和11所示.由圖10可知，不同油相折算速度下管道截面上油相速度分布特性相似，而水相速度差異較大.從圖11(a)～11(c)可知，靠近相界面處的水相部分在油相的攜帶作用下向下游流動，而靠近管壁底部的水相部分由于油相剪切不足以克服重力作用，致使部分水相回流至管道底部.隨著油相速度的增大，油相在油水相界面上的剪切力增大，油相對水相的攜帶作用增強，隨油流向下游流動的積水體積逐漸增多.

圖11 不同油速下，上傾管道截面中垂線速度分布

對比圖11(b)和11(d)還可以發(fā)現，油相速度越大，水相速度最小值(負值)位置逐漸向管壁偏移，當水相速度的最小值出現在管壁位置時，此時管內油相速度即為臨界完全攜積水油速，上述研究結果即上傾管道中油流取得臨界完全攜積水油速的數值判據.在本文的數值算例中，臨界完全攜積水油速為0.5 m/s.從圖11(e)和11(f)中還可以看出，在臨界完全攜積水流動狀態(tài)下，水相在管道中垂線上的速度從管壁到相界面近似呈線性增大，并在相界面處達到最大.

5 結論

本文首先分析了完全攜水油速下管內油攜水的流動型態(tài)，建立了一種高效的油攜水分層波狀流動的擬三維準穩(wěn)態(tài)數值模型，首次給出了臨界完全攜積水油速的數值判據.對上傾管道內油攜水流動過程進行了數值研究，得出如下主要結論.

(1) 在臨界完全攜積水油速下，上傾管內油攜水流動屬于油水兩相波狀分層流.

(2) 相比于水平管道內油水兩相流動的油相速度分布，在重力作用下，上傾管道截面上油相速度最大值位置向相界面方向偏移，水團在流動過程中主要受到重力壓降的影響，摩阻壓降相比于重力壓降可忽略.

(3) 管內油速較低時，靠近相界面處的水相部分在油相的攜帶作用下向下游流動，而靠近管壁底部的水相部分由于剪切作用不足以克服重力影響，水相在重力作用下回流至管道底部，隨著油相速度的增大，油相在油水相界面上的剪切力增大，對水相的攜帶作用增強，積水被攜帶出管道的體積增大.

(4) 油相速度越大，水相速度最小值(負值)位置逐漸向管壁偏移，當水相速度的最小值首次出現在管壁位置時，對應的管內油相速度為臨界完全攜積水油速.