韓宇波，常 暢，宋 琪，洪文鵬

(1.珠海橫琴能源發(fā)展有限公司，廣東 珠海 519015；2.東北電力大學(xué)能源與動力工程學(xué)院，吉林 吉林 132012)

目前大多比較先進(jìn)成熟的泄漏檢測定位方法和技術(shù)均是針對單一的長直管道，例如長輸油、輸天然氣管道[1-3].然而，對于拓?fù)浣Y(jié)構(gòu)復(fù)雜的多分支型管網(wǎng)來說，例如供熱管網(wǎng)、供水管網(wǎng)，這些方法將不能直接被利用.復(fù)雜供水管網(wǎng)廣泛存在破損、泄漏等問題.給社會生產(chǎn)和人民生活造成巨大經(jīng)濟(jì)損失以及諸多不便，因此，為了及時的發(fā)現(xiàn)，定位，并修補(bǔ)泄漏，消除水泄漏造成的潛在危險，本文針對供水管網(wǎng)泄漏檢測定位，提出了基于伴隨方程的反向?qū)ぴ捶椒?

1994年，靳世久[4]提出基于負(fù)壓波的管道泄漏檢測定位方法.1996年，楊開林、郭宗周[5]提出了基于恒定流動模擬的靜態(tài)泄漏檢測法和基于水力瞬態(tài)模擬的瞬態(tài)泄漏檢測法.2001年，Verde和Mpesha[6-7]分別提出了最小非線性觀測器和頻率響應(yīng)法來確定泄漏點(diǎn)的位置.2002年，Costa[8]運(yùn)用參數(shù)估計法來對管網(wǎng)進(jìn)行泄漏檢測和泄漏點(diǎn)定位.2003年，F(xiàn)errante[9]提出了基于小波的分析方法，根據(jù)上、下游傳感器的時間差進(jìn)行泄漏點(diǎn)定位.2007年，劉暢等[10]利用多元統(tǒng)計學(xué)上的聚類分析方法，對管網(wǎng)故障模擬的數(shù)據(jù)進(jìn)行聚類分析，并利用多元統(tǒng)計學(xué)上的判別分析方法，對泄漏情況進(jìn)行故障定位.2010年，楊紅英等[11]研究了基于Duffing振子的管道泄漏檢測方法.2014年，呂玉坤等[12]對于單點(diǎn)泄漏的枝狀供熱管網(wǎng)，引入插入虛擬節(jié)點(diǎn)的方法，來判斷管網(wǎng)是否泄漏與泄漏可能發(fā)生的位置.2015年，Wachla等[13]采用人工神經(jīng)模糊模型對供水管道泄漏進(jìn)行定位.2015年，燕山大學(xué)黃滿義等[14]針對單一參數(shù)研究輸油管道泄漏，引入混沌特性分析及遞歸圖方法對管道泄漏檢測與定位方法進(jìn)行分析.2016年，劉煒，劉宏昭[15]從時間序列角度出發(fā)，提出一種基于模糊集理論的輸油管道泄漏檢測定位方法.

Enting[16]將反問題用于大氣污染物擴(kuò)散研究.Bagtzoglou[17]將反演方法用于地下水污染定位研究.郭少東等[18]采用計算流體力學(xué)的方法求解伴隨方程和MCMC方法的室內(nèi)污染源反演結(jié)果的研究.Zhang[19]將源項(xiàng)反演方法用在室內(nèi)污染源定位研究.

Liu和Zhai[20-22]使用伴隨方法實(shí)現(xiàn)了室內(nèi)污染源的快速辨識，但他們并沒有將其作為最優(yōu)化方法使用，而是利用伴隨方法的特性，將在流場固定的基礎(chǔ)上實(shí)現(xiàn)了時間和對流項(xiàng)的取反，以反向方法的方式實(shí)現(xiàn)了污染源的辨識.

總結(jié)以上國內(nèi)外研究現(xiàn)狀發(fā)現(xiàn)，現(xiàn)有的管網(wǎng)漏損檢測定位方法并不能滿足日常管網(wǎng)漏損檢測的需求，在精度或者經(jīng)濟(jì)性上均存在一定的缺點(diǎn).目前，基于流場求解的反向?qū)ぴ醇夹g(shù)能夠?qū)崿F(xiàn)準(zhǔn)確地檢測和定位，在診斷空氣污染和水污染源頭的方向上取得了較好的效果，因此，本文將基于伴隨方程的反向?qū)ぴ醇夹g(shù)應(yīng)用于供水管網(wǎng)中，以更快的實(shí)現(xiàn)復(fù)雜管網(wǎng)的泄漏檢測和定位.

