劉兵洋 ，楊 魏 ，王福軍 ，劉竹青 ，黎耀軍

（1.中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，北京市 100083；2.北京市供水管網(wǎng)系統(tǒng)安全與節(jié)能工程技術(shù)研究中心，北京市 100083）

0 引言

水泵水輪機是抽水蓄能電站的核心部件[1]，其特有的S特性對機組的安全穩(wěn)定運行造成了一定的不良影響。研究表明，該特性由水輪機制動工況轉(zhuǎn)輪內(nèi)部不穩(wěn)定流動所造成的巨大損失所引起[2]。

近年來，國內(nèi)外研究人員對S特性區(qū)內(nèi)的不穩(wěn)定流動作了大量的研究。Widmer等[3]人通過試驗以及數(shù)值模擬，將S區(qū)內(nèi)的不穩(wěn)定流動劃分為 Stationary vortex（固定渦）、Unsteady vortex（非定常渦）以及Rotating stall（旋轉(zhuǎn)失速），并根據(jù)其各自壓力脈動以及頻率特征加以區(qū)分。ZHANG等[4]人在對水泵水輪機飛逸工況下無葉區(qū)壓力波動轉(zhuǎn)換機理的探索過程中，將S區(qū)內(nèi)轉(zhuǎn)輪內(nèi)產(chǎn)生的渦結(jié)構(gòu)分為FFVS（同主流方向傳播的渦）和BFVS （回流渦）。但它們均未進一步闡述它們之間的相互關(guān)系。

為了進一步理解水泵水輪機S特性區(qū)轉(zhuǎn)輪內(nèi)不穩(wěn)定流動產(chǎn)生的原因及演化機理，采用數(shù)值模擬的手段對水泵水輪機S特性附近工況轉(zhuǎn)輪內(nèi)的流場進行研究。

1 數(shù)值模擬方法

1.1 計算域與網(wǎng)格劃分

選用某模型水泵水輪機的全流道作為計算域，如圖1所示，部件包括蝸殼、固定導(dǎo)葉、活動導(dǎo)葉、轉(zhuǎn)輪以及尾水管，為了減小尾水渦帶對出口邊界的影響，對尾水管出口進行延長。水泵水輪機的模型參數(shù)如表1所示，其中，ns為水輪機工況比轉(zhuǎn)速，D1、D2分別為轉(zhuǎn)輪低壓邊與高壓邊直徑，Zb，Zgv和Zsv分別表示轉(zhuǎn)輪葉片數(shù)、活動導(dǎo)葉數(shù)和固定導(dǎo)葉數(shù)，αmax表示導(dǎo)葉最大開度，計算活動導(dǎo)葉開度為20.03°。

圖1 模型水泵水輪機計算域Figure 1 Computational domainof model pump turbine

表1 模型水泵水輪機幾何尺寸Table 1 Parameters of the model pump-turbine

網(wǎng)格部分運用ANSYS ICEM對各過流部件進行結(jié)構(gòu)化網(wǎng)格劃分，為了保證捕捉到轉(zhuǎn)輪內(nèi)的流動分離，對近壁區(qū)網(wǎng)格進行加密，以保證轉(zhuǎn)輪葉片表面平均y+<10。在保證轉(zhuǎn)輪第一層網(wǎng)格高度不變的情況下，分別劃分出網(wǎng)格節(jié)點數(shù)從201萬到550萬的5套網(wǎng)格，其水頭相對誤差如圖2所示，從圖中可以看出，當(dāng)網(wǎng)格節(jié)點數(shù)達到320萬時，水頭波動減小到3%以內(nèi)；此外，隨著節(jié)點數(shù)的增多，葉片表面平均y+逐漸保持在3左右。因此，考慮到計算精度及葉片表面y+，最終選擇網(wǎng)格節(jié)點數(shù)為320.67萬進行后續(xù)計算。

圖2 網(wǎng)格無關(guān)性分析Figure 2 Mesh independence analysis

1.2 湍流模型及邊界條件

