基于Godunov 格式的輸水管道阻塞數(shù)值模擬與識(shí)別

胡朝仲，楊旭，李赟杰，白光國(guó)，周領(lǐng)，4

（1.云南水利水電職業(yè)學(xué)院水利工程學(xué)院，云南昆明 650499；2.中水北方勘測(cè)設(shè)計(jì)研究有限責(zé)任公司，天津 300222；3.河海大學(xué)水利水電學(xué)院，江蘇南京 210098；4.河海大學(xué)長(zhǎng)江保護(hù)與綠色發(fā)展研究院，江蘇南京 210098）

0 引言

引調(diào)水工程和城市供水管網(wǎng)建設(shè)的發(fā)展，大大滿足了生活和生產(chǎn)用水，同時(shí)也暴露出了一些問(wèn)題，如管道泄漏、阻塞等。在實(shí)際輸水管道中，生物膜堆積、腐蝕和沉積過(guò)程會(huì)形成堵塞，且堵塞的大小和嚴(yán)重程度通常隨著時(shí)間的推移而增加［1］。管道阻塞會(huì)導(dǎo)致額外的能量損失和水質(zhì)問(wèn)題，因此管道阻塞的準(zhǔn)確模擬和精準(zhǔn)檢測(cè)十分重要。

以往的管道瞬變流數(shù)值模擬大都基于特征線法求解，該方法簡(jiǎn)單、高效，常用于各類瞬變流現(xiàn)象模擬［2，3］。然而，特征線法在求解復(fù)雜管道模型時(shí)需要插值近似或調(diào)整波速，容易產(chǎn)生數(shù)值耗散或色散。基于此，Zhao等［4］提出一種有限體積法Godunov 求解格式，該求解格式受庫(kù)朗特?cái)?shù)的影響較小。Zhou等［5，6］基于該求解格式，納入了不同類型的動(dòng)態(tài)摩阻，結(jié)果表明Brunone 摩阻模型結(jié)合Godunov 求解格式在低庫(kù)朗特?cái)?shù)工況條件下數(shù)值耗散小，能更好地計(jì)算復(fù)雜管道模型。管道阻塞分為局部阻塞和管段阻塞，局部阻塞即尺度小于最小波長(zhǎng)，可以視為一個(gè)不連續(xù)點(diǎn)，如一個(gè)未完全關(guān)閉的閥門(mén)；管段阻塞的尺度大于最小波長(zhǎng)，可以視為一段直徑較小的管段串聯(lián)在管道中［7］。對(duì)于管道阻塞，特別是管段阻塞，在數(shù)值模擬時(shí)可能涉及到庫(kù)朗特?cái)?shù)小于1的情況，此時(shí)特征線法求解格式會(huì)有較大誤差。

本文提出了一種基于Godunov 求解格式的管道阻塞模型。首先，給出含Brunone 動(dòng)態(tài)摩阻的控制方程和邊界條件，以此建立相應(yīng)的管道阻塞模型。然后，分別對(duì)局部阻塞和管段阻塞進(jìn)行模型驗(yàn)證，并根據(jù)波傳播理論分析不同阻塞工況壓力波動(dòng)曲線突變?cè)颍⒔o出阻塞參數(shù)識(shí)別方法。通過(guò)參數(shù)敏感性分析研究了局部阻塞和管段阻塞的影響因素和結(jié)果。

1 數(shù)學(xué)模型及求解

1.1 控制方程

管道瞬變流的基本控制方程為［2，8，9］：

式中：H為測(cè)壓管水頭；V為管道斷面的平均流速；a為波速；g為重力加速度；x為沿管道軸向長(zhǎng)度；t為時(shí)間；D為管道內(nèi)徑；ρ為水體密度；f為達(dá)西-維斯巴赫系數(shù)；SIGN（V）為符號(hào)函數(shù)；k3為Brunone摩阻系數(shù)［10］。k3通常采用Vardy等［11］所給的系數(shù)公式：

1.2 邊界條件

對(duì)于邊界條件，通過(guò)在邊界兩端各增加兩個(gè)虛擬控制體，可以得到對(duì)邊界控制體和內(nèi)部控制體的統(tǒng)一計(jì)算模擬，并達(dá)到二階精度要求［12］。結(jié)合黎曼問(wèn)題在邊界處的解析解、上下游邊界處的動(dòng)量和質(zhì)量方程可得到邊界處的未知量。

對(duì)于上游恒定水庫(kù)水位HR，流速可以表示為：