1 數(shù)學(xué)模型的建立

1.1 單相流瞬態(tài)伴隨方程的推導(dǎo)

建立管網(wǎng)流場瞬態(tài)方程，找出與流場瞬態(tài)壓力耦合的參數(shù)，保留這些參數(shù)，對流場瞬態(tài)方程進(jìn)行簡化.

不考慮摩擦、阻尼等因素的影響，一維半無限區(qū)間下壓力瞬態(tài)流動方程為

(1)

初始條件為

當(dāng)t=0時，p(x，t)=0，v(x，t)=0；

邊界條件為

當(dāng)x=0時，p(x，t)=p0(t)，v(x，t)=0；當(dāng)x→+∞時，p(x，t)=0；

公式中：p為壓力，Pa；v為液體流速，m3/s；a為壓力波速，m/s；ρ為液體密度，kg/m3.

(1)以管網(wǎng)泄漏對流場形成的壓力波影響為依據(jù)，選取適當(dāng)敏感度函數(shù).所選取的壓力波敏感函數(shù)

h(α，p)=p(x，t)δ(x-x*)δ(t-t*)，

(2)

公式中：p(x，t)為壓力波傳播函數(shù)；α為狄克拉函數(shù)，α與α為泄漏或堵塞故障發(fā)生的位置與時間；α為敏感度參數(shù)(過程變量)；p為壓力.

(2)將敏感度函數(shù)對壓強(qiáng)求導(dǎo)，并引入伴隨算子，推導(dǎo)管網(wǎng)流場瞬態(tài)伴隨方程.利用高斯散度定理處理，設(shè)置邊界條件，得到流場瞬態(tài)伴隨方程.

量化系統(tǒng)狀態(tài)量的度量目標(biāo)函數(shù)

(3)

公式中：h(α，p)為系統(tǒng)狀態(tài)函數(shù)；α為系統(tǒng)矢量參數(shù)(例如[v，D，a]；p為壓力，積分區(qū)間是給定的空間-時間區(qū)間，系統(tǒng)關(guān)于矢量參數(shù)αk的臨界靈敏度可以通過對參數(shù)αk求微分得到

(4)

(5)

初始條件為

當(dāng)t=0時，Φ(x，t)=0；當(dāng)t=0時，λ(x，t)=0；

邊界條件為

當(dāng)x→+∞時，φ(x，t)=0；

當(dāng)x=0時，λ(x，t)=0，φ(x，t)=0；當(dāng)x→+∞時，λ(x，t)=0.

(6)

初始條件為

φ*(x，T)=0，λ*(x，T)=0；

邊界條件為

φ*(0，t)=p0(t)，λ*(0，t)=0.

這樣伴隨方程(6)中未定義函數(shù)僅剩余h(α，p)項(xiàng)，它的形式還需要根據(jù)所要研究的具體工況條件確定，將方程(2)中選取的適當(dāng)函數(shù)兩邊取偏導(dǎo)可得

(7)

將方程(7)帶入到伴隨方程(6)中，經(jīng)過整理得

(8)

初始條件為

φ*(x，T)=0，λ*(x，T)=0；

邊界條件為

φ*(0，t)=p0(t)，λ*(0，t)=0.

考慮流動過程中摩擦阻力，最終得出管網(wǎng)流場瞬態(tài)伴隨方程

(9)

邊界條件為

φ*(0，t)=p0(t)，φ*(L，t)=0，λ*(0，t)=0，λ*(L，t)=0；

初始條件為

φ*(x，T)=0，λ*(x，T)=0.

公式中：t為時間變量；x為距離變量；L為距離常數(shù)；T為時間常數(shù).

(3)將伴隨方程進(jìn)行空間與時間的反向處理，應(yīng)用拉普拉斯變換和傅里葉變換求解析解，得到反向?qū)ぴ茨Ｐ?

做反向處理，τ=td-t，x為-x，可得反向伴隨方程

(10)

公式中：τ=td-t，td為檢測的時間；xd為檢測的位置.

1.2 模型驗(yàn)證

傳統(tǒng)的正向求解方法是通過對每個可能的泄漏位置進(jìn)行正向求解，判斷泄漏點(diǎn)的位置，但是這種方法因?yàn)樾枰啻吻蠼鈴亩兊牡托?因此，本文利用伴隨理論和數(shù)學(xué)方法結(jié)合流體力學(xué)中的單相流體控制方程推導(dǎo)伴隨方程，建立逆向模型，如圖1所示.

