盧坤銘，周領*，曹波，劉德有，王歡

(1. 河海大學水利水電學院，江蘇 南京 210098； 2. 國網新源水電有限公司新安江水力發(fā)電廠， 浙江 建德 311600)

在實際輸水系統(tǒng)中，管道布置常因地理條件限制及工程需要而呈現(xiàn)高低起伏狀態(tài).在這類管道系統(tǒng)中，氣團極易滯留于管道頂部，當系統(tǒng)啟動(閥門、水泵開啟)時，有壓水流將沖擊滯留氣團而產生瞬變流，從而易產生異常的危險壓力，極有可能造成管道爆裂，危害整個管道系統(tǒng)的安全運行[1-2].

目前，針對水流沖擊滯留氣團瞬變流現(xiàn)象，現(xiàn)有的模型大多為一維數(shù)學模型.作為該問題研究的先驅者，MARTIN[3]基于剛性水體理論，首次對水流沖擊滯留氣團的現(xiàn)象進行了研究，指出滯留氣團的存在可能引起異常的壓力波動.ZHOU等[4]對MARTIN的剛性模型進行了改進，充分考慮了沖擊水體長度的動態(tài)變化.CHAIKO等[5]考慮水體和管道的彈性，建立了3種彈性數(shù)學模型，并給出了各彈性模型的適用范圍.一些學者先后提出了一維“虛擬塞”“剛性塞”數(shù)學模型，與考慮全水體彈性數(shù)學模型相比，成功地避免了特征線法在水氣交界面動態(tài)追蹤時的復雜插值處理，并保持了同樣的計算精度；指出在氣團含量很小的情況下，剛性水體模型并不適用[6-8].隨后，ZHOU等[9]通過試驗和一維數(shù)值模擬，研究了2段滯留氣團和阻斷水體之間的關系，結果表明：與單個氣團相比，2個滯留氣團可能導致更加危險的峰值壓力，這與2個氣團的大小及阻斷水體長度密切相關.然而，有些學者的研究成果表明現(xiàn)有一維數(shù)學模型常常低估了氣團壓力峰值，這與一維模型的基本假定有關，包括：① 滯留氣團占據整個管道截面，水氣交界面與管中心線垂直，氣、水互不摻混；② 在穩(wěn)態(tài)和瞬態(tài)情況下，管道內的水流阻力特性不變；③ 管道內滯留氣團假定為完善氣體，滿足氣體熱力學多變過程方程[4,6-9].但試驗觀測發(fā)現(xiàn)，水-氣耦合過程呈現(xiàn)明顯的三維特征，整個瞬變過程中水-氣交界面自由變化，會出現(xiàn)氣團分離成2個或若干個氣團的情況，能量衰減加快，這將導致數(shù)學模型的計算結果與試驗有所偏差.

為了能精細模擬水氣耦合過程，近年來，論文作者研究團隊成功將計算流體動力學(CFD)方法引入對簡單管道內含滯留氣團瞬變流進行了二維、三維數(shù)值模擬，并通過試驗數(shù)據驗證CFD方法的可行性和準確性[10-11].水平管道情況模擬結果顯示，不同湍流模型對滯留氣團瞬態(tài)特性和壓力變化有較大影響[12]；同時，通過三維模擬手段揭示了豎直管道內水流沖擊滯留氣團瞬態(tài)過程中滯留氣團的熱力學特性.此外，MARTINS等[13]采用同樣的三維CFD方法研究了管道填充過程中滯留氣團的運動特性，指出三維模擬能更好地揭示水氣動態(tài)變化特性.

文中將采用三維CFD方法模擬起伏管道內水流沖擊滯留氣團的瞬變流.現(xiàn)有的一維模型存在諸多假定，且計算結果與試驗差別較大，很難采用一維模型來對該瞬變流進行深入的研究.對于水-氣耦合三維瞬變流，試驗手段很難準確、全面記錄水-氣的三維動態(tài)特性.三維CFD方法已被證實可用于簡單管道瞬變問題的研究，然而，起伏管道的情況更為復雜，需要進一步研究.文中將對比分析3種湍流模型(Standardk-ε、RNGk-ε和Realizablek-ε)對結果的影響，并將三維計算結果、一維計算結果和試驗結果進行對比；同時，研究瞬變過程中氣團壓力波動和氣水形態(tài)變化.

