何國璽，梁永圖*，李巖松，劉勝利，吳夢雨，謝成，席罡，李豐

1 中國石油大學(xué)(北京)城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室, 北京 102249

2 中國石化銷售華南分公司, 廣州 510620

3 中石化北海液化天然氣有限責(zé)任公司, 北海 536000

成品油管道打孔盜油量測算方法

何國璽1，梁永圖1*，李巖松1，劉勝利1，吳夢雨1，謝成2，席罡2，李豐3

1 中國石油大學(xué)(北京)城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室, 北京 102249

2 中國石化銷售華南分公司, 廣州 510620

3 中石化北海液化天然氣有限責(zé)任公司, 北海 536000

打孔盜油的事故危害很大，不僅會(huì)造成油品損失、影響管道正常運(yùn)行，還極易發(fā)生爆炸、泄漏、污染環(huán)境等二次事故。目前的研究著重于從管理和技術(shù)角度，對盜油事故進(jìn)行甄別，并進(jìn)行定位，但對盜油量的計(jì)算研究較少。分析了打孔盜油過程中整個(gè)管道系統(tǒng)的水力瞬變，以打孔盜油點(diǎn)、變徑點(diǎn)、批次界面和閥門等為邊界條件，建立了打孔盜油模型。根據(jù)盜油過程上下游壓力、上游流量和溫度等參數(shù)的變化，計(jì)算了整個(gè)盜油過程的瞬時(shí)盜油流量和累計(jì)盜油量。現(xiàn)場實(shí)驗(yàn)及盜油量測算結(jié)果表明：(1)孔口自由出流和短管出流的出流系數(shù)差別較大，因此計(jì)算盜油量時(shí)，應(yīng)考慮閥門的局部摩阻、盜油管線摩阻、盜油罐車內(nèi)的油品壓力等參數(shù)；(2)對盜油速率和累計(jì)盜油量影響從大到小的參數(shù)依次是：盜油點(diǎn)孔徑(閥門開度)、盜油管段上下游壓力變化、流量、溫度；(3)輸量越小、壓力越高、壓力波動(dòng)程度越大，模型計(jì)算結(jié)果精度越高。

成品油管道；打孔盜油；水熱力瞬變；盜油量測算

0 引言

打孔盜油時(shí)整條輸油管道和打孔盜油管、減壓閥、儲(chǔ)油罐構(gòu)成一個(gè)完整的水力系統(tǒng)，盜油事故發(fā)生后，管線和盜油設(shè)備整體處于水力瞬變狀態(tài)，盜油處產(chǎn)生減壓波同時(shí)向上下游傳播，管內(nèi)流量和壓力均會(huì)發(fā)生變化。由于盜油水力系統(tǒng)結(jié)構(gòu)復(fù)雜，水力瞬變過程復(fù)雜，目前完整準(zhǔn)確模擬整個(gè)過程的研究較少，也無計(jì)算此過程的盜油量的方法。雖有文獻(xiàn)指出泄漏的結(jié)果是上游流量增大、下游流量減小、上下游壓力都降低[1-3]，但打孔盜油過程與小孔泄漏過程具有一定的差異[4]。因打孔盜油過程存在人為控制因素，故其過程中的水力參數(shù)變化規(guī)律與小孔泄漏過程不完全相同。目前計(jì)算打孔盜油量主要有2種方法。一是根據(jù)上下游的流量數(shù)據(jù)，將流量測量值之差與盜油持續(xù)時(shí)間相乘即得到盜油量，該方法的誤差較大。二是將打孔盜油過程看作小孔泄漏，先測算出泄漏處的管內(nèi)壓力[5-9]，再根據(jù)小孔內(nèi)外壓差計(jì)算泄漏量[10-13]。這種方法的問題在于將盜油管線出流等效為小孔出流，計(jì)算誤差較大，而且無法得到準(zhǔn)確的小孔泄漏系數(shù)。另一方面由于管道內(nèi)外壓差對盜油量計(jì)算的影響較大，將打孔盜油當(dāng)作小孔泄漏時(shí)無法得到準(zhǔn)確的小孔外部壓力分布。還有學(xué)者利用數(shù)值模擬法[14-16]、統(tǒng)計(jì)學(xué)法[17-18]、事后間接法[19]、實(shí)驗(yàn)分析法[20-22]等計(jì)算管道在小孔泄漏過程的泄漏量。實(shí)際上，打孔盜油管路管內(nèi)(閥后)壓力明顯高于大氣壓，盜油出流過程中既有閥門節(jié)流摩阻，也有盜油管路摩阻，將油盜至罐車時(shí)，還有密閉罐車內(nèi)的油品蒸汽壓產(chǎn)生的背壓。因此用小孔泄漏的計(jì)算方法來計(jì)算打孔盜油的盜油量往往造成計(jì)算結(jié)果偏大。本文著重研究盜油點(diǎn)的水力特性，同時(shí)還考慮了管路、閥門、罐車等打孔盜油相關(guān)設(shè)備對摩阻、出流特性和背壓造成的影響。

1 打孔盜油模型

1.1 問題分析

成品油管道的穩(wěn)定運(yùn)行是指管內(nèi)油品的水熱力參數(shù)和外界環(huán)境參數(shù)均達(dá)到相對恒定時(shí)的狀態(tài)。打孔盜油發(fā)生后，管道運(yùn)行狀態(tài)變?yōu)榉欠€(wěn)定，產(chǎn)生的負(fù)壓波向管道兩端傳遞，因管道沿線摩阻的存在，產(chǎn)生的壓力波很快衰減，最后完全消失，達(dá)到新的穩(wěn)定狀態(tài)。盜油發(fā)生初期，壓力下降快，隨著時(shí)間增加，壓力逐漸趨于平穩(wěn)。打孔盜油后現(xiàn)場一般將采取關(guān)閉打孔盜油管段兩端閥門以及停運(yùn)相關(guān)機(jī)組等相應(yīng)措施，此時(shí)管道內(nèi)將發(fā)生復(fù)雜的瞬變流動(dòng)現(xiàn)象，打孔盜油點(diǎn)的盜油流量也隨之有規(guī)律地變化。

