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

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

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

0 引言

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

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

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

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

1.1 控制方程

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

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

1.2 邊界條件

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

對于上游恒定水庫水位HR，流速可以表示為：

對于下游閥門-水庫邊界條件，有：

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

1.3 Godunov求解格式

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

式中：u為求解向量；f（u）為控制體通量矢量；Aˉ為矢量f中各變量對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ì)算下一時步的向量U時，需涉及當(dāng)前時步的控制體兩邊界數(shù)值通量和源項(xiàng)。

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

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

2 模型驗(yàn)證

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

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

將模擬結(jié)果與Lee等［14］2015 年中的論文數(shù)據(jù)作對比，該系統(tǒng)上游為水庫，下游閥門，管道中間某處有一未完全關(guān)閉閥門。具體參數(shù)為：管道長度L=2 000 m，上游恒定水位HR=50 m，中間閥門位于距上游水庫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。閥門瞬間關(guān)閉，模型模擬結(jié)果對比如圖2 所示。圖2中，對縱坐標(biāo)壓力進(jìn)行標(biāo)準(zhǔn)化處理，即H*=（H-H0）/（aV0/g），其中H*為無量綱水頭，H0為穩(wěn)態(tài)壓力，V0為穩(wěn)態(tài)流速。

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

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

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

將模擬結(jié)果與Duan等［15］實(shí)驗(yàn)數(shù)據(jù)作對比，該系統(tǒng)上游為水庫，下游閥門，管道分為三段，中間段連接直徑較小的管道。具體參數(shù)為：管道總長度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，閥門關(guān)閉時間為0.01 s。模型與實(shí)驗(yàn)對比結(jié)果如圖3 所示，圖3（a）和（b）分別為穩(wěn)態(tài)摩阻與動態(tài)摩阻模型模擬結(jié)果。

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

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

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

3.1 局部阻塞

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

對管道初始流量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ù)也會增大，且增大的幅度與流量呈正比。這也證明了模擬的反射系數(shù)與解析解結(jié)果一致。此外，隨著流量的增加，壓力曲線的峰谷壓力也會增大，衰減速度更快。

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

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

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

3.2 管段阻塞

對管段阻塞進(jìn)行參數(shù)敏感性分析，系統(tǒng)上游為水庫，下游閥門，管道分為三段，中間段連接直徑較小的管道。具體參數(shù)為：管道總長度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，閥門瞬間關(guān)閉。分別對管段阻塞程度和阻塞長度進(jìn)行參數(shù)敏感性分析，結(jié)果如圖5所示。

對管段阻塞程度進(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)最大峰值壓力不再為第一峰值壓力。對于隨時間變化的壓力波動曲線，阻塞管段直徑減小使得系統(tǒng)壓力波動周期增大，這是由于布拉格共振效應(yīng)引起的［17，18］。

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

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

4 結(jié)論

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

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