呂玉坤，程 博，趙 鍇，彭 鑫，李 聰

(1．華北電力大學(xué) 能源與動力工程學(xué)院，河北 保定 071003;2．杭州和利時自動化有限公司，陜西 西安 710075;3．華能荊門熱電有限責(zé)任公司，湖北 荊門 448002)

0 引言

近年來，隨著計算機(jī)技術(shù)、信息技術(shù)、空間技術(shù)以及人工智能技術(shù)等的不斷發(fā)展[1，2]，熱力管網(wǎng)故障診斷方法研究工作取得了突破性的進(jìn)步。就總體而言，基于硬件的故障診斷技術(shù)雖然靈敏性好、定位準(zhǔn)確，但由于其成本和施工費(fèi)用高、實際管網(wǎng)發(fā)生泄漏情況無法預(yù)知，而且無法實時監(jiān)測每個測點(diǎn)，故該類技術(shù)一般不作為熱力管網(wǎng)故障診斷的主要手段[3，4];而基于軟件的故障診斷技術(shù)其主要優(yōu)點(diǎn)在于:可實時在線實現(xiàn)故障診斷、及時給出報警信號。因此，開展以軟件為主、軟硬件相結(jié)合的熱力管網(wǎng)故障診斷方法的研究工作，以實現(xiàn)優(yōu)勢互補(bǔ)、提高熱力管網(wǎng)管理的自動化水平，從而保證管網(wǎng)安全生產(chǎn)和經(jīng)濟(jì)運(yùn)行已成必然之勢。

筆者以國際上公認(rèn)的管道泄漏檢測與定位系統(tǒng)的9項性能指標(biāo)為依據(jù)，對各種故障診斷方法進(jìn)行綜合比較后發(fā)現(xiàn):基于狀態(tài)估計的實時模型法[5，6]具有可很好地解決管道內(nèi)流體特性參數(shù)隨時間變化的特點(diǎn)，但其對流量計的要求高、誤報警率高，且僅適用于小泄漏的管網(wǎng)故障診斷。有鑒于此，本文提出了適用于枝狀熱力管網(wǎng)單點(diǎn)泄漏的故障診斷方法——插入虛擬節(jié)點(diǎn)法。

插入虛擬節(jié)點(diǎn)法的基本思想是:假設(shè)管網(wǎng)某點(diǎn)泄漏后通過補(bǔ)水能夠維持水壓圖達(dá)到新的平衡，把泄漏點(diǎn)視為一個新的熱用戶，即將泄漏點(diǎn)看作新插入的虛擬輸出節(jié)點(diǎn)，則整個管網(wǎng)將仍滿足流量守恒方程。

為實現(xiàn)實時在線故障診斷，與插入虛擬節(jié)點(diǎn)法相應(yīng)地定義了一個管網(wǎng)報警閥的概念。以某市南市區(qū)的枝狀熱力管網(wǎng)從某正常運(yùn)行狀態(tài)過渡到故障狀態(tài)為研究對象，采用圖論和線性代數(shù)相結(jié)合的方法建立其熱力管網(wǎng)水力計算模型，并應(yīng)用環(huán)方程法求解之，以驗證所提方法和所設(shè)管網(wǎng)報警閥值域的合理性。

1 熱力管網(wǎng)水力模型與求解

圖1為某市南市區(qū)的枝狀熱力管網(wǎng)圖。供熱總面積為1210.38萬m2，單一熱源廠供熱，熱源廠出口壓力為0.5 MPa，供回水溫度分別為120℃和60℃，水的平均密度取983.3 kg/m3，管材為無縫鋼管，絕對粗糙度Δ=0.01 mm，供暖期為142天。

圖1 某市南市區(qū)的枝狀熱力管網(wǎng)圖Fig．1 Nanshi district of branched heat pipe network diagram