圖1 某真實(shí)打孔盜油發(fā)生過程示意圖Fig. 1 Illustration of a real event of oil stolen by drilling hole

結(jié)合上述分析可知，影響打孔盜油過程的因素有很多，如圖1所示，包括：環(huán)境參數(shù)，如地形、環(huán)境溫度；管道參數(shù)，如管徑、管長、粗糙度、導(dǎo)熱系數(shù)等；油品性質(zhì)，如密度、黏度、比熱容、導(dǎo)熱系數(shù)、體積壓縮系數(shù)、溫度膨脹系數(shù)等；流動(dòng)狀態(tài)；其他邊界條件，如打孔盜油孔特性、混油界面位置、變徑點(diǎn)、閥門、泵、液柱分離等。此外，非穩(wěn)定狀態(tài)下的打孔盜油還與管道運(yùn)行參數(shù)有關(guān)，是一個(gè)動(dòng)態(tài)過程。因此，不同時(shí)間、不同位置處油品的流量、壓力存在很大差異。當(dāng)管內(nèi)油品流量發(fā)生變化時(shí)，油品的摩擦熱會(huì)改變，進(jìn)而引起油品溫度改變，尤其在流量較大時(shí)更是如此。油品溫度改變使得其黏度及水力摩阻系數(shù)發(fā)生變化，進(jìn)而使流量和壓力等發(fā)生變化，流量變化對全線溫度又將產(chǎn)生一定影響。另外，成品油管道順序輸送多種油品，前后行油品的物性相差很大，而水熱力動(dòng)態(tài)仿真模型中的參數(shù)受油品物性的影響，故對順序輸送管道進(jìn)行動(dòng)態(tài)仿真之前必須結(jié)合混油界面跟蹤，明確管道中油品的種類、數(shù)量等相關(guān)基礎(chǔ)參數(shù)。

依據(jù)油品在管道內(nèi)的流動(dòng)狀態(tài)和流體相態(tài)將打孔盜油過程劃分為3個(gè)階段：(1)液態(tài)油品瞬變漏失過程(因打孔盜油引起的瞬變)；(2)液態(tài)油品穩(wěn)定漏失過程(打孔盜油后達(dá)到新的穩(wěn)態(tài))；(3)液態(tài)油品瞬變漏失(執(zhí)行相應(yīng)的停輸方案后引起的瞬變)。對上述各階段分別建立與之對應(yīng)的數(shù)學(xué)模型并求解，得到各階段相應(yīng)時(shí)段內(nèi)打孔盜油點(diǎn)處的壓力、溫度、流速、相態(tài)變化，從而實(shí)時(shí)跟蹤記錄油品的打孔盜油過程。

1.2 物理模型

成品油管道的非穩(wěn)態(tài)工況是一個(gè)管道內(nèi)油品的流動(dòng)、傳熱、熱力及水力相耦合過程。計(jì)算打孔盜油量的數(shù)學(xué)模型包含連續(xù)性方程、動(dòng)量守恒方程、能量守恒方程以及相關(guān)的初始條件和邊界條件。為得到打孔盜油過程中油管段兩端的壓力(作為邊界條件)，選擇上、下游泵站間的管道作為計(jì)算的基本單元。為保證模型的通用性，打孔盜油管段包含的物理單元有閥門、變徑點(diǎn)、混油界面、打孔盜油點(diǎn)、上游泵站出口和下游泵站入口，如圖2所示。

對于該耦合問題，做如下假設(shè)：(1)流動(dòng)為沿管道軸向的一維流動(dòng)，管截面流速為均值；(2)截面上的油品溫度均勻分布，僅與軸向位置和運(yùn)行時(shí)間長短有關(guān)；(3)忽略軸向?qū)幔?4)不考慮管道外部土壤溫度的變化，以平均地溫計(jì)算；(5)總傳熱系數(shù)是一個(gè)常數(shù)；(6)不考慮混油，假設(shè)相鄰油品接觸面為管道截面大??；(7)不考慮打孔盜油點(diǎn)外界環(huán)境對打孔盜油的影響；(8)不考慮油品經(jīng)過閥門、變徑點(diǎn)等引起的溫度變化；(9)流體為不可壓縮流體；(10)盜油罐內(nèi)背壓保持不變。

基于以上假設(shè)，描述管內(nèi)油流的連續(xù)性方程、動(dòng)量方程和能量方程如下：

圖2 通用物理模型Fig. 2 General physical model

根據(jù)上述基本方程可以得到管內(nèi)油流的換熱方程：

式中，t為時(shí)間，s；ρ為管截面上油品的平均密度，kg/m3；x為距盜油管段起點(diǎn)的距離，m；v為管內(nèi)油品的平均流速，m/s；g為重力加速度，m/s2；A為盜油管段截面積，m2；θ是盜油管段與水平方向的夾角，rad；p為油品在管截面上的壓力，Pa；D為管內(nèi)徑，m；λ為達(dá)西摩阻系數(shù)；e為油品比內(nèi)能，J/kg；h為油品比焓，J/kg；s為相鄰計(jì)算節(jié)點(diǎn)間的高程差，m；q為油品與單位面積管壁單位時(shí)間內(nèi)的熱流密度，W/m2；αp為油品的膨脹系數(shù)，1/℃；T為管內(nèi)油品溫度，℃；c為油品的比熱容，J/(kg·℃)；ao,g為壓力波在不同介質(zhì)中的傳播速度，m/s；ρo為不同批次的油品密度kg/m3；ko為對應(yīng)批次油品的彈性模量，Pa；Dg為不同管段直徑，m；Eg為對應(yīng)管段的楊氏彈性模量，Pa；δ為對應(yīng)管段的壁厚，m。

埋地成品油管道的傳熱過程由3部分組成，分別為油品至管壁的傳熱、鋼管壁至瀝青層的傳熱及管壁至外部土壤的傳熱。當(dāng)管道運(yùn)行達(dá)到穩(wěn)定狀態(tài)后，在管內(nèi)外建立了穩(wěn)定的溫度場，各部分在同一時(shí)間內(nèi)傳遞的熱量相等。

油流與外界環(huán)境間的熱流密度為：

