考慮動態(tài)摩阻的管道泄漏水力瞬變數(shù)值模擬

李赟杰，周 領(lǐng)，歐傳奇，劉德有

(1.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098； 2.國際小水電中心，浙江 杭州 310002)

管道漏損廣泛存在于流域輸水、城市供水等工程中，泄漏不僅造成巨額的經(jīng)濟損失，而且可能造成次生災(zāi)害，如影響水質(zhì)、損壞建筑基礎(chǔ)等，因此尋求準(zhǔn)確有效的泄漏檢測方法是一個重要的研究課題。目前，泄漏檢測方法主要有目視觀測法、探地雷達(dá)檢漏法、示蹤劑檢漏法、振動聲學(xué)檢漏法、光纖傳感技術(shù)法、瞬態(tài)檢漏法[1]以及穩(wěn)態(tài)水力模型法等，其中，瞬態(tài)檢漏法是泄漏檢測的主要手段[2-3]。泄漏主要通過增加阻尼率和在產(chǎn)生的壓力中創(chuàng)建新的反射信號兩個方面影響瞬態(tài)信號，識別和量化這些影響是瞬態(tài)泄漏檢測技術(shù)的核心[4]。

Pudar等[1]在1992年第一次提出基于壓力測量的檢漏程序和逆問題的解決方法，即已知壓力和流量求管道粗糙度和泄漏情況。誤差分析表明系統(tǒng)摩擦系數(shù)估計的誤差帶來的影響比水頭測量誤差的影響大得多，摩擦系數(shù)的精確性是逆瞬態(tài)法的關(guān)鍵[5-7]。目前，為了提高泄漏檢測的精確度，除了較大泄漏問題，水力瞬變模型中一般均需要考慮動態(tài)摩阻的影響[2,8-10]。Zielke[11]推導(dǎo)出適用于管道瞬變層流的卷積加權(quán)模型，該模型考慮了管道瞬時平均流速和加權(quán)歷史流速變化，被認(rèn)為是精確動態(tài)摩阻模型；Vardy等[12-14]在Zielke[11]的基礎(chǔ)上建立了光滑管和粗糙管內(nèi)流動的瞬態(tài)摩擦加權(quán)函數(shù)模型，將卷積加權(quán)模型推廣到阻力平方區(qū)；Brunone等[15]在Daily等[16]研究基礎(chǔ)上考慮對流加速度，引入動態(tài)摩阻系數(shù)，提出了適用于層流和紊流的動態(tài)摩阻模型，該模型能很好地平衡運算量和計算精度的問題；Bergant等[17]和Vitkovsky等[18]對Brunone模型分別進行了公式改進和近似求解，使Brunone模型更加簡潔、適用。郭新蕾等[3]根據(jù)實際工程需要研究了Brunone模型的離散網(wǎng)格和求解方法，將其應(yīng)用于管道泄漏模擬，壓力波形對比表明該模型比傳統(tǒng)的準(zhǔn)恒定泄漏摩阻模型更加精確。然而，加權(quán)函數(shù)類動態(tài)摩阻模型很少用于管道泄漏檢測的研究。

本文將瞬態(tài)加速類(Brunone模型)和加權(quán)函數(shù)類(Zielke模型和Vardy模型)兩類動態(tài)摩阻模型用于水力瞬變泄漏建模和模擬，并將計算結(jié)果與已有的試驗數(shù)據(jù)和模擬結(jié)果進行對比，采用3種模型進行了算例計算分析，通過對比不同瞬變流模型泄漏過程中的壓力和流量變化情況，分析泄漏畸變點突變原因、Brunone模型計算誤差和泄漏中部分參數(shù)的敏感性。

1 動態(tài)摩阻模型

常見的動態(tài)摩阻模型主要有以下兩類[10,19]：①基于加權(quán)函數(shù)考慮歷史速度和歷史加速度對當(dāng)前流態(tài)影響的加權(quán)函數(shù)類模型，如Zielke模型和Vardy模型；②基于瞬時加速度的模型，如Brunone模型及其改進型。本文將Zielke模型、Vardy模型和Brunone模型用于管道泄漏檢測的水力瞬態(tài)建模。

1.1 Zielke模型

Zielke[11]在考慮管道瞬時平均流速和加權(quán)歷史流速變化的情況下，推導(dǎo)出了適用于層流的動態(tài)摩阻項JU：

(1)

1.2 Vardy模型

