基于動(dòng)態(tài)摩阻模型的大落差輸油管道不穩(wěn)定流動(dòng)模擬

李旺，馮亮, 黃立宇, 廖代兵, 陳莎, 賈文龍, 孫璐璐

(1.國(guó)家管網(wǎng)集團(tuán)西南管道有限責(zé)任公司，成都 610000；2.西南石油大學(xué)石油與天然氣工程學(xué)院，成都 610000)

中國(guó)西南地區(qū)的中緬原油管道、蘭成渝成品油管道、蘭成原油管道都是大落差輸油管道，管道高低點(diǎn)的局部落差高達(dá)500～1 500 m，壓力差值高達(dá)4～8 MPa[1]。特別是在管道停泵、關(guān)閥等不穩(wěn)定流動(dòng)工況下，管道局部高點(diǎn)處容易出現(xiàn)低壓，使油品中的輕質(zhì)組分揮發(fā)，從而在高點(diǎn)處出現(xiàn)氣泡聚集，形成液柱分離。過(guò)低的壓力容易導(dǎo)致管道壓癟，而氣泡破裂后容易形成彌合水擊導(dǎo)致管道內(nèi)出現(xiàn)超壓和爆管。準(zhǔn)確分析大落差輸油管道在不穩(wěn)定工況下的壓力、流量變化是合理預(yù)判彌合水擊現(xiàn)象的基礎(chǔ)和前提。

在管道不穩(wěn)定流動(dòng)模擬[2]中，水力摩阻系數(shù)是計(jì)算管道壓力及流量損耗的重要參數(shù)。傳統(tǒng)的摩阻處理方法多采用擬穩(wěn)態(tài)摩阻假設(shè)，即管壁切應(yīng)力采用Darcy-Weisbach摩阻公式計(jì)算[3]。然而，在管道的不穩(wěn)定流動(dòng)過(guò)程中，油品在管道內(nèi)始終處于不斷加速、減速狀態(tài)，流體的慣性會(huì)和管壁之間產(chǎn)生相比于穩(wěn)態(tài)流動(dòng)更大的摩阻損失(即動(dòng)態(tài)摩阻)[4]。實(shí)驗(yàn)和理論研究表明，忽略不穩(wěn)定流動(dòng)中流體慣性影響(即動(dòng)態(tài)摩阻)，會(huì)低估水擊摩阻值，不能準(zhǔn)確預(yù)測(cè)水擊壓力波的衰減過(guò)程和波形畸變[5-7]。

動(dòng)態(tài)摩阻系數(shù)與多因素有關(guān)，計(jì)算較為復(fù)雜，許多學(xué)者對(duì)此進(jìn)行了實(shí)驗(yàn)和理論研究[8]，并提出了多種模型。Zielke[9]基于Navier-Stokes方程，提出了層流瞬變流的管壁切應(yīng)力計(jì)算模型，該模型計(jì)算簡(jiǎn)單，不依賴經(jīng)驗(yàn)參數(shù)，對(duì)計(jì)算與頻率相關(guān)的摩阻損失精度較高，但需要考慮每一時(shí)步的流速值，計(jì)算時(shí)間長(zhǎng)、內(nèi)存消耗大。Trikha[10]采用加權(quán)函數(shù)改進(jìn)了Zielke模型，使得計(jì)算只需前一步的值，簡(jiǎn)化了計(jì)算過(guò)程。Martin等[11]將管道分為近壁層流區(qū)和核心區(qū)，推導(dǎo)出了適用于水力光滑管中任意雷諾數(shù)湍流非恒定摩阻計(jì)算模型。Brunone等[12]、Vardy等[13-14]考慮瞬時(shí)當(dāng)?shù)丶铀俣群蛯?duì)流加速度，得到了非恒定摩阻計(jì)算模型，更真實(shí)地模擬了不穩(wěn)定流摩阻。近年來(lái)，有學(xué)者注意到了動(dòng)態(tài)摩阻對(duì)水利工程的影響[15-17]，研究了動(dòng)態(tài)摩阻對(duì)輸水管道泄漏過(guò)程壓力和流量畸變的影響[18-20]，但國(guó)內(nèi)外對(duì)于動(dòng)態(tài)摩阻對(duì)輸油管道的不穩(wěn)定流動(dòng)過(guò)程的影響研究相對(duì)較少。