其中，T為管內(nèi)油品平均溫度，℃；T0為外界環(huán)境溫度，℃；q為熱流密度，W/m2；K為總傳熱系數(shù)，W/(m2·℃)。

根據(jù)特征方程一般形式，建立式(1)、式(2)和式(3)的特征方程：

各管段內(nèi)部節(jié)點(diǎn)的特點(diǎn)是壓力和流量均相同，若采用工程上常用的H(x，t)、Q(x，t)代替p(x，t)、v(x，t)，由于：

1.3 邊界條件

1.3.1 通用邊界條件

管道的首末端邊界為外部邊界，其它都為內(nèi)部邊界。壓力波經(jīng)過邊界時(shí)會(huì)發(fā)生反射，反射后的壓力波與邊界類型有關(guān)；非穩(wěn)態(tài)過程一般起始于邊界的擾動(dòng)，隨著擾動(dòng)的結(jié)束而結(jié)束。

(1)上下游邊界

依據(jù)建立的通用物理模型，因上游只有C-特征線，下游只有熱力特征線v和C+特征線，若要確定上下游邊界處的水熱力參數(shù)，需知上游端溫度、壓力或流量與下游端壓力或流量。上下游邊界示意如圖3所示，其中j為時(shí)間點(diǎn)，i為管段，0、1、N、N-1為節(jié)點(diǎn)，具體描述如下：

圖3 上下游邊界節(jié)點(diǎn)Fig. 3 The upstream and downstream boundary nodes

(2)閥門邊界

閥門邊界又稱為擾動(dòng)邊界，它本身的特性是隨時(shí)間改變的。從水力學(xué)可知，閥門的阻力特性為：

其中，ξ為阻力因數(shù)，取決于閥的結(jié)構(gòu)、開度和口徑；K為集合系數(shù)，等于；w為閥門通道的面積，m2。

當(dāng)閥門型號(hào)確定之后，即確定了w，可將集合系數(shù)K稱為閥的阻力系數(shù)。阻力系數(shù)與閥門開度的關(guān)系稱為閥的阻力特性，是閥本身的特性。閥門邊界的特點(diǎn)是，閥門上下游的壓力不等，但流量相同。閥門邊界示意如圖4所示，其中分別為N、0節(jié)點(diǎn)對應(yīng)的參數(shù)值，具體描述如下：

(3)變徑點(diǎn)邊界

因變徑點(diǎn)兩側(cè)管徑不同，水力特征線不同，壓力波在變徑點(diǎn)處發(fā)生反射。變徑點(diǎn)邊界的特點(diǎn)是，邊界上下游的壓力和流量均相同。變徑點(diǎn)邊界示意如圖5所示，具體描述如下：

(4)混油邊界

把混油段中央放置于相應(yīng)的差分節(jié)點(diǎn)上，以它作為內(nèi)部邊界，認(rèn)為節(jié)點(diǎn)上游管段內(nèi)全是后行油品，下游管段全是前行油品。當(dāng)混油界面移動(dòng)到下一節(jié)點(diǎn)時(shí)，混油界面邊界相應(yīng)的移至該節(jié)點(diǎn)。因節(jié)點(diǎn)之間的長度(即距步)遠(yuǎn)大于每一時(shí)步混油界面移動(dòng)的距離，故混油界面由當(dāng)前節(jié)點(diǎn)移動(dòng)至下一節(jié)點(diǎn)需要經(jīng)歷很長一段時(shí)間，在此不考慮混油界面移動(dòng)至節(jié)點(diǎn)之間時(shí)對計(jì)算帶來的影響。混油邊界的特點(diǎn)是，該點(diǎn)壓力和流量均相同?；煊瓦吔缡疽馊鐖D6所示，具體描述如下：

1.3.2 打孔盜油邊界條件

打孔盜油與小孔泄漏的邊界條件不同之處在于：(1)出流系數(shù)不同，閥門和管線的存在使末端出流不同于小孔出流；(2)局部摩阻不同，閥門和管線的存在會(huì)引入新的摩阻，而小孔泄漏則沒有；(3)環(huán)境壓力不同。埋地管道周圍土壤的壓力與管內(nèi)油品蒸氣壓變化規(guī)律不同。

假定盜油設(shè)備出流時(shí)具有閥門出流的性質(zhì)，則有：

將式(18)和式(19)代入式(17)中有計(jì)算盜油流量公式：

其中，qP為盜油流量，m3/s；C0為流量系數(shù)；為閥門局部摩阻，m；為盜油管摩阻，m，dL為盜油設(shè)備管的直徑，m；LL為盜油設(shè)備管長度，m；zL為盜油點(diǎn)高程，m；pL為盜油罐內(nèi)油品蒸汽壓，Pa。盜油點(diǎn)邊界的特點(diǎn)是：具有公共的節(jié)點(diǎn)壓頭HL，進(jìn)出盜油點(diǎn)流量的代數(shù)和為0。盜油點(diǎn)邊界示意如圖7所示，具體描述如下：

蝶閥的流量系數(shù)按下式計(jì)算，

圖4 閥門邊界節(jié)點(diǎn)Fig. 4 The boundary nodes of valves

圖5 變徑點(diǎn)邊界節(jié)點(diǎn)計(jì)算方法Fig. 5 The boundary nodes of diameter changing points

圖6 混油邊界節(jié)點(diǎn)計(jì)算方法Fig. 6 The boundary nodes of mixed oil interfaces

圖7 盜油邊界節(jié)點(diǎn)計(jì)算方法Fig. 7 The boundary nodes of pipeline damage point

α為出流系數(shù)；Cd為流速系數(shù)；AL為等效圓孔的截面積，m2；DL為等效圓孔的直徑，m。根據(jù)負(fù)壓波傳播規(guī)律式(23)[23-24]和慣性水擊壓力公式(24)有

式(26)是關(guān)于DL的隱函數(shù)，需要迭代求解。其中，Δpu為盜油管段入口的壓力變化，Pa；AG為管道截面積，m2；Δx為劃分管段距步，m；Li為被盜油管線劃分出的第i段，m；Di為第i段的管徑，m；λi為第i段的摩阻系數(shù)；ai為第i段的波速，m/s；vi為第i段的平均流速，m/s。