圖1 泄漏位置檢測原理圖

從檢測點(diǎn)開始傳播的壓力變化曲線，如圖2所示.從圖2中可以看到明顯的突升變化過程，該曲線是通過在假定的壓力測定點(diǎn)出檢測到負(fù)壓波信號，并設(shè)置為xd=0，即為原點(diǎn)，得到反向傳播時間，τ=5 000 m/1 300 m/s，且泄漏點(diǎn)出的壓力變化為突升25 kPa.

圖2 從檢測點(diǎn)開始傳播的壓力變化曲線圖3 檢測點(diǎn)開始傳播的壓力曲線

假設(shè)在xd=0監(jiān)測點(diǎn)處開始檢測壓力波，如圖3所示，模擬了泄漏時間τ=5.36 s時，壓力波反向傳播變化曲線，壓力值從25 kPa突變到50 kPa，泄漏突變的位置在x=6 968 m處.

綜上所述，基于壓力輸運(yùn)控制方程，應(yīng)用伴隨理論及一些數(shù)學(xué)方法推導(dǎo)出相應(yīng)的管網(wǎng)瞬態(tài)伴隨方程，通過MATLAB模擬了檢測點(diǎn)處逆向壓力波變化規(guī)律，明確了由壓力輸運(yùn)方程推導(dǎo)出來反向伴隨方程可以應(yīng)用于管道泄漏檢測與定位.

2 實(shí)驗(yàn)結(jié)構(gòu)及數(shù)值模擬計算

2.1 實(shí)驗(yàn)結(jié)構(gòu)

以實(shí)驗(yàn)為基準(zhǔn)驗(yàn)證使用Open FOAM進(jìn)行數(shù)值模擬的結(jié)果，實(shí)驗(yàn)中對管道泄漏前后以及泄漏瞬間的壓力進(jìn)行了測量.實(shí)驗(yàn)系統(tǒng)圖如圖4所示，主要由一個實(shí)驗(yàn)段，五個壓力變送器，兩個電磁流量計，兩個過濾器，一個水箱，一個水泵，一個數(shù)據(jù)采集儀和一個電腦組成，整體分為復(fù)雜管網(wǎng)微縮實(shí)驗(yàn)系統(tǒng)、信號傳輸系統(tǒng)和數(shù)據(jù)處理系統(tǒng)三個基礎(chǔ)部分.

圖4 實(shí)驗(yàn)系統(tǒng)圖

由于實(shí)際供水管網(wǎng)雷諾數(shù)大致范圍在5×104到4×105之間，且管網(wǎng)穩(wěn)定運(yùn)行時管內(nèi)流態(tài)為阻力平方區(qū)，結(jié)合模擬供水管網(wǎng)、運(yùn)行參數(shù)及實(shí)驗(yàn)室空間，根據(jù)相似理論進(jìn)行了管網(wǎng)模型設(shè)計.選擇水為管網(wǎng)運(yùn)行介質(zhì)，為了滿足可視性，使用強(qiáng)度大、安全性高、安裝簡單等要求，選擇有機(jī)玻璃為管網(wǎng)主體管道，U-PVC為閥門和四通材料.

管網(wǎng)主體為4×4復(fù)雜環(huán)狀管網(wǎng)，如圖5所示.管材為內(nèi)徑12 mm外徑20 mm的有機(jī)玻璃，泄漏內(nèi)徑4 mm，外徑6 mm.管網(wǎng)各環(huán)管段尺寸相同，長均為0.85 m，寬均為0.38 m，管段中共有25個節(jié)點(diǎn)，其中5個節(jié)點(diǎn)為壓力測點(diǎn).共設(shè)置20個泄漏點(diǎn)，其中10個位于管段節(jié)點(diǎn)，另外10個為管段中點(diǎn)，分別對節(jié)點(diǎn)和管段中泄漏進(jìn)行研究.

圖5 實(shí)驗(yàn)臺壓力測點(diǎn)和泄漏點(diǎn)的分布

2.2 數(shù)值模擬計算

2.2.1 網(wǎng)格的劃分

