黏彈性管道水柱分離彌合水錘有限體積法模型1)

王 寧 周 領 , 李赟杰 潘天文

* (河海大學水利水電學院,南京 210098)

? (河海大學長江保護與綠色發展研究院,南京 210098)

** (長江勘測規劃設計研究有限責任公司,武漢 430010)

引言

水泵的啟停、閥門開關以及水庫水位波動等情況都可能引起管道負壓,導致沿線某些位置的壓強低于水體的汽化壓強,此時管道中會產生蒸汽空穴,其經歷生成、增長、聚集、壓縮至潰滅過程被稱為水柱分離彌合現象.空穴潰滅水柱彌合過程會生成高速射流、沖擊波等[1],空穴的大部分能量會在此時轉化為沖擊波的波能[2],分離的水柱重新彌合將產生比正常瞬變流水錘大許多的升壓,因此對水柱分離彌合水錘的研究十分重要[3].近年來,黏彈性管材(如PE,HDPE,PVC 等)越來越多地應用到輸水線路和供水系統中,準確地模擬黏彈性管道中水力瞬變現象可為基于瞬變流理論管道泄漏堵塞檢測等前沿技術提供數據,為新型管材管道設計提供依據,降低工程造價并且有效預防水錘危害.

針對長距離輸水管道中水柱分離彌合水錘現象,國內外學者進行了大量模擬研究.Streeter[4]首次提出離散蒸汽空穴模型(discrete vapor cavity model,DVCM),并用于預測水柱分離過程中瞬變壓力,此后許多學者[5-8]進行了研究和討論,認為該模型在計算網格數增加時會產生不真實的壓力峰值.Wylie[9]基于理想氣體方程,考慮計算區段自由氣體總質量來模擬均勻混合物中分布在液體中的自由氣體,提出在管道任意位置空穴體積隨壓力而變化的離散氣體空穴模型(discrete gaseous cavity model,DGCM),并認為該模型比DVCM 模擬更準確.Soares 等[10-11]將DGCM 模型用于高密度聚乙烯管道HDPE 中水柱分離彌合瞬變流現象的研究,對瞬變過程中單相流和空化流進行試驗驗證,并對黏彈性管道瞬變流模型的元件松弛時間及蠕變柔量等參數進行了校核.相關研究工作主要結合特征線法(method of characteristics,MOC)進行計算,有限體積法Godunov 格式因對間斷、水躍和激波等結構具有良好的捕捉能力和效率,且容易擴展到高階精度,近年來被廣泛應用于水力瞬變的模擬,文獻[12-13]基于有限體積法(finite volume method,FVM)二階Godunov 格式對瞬變流問題進行模擬,并指出在相同計算精度下,與有限差分法中特征線法相比,FVM二階Godunov 格式具有計算高效、穩定的優點對于水力瞬變過程數值模擬中常見庫朗數小于1 的情況,更能體現FVM 的優越性.Zhou 等[14]通過引入壓力修正系數,導出彈性管道中有限體積法的DGCM模型,成功消除了該模型在MOC 格式中因插值而出現的虛假數值振蕩.

本文采用有限體積法二階Godunov 格式對離散氣體空穴模型進行研究分析,在傳統的彈性管道水柱分離彌合水錘模型基礎上考慮黏彈性效應的影響.將所建FVM 黏彈性瞬變流模型與試驗結果以及僅考慮管道彈性的瞬變流模型進行比較,并對動態摩阻、庫朗數以及初始含氣率等影響因素進行了分析.

1 數學模型

有限體積法求解DGCM 模型需假定[13]自由氣體集中在控制體中心,形成離散氣體空穴;相鄰空穴之間為純液體,波速恒定;為計算氣體空穴體積,將每一個控制體的體積平分為兩等份,兩個等分體內的壓力及流速均由FVM 求解.根據上述基本假定,管道內空穴瞬變流的數學模型由純水體和離散氣體空穴兩部分控制方程組成.

1.1 控制方程