通過將實際熱力管網(wǎng)簡化和抽象為管線和節(jié)點(diǎn)兩類組成元素，體現(xiàn)出它們之間的水力特性和拓?fù)潢P(guān)系，以便于圖形、數(shù)據(jù)的表達(dá)和分析，為后續(xù)的水力計算、工況分析和故障診斷等工作奠定基礎(chǔ)。熱力管網(wǎng)數(shù)學(xué)模型本質(zhì)上就是由管線和節(jié)點(diǎn)組成的一個有序的集合體。因此，可采用圖論和線性代數(shù)相結(jié)合的方法，系統(tǒng)地建立其熱力管網(wǎng)數(shù)學(xué)模型并求解之[7]。

1.1 熱力管網(wǎng)水力模型的建立

根據(jù)實際工程的要求，為便于進(jìn)行管網(wǎng)網(wǎng)絡(luò)方程的建立，做出以下5個假設(shè):

(1)管網(wǎng)中工質(zhì)的流動為不可壓縮流動;

(2)準(zhǔn)穩(wěn)態(tài)，熱力管網(wǎng)特性參數(shù)(流量、壓力等)不隨時間變化;

(3)流動處于阻力平方區(qū)，ΔH=Sq2;

(4)動壓忽略不計;

(5)熱力管網(wǎng)的管線數(shù)為M、節(jié)點(diǎn)數(shù)為N，已知壓力節(jié)點(diǎn)數(shù)為m。

基于管網(wǎng)水力計算基本方程的理論，應(yīng)用圖論理論和線性代數(shù)的思想建立管線與節(jié)點(diǎn)間的關(guān)系式。熱力管網(wǎng)穩(wěn)態(tài)運(yùn)行時，均滿足節(jié)點(diǎn)流量方程、管線能量方程、壓降方程以及環(huán)方程[8]。

(1)節(jié)點(diǎn)流量方程。根據(jù)質(zhì)量守恒定律可得:

將N個節(jié)點(diǎn)流量連續(xù)性方程相加可得到節(jié)點(diǎn)流量守恒方程式:

(2)管線能量方程。根據(jù)能量守恒定律可得:

(3)壓降方程?？紤]管線上加壓泵的設(shè)置和水平位置的變化等情況，管線壓降方程的表達(dá)式為

若不考慮上述因素影響，由假設(shè)(3)可知，對任意管線的管線壓降公式可表示為

式中:Si為管線的阻抗系數(shù)，取決于管材、管長和管徑。

(4)環(huán)方程。環(huán)方程的物理意義是指每個閉合回路中各管道的壓力損失總和為零，也稱回路方程;表達(dá)式為

由于熱力管網(wǎng)供、回水管一般對稱布置，且同一位置的供、回水管的管徑、流量和阻力均相等。為了便于編程和計算，將熱力管網(wǎng)的供、回水管分離，建模時只考慮供水管或回水管的水力特性，對于已知的熱源輸出流量和熱力站設(shè)計流量均當(dāng)作節(jié)點(diǎn)流量處理。

根據(jù)熱力管網(wǎng)水力計算的基本方程，對于一個管線數(shù)為M、節(jié)點(diǎn)數(shù)為N的熱力管網(wǎng)，可建立熱力管網(wǎng)的通用數(shù)學(xué)模型如式(7):

式(1―7)中:A為N×M階關(guān)聯(lián)矩陣，表示管線和節(jié)點(diǎn)間的拓?fù)潢P(guān)系;AT為管網(wǎng)關(guān)聯(lián)矩陣A的轉(zhuǎn)置矩陣;A0為管網(wǎng)A的降階關(guān)聯(lián)矩陣，為(N－m)×M 階矩陣;q 為管線流量向量，q=(q1，q2，…qM)T;Q 為 節(jié) 點(diǎn) 流 量 向 量，Q = (Q1，Q2，…QN-m)T;H 為 節(jié) 點(diǎn) 壓 力 向 量，H =(H1，H2，…HN)T;ΔH為管段壓力損失或壓降向量，ΔH=(ΔH1，ΔH2，…ΔHM)T;S 為 N 階對角阻力系數(shù)矩陣，S=diag(S1，S2，…SM);Z為管線上兩節(jié)點(diǎn)的勢能差向量，Z=(Z1，Z2，… ZM)T;Hp為管線中加壓泵的揚(yáng)程向量，Hp=(Hp1，Hp2，… HpM)T;B 為(M－R)×M階基本回路矩陣，R為關(guān)聯(lián)矩陣的秩。

