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

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

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

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

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

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

1 模型理論

1.1 控制方程

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

連續(xù)性方程：

動量方程：

溫度方程[6]：

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

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

1.2 封閉關(guān)系

剪切應(yīng)力由摩阻系數(shù)與動能的乘積表示。

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

2 數(shù)值求解

2.1 數(shù)值離散

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

圖1 交錯網(wǎng)格示意圖Fig.1 Schematic diagram of staggered grid

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

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

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

將修正前后的速度與壓力代入動量方程中相減，可推導(dǎo)出速度修正量與壓力修正量之間的關(guān)系。

考慮氣體壓縮性，根據(jù)狀態(tài)方程推導(dǎo)出密度和壓力的關(guān)系式。

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

2.2 CFL 條件

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

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

2.3 收斂判據(jù)

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

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

2.4 算法步驟

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

算法步驟如下：

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

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

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

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

3 算例分析

3.1 基礎(chǔ)數(shù)據(jù)

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

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

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

3.2 計算結(jié)果

利用一維瞬態(tài)非等溫可壓縮氣體單流體模型以時間推進的方式模擬穩(wěn)定流動，此過程中每一時間步下的模擬均收斂。對于最終穩(wěn)定的流場，對比了本研究與OLGA[22]對模擬管道相同運行工況下五個組分（五個組分計算得到管道中最大持液率分別為0，1%，5%，10%，20%）的流動參數(shù)計算結(jié)果（表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 五種持液率組分在不同計算方法下的水力熱力參數(shù)計算結(jié)果Tab.3 Results of hydraulic and thermal parameters of five liquid holdup components under different calculation methods

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

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

4 結(jié)論

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