Vardy等[12]建立了光滑管內(nèi)流動的瞬態(tài)摩擦加權(quán)函數(shù)模型，并假定湍流黏度在均勻流速核心周圍的剪切層內(nèi)呈線性變化，對于高雷諾數(shù)的流動，給出了水力光滑管道動態(tài)摩阻模型：

(2)

1.3 Brunone模型及其改進型

(3)

式中：x為沿管道中心線方向的距離；a為波速；k3為Brunone動態(tài)摩阻系數(shù)，取值的適當(dāng)性對模型精確程度有很大影響，k3的確定方法有兩種：一是由試驗確定；二是采用Vardy等[12]所給的系數(shù)公式：

(4)

Bergant等[17]發(fā)現(xiàn)Brunone模型動態(tài)摩阻項計算公式具有局限性：對于某些加速或減速條件，該公式對流項符號錯誤，因此他們在其基礎(chǔ)上提出了改進模型：

(5)

Vitkovsky等[18]隨后給出了該模型的顯性和隱性近似解法。

2 考慮動態(tài)摩阻的管道泄漏數(shù)學(xué)模型

管道泄漏的數(shù)學(xué)模型由控制方程、動態(tài)摩阻和邊界條件3部分組成。本文采用特征線法對管道泄漏模型進行求解。動態(tài)摩阻的求解一般包含隱式和顯式求解兩種格式[3]：隱式求解格式中，將穩(wěn)態(tài)摩阻項和動態(tài)摩阻項共同作為管道流中摩擦引起的水頭損失代入雙曲型微分方程，這種方法改變了特征線網(wǎng)格，運算復(fù)雜且求解時間長，優(yōu)點是更為精確；顯式求解格式中，將動態(tài)摩阻作為附加項單獨計算，可以較好地代入各種模型，模型簡單且計算時間短。綜合考慮模型精度和運算效率，本文采用第二種方法。

2.1 控制方程

管道瞬變流的連續(xù)性方程和動量守恒方程[20-22]如下：

(6)

(7)

式中：H為測壓管水頭；J為管道水體運動摩擦引起的單位長度的水頭損失；f為Darcy-Weisbach摩擦系數(shù)，其值大小依賴于雷諾數(shù)和管道相對粗糙度。

如圖1所示，將線性雙曲型微分方程(6)(7)變換為常微分方程，采用一階有限差分近似，分別沿著正負(fù)特征線進行積分，得到如下C+、C-特征線方程：

(8)

(9)

式中：vP、HP分別為在t時刻節(jié)點P處的流速和水頭；vR、HR分別為在t-Δt時刻節(jié)點R處的流速和水頭；vS、HS分別為在t-Δt時刻節(jié)點S處的流速和水頭；JR、JS分別為在t-Δt時刻節(jié)點R、S處的水頭損失；Δt為計算時間間隔。

圖1 特征線積分網(wǎng)格

2.2 邊界條件

上游端為已知水位的水庫，下游為一閥門，泄漏孔口邊界條件如圖2所示。

圖2 泄漏孔邊界條件

假定泄漏發(fā)生在第i個節(jié)點處，泄漏孔口外部為大氣壓，正負(fù)特征線測壓管水頭H1、H2與泄漏處測壓管水頭相等。由孔口出流公式得泄漏孔口的流量為

(10)

式中：Cd為泄漏孔口的流量系數(shù)，可由試驗確定具體值，取值范圍一般在0.60～0.61之間[23]；Ag為泄漏孔口面積；ΔHL為泄漏點的壓力水頭，ΔHL=HL-zL，其中HL為管道泄漏處的測壓管水頭，zL為管道泄漏點高程。

由連續(xù)性方程可得：

Q1=Q2+QL

(11)

式中：Q1、Q2分別為泄漏節(jié)點處的正向特征線和負(fù)向特征線對應(yīng)的流量。將式(10)(11)與式(8)(9)聯(lián)立求解，可得流量Q1、Q2、QL和測壓管水頭H1、H2、HL。

3 模型驗證

3.1 無泄漏純水錘情況

不考慮泄漏情況，將本文動態(tài)摩阻模型(Zielke模型、Vardy模型和Brunone模型)的模擬結(jié)果與Bergant等[17]的實測數(shù)據(jù)進行對比。試驗裝置基本參數(shù)如下：管道長L=37.23 m，直徑D=0.022 1 m，上游水箱恒定水頭h=32 m，波速a=1 319 m/s，管道初始流速v0=0.3 m/s，雷諾數(shù)Re=5 600。閥門瞬間關(guān)閉，模擬的閥門處壓力(用測壓管水頭Hv表示)如圖3所示。由圖3可知，與不考慮動態(tài)摩阻的傳統(tǒng)摩阻模型相比，動態(tài)摩阻模型能精準(zhǔn)描述瞬變管道壓力幅值和衰減，隨著周期數(shù)增加，動態(tài)摩阻模型的模擬值與實測值仍能吻合較好，而傳統(tǒng)摩阻模型模擬的壓力波在兩個周期后誤差逐漸增大，相位偏移加大。