對(duì)于下游閥門(mén)-水庫(kù)邊界條件，有：

對(duì)于管道中閥門(mén)或串聯(lián)節(jié)點(diǎn)邊界條件，主要考慮連接邊界處的水頭損失，即：

1.3 Godunov求解格式

對(duì)于有限體積法Godunov 求解格式，首先將瞬變流基本方程以黎曼問(wèn)題的求解格式表示為矩陣形式［4-6］：

式中：u為求解向量；f（u）為控制體通量矢量；Aˉ為矢量f中各變量對(duì)u中各求解變量的偏導(dǎo)數(shù)矩陣；s（u）為源項(xiàng)，代表了管道的摩阻項(xiàng)。

圖1 沿管道軸向的計(jì)算網(wǎng)格

式中：fi-1/2表示i-1/2界面的數(shù)值通量；fi+1/2表示i+1/2界面的數(shù)值通量。因此，在計(jì)算下一時(shí)步的向量U時(shí)，需涉及當(dāng)前時(shí)步的控制體兩邊界數(shù)值通量和源項(xiàng)。

針對(duì)Godunov 格式，控制體邊界數(shù)值通量可由局部黎曼問(wèn)題得到，可描述為初值問(wèn)題：

對(duì)于和，采用分段多項(xiàng)式重構(gòu)近似，其精度取決于多項(xiàng)式的形式。對(duì)于二階Godunov 近似，通過(guò)引入MUSCLHancock 格式進(jìn)行線性插值重構(gòu)。因?yàn)槎A格式會(huì)在高梯度附近產(chǎn)生假振蕩，采用MINMOD 斜率限制器削減這種數(shù)值震蕩［6，13］。因此，可以得到控制體積各界面的通量向量。源項(xiàng)采用顯式的二階龍格-庫(kù)塔離散化方法，考慮到Courant-Friedrichs-Lewy（CFL）定律和龍格庫(kù)塔離散穩(wěn)定性條件，可得到源項(xiàng)的允許時(shí)間步長(zhǎng)［6］。

2 模型驗(yàn)證

針對(duì)局部阻塞、管段阻塞兩種堵塞類型對(duì)本文所建模型進(jìn)行驗(yàn)證、討論和分析。

2.1 局部阻塞模型驗(yàn)證

將模擬結(jié)果與Lee等［14］2015 年中的論文數(shù)據(jù)作對(duì)比，該系統(tǒng)上游為水庫(kù)，下游閥門(mén)，管道中間某處有一未完全關(guān)閉閥門(mén)。具體參數(shù)為：管道長(zhǎng)度L=2 000 m，上游恒定水位HR=50 m，中間閥門(mén)位于距上游水庫(kù)0.3L處，穩(wěn)態(tài)損失系數(shù)K=1 890，波速a=1 000 m/s，管道直徑D=0.5 m，初始流量Q0=0.02 m3/s，穩(wěn)態(tài)摩阻系數(shù)f=0.02。閥門(mén)瞬間關(guān)閉，模型模擬結(jié)果對(duì)比如圖2 所示。圖2中，對(duì)縱坐標(biāo)壓力進(jìn)行標(biāo)準(zhǔn)化處理，即H*=（H-H0）/（aV0/g），其中H*為無(wú)量綱水頭，H0為穩(wěn)態(tài)壓力，V0為穩(wěn)態(tài)流速。

由圖2 可知，本文提出的基于Godunov 求解格式的管道局部阻塞模型與Lee 等2015 年模型模擬結(jié)果基本重合，驗(yàn)證了局部阻塞模型的正確性。此外，根據(jù)模擬結(jié)果可知，對(duì)于管道局部阻塞，其壓力波動(dòng)曲線在第一峰值處會(huì)突然增加，這是由于壓力波從末端閥門(mén)傳播到阻塞處時(shí)會(huì)發(fā)生反射，反射的正波與末端閥門(mén)突然關(guān)閉產(chǎn)生的主波疊加，導(dǎo)致反射波傳播到閥門(mén)末端時(shí)壓力驟增，由此可以根據(jù)壓力驟增的時(shí)間判斷局部阻塞發(fā)生位置。在本算例中，第一峰值壓力驟增的時(shí)刻為t=2.8 s，根據(jù)波速a=1 000 m/s，可以得到阻塞點(diǎn)距下由閥門(mén)位置xl=at/2=1 400 m，與實(shí)際阻塞點(diǎn)位置一致。

