燈泡貫流式水輪機尾水管流動特性分析

敏 政, 李秋桐, 李琪飛, 張毅鵬, 李正貴，鄧育軒(. 蘭州理工大學 能源與動力工程學院，蘭州 730050； . 南昌工程學院 機械與電氣學院，南昌 330099)

0 引 言

貫流式水輪機是適用于低水頭，大流量的水電機組[1]，由于河流上游流量的不穩(wěn)定，機組需要經(jīng)常改變運行工況[2]，特別是在非協(xié)聯(lián)工況下，水輪機內(nèi)部流態(tài)會急劇惡化，對過流部件的空化性能會產(chǎn)生很大的影響[3]。葉片和導葉之間漩渦與尾水管渦帶是引起壓力脈動的重要原因之一，其中尾水管渦帶所引起的低頻振動引起眾多學者興趣。當水輪機在偏離于最優(yōu)工況點運行時，尾水管中的會產(chǎn)生空腔渦帶，且其渦帶形態(tài)與此時水輪機負荷的大小有著密切的關系，在低負荷時出現(xiàn)螺旋狀偏心渦帶，高負荷時產(chǎn)生柱狀渦帶。尾水管內(nèi)流場的特性會因為空化渦帶的形態(tài)變化而隨之改變。楊靜[4]采用LES模型對混流式水輪機尾水管內(nèi)渦帶引起的壓力脈動進行了研究。夏林生[5]等對燈泡貫流式水輪機在飛逸狀態(tài)下進行了全流道的數(shù)值模擬，計算了此過渡過程中尾水管內(nèi)渦帶的演變規(guī)律。季斌等[6]通過計算小流量工況下過流部件之間的動靜干涉作用以及渦帶的演化過程。錢忠東等[7]采用大渦模擬湍流模型，對燈泡貫流式水輪機額定工況下和偏工況下進行了計算，分析不同流量下尾水管內(nèi)產(chǎn)生的渦帶對壓力脈動的影響，并與實驗數(shù)據(jù)進行了比較。王磊等[8]基于氣泡兩相流的方法，對白鶴灘模型水輪機進行定長和非定常計算，預估了水輪機轉(zhuǎn)輪和尾水管內(nèi)空化的發(fā)生情況，并與實驗結(jié)果進行了比較。由于燈泡貫流式水輪機在國內(nèi)研究起步較晚，因此相關的一些研究多見于國外[9-11]。

本文采用汽液混輸相均質(zhì)假設和基于組分輸運方程的空化模型，通過求解汽液混相均質(zhì)流的雷諾平均N-S方程以及汽相組分輸運方程，并考慮不可凝結(jié)汽相的影響，計算了貫流式水輪機尾水管內(nèi)部的空化流場和壓力脈動變化特點。

1 基本參數(shù)和計算工況

1.1 基本幾何參數(shù)

根據(jù)提供的基本參數(shù)，以Pro/e為三維建模軟件，圖1為建立從進口到尾水管出口的全流道幾何模型，圖2為轉(zhuǎn)輪的幾何模型，表1列出了水輪機基本參數(shù)。

圖1 水輪機全流道三維實體圖Fig.1 The geometrical model of tubular-turbine

圖2 葉片示意圖Fig.2 The sketch of blades

表1 水輪機主要參數(shù)Tab.1 Main characteristics of the tubular turbine

1.2 網(wǎng)格參數(shù)

對幾何模型進行網(wǎng)格劃分，采用ICEM網(wǎng)格劃分軟件對進口區(qū)域，尾水管和流道進口劃分為六面體結(jié)構(gòu)化網(wǎng)格，對轉(zhuǎn)輪和導葉劃分為四面體非結(jié)構(gòu)化網(wǎng)格。為了減小網(wǎng)格數(shù)對計算結(jié)果的影響，對設計工況點進行網(wǎng)格無關性假設，最終確定的網(wǎng)格數(shù)為10 065 610，節(jié)點數(shù)為2 788 414，計算所用的網(wǎng)格如圖3所示，網(wǎng)格無關性假設，如圖4所示。

圖3 水輪機全流道網(wǎng)格圖Fig.3 Mesh of turbine

圖4 網(wǎng)格無關性假設Fig.4 Verifying grid independence

2 數(shù)值計算方法

2.1 數(shù)學模型