基于以上研究，在管道不穩(wěn)定流動(dòng)連續(xù)性方程、動(dòng)量方程和能量方程基礎(chǔ)上，現(xiàn)引入Brunone-Vitkovsky動(dòng)態(tài)摩阻模型[21]，結(jié)合特征線方法[22]和Newton-Raphson方法，對(duì)大落差輸油管道不穩(wěn)定流動(dòng)工況進(jìn)行模擬，分析大落差管道水力瞬變特征。

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

原油在長(zhǎng)輸管道內(nèi)的流動(dòng)可視為一維流動(dòng)。根據(jù)流體流動(dòng)的質(zhì)量守恒、動(dòng)量守恒、和能量守恒原理，建立管道內(nèi)的一維不穩(wěn)定流動(dòng)模型為

(1)

式(1)中：U為自變量列向量，U=[P,V,T]T；A為系數(shù)矩陣；B為列向量。穩(wěn)態(tài)流動(dòng)條件下，原油單位面積管內(nèi)壁的摩擦力τ可按Darcy-Weisbach[13]摩阻公式計(jì)算。對(duì)于不穩(wěn)定流動(dòng)，τ按Brunone-Vitkovsky動(dòng)態(tài)摩阻公式計(jì)算[9-11]，公式為

(2)

A=

(3)

B=

(4)

式中：kt、kx分別為當(dāng)?shù)丶铀俣群瓦w移加速度系數(shù)；sign(V)為符號(hào)函數(shù)，當(dāng)V≥0時(shí)，sign(V)=1；當(dāng)V<0時(shí)，sign(V)=-1；ρ為介質(zhì)密度，kg/m3；P為管內(nèi)壓力，Pa；θ為管道與水平面間的傾角，rad；t為時(shí)間變量，s；V為介質(zhì)流速，m/s；d為管道內(nèi)徑，m；g為重力加速度，m/s2；a為流體波速，m/s；τ為流體與單位面積管內(nèi)壁的摩擦力，Pa；x為管道位置變量，m；K為管道的總熱傳遞系數(shù)，W/(m2·K)；Cv為介質(zhì)比熱容，J/(kg·K)；T為介質(zhì)溫度，K；T0為地溫，K；f為水力摩阻系數(shù)；D為管道外徑，m。

2 模型求解

2.1 特征線方程

基本方程[式(1)]為擬線性雙曲型偏微分方程組，它有3條特征線，可以沿特征線將其轉(zhuǎn)化為全微分方程，然后對(duì)全微分方程積分便可得易于數(shù)值處理的全微分方程。從物理意義上講，特征線即為壓力波和溫度波的傳播路徑。為此，建立了基于特征線法的模型求解格式。根據(jù)線性代數(shù)理論，系數(shù)矩陣A的3個(gè)特征值表達(dá)式為

(5)

式(5)中：

(6)

(7)

如考慮kx=kt=k3，且a>>V，則有

(8)

(9)

k3按式(10)取值為

(10)

式(10)中：對(duì)于層流，C=0.004 76；對(duì)于紊流，有

(11)

由Li(A-λiE)=0可得每個(gè)λi的相應(yīng)特征向量Li為

(12)

由此可得3個(gè)特征值對(duì)應(yīng)特征線方程及其相容性方程如式(13)～式(15)所示。

前向特征方程C+：

(13)

后向特征方程C-：

(14)

將連續(xù)性方程帶入溫度特征方程，可得溫度特征方程C為

(15)

2.2 有限差分方程

管道瞬變流動(dòng)計(jì)算要求給出管道沿線的壓力、流量和溫度隨時(shí)間變化的規(guī)律，為此將管道沿時(shí)間和長(zhǎng)度方向離散成矩形網(wǎng)格，如圖1所示。為保證解的穩(wěn)定性，網(wǎng)格長(zhǎng)度必須滿足Courant條件[13]，即

圖1 管道不穩(wěn)定流動(dòng)差分網(wǎng)格示意圖

(16)

對(duì)于上述特征方程中的相容性方程沿各自的特征線進(jìn)行積分，對(duì)摩阻項(xiàng)采用二階近似來(lái)防止解的不穩(wěn)定?？梢缘玫紺+、C-、C相對(duì)應(yīng)的差分方程為

(17)

(18)

(19)

式中：Q為管道內(nèi)輸量，m3/s；A為管道截面積，m2；下角標(biāo)參數(shù)N、S、R、H見(jiàn)圖1；其他參數(shù)見(jiàn)第1節(jié)。

