盧坤銘，周領(lǐng)，劉靜

(河海大學水利水電學院，江蘇 南京 210098)

實際工程中，輸水管道的高程會根據(jù)地形條件發(fā)生變化，往往會出現(xiàn)氣團滯留于管道中來不及排出.水流沖擊滯留氣團現(xiàn)象，不但會影響系統(tǒng)正常運行，甚至可能引發(fā)爆管事故，從而造成極其嚴重的經(jīng)濟損失，甚至導致人員傷亡[1-3].

目前，針對水流沖擊滯留氣團瞬變流現(xiàn)象，現(xiàn)有的模型大多為一維數(shù)學模型，且主要研究對象為單個氣團，對于起伏管道內(nèi)水流沖擊多段氣團的瞬變流研究，相關(guān)成果很少.MARTIN[4]首次對水流沖擊多段滯留氣團現(xiàn)象進行了理論研究，基于剛性水體理論，建立了相應(yīng)的數(shù)學模型，但其模型沒有考慮水氣交界面的運動.劉德有等[5]建立了起伏變特性管道系統(tǒng)中水流沖擊單個滯留氣團的剛性數(shù)學模型.剛性模型具有簡單、快捷的優(yōu)點，但其應(yīng)用具有一定的局限性，如：在復雜的管網(wǎng)系統(tǒng)中，剛性模型處理起來尤為煩瑣和復雜，尤其是出現(xiàn)液柱分離的情況.EPSTEIN[6]將積分法用于水錘基本方程的求解，建立水流沖擊兩段封閉滯留氣團的數(shù)學模型，由于沒有相應(yīng)的多段氣團的彈性數(shù)學模型，該方法僅與剛性數(shù)學模型進行了比較，結(jié)果顯示當氣團長度很小時，剛性模型預(yù)測的氣團壓力過高.ZHOU等[7-9]研究了水流沖擊滯留氣團的瞬變壓力及參數(shù)變化規(guī)律，發(fā)現(xiàn)水氣耦合作用機理復雜，且可能引起10倍于入口壓力的危險峰值壓力；同時，提出了“剛性塞”、“虛擬塞”等方法，建立了彈性水體模型，成功避免了特征線法在水氣交界面動態(tài)追蹤時的復雜插值處理，并保持了同樣的計算精度；指出在氣團含量很小的情況下，剛性水體模型的不適用性.

現(xiàn)有的一維數(shù)學模型存在諸多假定，其計算精度往往不能滿足實際需求.對于水氣耦合兩相瞬變流而言，整個過程往往很難通過試驗手段準確記錄.文中采用三維CFD方法模擬起伏管道內(nèi)含多段滯留氣團的瞬變流，應(yīng)用Standardk-ε湍流模型，研究兩段滯留氣團的相互作用規(guī)律，并將三維計算結(jié)果、一維模型計算結(jié)果與試驗結(jié)果進行對比分析；同時，研究整個過程中水-氣作用規(guī)律.

1 試驗裝置

整個系統(tǒng)從上游至下游依次由蓄水池、潛水泵、氣罐(壓力罐)、電磁流量計、球閥、進氣孔口、泄水閥、封閉末端組成.水泵與氣罐之間通過不銹鋼鋼管連接，壓力罐至下游由一段1 m長的不銹鋼和多段起伏的有機玻璃透明管道組成.從氣罐出口至管道末端為水流沖擊滯留氣團的試驗研究管道，其總長為10.97 m，其中，有機玻璃管道內(nèi)徑4 cm，壁厚1 cm，具體如圖1所示.

圖1 含多個滯留氣團試驗裝置布置圖