根據(jù)實(shí)驗(yàn)段尺寸，建立模擬模型，忽略壁厚，管道直徑12 mm，泄漏孔直徑4 mm，整體采用三維六面體結(jié)構(gòu)化網(wǎng)格，管道交叉處的與泄漏口處分別采用Y-Block和O-block，實(shí)驗(yàn)段為4×4管網(wǎng)，網(wǎng)格數(shù)達(dá)9 061 188，為了便于計算，將其簡化為2×2管網(wǎng)，網(wǎng)格數(shù)減少到4 420 092，計算效率大大提高.

圖6 4*4模型網(wǎng)格劃分圖7 2*2模型網(wǎng)格劃分

本文選擇三維六面體結(jié)構(gòu)化網(wǎng)格，如表1所示.表1做了網(wǎng)格無關(guān)性驗(yàn)證，本次模擬采用的網(wǎng)格數(shù)為4 420 092.

表1 網(wǎng)格無關(guān)性驗(yàn)證

由表1可知，以泄漏點(diǎn)1為例，當(dāng)網(wǎng)格數(shù)達(dá)到400萬，監(jiān)測點(diǎn)壓力的壓降幾乎不變，為了提高管網(wǎng)泄漏模擬計算的精度和效率，對管道相貫處，相貫處的泄漏點(diǎn)，以及管道進(jìn)出口等進(jìn)行了網(wǎng)格加密處理.

2.2.2 復(fù)雜管網(wǎng)泄漏檢測及相關(guān)性分析

同一泄漏源產(chǎn)生的壓力波，其傳播到不同傳感器所需的時間不同，因此不同傳感器上所接收到的信號記錄也不會相同，但是當(dāng)把兩個傳感器記錄其中之一給予一定的時移后，會發(fā)現(xiàn)此時兩個信號序列之間具有很大的相似性.根據(jù)壓力波傳到兩個監(jiān)測點(diǎn)的時間差，壓力波波速即可確定泄漏點(diǎn)的位置.

為提高時延估計的精度，采用了希爾伯特變換和二次相關(guān)法相結(jié)合的新時延估計方法，其基本思路在于結(jié)合兩種分析方法，將其與相關(guān)函數(shù)的絕對值相減，相關(guān)峰值點(diǎn)被保留，其附近的點(diǎn)的相關(guān)值相應(yīng)的減小，主峰值點(diǎn)被銳化，精度大大提高，在如上所述的二次相關(guān)法的基礎(chǔ)上，加入了希爾伯特變換，并結(jié)合了二次相關(guān)方法，原理如圖8所示.

圖8 新時延估計法原理圖

銳化了主峰值點(diǎn)后，利用互相關(guān)系數(shù)法判斷信號相似程度，其原理是兩列信號的相關(guān)系數(shù)在泄漏前通常在零值附近，出現(xiàn)泄漏后，信號的相關(guān)系數(shù)在零附近有所偏移，信號峰值點(diǎn)所對應(yīng)的時間就是相關(guān)系數(shù)最大時所對應(yīng)的時間點(diǎn).假定兩列不同隨機(jī)信號s1(t)、s2(t)，互相關(guān)函數(shù)為

(11)

公式中：RS1S2(τ)可以計算出s1(t)和s2(t)信號的相似程度，它是信號的相關(guān)函數(shù)，其中，τ為時移.

以泄漏孔1、泄漏孔2、泄漏孔3為例，壓力波的傳播速度1 483.2 m/s，設(shè)置時間步長為5e-7s，采用PISO算法，自定義求解器，管網(wǎng)內(nèi)流體流動穩(wěn)定到2 s，出現(xiàn)泄漏.

簡化4×4管網(wǎng)之后的管網(wǎng)系統(tǒng)圖，如圖9所示.設(shè)置兩個壓力測點(diǎn)，測點(diǎn)1(左上方測點(diǎn))和測點(diǎn)2(右下方測點(diǎn))，監(jiān)測泄漏前后的壓力波信號.模擬泄漏前和整個過程測點(diǎn)1和測點(diǎn)2通過新時延估計法和互相關(guān)分析法處理后得到的信號相關(guān)性圖，如圖10～圖12所示.相關(guān)系數(shù)的最大值在泄漏前通常在零值附近，出現(xiàn)泄漏后，信號相關(guān)系數(shù)的最大值在零附近有所偏移，通過相關(guān)性對比圖可知，管道泄漏的出現(xiàn)，最大相關(guān)系數(shù)的值下降，我們通過對比管道泄漏前后的閾值，判斷管道是否出現(xiàn)了泄漏.

圖9 管網(wǎng)布置示意圖

圖10 不同測點(diǎn)，1點(diǎn)泄漏前后壓力及相關(guān)性對比圖

