水流沖擊多段滯留氣團(tuán)的三維數(shù)值模擬

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

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

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

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

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

1 試驗(yàn)裝置

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

圖1 含多個(gè)滯留氣團(tuán)試驗(yàn)裝置布置圖

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

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

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

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

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

2 CFD數(shù)值計(jì)算

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

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

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

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

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

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

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

容積比率方程為

(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)

(5)

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

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

2.4 求解控制及邊界條件

求解控制：瞬態(tài)計(jì)算中，采用有限體積法(FVM)對(duì)網(wǎng)格進(jìn)行離散，時(shí)間步長(zhǎng)為1.0×10-4s，所有參數(shù)收斂標(biāo)準(zhǔn)均為1.0×10-6，每一時(shí)間步最大迭代步數(shù)為100步，經(jīng)后期觀測(cè)發(fā)現(xiàn)迭代60次以內(nèi)均能收斂.速度、壓力和密度的耦合方式選擇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ì)算時(shí)長(zhǎng)5 s.

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

3 計(jì)算結(jié)果

3.1 瞬變壓力三維CFD計(jì)算與試驗(yàn)結(jié)果對(duì)比

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

表2 氣團(tuán)最大壓力值的CFD值和試驗(yàn)值

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

圖3 2個(gè)氣團(tuán)計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比

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

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

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

圖4 一維-CFD-試驗(yàn)結(jié)果對(duì)比

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

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

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

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

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

4 結(jié) 論

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

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

3) 一維彈性模型計(jì)算結(jié)果誤差大于CFD結(jié)果，這是因?yàn)橐痪S模型無(wú)法有效地模擬水氣交界面變化過(guò)程、阻斷水體長(zhǎng)度的動(dòng)態(tài)變化.

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