1.2 熱力管網(wǎng)水力模型的求解

當(dāng)已知管徑時，熱力管網(wǎng)模型的方程組含有2M個未知數(shù)，可列出的方程數(shù)為

因此，方程個數(shù)與未知數(shù)的數(shù)目相等，該模型可以求解。但是，由于非線性方程組直接求解困難，本文中采用環(huán)方程法求解熱力管網(wǎng)模型。

環(huán)方程法的基本思想是在滿足節(jié)點(diǎn)流量方程的條件下，構(gòu)造各環(huán)的能量方程組，通過求解各環(huán)校正流量的方法間接求出各管線流量，即Hardy Cross法，采用的環(huán)方程法需先擬定各管段流量初值[9]。

一般熱力管網(wǎng)最初的設(shè)計流量難以滿足環(huán)能量方程，因而第i環(huán)的環(huán)能量方程可表示為

式中:ΔHi為第i環(huán)的壓降閉合差。

為了滿足各環(huán)壓降閉合差的允許精度和節(jié)點(diǎn)流量方程的要求，引入第i環(huán)校正流量δqi，此時第i環(huán)的環(huán)能量方程為

式中:δqi的正負(fù)號與 bij一致;δqk為第 i環(huán)第 j管段鄰環(huán)校正流量，正負(fù)號與bij相反。

采用麥克勞林級數(shù)展開(10)式，忽略 δqi，δqk二階多項式，簡化后表達(dá)式為

由式(11)構(gòu)造環(huán)能量方程組，采用高斯－賽德爾法求解關(guān)于δHi的線性方程組，其方程組表達(dá)式為

[B·Λ·BT]·δq=1/2·B·Λ·q(12)式中:Λ為管線特性的對角矩陣，Λ=diag(S1|q1|，S2|q2|，…SM|qM|);δq 為管網(wǎng)環(huán)的校正流量，δq=(δq1，δq2，… δqU)T。其求解步驟如圖2 所示。

2 報警閥 ε0和 εF定義

為實現(xiàn)實時在線故障診斷、定位熱力管網(wǎng)泄漏點(diǎn)、減少管網(wǎng)診斷系統(tǒng)誤報率[10]，需設(shè)置管線報警閥ε0和泄漏點(diǎn)報警閥εF。

2.1 管線報警閥ε0定義

熱力管網(wǎng)正常運(yùn)行時，假設(shè)管線i兩端的節(jié)點(diǎn)壓力為 Hi，Hi+1，則

假設(shè)泄漏發(fā)生在管線i上，則其兩端節(jié)點(diǎn)壓力將變?yōu)镠'i，H'i+1，泄漏量為QF，若將泄漏點(diǎn)視為插入輸出節(jié)點(diǎn)，此時管線i仍將滿足壓降方程，即

泄漏前后，管線i兩端的節(jié)點(diǎn)壓力差值為

基于建模假設(shè)(3)，沿程損失系數(shù)λ僅與管壁相對粗糙度Δ/d有關(guān)。對于已建成的熱力管網(wǎng)，管壁絕對粗糙度Δ和管徑d為常數(shù)，可令管線阻抗 S=0.8106λρ0L/d5，忽略 QF二階函數(shù)，可得:

則管線報警閥ε0可定義為

