楊衛(wèi)彬，王煥茂，張 韜

(哈爾濱電機廠有限責任公司，黑龍江 哈爾濱 150040)

當水泵水輪機機組啟動或停機時，駝峰區(qū)會引發(fā)機組的振動、噪聲等不穩(wěn)定現(xiàn)象，從而直接影響機組的運行[1]，因此需要在水力開發(fā)階段對駝峰區(qū)壓力脈動予以特別關(guān)注。李德友等對水泵水輪機不同導(dǎo)葉開口的駝峰特性進行了研究[2]；陳順義等人對水泵水輪機穩(wěn)定性判斷進行過研究[3]；國內(nèi)還有很多學者對水泵水輪機導(dǎo)葉開度等方面進行過研究[4-8]；但對駝峰區(qū)多個工況點進行壓力脈動數(shù)值模擬分析，同時與試驗值進行比較和分析的還很少。本文采用專業(yè)的旋轉(zhuǎn)機械領(lǐng)域設(shè)計與分析軟件Ansys CFX軟件，通過非定常數(shù)值分析的方式，采用DES數(shù)值模擬方法對駝峰區(qū)多個不同工況的壓力脈動頻譜特性和幅值進行了分析研究，同時對駝峰區(qū)流態(tài)進行了分析研究。

1 湍流模型與計算方法的選擇

水力機械中的湍流運動可用連續(xù)方程和Navier-Stokes方程來描述：

(1)

(2)

式中：ρ為水的密度；u為速度矢量；F為質(zhì)量力；p為壓力；μ為動力粘度。

式(2)左邊第一項為由于運行非定常性而引起的局部慣性力，第二項為由于運動的非均勻性而引起的變位慣性力，右邊三項分別為質(zhì)量力、壓力和粘性力。通?？蓪⑺靼床豢蓧嚎s處理，故方程可以簡化。對于湍流的處理方式，式(2)中的最后一項是決定流場數(shù)值模擬精度與效率的關(guān)鍵。各種k-ε雙方程模型及雷諾應(yīng)力模型是雷諾時均法的幾種不同形式，本文選擇RNGk-ε雙方程模型和DES模擬方法同時對混流式水泵水輪機駝峰區(qū)的非穩(wěn)定流動進行計算。

計算的水力模型采用9葉片混流式水泵水輪機，其高壓邊直徑為470 mm，活動導(dǎo)葉和固定導(dǎo)葉各20個。計算域包括蝸殼、固定導(dǎo)葉、活動導(dǎo)葉、轉(zhuǎn)輪和尾水管，采用非結(jié)構(gòu)化網(wǎng)格，網(wǎng)格數(shù)量為323萬，水泵水輪機全流道數(shù)值模擬的計算域如圖1所示。

圖1 水泵水輪機全流道數(shù)值模擬計算域圖

使用有限體積法對控制方程進行離散，尾水管采用進口質(zhì)量流量，蝸殼采用出口相對靜壓P=0 Pa；壁面采用無滑移邊界條件處理，并采用壁面函數(shù)對近壁區(qū)域進行流動模擬；動靜交界面采用interface模擬。

通過定常計算得到初始流場，非定常計算時以定常計算的結(jié)果作為初始條件。選取駝峰區(qū)的多個工況點進行數(shù)值模擬，計算總時長為10個旋轉(zhuǎn)周期，時間步設(shè)為旋轉(zhuǎn)周期的1/512，步長是0.000 109 5 s，以10-4作為收斂殘差計算。

2 仿真計算結(jié)果

測點如圖2布置，包含無葉區(qū)壓力脈動測點P1和錐管壓力脈動測點P2。

圖2 壓力脈動監(jiān)測點分布圖

圖3、圖4分別為試驗的壓力脈動幅值與數(shù)值模擬的壓力脈動幅值的對比圖。從圖3、圖4中可見，無論是試驗還是數(shù)值模擬，壓力脈動幅值都是隨著流量的增大而減小，兩者趨勢具有很好的一致性。同時需要注意到，試驗與數(shù)值模擬的壓力脈動幅值在同一工況點相差較大，原因還要進一步的分析。

表1為數(shù)值模擬的無葉區(qū)壓力脈動主頻，從中可以看出，無葉區(qū)壓力脈動主頻隨著流量的增大而增大，但最小流量的兩個工況點主頻又略大，可依舊是低頻脈動。

圖3 P1測點壓力脈動幅值

圖4 P2測點壓力脈動幅值

表1 數(shù)值模擬的無葉區(qū)壓力脈動主頻表

圖5和圖6是Z=0.03 m平面內(nèi)的速度矢量分布圖，可以看出流量越小，固定導(dǎo)葉、活動導(dǎo)葉和轉(zhuǎn)輪間的渦流越強烈，無葉區(qū)壓力脈動由原來的轉(zhuǎn)輪旋轉(zhuǎn)引起，變成由渦流引起。由此可見，水泵水輪機駝峰區(qū)內(nèi)部流動的主要特征就是渦流運動。

3 采用DES數(shù)值方法模擬結(jié)果

分離渦模型(DES)方法綜合了非定常雷諾平均方法(URANS)和大渦模擬方法(LES)的優(yōu)點，在近壁區(qū)采用URANS模擬，遠離近壁處采用LES模擬，目前在大攻角翼型等分離流動中得到廣泛應(yīng)用[9]。

圖5 Q/Qid=0.8，Z=0.03 m平面速度矢量分布