試驗中，將氣罐出口和管道入口交界面處定為x= 0，水平管道中心線定為z= 0；圖1中給出了管道彎曲段最高點和最低點處的沿線長度x和高程位置z.球閥距上游入口距離為x= 2.236 m；共8個壓力傳感器，其安裝位置分別如下：PT-1#為(x=10.82 m,z=0.632 m)，PT-2#為(x=9.15 m,z=0.075 m)，PT-3#為(x=7.60 m,z=0.82 m)，PT-4#為(x=6.68 m,z=0.24 m)，PT-5#為(x=5.75 m,z=0.32 m)，PT-6#為(x=4.33 m,z=0.96 m)，PT-7#為(x=1.707 m,z=0 m)，PT-8#為(x=1.337 m,z=0 m)；在彎曲管道最頂部和最底部分別安裝3個進氣孔、3個排水閥，僅用于調(diào)節(jié)初始狀態(tài)下氣團段數(shù)和長度，但在瞬變過程中均處于關(guān)閉狀態(tài).各壓力傳感器性能一致，測量范圍：0～1 MPa，線性誤差和滯后誤差范圍為0.03%～0.30%(廠家提供).電磁流量計性能參數(shù)，公稱通徑：DN40，流量范圍：0～25 m3/h.數(shù)據(jù)采集系統(tǒng)為美國國家儀器(NI)有限公司產(chǎn)品：PCI-6221，采樣率250 kS/s.

試驗初始條件：上游球閥完全關(guān)閉，球閥上游為有壓水體與氣罐相通，球閥下游為起伏管段，初始氣團的個數(shù)和長度通過3個進氣孔、4個排水閥進行調(diào)節(jié)實現(xiàn)的，管道下游末端始終處于封閉狀態(tài).水流沖擊滯留氣團瞬變過程是通過瞬間開啟球閥實現(xiàn)的，整個瞬變過程中氣團始終滯留于管道內(nèi).

表1為3種試驗工況的初始狀態(tài)，各物理參數(shù)：xf為起伏管道系統(tǒng)中，閥門開啟后t時刻沖擊水體的長度，m；x1u為第1 段阻斷水體的上游端沿線位置，m；x1d為第1 段阻斷水體的下游端沿線位置，m；x2u為第2 段阻斷水體的上游端沿線位置，m；x2d為第2 段阻斷水體的下游端沿線位置，m；La0,1為第1段氣團的長度，m；La0,2為第2段氣團的長度，m；La0為氣團總長度.試驗針對含2個初始滯留氣團情況，進行了3種工況試驗入口壓力均為Pr= 0.16 MPa，初始沖擊水體長度均為xf= 3.3 m，具體初始狀態(tài)值見表1.為保證試驗精確性，嚴格、精確控制試驗的初始條件和邊界條件；每種工況重復至少3次，直至試驗結(jié)果基本一致.

表1 3種試驗工況的初始狀態(tài)值

為了便于數(shù)值模擬計算比較分析，試驗通過高速攝像機的拍攝記錄顯示，試驗中手動球閥的開啟時間(從全關(guān)至全開)為0.07～0.09 s，通過純水錘試驗得水錘波速約為400 m/s.

2 CFD數(shù)值計算

2.1 三維物理模型及網(wǎng)格劃分

通過Design Modeler對試驗系統(tǒng)進行三維物理建模，導入ICEM劃分網(wǎng)格，如圖2所示.為了獲得更加準確的計算結(jié)果，對彎管處、水箱下部、近壁面處進行加密處理，均采用六面體網(wǎng)格.整個管道的網(wǎng)格數(shù)為33.2萬，最小網(wǎng)格質(zhì)量為0.350，平均網(wǎng)格質(zhì)量為0.872.

圖2 三維模型及網(wǎng)格示意圖

2.2 網(wǎng)格無關(guān)性驗證

針對工況1，采用3種不同數(shù)量網(wǎng)格進行數(shù)值模擬計算并與試驗結(jié)果對比分析驗證，得出計算耗時與壓力誤差.通過改變網(wǎng)格節(jié)點數(shù)獲取不同網(wǎng)格數(shù)量，分別為26.2萬、33.2萬、53.8萬，對應(yīng)的壓力峰值與試驗結(jié)果的誤差分別為0.70%,0.45%和0.20%.基于模擬精度與計算時間的考慮，最終選用網(wǎng)格數(shù)量為33.2萬的網(wǎng)格.

2.3 基本控制方程及湍流模型

基本控制方程包括容積比率方程、動量方程、能量方程以及水體狀態(tài)方程.

容積比率方程為

(1)

動量方程為

(2)

能量方程為

(3)

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

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

(4)

(5)

其中：a為水錘波速；D和E為管徑和管壁厚度；E為彈性管道楊氏模量.