圖3 不同摩阻模型模擬的閥門處測壓管水頭與實測值對比

3.2 存在泄漏情況

考慮管道泄漏情況，將本文模擬結(jié)果與Lee等[24]模擬的結(jié)果進行對比，其模擬參數(shù)為：L=2 000 m，波速a=1 000 m/s，D=0.5 m，h=50 m，初始管道流量Q0=0.02 m3/s，Re=50 000，f=0.02，泄漏孔口距上游水庫300 m，初始泄漏量QL=0.2Q0。不考慮動態(tài)摩阻的影響，模擬結(jié)果見圖4。從圖4可以看出，兩次模擬的壓力幅值和衰減幾乎完全一致，從而證明了本文模型的正確性。

圖4 泄漏情況下閥門處壓力對比

4 算例分析

如圖5所示，選擇水庫-管道-閥門(RPV)模型進行管道泄漏瞬變模擬。假定管道為一水平直管，模型參數(shù)為：L=2 000 m，D=0.2 m，上游水庫水位恒定為60 m，a=1 000 m/s，Q0=0.02 m3/s，Re=85 000，f=0.03。

圖5 泄漏管道概化物理模型

4.1 畸變點分析

假設(shè)泄漏孔外界為大氣壓強，泄漏位置(距上游水庫的距離)XL=1 000 m，QL=0.1Q0。節(jié)點數(shù)N=100，t=0時末端閥門瞬間關(guān)閉，計算時間為32 s，分別模擬了該泄漏情況下的閥門處壓力Hv和管道入口處流量Q，結(jié)果如圖6所示。

圖6 有泄漏時不同摩阻模型模擬的閥門處壓力和管道入口處流量

由圖6可知，管道出現(xiàn)泄漏時，管中壓力和流量會出現(xiàn)畸變。相比于流量，壓力畸變更加容易觀察，因此使用壓力畸變信號能更好地探尋泄漏點。

以泄漏時閥門處的壓力變化為例，閥門關(guān)閉前，管道中水體處于恒定流狀態(tài)，水庫水位保持恒定，壓強為p0，當(dāng)閥門突然關(guān)閉時，緊靠閥門處的水體流速瞬間變?yōu)榱?，壓強立刻升高Δp，上游水體碰到靜止不動的水體時就像碰到完全關(guān)閉的閥門一樣，類似狀態(tài)以波速a向上游傳播，這種波稱為主壓力波，主壓力波對管道壓力和流量變化起主要作用。當(dāng)主壓力波到達(dá)泄漏點時，泄漏處的壓力發(fā)生突然變化，在泄漏點形成向泄漏兩側(cè)傳播的次壓力波，次壓力波也會影響管道壓力變化。圖7給出了管道泄漏時該工況下不同時刻的壓力波運動圖形。

由圖7(a)可知t=2 s時閥門處壓力畸變的原因為：主壓力波從閥門向上游傳播，經(jīng)過泄漏點時形成新的次壓力波。主壓力波從閥門傳遞至泄漏點前引起水體壓力增大，當(dāng)這種增壓波遇到泄漏點(假定為大氣壓)時，由于壓差作用會向閥門處反射一個引起水體壓力減小的減壓信號，在t=2 s時到達(dá)閥門，引起閥門處壓力突降。

圖7 泄漏時不同時刻管道的壓力波運動

由圖7(b)可知t=4 s時閥門處壓力驟降的原因為：增壓波從閥門向水庫傳播，引起逐段的水體壓力增大，傳播至泄漏點時，不僅有向閥門處的減壓波，也會產(chǎn)生向水庫方向的次壓波，這種波也為減壓波。由于主壓力波起主導(dǎo)作用，主增壓波和次減壓波疊加仍然會使泄漏點上游的管道壓力增加，傳播至水庫時，由于壓力差Δp的作用，水庫向閥門反射一個使下游壓力降低的波，在t= 3 s時又到達(dá)泄漏處產(chǎn)生新的次波，此時仍有從泄漏點反射的其他次波，同樣的疊加作用使得泄漏點下游管道壓力降低，t= 4 s時到達(dá)閥門導(dǎo)致閥門處壓力驟降。

