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

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

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

1 湍流模型與計(jì)算方法的選擇

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

(1)

(2)

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

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

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

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

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

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

2 仿真計(jì)算結(jié)果

測(cè)點(diǎn)如圖2布置，包含無葉區(qū)壓力脈動(dòng)測(cè)點(diǎn)P1和錐管壓力脈動(dòng)測(cè)點(diǎn)P2。

圖2 壓力脈動(dòng)監(jiān)測(cè)點(diǎn)分布圖

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

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

圖3 P1測(cè)點(diǎn)壓力脈動(dòng)幅值

圖4 P2測(cè)點(diǎn)壓力脈動(dòng)幅值

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

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

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

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

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

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

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

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

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

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

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

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

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

圖12至圖15是轉(zhuǎn)輪一個(gè)旋轉(zhuǎn)周期內(nèi)不同時(shí)刻導(dǎo)葉中心線平面以上30 mm處壓力分布圖。由圖中可以看出，不論轉(zhuǎn)輪、雙列葉柵還是蝸殼內(nèi)壓力分布都很不均勻，在轉(zhuǎn)輪出口壓力面出現(xiàn)局部的高壓區(qū)，該局部高壓區(qū)在不同時(shí)刻的大小和位置都隨轉(zhuǎn)輪的旋轉(zhuǎn)而變化，其壓力值大小甚至有3～4倍的增加(如圖15所示)，這將導(dǎo)致駝峰區(qū)工況無葉區(qū)較高的壓力脈動(dòng)幅值。圖16至圖19是轉(zhuǎn)輪一個(gè)旋轉(zhuǎn)周期內(nèi)不同時(shí)刻錐管截面壓力分布圖，可以看到不同時(shí)刻錐管截面壓力分布都存在一個(gè)環(huán)狀低壓區(qū)，該環(huán)狀低壓區(qū)是由錐管內(nèi)的環(huán)狀旋流造成，其方向與轉(zhuǎn)輪旋轉(zhuǎn)方向相同；采用DES方法數(shù)值模擬所得到的錐管內(nèi)部流態(tài)與模型試驗(yàn)中經(jīng)常觀測(cè)到的駝峰區(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-ε湍流模型對(duì)水泵水輪機(jī)駝峰區(qū)進(jìn)行非定常數(shù)值模擬，可以對(duì)駝峰區(qū)壓力脈動(dòng)頻譜特性和壓力脈動(dòng)幅值進(jìn)行初步的預(yù)測(cè)與評(píng)估。DES數(shù)值模擬方法可以對(duì)駝峰區(qū)壓力脈動(dòng)進(jìn)行更為準(zhǔn)確的預(yù)測(cè)，在工程應(yīng)用領(lǐng)域，其預(yù)測(cè)精度是足夠的。

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