由Rayleigh-Plesset方程得到的空泡特性和考慮各種影響的空化模型,水輪機中空化流動計算的連續(xù)性控制方程描述如下。

(1)混合流體的連續(xù)性方程：

(ρ)+▽·(ρv)=0

(2)汽相的連續(xù)性方程：

(ρf)+▽·(ρfv)=Re-Rc

蒸汽生成率 ：

蒸汽凝結(jié)率：

(3)混合流體的動量方程：

▽ (ρvv)=

式中：ρ為水流相和空泡相形成的混合流體質(zhì)量密度；v為混合流體的速度矢量；f為空泡相的質(zhì)量組分；Re為水蒸氣的生成率；Rc為水蒸氣的凝結(jié)率；P為靜壓力；μ為分子黏性系數(shù)；μt為湍流黏性系數(shù)；g為重力加速度矢量。

由于采用了混合流體處理水輪機內(nèi)部的空化流動，因此湍流模型中的和在形式上與單相流動相同，但其中的混合變量均為混合流體的平均量。計算中的物性參數(shù)取工作水溫為25時的值，根據(jù)水輪機的實際運行情況，其中汽化壓力=3 540，水-空泡表面張力-0.071 7 NM，假定水中不凝結(jié)性氣體的質(zhì)量組分為1×10-5。

2.2 離散和湍流模型

對計算區(qū)域采用混合網(wǎng)格，采用有限體積法對上述方程進行離散求解。時間項采用全隱世離散，擴散項和離散相采用中心差分格式，動量方程的對流項采用二階迎風差分格式，變量儲存在控制體中心，壓力和速度采用SIMPLEC算法進行耦合求解。由于水輪機在偏工況運行時，存在漩渦，二次流等復雜流態(tài)，本文定常計算使用RNG k-ε模型代替標準的k-ε模型，非定常計算采用SST湍流模型，空化模型選用基于輸運方程模型的Singhal完全空化模型[9,10]。

2.3 邊界條件

為了模擬真實情況下水輪機內(nèi)部的流動情況，在計算流動區(qū)域進口面上根據(jù)流量給定速度條件，并假定速度垂直于進口面；出口邊界條件：自由出流；壁面邊界條件：采用無滑移邊界條件。初始流場條件：設置導葉和槳葉開度為某一固定值，通過三維定長計算得到流場結(jié)果，以此作為非定常計算的初始結(jié)果。為了預測尾水管的壓力脈動，全文進行全流場的非定常湍流計算，考慮后期進行頻譜分析的需要，計算的時間步長定位轉(zhuǎn)輪旋轉(zhuǎn)一圈的1/360，轉(zhuǎn)速n=68.18 r/min，得到非定常的時間步長為0.002 444 5 s。

3 計算監(jiān)測點位置與試驗

3.1 計算監(jiān)測點

從導葉段到尾水管一共取了6個截面，截面A位于導葉出口處，截面B位于轉(zhuǎn)輪出口處，截面C位于輪轂尾部，截面D為尾水管進口處，截面E，截面F位于尾水管中后部，在每個截面上均勻取四個監(jiān)測點，共有24個監(jiān)測點。如圖5所示。

圖5 監(jiān)測點位置分布圖Fig.5 Positions of sampling points

3.2 模型試驗

為了保證幾何模型建立的準確性，對所計算的結(jié)果與模型試驗進行了對比。本次的模型試驗由東方電機股份有限公司承擔，如圖6所示，試驗臺參數(shù)均滿足ICE所規(guī)定的要求[11]，圖7為在額定水頭下數(shù)值計算和模型試驗的效率變化曲線。

圖6 試驗現(xiàn)場Fig.6 Test site

通過圖7線可以看出，效率的變化趨勢基本一致。在較大空化數(shù)下，空化并沒有對水輪機的外特性產(chǎn)生明顯影響，效率基本保持不變。然而當空化數(shù)降低至初生空化附近時，效率曲線出現(xiàn)小幅度下降，其原因可能是由于少量的空化氣泡出現(xiàn)在葉片出口邊附近時，轉(zhuǎn)輪出口速度矢量的方向發(fā)生了改變，尾水管進口的軸向速度有所增加，因此效率略微下降。

圖7 效率變化曲線Fig.7 The curves of efficiency

