長輸管道順序輸送混油高效數值模擬方法研究

王軍防 陳雪嬌 沈允 陳宇杰

1中國石化管道儲運有限公司科技研發中心

2北京航天試驗技術研究所

3中國石油北京油氣調控中心

4中國石油大學（北京）

隨著社會發展，我國對石油資源的需求量不斷增加，尤其是沿海發達地區，但我國的石化資源主要分布在西部和東北地區，資源產地與經濟發達地區存在空間的差異。因此，管道輸送成為油品輸送的主要方式，其具有效率高、運輸能力大、連續性高、費用低等優越性[1]。但在成品油順序輸送的過程中，受到油品輸送順序、流速、地形、管道尺寸等多種因素的影響，使得相鄰油品之間相互摻混，產生大量混油，造成極大經濟損失。混油規律復雜，目前的技術很難對混油實現精確切割，而精確切割的關鍵是對混油段的跟蹤[2]。

目前，混油研究方法主要包括實驗[3-6]和數值模擬[7-9]。實驗方法所需成本高，且采集數據困難，難以有效探究混油變化規律；數值模擬方法則成本低、操作簡單，可運用于求解復雜流動問題。數值模擬方法廣泛應用于成品油順序輸送管道混油濃度分布以及混油量的計算中，常用的混油數值模擬模型包括三維模型、二維模型和一維模型。三維模型計算精確，但計算量大，只適用于很短的管道[10]；一維模型簡單，計算速度快，但計算精度較低，適用于雷諾數很大的情況[11]，是目前現場最常用的模型；二維模型可兼顧計算效率和精度，適用于較長管道的混油模擬[12]。以上三種模型的計算量均隨著管道長度的增加不斷增大，若采用常見的二維對流擴散方程求解一條100 km 的順序輸送管道混油量，所需的計算時間甚至長達一個月。低效的計算方法限制了混油數值模擬研究的發展，并且無法滿足現場對混油量快速預測的需求。因此，亟待高效方法實現混油的快速求解，以滿足長距離管道混油的研究與預測。

本文從控制方程求解方法和計算區域兩個方面開展研究來實現高效的混油數值模擬。數值解法方面，常用的對流擴散方程求解方法包括直接求解和迭代求解。當求解區域很大時，直接解法的存儲量和計算量均很大，而迭代解法的計算量和所需的存儲空間均比直接法要少，更適合于長距離管道混油研究。常見的迭代求解方法包括雅克比迭代法、高斯-賽德爾迭代法、多重網格方法和預條件共軛梯度法，可通過選用高效的求解方法提高計算效率。計算區域方面，考慮到混油界面前方和后方均是純油品，可跟蹤混油區域，從減少計算區域角度入手提高計算效率。本文針對較精確的二維對流擴散模型，為滿足混油高效、精確的數值模擬需求，進行以下方面研究：①對比一維和二維對流擴散模型求解結果和求解速度的差異；②基于二維模型，在不同管長條件下，對比高斯-賽德爾迭代法、多重網格方法和預條件共軛梯度法求解效率，選取用于混油數值模擬的高效求解方法；③基于高效求解方法，采用追蹤方法跟蹤混油段起始界面，探究混油段追蹤方法在混油高效模擬中的適用性；④進行一維和二維模型的現場實例對比驗證，探究二維模型能否實現快速、準確的混油數值模擬。

1 模型的建立

1.1 物理模型

直接對實際的三維管道進行混油模擬，計算量巨大。若忽略重力對混油的影響，則順序輸送管道混油的濃度分布可用一個過管道中心軸線的二維平面來表示，三維問題就可簡化為平面二維問題進行處理（圖1）。而二維平面上的成品油濃度上下對稱分布，計算區域進一步簡化為x軸上半部分，如果管道徑向的流速和油品濃度分布均勻，則可近似化簡成一維問題進行處理。

圖1 管道簡化示意圖Fig.1 Schematic diagram of simplified pipeline

1.2 數學模型

1.2.1 一維模型

正常輸送情況下，物性不同的兩種油品之間的摻混過程主要包括對流擴散和分子擴散兩部分。其中對流擴散主要是由湍流脈動引起的，而分子擴散是由徑向濃度梯度引起的。當雷諾數足夠大時，對流擴散遠大于分子擴散，湍流脈動使得管道界面流速和混油濃度分布趨于均勻，此時引入有效擴散系數，借助費克擴散定律和質量守恒定律可得一維簡化模型。

式中：c為前行油品截面處油品平均濃度（體積分數，下同）；um為油品的平均流速，m/s；K為有效擴散系數，可采用式（2）計算，m2/s；t為時間，s；x為軸向坐標位置，m。