實(shí)際現(xiàn)場應(yīng)用中，盜油點(diǎn)位置、基本管道數(shù)據(jù)和輸送油品物性均為已知條件，而上下游的壓力、上游的溫度、上游流量可從SCADA系統(tǒng)中獲得，盜油時(shí)間可由壓力變化趨勢預(yù)測，也可由其他方式設(shè)定如發(fā)現(xiàn)盜油后手動(dòng)關(guān)閉盜油閥門。由這些數(shù)據(jù)便可數(shù)值模擬出盜油發(fā)生時(shí)間段內(nèi)管道中每一點(diǎn)處的壓力、流量值。由于SCADA系統(tǒng)記錄數(shù)據(jù)的頻率一般為1 s，因此可以將1 s時(shí)間內(nèi)的油品的流動(dòng)看作準(zhǔn)穩(wěn)態(tài)流動(dòng)。取120 s內(nèi)管道穩(wěn)態(tài)運(yùn)行的數(shù)據(jù)，便可模擬出盜油事故發(fā)生前管道的運(yùn)行狀態(tài)。接著盜油發(fā)生，將管道兩端的頻率為1 s的壓力值和壓降值、上游流量值、上游溫度值作為邊界條件，解方程(1)-(4)便可得到盜油點(diǎn)的壓力值(轉(zhuǎn)換為水頭后為HL)和壓降值ΔpL，將此代入公式(26)中便可迭代計(jì)算DL。由于每1 s的上下游的壓力值和壓降值、上游流量值、上游溫度值都不同，因此計(jì)算得到的盜油點(diǎn)的壓力值和壓降值也不同，用此計(jì)算得到的DL也不同，因此取從壓力開始下降到壓力下降到最低并開始恢復(fù)這段時(shí)間內(nèi)的數(shù)據(jù)計(jì)算得到多個(gè)DL，取平均得到最終的DL。得到DL后，再由式(22)得到盜油孔參數(shù)，最后再由式(20)計(jì)算得到瞬時(shí)盜油流量qP。

2 模型求解

成品油順序輸送非穩(wěn)定流動(dòng)的基本方程為一組擬線性雙曲型偏微分方程，非穩(wěn)定打孔盜油過程涉及的問題為快瞬變流動(dòng)問題。采用特征線法結(jié)合有限差分求解水熱力耦合瞬變流動(dòng)問題，并對算法的穩(wěn)定性和計(jì)算精度進(jìn)行檢驗(yàn)。求解步驟如下(見圖8)：

(1)首先對打孔盜油管段進(jìn)行穩(wěn)態(tài)計(jì)算，確定管道全線的壓力、流量和溫度；

(2)當(dāng)邊界條件發(fā)生變化時(shí)，按照上述邊界條件處理方法，以前一時(shí)步的溫度計(jì)算油品物性，確定管道各節(jié)點(diǎn)當(dāng)前時(shí)步的壓力和流量；

(3)采用v特征方程計(jì)算各節(jié)點(diǎn)油品的溫度，由于水力瞬變的空間步長遠(yuǎn)大于熱力瞬變的空間步長，故在進(jìn)行熱力瞬變計(jì)算時(shí)，需要對一些點(diǎn)的壓力及流量進(jìn)行插值；

(4)由(3)中的計(jì)算結(jié)果，確定油品密度ρ、黏度μ和列賓宗摩阻系數(shù)f；

(5)重復(fù)步驟(2)、(3)和(4)，直到達(dá)到設(shè)定的計(jì)算精度；

(6)對以后每一時(shí)步重復(fù)上述計(jì)算，直至達(dá)到設(shè)定的計(jì)算時(shí)間。

3 算例分析及討論

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

打孔盜油管段全線的管道參數(shù)，包括管外(內(nèi))徑、壁厚、絕對粗糙度、變徑點(diǎn)位置、閥門位置、管道長度、管道的楊氏彈性模量、管材的線性膨脹系數(shù)、管材密度、比熱容及全線高程。打孔盜油之前，管道穩(wěn)定運(yùn)行條件下全線的油品分布、上下游壓力、流量和溫度；打孔盜油及停輸過程中，打孔盜油管段上下游壓力、流量和溫度隨時(shí)間的變化。停輸過程中閥門的操作過程(關(guān)閥時(shí)間，包括起止時(shí)間)、閥門特性(關(guān)閥過程中對應(yīng)的閥門開度)。管道的外界環(huán)境溫度、油品與外界換熱的總傳熱系數(shù)或通過歷史數(shù)據(jù)反算得到(需已知管道上下游流量、壓力和溫度的歷史數(shù)據(jù))。油品密度、比熱容、黏度與溫度的關(guān)系分別由式(27)-(29)確定，除此之外還應(yīng)給出油品的體積膨脹系數(shù)和壓縮系數(shù)。

圖8 模型求解流程圖Fig. 8 Flow chart of model solving

其中，為油品相對密度，kg/m3；為油品在20 ℃的相對密度kg/m3；ξ為溫度系數(shù)；cy為油品比熱容，kJ/(kg·℃)；T為油品溫度，℃；v1、v2分別為溫度為T1、T2時(shí)油品的黏度，m2/s；u為黏溫指數(shù)，1/℃。表1給出了油品的基礎(chǔ)物性。

表1 20 ℃油品物性Table 1 Physical properties of oil at 20 ℃

3.2 實(shí)際成品油管道打孔盜油分析

3.2.1 管道分輸實(shí)驗(yàn)

(1)實(shí)驗(yàn)過程及實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)依托某成品油管道某管段“YL-GG-LT”，以GG站分輸下載作業(yè)模擬打孔盜油過程，實(shí)驗(yàn)管段全長149.7 km，分輸點(diǎn)位于85.2 km處，沿線地形如圖9所示。

GG站開始分輸后，YL站出站壓力和LT站的進(jìn)站壓力隨時(shí)間的變化趨勢如圖10(a)所示。

圖9 YL-LT管道縱斷面圖Fig. 9 Pipeline route pro fi le between YL and LT

圖10 分輸實(shí)驗(yàn)數(shù)據(jù)及計(jì)算結(jié)果Fig. 10 Experimental data and calculation results of downloading