1 試驗研究

針對起伏管道內水流沖擊滯留氣團的瞬變流現(xiàn)象，設計搭建了如圖1所示的試驗系統(tǒng)，分析多個滯留氣團的相互作用規(guī)律[7].整個系統(tǒng)主要由潛水泵、氣壓罐、電磁流量計、球閥、管道、空氣閥、排水閥等組成.整個管道總長10.97 m，由內徑D為4 cm的透明有機玻璃管道(長9.883 m，壁厚e為1 cm)及鋼管(長1.087 m，壁厚0.5 cm)組成.起伏段采用透明的有機玻璃，以便于觀察水-氣交界面的變化及水-氣兩相分布規(guī)律.進氣孔直徑3.2 mm，排水閥直徑4.2 mm.圖1中給出了管道彎曲段最高點和最低點處的沿線長度x和高程位置z.管道上安裝了1個球閥及8個壓力傳感器(PT1, PT2, PT3, PT4, PT5, PT6, PT7, PT8)，其安裝位置見表1.

圖1 試驗裝置布置圖

Fig.1 Test equipment arrangement diagram

氣壓罐直徑0.75 m，高2.075 m，有效容積750 L，可以提供0～1 MPa的基本恒定的水壓力.電磁流量計用于測量水流量，其流量范圍為0～25 m3/h；壓力表的最大量程為1 MPa，用于讀取壓力罐的壓力值；壓力傳感器的測量范圍為0～5 MPa，以測量不同位置處的壓力值；數(shù)據采集系統(tǒng)將相應的電壓信號轉換為壓力和流量值，文中采集頻率為1 000 Hz.

表1 球閥及壓力傳感器位置

文中將針對起伏管道內初始單個滯留氣團的瞬變過程進行三維模擬，模擬工況參數(shù)見表2，表中pr，xf，xu，La0分別為入口壓力、初始沖擊水體下游斷面位置、阻斷水體上游斷面位置、氣團長度.初始時所有閥門關閉，球閥快速開啟以實現(xiàn)水流沖擊氣團瞬變過程，球閥開啟時間0.1 s.

表2 模擬工況參數(shù)

2 數(shù)值模型

2.1 基本控制方程

采用VOF模型[14]追蹤水氣交界面，三維模型的基本控制方程[10, 15]如下.

容積比率方程為

(1)

動量方程為

(2)

能量方程為

(3)

式中:ρw為水體密度；αw為水體體積分數(shù)；t為時間；v為流體速度；ρ為水-氣混合物密度；p為靜壓；μ為水-氣混合物黏度；μw，μa分別為水、氣相流體黏度；g為重力加速度；F為外部體積力；E為總能；T為溫度；keff為有效傳熱系數(shù).

模擬中，氣體、水體均可壓.假定氣體遵循理想氣體定律，水體狀態(tài)方程為

(4)

2.2 湍流模型

