劉光,吳國鴻,吳如坤

(1. 海軍潛艇學(xué)院,山東 青島 266042; 2. 清華大學(xué)能源與動力工程系流體機械及工程研究所,北京100084)

壓力輸水管道是灌溉管網(wǎng)、民用輸調(diào)水工程、水泵和水電站等現(xiàn)代管道水力系統(tǒng)中的重要組成部分,在潛艇脫險系統(tǒng)及岸港模擬訓(xùn)練系統(tǒng)中也有廣泛應(yīng)用,其穩(wěn)定運行對可靠、安全地輸送液體至關(guān)重要[1-2].在壓力管道中,瞬變流會造成沿管道傳播的壓力波和水錘現(xiàn)象,是系統(tǒng)可變載荷的來源,在特定的情況下會造成水力系統(tǒng)的嚴重故障[3].因此,在各類水力系統(tǒng)的設(shè)計、制造和運行階段,可靠的瞬變流動分析方法和抵消其破壞性影響的手段至關(guān)重要.盡管經(jīng)過多年的研究,管道中的瞬變流動及相關(guān)的水錘現(xiàn)象還是沒能被深入認識和數(shù)學(xué)定義.

近年來,管道瞬變流動的數(shù)值模擬研究已取得了很大進展,在空間尺度上,由1D模擬發(fā)展到(準)2D模擬,乃至3D模擬.其中2D模擬計算考慮了管道截面的速度分布,相比1D模擬能夠獲得更多的流動信息,但由于其采取了軸對稱假設(shè),從而無法考慮動量方程中的徑向速度分量以及導(dǎo)數(shù),并且其假設(shè)橫截面壓力為常數(shù),這與實際情況并不符合[4].3D模擬方法可以同時考慮管道內(nèi)3個方向的速度分量,獲得管道系統(tǒng)的瞬態(tài)特性以及詳細的流場信息,因此被廣泛應(yīng)用于包括管道非對稱流動特性在內(nèi)的整體流場信息研究.但其缺點是對于計算資源的要求很高,全流道計算時間很長,難以滿足大多數(shù)工程實際的應(yīng)用需求[5].相比于2D及3D模擬,1D模擬以其精度高、處理簡單、消耗計算資源少等特點,得到了廣泛應(yīng)用.

目前大多數(shù)相關(guān)文獻中1D瞬變流模型中采用的是穩(wěn)態(tài)假說,即假定管道瞬態(tài)流動中的摩擦損失與穩(wěn)定流動條件下確定的摩擦損失相等.該假說被認為是造成試驗測試結(jié)果與水錘過程中壓力變化計算結(jié)果之間差異的最重要原因之一.為了獲得可靠的水錘模擬結(jié)果,需要充分考慮對瞬變流動中摩擦損失的影響.目前已有的摩阻模型包括基于層流和湍流流動開發(fā)的模型[6-8],這些模型雖然能夠獲得較準確的結(jié)果,但對計算的要求非常高.在進一步的研究中,學(xué)者們[9-11]提出了考慮水錘過程摩擦力瞬時效應(yīng)的Zielke和Vardy-Brown近似模型(等效模型).另一方面,BRUNONE等[12]提出了基于瞬時加速度的MIAB (modified instantaneous acceleration based model, MIAB)模型,該模型在預(yù)測水錘波衰減時得到了較準確的結(jié)果,尤其對于高雷諾數(shù)瞬變流動.

在試驗方面,一些學(xué)者[13-14]通過建立水力管道系統(tǒng),對管道瞬變流動的壓力瞬變特性進行了研究.還有一些學(xué)者[15-19]采用PIV,LDV等先進測試技術(shù)探究了瞬變流動局部流場和壁面剪切應(yīng)力等特性.除此之外,EGGELS等[20]介紹了大量非空化條件下管道中的瞬變流動研究,為摩阻模型的對比驗證提供了重要參考.此外,MARTINS等[21]建立了典型水箱-管道-閥門系統(tǒng),并包含確定閥門關(guān)閉規(guī)律下不同雷諾數(shù)工況的試驗數(shù)據(jù),也能用于對不同摩阻模型的驗證.

文中在1D恒定摩阻模型研究的基礎(chǔ)上,引入包含恒定摩阻和附加瞬時摩阻的改進非恒定MIAB摩阻模型,并結(jié)合壓力管道的特征線法,對壓力管道的瞬變流動開展理論與數(shù)值研究,并將不同模型計算得到的壓力演變與試驗測量數(shù)據(jù)進行對比,以進一步驗證不同摩阻模型在管道瞬變流動計算中的適用性.

1 瞬變流摩阻模型

為了對比管道瞬變流模型的適用性和預(yù)測精度,采用1D恒定摩阻模型(1D-SF)和1D非恒定摩阻模型(1D-UF)開展數(shù)值計算.

1.1 瞬變流動控制方程

管道瞬變流動的1D控制方程含如下連續(xù)性方程和運動方程

(1)

(2)

式中:v為平均流速,m/s;H為壓力水頭,m;a為管中壓力波速,m/s;g為重力加速度,m2/s;φ為管道傾角,(°);J為水力損失,m.

