盧坤銘，周領(lǐng)*，曹波，劉德有，王歡

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

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

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

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

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

1 試驗(yàn)研究

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

圖1 試驗(yàn)裝置布置圖

Fig.1 Test equipment arrangement diagram

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

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

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

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

2 數(shù)值模型

2.1 基本控制方程

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

容積比率方程為

(1)

動(dòng)量方程為

(2)

能量方程為

(3)

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

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

(4)

2.2 湍流模型

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

2.3 離散格式

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

2.4 幾何模型及邊界條件

圖2為3D模擬計(jì)算域及網(wǎng)格分布圖.計(jì)算域中所有網(wǎng)格采用六面體結(jié)構(gòu)化網(wǎng)格.

圖2 計(jì)算域及網(wǎng)格分布

3 結(jié)果與分析

3.1 網(wǎng)格無關(guān)性分析

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

表3 網(wǎng)格信息

圖3 壓力計(jì)算與試驗(yàn)結(jié)果

3.2 湍流模型的影響

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

圖4 不同湍流模型下，PT1#處壓力計(jì)算結(jié)果與試驗(yàn)結(jié)果對比

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

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

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

結(jié)果顯示，與一維模型相比，三維CFD計(jì)算瞬變壓力結(jié)果與試驗(yàn)結(jié)果更為接近.

圖5 計(jì)算結(jié)果與試驗(yàn)結(jié)果對比

3.4 三維模擬誤差分析

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

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

圖6為工況2(已證實(shí)該工況下，三維計(jì)算結(jié)果的準(zhǔn)確性)下水-氣兩相分布變化云圖.

圖6 不同時(shí)刻水-氣動(dòng)態(tài)特性

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

4) 與一維模型相比，三維CFD計(jì)算結(jié)果與試驗(yàn)結(jié)果更為接近.一維模型存在較大計(jì)算誤差，是因?yàn)槠浔旧砑俣ㄔ斐傻模@些假定均與實(shí)際情況有較大差別.

5) 對于初始?xì)鈭F(tuán)較短的情況，三維CFD計(jì)算結(jié)果明顯好于一維模型結(jié)果，但與試驗(yàn)結(jié)果仍有差別.原因可能是，隨著氣團(tuán)初始長度減小，水氣耦合、水氣摻混的程度更加劇烈，涉及復(fù)雜的熱力學(xué)問題.