式中：ρ為油品密度，kg/m3；τw為壁面剪切應力[13]，N/m2；γ˙w,t為壁面剪應變，s-1；μ為油品的動力黏度，Pa·s；υ為油品的流速，m/s；Re為雷諾數；d為管道直徑，m。

1.2.2 二維模型

成品油順序輸送的正常工況一般為圓管充分發展的恒定湍流，主要包括層流底層、緩沖層和湍流核心區。二維模型充分考慮徑向不同流態區域流速和混油擴散系數的差異，針對不同的區域分別進行處理。在圓柱坐標下，針對圖1 中2πr、dr、dx的微元體，同樣借助費克擴散定律和質量守恒定律建立二維對流擴散方程[14]。

式中：u為管道截面流速分布，m/s；r為到管道中心的距離，m；c為前行油品濃度；D為擴散系數，m2/s。

不同流態區域的速度分布和擴散系數采用經驗公式計算[12]。

1.2.3 邊界條件和初始條件

管道中充滿前行油品，從0 時刻開始，后行油品進入管道，二維混油數值模擬在管道壁面和中心處沿徑向的濃度梯度為0，邊界條件見式（6），一維混油數值模擬邊界條件見式（7）。值得注意的是，對于二維模型，壁面處采用無滑移邊界條件，無法有效考慮壁面對油品的黏附作用。因此，無法考慮油品輸送順序對混油帶來的影響。

式中：c0為管道起點前行油品初始濃度；L為管道全長，m。

1.3 網格劃分

一維模擬直接采用一維均分網格，二維模擬管壁處采用無滑移邊界條件，在壁面的作用下，存在較長的混油尾跡，對混油的產生和發展具有重要的影響。為了精確捕捉近壁面區的信息，需要在近壁面區采用密集的網格，而研究表明，湍流核心區可以采用尺度較大的網格而計算精度并不會明顯降低[15]。因此在壁面法向采用式（8）和（9）進行網格劃分。

式中：yj為y方向第j個節點的位置；N為y方向總的節點個數；a為網格變化的調節參數；?j為中間變量。

1.4 方程求解

求解混油濃度的控制方程式（1）屬于對流占優的擴散方程，一般采用兩步法[16]進行求解，見式（10）、式（11）。其基本思路是：將對流擴散方程分解成一個純對流方程和一個純擴散方程，在前半個時步采用特征線法求解純對流方程，在后半時步采用隱式有限容積方法求解擴散方程，如此循環直到求解完成。值得注意的是，在進行純對流方程求解時，需采用高階格式，防止假擴散帶來的數值誤差[16]。

特征線法直接求解純對流方程速度較快，而純擴散方程則一般采用迭代法求解，速度較慢。目前，主要采用高斯賽德爾迭代法、多重網格方法和預條件共軛梯度法求解純擴散方程。

2 模擬結果

為高效準確地進行混油數值研究，首先對比一維和二維模型的求解結果和求解效率，其次從高效求解方法和減少計算區域兩個角度分別提升二維混油模型計算效率，最后將高效求解方法運用到二維模型中，進行現場實例驗證分析，并和一維模型計算效率相對比。

2.1 網格無關解驗證

在進行具體的混油數值模擬研究之前需確定一套合適的網格，為后續工作的開展提供支持。本文以規格為355 mm×6.4 mm、長度為1 km 的管道為例進行網格無關解驗證?；煊烷L度定義為前行油品濃度為1%～99%的管段長度。時間步長（一維模型中umax=um，本文所有時間步長都按照該方法計算），當相鄰兩套網格計算得到的混油長度之差小于3%時，即認為得到網格無關解。二維模型網格信息見表1，一維模型網格信息見表2，油品物性以呼包鄂管道為參考，具體參數見表3（本文所涉及的計算均采用該油品的物性參數）。

表1 二維模型網格無關解驗證的網格信息Tab.1 Grid information for grid independent solution verification of 2D model

表2 一維模型網格無關解驗證的網格信息Tab.2 Grid information for grid independent solution verification of 1D model

表3 呼包鄂管道油品物性Tab.3 Physical properties of oil products in Hu-Bao-E Pipeline

在表1 和表2 中的三種網格條件下，二維模型和一維模型在終點處混油分布如圖2 所示。兩種模型在三套網格下的計算結果近似，其中二維模型1 000×20 的網格和2 000×40 的網格所計算的混油長度為57.5 m 和57.75 m，相對差異為0.43%；一維模型1 000 個網格和2 000 個網格的計算的混油長度為58.06 m 和59.03 m，結果相差1.64%。所以可選擇1 000×20 的網格作為二維模型的計算網格，1 000 的網格作為一維模型的計算網格，基于相應的網格比例開展具體混油數值模擬研究。

