王軍防 陳雪嬌 沈允 陳宇杰

1中國石化管道儲(chǔ)運(yùn)有限公司科技研發(fā)中心

2北京航天試驗(yàn)技術(shù)研究所

3中國石油北京油氣調(diào)控中心

4中國石油大學(xué)（北京）

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

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

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

1 模型的建立

1.1 物理模型

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

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

1.2 數(shù)學(xué)模型

1.2.1 一維模型

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

式中：c為前行油品截面處油品平均濃度（體積分?jǐn)?shù)，下同）；um為油品的平均流速，m/s；K為有效擴(kuò)散系數(shù)，可采用式（2）計(jì)算，m2/s；t為時(shí)間，s；x為軸向坐標(biāo)位置，m。

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

1.2.2 二維模型

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

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

不同流態(tài)區(qū)域的速度分布和擴(kuò)散系數(shù)采用經(jīng)驗(yàn)公式計(jì)算[12]。

1.2.3 邊界條件和初始條件

管道中充滿前行油品，從0 時(shí)刻開始，后行油品進(jìn)入管道，二維混油數(shù)值模擬在管道壁面和中心處沿徑向的濃度梯度為0，邊界條件見式（6），一維混油數(shù)值模擬邊界條件見式（7）。值得注意的是，對于二維模型，壁面處采用無滑移邊界條件，無法有效考慮壁面對油品的黏附作用。因此，無法考慮油品輸送順序?qū)煊蛶淼挠绊憽?/p>

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

1.3 網(wǎng)格劃分

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

式中：yj為y方向第j個(gè)節(jié)點(diǎn)的位置；N為y方向總的節(jié)點(diǎn)個(gè)數(shù)；a為網(wǎng)格變化的調(diào)節(jié)參數(shù)；?j為中間變量。

1.4 方程求解

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

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

2 模擬結(jié)果

為高效準(zhǔn)確地進(jìn)行混油數(shù)值研究，首先對比一維和二維模型的求解結(jié)果和求解效率，其次從高效求解方法和減少計(jì)算區(qū)域兩個(gè)角度分別提升二維混油模型計(jì)算效率，最后將高效求解方法運(yùn)用到二維模型中，進(jìn)行現(xiàn)場實(shí)例驗(yàn)證分析，并和一維模型計(jì)算效率相對比。

2.1 網(wǎng)格無關(guān)解驗(yàn)證

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

表1 二維模型網(wǎng)格無關(guān)解驗(yàn)證的網(wǎng)格信息Tab.1 Grid information for grid independent solution verification of 2D model

表2 一維模型網(wǎng)格無關(guān)解驗(yàn)證的網(wǎng)格信息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 中的三種網(wǎng)格條件下，二維模型和一維模型在終點(diǎn)處混油分布如圖2 所示。兩種模型在三套網(wǎng)格下的計(jì)算結(jié)果近似，其中二維模型1 000×20 的網(wǎng)格和2 000×40 的網(wǎng)格所計(jì)算的混油長度為57.5 m 和57.75 m，相對差異為0.43%；一維模型1 000 個(gè)網(wǎng)格和2 000 個(gè)網(wǎng)格的計(jì)算的混油長度為58.06 m 和59.03 m，結(jié)果相差1.64%。所以可選擇1 000×20 的網(wǎng)格作為二維模型的計(jì)算網(wǎng)格，1 000 的網(wǎng)格作為一維模型的計(jì)算網(wǎng)格，基于相應(yīng)的網(wǎng)格比例開展具體混油數(shù)值模擬研究。

圖2 一維和二維網(wǎng)格無關(guān)解驗(yàn)證Fig.2 Verification of 1D and 2D models gird independent solutions

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

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

2.3 混油高效求解方法探究

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

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

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

表4 不同求解方法計(jì)算結(jié)果Tab.4 Calculation results of different solution methods

圖5 不同求解方法計(jì)算耗時(shí)曲線Fig.5 Calculation time curve of different sloution methods

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

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

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

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

2.4 混油段跟蹤方法

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

式中：l為管道長度，m；t為計(jì)算時(shí)間，s。

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

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

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

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

圖8 混油段跟蹤方法計(jì)算結(jié)果精度驗(yàn)證Fig.8 Accuracy verification of calculation results of contamination segment tracking method

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

表5 混油段跟蹤方法計(jì)算結(jié)果Tab.5 Calculation results of contamination segment tracking method

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

從圖9 中可以看出，混油段跟蹤方法的運(yùn)用實(shí)現(xiàn)了進(jìn)一步加速，且隨著管道長度的增加，加速比從15.67 增加到48.77，10 km 的管道可實(shí)現(xiàn)48.77 倍的加速，這是因?yàn)榛煊投胃櫡椒O大縮小了計(jì)算區(qū)域，管道越長，混油段占管道全長的比例越小，使得加速效果隨計(jì)算區(qū)域的變大而大幅度提高。

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

2.5 現(xiàn)場實(shí)例運(yùn)用

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

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

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

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

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

表7 模型驗(yàn)證結(jié)果Tab.7 Validation results of model

3 結(jié)論

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

（1）預(yù)條件共軛梯度法可實(shí)現(xiàn)混油純擴(kuò)散過程的高效計(jì)算。在不同管道長度下，預(yù)條件共軛梯度法的求解效率均高于傳統(tǒng)的高斯賽德爾迭代法，針對0.3～5 km 的管道，實(shí)現(xiàn)了6.14～7.41 倍的加速，隨著管道長度的增長（網(wǎng)格量變大），預(yù)條件共軛梯度法的優(yōu)勢更加明顯。但多重網(wǎng)格方法由于擴(kuò)散系數(shù)隨空間變化，導(dǎo)致求解效率比傳統(tǒng)方法更低。

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

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