采用商業(yè)軟件ANSYS CFX進行流場計算，定常計算采用SST k-ω湍流模型，非定常計算采用基于SST k-ω的DES（Detached Eddy Simulation）湍流模型。該模型結(jié)合了RANS和LES的優(yōu)點，其主要思想是通過建立函數(shù)來改變不同位置處的湍流長度尺度，最終達到在壁面附近求解雷諾平均N-S方程，在遠離壁面的主流區(qū)采用大渦模擬的算法。這樣既能在近壁面發(fā)揮雷諾平均法計算量小的優(yōu)點，也能避免高雷諾數(shù)下大渦模擬對網(wǎng)格要求太高的缺點[5]?，F(xiàn)有研究表明，DES湍流模型在流體機械中存在大分離的流動中能夠取得很好的結(jié)果[6]。

邊界條件：在蝸殼進口設(shè)置質(zhì)量流量邊界，尾水管出口設(shè)置為靜壓邊界。固體壁面設(shè)置為無滑移壁面條件。動靜交界面在定常計算中采用Frozen Rotor，在非定常設(shè)置中動靜交接壁面設(shè)置為Transient Rotor Stator，即滑移網(wǎng)格壁面。設(shè)置轉(zhuǎn)輪每旋轉(zhuǎn)1°的時間為一個時間步，每步收斂殘差標(biāo)準(zhǔn)RMS=10-5，最大迭代步數(shù)為20步，共計算15個周期。

1.3 監(jiān)測點設(shè)置

為了分析流場內(nèi)的壓力脈動特征，在流場的特定位置布置監(jiān)測點。其中，在活動導(dǎo)葉中部以及固定導(dǎo)葉中部各布置一個監(jiān)測點，活動導(dǎo)葉與轉(zhuǎn)輪之間的無葉區(qū)沿周向布置20個壓力監(jiān)測點，即圓周方向每18°布置一個，此外，在轉(zhuǎn)輪通道中沿徑向依次布置3個壓力監(jiān)測點，監(jiān)測點位置如圖3所示。

圖3 監(jiān)測點位置示意圖Figure 3 Schematic diagram of monitoring points

2 結(jié)果與分析

為了得到模型水泵水輪機內(nèi)部壓力脈動以及渦結(jié)構(gòu)的演化，對20.03°導(dǎo)葉開度的10個工況點進行模擬計算。圖4是試驗與數(shù)值模擬結(jié)果對比圖，其中，n11為單位轉(zhuǎn)速，Q11為單位流量，按照式（1）、式（2）定義。

圖4 S特性曲線試驗與模擬比較Figure 4 Comparisons between the experimental and numerical characteristic curves

式中：n——轉(zhuǎn)輪轉(zhuǎn)速，r/min；

Q——質(zhì)量流量，kg/s；

D1——轉(zhuǎn)輪低壓邊直徑，m；

H——水輪機水頭，m。

在圖4中，OP.1～OP.7位于水輪機工況，OP.8～OP.10位于水輪機制動工況。OP.7和OP.8分別位于飛逸工況兩側(cè)，在水輪機制動工況OP.8水輪機特性曲線表現(xiàn)為正斜率?？梢钥闯?，由CFD模擬出的結(jié)果與實驗結(jié)果基本吻合，單位轉(zhuǎn)速n11和單位流量Q11的誤差均在3%以內(nèi)。由此可以得出，本研究采用的數(shù)值模擬方法是可靠的。

2.1 壓力脈動特征

為了獲得水泵水輪機在水輪機工況以及水輪機制動工況的主要流動特征，對工況點OP.2、OP.5、OP.7和OP.8的壓力脈動特征進行分析。4個工況點的詳細參數(shù)及所處工況如表2所示。

表2 工況點參數(shù)（試驗）Table 2 Operating parameters （experiment）

對采集到的壓力信號進行無量綱化處理，用ΔH/H表示各監(jiān)測點壓力脈動峰值，其定義如式（3）所示。

式中：pi——瞬時壓力值，Pa；

p1——蝸殼進口壓力，Pa；

ρ——密度，一般取1000kg/m3；

