基于Runge-Kutta法的輸氣管道泄漏定位技術(shù)研究

吳波 楊旭東 張小芳 司宗慶 牟燕 梁昌晶

1中國(guó)石油華北油田公司第三采油廠

2中國(guó)石油華北油田公司第一采油廠

3中國(guó)石油華北油田公司二連分公司

輸氣管道的泄漏定位技術(shù)是管道完整性管理的重要組成部分，管道防腐層的缺失、陰極保護(hù)的失效、土壤介質(zhì)對(duì)管道外壁的腐蝕，以及管輸介質(zhì)對(duì)管道內(nèi)壁產(chǎn)生的腐蝕，均會(huì)引起管道不同程度的泄漏[1]。據(jù)統(tǒng)計(jì)，1995—2005 年，中國(guó)石油管道分公司每年大約產(chǎn)生12～15 起重大輸氣管道泄漏事故，其中2004 年因腐蝕引起的管道泄漏事故就多達(dá)20余起[2-3]。目前，已有多種針對(duì)管道泄漏的檢測(cè)與定位技術(shù)，但每種檢測(cè)方法的原理和手段差別較大，同時(shí)由于管輸天然氣組分的多樣性，敷設(shè)方式的差異性（如埋地、直敷、架空、穿跨越等），以及泄漏方式的多元性（如滲透、穿孔、斷裂等），使得各種檢測(cè)方法在定位精度和檢測(cè)時(shí)間上都存在一定的局限性[4-5]。

韓玲娟等[6]采用分布式光纖傳感器對(duì)水下輸氣管道的泄漏進(jìn)行了有效定位，采用小波神經(jīng)網(wǎng)絡(luò)對(duì)信號(hào)進(jìn)行了有效提取，用聚類分析判斷管道泄漏的嚴(yán)重程度，并采用最小二乘法擬合信號(hào)參數(shù)進(jìn)行定位，該方法與傳統(tǒng)定位方法相比誤報(bào)率下降了6.6%；焦敬品等[7]通過對(duì)20種不同的泄漏信號(hào)進(jìn)行識(shí)別，運(yùn)用BP 神經(jīng)網(wǎng)絡(luò)進(jìn)行泄漏點(diǎn)的檢測(cè)和定位，整體識(shí)別率≥92.5%；田野[8]設(shè)計(jì)了一套基于次生波原理的泄漏監(jiān)測(cè)系統(tǒng)，結(jié)合了比例-積分-微分（PID）整定調(diào)節(jié)和自適應(yīng)算法，對(duì)輸氣管道的泄漏定位誤差≤20 m；李鳳等[9]將負(fù)壓波和聲波技術(shù)進(jìn)行了耦合，通過前處理軟件（ICEM）對(duì)管道進(jìn)行了網(wǎng)格劃分，應(yīng)用計(jì)算流體力學(xué)（CFD）軟件進(jìn)行了流場(chǎng)變化分析，有效解決了單一信號(hào)噪聲大、魯棒性差的缺點(diǎn)。綜上所述，管道的定位檢測(cè)需要在首末點(diǎn)設(shè)置大量的傳感器，并且采用多種軟件進(jìn)行模型構(gòu)建和數(shù)據(jù)模擬，工作量大，一般操作人員較難掌握。

我國(guó)大部分輸氣管道由于管道沿程高差的區(qū)別和地理跨度的差異，屬于非等溫輸氣管道，隨著計(jì)算機(jī)科學(xué)和數(shù)值分析理論的不斷發(fā)展，通過現(xiàn)有傳感器測(cè)量的溫度、壓力和流量數(shù)據(jù)，采用改進(jìn)的變步長(zhǎng)龍格-庫(kù)塔（Runge-Kutta）法對(duì)非等溫輸氣管道進(jìn)行檢測(cè)和定位，同時(shí)通過實(shí)例進(jìn)行仿真計(jì)算，并與傳統(tǒng)的負(fù)壓波檢測(cè)結(jié)果進(jìn)行對(duì)比，以期可以快速、準(zhǔn)確地對(duì)泄漏點(diǎn)進(jìn)行定位。

1 模型建立

1.1 輸氣管道模型

將輸氣管道中介質(zhì)的流動(dòng)視為穩(wěn)態(tài)流動(dòng)，而穩(wěn)態(tài)流動(dòng)與時(shí)間變量無關(guān)，根據(jù)流體力學(xué)和熱力學(xué)相關(guān)公式建立沿管輸方向的穩(wěn)態(tài)流動(dòng)方程[10-11]。其中連續(xù)性方程、動(dòng)量方程、能量方程和狀態(tài)方程分別為