1.1.1 水錘模型

當黏彈性管道中僅有液體時的連續方程和動量方程

式中,H為測壓管水頭,m;x為沿管道軸向距離,m;g為重力加速度,m/s2;V為斷面平均流速,m/s;t為計算時間,s;a為波速,m/s;J為管道摩阻引起的水頭損失;εr為管壁的滯后應變.

采用FVM 法求解,式(1)和式(2)瞬變流控制方程可以寫成帶源項的守恒性雙曲方程形式[15]

式中,U為待求的未知向量;F為通量向量;A為系數矩陣;S為源項,包含了管道的黏彈性和摩阻項.

1.1.2 離散氣體空穴模型

DGCM 將自由氣體的質量集中在控制體中心,可模擬自由氣體在均質混合流的分布.根據理想氣體方程,控制體中小體積氣穴隨著壓力的變化而進行等溫地[16]膨脹和收縮.氣穴變化可表示為

式中,Qu為計算步時內平均流入截面的流量;Q為計算步時內平均流出截面的流量;?g為氣體體積.

采用FVM 法Godunov 格式計算兩個等分單元內測壓管水頭及流量,再根據式(7)和式(8)計算空穴處測壓管水頭和氣穴體積

式中,ρ為水體的密度,kg/m3;z為計算節點高程,m;為第j個控制單元空穴上游側和下游側在計算時刻的測壓管水頭;為第j個控制單元平均測壓管水頭;Hv為液體的蒸汽壓頭.

當兩個等分單元內壓強均大于水體汽化壓強時,氣體空穴處測壓管水頭取兩部分均值,按式(7)和式(8)計算.此時離散氣體空穴會影響等分單元內的流量變化,引入壓力修正系數C_ap來考慮自由氣體的影響并作出調整[17].當任一等分單元內壓強降低到或低于汽化壓強,則設定控制體保持液體汽化壓強不變,此時,兩個等分單元內壓力均等于氣體空穴內壓力.當計算空穴體積小于0 時,空穴潰滅,繼續純液相水錘計算.

1.1.3 黏彈性效應和管道摩阻

與彈性材料不同,黏彈性管材在受力時不遵循胡克定律,產生的應變可由波爾茲疊加原理[18]描述:對于較小的應變,總應變等于相互獨立的一系列應力作用下產生應變的線性相加.式(1)中遲滯應變率?εr/?t表示水力瞬變過程中管壁黏彈性形變隨時間變化,用廣義Kelvin-Voigt (K-V)模型表達,該模型將彈簧和阻尼并聯看作一組元件,而管壁的應變可等效為若干組元件串聯收到加載后的形變變化,廣義K-V 模型示意圖如圖1 所示.

圖1 廣義K-V 模型示意圖Fig.1 Schematic diagram of generalized K-V model

蠕變函數表達式及數值求解[19]

式中,Jk為第k個元件中彈簧的蠕變柔量,m2/s;Tk為第k個元件中阻尼的松弛時間,s.

式中,D為管道直徑,m.

式(2)中的管道摩阻J由穩態摩阻Js和動態摩阻Ju兩部分組成,其中穩態摩阻為恒定狀態計算值,動態摩阻通過計算流動瞬時加速度和加權函數W的卷積得到[20].

穩態摩阻

動態摩阻

式中,ν 為運動黏滯系數,m2/s;u為在0～t之間積分的時間變量;W為關于無量綱時間 τ 的函數.

式(13)離散[21]成

式中,yi為歷史速度對當前時刻壁面剪切應力的影響;Δτ=4νΔt/D2;mi,ni為試驗擬合得到的系數.

1.2 FVM 二階Godunov 格式求解

1.2.1 數值通量計算及黎曼求解

采用有限體積法,將水錘計算區域離散為N個控制體,計算單元長度為Δx,計算時間間隔為Δt.計算區域局部離散網格系統示意圖如圖2 所示.

圖2 離散網格系統Fig.2 Discrete grid system