g——常量，一般取9.8N/kg；

H——各工況水頭模擬值。

圖5顯示了工況點OP.5和OP.8的固定導(dǎo)葉、活動導(dǎo)葉、無葉區(qū)、轉(zhuǎn)輪進口壓力脈動峰值。從圖5（a）可以看出，隨著流量減小，在水輪機工況也會出現(xiàn)較大壓力波動，其中，轉(zhuǎn)輪內(nèi)的壓力脈動幅值最大，無葉區(qū)次之。在水輪機制動工況，如圖5（b）所示，除轉(zhuǎn)輪區(qū)域外，其余流道各處的壓力脈動幅值均不同程度增加，無葉區(qū)增幅最大，受轉(zhuǎn)輪及無葉區(qū)影響，活動導(dǎo)葉內(nèi)壓力脈動幅度也逐漸增加?？梢钥闯觯S著工況點偏離最優(yōu)工況，水泵水輪機的流動不穩(wěn)定現(xiàn)象首先集中在轉(zhuǎn)輪中，隨著工況點進入水輪機制動工況，無葉區(qū)及活動導(dǎo)葉內(nèi)不穩(wěn)定流動現(xiàn)象隨之加劇。

圖5 不同工況壓力脈動時域圖Figure 5 Pressure fluctuation under different operating mode

圖6 表示不同工況無葉區(qū)的壓力信號時域圖以及頻域圖。從時域圖可以看出，隨著流量減小，無葉區(qū)的壓力脈動幅值逐漸增加，也變得更加無序。頻率特性表明，在所有工況均出現(xiàn)了9fn的葉片頻率（fn為葉片旋轉(zhuǎn)頻率）及其他由動靜干涉引起的倍頻。在水輪機制動工況OP.8出現(xiàn)了0.57fn的低頻成分，且表現(xiàn)為主頻，這說明在水輪機制動工況，除了動靜干涉外，其他因素也引起了無葉區(qū)空間的壓力波動。為了進一步研究該波動形成原因，需要進一步對轉(zhuǎn)輪內(nèi)的流動特征進行分析。

圖6 不同工況無葉區(qū)壓力脈動特征Figure 6 Pressure fluctuation in vaneless space under different mode

2.2 轉(zhuǎn)輪內(nèi)流態(tài)分析

現(xiàn)有研究表明[7]，在較大導(dǎo)葉開度的S特性區(qū)內(nèi)，轉(zhuǎn)輪內(nèi)部會出現(xiàn)非定常渦以及旋轉(zhuǎn)失速[8]，本節(jié)基于對這兩個不穩(wěn)定特征的捕捉，進而描述其主要特征及主要影響。

2.2.1 非定常渦

圖7 為工況點OP.5、OP.7和OP.8的轉(zhuǎn)輪進口的壓力脈動特征。3個工況點的主頻及振幅如表3所示。圖8～圖10是對應(yīng)3個工況為0.9、0.5以及0.1葉輪展向位置處的流線圖。流線圖表明，在工況點OP.5，渦結(jié)構(gòu)主要靠近hub側(cè)，即0.1倍葉輪展向位置；且在所有葉片通道中都存在有渦結(jié)構(gòu)，渦結(jié)構(gòu)并未完全阻塞流道，也未出現(xiàn)明顯回流，該渦結(jié)構(gòu)未對上游流場產(chǎn)生明顯影響。結(jié)合頻域圖可以得出，0.2fn的壓力波動主要由此刻轉(zhuǎn)輪內(nèi)的渦結(jié)構(gòu)演化所造成，該渦結(jié)構(gòu)與Widmer等人[3]所描述的非定常渦特征相同。

表3 不同工況頻域信息Table 3 Frequency domain under different mode

圖7 不同工況轉(zhuǎn)輪進口頻域圖Figure 7 Frequency domain on runner inlet under different conditions

圖8 工況OP.5 各展向流線圖Figure 8 Streamline distributions under OP.5

圖9 工況OP.7 各展向流線圖Figure 9 Streamline distributions under OP.7