圖2 局部阻塞模型驗(yàn)證對(duì)比圖

2.2 管段阻塞模型驗(yàn)證

將模擬結(jié)果與Duan等［15］實(shí)驗(yàn)數(shù)據(jù)作對(duì)比，該系統(tǒng)上游為水庫(kù)，下游閥門(mén)，管道分為三段，中間段連接直徑較小的管道。具體參數(shù)為：管道總長(zhǎng)度L=41.5 m，其中L1=8.8 m，L2=12.1 m，L3=20.6 m；上游恒定水位HR=38.2 m，波速a=1 180 m/s，管道分別為直徑D1=0.072 4 m，D2=0.022 25 m，D3=0.072 4 m，初始流量Q0=9.7×10-5m3/s，穩(wěn)態(tài)摩阻系數(shù)f=0.045 9，閥門(mén)關(guān)閉時(shí)間為0.01 s。模型與實(shí)驗(yàn)對(duì)比結(jié)果如圖3 所示，圖3（a）和（b）分別為穩(wěn)態(tài)摩阻與動(dòng)態(tài)摩阻模型模擬結(jié)果。

由圖3 可知，提出的管段阻塞模型能較好地復(fù)現(xiàn)壓力波動(dòng)曲線周期和衰減。圖3（a）所示，僅考慮穩(wěn)態(tài)摩阻時(shí)，在第一峰值處與實(shí)驗(yàn)結(jié)果基本重合，但隨著時(shí)間變化，最大與最小壓力值和實(shí)驗(yàn)結(jié)果偏差越大。圖3（b）所示，考慮動(dòng)態(tài)摩阻模型，隨著時(shí)間變化仍能較好的模擬實(shí)驗(yàn)最大最小壓力值，僅會(huì)發(fā)生微小的相位偏移。與局部阻塞一致，末端閥門(mén)產(chǎn)生的壓力波在管段阻塞處也會(huì)發(fā)生反射。當(dāng)管徑突然變小時(shí)，反射的壓力波為正波，管徑突然變大時(shí)，反射的壓力波為負(fù)波。當(dāng)正波或負(fù)波傳播到末端閥門(mén)時(shí)，會(huì)產(chǎn)生相應(yīng)的壓力增加或減小，因此在本結(jié)果前半個(gè)周期中會(huì)出現(xiàn)壓力波先增加后減小的現(xiàn)象。同理，據(jù)此可以計(jì)算管段阻塞的位置和長(zhǎng)度。以本工況為例，前半個(gè)周期內(nèi)壓力波驟增和驟減的時(shí)間分別為t1=0.035 s，t2=0.055 s，根據(jù)波速a=1 180 m/s，可以得到阻塞管段末端距下游閥門(mén)20.65 m，阻塞長(zhǎng)度為11.8 m，這與實(shí)際阻塞位置和長(zhǎng)度基本一致。

圖3 管段阻塞模型驗(yàn)證對(duì)比圖

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

3.1 局部阻塞

采用2.1 節(jié)局部阻塞工況分別對(duì)管道初始流量Q0、局部阻塞損失系數(shù)K和局部阻塞位置xl進(jìn)行參數(shù)敏感性分析，結(jié)果如圖4所示。

對(duì)管道初始流量Q0參數(shù)敏感性分析，分別采用0.02、0.04和0.08 m3/s，結(jié)果如圖4（a）所示。圖4（a）表明，阻塞大小不變的情況下，隨著初始流量的增加，在第一峰值由阻塞引起的反射壓力幅值隨之增大，即反射系數(shù)隨初始流量增大而增大。Yan［16］理論推導(dǎo)了管道局部阻塞引起的反射系數(shù)解析解，反射系數(shù)r可以表示為：

式中：A為管道橫截面積。

由解析解可知：其他條件不變的情況下，隨著管道初始流量的增大，反射系數(shù)也會(huì)增大，且增大的幅度與流量呈正比。這也證明了模擬的反射系數(shù)與解析解結(jié)果一致。此外，隨著流量的增加，壓力曲線的峰谷壓力也會(huì)增大，衰減速度更快。

對(duì)局部阻塞損失系數(shù)K進(jìn)行參數(shù)敏感性分析，K值分別取1 890、3 000 和5 000，結(jié)果如圖4（b）所示。所得規(guī)律與管道初始流量類似，具體表現(xiàn)為：在其他條件不變的情況下，隨著K值的增大（即局部阻塞程度增加），峰值反射增大，且反射系數(shù)與K呈正比（與解析解一致）；峰谷壓力增大，衰減速度加快。