式中:ξ為衡量熱力管網(wǎng)管線泄漏報警率的特征因子，稱為靈敏度因子。ξ取值減小，報警閥ε0隨之減小，這說明熱力管網(wǎng)泄漏診斷的靈敏度升高，可導(dǎo)致熱力管網(wǎng)故障診斷的誤報率增加，反之亦然。因此，靈敏度因子ξ的選定是至關(guān)重要的。在下述算例中，靈敏度因子ξ=6，取自熱力管網(wǎng)故障診斷系統(tǒng)的調(diào)試經(jīng)驗值。

2.2 泄漏點(diǎn)報警閥εF的定義

城市熱力管網(wǎng)一般采用中壓供暖，管線壓力波動范圍為3% ～5%。若管線出現(xiàn)泄漏點(diǎn)，則壓力變化將超過該范圍。故泄漏點(diǎn)報警閥εF可定義為管線i正常工作時虛擬泄漏點(diǎn)處壓力值的4%，即

在定位泄漏點(diǎn)的過程中，采用黃金分割法插入泄漏點(diǎn)，所插入泄漏點(diǎn)的定位邏輯圖如圖3所示。

圖3 泄漏管線泄漏點(diǎn)的定位邏輯圖Fig．3 Logic diagram of leakage of pipeline leak location

3 算例及分析

3.1 管網(wǎng)正常運(yùn)行狀態(tài)的水力計算結(jié)果

熱力管網(wǎng)正常運(yùn)行時，所有節(jié)點(diǎn)的流量代數(shù)和為零，節(jié)點(diǎn)流量和壓力數(shù)值如表1所示，管線計算流量和壓降數(shù)值如表2所示。

表1 管網(wǎng)正常運(yùn)行狀態(tài)時節(jié)點(diǎn)參數(shù)值Tab.1 Node parameter value of network normal operation

表2 管網(wǎng)正常運(yùn)行狀態(tài)時管線參數(shù)值Tab.2 Pipeline parameter value of network normal operation

3.2 管網(wǎng)泄漏故障狀態(tài)的水力計算結(jié)果及分析

假設(shè)熱力管網(wǎng)管線7上坐標(biāo)為(294.85，904.44)點(diǎn)發(fā)生泄漏，泄漏量為5 m3/h，此時熱力管網(wǎng)所有節(jié)點(diǎn)的流量代數(shù)和不為零，管網(wǎng)節(jié)點(diǎn)流量和壓力數(shù)值見表3，管線計算流量和壓降數(shù)值見表4，熱力管網(wǎng)各管線報警閥ε0的取值與管線壓降實際值和計算值的差值ε的比較結(jié)果見表5。由上述計算結(jié)果可知:

(1)管線7，13，18的壓降實際值和計算值的差值ε均大于報警閥。由圖1可以看出:管線13，18均位于管線7的上游，而管線泄漏搜索從最末端開始，且管線7單位長度壓力損失變化率遠(yuǎn)大于管線13，18的壓力損失變化率;依據(jù)判定標(biāo)準(zhǔn)，確定管線7上存在泄漏點(diǎn);管線13，18出現(xiàn)的壓力波動和供水量變化主要是由于下游泄漏點(diǎn)泄漏對上游管路的影響造成的，也不排除與計算時造成的舍入誤差有關(guān)。

表3 管網(wǎng)泄漏故障狀態(tài)時節(jié)點(diǎn)參數(shù)值Tab.3 Node parameter value of network leakage fault

表4 管網(wǎng)泄漏故障狀態(tài)時管線參數(shù)值Tab.4 Pipeline parameter value of network leakage fault

表5 ε0與ε比較結(jié)果Table．5 Comparison results of ε0and ε

(2)管線7的上游主干管線13出現(xiàn)壓降實際值大于計算值的管線，且差值ε大于報警閥ε0，這是由于管網(wǎng)發(fā)生泄漏時，管線13較長，其累積壓降大，因此出現(xiàn)報警;同時由于管線7的泄漏僅使上游支管流量略有增加，使得管線7上游各支管壓降實際值近似等于計算值。