式中：ρ為管輸氣體密度，kg/m3；ω為氣體流速，m/s；x為管道長(zhǎng)度，m；p為管輸氣體壓力，Pa；g為重力加速度，取9.8 m/s2；θ為管道傾角，rand；λ為水力摩阻系數(shù)；d為管道內(nèi)徑，m；h為管輸氣體的焓，J/kg；s為管道的位置高度，m；Q為單位質(zhì)量氣體向外界釋放的能量，J/kg；T為管輸氣體溫度，K。

在已知?dú)怏w組分及溫度、壓力的前提下，可以由實(shí)際氣體狀態(tài)方程對(duì)氣體密度和壓縮系數(shù)進(jìn)行求解，比較實(shí)用的狀態(tài)方程有RK、SRK、BWRS 和PR方程，在此采用BWRS方程進(jìn)行計(jì)算[12]。

1.2 Runge-Kutta法

該算法是一種在工程上被廣泛應(yīng)用的高精度單步算法，主要原理是通過對(duì)常微分方程組進(jìn)行離散化處理，求解該方程的初值；隨后選擇不同的步長(zhǎng)，在之間計(jì)算多個(gè)點(diǎn)的斜率值，將其加權(quán)平均作為平均斜率。為了便于計(jì)算，通常采用較為經(jīng)典的四階Runge-Kutta 進(jìn)行計(jì)算[13]，通過與泰勒級(jí)數(shù)展開式進(jìn)行對(duì)比后求解未知數(shù)，簡(jiǎn)化后的公式如下

式中：xi、yi分別為自變量和因變量；yi+1為下一個(gè)步長(zhǎng)的因變量；h為步長(zhǎng)系數(shù)；K1、K2、K3、K4為待定系數(shù)，表示輸出變量的一階倒數(shù)。