圖6 Q/Qid=0.70，Z=0.03 m平面速度矢量分布

計算域中蝸殼和尾水管部分采用四面體單元，導(dǎo)葉和轉(zhuǎn)輪區(qū)域采用棱柱形和四面體相結(jié)合的混合單元。為了滿足DES計算對近壁網(wǎng)格的要求，導(dǎo)葉和轉(zhuǎn)輪的近壁處進行了加密，第一個節(jié)點的位置距壁面的距離為0.001 mm，網(wǎng)格高度比設(shè)為1.35，y+<5，網(wǎng)格單元數(shù)合計926.1萬。計算選用水泵水輪機泵工況某導(dǎo)葉開度下駝峰區(qū)的1個工況點Q/Qid=0.64，對應(yīng)單位流量Q11=0.438 L/s。

由圖7、圖8可以看出，無葉區(qū)壓力脈動DES數(shù)值模擬頻譜圖與試驗頻譜圖相比，壓力脈動幅值比較大的部分均集中在100 Hz以下，低頻部分壓力脈動幅值均較大，反映出DES數(shù)值模擬能夠?qū)o葉區(qū)渦流動進行比較準確的模擬；由圖9、圖10可以看出，錐管壓力脈動頻譜圖則反映出DES數(shù)值模擬與試驗均具有很高的一致性，二者所得的主頻幾乎是一致的，表明DES數(shù)值模擬對錐管內(nèi)部渦流的模擬是很準確的。圖11為采用RNGk-ε湍流模型、DES湍流模型進行非定常數(shù)值模擬所得無葉區(qū)壓力脈動幅值與模型試驗測量值的比較?？梢钥闯?，與試驗值相比，DES數(shù)值模擬對無葉區(qū)壓力脈動幅值的預(yù)測比采用RNGk-ε湍流模型更為準確。由以上對壓力脈動分析可以看出，采用DES湍流模型對水泵水輪機駝峰區(qū)進行的非定常模擬是可以對駝峰區(qū)壓力脈動進行準確的預(yù)測，在此基礎(chǔ)上對駝峰區(qū)內(nèi)部流動進行分析是可以得到可信的結(jié)果。

圖7 無葉區(qū)壓力脈動模型試驗頻譜圖(對應(yīng)紅色曲線)

圖8 無葉區(qū)壓力脈動DES數(shù)值模擬頻譜圖

圖9 錐管壓力脈動模型試驗頻譜圖(對應(yīng)紅色曲線)

圖10 錐管壓力脈動DES數(shù)值模擬頻譜圖

圖11 采用不同湍流模型所得無葉區(qū)壓力脈動幅值與試驗值比較圖

圖12至圖15是轉(zhuǎn)輪一個旋轉(zhuǎn)周期內(nèi)不同時刻導(dǎo)葉中心線平面以上30 mm處壓力分布圖。由圖中可以看出，不論轉(zhuǎn)輪、雙列葉柵還是蝸殼內(nèi)壓力分布都很不均勻，在轉(zhuǎn)輪出口壓力面出現(xiàn)局部的高壓區(qū)，該局部高壓區(qū)在不同時刻的大小和位置都隨轉(zhuǎn)輪的旋轉(zhuǎn)而變化，其壓力值大小甚至有3～4倍的增加(如圖15所示)，這將導(dǎo)致駝峰區(qū)工況無葉區(qū)較高的壓力脈動幅值。圖16至圖19是轉(zhuǎn)輪一個旋轉(zhuǎn)周期內(nèi)不同時刻錐管截面壓力分布圖，可以看到不同時刻錐管截面壓力分布都存在一個環(huán)狀低壓區(qū)，該環(huán)狀低壓區(qū)是由錐管內(nèi)的環(huán)狀旋流造成，其方向與轉(zhuǎn)輪旋轉(zhuǎn)方向相同；采用DES方法數(shù)值模擬所得到的錐管內(nèi)部流態(tài)與模型試驗中經(jīng)常觀測到的駝峰區(qū)錐管內(nèi)部流態(tài)具有很好的一致性。

圖12 t=0，導(dǎo)葉中心線平面以上30 mm處壓力分布圖

圖13 t=1/4T，導(dǎo)葉中心線平面以上30 mm處壓力分布圖

圖14 t=1/2T，導(dǎo)葉中心線平面以上30 mm處壓力分布圖

圖15 t=3/4T，導(dǎo)葉中心線平面以上30 mm處壓力分布圖

圖16 t=0，錐管截面壓力分布圖

圖17 t=1/4T，錐管截面壓力分布圖

圖18 t=1/2T，錐管截面壓力分布圖

圖19 t=3/4T，錐管截面壓力分布圖

4 結(jié) 語

采用RNGk-ε湍流模型對水泵水輪機駝峰區(qū)進行非定常數(shù)值模擬，可以對駝峰區(qū)壓力脈動頻譜特性和壓力脈動幅值進行初步的預(yù)測與評估。DES數(shù)值模擬方法可以對駝峰區(qū)壓力脈動進行更為準確的預(yù)測，在工程應(yīng)用領(lǐng)域，其預(yù)測精度是足夠的。

在水泵水輪機駝峰區(qū)，隨著流量的不斷減小，壓力脈動的主頻也隨之降低，渦流運動將在水泵水輪機內(nèi)部流動中起著越來越重要的作用，大量的渦流運動和流道內(nèi)局部壓力的劇烈變化是水泵水輪機駝峰區(qū)壓力脈動的主要根源。