4.2 Brunone模型計算誤差分析

Zielke模型主要適用于層流，也可用于雷諾數(shù)不高(Re<106)的湍流；Vardy水力光滑管道動態(tài)摩阻模型[12]適合于湍流情況，Vardy改進模型[13]給出了雷諾數(shù)在2×103～108情況下更為精確的系數(shù)計算公式；Brunone模型在層流和紊流時均可適用。Zielke模型和Vardy模型通常被認(rèn)為是精確動態(tài)摩阻模型，Brunone模型有簡化，存在一定的簡化誤差[25]。在XL=1 000 m，QL=0.1Q0工況下得到壓力如圖8所示，其中圖8(a)為3種動態(tài)摩阻模型模擬的閥門處實際壓力對比，圖8(b)是以Vardy模型為標(biāo)準(zhǔn)，3種模型模擬的閥門處的壓力分別減去Vardy模型的模擬值所得到的壓力對比曲線。

圖8 不同動態(tài)摩阻模型模擬的泄漏時閥門處壓力對比曲線

由圖8(a)可知，Brunone模型在前半個周期內(nèi)所模擬的壓力幅值和衰減與Zielke模型、Vardy模型基本一致，隨著時間增長，因泄漏導(dǎo)致的壓力畸變值雖趨于一致，但其整體壓力衰減值小于Zielke模型和Vardy模型，兩個周期后的壓力曲線與Zielke模型和Vardy模型產(chǎn)生明顯差異。圖8(b)可以更好地印證圖8(a)中所觀察到的現(xiàn)象，Brunone模型在t=4 s時開始與Zielke模型和Vardy模型產(chǎn)生明顯偏離，偏離值隨時間增長而增大，并每隔半個周期出現(xiàn)一次較大的異變點，Zielke模型和Vardy模型在圖8(a)(b)中均無明顯差異?？赡艿脑蚴呛喕腂runone模型有一定的模型誤差，Brunone模型的系數(shù)k3對模型的精確程度有較大影響，需要現(xiàn)場率定。

4.3 參數(shù)敏感性分析

圖9 不同泄漏量影響下的管道壓力、流量變化

以Vardy模型為例，考慮泄漏模型在不同泄漏量和不同泄漏位置情況下的管道壓力、流量變化。模型參數(shù)不變，泄漏位置距上游水庫1 000 m，初始泄漏量QL分別取0.01Q0、0.05Q0、0.1Q0和0.2Q0，圖9比較了4種工況下閥門處壓力Hv、泄漏處壓力HL、管道入口處流量Q和泄漏處流量Q1。由圖9可見，初始泄漏量越大，閥門起始壓力越小，這可由水錘基本公式驗證。隨著泄漏量增加，壓力和流量阻尼增大，收斂速度加快，壓力畸變點隨著泄漏量增加變得越來越清晰，畸變量與初始泄漏量呈正相關(guān)，與時間呈負(fù)相關(guān)。

為了分析不同泄漏位置情況下的管道壓力、流量變化，固定QL=0.1Q0不變，泄漏位置XL分別取距上游水庫500 m、1 000 m、1 500 m，模擬結(jié)果如圖10所示。由圖10可知，不同泄漏位置會導(dǎo)致壓力和流量畸變點位置變化，以閥門壓力第一周期為例，泄漏點越靠近閥門，畸變越早，畸變時間基本與泄漏距閥門距離成正比，據(jù)此可根據(jù)畸變時間反推泄漏位置，其他壓力和流量畸變情況類似，可由壓力波運動情況推得；泄漏量相同，泄漏位置不同時壓力和流量的幅值不同，阻尼不同，泄漏位置越靠近閥門，壓力和流量曲線幅值越小，衰減越快。

圖10 不同泄漏位置影響下的管道壓力、流量變化

5 結(jié) 論

a. 動態(tài)摩阻模型比傳統(tǒng)準(zhǔn)恒定摩阻模型更加精確，能更好地模擬管道泄漏時管中各點的壓力和流量變化。與Zielke和Vardy動態(tài)摩阻模型相比，Brunone模型存在一定的簡化誤差，模型系數(shù)須現(xiàn)場率定。

b. 與流量畸變信號相比，在管道中用壓力畸變信號探尋泄漏點更加方便快捷。泄漏參數(shù)敏感性分析發(fā)現(xiàn)，管道泄漏量與泄漏位置均對管道壓力和流量信號產(chǎn)生影響。