將管道分成n段，每段管道長(zhǎng)度定為Δx（代表每段管道的微元）。Runge-Kutta法的關(guān)鍵是對(duì)于Δx的合理選擇，單從每一個(gè)截?cái)鄟砜?，步長(zhǎng)越小，截?cái)嗾`差越小，但隨著步長(zhǎng)的不斷減少，計(jì)算時(shí)間也在不斷延長(zhǎng)，一定范圍內(nèi)的求解步數(shù)也會(huì)增加。不同的迭代計(jì)算會(huì)導(dǎo)致舍入誤差的大量累計(jì)，在計(jì)算效率和精度上有很大缺陷，因此首先設(shè)置較大的計(jì)算步長(zhǎng)，在整個(gè)管段內(nèi)進(jìn)行泄漏位置的初選，隨后逐步縮小泄漏定位的搜索范圍，直到滿足定位精度的要求，一般取0.5%～1%的初始步長(zhǎng)進(jìn)行運(yùn)算。

在經(jīng)典四階Runge-Kutta 算法的基礎(chǔ)上采用精細(xì)Runge-Kutta 法計(jì)算，該方法融合了精細(xì)積分法和Runge-Kutta 法的優(yōu)點(diǎn)，可以將特解的數(shù)值代入Runge-Kutta公式中積分，并利用Runge-Kutta法對(duì)數(shù)值積分點(diǎn)的未知參數(shù)進(jìn)行預(yù)估，得到更精細(xì)的算法。根據(jù)現(xiàn)場(chǎng)傳感器采集到的管道起點(diǎn)溫度Ti、壓力pi和流量mi，計(jì)算從起點(diǎn)到末點(diǎn)管道各截面的溫度、壓力和流量數(shù)值。首先計(jì)算第一個(gè)微元Δx1的斜率k1，從第一個(gè)積分點(diǎn)出發(fā)，前進(jìn)半個(gè)步長(zhǎng)，得到第二個(gè)積分點(diǎn)的預(yù)測(cè)斜率k2；然后以k2為基礎(chǔ)，再次從第一個(gè)積分點(diǎn)出發(fā)，前進(jìn)半個(gè)步長(zhǎng)，重新計(jì)算第二個(gè)積分點(diǎn)的預(yù)測(cè)斜率k′2；將兩次計(jì)算的k2和k′2進(jìn)行加權(quán)平均，作為第二個(gè)積分點(diǎn)的近似斜率；以此類推，直到計(jì)算出最后一個(gè)積分點(diǎn)也就是管道末點(diǎn)的溫度、壓力和流量數(shù)值。其中以第一個(gè)微元管道Δx1為例，運(yùn)用Runge-Kutta法求解公式（1）～（4）來計(jì)算管道截面的參數(shù)，公式如下

以此類推，得到管道末點(diǎn)的溫度Ti+1、壓力pi+1和流量mi+1的計(jì)算公式如下

從管道起點(diǎn)開始，運(yùn)用Matlab 對(duì)公式（7）進(jìn)行反復(fù)的迭代計(jì)算，即可完成從管道起點(diǎn)到末點(diǎn)不同截面的溫度、壓力和流量的計(jì)算；從末點(diǎn)到起點(diǎn)的參數(shù)計(jì)算也類似，只是將公式（2）～（3）的摩擦阻力項(xiàng)和熱損失項(xiàng)進(jìn)行調(diào)整，具體見式（8）～（9），這種模型非常適合非等溫輸氣管道的模擬和仿真。

2 泄漏定位算法

通過現(xiàn)場(chǎng)傳感器采集到的起點(diǎn)溫度Ti、壓力pi和流量mi，依據(jù)Runge-Kutta法計(jì)算末點(diǎn)管道的溫度、壓力和流量，同樣通過采集到的末點(diǎn)溫度Tn、壓力pn和流量mn，反算起點(diǎn)管道的溫度、壓力和流量，得到圖1。

圖1 泄漏定位原理Fig.1 Leak location principle

假設(shè)在某一時(shí)刻t，由管道末點(diǎn)計(jì)算到管道起點(diǎn)的參數(shù)為Tnt、pnt、mnt，將判斷泄漏事故發(fā)生的壓力閾值定義為εt，則判斷泄漏的邊界條件為

當(dāng)公式（10）成立，表示管內(nèi)無泄漏點(diǎn)；當(dāng)存在泄漏點(diǎn)時(shí)，管道內(nèi)各點(diǎn)的相關(guān)參數(shù)會(huì)偏離計(jì)算結(jié)果，式（10）不再成立。當(dāng)pnt＜pn時(shí)，管道內(nèi)存在泄漏點(diǎn)。在管道的運(yùn)行過程中，難免會(huì)出現(xiàn)關(guān)停壓縮機(jī)或停輸再啟動(dòng)等工況，此時(shí)的監(jiān)測(cè)數(shù)據(jù)會(huì)有較大波動(dòng)，容易產(chǎn)生誤報(bào)，因此不考慮改變操作引起的擾動(dòng)在時(shí)間內(nèi)的報(bào)警，t1為泄漏點(diǎn)傳播到起點(diǎn)的時(shí)間，t2為泄漏點(diǎn)傳播到末點(diǎn)的時(shí)間，t3為擾動(dòng)的熵變傳播到末點(diǎn)的時(shí)間，n為劃分管道的段數(shù)，Δt為報(bào)警閾值。

圖1 中L1代表未發(fā)生泄漏時(shí)的壓降，L2代表泄漏發(fā)生后從起點(diǎn)計(jì)算到末點(diǎn)的壓降，L3代表泄漏發(fā)生后從末點(diǎn)計(jì)算到起點(diǎn)的壓降。當(dāng)發(fā)生泄漏時(shí)，泄漏點(diǎn)處氣體的密度和壓縮系數(shù)會(huì)發(fā)生變化，壓力也會(huì)發(fā)生瞬變。對(duì)于L2，泄漏點(diǎn)之前的仿真結(jié)果無限接近L1，泄漏點(diǎn)之后的仿真結(jié)果偏離L1，且距離x0越遠(yuǎn)，偏離程度越嚴(yán)重；對(duì)于L3，泄漏點(diǎn)之后的仿真結(jié)果接近L1，泄漏點(diǎn)之前的仿真結(jié)果偏離L1，L2、L3的交點(diǎn)為x0即為泄漏點(diǎn)。泄漏定位流程圖見圖2。

3 實(shí)例計(jì)算

華北油田采油二廠從別古莊到霸一聯(lián)的輸氣管道全長(zhǎng)2.7 km，管道規(guī)格為Φ213 mm×5.5 mm，管道輸量200 m3/min，實(shí)際運(yùn)行壓力0.4～0.8 MPa，管輸介質(zhì)成分CH4占98%，另外還有少量的C2H6、C3H8等。利用管線全面清管期間，在一些特定位置（充分考慮安全、環(huán)保的要求）人工制造幾處泄漏點(diǎn)，通過球閥控制泄漏過程。為了避免放氣時(shí)間與信號(hào)采集不同步，先進(jìn)行閥門放氣，一段時(shí)間后關(guān)閉閥門，設(shè)置首站壓力0.65 MPa。圖3 和圖4 分別為泄漏點(diǎn)位于40.5 m時(shí)，發(fā)生泄漏后管道起點(diǎn)和末點(diǎn)的壓力隨時(shí)間變化的曲線。

圖2 泄漏定位流程圖Fig.2 Flow chart of leak location

圖3 起點(diǎn)壓力隨時(shí)間變化曲線Fig.3 Change curve of start point pressure with time

圖4 末點(diǎn)壓力隨時(shí)間變化曲線Fig.4 Change curve of final point pressure change with time

由圖3 和圖4 可知，在15 s 處管道發(fā)生泄漏，15～25 s 的時(shí)間段內(nèi)首末點(diǎn)壓力急劇下降，而后趨于平穩(wěn)，可見整個(gè)泄漏過程經(jīng)歷了穩(wěn)態(tài)—快瞬態(tài)—慢瞬態(tài)—穩(wěn)態(tài)四個(gè)過程，流體從恒定流狀態(tài)向瞬變流狀態(tài)轉(zhuǎn)變，穩(wěn)定之后的壓力值低于泄漏前的壓力值。為了更好地驗(yàn)證改進(jìn)后的Runge-Kutta 法對(duì)泄漏定位預(yù)測(cè)的準(zhǔn)確性，控制閥門開度，以不同的泄漏率（即泄漏孔徑與管徑的比值）對(duì)模型進(jìn)行驗(yàn)證，同時(shí)與傳統(tǒng)的負(fù)壓超聲波定位精度進(jìn)行對(duì)比，結(jié)果如表1所示。

表1 不同泄漏率下泄漏定位結(jié)果對(duì)比Tab.1 Comparison of leakage location results under different leakage rates

由表1可知，在泄漏位置一定的情況下，隨著泄漏孔徑和泄漏率的增加，Runge-Kutta 法和負(fù)壓波定位的定位誤差都在逐漸變小，通過結(jié)果對(duì)比，Runge-Kutta 法的定位誤差為0.4～2.5 m，小于傳統(tǒng)的負(fù)壓波定位誤差（5～12 m）。

地形地貌及管輸介質(zhì)成分的不同決定了實(shí)際工況下的輸氣管道大多數(shù)為非等溫管道，為考察首末段溫差對(duì)定位精度的影響，在泄漏位置和泄漏率一定的前提下，通過提高首站和末站的溫度，在泄漏率5%的條件下用Runge-Kutta法和負(fù)壓波對(duì)泄漏進(jìn)行定位，結(jié)果如表2所示。

由表2可知，隨著首末點(diǎn)溫差的不斷增加，兩種方法的定位誤差也隨之增大，但Runge-Kutta 法在定位精度上明顯優(yōu)于負(fù)壓波定位，因此在溫差較大的情況下，Runge-Kutta 法可用于非等溫輸氣管道的泄漏檢測(cè)和定位。

管道輸量的大小對(duì)于不同泄漏孔徑也有一定影響?？刂剖渍玖髁繛?50～300 m3/min，在泄漏率5%條件下，首站壓力設(shè)置為0.5～0.7 MPa，對(duì)326 m處的泄漏點(diǎn)進(jìn)行定位，結(jié)果如表3所示。

表2 首末點(diǎn)溫差變化對(duì)定位精度的影響Tab.2 Influence of temperature difference between start and final point on positioning accuracy

表3 不同管道輸量對(duì)定位精度的影響Tab.3 Influence of different pipeline volumes on positioning accuracy

由表3可知，隨著管道輸量的不斷增加，定位誤差整體減小。這是因?yàn)檩斄吭酱螅艿纼?nèi)的氣體流速越大，在泄漏孔處的流出速度越大，壓力波動(dòng)越明顯。

綜上所述，對(duì)比不同情況下的定位誤差，對(duì)定位精度的影響因素排序由大至小是管道輸量、泄漏率、首末點(diǎn)溫差。

4 結(jié)束語(yǔ)

（1）通過采集到的首末點(diǎn)溫度、壓力和流量數(shù)據(jù)，應(yīng)用改進(jìn)后的Runge-Kutta 法可以計(jì)算管道各點(diǎn)的壓力分布，進(jìn)而對(duì)泄漏點(diǎn)進(jìn)行有效定位，該方法原理簡(jiǎn)單，對(duì)現(xiàn)場(chǎng)操作人員的專業(yè)要求不高，易于實(shí)現(xiàn)，與傳統(tǒng)的負(fù)壓波相比具有一定的優(yōu)越性。

（2）泄漏定位技術(shù)涉及多學(xué)科的融合和交叉，不應(yīng)只注重現(xiàn)場(chǎng)傳感器的數(shù)據(jù)，隨著機(jī)器學(xué)習(xí)、數(shù)值分析的發(fā)展，今后應(yīng)加強(qiáng)各學(xué)科之間的融合，提高定位精度。

（2）該方法只適用于干氣管道，對(duì)于含液較多的管道由于氣液流速不一致、氣液滑脫等現(xiàn)象并不適用，同時(shí)對(duì)于同一時(shí)刻有多處泄漏點(diǎn)的數(shù)據(jù)擾動(dòng)考慮不足，今后應(yīng)深入研究，加強(qiáng)模型的完善以提高定位精度。