胡朝仲，楊旭，李赟杰，白光國，周領，4

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

0 引言

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

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

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

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

1.1 控制方程

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

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

1.2 邊界條件

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

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

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

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

1.3 Godunov求解格式

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

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

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

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

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

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

2 模型驗證

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

2.1 局部阻塞模型驗證

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

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

圖2 局部阻塞模型驗證對比圖

2.2 管段阻塞模型驗證

將模擬結果與Duan等［15］實驗數(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，閥門關閉時間為0.01 s。模型與實驗對比結果如圖3 所示，圖3（a）和（b）分別為穩(wěn)態(tài)摩阻與動態(tài)摩阻模型模擬結果。

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

圖3 管段阻塞模型驗證對比圖

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

3.1 局部阻塞

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

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

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

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

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

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

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

3.2 管段阻塞

對管段阻塞進行參數(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，閥門瞬間關閉。分別對管段阻塞程度和阻塞長度進行參數(shù)敏感性分析，結果如圖5所示。

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

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

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

4 結論

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

同時應當注意到，在計算阻塞位置或長度時，壓力變化點的時間不精確會導致誤差，利用一些信號處理技術如小波變換等可能會提高計算精度。