利用SCADA等現(xiàn)場工況數(shù)據(jù)，泄漏模型計(jì)算的分輸量為5.51 m3(圖10(d))，現(xiàn)場用計(jì)量設(shè)備測定的實(shí)際分輸量為4.97 m3，相對誤差為10.9 %?，F(xiàn)場打孔盜油設(shè)備安裝區(qū)域(分輸實(shí)驗(yàn)中表現(xiàn)為站內(nèi)分輸閥門)開口處不規(guī)則，文獻(xiàn)中尚無將不規(guī)則區(qū)域換算為等效圓孔的方法，本算法和模型將打孔盜油閥門等效為圓孔，計(jì)算得到等效圓孔孔徑為25.4 mm。

由圖10(a)和(b)可以看出，前120 s是穩(wěn)態(tài)數(shù)據(jù)，因此(a)中上游站出站壓力和下游站進(jìn)站壓力以及(c)中盜油點(diǎn)處壓力以準(zhǔn)穩(wěn)態(tài)形式略微升高。120 s后發(fā)生盜油，(c)中盜油點(diǎn)處壓力迅速降低，同時(shí)(c)中盜油點(diǎn)處管內(nèi)流量也迅速降低，形成減壓波分別向上下游傳播。(a)中上下游距離盜油點(diǎn)處距離不同，因此減壓波到達(dá)時(shí)間不同，壓力下降的時(shí)間也不同。盜油發(fā)生后，盜油點(diǎn)處壓力持續(xù)下降，但流量開始從最低點(diǎn)慢慢回升。整個(gè)盜油發(fā)生過程中，上下游和盜油點(diǎn)處壓力持續(xù)降低，但盜油點(diǎn)處管內(nèi)流量從最低點(diǎn)回升至一定值后又逐漸下降。(b)中上游泵站出口溫度和流量始終處于波動(dòng)狀態(tài)，其中溫度變化幅度小，規(guī)律性不強(qiáng)，而流量雖然波動(dòng)但平均值是升高的。(d)中盜油速率在盜油發(fā)生時(shí)最大，之后持續(xù)降低，而累積盜油量一直持續(xù)增加。

3.2.2 現(xiàn)場打孔實(shí)驗(yàn)

(1)打孔實(shí)驗(yàn)1

打孔實(shí)驗(yàn)依托某成品油管道“GY-AS”管段改線工程。實(shí)驗(yàn)將GY與AS之間改線點(diǎn)作為泄放點(diǎn)，實(shí)驗(yàn)管段全長93.9 km，泄放點(diǎn)位于90.5 km處，沿線地形如圖11所示。打孔盜油實(shí)驗(yàn)位置為圖11中的打孔盜油實(shí)驗(yàn)點(diǎn)1。實(shí)驗(yàn)通過軟管將油品泄放至配有流量泵的油罐車內(nèi)以模擬打孔盜油過程(如圖12)，現(xiàn)場通過流量泵計(jì)量油品泄放量。實(shí)驗(yàn)包括不同閥門開度下的3組實(shí)驗(yàn)。

圖11 GY-AS管道縱斷面圖Fig. 11 Route pro fi le of pipeline GY-AS

圖12 GY-AS打孔盜油測試實(shí)驗(yàn)現(xiàn)場圖(左圖為實(shí)驗(yàn)1，右圖為實(shí)驗(yàn)2)Fig. 12 Experiment site of oil stolen by drilling hole

第2組實(shí)驗(yàn)中測試數(shù)據(jù)如圖13(a)和(b)所示，第2組實(shí)驗(yàn)中模型計(jì)算結(jié)果如圖13(c)和(d)所示。

表2中實(shí)驗(yàn)1、2、3閥門開度依次減小，平均盜油速率也依次減小，表明盜油速率與閥門開度成正相關(guān)關(guān)系。第2組實(shí)驗(yàn)過程中，人為干擾較少，因此實(shí)驗(yàn)效果較為理想。模型測算誤差均控制在15%以內(nèi)，說明模型能夠較為精確地測算管道發(fā)生打孔盜油后的盜油量。

(2)打孔實(shí)驗(yàn)2

此次實(shí)驗(yàn)依然基于某成品油管道“GY-AS”管段改線工程。實(shí)驗(yàn)管段全長93.889 km，泄放點(diǎn)位于84.2 km處(如圖11)。打孔盜油實(shí)驗(yàn)位置為圖11中的打孔盜油實(shí)驗(yàn)點(diǎn)2。本次實(shí)驗(yàn)共進(jìn)行了3組實(shí)驗(yàn)。

第4組具體實(shí)驗(yàn)數(shù)據(jù)如圖14中(a)和(b)所示，模型計(jì)算結(jié)果如圖14中(c)和(d)所示。

表3中，實(shí)驗(yàn)4、5、6的閥門開度幾乎相同，結(jié)合圖14(c)、圖15(a)和圖16(a)可以看出，實(shí)驗(yàn)4、5、6盜油處管內(nèi)壓降分別為0.150、0.075、0.095 MPa，平均盜油速率分別為14.70 、6.84 、10.59 m3/h，表明盜油處管內(nèi)壓降越大時(shí)，平均盜油速率也越大，盜油速率與盜油處管內(nèi)壓力成正相關(guān)關(guān)系。同理，盜油處管內(nèi)流量分別為416.08、422.43、417.77 m3/h，盜油速率與盜油處管內(nèi)流量成負(fù)相關(guān)關(guān)系。因此結(jié)合打孔實(shí)驗(yàn)1可以認(rèn)為，對盜油速率影響從大到小的參數(shù)依次是：打孔盜油點(diǎn)孔徑(閥門開度)、打孔盜油處管內(nèi)油品壓力、流量和溫度。

圖13 第2組實(shí)驗(yàn)的數(shù)據(jù)及計(jì)算結(jié)果Fig. 13 Experimental data and calculation results of experiment No.2

表2 GY-AS打孔盜油實(shí)驗(yàn)及計(jì)算結(jié)果數(shù)據(jù)Table 2 Calculation results of experiment No.1, 2 and 3