的特征線與上一時(shí)間層網(wǎng)格線的交點(diǎn)。S、N和R點(diǎn)的值要通過(guò)對(duì)上一時(shí)步的值進(jìn)行線性插值而得到。上述方程僅能用于計(jì)算管道內(nèi)節(jié)點(diǎn)的參數(shù)，并且以穩(wěn)態(tài)時(shí)的水力參數(shù)作為初始值。對(duì)于上游邊界點(diǎn)，需要將C-上的特征線方程與邊界條件聯(lián)解得到。對(duì)于下游邊界，需要將C+和C上的特征線方程與下游邊界條件聯(lián)解得到。對(duì)于中間泵站、閥門等非管元件，需要將C+、C-、C上的特征線方程與泵的特性方程，閥門方程聯(lián)立求解得到。上述由邊界條件、以及有限差分方程構(gòu)成的非線性方程組可采用Newton-Raphson迭代法進(jìn)行求解。

3 算例分析

3.1 算例描述

某大落差加熱輸油管道全長(zhǎng)438 km，管道沿線高程及水熱力分布圖如圖2所示。管道規(guī)格為Φ273.1 mm×8.7 mm，不設(shè)保溫層，管道總傳熱系數(shù)取2.6 W/(m2·K)，環(huán)境溫度為30 ℃，全線共設(shè)置4個(gè)熱泵站，各熱泵站高程、里程位置如表1所示。管道最高點(diǎn)為P(100.94，3.424)km，沿線最大落差為549 m。穩(wěn)定流動(dòng)時(shí)管道內(nèi)流量為418.6 m3/h；為減少油品密度、黏度等參數(shù)對(duì)于水力計(jì)算結(jié)果的影響，取輸送油品相對(duì)密度為0.82，動(dòng)力黏度為120 mPa·s。

表1 站場(chǎng)布置及進(jìn)出站壓頭、溫度參數(shù)

圖2 管道沿線高程及水熱力分布圖

3.2 計(jì)算過(guò)程

基于考慮動(dòng)態(tài)摩阻的輸油管道不穩(wěn)定流動(dòng)模型，采用C#語(yǔ)言編制了輸油管道不穩(wěn)定流動(dòng)計(jì)算程序，分析大落差管道的不穩(wěn)定流動(dòng)特征。設(shè)定管道起點(diǎn)、終點(diǎn)邊界條件分別為壓力溫度和流量邊界條件，結(jié)合泵站的進(jìn)出站壓頭，計(jì)算得到沿線壓力、溫度和流量分布。以此穩(wěn)態(tài)條件作為初始值，分析No.3泵站停運(yùn)的瞬態(tài)工況。停泵運(yùn)行工況具體設(shè)置如下：分析No.3泵站在1 s時(shí)兩臺(tái)泵同時(shí)停運(yùn)，同時(shí)當(dāng)No.3站上游壓力大于下游時(shí)旁通閥門打開(kāi)工況下的權(quán)限壓力和流量變化。分別按照擬穩(wěn)態(tài)摩阻和動(dòng)態(tài)摩阻方法計(jì)算管道內(nèi)的壓力和流量變化，對(duì)比高點(diǎn)P與No.4站的進(jìn)站壓力及流量的差異，結(jié)果如圖3～圖6所示。

圖3 高點(diǎn)P處的壓力變化

圖4 高點(diǎn)P處的流量變化

圖5 No.4站進(jìn)站壓力變化

圖6 No.4站進(jìn)站流量變化

3.3 結(jié)果分析

圖3和圖4的計(jì)算結(jié)果表明，在No.3站兩臺(tái)泵同時(shí)停運(yùn)后，其增壓波經(jīng)過(guò)87 s傳遞到最高點(diǎn)P(圖3)，在此之前由兩種方法計(jì)算P點(diǎn)的壓力、流量一致，表明在穩(wěn)定流動(dòng)條件下動(dòng)態(tài)摩阻對(duì)管道水力參數(shù)的影響可忽略不計(jì)；87 s后，由兩種方法計(jì)算的P點(diǎn)壓力、流量開(kāi)始出現(xiàn)差異，最大差異出現(xiàn)在208 s時(shí)刻，此時(shí)P點(diǎn)波動(dòng)最劇烈(流量、壓力變化最快)，壓力差0.024 4 MPa(圖3)，流量差6.01 m3/h(圖4)。基于穩(wěn)態(tài)摩阻和動(dòng)態(tài)摩阻計(jì)算的壓力、流量之間的最大相對(duì)偏差分別為2.06%、2.11%，且基于穩(wěn)態(tài)摩阻計(jì)算的壓力、流量值偏高。