對(duì)局部阻塞位置xl進(jìn)行參數(shù)敏感性分析，局部阻塞位置xl分別距上游水庫(kù)0.3L、0.5L和0.7L，結(jié)果如圖4（c）所示。由圖4（c）可知，不同局部阻塞位置在前半個(gè)周期內(nèi)表現(xiàn)為壓力上升時(shí)間不同，即發(fā)生反射的時(shí)間不同。因此，可根據(jù)壓力波反射的時(shí)間估計(jì)局部阻塞位置，具體見(jiàn)3.1節(jié)。

圖4 局部阻塞參數(shù)敏感性分析

3.2 管段阻塞

對(duì)管段阻塞進(jìn)行參數(shù)敏感性分析，系統(tǒng)上游為水庫(kù)，下游閥門(mén)，管道分為三段，中間段連接直徑較小的管道。具體參數(shù)為：管道總長(zhǎng)度L=300 m，其中L1、L2、L3均為100 m；上游恒定水位HR=50 m，波速a=1 000 m/s，管道直徑分別為D1=0.2 m，D2=0.14 m，D3=0.2 m，初始流量Q0=0.01 m3/s，穩(wěn)態(tài)摩阻系數(shù)f=0.02，閥門(mén)瞬間關(guān)閉。分別對(duì)管段阻塞程度和阻塞長(zhǎng)度進(jìn)行參數(shù)敏感性分析，結(jié)果如圖5所示。

對(duì)管段阻塞程度進(jìn)行參數(shù)敏感性分析，保持其他參數(shù)不變，中間管道直徑D2分別選取0.14、0.16、0.18和0.2 m 進(jìn)行計(jì)算分析，結(jié)果如圖5（a）所示。圖5（a）表明，隨著阻塞管段直徑的減小，阻塞反射幅值增加，表現(xiàn)為第一壓力峰值增大，且系統(tǒng)最大峰值壓力不再為第一峰值壓力。對(duì)于隨時(shí)間變化的壓力波動(dòng)曲線，阻塞管段直徑減小使得系統(tǒng)壓力波動(dòng)周期增大，這是由于布拉格共振效應(yīng)引起的［17，18］。

對(duì)管段阻塞長(zhǎng)度進(jìn)行參數(shù)敏感性分析，保持其他參數(shù)不變，第三段管道長(zhǎng)度恒定100 m，第一段管道長(zhǎng)度和中間管道長(zhǎng)度之和保持200 m 不變，中間管道長(zhǎng)度L2分別選取20、60 和100 m 進(jìn)行計(jì)算分析，結(jié)果如圖5（b）所示。圖5（b）表明，隨著管段阻塞長(zhǎng)度增加，第一峰值壓力下降時(shí)間延長(zhǎng)，這是因?yàn)樗矐B(tài)波穿過(guò)阻塞管段時(shí)間變短。波動(dòng)周期隨管段阻塞長(zhǎng)度增加而逐漸增大，與阻塞管段直徑變化引起的相移一致，這也是由布拉格共振效應(yīng)引起的。

圖5 管段阻塞參數(shù)敏感性分析

4 結(jié)論

本文首先提出了基于Godunov 求解格式的管道阻塞模型，并針對(duì)局部阻塞和管段阻塞分別進(jìn)行模型驗(yàn)證和結(jié)果分析，主要結(jié)論有：①本文提出的基于Godunov 求解格式的管道阻塞模型精確度高，能較好的模擬管道阻塞壓力軌跡；②根據(jù)壓力波傳播理論，可根據(jù)壓力的驟增或驟減計(jì)算管道中局部阻塞的位置以及管段阻塞的位置和長(zhǎng)度；③局部阻塞壓力峰值大小受管道初始流量、局部阻塞損失系數(shù)的影響，壓力峰值出現(xiàn)的時(shí)間受局部阻塞位置的影響；④管段阻塞程度影響壓力峰值和波動(dòng)周期大小，管段阻塞長(zhǎng)度影響壓力峰值出現(xiàn)的時(shí)間和波動(dòng)周期。

同時(shí)應(yīng)當(dāng)注意到，在計(jì)算阻塞位置或長(zhǎng)度時(shí)，壓力變化點(diǎn)的時(shí)間不精確會(huì)導(dǎo)致誤差，利用一些信號(hào)處理技術(shù)如小波變換等可能會(huì)提高計(jì)算精度。