若將實(shí)驗(yàn)4的數(shù)據(jù)當(dāng)作小孔泄漏計(jì)算泄漏量，從圖17可以看出計(jì)算所得泄漏孔為22.63 mm和泄漏量為16.02 m3，比盜油孔(16.49 mm)和盜油量(4.85 m3)大很多。原因在于孔口自由出流與短管出流的出流系數(shù)差別較大，輸油主管內(nèi)外的壓力分布及壓差也不同，將盜油設(shè)備內(nèi)的壓力損失當(dāng)作小孔出流壓力損失則必然導(dǎo)致計(jì)算所得出流速率偏大，因此打孔盜油過程需要把閥門的局部摩阻、盜油管線摩阻、盜油罐車內(nèi)的油品壓力等考慮在內(nèi)。

研究發(fā)現(xiàn)，本文模型考慮“帶閥管路”模塊后，用于模擬打孔盜油過程，與實(shí)際情況吻合程度較高，模型計(jì)算誤差都在15%以內(nèi)。實(shí)驗(yàn)過程中利用長約100 m的橡膠軟管連接管道平衡孔和油罐車，其產(chǎn)生的摩阻不容忽略，且由于地面不平整，軟管無法避免盤曲現(xiàn)象，增加了局部摩阻，而模型“帶閥管路”模塊算法中計(jì)算的是等長鋼管的摩阻，比實(shí)際軟管摩阻小，因此導(dǎo)致計(jì)算結(jié)果偏大。

圖14 實(shí)驗(yàn)4的數(shù)據(jù)及計(jì)算結(jié)果Fig. 14 Experimental data and calculation results of experiment No. 4

表3 GY-AS打孔盜油實(shí)驗(yàn)及計(jì)算結(jié)果數(shù)據(jù)Table 3 Calculation results of experiment No. 4, 5 and 6

3.3 基于SPS的打孔盜油模型有效性分析

由于大落差管段無法通過現(xiàn)場實(shí)驗(yàn)驗(yàn)證打孔盜油量估算模型的有效性，本文結(jié)合平原和起伏地形條件下成品油管道的盜油情況，基于SPS軟件模擬打孔盜油工況，并對模擬結(jié)果與模型計(jì)算結(jié)果進(jìn)行了對比。SPS軟件建模如圖18所示。

圖15 實(shí)驗(yàn)5的計(jì)算結(jié)果Fig. 15 Calculation results of experiment No. 5

圖16 實(shí)驗(yàn)6的計(jì)算結(jié)果Fig. 16 Calculation results of experiment No. 6

圖17 實(shí)驗(yàn)4當(dāng)作小孔泄漏而計(jì)算泄漏量Fig. 17 Calculation results of experiment No. 4 regarding as small hole leaking process

3.3.1 平原地區(qū)管道打孔盜油

選取某成品油管段相關(guān)數(shù)據(jù)，用分輸工況模擬打孔盜油，利用SPS建模，上下游采取壓力控制，分輸點(diǎn)采取壓力控制。管段全長92 km，打孔盜油點(diǎn)距離首站36 km。取打孔盜油時(shí)間點(diǎn)前120 s作為穩(wěn)態(tài)數(shù)據(jù)反算管道摩阻系數(shù)、管道與環(huán)境的換熱系數(shù)等參數(shù)。管段地形如圖19所示，地形較為平緩，管道沿線最大高程差小于40 m，屬于典型的平原地區(qū)。

設(shè)定5種不同流量下的運(yùn)行工況，分別為150、300、600、1 000、1 300 m3/h，模擬分輸工況1 565 s。圖20為不同流量工況下的模型計(jì)算結(jié)果。

從圖20可以看出，對于平原地區(qū)管道，在設(shè)定的流量變化范圍中，模型計(jì)算結(jié)果與SPS計(jì)算結(jié)果的相對偏差依次為：1.69%、0.38%、4.43%、7.65%、8.92%，均在15%以內(nèi)，最大相對偏差僅為8.92%，特別是在流量較小的情況下，模型計(jì)算結(jié)果與SPS模擬結(jié)果更為接近，相對偏差僅為0.38%。因此，在打孔盜油速率較小時(shí)，模型的工況模型計(jì)算精確度更高。對于平原地區(qū)管道，相對誤差隨著輸量增大而增大。

3.3.2 起伏地區(qū)管道打孔盜油

選取某起伏地區(qū)管段相關(guān)數(shù)據(jù)(圖21)，用分輸工況模擬打孔盜油工況，利用SPS建模，上下游采取壓力控制，分輸點(diǎn)采取壓力控制。管段全長215.404 km，盜油點(diǎn)距離首站100.09 km，取打孔盜油時(shí)間點(diǎn)前120 s作為穩(wěn)態(tài)反算數(shù)據(jù)。

設(shè)定5種不同流量下的運(yùn)行工況，分別為195，300，450，600，750 m3/h模擬時(shí)長為1 565s的分輸工況。圖22為不同流量工況下，模型計(jì)算結(jié)果。

從圖22可以看出，對于存在翻越點(diǎn)的管道，在設(shè)定的流量范圍內(nèi)進(jìn)行模擬計(jì)算時(shí)，模型計(jì)算與SPS模擬相對誤差依次為9.77%、14.56%、7.69%、1.77%、0.16%，均在15%以內(nèi)；尤其在流量較大時(shí)，模型計(jì)算結(jié)果與SPS模擬結(jié)果較為接近，即在打孔盜油速率較大的情況下，模型的工況模型計(jì)算精確度更高。此外，結(jié)合相對誤差分析可知，對于起伏地區(qū)管道，相對誤差隨著輸量增大而減小。

圖18 SPS建模示意圖Fig. 18 Illustration of SPS model

圖19 某平原地區(qū)管道里程及高程數(shù)據(jù)Fig. 19 Pipeline route pro fi le in plain area

圖20 平原地區(qū)管道盜油量測算Fig.20 The calculation result of stolen oil volume of pipeline in fl atlands

4 結(jié)論

針對成品油管道打孔盜油量預(yù)測問題建立了與之對應(yīng)的水熱力耦合模型，利用特征線法、有限差分法對模型及參數(shù)進(jìn)行了處理，考慮各類邊界條件，對非穩(wěn)態(tài)打孔盜油過程進(jìn)行了數(shù)值模擬研究，得出了打孔盜油點(diǎn)孔徑預(yù)測的方法和打孔盜油量的估算方法。通過理論研究和數(shù)值計(jì)算，主要得到以下結(jié)論：