文中對比了Standardk-ε、RNGk-ε、Realizablek-ε這3種湍流模型計算結果.3種模型都適用于湍流，分子之間的黏性可以忽略，區(qū)別在于計算湍流黏性的方法不同，`以及湍流的產生和擴散的因素不同.

2.3 離散格式

采用的離散方法為有限體積法(FVM)，速度、壓力和密度的耦合方法選擇PISO算法，壓力項采用PRESTO格式，動量項、能量項采用二階迎風格式，湍動能項、湍流耗散項采用一階迎風格式，體積分數(shù)采用Geo-Reconstruct格式進行離散，時間差分采用一階隱式格式，時間步長為0.000 1 s，計算時長為5 s.

2.4 幾何模型及邊界條件

圖2為3D模擬計算域及網格分布圖.計算域中所有網格采用六面體結構化網格.

圖2 計算域及網格分布

3 結果與分析

3.1 網格無關性分析

采用3套網格針對工況1進行對比分析.表3為網格信息，表中n，T，ε分別為網格數(shù)量、計算時間、最大壓力誤差.圖3為壓力計算與試驗結果.隨著網格數(shù)的增加，計算耗時顯著延長，系統(tǒng)最大壓力的精度有所提升，而3種網格下計算的壓力變化曲線的區(qū)別并不明顯.考慮到計算效率，最終采用第1套網格來進行相關模擬分析.

表3 網格信息

圖3 壓力計算與試驗結果

3.2 湍流模型的影響

圖4為工況1及工況2時，不同湍流模型下，壓力波動計算結果與試驗的對比結果.可以看出：① 對工況1(pr=0.08 MPa)，壓力幅值大小、波動周期都與試驗數(shù)據十分吻合，3種湍流模型都有效地模擬了氣團壓力變化過程.對于第1壓力峰值，Standardk-ε湍流模型的結果略高于試驗值，誤差為0.7%，RNGk-ε湍流模型的結果略低于試驗值，誤差為-0.2%；Realizablek-ε模型與試驗結果的誤差最大，為4.8%.在第2壓力峰值處，RNGk-ε,Realizablek-ε這2種湍流模型與試驗結果的誤差較大，其原因為在較低入口壓力下，在近壁面區(qū)湍流的發(fā)展并不充分，使考慮了高應變率和旋流的修正模型產生了較大的偏差.對時間響應，Standardk-ε湍流模型的時間響應與試驗結果幾乎同步，而2種修正模型則略早于試驗結果.② 對工況2(pr=0.16 MPa)，對于第1壓力峰值，Standardk-ε湍流模型的結果與試驗非常吻合，誤差僅為0.1%.而另外2種湍流模型的誤差則較大，分別為4.7%和9.0%，可能的原因是在較高壓入口壓力下，整個系統(tǒng)湍流充分發(fā)展，能量損耗主要以水體與管壁的熱傳導，以及水-氣的熱交換為主，修正模型的湍動黏度計算則考慮了流動中的旋流和曲率的影響.對于時間響應，該工況與pr=0.08 MPa工況存在相似的規(guī)律.此外，工況3,4可得類似結論，由于篇幅限制，在此不再贅述.

圖4 不同湍流模型下，PT1#處壓力計算結果與試驗結果對比

Fig.4 Comparisons between calculated and experimental results of pressure at PT1# with different turbulence models

3.3 現(xiàn)有一維模型計算結果

圖5為4種試驗工況下，采用Standardk-ε湍流模型三維模擬得出的壓力波動曲線.與現(xiàn)有一維模型計算結果、試驗觀測數(shù)據的對比結果.其中，一維模型是采用本課題組已發(fā)表的彈性水體數(shù)學模型[10]；管道恒定摩阻取f=0.1.

結果顯示，與一維模型相比，三維CFD計算瞬變壓力結果與試驗結果更為接近.

圖5 計算結果與試驗結果對比

3.4 三維模擬誤差分析

工況1和工況2時，三維CFD方法能較為準確地預測出瞬變壓力，如圖5所示.然而，對于初始氣團較短的工況3，4，雖然三維CFD計算結果明顯好于一維模型結果，但與試驗結果仍有差別.可能的原因是，隨著氣團初始長度減小，水氣耦合變得更加劇烈，如工況4的壓力峰值和波動頻率遠高于工況2，這意味著水氣摻混的程度更加劇烈，將涉及復雜的熱力學問題.另外，試驗發(fā)現(xiàn)，工況4中，壓力峰值較大，出現(xiàn)輕微的管道振動.然而，文中三維模擬沒有考慮水-氣-管壁復雜的熱傳遞和熱交換問題，以及流固耦合現(xiàn)象.另外，數(shù)值模擬很難準確設定和考慮實際管道材料、流體特性等.

3.5 水-氣動態(tài)特性

圖6為工況2(已證實該工況下，三維計算結果的準確性)下水-氣兩相分布變化云圖.

圖6 不同時刻水-氣動態(tài)特性

結合圖5b，起伏管道內水流沖擊滯留氣團過程中，水-氣三維動態(tài)特性如下：

1)t=0 s, 初始時刻，上游管道、閥門處都充滿了有壓水體，沖擊水體至起伏管道末端為氣體，壓強為大氣壓力，水-氣交界面位于6#管道處，且?guī)缀醮怪庇诠艿乐行木€.

2)t=0.30 s，迅速打開閥門后，在有壓水流的快速沖擊下，滯留氣團首次被壓縮，其壓力迅速增大，水-氣交界面在1.5 s內沿管道迅速運動至4#和3#管道的彎曲處.水體和滯留氣團的交界面十分明顯，但由于重力的作用，很明顯沒有垂直于管道中心線.

3)t為0.30～0.60 s時，沖擊水體從6#管道運動至6#和5#管道的彎曲處，在重力的作用下，管道下方的水體運動速度更快，比管道上方的水體提前到達5#和4#管道的彎曲處，氣團由于密度小被沖擊水體擠壓在管道的上方，水-氣交界面擴大為“楔形”，與管中心線大致呈-45°夾角，面積逐漸增大.

4)t=0.95 s，在高速水流的沖擊作用下，5#和4#管道的彎曲處填充了大量水體，滯留氣團繼續(xù)被壓縮，但存在少部分的氣團被分離出來，水-氣互相摻混，部分氣體溶解在水體中，形成氣泡流存在管道彎曲處.在滯留氣團的阻力影響以及重力作用下，管道下方的水體運動速度減慢，但仍大于管道上方的水體運動速度，導致水-氣交界面繼續(xù)以“楔形”存在，與管道中心線的夾角為45°.相比于t=0.6 s，滯留氣團的壓力繼續(xù)增大，對水體的緩沖作用逐漸增強，導致水-氣交界面的面積明顯減小.

5)t=1.60 s，氣團開始膨脹，擴大了水-氣交界面的面積，由于慣性和重力作用，管道下方部分水體沿著3#管道運動至3#和2#管道的彎曲處，同時氣團壓力開始減小.

6)t=2.15 s，隨著氣團壓力達到最小值，滯留氣團的體積膨脹至最大，水體回流至4#管道，水-氣交界面垂直于管道中心線.

4 結 論

1) 針對起伏管道內水流沖擊滯留氣團的瞬變流現(xiàn)象，所采用的三維CFD模型能夠有效模擬瞬變過程中壓力波動、氣-水界面變化，且與試驗結果吻合較好.三維CFD方法可以對整個流場實現(xiàn)可視化處理，這是試驗及一維模型都難以達到的.

2) 湍流模型是三維CFD建模的核心部分之一，影響管路系統(tǒng)的阻力、能量衰減的核心因素.對比發(fā)現(xiàn)，Standardk-ε湍流模型能夠更好地模擬起伏管道內水流沖擊滯留氣團瞬變過程的壓力峰值和波動周期.

3) 三維CFD模擬結果與試驗觀察結果顯示，隨著氣團的壓縮和膨脹，水-氣兩相互相摻混，交界面并不垂直于管道中心線，而是與管道中心線呈一定的夾角，水-氣交界面自由變化，以及水氣摻混等復雜兩相流態(tài).

4) 與一維模型相比，三維CFD計算結果與試驗結果更為接近.一維模型存在較大計算誤差，是因為其本身假定造成的，這些假定均與實際情況有較大差別.

5) 對于初始氣團較短的情況，三維CFD計算結果明顯好于一維模型結果，但與試驗結果仍有差別.原因可能是，隨著氣團初始長度減小，水氣耦合、水氣摻混的程度更加劇烈，涉及復雜的熱力學問題.