圖11 不同測點(diǎn)，2點(diǎn)泄漏前后壓力及相關(guān)性對比圖

圖12 不同測點(diǎn)，3點(diǎn)泄漏前后壓力及相關(guān)性對比圖

從上圖10至圖12中可知，管網(wǎng)出現(xiàn)泄漏后，壓力瞬間增大，與實(shí)際實(shí)驗(yàn)數(shù)據(jù)相反，這是由于模擬過程添加了逆向伴隨模型，大大加快了計算效率.管網(wǎng)流動穩(wěn)定后，2 s時出現(xiàn)泄漏，壓力波迅速向管道上下游進(jìn)行傳播，在泄漏后的0.2 s內(nèi)壓力變化明顯，出現(xiàn)上升的趨勢，管道內(nèi)部的沿程壓力也升高，略微波動后穩(wěn)定于略微高于管道正常壓力的值.管道入口我們采用速度進(jìn)口，入口速度一定即流量一定，泄漏后短時間內(nèi)，泄漏點(diǎn)上游的流速下降后又恢復(fù)到了泄漏前的穩(wěn)定值，左上方測點(diǎn)距離泄漏孔較近，右下方測點(diǎn)較遠(yuǎn)，其接收到壓力波信號時間大于左上方測點(diǎn)，壓力值大于左上方測點(diǎn)，總體趨勢相同，壓力變化范圍在10 kPa～28 kPa之間，泄漏后相關(guān)系數(shù)減小，最大相關(guān)系數(shù)對應(yīng)的時間有延遲，符合相關(guān)性驗(yàn)證.

管道內(nèi)流體流動2 s穩(wěn)定后，1，2兩點(diǎn)同時發(fā)生泄漏的模擬圖，相對于一點(diǎn)泄漏，兩點(diǎn)泄漏的峰值差更小，壓力變化范圍也較小，但總體趨勢一致，符合逆向模型得出的結(jié)果，可以利用上述方法對管網(wǎng)實(shí)現(xiàn)在線監(jiān)測，如圖13所示.

圖13 不同測點(diǎn)，1，2兩點(diǎn)同時泄漏前后壓力對比圖

綜上所述，管道突然出現(xiàn)泄漏后，相比于泄漏前，管道首末兩端的壓力差更大，壓降作用也更明顯.加入了反向模型的管網(wǎng)，泄漏后管網(wǎng)各個監(jiān)測點(diǎn)的壓力迅速上升后略有降低，之后維持在一個較高的值，比較不同泄漏位置的壓力值，泄漏位置靠近管道進(jìn)口，下游壓力降低幅度越高，越靠近管道出口，上游壓力下降幅度越小，根據(jù)管道沿程各監(jiān)測點(diǎn)收集到的壓力值梯度，可知泄漏情況的發(fā)生.

2.2.3 復(fù)雜管網(wǎng)泄漏位置的確定

管網(wǎng)泄漏檢測及定位流程圖，如圖14所示.當(dāng)管網(wǎng)出現(xiàn)泄漏時，我們首先通過上述方法進(jìn)行泄漏檢測，確定可能發(fā)生泄漏的位置，判斷各部分泄漏量的多少，然后利用全局搜索方式，最終確定泄漏點(diǎn).

結(jié)合模擬結(jié)果，我們分析了兩列信號s1(t)和s2(t)的相關(guān)性(以1點(diǎn)泄漏采集為例)，計算可知壓力波傳到兩個監(jiān)測點(diǎn)的時間差為△t=1.15×10-3s，根據(jù)已知的壓力波波速信息就能計算出泄漏孔距離兩個測點(diǎn)的距離，兩者距離差約五倍，由于管網(wǎng)線是同程的，可知可能出現(xiàn)泄漏的位置如下圖所示的A-D，根據(jù)圖15所示，測點(diǎn)1測得的壓力值較大，且有所延遲，說明泄漏點(diǎn)在測點(diǎn)2附近，泄漏位置則是泄漏點(diǎn)A或B.

圖15 泄漏可能出現(xiàn)位置示意圖

因此，為了識別泄漏位置，提出了一種搜索算法.為了排出其他泄漏點(diǎn)，我們采取了新的方案，該方案包括兩個特點(diǎn)：它包括環(huán)形網(wǎng)絡(luò)；測量點(diǎn)的排列可以符合試驗(yàn)的要求.這個管道泄漏產(chǎn)生的負(fù)壓波可以到達(dá)壓力測量點(diǎn)通過不同的方式，我們可以利用這個特性實(shí)驗(yàn)誤差分析.具體地：