圖21 某起伏地區(qū)管道里程及高程數(shù)據(jù)Fig. 21 Pipeline route pro fi le in mountain area

圖22 起伏地區(qū)管道盜油量測算Fig. 22 Calculation results of the stolen oil amount of pipeline in mountain area

(1)構(gòu)建和求解描述成品油管道的非穩(wěn)態(tài)水熱力瞬變過程的數(shù)學(xué)模型，可以得出瞬時(shí)打孔盜油量和累計(jì)打孔盜油量隨時(shí)間的變化規(guī)律。穩(wěn)態(tài)和非穩(wěn)態(tài)的水熱力耦合程度不同，在進(jìn)行打孔盜油量計(jì)算之前，應(yīng)使管道運(yùn)行參數(shù)達(dá)到真正的穩(wěn)定狀態(tài)。

(2)打孔盜油過程與小孔泄漏過程不完全相同。原因在于孔口自由出流和短管出流的出流系數(shù)差別較大，出流處管內(nèi)外壓差不同，因此打孔盜油過程需要考慮閥門的局部摩阻、盜油管線摩阻和盜油罐車內(nèi)的油品壓力。

(3)現(xiàn)場打孔盜油設(shè)備安裝區(qū)域不規(guī)則，目前尚無將不規(guī)則區(qū)域換算為等效圓孔的方法。本算法和模型將打孔盜油閥門等效為圓孔，模型可根據(jù)上下游壓力變化反算出此當(dāng)量孔徑；反算孔徑后，模型可進(jìn)一步計(jì)算非穩(wěn)態(tài)和穩(wěn)態(tài)時(shí)的瞬時(shí)打孔盜油速率和累計(jì)盜油量。

(4)穩(wěn)態(tài)和非穩(wěn)態(tài)打孔盜油過程中，對打孔盜油速率影響從大到小的參數(shù)依次是：打孔盜油點(diǎn)等效孔徑(閥門開度)、打孔盜油處管內(nèi)油品壓降、流量和溫度。

(5)本文模型既可用于平原地區(qū)管段打孔盜油量測算，也可用于起伏地區(qū)管段打孔盜油量測算。分析發(fā)現(xiàn)，輸量越小、壓力越高、打孔盜油速率越大、壓力波動(dòng)程度越大，模型計(jì)算結(jié)果精度越高。輸量越大、壓力越低、打孔盜油速率越小、壓力波動(dòng)程度越小，模型計(jì)算結(jié)果精度越低。

[1] WANG S X, JOHN J C. Leak detection for gas and liquid pipelines by transient modeling[C]. International Oil & Gas Conference and Exhibition in China. Society of Petroleum Engineers, 2006.

[2] ZHANG J, HOFFMAN A, MURPHY K, et al. Review of pipeline leak detection technologies[C]. PSIG Annual Meeting. Pipeline Simulation Interest Group, 2013.

[3] WOODARD J L, MUDAN K S. Liquid and gas discharge rates through holes in process vessels[J]. Journal of Loss Prevention in the Process Industries, 1991, 4(3): 161-5.

[4] 孫天宇. 打孔盜油支管對輸送成品油管內(nèi)壓力分布的影響[J]. 科學(xué)技術(shù)與工程, 2012, 20(8): 1906-1908.[SUN T Y. Study on the in fl uence of slotting and stealing branch pipe on the pressure fi eld distribution in oil product transportation pipeline[J]. Science Technology and Engineering, 2012, 20(8): 1906-1908.]

[5] 李大全, 姚安林.成品油管道泄漏擴(kuò)散規(guī)律分析[J]. 油氣儲(chǔ)運(yùn), 2006, 25(8): 18-24.[LI D Q,YAO A L. The study on oil-diffusing law of product pipeline leaking[J]. Oil & Gas Storage and Transportation, 2006, 25(8): 18-24.]

[6] NAZMUL Z S. Dynamic seals lower life costs of wastewater pumps[J]. Sealing Technology, 2004, 8(8): 11-12.

[7] 聶平, 鄧松圣, 文俊. 管道泄漏水力瞬變模擬研究[J]. 油氣儲(chǔ)運(yùn), 2006, 25(7): 19-22.[NIE P, DENG S S, WEN J. Study on the leak in the hydraulic transient simulation of the pipeline[J]. Oil & Gas Storage and Transportation, 2006, 25(7): 19-22.]

[8] 劉恩斌, 彭善碧, 李長俊, 等. 基于瞬態(tài)模型的油氣管道泄漏檢測[J]. 天然氣工業(yè), 2005, 25(6): 102-104.[LIU E B, PENG S B, LI C J, et al. Leakage detection of oil and gas pipeline by transient model[J]. Natural Gas Industry, 2005, 25(6): 102-104.]

[9] 孫良, 王建林. 基于泄漏瞬變模型的管道泄漏檢測與定位方法[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào), 2012, 20(1): 159-168.[SUN L, WANG J L. Release transient modeling based leak detection and location method for pipelines[J]. Journal of Basic Science and Engineering, 2012, 20(1): 159-168.]

[10] SALUFU S O, IBUKUN S. Prediction of rate and volume of oil spill in horizontal and inclined pipelines[C]. Nigeria Annual International Conference and Exhibition. 2013.

[11] OBIBUIKE U J, OBAH B, OGWO O U J, et al. A novel approach to estimation of leak volume in an oil pipeline[J]. Petroleum & Coal, 2015, 57(2): 154-168.

[12] MODE A W, AMOBI J O, SALUFU S O. A model for predicting rate and volume of oil spill in horizontal and vertical pipelines[J]. Journal of Environment and Earth Science, 2013, 3(9): 12-19.

[13] 王新穎, 邵輝, 王宏鑫. 基于量綱分析的輸液管道泄漏量計(jì)算模型研究[J]. 工業(yè)安全與環(huán)保, 2010, 36(4): 43-44.[WANG X Y, SHAO H, WANG H X. Research on the calculation model of pipeline leakage rate based dimension analysis[J]. Industrial Safety and Environmental Protection, 2010, 36(4): 43-44.]