隨著出口壓力的進一步降低，較大范圍的空化氣泡產(chǎn)生在葉片區(qū)域及尾水管內(nèi)，嚴重阻塞流道致使流量及輸出功率降低，進而導致效率曲線明顯下降，效率下降3%時的空化點即為臨界空化點。然而受到計算精度及非定常計算的數(shù)值模型等因素的影響，數(shù)值計算的結(jié)果略高于實驗得到的結(jié)果。

4 計算結(jié)果與分析

4.1 定常計算結(jié)果

圖8為在不同空化系數(shù)下流場內(nèi)壓力和流線分布。從圖8中可以看出，隨著空化系數(shù)的降低，流道內(nèi)的流線密集程度逐漸在降低，說明流道內(nèi)的流體流通性較好，沒有發(fā)生邊界層分離，二次流和漩渦等干擾正常流動的因素較少，這時流場比較穩(wěn)定，水力損失少。當空化系數(shù)繼續(xù)降低后，轉(zhuǎn)輪出口處壓力持續(xù)下降，在尾水管內(nèi)出現(xiàn)的“柱狀”低壓區(qū)面積增大，此時流道內(nèi)空化發(fā)生加劇，含有氣體的空泡附著在葉片表面，降低了葉輪的做功能力，并且氣泡充滿流道阻礙了水流的正常流動，使流道內(nèi)流場結(jié)構(gòu)產(chǎn)生了巨大的改變，致使水輪機的效率繼續(xù)降低，這與圖7中所顯示的效率變化是相符的。通過定常計算，數(shù)值模擬的效率和模型試驗的結(jié)果相差不大，說明了本次計算的準確性和數(shù)值計算的可靠性。

圖8 不同空化系數(shù)下流道內(nèi)壓力和流線分布Fig.8 Pressure and streamline contour in different cavitation coefficient

圖9 尾水管渦帶形態(tài)Fig.9 Pattern of vortex rope in draft tube

4.2 非定常計算結(jié)果

圖9為在空化系數(shù)σ=1.2尾水管內(nèi)渦帶隨時間的變化圖。從圖9中可以看出，在水輪機尾水管內(nèi)的空化渦帶形態(tài)呈現(xiàn)為接近“柱狀”，渦帶長度以及直徑均隨著時間的變化而演變。在轉(zhuǎn)輪轉(zhuǎn)動一個周期剛開始時，在尾水管內(nèi)出現(xiàn)體積較小的柱狀空化渦帶。隨著時間的向后推移，空化區(qū)域逐漸增大，空化渦帶體積變化劇烈，在渦帶體積增大的過程中，尾部會出現(xiàn)脫流；當?shù)竭_0.6t時，渦帶的尾部直徑會出現(xiàn)很明顯的突然減小，這種現(xiàn)象稱之為“環(huán)形水躍”，Dorfler[12]在使用兩相流對空化渦帶進行數(shù)值模擬中也發(fā)現(xiàn)了這個現(xiàn)象，并且這種現(xiàn)象在單相流的模擬中不會出現(xiàn)。這與Alligne S[13]等所研究的現(xiàn)象是一致的。當轉(zhuǎn)輪周期繼續(xù)增加時，尾水管內(nèi)的空化渦帶直徑及長度仍會繼續(xù)增大，此時尾水管內(nèi)空化已經(jīng)發(fā)展至非常嚴重的程度，大量的氣泡會堵塞流道。環(huán)形水躍現(xiàn)象一般解釋為：柱狀空化渦帶主要通過體積波動來影響流場，但空化渦帶本身也是在不斷旋轉(zhuǎn)的，受尾水管擴散段及下游流場的影響，空化渦帶的尾部附近流場壓力會越來越高，空泡體積濃度越來越低，因而無法繼續(xù)形成直徑較大的空化渦帶，所以會出現(xiàn)突然的水躍。隨著空化渦帶體積波動過程，渦帶環(huán)形水躍以下直徑較小的部分會出現(xiàn)頻繁的空泡增長和脫落現(xiàn)象。

4.3 壓力脈動計算結(jié)果分析