圖10 工況OP.8 各展向流線圖Figure 10 Streamline distributions under OP.8

2.2.2 旋轉(zhuǎn)失速

與工況OP.5相比，工況OP.7和工況點OP.8均出現(xiàn)0.57fn的頻率成分，且對應(yīng)壓力脈動幅值更高，在工況OP.8表現(xiàn)為主頻。通過流線圖9和圖10可以看出，在兩個工況點轉(zhuǎn)輪通道內(nèi)的渦結(jié)構(gòu)都呈現(xiàn)出明顯的非對稱性，僅存在于連續(xù)的幾個流道當(dāng)中，且在部分流道造成阻塞以及回流。尤其在工況點OP.8中，葉道內(nèi)渦的非對稱分布更加明顯，流道被嚴重阻塞，回流增加，并對無葉區(qū)及上游流場產(chǎn)生影響，即轉(zhuǎn)輪內(nèi)出現(xiàn)旋轉(zhuǎn)失速。

此外，在工況點OP.7和OP.8，渦結(jié)構(gòu)在靠近hub側(cè)（0.1倍轉(zhuǎn)輪展向）最明顯，因此，用hub側(cè)的流線圖進一步分析旋轉(zhuǎn)失速傳播特征。圖11是工況點OP.8從時刻t1開始每隔一個轉(zhuǎn)輪周期T的流線圖（hub側(cè)），紅線表示失速團所在流道。從圖中可以看出，渦結(jié)構(gòu)在葉片通道中沿旋轉(zhuǎn)方向，以小于葉輪轉(zhuǎn)速的速度傳播。

圖11 工況點OP.8流線圖（0.1span）Figure 11 Streamline distributions under OP.8（0.1span）

圖12 是工況點OP.8無葉區(qū)監(jiān)測點的壓力脈動時域圖，從圖中可以很明顯看出，在無葉區(qū)的壓力脈動在該工況呈現(xiàn)出周期性，在7個周期內(nèi)出現(xiàn)了約4次壓力波動，傳播周期約為1.75倍的轉(zhuǎn)輪旋轉(zhuǎn)周期，此傳播周期與該工況無葉區(qū)0.57fn特征頻率相佐證。

圖12 工況點OP.8無葉區(qū)壓力脈動時域圖Figure 12 Time history of pressure fluctuation in vaneless space under OP.8

綜上可知，沿20.03°等開度線，從水輪機工況到水輪機制動工況，轉(zhuǎn)輪內(nèi)先后出現(xiàn)了非定常渦以及旋轉(zhuǎn)失速。其中，非定常渦沿周向分布在所有流道，旋轉(zhuǎn)失速沿轉(zhuǎn)輪周向非對稱分布，存在于連續(xù)的幾個流道當(dāng)中，并以亞同步轉(zhuǎn)速在流道中傳播，失速團在正斜率工況（OP.8）的葉片通道中發(fā)展完全。

2.3 轉(zhuǎn)輪內(nèi)不穩(wěn)定流動轉(zhuǎn)化過程

圖13是工況點OP.5和OP.8某時刻各葉片通道流量圖，其對應(yīng)的x軸表示轉(zhuǎn)輪的9個葉片通道，y軸表示該時刻瞬時流量Q與葉片平均流量（Qav）的比值。可以看出在工況點OP.5，受非定常渦沿周向分布，各葉片通道的流量均出現(xiàn)了不同程度的波動，但未造成明顯的通道阻塞現(xiàn)象。在工況點OP.7，部分流道發(fā)生阻塞，但阻塞的流道位置隨機。然而，在工況點OP.8，受旋轉(zhuǎn)失速影響，流道中的渦結(jié)構(gòu)已經(jīng)進入穩(wěn)定的周期旋轉(zhuǎn)模式，通道阻塞現(xiàn)象表現(xiàn)出連續(xù)性，即在相鄰的3個流道（流道9、1、2）流量均明顯偏離平均流量。圖14是OP.8對應(yīng)時刻的流線圖，可以觀察到在對應(yīng)的3個流道均出現(xiàn)了發(fā)展完全的失速團。

