基于一維氣體單流體模型的天然氣-凝析液管道流體流動預測

張璐瑤 路宏 陳宏舉， 周良勝 李鵬程 王瑋 宮敬

1中國石油大學（北京）城市油氣輸配技術北京市重點實驗室·油氣管道輸送安全國家工程實驗室

2中海石油（中國）有限公司北京研究中心

隨著油氣勘探開發從陸地走向深海、沙漠等地處偏遠及環境復雜的地區，惡劣的生產環境促使油氣混輸技術應用更加廣泛。穩定并準確地獲取油氣管道相關運行參數（壓力、流速、溫度等），對油氣集輸管道的設計及安全運行管理具有重要意義。

國內外學者對氣液兩相流進行了大量的研究[1-8]，但是氣液兩相流的數值模擬仍不成熟。對于存在相變的天然氣-凝析液管道，因其涉及相間質量傳遞、動量和能量傳遞，通常采用均相流模型[9]、漂移流模型[10]和雙流體模型[11]進行描述和流動預測。均相流的簡化過于理想，計算結果與實際相差較大。漂移流模型中含有相關式，導致適用范圍較窄。雙流體模型是在充分考慮兩相流各相流動主要特征的基礎上建立的，相對更加穩健且使用廣泛，但學者指出采用雙流體模型時模型的雙曲性[12]非常重要，模型不適定將產生不符合流動規律的預測結果。然而，雙流體模型并不是無條件雙曲的，為保證方程的雙曲性，一些學者在控制方程中增加壓力松弛輸運方程[13]、虛擬質量力[14]、界面壓差[15]等，但因沒有明確的物理意義而存在爭議。

考慮到上述情況，針對低持液率天然氣-凝析液管道內液體和氣體的自身流動特點，在一維瞬態雙流體模型基礎上，建立了一維瞬態氣體單流體簡化模型。該模型可以避免計算過程中由于相的出現與消失而帶來的計算不穩定問題，并進一步補充封閉關系以及熱力計算，討論了一維瞬態非等溫可壓縮氣體單流體模型用于低持液率天然氣-凝析液管道內流體流動預測的適用性。

1 模型理論

1.1 控制方程

流體流動時假設忽略各參數在截面分布的不均勻性[16]?；谫|量守恒定律、動量守恒以及能量守恒定律，對截面上的流動參數進行積分，得到一維瞬態非等溫可壓縮氣體單流體模型[17]。

連續性方程：

動量方程：

溫度方程[6]：

定壓比熱和JT系數分別定義為[6]：

式中：ρg為氣相的密度，kg/m3；ug為氣相的速度，m/s；p為壓力，Pa；g為重力加速度，m/s2；θ為管道傾角；τg為剪切應力，N/m2；Sg為天然氣潤濕周長（即為管道截面周長），m；A為橫截面積，m2；cpg為定壓比熱容，J/(kg·K)；αhg為JT效應系數；U為流體與外界總換熱系數，W/(m2·K)；D為管徑，m；Te為環境溫度，K。

1.2 封閉關系

剪切應力由摩阻系數與動能的乘積表示。

氣體與管壁的摩阻系數根據ISSA[1]推薦，層流時采用Hagen-Poiseulle關系式，湍流時采用Taitel-Dukler[18]關系式。

2 數值求解

2.1 數值離散

在節點設置上采用內節點法，采用交錯網格（圖1），實線表示主網格，虛線表示交錯網格。在交錯網格中，標量設置在主網格上，如壓力、密度、溫度；矢量設置在交錯網格上，如速度[19]。對流項的離散采用一階迎風格式，時間項的離散采用一階隱式格式。

圖1 交錯網格示意圖Fig.1 Schematic diagram of staggered grid

將動量方程在速度控制容積上積分得到

將溫度方程在主控制容積上積分得到

公式中上標n和n+1 分別表示當前時層和下一時層。采用SIMPLE 壓力修正算法[20]，基本過程為：假設初始壓力為p，由動量方程可求得相應的速度u。這時的速度并不能滿足連續性方程，因此需要修正，設修正量為u′。由動量方程可知，若速度變化，壓力也會隨之改變，因此設對應于速度修正量u′的壓力修正量為p′，則新的速度與壓力為

將修正前后的速度與壓力代入動量方程中相減，可推導出速度修正量與壓力修正量之間的關系。

考慮氣體壓縮性，根據狀態方程推導出密度和壓力的關系式。

對于壓力方程的構造，將式（15）、（16）和（17）代入式（1）連續性方程中整理后可得如下形式

2.2 CFL 條件

為了保證數值計算的穩定性，在確定網格步長后，通常采用CFL條件來確定時間步長[21]，其描述了流體流動的距離在一個時步內不得超過一個網格的長度。設流體流速為u，網格步長為Δx，則時間步長Δt可表示為

其中CFL≤1。隨著網格的細化，時間步長會相應減小，故在模擬中，需檢查所有網格的CFL條件，選取最小的時間步長，作為當前的時間步長。

2.3 收斂判據

壓力修正算法的基本思想是速度與壓力相互修正，最終達到同時滿足動量方程和質量守恒方程的目的。而壓力修正方程中的源項b剛好是質量守恒方程的離散形式。如果上一迭代層次的流場滿足質量守恒，則b=0。因此b反映了每一個控制容積是否滿足質量守恒，可以作為一個層次迭代是否收斂的一個指標[19]，即