圖2 一維和二維網格無關解驗證Fig.2 Verification of 1D and 2D models gird independent solutions

2.2 一維模型和二維模型比較

目前常用的混油理論模型包括一維對流擴散模型和二維對流擴散模型。本文以3 km 管道為例，對比兩種模型的計算結果和效率。一維模型和二維模型在終點處的濃度分布如圖3 所示，3 km 處前行油品濃度均為99%，對應的混油長度為103 m 和101.5 m，相對差異為1.46%。所以兩種模型從計算結果角度來看，差異不大。但從計算角度來看，一維模型的計算速度遠快于二維模型。盡管如此，一維模型只能使用于高雷諾數，且一維模型無法描述混油濃度在管道中分布的情況，這無疑限制了一維模型的使用。而二維模型與一維模型相反，有很強的適用性，可描述混油濃度分布情況（圖4），但計算效率低。如果能大幅度提高二維模型計算效率，則可為高效精確的混油數值模擬研究提供參考。

2.3 混油高效求解方法探究

圖3 一維和二維模型計算結果對比Fig.3 Comparison of calculation results of 1D and 2D models

圖4 二維混油分布云圖Fig.4 2D contamination distribution cloud chart

為了探究適合于二維混油純擴散方程的高效求解方法，在保持網格尺寸不變的情況下，通過改變管道長度來對比高斯賽德爾迭代法、多重網格方法和預條件共軛梯度法在不同管段長度下的計算效率。計算效率的判定標準為求解時間，從0 時刻后行油品進入管道開始計時，當終點處前行油品的管道截面處油品平均濃度小于99%時結束計時（表4、圖5）。

表4 不同求解方法計算結果Tab.4 Calculation results of different solution methods

圖5 不同求解方法計算耗時曲線Fig.5 Calculation time curve of different sloution methods

從表1 和圖5 中可以看出，當管道長度相同時，與傳統的高斯賽德爾迭代法相比，預條件共軛梯度法的求解速度更快，在求解過程中只需迭代一次便可得到收斂解。而多重網格方法反而變慢，其主要原因是擴散系數在不同的流態區域大小不同，且差異較大。在多重網格實施過程中，粗網格上的物性參數可能是從不同區域限定得到的，反而導致誤差被放大，最終導致密網格上的收斂速度減緩。

從圖5 中無法直觀看出不同求解方法的加速效果隨管道長度的變化，因此，以傳統的高斯賽德爾迭代法為基準來計算多重網格方法和預條件共軛梯度法與其的加速比。如圖6 所示，隨著管道長度增加，預條件共軛梯度法相對高斯賽德爾迭代法的加速比從6.14 增加到7.41，加速效果越來越好。除此之外，多重網格方法的計算效率也在不斷接近高斯賽德爾迭代法。

圖6 不同求解方法計算耗時加速比Fig.6 Calculation time speed-up ratio of different solution methods

綜上可以看出，預條件共軛梯度法求解混油純擴散方程的效率最高，并且加速效果隨著管道長度的增加而略有提高。

2.4 混油段跟蹤方法

雖然預條件共軛梯度法的求解效率相對于傳統的高斯賽德爾方法已有較大幅度提高，但仍無法滿足長距離管道混油數值模擬研究。從圖5 可以看出，當管道長度增加后，其計算時間顯著增加，采用擬合方法對共軛梯度法計算時間與計算管道長度的關系進行擬合，得

式中：l為管道長度，m；t為計算時間，s。

按照該擬合式，預估300 km 的管道至少需要超過10 d 的模擬時間，而現場混油段只需不到3.5 d的時間便可到達終點，連現場都無法做到混油提前預測。因此，需要對混油模擬過程進行進一步加速。

混油只存在于不同油品的交界面上，而管道的絕大部分區域都是純油品，在進行混油模擬過程中對純油品區域的計算無任何意義。因此，基于預條件共軛梯度法，通過采用混油段跟蹤方法減小計算區域大小，從而實現混油模擬過程的進一步加速，其基本操作過程如圖7 所示。本文選取的混油段起始位置分別為前行油品濃度為10-10和1-10-10的位置。值得注意的是，在保存混油段起始位置的時候，需增加少量網格（本文選取10 個網格），以保證混油界面在求解過程中不會超出求解區域。

圖7 混油段跟蹤思路Fig.7 Tracking idea of contamination segment