圖13 不同工況各葉片通道流量分布圖Figure 13 Flow distribution of each blade channel under different conditions

圖14 工況點OP.8流線圖（0.5span）Figure 14 Streamline distributions under OP.8（0.5span）

對轉(zhuǎn)輪內(nèi)部流場進行進一步分析發(fā)現(xiàn)，轉(zhuǎn)輪中渦結(jié)構(gòu)的轉(zhuǎn)化與葉片通道的流量變化密切相關(guān)。即隨著流量的減小，葉片通道的阻塞驅(qū)動了非定常渦逐漸表現(xiàn)出穩(wěn)定的周期性旋轉(zhuǎn)模式。如圖15所示，由于流量減小，入流角與葉片安放角存在一定沖角，導(dǎo)致水流沖擊葉片面之后向吸力面偏移，在吸力面形成擾流并最終在流道中產(chǎn)生非定常渦結(jié)構(gòu)，在所有流道均出現(xiàn)了幾乎相同的現(xiàn)象。隨著流量進一步減小，渦結(jié)構(gòu)進一步發(fā)展，部分流道出現(xiàn)阻塞現(xiàn)象但分布隨機。如圖16所示，某一流道的阻塞會造成流體向相鄰流道流動，當(dāng)受阻塞流體流入與旋轉(zhuǎn)方向相同的相鄰流道時，會沖擊葉片的吸力面，從而抑制該流道吸力面的渦結(jié)構(gòu)。然而，在與旋轉(zhuǎn)方向相反的相鄰流道，受阻塞流體的流入會加劇水流沖擊葉片的壓力面，進而在葉片背面產(chǎn)生更強的渦結(jié)構(gòu)。也就是在這樣的相互作用下，隨著流量的降低，不穩(wěn)定渦結(jié)構(gòu)逐漸展現(xiàn)出穩(wěn)定的周向傳播模式，即旋轉(zhuǎn)失速。

圖15 非定常渦形成機理示意圖Figure 15 Formation mechanism of unsteady vortex

圖16 旋轉(zhuǎn)失速形成機理示意圖Figure 16 Formation mechanism of rotating stall

旋轉(zhuǎn)失速由轉(zhuǎn)輪中的非定常渦發(fā)展而來，隨著流量的降低，轉(zhuǎn)輪通道中形成非定常渦結(jié)構(gòu)，當(dāng)某一流道的渦結(jié)構(gòu)增加導(dǎo)致通道阻塞時，會逐漸驅(qū)動不穩(wěn)定渦沿與轉(zhuǎn)輪旋轉(zhuǎn)相反的方向傳播，并最終表現(xiàn)出穩(wěn)定的周期性特征同時阻塞流道，即發(fā)展為旋轉(zhuǎn)失速。

3 結(jié)論

本文以模型水泵水輪機為研究對象，通過數(shù)值模擬方法，分析了水泵水輪機S特性區(qū)附近工況轉(zhuǎn)輪內(nèi)的主要不穩(wěn)定特征及其演化機理。主要結(jié)論如下：

沿等開度線，轉(zhuǎn)輪中依次形成非定常渦和旋轉(zhuǎn)失速。其中，旋轉(zhuǎn)失速表現(xiàn)出穩(wěn)定的周期性特征，失速團造成通道嚴重阻塞，并以亞同步轉(zhuǎn)速在葉片通道內(nèi)傳播。分析表明：旋轉(zhuǎn)失速由轉(zhuǎn)輪內(nèi)的不穩(wěn)定渦結(jié)構(gòu)發(fā)展而來，渦結(jié)構(gòu)發(fā)展造成的部分葉片通道阻塞促使了這一轉(zhuǎn)化過程。

在本次模擬中，旋轉(zhuǎn)失速在特性曲線“正斜率”處充分發(fā)展，因此，進一步明確了旋轉(zhuǎn)失速與“正斜率”產(chǎn)生的關(guān)系，對削弱機組“S特性”有重要參考價值。