從圖中可以看出,任一計算單元i邊界界面分別為i-1/2 和i+1/2.變量(H,V)取控制體平均值,單元體的壓力、流速將通過界面i-1/2 和i+1/2 處數值通量計算,即

式中,n表示t時刻;n+1 表示t+Δt時刻;Fi-1/2,Fi+1/2分別表示界面i-1/2,i+1/2 的數值通量,在二階Godunov 格式[22-23]中根據計算單元i界面處的局部黎曼問題求解,即

本文采用MUSCL-Hancock 法[24]進行線性重構得到二階精度Godunov 格式,為避免線性插值中數據重構時出現虛假震蕩現象,引入斜率限制器MINMOD 函數[25]

黎曼解以控制體中間區域的待求變量值來近似代替精確解

1.2.2 邊界構建方法

邊界處理影響數值模擬的計算精度及效率,本文采用虛擬單元法構建邊界,基本思想與文獻[26]相似,將邊界信息直接通過計算區域體內虛擬點加以反映,可直接與已有的流體動力學求解器相結合.虛擬單元法關鍵在于重構虛擬點的信息,使邊界信息刻畫更加準確[27],最常見的處理方式是虛擬點沿著邊界法向鏡像外插[28].由二階Godunov 格式可知,計算控制體變量時需與之相鄰上下游各兩個單元數值,通過在兩端邊界各增加兩個虛擬控制體-1,0 和N+1,N+2,可以對邊界控制體和內部控制體進行統一計算模擬,并達到二階精度要求.其中虛擬單元0 和N+1 滿足鏡像外插,其數據與單元1 和N相等.

根據黎曼不變量理論[29]在計算域邊界處的解析解,可推得與正負特征線相關的黎曼不變量方程.對于上游水庫邊界,滿足虛擬單元U-1=U0=U1/2,代入=HR可得

對于下游閥門邊界,虛擬單元的數據需滿足UN+1=UN+2=UN+1/2

式中,CV=(V0τ)2/(2ΔH0),V0和ΔH0分別為初始時刻閥門處的流量和水頭損失,τ為閥門開度

1.2.3 時間積分與穩定性約束

將源項引入非線性系統求解,考慮黏彈性項和動態摩阻項,求解精度取決于對源項的積分方式.對于二階求解精度,采用二階龍格-庫塔法進行時間離散,其顯式求解過程如下

穩定性約束條件主要包括對流項的CFL(Courant-Friedrich-Lewy)條件,以及二階龍格-庫塔離散的約束,具體時間步控制方法如下

由此可以得到

因此,計算所允許的最大時間步長需滿足條件Δt≤min(Δtmax,CFL,Δtmax,s).

2 算例驗證

為驗證所建數學模型,以里斯本大學Soares AK 等[10]在葡萄牙搭建的試驗臺上測得的試驗數據為依據,將本文研究的黏彈性管道中有限體積法求解純水錘模型(FVM VE)及離散氣體空穴模型(FVM-DGCM VE)與常用的MOC 算法結果、傳統彈性管道模型計算結果、試驗數據進行對比分析.本文引用Soares 等[10]搭建的試驗裝置圖如圖3所示.

圖3 黏彈性管道試驗裝置圖[10]Fig.3 Diagram of viscoelastic pipeline test device[10]

試驗系統由水泵、空氣罐、上下游閥門、壓力傳感器和HDPE 管道組成.管道總長度為203.00 m,內徑為44.0 mm.水由水泵加壓流入空氣罐,保持上游壓力恒定,流經上游球閥后進入水平放置的管道系統,由回水管流出.選擇位于上游閥門、管道中點、下游閥門處的3 個測壓點試驗結果供參考校核.為分析系統中的壓力瞬變,分別進行快速關閉下游閥門試驗對純水錘驗證分析、快速關閉上游閥門試驗對水柱分離彌合水錘驗證分析.

2.1 純水錘驗證分析