圖5和圖6的計(jì)算結(jié)果表明，No.3站在1 s后兩臺(tái)泵同時(shí)停運(yùn)，其減壓波經(jīng)過(guò)162 s傳遞到No.4站，之前由兩種方法計(jì)算No.4站進(jìn)站壓力、流量也是一致的。但是，162 s后，由兩種方法計(jì)算的No.4站壓力、流量開(kāi)始出現(xiàn)差異，最大差異出現(xiàn)在300 s時(shí)刻，此時(shí)No.4站處波動(dòng)也是最劇烈的(流量、壓力變化最快)，壓力差0.027 6 MPa(圖5)，流量差1.49 m3/h(圖6)?；诜€(wěn)態(tài)摩阻和動(dòng)態(tài)摩阻計(jì)算的壓力、流量之間的最大相對(duì)偏差分別為3.24%、0.43%，且基于穩(wěn)態(tài)摩阻計(jì)算的壓力、流量值偏高。

模擬結(jié)果表明，當(dāng)管道運(yùn)行過(guò)程中發(fā)生不穩(wěn)定流動(dòng)時(shí)，由兩種方法計(jì)算的壓力和流量會(huì)出現(xiàn)較大的差異。這是因?yàn)椴环€(wěn)定流動(dòng)過(guò)程中的動(dòng)態(tài)摩阻是由擬穩(wěn)態(tài)摩阻項(xiàng)和附加摩阻項(xiàng)之和進(jìn)行表征的，動(dòng)態(tài)摩阻計(jì)算方法考慮了附加摩阻項(xiàng)對(duì)管道不穩(wěn)定流動(dòng)過(guò)程的影響，使計(jì)算結(jié)果更為準(zhǔn)確的反應(yīng)管道不穩(wěn)定流動(dòng)過(guò)程，這與式(4)是吻合的。由于基于動(dòng)態(tài)摩阻計(jì)算的管道壓力低于穩(wěn)態(tài)摩阻值，因此能夠更加準(zhǔn)確地預(yù)測(cè)在大落差管道高點(diǎn)可能出現(xiàn)的負(fù)壓，進(jìn)而分析其可能引發(fā)的彌合水擊等現(xiàn)象，為采取更加切實(shí)可行的管道安全保障措施提供理論和技術(shù)支撐。

4 結(jié)論

(1)考慮輸油管道不穩(wěn)定流動(dòng)過(guò)程中由流體不穩(wěn)定流動(dòng)帶來(lái)的流體慣性影響，基于一維流動(dòng)的質(zhì)量、動(dòng)量、能量守恒方程，以及由穩(wěn)定流摩阻和附加項(xiàng)摩阻構(gòu)成的Brunone-Vitkovsky動(dòng)態(tài)摩阻模型，建立了大落差管道不穩(wěn)定流動(dòng)模型。

(2)采用特征線法、有限差分法，求解輸油管道不穩(wěn)定流動(dòng)數(shù)學(xué)模型，分析了某大落差輸油管道中間泵站停運(yùn)后站前局部高點(diǎn)與下游泵站進(jìn)站的壓力、流量變化趨勢(shì)，對(duì)比穩(wěn)態(tài)摩阻模型與瞬態(tài)模型的結(jié)果差異。

(3)瞬態(tài)摩阻模型計(jì)算局部高點(diǎn)處的壓力、流量較穩(wěn)態(tài)摩阻模型的計(jì)算值之間的最大偏差分別為2.06%、2.11%，下游泵站進(jìn)站壓力、流量的最大偏差分別為3.24%、0.43%，且基于穩(wěn)態(tài)摩阻計(jì)算的壓力、流量值偏高。基于動(dòng)態(tài)摩阻模型對(duì)大落差管道不穩(wěn)定流動(dòng)工況進(jìn)行模擬，對(duì)于提高管道壓力、流量等工藝參數(shù)的預(yù)測(cè)精度，制訂更加有效的不穩(wěn)定流動(dòng)安全防護(hù)措施具有重要意義。