壓力波速a計算式[5]為

(3)

式中:K為流體彈性模量;E為管道彈性模量;e為管道壁厚,且在模型中,同時考慮了管道和流體的可壓縮性;ρ為密度,kg/m3.

在1D恒定摩阻模型中,只考慮管道的穩(wěn)態(tài)摩阻Jq,即

(4)

式中:D為管道直徑,m;f為管道的達西-維斯巴赫摩阻系數(shù).

而在1D非恒定摩阻模型中,同時考慮了管道的穩(wěn)態(tài)摩阻Jq和非恒定摩阻Ju,

J=Jq+Ju,

(5)

其中,非恒定摩阻Ju為

(6)

式中:k為Brunone摩阻系數(shù);sign為符號函數(shù),采用計算式為

(7)

1.2 特征線法

基于特征線法,將一維恒定摩阻模型的偏微分方程(1),(2)轉(zhuǎn)化為沿著正負特征線的常微分方程(8),(9),即

(8)

(9)

一維非恒定摩阻模型的相應(yīng)常微分方程表示為

(10)

(11)

將1D水錘方程組(8),(9) 和(10),(11)進一步轉(zhuǎn)化為差分方程組,并利用C++編程求解,得到有壓管道1D瞬變過程中管道任意斷面的水頭H和流量Q隨時間變化的規(guī)律.

1.3 計算模型和工況

利用MARTINS等[21]的試驗數(shù)據(jù)對上述數(shù)值方法的準確性進行驗證.如圖1所示,試驗系統(tǒng)包括1個具有恒定水頭的上游水箱、1根銅管和1個位于管道下游端部的球閥.管道中的瞬變流是由閥門的快速關(guān)閉引起的,關(guān)閥時長為14 ms.使用位于閥門上游(PT1)和整個管道中部(PT2)的2個壓力傳感器測量壓力變化,數(shù)據(jù)采集頻率為3 kHz.圖中的X,Y和Z坐標分別表示流動的軸向、垂直和水平方向.

圖1 幾何模型示意圖

管道系統(tǒng)主要參數(shù):管長l=15.22 m,壁厚e=1 mm,管道直徑D=20 mm.水錘波速為1 250 m/s.對試驗中3種不同初始流動工況進行計算和分析.主要的參數(shù)如表1所示,表中v0為初始流速,Re為初始雷諾數(shù),H0為上游水箱水頭.

表1 測試工況初始流動參數(shù)

2 試驗結(jié)果和討論

2.1 1D摩阻模型適用性研究

由Case1工況的計算得到PT1測壓點的瞬時壓力演變與試驗數(shù)據(jù)對比,如圖2 所示.

圖2 Case1工況PT1壓力隨時間變化

從圖2可以看出,在壓力演變的第1個周期,1D-SF和1D-UF模型預(yù)測的振幅和相位都與試驗數(shù)據(jù)接近.這說明2個模型對于閥門關(guān)閉引起的壓力升高以及此后由上游水箱返回的降壓波的計算具有較好的精度.進一步提取0～0.050 s的數(shù)據(jù),如圖3所示,可以看出第1個壓力波周期內(nèi)2種模型的預(yù)測曲線幾乎重合,并且與試驗測試值相差很小.

圖3 Case1工況PT1壓力隨時間變化(0～0.050 s)

在第1個壓力波周期后,隨著水錘波在管道內(nèi)的傳播,試驗測試得到的壓力波逐漸衰減,壓力波幅值不斷減小.然而,1D-SF 模型僅僅精確預(yù)測了第1個壓力峰值,在第1個壓力峰值后,水錘波的幅值幾乎沒有減小,模型預(yù)測結(jié)果與試驗數(shù)據(jù)的誤差在之后的4個水錘波周期內(nèi)不斷增大.與之相比,1D-UF 模型預(yù)測得到的壓力幅值與試驗值符合較好,基本預(yù)測到了第1個周期之后的壓力波衰減.在低流速的Case1工況,1D-SF模型和1D-UF模型均準確預(yù)測到了第1個周期內(nèi)的壓力波動,但1D-SF模型幾乎沒有預(yù)測到第1個周期之后的壓力耗散,而1D-UF模型對壓力耗散的預(yù)測也和試驗值較吻合.1D-SF 模型和1D-UF 模型在0.220 s 附近的壓力峰值相對誤差分別為10.3%和4.7%,并且誤差值會隨著時間增長逐漸增大.

為了驗證2種模型在瞬變壓力的預(yù)測效果在相對較高流速工況下是否相同,進一步對Case2和Case3工況的壓力演變進行了分析.

圖4和圖5分別為Case2和Case3工況的計算結(jié)果與試驗結(jié)果對比.對于Case2和Case3,壓力波的整體演變過程與Case1基本相似.然而由于管道內(nèi)初始流速更高,因此壓力峰值也更大.1D-SF模型對壓力峰值預(yù)測的相對誤差在Case2和Case3工況均大于1D-UF模型,這說明在高流速工況下,1D-UF模型對于管道瞬變流動的預(yù)測仍然與試驗值吻合較好.