湍流模型：通過對3種湍流模型(Standardk-ε，RNGk-ε和Realizablek-ε)進行適用性分析[10]，Standardk-ε模型能更準確地模擬流態(tài)變化過程.故文中采用Standardk-ε湍流模型進行數(shù)值模擬研究.

2.4 求解控制及邊界條件

求解控制：瞬態(tài)計算中，采用有限體積法(FVM)對網(wǎng)格進行離散，時間步長為1.0×10-4s，所有參數(shù)收斂標準均為1.0×10-6，每一時間步最大迭代步數(shù)為100步，經(jīng)后期觀測發(fā)現(xiàn)迭代60次以內(nèi)均能收斂.速度、壓力和密度的耦合方式選擇PISO方法，壓力項采用PRESTO格式，動量項、能量項采用二階迎風格式，湍動能項、湍流耗散項采用一階迎風格式，體積分數(shù)采用Geo-Reconstruct格式.計算時長5 s.

邊界條件：壓力罐下部Patch為有壓水體，阻斷水體及氣團均為大氣壓.閥門區(qū)域采用動網(wǎng)格模型，開閥時間為0.08 s.上游管道出口與閥門進口、閥門出口與起伏管道進口設(shè)置為2對“interface”邊界條件，其他均為默認邊界條件.

3 計算結(jié)果

3.1 瞬變壓力三維CFD計算與試驗結(jié)果對比

圖3為不同氣團長度下，PT-1#和PT-4#壓力傳感器處壓力波動計算結(jié)果與試驗結(jié)果.表2給出了PT-1#和PT-4#壓力傳感器處最大壓力(Ha1-max,Ha2-max)的CFD值pCFD和試驗值pE.與管道內(nèi)含一段滯留氣團結(jié)果相比，管道中含兩段滯留氣團增加了水力瞬變過程的復雜性.試驗結(jié)果表明：① 系統(tǒng)的最大壓力始終出現(xiàn)在靠近管道末端的滯留氣團的活動區(qū)域；② 位于管道末端的第2段阻斷水體(僅工況3中)，始終保持靜止，內(nèi)部壓力處處相等，其與氣團壓力變化基本一致；③ 在水流沖擊過程中，由于兩段滯留氣團間相互作用，導致氣團的壓力變化曲線出現(xiàn)峰值波動，其發(fā)生位置和劇烈程度與初始兩段氣團相對體積大小有密切關(guān)系.

表2 氣團最大壓力值的CFD值和試驗值

圖3和表2表明，2個氣團的數(shù)值模擬結(jié)果與試驗數(shù)據(jù)十分吻合，三維數(shù)值模擬可以很好地研究水流沖擊多段滯留氣團的瞬態(tài)壓力.此外，水流沖擊兩段滯留氣團壓力波動曲線中，兩氣團的瞬變壓力并非同步變化，可能呈現(xiàn)兩氣團峰值壓力交替出現(xiàn)的情況，這與初始氣團長度有一定的關(guān)系.如圖3a所示，第2段滯留氣團的壓力峰值大于第1段滯留氣團，初始的兩段滯留氣團的壓力變化過程中，兩者的峰值壓力交替出現(xiàn).如圖3b所示，與工況1相比，工況2中初始第1段滯留氣團長度增大，其壓力變化變得平緩，而第2段滯留氣團長度減小，壓力變化較為劇烈，2種情況的壓力曲線的總體趨勢是一致的.圖3c顯示，與工況1相比，第1段氣團長度保持不變，第2段滯留氣團縮短91%，導致系統(tǒng)壓力峰值增加23%，壓力峰值響應(yīng)時間縮短58%.3種工況峰值壓力和初始氣團總體含量對比表明，同一入口壓力下，隨著初始氣團總體含量的增加，系統(tǒng)壓力峰值相應(yīng)減小.

圖3 2個氣團計算結(jié)果與試驗結(jié)果對比

3.2 現(xiàn)有一維模型計算結(jié)果對比

圖4為3種試驗工況下，采用Standardk-ε湍流模型三維模擬得出的系統(tǒng)最大壓力與一維模型計算結(jié)果、試驗觀測數(shù)據(jù)的對比結(jié)果.其中，一維模型采用本課題組已發(fā)表的彈性水體數(shù)學模型[10].