圖10為圖5中所示截面中監(jiān)測點1處的壓力脈動時域圖。由圖10中可以看出，所檢測截面上壓力的相位和脈動周期均表現(xiàn)出了相似的變化規(guī)律，但是流體經(jīng)過轉(zhuǎn)輪做功后使這些位置表現(xiàn)出不同的壓力脈動特性。從導葉出口處到尾水管內(nèi)所監(jiān)測的截面上可以看到，在導葉出口處，由于并未發(fā)生空化，所以此處壓力值與其他部分相比較高，壓力波動平穩(wěn)；在轉(zhuǎn)輪出口處因為發(fā)生了空化，此處的壓力波動劇烈，在一個大的波動后出現(xiàn)若干個較小的脈動區(qū)域，隨著時間的向后推移，這種變化會不斷地持續(xù)下去，壓力在數(shù)值上呈現(xiàn)出劇烈下降趨勢，這與流場中所展示的壓力變化是相吻合的；隨著在輪轂末端區(qū)域出現(xiàn)明顯的柱狀空腔渦帶，沿著水流方向尾水管兩側(cè)呈現(xiàn)出對稱性延生分布，且渦帶區(qū)域有著不斷擴大的趨勢，表現(xiàn)出了極強的不穩(wěn)定性，因而壓力進一步降低波動更加劇烈；在尾水管中段，渦帶消失，空化渦帶由于自身的旋轉(zhuǎn)作用的減弱而對內(nèi)部流場的壓力影響同樣愈見下降，此時壓力數(shù)值回升，所以表現(xiàn)為監(jiān)測面E和F上壓力脈動的相位差也愈來越小，波動趨于平緩，這對于尾水管的壓力回復性能提高是有利的。當空化系數(shù)降低時，整個流場內(nèi)壓力波動加劇，并且壓力數(shù)值大幅降低，這是由于空化發(fā)生更加劇烈所導致的。

水輪機在穩(wěn)定運行過程中，機組的轉(zhuǎn)頻為fn=1.123 6 Hz，對計算數(shù)據(jù)進行快速傅里葉變換得到頻域圖。

圖10 監(jiān)測斷面壓力脈動時域圖和頻域圖Fig.10 Frequency spectra of pressure fluctuations at different sections

圖10中右側(cè)為在3個空化系數(shù)下相應截面監(jiān)測點1的壓力頻譜，3個工況下的結(jié)果比較類似，轉(zhuǎn)輪出口處、尾水管內(nèi)轉(zhuǎn)頻為4fn，為葉片通過頻率，此時頻率下壓力脈動的幅值較小；輪轂末端和尾水管進口處的轉(zhuǎn)頻為0.8fn，屬于低頻壓力脈動，這時壓力脈動幅值相比較大，這是因為柱狀空化渦帶工況尾水管進口處流場出現(xiàn)“環(huán)形水躍”現(xiàn)象，導致渦帶尾部受到水流產(chǎn)生的脫流現(xiàn)象而導致空泡不斷發(fā)生潰滅的影響，這時流場內(nèi)部出現(xiàn)十分明顯的流量不均勻性，這就是導致壓力脈動幅值在此處高于其他地方的原因。

不同頻率的壓力脈動同時向著上游和下游同時傳播開來，與此同時幅值不斷降低，說明幅值會隨著距離的增加而降低，最小幅值出現(xiàn)在尾水管中段，這是由于此處渦帶已經(jīng)逐漸消失，受到渦帶影響很小，流體運動平穩(wěn)。從圖10中還可以看出，除了空化渦帶所占據(jù)的主頻以外，流場中還存在一個空化渦帶所引起的頻率為5fn的次頻，從B、C、D 3個截面可以明顯看出其對于流場的影響，這和圖9中出現(xiàn)的小幅動的壓力波動相對應。隨著空化系數(shù)的降低，壓力脈動幅值增大，這是由于流道內(nèi)空化發(fā)生加劇，形成各種復雜的空化漩渦，特別是在輪轂末端和尾水管進口處不斷增強空腔渦帶，誘發(fā)了大量的低頻高幅壓力脈動，易引發(fā)由于空化渦帶產(chǎn)生的噪聲。導葉和轉(zhuǎn)輪之間截面A中監(jiān)測點壓力脈動主頻為fn=4.54 Hz，接近轉(zhuǎn)輪的葉頻，這主要是由于轉(zhuǎn)輪在旋轉(zhuǎn)過程中切割活動導葉所形成的射流-尾跡產(chǎn)生強烈的動靜干涉作用，在此作用下，造成了頻率較高的低幅壓力脈動，因而很容易引起機組的振動。此時壓力脈動的最大幅值ΔH/H達到6.7%，幅值沿著下游方向減小，動靜干涉作用隨著逐漸遠離葉輪而影響逐漸降低。