步驟1：搜索最接近泄漏位置的節(jié)點(diǎn)

(12)

公式中：tj、tk為在tj和tk時刻，在節(jié)點(diǎn)j和k檢測到泄漏瞬變；Si為較小的殘值表示泄漏發(fā)生在i節(jié)點(diǎn)的概率較高；τik為讓我們把節(jié)點(diǎn)i到k的最短旅行時間表示出來.

假如泄漏在節(jié)點(diǎn)i發(fā)生，當(dāng)i=1……N(N=管網(wǎng)中的節(jié)點(diǎn)數(shù)目)，從而

(tj-tk)-(τij-τik)=0.

(13)

然而，由于時間、測量和其他誤差，這個等式的左邊永遠(yuǎn)不會為零.

步驟2：沿著連接到最近節(jié)點(diǎn)的管道段搜索泄漏位置.

在這一步驟中，沿著步驟1確定的與節(jié)點(diǎn)ni相連的管道部分(即沿著圖模型的邊緣)放置一組新的虛擬節(jié)點(diǎn).

第一步是對網(wǎng)絡(luò)中所有節(jié)點(diǎn)進(jìn)行粗略的全局搜索；

第二步在最近的節(jié)點(diǎn)估計附近執(zhí)行本地搜索，以確定沿著管道段的最可能泄漏位置.

此外，管道泄漏產(chǎn)生的壓力波傳播將通過不同的路徑到達(dá)測量點(diǎn)，我們可以以A點(diǎn)泄漏時，M1和M2點(diǎn)為例.

圖16 A點(diǎn)泄漏時，傳遞到測點(diǎn)M2的路徑

圖17 A點(diǎn)泄漏時，傳遞到測點(diǎn)M1的路徑

如圖16和17所示，假設(shè)A點(diǎn)泄漏，則最符合本次模擬的路線是圖16中的(a)和圖17中的(b).

表2為根據(jù)圖14所示的方法計算出的的每個節(jié)點(diǎn)的Ei值，由表中可以看出節(jié)點(diǎn)3的Ei值最小，節(jié)點(diǎn)4的Ei值次之，故離泄漏位置最近的節(jié)點(diǎn)為3，并且泄漏位置介于節(jié)點(diǎn)3和4之間，泄漏點(diǎn)為點(diǎn)A.根據(jù)△t=1.15×10-3s，已知壓力波的波速，我們可以求出泄漏點(diǎn)最終確定泄漏位置414 mm，實(shí)際泄漏點(diǎn)在431 mm處(泄漏孔1)，在誤差允許的范圍內(nèi)，所以該方法可行.

表2 不同節(jié)點(diǎn)Ei計算結(jié)果

同時，模擬了不同泄漏孔大小和不同進(jìn)口壓力下的5個泄漏點(diǎn)的定位情況，分別進(jìn)行了30次和20次模擬計算，為了便于分析，取以測點(diǎn)2為起點(diǎn)，傳播路徑最短的一組數(shù)據(jù)，現(xiàn)以泄漏孔3和泄漏孔5為例，定位情況如表3、表4所示，誤差均在允許的范圍內(nèi)，說明該方法可行.

表3 以泄漏孔3為例的不同泄漏孔的定位情況

表4 以泄漏孔5為例的不同進(jìn)口壓力的定位情況

2.2.4 不同方法及算法對比

2.2.4.1 不同方法對比

由表5可知，近些年現(xiàn)有的檢測方法如ITA(Inverse Transient Analysis)，Model Falsification & EPANET，PSO-SVM(Particle Swarm Optimization-Support Vector Machine)等方法的定位能力較高，均在90%以上.本文提到的于反向?qū)ぴ捶ǘㄎ坏臏?zhǔn)確度在92%～96%之間，在復(fù)雜管網(wǎng)的泄漏定位中有一定的應(yīng)用價值.

表5 不同泄漏檢測及定位方法對比

2.2.4.2 算法對比

如圖18所示，利用表1中經(jīng)過網(wǎng)格無關(guān)性驗(yàn)證的網(wǎng)格數(shù)為442 009的管網(wǎng)模型，在相同的邊界條件和參數(shù)的設(shè)置下，對比了在50 000、100 000、150 000、200 000計算步數(shù)下，SIMPLE、PISO及加入伴隨方程的PISO算法所需的迭代時間，隨著計算步數(shù)的增多，三種方式的計算時間均增長.結(jié)果表明，SIMPLE算法應(yīng)用于管網(wǎng)泄漏模擬中，計算速度最慢，伴隨方程的PISO算法計算速度最快.