工況1 為純水體水力瞬變試驗[10],系統中上游空氣罐水頭為35.6 m,流量為2.72 L/s,波速為315.0 m/s,雷諾數為80 000.計算時Cr=1.0,N=32.管道材料的蠕變函數常通過蠕變測試確定[19],但當黏彈性材料集成在管道系統中時,蠕變還取決于管道的軸向和環向約束,因此需基于實測瞬態壓力數據進行校核.Soares 等[11]提出的系統黏彈性參數基于遺傳算法和Levenberg 算法兩步逆模型確定.經校核發現當取3 組Kelvin-Voigt 單元可以較好擬合蠕變曲線,各元件松弛時間及蠕變柔量的取值為τ1=0.018 s,J1=0.256 GPa-1;τ2=0.50 s,J2=0.256 GPa-1;τ3=3.0 s,J3=0.256 GPa-1.為驗證模型的有效性,選擇工況1 的試驗數據與本文模型計算值、MOC 模擬結果進行對比,結果如圖4 所示.

圖4 工況1 下游閥門處計算值與試驗結果Fig.4 Calculation value and test result at downstream valve in Case 1

采用工況1 研究不同Cr的影響,圖5 分別給出閥門處MOC VE 及FVM VE 模型壓力計算的結果.

圖5 庫朗數對計算結果的影響Fig.5 Influence of Cr on the calculation result

由圖4 可知,工況1 數值模擬的最大峰值比試驗測得數據略高,可能由試驗關閥時間比數值模擬采用的時間稍大造成.FVM VE 模型與MOC VE 模擬的閥前水錘結果基本一致,總體上與試驗壓力波動曲線在極值、相位上能較好地吻合,根據圖4 結果可認為,本文模型能準確模擬試驗壓力波動整個過程.

圖5 中,計算中庫朗數Cr分別取1.0,0.5,0.1.圖5(a)結果表明:在MOC 求解格式下,隨著Cr的減小,計算壓力波動結果出現了明顯的數值耗散,這是由于傳統MOC 計算方法只有一階計算精度.圖5(b)結果顯示:隨著Cr減小,壓力波動只產生了輕微的數值衰減的現象,有限體積法二階Godunov求解格式具有二階精度,本文模型在Cr小于1.0 時仍然能保證數值模擬結果的穩定性和準確性.

在輸水管道中管道摩阻和黏彈性效應均可引起瞬變壓力衰減,考慮動態摩阻UF 的有限體積法黏彈性模型FVM VE+UF 與考慮穩態摩阻的FVM VE 模型計算結果如圖6 所示.

圖6 動態摩阻對計算壓力的影響Fig.6 Influence of dynamic friction on calculated pressure

由圖6 可知,在黏彈性輸水管道瞬態壓力計算中,動態摩阻的影響不再明顯.這是因為實際情況中影響管道因素還有流體的軸向運動、管道錨固約束、流體慣性效應等,因此校核得到的蠕變函數不一定僅代表HDPE 管壁準確黏彈性效應[30],而是包含動態摩阻等阻尼效應的組合.在瞬變流過程中,由于黏彈性效應造成的管道壓力衰減太快[31-32],瞬變流持續時間短,因此關閥完成后管壁摩擦做功時間短,做功較少[33].綜上所述,在黏彈性管道水力瞬變壓力衰減過程中,黏彈性效應相比于非恒定摩阻起主要作用,在后續計算中將不再考慮動態摩阻的影響.

2.2 水柱分離彌合水錘驗證分析

工況2,3 為水柱分離彌合試驗[10],系統中上游空氣罐水頭為30.6 m,波速287.0 m/s.計算時Cr=1.0,N=32,松弛時間及蠕變柔量取值τ1=0.10 s,J1=0.60 GPa-1;τ2=0.50 s,J2=0.35 GPa-1;τ3=3.0 s,J3=0.50 GPa-1.其中工況2 流量為3.00 L/s,雷諾數為87 000;工況3 流量為4.00 L/s,雷諾數為1.2×105.選擇工況2 驗證水柱分離彌合水錘現象.試驗初始流態沒有出現分布式氣泡,故采用建議較小的初始含氣率.