(3)管線7的下游管線的壓降實際值近似等于計算值，這是供水量的補(bǔ)償導(dǎo)致下游管線泄漏前后壓降變化不大的結(jié)果。

作者所開發(fā)的熱力管網(wǎng)故障診斷系統(tǒng)按圖3所示邏輯可定位泄漏管線泄漏點(diǎn)，即xF=73.64，其所在的管線7長度L7=180 m。將其轉(zhuǎn)換成地理圖集上的坐標(biāo)為(x1，y1)=(292.97，910.83)，假設(shè)的泄漏點(diǎn)坐標(biāo)為(x2，y2)=(294.85，904.44)，則該系統(tǒng)定位的誤差可按式(18)計算:

代入數(shù)據(jù)得δ=4.82%

可見系統(tǒng)泄漏點(diǎn)的定位誤差在工程要求范圍內(nèi)。因此，本文所提出的枝狀熱力管網(wǎng)單點(diǎn)泄漏的故障診斷方法是合理的，具有一定應(yīng)用價值[11，12]。

4 結(jié)論

(1)為某市熱力管網(wǎng)建立了熱力管網(wǎng)水力計算模型，以此為基礎(chǔ)開發(fā)的熱力管網(wǎng)故障診斷系統(tǒng)可確定泄漏點(diǎn)的大致位置，其定位誤差在工程要求的范圍內(nèi)。因此，文中所提出的插入虛擬節(jié)點(diǎn)法及報警閥ε0，εF的值域的設(shè)定是合理的，具有一定的應(yīng)用價值。

(2)泄漏管線對上游管線誤報警有一定的影響，但也不排除與計算時造成的舍入誤差有關(guān)，此點(diǎn)尚需進(jìn)一步研究。

(3)本文的研究成果雖然僅適用于單點(diǎn)泄漏的枝狀熱力管網(wǎng)故障診斷，但對多泄漏點(diǎn)的熱力管網(wǎng)可找出離末端最近的泄漏管線，具有一定的參考價值。

[1]Alan T．Murray，Tong Daoqin．GIS and spatial analysis in the media[J]．Applied geography，2009，(29):250－259．

[2]Lopes dos Santos P，Azevedo-Perdicou，Ramos J A，et al．An LPV Modeling and Identification Approach to Leakage Detection in High Pressure Natural Gas Transportation Networks[J]．IEEE Control Systems Society，2011，19(1):77－79．

[3]楊杰，王桂增．輸氣管道泄漏診斷技術(shù)綜述[J]．化工自動化及儀表，2004，31(11):1－5．

[4]王瑛周，李春剛．SCADA技術(shù)在城市供水系統(tǒng)中的應(yīng)用[J]．經(jīng)濟(jì)技術(shù)協(xié)作信息，2010，11(1030):198－199．

[5]陶建科，胡永清．基于GIS和SCADA技術(shù)的供水管網(wǎng)管理信息系統(tǒng)[J]．供水技術(shù)，2008，2(2):29－32．

[6]Zhang Yu，Liu Jian，Zeng Zhoumo．A combined kalman filter-Discrete wavelet transform method for leakage detection of crude oil pipelines[C]．Shanghai:ICEMI，2009:1086 －1090．

[7]卜月華．圖論及其應(yīng)用[M]．南京:東南大學(xué)出版社，2002．

[8]秦芳芳．供熱管網(wǎng)水力計算模型研究[D]．保定:華北電力大學(xué)，2008．

[9]劉賢榮，倪鵬．外部多水源給水管網(wǎng)水力計算程序的應(yīng)用[J]．華僑大學(xué)學(xué)報，2005，26(2):15－19．

[10]畢鋒東，李樹杰，陳久會，等．人工智能型管道泄漏監(jiān)測系統(tǒng)在克—烏成品油管道的應(yīng)用[J]．管道技術(shù)與設(shè)備，2005，22(6):22－24．