在進行不同管道長度下計算效率的對比之前，需對混油段追蹤方法計算結果的正確性進行驗證，以確保該方法不會給計算結果增加誤差。以355 mm×6.4 mm、長度為1 km 的管道，流速為1 m/s 工況為例計算，結果如圖8 所示。從圖8 可以看出，兩條曲線完全重合，并且在浮點數類型自動保存的6 位小數的范圍內誤差為0。因此，混油段跟蹤方法不會給模擬結果帶來額外的誤差。

圖8 混油段跟蹤方法計算結果精度驗證Fig.8 Accuracy verification of calculation results of contamination segment tracking method

由于常規方法求解效率較低，無法應用于較長管道，因此，采用表4 中對應的不同管道長度進行具體的混油段追蹤方法應用分析。利用混油段跟蹤方法計算前幾組算例，所需的時間較短，比較意義不大，因此將對后幾組數據進行重點比較。采用混油段跟蹤方法計算的結果如表5，圖9 所示。

表5 混油段跟蹤方法計算結果Tab.5 Calculation results of contamination segment tracking method

圖9 混油段跟蹤加速結果Fig.8 Acceleration results of contamination segment tracking method

從圖9 中可以看出，混油段跟蹤方法的運用實現了進一步加速，且隨著管道長度的增加，加速比從15.67 增加到48.77，10 km 的管道可實現48.77 倍的加速，這是因為混油段跟蹤方法極大縮小了計算區域，管道越長，混油段占管道全長的比例越小，使得加速效果隨計算區域的變大而大幅度提高。

綜上可以看出，混油段跟蹤方法可顯著提高混油模擬效率，并且只跟蹤混油段，計算區域隨管道全長變化很小，可適用于長距離管道的混油模擬計算。

2.5 現場實例運用

根據前文分析結果，采用呼包鄂成品油順序輸送管道實際混油驗證程序的正確性和適用性。由于無法獲取管道終點處密度變化的數據，因此直接采用管道處終點的混油量進行對比驗證。油品物性參數、管道參數如表3、表6 所示，油品流速采用平均流速。

表6 呼包鄂管道參數Tab.6 Parameters of Hu-Bao-E Pipeline

根據已知現場參數進行數值模擬，計算結果如圖10和表7所示。鄂爾多斯進站的混油量為90.32 m3，與現場數據95.00 m3之間的相對誤差為4.93%。由于模擬過程中未考慮局部管件、重力、地形、流速變化等因素對混油量的影響，模擬結果應比實際過程偏小一些，因此本文模擬計算結果在合理范圍內。另一方面，現場混油界面從起點到終點所需時間約為275 005.0 s，而模擬耗時為3 293.5 s，為現場所需時間的1.12%，可實現現場混油量的快速超前預測。

在相同的管長情況下，如果不采用混油段追蹤方法，只采用預條件共軛梯度法求解一維模型，所需的計算時間為3 592.7 s，比二維模型采用最終方法計算所需時間還長，所以預條件共軛梯度法+混油段跟蹤方法可實現快速、精確的混油數值模擬。

圖10 呼包鄂管道混油模擬結果Fig.10 Simulation results of contaminated oil in Hu-Bao-E pipeline

表7 模型驗證結果Tab.7 Validation results of model

3 結論

從優選求解方法和減小計算區域兩方面研究提高混油數值模擬的計算效率。對比得到預條件共軛梯度法可實現混油代數方程組的高效求解；通過實現混油段的跟蹤模擬，進一步大幅度提高了混油數值模擬的計算效率。將本文研究成果成功運用于現場管道，得到具體結論如下：

（1）預條件共軛梯度法可實現混油純擴散過程的高效計算。在不同管道長度下，預條件共軛梯度法的求解效率均高于傳統的高斯賽德爾迭代法，針對0.3～5 km 的管道，實現了6.14～7.41 倍的加速，隨著管道長度的增長（網格量變大），預條件共軛梯度法的優勢更加明顯。但多重網格方法由于擴散系數隨空間變化，導致求解效率比傳統方法更低。

（2）混油段跟蹤方法在保證精度的前提下可實現混油數值模擬效率的顯著提升。在求解過程中，純油品區域不進行求解而只求解混油段，因此管道越長，混油段所占的比例越小，加速比越大，10 km的管道實現了48.77 倍的加速。更為重要的是，該方法不受管道長度的限制，計算量只依賴于混油段長度，可適用于長距離管道。

（3）研究成果可適用于現場對混油的提前預測。通過將研究成果運用于呼包鄂成品油順序輸送管道，所得混油結果的相對誤差為4.93%，模擬耗時為現場混油界面從起點到終點所需時間的1.12%。因此，本文研究成果可為現場管道混油分布及混油量的高效模擬提供參考。