3種工況均表明三維CFD方法能夠準確地模擬整個壓力波動過程，其壓力峰值與試驗結(jié)果更接近.一維彈性水體模型可以簡單模擬整個瞬變過程的壓力變化趨勢，與試驗結(jié)果存在較大誤差，其壓力峰值均大于試驗值.一維數(shù)值模擬結(jié)果表明，工況2在第一個波動周期后存在明顯的壓力振蕩，說明第2段氣團受到阻斷水體的反復壓縮.彈性水體模型假定：① 沖擊水體和阻斷水體間被氣團完全隔開；② 滯留氣團占據(jù)整個管道截面，水氣交界面垂直于管道中心線，水、氣互不混摻；③ 管道內(nèi)滯留氣團假定為完善氣團，滿足氣體熱力學多變過程方程.試驗以及三維數(shù)值模擬過程中發(fā)現(xiàn)，水體運動至管道彎曲處，可能將氣團分成若干部分，能量衰減加快，阻斷水體長度發(fā)生變化，導致一維數(shù)學模型的計算結(jié)果與試驗有所偏差，無法準確模擬整個壓力過程.三維CFD結(jié)果明顯好于一維模型結(jié)果，但與試驗結(jié)果仍有差別，表現(xiàn)為從第二個周期開始提前于試驗結(jié)果，可能的原因是很難真實地模擬水-氣-管壁的熱傳遞.

圖4 一維-CFD-試驗結(jié)果對比

3.3 動態(tài)特性分析

圖5為工況1下水-氣兩相分布圖.

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

結(jié)合試驗觀察結(jié)果，水流沖擊管道內(nèi)含多段滯留氣團的瞬變流過程可描述如下：初始時，上游閥門完全關(guān)閉，閥門上游段為高壓段，與壓力罐相連，閥門下游段管道系統(tǒng)中滯留兩段氣團和阻斷水體.當上游閥門瞬間打開之后，管道系統(tǒng)開始充水，靠近上游的第一段滯留氣團在沖擊水體的作用下開始壓縮，其壓力逐漸增大；第一段滯留氣團壓力增大到一定值時，開始推動阻斷水體向下游運動.與此類似，第二段滯留氣團也受到壓縮，由于管道系統(tǒng)末端完全封閉，管道下部阻斷水體撞擊在管道末端，氣團壓力迅速增大并反向壓縮阻斷水體，隨著水氣不斷混摻，最終阻斷水體被分隔成兩段，滯留氣團則變?yōu)槿?

數(shù)值模擬結(jié)果發(fā)現(xiàn)，隨著氣團的壓縮-膨脹-變形，水氣交界面自由變化，并不垂直于管道中心線，同時，沖擊水體和阻斷水體互相摻混，當水體運動到管道彎曲處時會產(chǎn)生新的阻斷水體，其長度發(fā)生變化，將氣團分成若干部分.這與試驗結(jié)果保持一致，并進一步解釋三維CFD數(shù)值模擬結(jié)果優(yōu)于一維數(shù)學模型.

4 結(jié) 論

1) 三維CFD模型能夠準確地對起伏管道內(nèi)水流沖擊多段滯留氣團的瞬變流現(xiàn)象進行可視化分析與處理，且與試驗結(jié)果吻合較好.

2) 三維CFD模擬結(jié)果與試驗觀察結(jié)果顯示，系統(tǒng)的最大壓力始終出現(xiàn)靠近管道末端的滯留氣團的活動區(qū)域.在水流沖擊過程中，由于兩段滯留氣團間相互作用，導致氣團的壓力變化曲線出現(xiàn)峰值波動，其發(fā)生位置和劇烈程度與兩段氣團相對體積大小有密切關(guān)系.

3) 一維彈性模型計算結(jié)果誤差大于CFD結(jié)果，這是因為一維模型無法有效地模擬水氣交界面變化過程、阻斷水體長度的動態(tài)變化.

4) 對水流沖擊多段滯留氣團進行動態(tài)模擬并分析誤差.整個過程中，水氣兩相互相摻混，水氣交界面自由變化，阻斷水體長度時刻發(fā)生變化，當水體運動到管道彎曲處時會產(chǎn)生新的阻斷水體，將氣團分成若干部分.