[14] REED M, EMILSEN M H, HETLAND B, et al. Numerical model for estimation of pipeline oil spill volumes[J]. Environmental Modelling & Software, 2006, 21(2): 178-189.

[15] 劉恩斌, 彭善碧, 李長俊, 等. 基于瞬態(tài)模型的油氣管道泄漏檢測[J]. 天然氣工業(yè), 2005, 25(6): 102-104.[LIU E B, PENG S B, LI C J, et al. Leakage detection of oil and gas pipeline by transient model[J]. Natural Gas Industry, 2005, 25(6): 102-104.]

[16] ZHU H, LIN P, PAN Q. A CFD (Computational Fluid Dynamic) simulation for oil leakage from damaged submarine pipeline[J]. Energy, 2014, 64(1): 887-899.

[17] KIM B I, SHARMA M P, HARRIS H G. A statistical approach for predicting volume of oil spill during pipeline operations[C]. SPE Annual Technical Conference and Exhibition. 1991.

[18] RYAN K. GIS Method for calculating maximum potential spill volume due to natural landforms[M]. Saint Mary’s University of Minnesota Central Services Press, Winona, MN, 2008(4): 100-110.

[19] 劉國華, 趙立奇. 輸油管道泄漏損失量的評估方法[J]. 油氣儲(chǔ)運(yùn), 2013, 32(6): 672-674.[LIU G H, ZHAO L Q. Assessment of leakage volume in oil pipeline[J]. Oil & Gas Storage and Transportation, 2013, 32(6): 672-674.]

[20] 邵輝, 張蓉愛, 王新穎, 等. 氣、液管道泄漏實(shí)驗(yàn)裝置的設(shè)計(jì)[J]. 江蘇工業(yè)學(xué)院學(xué)報(bào), 2008, 20(2): 45-48.[SHAO H, ZHANG R A, WANG X Y. Study of experimental equipment deisigns for leakage of gas and liquid pipelines[J]. Journal of Jiangsu Polytechnic University, 2008, 20(2): 45-48.]

[21] JASPER A A, OVE G T, TORLEIV B. Experimental study of oil pipeline leak processes[J]. Journal of Environmental Protection, 2012, 18(03): 597-604.

[22] 付建民, 趙振洋, 陳國明, 等. 液相管道流量與壓力對小孔泄漏速率的影響[J]. 石油學(xué)報(bào), 2016, 37(2): 257-265.[FU J M, ZHAO Z Y, CHEN G M, et al. Influences of liquid pipeline flow and pressure on small-hole leakage rate[J]. Acta Petrolei Sinica, 2016, 37(2):257-265.]

[23] 孫良, 王建林, 趙利強(qiáng). 負(fù)壓波法在液體管道上的可檢測泄漏率分析[J]. 石油學(xué)報(bào), 2010, 31(04): 654-658.[SUN L, WANG J L, ZHAO L Q. Analysis on detectable leakage ratio of liquid pipeline by negative pressure wave method[J]. ACTA PETROLEI SINICA, 2010, 31(04): 654-658.]

[24] 孫良, 王建林. 基于泄漏瞬變模型的管道泄漏檢測與定位方法[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào), 2012, 20(1): 159-168.[SUN L, WANG J L. Release transient modeling based leak detection and location method for pipelines[J]. TOURNAL OF BASIC SCIENCE AND ENGINEERING, 2012, 20(1): 159-168.]

Calculation methods for the amount of oil stolen by drilling holes in product pipelines

HE Guoxi1, LIANG Yongtu1, LI Yansong1, LIU Shengli1, WU Mengyu1, XIE Cheng2, XI Gang2, LI Feng3
1 Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum-Beijing, Beijing 102249, China
2 SINOPEC Sales Company South China Branch, Guangzhou 510620, China
3 Sinopec Beihai Lique fi ed Natural Gas Limited Liability Company, Beihai 536000, China

Stealing oil by drilling hole brings great harm, including the oil loss, the negative influence on normal pipeline operation and risks of potential explosion, leakage and environmental pollution. At present, many researchers aim to identify the problems and have put forward some preventive measures, but few focus on the calculation of the loss. The transient hydraulic properties of the whole pipeline system during the process of oil theft by hole drilling is analyzed. Considering boundary conditions of the oil-stealing point, diameter-changing point, batch interface and valves, the model is established. Based on the pressure, fl ow rate and temperature, the transient and accumulated volume of stolen oil are given. The results show that: (1) the ori fi ce free fl ow is different from short-pipe out fl ow, therefore the valves friction, the friction of the short oil-stealing pipe and the inside pressure of the oil-stealing tankers should be taken into in calculation; (2) The parameters in fl uencing the amount of oil stolen can be ranked from large to small, namely the size of ori fi ce, the pressure, fl ow rate and temperature; (3) The model has relatively better accuracy when the fl ow rate is lower, the pressure is higher and the pressure fl uctuation is higher.

products pipeline; oil stolen by hole drilling; hydraulic and thermal transient process; calculation for the amount of oil stolen

10.3969/j.issn.2096-1693.2017.01.009

(編輯 馬桂霞)

*通信作者, liangyt21st@163.com

2016-11-09

國家自然科學(xué)基金項(xiàng)目“成品油管道批次輸送過程中的復(fù)雜傳熱傳質(zhì)機(jī)理研究”(51474228)、北京市科學(xué)研究與研究生培養(yǎng)共建項(xiàng)目科研項(xiàng)目“成品油管道泄漏量預(yù)測及停輸方案的優(yōu)化研究”(ZX20150440)聯(lián)合資助

何國璽, 梁永圖, 李巖松, 劉勝利, 吳夢雨, 謝成, 席罡, 李豐. 成品油管道打孔盜油量測算方法. 石油科學(xué)通報(bào), 2017, 01: 86-101

HE Guoxi, LIANG Yongtu, LI Yansong, LIU Shengli, WU Mengyu, XIE Cheng, XI Gang, LI Feng. Calculation methods for the amount of oil stolen by drilling holes in product pipelines. Petroleum Science Bulletin, 2017, 01: 86-101. doi: 10.3969/ j.issn.2096-1693.2017.01.009