允許誤差數值ε通常取10-3～10-6，本文取ε為10-5。

2.4 算法步驟

假定管道內流體組分是恒定的，但流體特性（如黏度、密度、氣液相體積分數等）會隨溫度和壓力變化。在此假設下，首先計算出一定壓力和溫度范圍下流體特性參數表，計算過程中每一時步下參數根據這個表插值更新。當管道內存在液相時，則利用氣液相混合物性參數（密度、黏度等）來代替氣體單流體的物性參數，并將氣液混合物質量流量等效為氣相單流體質量流量。

算法步驟如下：

（1）由初始或上一步的壓力pγ、溫度Tγ更新流體參數，計算動量方程的系數與源項。

（3）求解壓力修正方程（19）得到壓力修正量p′，并根據式（17）得到速度修正量。

（4）根據式（15）、（16）得到修正后的壓力pγ+1和速度。

（6）判斷是否收斂，如果不收斂，將當前變量值作為下一層次迭代的初始值返回第一步；如果收斂，利用修正后的壓力、速度和密度求解溫度方程（3），得到Tγ+1，然后進入下一時步計算。

3 算例分析

3.1 基礎數據

某模擬海底管道總長77.5 km，規格Φ558.8 mm×23.8 mm，其地形數據及管道周圍環境溫度見表1，天然氣組分見表2。該算例邊界條件：入口定質量流量為119.98 kg/s，入口溫度為34 ℃，管壁絕對粗糙度為0.02 mm，管道總換熱系數為45 W/(m2·K)，出口定壓為10 MPa。

表1 模擬海底管道高程及環境溫度Tab.1 Elevation and ambient temperature of simulated submarine pipeline

表2 模擬海底管道天然氣組分Tab.2 Natural gas composition of simulated submarine pipeline 摩爾分數/%

3.2 計算結果

利用一維瞬態非等溫可壓縮氣體單流體模型以時間推進的方式模擬穩定流動，此過程中每一時間步下的模擬均收斂。對于最終穩定的流場，對比了本研究與OLGA[22]對模擬管道相同運行工況下五個組分（五個組分計算得到管道中最大持液率分別為0，1%，5%，10%，20%）的流動參數計算結果（表3）。計算所得氣體流速、壓力、溫度沿線分布對比見圖2～圖6。

圖2 本研究與OLGA軟件計算0持液率組分沿線氣速、壓力、溫度變化Fig.2 Gas velocity，pressure，and temperature changes along the line calculated by this research and the OLGA software at 0 liquid holdup

圖3 本研究與OLGA軟件計算1%持液率組分沿線氣速、壓力、溫度變化Fig.3 Gas velocity，pressure，and temperature changes along the line calculated by this research and the OLGA software at 1%liquid holdup

圖4 本研究與OLGA軟件計算5%持液率組分沿線氣速、壓力、溫度變化Fig.4 Gas velocity，pressure，and temperature changes along the line calculated by this research and the OLGA software at 5%liquid holdup

圖5 本研究與OLGA軟件計算10%持液率組分沿線氣速、壓力、溫度變化Fig.5 Gas velocity，pressure，and temperature changes along the line calculated by this research and the OLGA software at 10%liquid holdup

圖6 本研究與OLGA軟件計算20%持液率組分沿線氣速、壓力、溫度變化Fig.6 Gas velocity，pressure，and temperature changes along the line calculated by this research and the OLGA software at 20%liquid holdup

表3 五種持液率組分在不同計算方法下的水力熱力參數計算結果Tab.3 Results of hydraulic and thermal parameters of five liquid holdup components under different calculation methods

由對比結果可知：在五個算例中，本研究計算的氣速、溫度與OLGA結果均吻合較好。隨著持液率增加，氣速逐漸降低。因管道總換熱系數較大，入口環境溫度較低，管內流體溫降主要出現在入口0～15 km位置，隨后伴隨周圍環境溫度的升高，管內流體溫度逐漸回升；在立管段由于壓力降低帶來的節流效應將導致流體溫度有明顯下降。由圖7可以看出，在不考慮JT 效應時，立管段的流體溫度沒有下降，與實際不符，而考慮JT 效應時立管段的流體溫度出現明顯下降，隨著持液率的增加，本文計算出的終點溫度逐漸下降。由于液相密度明顯大于氣相，故管道的入口壓力隨持液率的增加而增大。本研究計算的壓力結果與OLGA結果趨勢保持一致，但數值上略有差異。當持液率達到10%時，本研究壓力計算結果與OLGA 的相對誤差為4.65%，驗證了采用一維瞬態可壓縮氣體單流體模型預測低持液率天然氣-凝析液管道內流體的流動特性是可行的。本研究算例中，當持液率達到20%時，管道入口壓力預測的相對偏差進一步增大。

圖7 0持液率下不考慮JT效應與考慮JT效應的溫度結果對比Fig.7 Comparison of temperature results without and with considering the JT effect at 0 liquid holdup

4 結論

本文采用SIMPLE 壓力修正算法，開展了一維瞬態非等溫可壓縮氣體單流體模型對低持液率天然氣-凝析液管道入口壓力、流速及溫度的預測分析。通過某模擬海管五個不同持液率（0，1%，5%，10%，20%）算例的對比分析，驗證了一維瞬態非等溫可壓縮氣體單流體模型對低持液率天然氣-凝析液管道內流體的流動預測具有適用性，可指導低持液率天然氣-凝析液管道的工程計算。