圖4 Case2工況PT1壓力隨時間變化

圖5 Case3工況PT1壓力隨時間變化

2.2 管道中點壓力瞬變特性

利用1D-UF模型,對于管道中點的壓力瞬變特性進行了分析.圖6—8分別為由1D-SF模型和1D-UF模型預(yù)測得到的Case1,Case2和Case3工況管道中點壓力演變.

圖6 Case1工況PT2壓力隨時間變化

圖7 Case2工況PT2 壓力隨時間變化

圖8 Case3工況PT2壓力隨時間變化

在Case1工況,根據(jù)壓力波的演變特性,可以將1個完整壓力波周期分為4個階段.

第1階段:閥門快速關(guān)閉,離閥門較近的上游處流動首先停止,在該處引起了升壓波,該升壓波向上游傳播,其傳播到的地方壓力增加,當(dāng)t=0.006 s時,升壓波傳到管道中點處,導(dǎo)致此處壓力開始急劇升高.直至t=0.012 s時,升壓波傳至上游水箱,使得整個管道內(nèi)為增壓狀態(tài).

第2階段:第1個階段的升壓波傳至上游水箱時,由于水箱內(nèi)的壓力保持恒定,而管道內(nèi)的液體壓力大于上游水箱,從而產(chǎn)生壓差,在壓差作用下,使得上游水箱處產(chǎn)生與原先流速相反的水流,在其作用下,管路內(nèi)壓強恢復(fù)至原來壓力值.該過程即壓力波的全負反射,升壓波在水箱處轉(zhuǎn)換為降壓波,并向下游傳播,在t=0.024 s時,降壓波傳至管道中點,使得管道中點處壓力快速下降至原先壓力值.直至t=0.030 s時,降壓波傳至閥門,整個管道內(nèi)壓力恢復(fù)至原始值.

第3階段:第2個階段的降壓波傳播至閥門處,此時管道內(nèi)的流速與初始流速方向相反,而閥門已完全關(guān)閉,閥門處微小流段速度變化為0,在此處引起了降壓波,即第2階段的降壓波在閥門處被反射,從而改變方向,向上游傳播,在t=0.036 s時,該降壓波傳播至管道中點,造成此處壓力驟降.到t=0.042 s時,降壓波再次傳至上游水箱處,此時整個管道均為降壓狀態(tài).

第4階段:第3個階段的降壓波傳播至上游水箱處,水箱內(nèi)的壓力大于管道內(nèi)壓力,在該壓差影響下,上游水箱處產(chǎn)生和原先流動方向相同的流速,管路內(nèi)的壓力恢復(fù)至原先壓力值,即降壓波在此處發(fā)生全負反射,轉(zhuǎn)換為升壓波,該升壓波繼續(xù)向下游傳播,在t=0.048 s時,升壓波傳至管道中點,而在t=0.054 s時,升壓波傳至閥門處,導(dǎo)致整個管道內(nèi)的壓力恢復(fù)至與初始值接近的狀態(tài).

以上即為1個完整壓力波周期內(nèi)4個階段壓力波的演變特性,在此之后,壓力波在管道內(nèi)不斷循環(huán),但由于受到管道摩阻和瞬時摩阻的作用,壓力幅值不斷減小,壓力波動不斷減弱,最終在管道內(nèi)消失.在Case2和Case3工況,管道中點處的壓力波瞬變特性與Case1基本相似,只是由于初始流速更高,壓力幅值更大.

3 結(jié) 論

1) 恒定摩阻模型和非恒定摩阻模型在壓力波第1個周期內(nèi)的壓力變化預(yù)測值均與試驗值較為接近.在考慮到閥門關(guān)閉過程壓力升高的情況下,2種模型對于第1個壓力波周期內(nèi)幅值和相位的預(yù)測值均與試驗值吻合較好.

2) 在第1個壓力波周期以后,MIAB模型能很好地預(yù)測管道內(nèi)不同位置、不同初始流速工況下的壓力波幅值和相位,在低流速工況,預(yù)測閥門處壓力峰值與試驗值誤差為4.7%,相比于恒定摩阻模型10.3%的誤差,預(yù)測精確度大幅提高.2種模型的計算精度主要體現(xiàn)在對于第1個壓力波周期后壓力波衰減的預(yù)測上.

3) 根據(jù)壓力波傳播特性,1個完整水錘波周期可以分為4個階段,并且壓力波在這4個階段內(nèi)不同位置處的傳播和反射特性存在明顯不同.在水箱處,由于水箱和管道內(nèi)的壓強差,導(dǎo)致壓力波在此處形成全負反射,即第1階段向上游傳播的升壓波轉(zhuǎn)換為方向相反的降壓波,而第3階段的降壓波轉(zhuǎn)換為方向相反的升壓波.而在閥門處,由于閥門完全關(guān)閉后,緊靠閥門的微流段速度變化為0,因此第2階段的降壓波在閥門處僅僅改變方向,轉(zhuǎn)換為向上游傳播的降壓波.