圖18 不同計算方法計算速度對比

我們利用了進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證的管網(wǎng)模型了，進(jìn)行了一段時間的計算，對比了SIMPLE算法，PISO算法和本文加入了伴隨方程的PISO算法，在相同的時間步長下，本文的計算速度明顯加快.

2.3 實(shí)驗(yàn)驗(yàn)證

在此部分通過對不同漏點(diǎn)進(jìn)行實(shí)驗(yàn)，分析管網(wǎng)中的壓力變化規(guī)律，將實(shí)驗(yàn)數(shù)據(jù)同模擬結(jié)果進(jìn)行對比分析.以泄漏點(diǎn)1為例，模擬和實(shí)驗(yàn)壓力變化規(guī)律同模擬結(jié)果符合度較高.

如圖19所示為以泄漏孔1為例的泄漏孔處的壓力變化曲線，圖中藍(lán)色線條代表壓力變化的總體趨勢，黑色線條分別為試驗(yàn)和模擬條件下壓力變化趨勢.從壓力變化曲線中可以發(fā)現(xiàn)兩者變化規(guī)律相同，實(shí)驗(yàn)數(shù)據(jù)在泄漏前壓力值保持在一個高位水平，泄漏瞬間壓力值突降，在到達(dá)最底端后，恢復(fù)值一個較低的壓力值，并逐漸穩(wěn)定.而模擬數(shù)據(jù)與實(shí)驗(yàn)變化規(guī)律相反，這是由于添加了逆向模型的緣故.對比實(shí)驗(yàn)和模擬數(shù)據(jù)，模擬過程阻力小于實(shí)驗(yàn)，無噪聲的影響，因此模擬壓力的最大值高于實(shí)驗(yàn)，而泄漏前后的壓降大小基本相同.總的來說，從橫軸觀察泄漏瞬間對應(yīng)的時間差值可以發(fā)現(xiàn)，泄漏是在一瞬間發(fā)生的壓力變化時間極端，這與模擬的結(jié)果相吻合，并且正向壓力傳播曲線與模擬所得反向壓力傳播曲線呈現(xiàn)的壓力變化趨勢相吻合.

圖19 實(shí)驗(yàn)?zāi)M對比圖

接下來為了進(jìn)一步驗(yàn)證方法的合理性，對與模擬工況對應(yīng)的多組實(shí)驗(yàn)數(shù)據(jù)做了處理，將50組實(shí)驗(yàn)和模擬的工況進(jìn)行了計算，對比了各自的定位準(zhǔn)確度，如圖20所示.

圖20 不同泄漏孔直徑下實(shí)驗(yàn)?zāi)M定位準(zhǔn)確度對比圖圖21 不同泄漏點(diǎn)管道泄漏定位準(zhǔn)確度變化曲線圖

不同泄漏孔直徑下實(shí)驗(yàn)?zāi)M定位準(zhǔn)確度對比圖，如圖20所示.可知隨著泄漏孔直徑增大實(shí)驗(yàn)和模擬點(diǎn)位準(zhǔn)確度均增大，但當(dāng)泄漏孔直徑達(dá)到12 mm時，再增大泄漏孔，由于泄漏孔徑很大，泄漏產(chǎn)生的壓力波經(jīng)過兩個測點(diǎn)時間很短，時間差計算不準(zhǔn)確，進(jìn)而導(dǎo)致定位位置計算不準(zhǔn)確，準(zhǔn)確度有所下降.總的來說，實(shí)驗(yàn)結(jié)果與模擬結(jié)果相近，定位準(zhǔn)確度在2%左右波動.

不同泄漏點(diǎn)管道泄漏定位準(zhǔn)確度變化曲線圖，如圖21所示.實(shí)驗(yàn)管網(wǎng)的形狀以及五個泄漏點(diǎn)的位置一定，泄漏點(diǎn)1和泄漏點(diǎn)5距離兩個壓力測點(diǎn)路徑長度一致，定位準(zhǔn)確度基本一致，泄漏點(diǎn)2距離測1距離相對較遠(yuǎn)，誤差比1、5大，泄漏點(diǎn)3的定位準(zhǔn)確度最低，3點(diǎn)泄漏信號傳遞到兩個壓力測點(diǎn)，無論怎樣傳遞都要兩次經(jīng)過四通，并且距離兩個測點(diǎn)的長度較長，定位準(zhǔn)確度比較低.