已有學者[17,34]討論采用MOC-DGCM 模型模擬時初始含氣率α0的變化對計算壓力產生的影響,且表明較小初始含氣率α0≤10-7時能較好地模擬試驗結果.本文FVM-DGCM VE 模型基于上述理論,對比α0分別為10-7,10-8,10-9時閥門處計算壓力,如圖7 所示,結果表明:當初始含氣率低于或等于10-7時,計算結果基本一致,因此后續計算中,α0均取為10-7.

圖7 FVM-DGCM VE 中α0 對上游閥門壓力的影響Fig.7 Effects of α0 in FVM-DGCM VE on pressure at upstream valve

圖8 是為關閉上游閥門時閥門處試驗測得的壓力-時間曲線與數值計算結果對比,模擬考慮壓力降低到蒸汽壓時汽化發生的水柱分離現象.

如圖8 所示,FVM-DGCM VE 模型模擬的閥后水錘與傳統算法MOC-DGCM VE 計算結果在每個周期內極值接近,隨著瞬變流的進行,MOC-DGCM VE 模型中的相位差不斷累積,使得模擬結果相位提前于試驗壓力曲線現象逐漸明顯,綜上所述,FVMDGCM VE 能更為準確地模擬出試驗的壓力峰值、波動周期.

圖8 工況2 上游閥門處計算值與試驗結果Fig.8 Calculation value and test result at upstream valve in Case 2

圖9 中三條曲線分別是試驗實測曲線、考慮管道黏彈性以及僅考慮管道彈性時的水柱分離彌合現象模擬結果.

圖9 工況3 上游閥門處計算值與試驗結果Fig.9 Calculation value and test result at the upstream valve in Case 3

由圖9 可以看出,在相位變化上,考慮管道黏彈性的水柱分離模型模擬結果與試驗結果吻合較好,而彈性管道模型在后續周期中與試驗結果相差較大.在壓力波峰值上FVM-DGCM VE 計算得到的相對誤差值為3.34%,模擬結果更接近試驗結果,而彈性管道模型結果相對誤差值為25.15%,振幅值偏大.從壓力波振幅值和相位變化綜合來看,考慮管道黏彈性的水柱分離彌合瞬變流模型的模擬結果更接近實測數據.

3 結論

本文基于有限體積法二階Godunov 求解格式,針對水柱分離彌合水錘現象,建立了考慮管道黏彈性的離散氣體空穴模型進行研究分析.將本文模型計算結果與試驗結果以及彈性管道瞬變流模型計算值進行比較,分析了動態摩阻、庫朗數及初始含氣率等因素的影響,主要得出以下結論.

(1)本文建立的考慮黏彈性效應的有限體積法求解水柱分離彌合水錘模型能夠準確模擬出純水錘和水柱分離彌合水錘兩種情況的瞬變壓力,均與試驗數據高度吻合.

(2)與已有的特征線求解模型相比,本文建立的有限體積法求解水柱分離彌合水錘模型更具穩定性和準確性.當庫朗數Cr小于1 時,特征線模型會產生嚴重的數值耗散現象,且MOC-DGCM VE 模型模擬結果相位提前于試驗壓力曲線,本文模型能更準確地模擬試驗壓力波動曲線.

(3)在黏彈性管道水力瞬變壓力衰減過程中,黏彈性效應相比于非恒定摩阻起主要作用.在水柱分離彌合數值模擬過程中,當初始含氣率低于或等于10-7時,計算結果基本保持一致,因此對于較小的初始含氣率,α0可取為10-7.

(4)對于黏彈性管道中的瞬變流現象,考慮管道黏彈性可顯著提高模擬準確性.黏彈性管道瞬變流模型模擬結果能較好地與實測數據相吻合,而彈性管道瞬變流模型模擬的結果在壓力波的振幅、相位變化等方面都與試驗結果有較大的差異.