5 結(jié) 論

(1)隨著時間的變化，空腔渦帶的體積逐漸增大，渦帶的長度增幅遠大于直徑的增幅，渦帶分布均勻，位于尾水管中部；在渦帶的尾部由于不斷出現(xiàn)的小的空泡脫流及潰滅，從而導致了渦帶尾部出現(xiàn)“環(huán)形水躍”現(xiàn)象。

(2)渦帶出現(xiàn)時，尾水管內(nèi)所取的截面監(jiān)測點位置處壓力脈動的相位差并不是完全相同的，尾水管內(nèi)壓力脈動以低頻為主，渦帶主頻為轉(zhuǎn)頻，其余頻率由渦帶公轉(zhuǎn)產(chǎn)生以2～5倍主頻的次頻；從轉(zhuǎn)輪出口到尾水管進口處壓力脈動變化最為劇烈，這與此處產(chǎn)生空腔渦帶有關；隨著空化系數(shù)的逐漸降低，監(jiān)測點處壓力脈動變化更加劇烈。

(3)壓力脈動對空化系數(shù)的變化很敏感，空化系數(shù)越低，壓力脈動越嚴重，會降低機組運行中的穩(wěn)定性，使水輪機的水力效率下降。

□

[1] 程良駿.水輪機[M].北京:機械工業(yè)出版社,1981.

[2] 于 波,肖惠民.水輪機原理與運行[M].北京:中國電力出版社,2008.

[3] 揚 波，李培根.貫流式水輪機轉(zhuǎn)輪試驗研究[J].東方電機,2001,(1):41-47.

[4] 楊 靜, 混流式水輪機尾水管空化流場研究[D].北京:中國農(nóng)業(yè)大學,2013.

[5] 夏林生,程永光，張曉曦.燈泡式水輪機飛逸過渡過程3維CFD模擬[J].四川大學學報(工程科學版),2014,46(5):35-41.

[6] 季 斌,羅先武，西道宏.混流式水輪機渦帶工況下兩級動靜干涉及其壓力脈動傳播特性分析[J].水力發(fā)電學報,2014,31(1):191-196.

[7] 錢忠東,魏 巍,馮曉波.燈泡貫流式水輪機全流道壓力脈動數(shù)值模擬[J].水力發(fā)電學報,2014,33(4):242-249.

[8] 王 磊,婁 瑜，王照福.混流式模型水輪機空化流動分析與試驗研究[J].排灌機械工程學報，2014,32(9):771-775.

[9] Necker J,Aschenbrenner T. Model test and CFD calculation of a cavitating bulb turbine [C]∥ Proceeding of the 25th IAHR Symposium on Hydraulic Machinery and Systems. Timisoara:IPO Publishing,2010.

[10] Yang W, Wu Y L,Liu S H. An optimization method on runner blades in bulb turbine based on CFD analysis [J].Science China Technological Sciences,2011,54(2):338-344.

[11] Thaithacha Sudsuansee, Udomliat Nontakaew. Simulation of leading edge cavitation on bulb turbine [J]. Songklanakarin Journal of Science and Technology, 2011,(1):51-60.

[12] Zwart P,Gerber A G,Belamri T.A two-phase model for predicting cavitation dynamics[C]∥ Proceedings of ICMF 2004 International Conference on Multiphase Flow,Yokohama,Japan,2004:1-11.

[13] Mejri I,Bakir F,Rey R.Comparison of computational results obtained from a homogeneous cavitation model with experimental investigations of three Inducers[J].Journal of Fluids Engineering,2006,128(6):1 308-1 323.

[14] Liu Shuhong,Li Shengkai,Wu Yulin. Pressure fluctuation Prediction of a model Kaplan turbine by unsteady turbulent flow simulation[J]. Journal of Fluid Engineering,2009,131(10):101012.1-101012.9.

[15] Dorfler P K, Keller M,Braun O. Francis full-load surge mechanism identified by unsteady 2-phaseCFD[C]∥ Proceedings of 25th IAHR Symposium on Hydraulic Machinery and Systems, Timisoara,Romania,2010.

[16] Alligne S, Maruzewski P, Dinh T,et al. Prediction of a Francis turbine prototype full load instability from investigations on a reduced scale model[C]∥ Proceedings of the 25th IAHR Symposium on Hydraulic Machinery and Systems,Timisoara,Romania,2010.