總體來說，實(shí)驗(yàn)結(jié)果與模擬結(jié)果基本吻合，均在誤差允許的范圍內(nèi)，該方法可行.另外本研究是基于實(shí)際管網(wǎng)進(jìn)行的簡化，實(shí)際地下供熱管網(wǎng)沒有四通，多用二個臨近三通來代替，壓力波的反射情況完全不同，給復(fù)雜實(shí)際管網(wǎng)檢測中帶來一定的干擾，但是這種干擾只是增加了一部分計算過程中的局部阻力，壓力波傳遞情況稍有改變，對整體壓力變化趨勢以及泄漏檢測和定位精度影響較小.

3 結(jié) 論

本文的主要結(jié)合反向?qū)ぴ丛砑鞍殡S方法，推導(dǎo)瞬態(tài)流場反向伴隨方程.借助MATLAB軟件和開源有限元軟件OpenFOAM，構(gòu)建單管和復(fù)雜管網(wǎng)流體流動數(shù)值模擬模型，在OpenFOAM自定義求解器，將反向?qū)ぴ磫栴}轉(zhuǎn)化為伴隨方程的求解問題.根據(jù)得到的圖形分析不同測點(diǎn)的壓力等參數(shù)的變化，結(jié)合新時延估計方法和互相關(guān)方法，對管網(wǎng)的運(yùn)行狀態(tài)進(jìn)行分析檢測，完成對伴隨方程的驗(yàn)證和求解，確定泄漏源的可能泄漏位置，最終利用全局搜索方法完成了泄漏源的定位.得出以下結(jié)論：

(1)基于正向壓力輸運(yùn)模型，引入目標(biāo)函數(shù)，將伴隨理論與靈敏度分析方法相結(jié)合推導(dǎo)了正向伴隨方程.使敏感度函數(shù)對壓強(qiáng)求導(dǎo)，并引入伴隨算子，推導(dǎo)了瞬態(tài)流場反向伴隨方程，求得了其解析解，并用MATLAB軟件，完成了伴隨模型的驗(yàn)證，結(jié)果證明反向?qū)ぴ捶椒捎糜诠芫W(wǎng)泄漏檢測.

(2)在OpenFOAM中完成了復(fù)雜管網(wǎng)流體流動數(shù)值模擬模型的構(gòu)建和伴隨方程的求解，并監(jiān)測了泄漏孔泄漏瞬間壓力變化，運(yùn)用新時延估計方法銳化了峰值點(diǎn)，最后用互相關(guān)法進(jìn)行了相關(guān)性分析，發(fā)現(xiàn)管道泄漏后，信號相關(guān)系數(shù)的最大值在零附近有所偏移，通過相關(guān)性對比圖可知，管道泄漏的出現(xiàn)，最大相關(guān)系數(shù)的值下降，以泄漏點(diǎn)1為例證明了方法的可行性.

(3)比較了不同泄漏孔大小及不同進(jìn)口壓力下，模擬泄漏定位的準(zhǔn)確度，結(jié)果表明，運(yùn)用基于反向?qū)ぴ捶ǘㄎ坏臏?zhǔn)確度在92%～96%之間，說明此方法可以應(yīng)用在復(fù)雜管網(wǎng)泄漏檢測及定位中.

(4)本文在不同計算步數(shù)時，對比了單一的SIMPLE、PISO及加入伴隨方程的PISO算法所需的迭代時間，結(jié)果表明，結(jié)合伴隨方程的PISO算法的計算速度要明顯優(yōu)于單一的SIMPLE和PISO算法，以此證明了伴隨方法的優(yōu)越性.

在本文中，基于復(fù)雜管網(wǎng)泄漏檢測及定位模型，未考慮同樣引起泄漏的堵塞、變形及污染等問題，應(yīng)定義相應(yīng)的求解器，設(shè)置相應(yīng)的參數(shù)及邊界條件，盡可能將定位方法轉(zhuǎn)化成OpenFOAM語言，便可通過模擬體系，直接得到不同情況下的診斷結(jié)果.此外，實(shí)驗(yàn)系統(tǒng)存在一定的不穩(wěn)定性和局限性，與實(shí)際復(fù)雜管網(wǎng)還有一定差距，只是在宏觀上完成了伴隨方程的求解和驗(yàn)證，需進(jìn)一步結(jié)合實(shí)際工況強(qiáng)化實(shí)驗(yàn)和模擬.