混流式水輪機偏工況運行的大渦模擬方法驗證

蘇文濤，鄭智穎，李小斌，李鳳臣，蘭朝鳳，2

(1.哈爾濱工業大學能源科學與工程學院，150001哈爾濱;2.哈爾濱理工大學電氣與電子工程學院，150080哈爾濱)

混流式水輪機在偏工況運行時，通常會出現內部流動不穩定，機組振動和噪聲增大現象，除了固有的機械因素外，主要是由于流致振動造成，其中包括30%～60%額定負荷時尾水管渦帶誘發大幅壓力脈動，50%額定負荷時輪葉片進口Kármán渦脫流，這樣不僅降低水輪機運行效率，嚴重時還會造成轉輪相鄰葉片中間區域出現葉道渦，并伴隨空化現象［1］.偏工況運行中，在轉輪區主要出現泡狀或片狀的翼型空化，空泡潰滅即為轉輪破壞的主要因素.當轉輪中至少有3個葉片表面出現氣泡時，表明已經發生空化，轉輪中空化程度主要由空化數σ來描述.實驗證明，水輪機壓力脈動的強度隨空化數的減小而迅速提高，而運行效率則隨之下降［2］.針對混流式水輪機中內部的不穩定流動和空化現象，眾多學者展開了實驗研究和數值模擬研究.實驗研究主要針對水輪機流道的各個部分，采用粒子圖像測速儀(PIV)來進行流場的測量和分析.Iliescu等［3］采用 PIV技術對混流式水輪機尾水管內汽、液兩相流動進行測量，研究空化尾水管渦帶結構.Palafox 等［4］和 Uzol等［5］分別采用PIV實驗研究了軸流式水輪機葉頂間隙內的湍流流動，該軸流式水輪機比轉速較低，且葉頂泄漏較為復雜，在下游將形成二次流.Katz等［6-7］基于對軸流式水輪機的PIV測量結果，引入湍動能組分平衡的觀點，分析了葉頂間隙流動及其尾跡的流動特征.針對偏工況運行時的混流式水輪機尾水管內流動，李丹等［8］、王軍等［9］和陳次昌等［10］對其二維流場和渦帶運動特征展開了研究，分析渦帶的周期性運動特征.實驗研究已經揭示了水輪機內復雜流場的形成原因及其與流場邊界之間的作用結果;然而，由于三維流動伴隨空化這一現象的復雜性，目前在眾多發生流動空化的場合中，實驗研究仍存在很多困難.如導葉之間與轉輪葉片之間的流場測試十分復雜，對于具有復雜三維流場的水輪機，采用數值模擬方法研究空化非定常流動成為觀測空化現象發展的重要技術.

對大型混流式水輪機機組內部空化流動數值研究的主要目的是從宏觀和微觀兩方面揭示流動穩定性問題.目前，對于葉片式流體機械內部流場的數值計算，仍以雷諾平均數值模擬(RANS)方法為主，而該方法中雷諾應力基于線性模型計算，使某些包含細節信息的高頻分量被濾除，數值模擬不能體現真實的多尺度湍流流動.大渦模擬(LES)的思想是將小尺度的渦從流動中分離出來，通過小尺度渦和大尺度渦之間的相互影響來描述不同尺度渦結構的相互作用.數學表達上，可在大渦流場的運動控制方程中添加相應的附加應力項來體現小渦的影響，該附加應力項即為亞格子(SGS)應力項.一般認為湍流中小尺度渦的結構趨于各向同性，基于不同渦特征即可構造不同的SGS模型.在水輪機數值模擬方面，Zhang等［11］對某一混流式水輪機進行了三維LES，并將計算結果與壓感片實驗獲得的結果作了比較，得到了令人滿意的結果，但其研究僅限于葉片部位;Liu等［12］采用LES研究了混流式水輪機尾水管渦帶的形態，通過對比不同負荷下的渦帶形態，發現越偏離設計工況，漩渦強度越大;黃劍峰等［13］使用LES動態亞格子湍流模型和先進的滑移網格技術對混流式水輪機全流道動靜干涉的三維非定常湍流場進行了動態LES，計算結果表明LES對復雜流場的仿真效果更加準確.針對上游尾跡對入口流動的擾動，Wang等［14］對某一混流式水輪機的單流道模型進行了LES模擬，發現強尾跡將造成葉片前緣處的大面積流動分離，并沿展向形成主渦.Ma等［15］采用混流式水輪機單流道模型，驗證了LES中SGS模型對計算結果的影響，并和實驗進行了對比，發現適應壁面的渦黏模型和實驗結果符合最好.而最近，RANS-LES混合湍流模型也被逐漸應用在水輪機模擬中［16-17］，較好的模擬結果體現了LES自身特性帶來的優勢.

由于水輪機全流道的幾何復雜性，上、下游部件之間流動干擾，以及旋轉部件與靜止部件之間的流動干涉，使全流道計算及分析相對困難;另外，由于空化的產生是由于液體在低壓區發生相變，如何準確地捕捉到空化區域也成為了必要問題.目前，對于水輪機的數值研究，很少進行全特性驗證.本文用LES方法，并結合空化模型，對某一大型混流式水輪機在不同活動導葉開度下的特性及在不同工況點的空化流場進行全通道三維數值計算，對比分析單相流與空化兩相流計算結果，并對外特性和內流場進行實驗結果對比，驗證針對偏工況運行的混流式水輪機的基于空化模型的大渦模擬方法.

1 計算模型與數值模擬方法

1.1 水輪機模型與幾何參數

對某大型混流式水輪機的模型進行數值模擬，幾何模型如圖1所示.從流動入口到流動出口依次分為蝸殼、固定導葉、活動導葉、轉輪及尾水管5個通道部件，其中入口段前面增加了直管段部分，用來穩定蝸殼的入口流動狀態.該模型的具體幾何參數參見表1.

表1 模型水輪機基本幾何參數

圖1 混流式水輪機幾何模型

1.2 控制方程

采用LES方法(Smagorinsky-Lilly模型)對混流式水輪機內部流動進行數值模擬計算.對不可壓縮牛頓流體流動，用過濾算子對Navier-Stokes方程進行變量過濾，可得基于LES的動量守恒方程為

式中:ui為i方向流動速度，xi為笛卡爾坐標系下x、y、z3 個分量，p為壓力，ρ為流體的密度，υ為流體的運動黏度，τij為SGS應力項.

采用Smagorinsky-Lilly模型有

其中CS為Smagorinsky系數，Δ為LES過濾尺度，Sij為可解尺度的變形率張量，δij為克羅內克符號.Smagorinsky-Lilly模型中，SGS應力可表述為

式中，渦黏性μt為亞格子尺度的湍動黏度.

Smagorinsky常數CS可根據Kolmogorov常數CK獲得:

通常為減小SGS應力的擴散，提高計算過程收斂性，特別是在近壁區，CS取值可適當減小.考慮到全通道中壁面對流動的強烈作用，在近壁區可以按照 Van Driest模型來修正CS值［18］，即

其中y+為網格節點到壁面的量綱一的距離，參數A+=25，CS0=0.1.

為更準確模擬空化流動的發生情況，計算中除使用單相流動模型之外，還引入了描述空泡運動的空化模型.商業CFD軟件ANSYS中目前已經集成了基于Rayleigh-Plesset方程的幾種空化模型，如 Singhal模型、Schnerr-Sauer模型和 Zwart-Gerber-Belamri(ZGB)模型［19］，這幾種模型均描述了空化過程中的相變，考慮了流場中空泡的密度變化與其運動特征的聯系.薛瑞等［20］和劉厚林等［21］針對不同的湍流場，對3種空化模型進行了比較研究，發現ZGB模型的計算收斂性較佳，計算穩定.故文中計算選用ZGB空化模型.

1.3 網格劃分與邊界條件

根據水輪機的幾何模型，依次按照蝸殼、固定導葉、活動導葉、轉輪、尾水管5個部分分別劃分網格，然后組合為模型水輪機的全流道網格，其中在相鄰部件之間，均設置了網格交界面(interface).網格劃分的最終結果如圖2所示，針對具有翼型結構的導葉和轉輪部件，在其近壁面處均進行了邊界層網格加密，以滿足壁面流動計算的要求.

圖2 網格劃分細節圖

對計算結果進行網格無關性驗證，選定設計工況對不同節點數的網格進行了模擬驗證，通常無關性驗證以模型實驗效率為標準，以核對數值計算得到的水力效率.該計算工況下，活動導葉開度a=18 mm、單位轉速n11=60.6 r/min(即最優工況)，模型實驗效率為η=94.8%.計算發現，在網格數超過一定數量后，隨網格節點數目的增加，計算得到的水輪機水力效率值趨于穩定，并和模型實驗的效率相吻合，故當計算域網格數量達到829萬時即可滿足要求.該工況下葉片表面的y+值見圖2，可以看出，y+值符合大渦模擬的計算要求.水輪機各個部件的網格參數如表2所示.

表2 計算域網格結點數及單元數

邊界條件設置中，流動入口設置為計算域進口(蝸殼前直管段);采用總壓進口條件，該值通過模型實驗中所采用的水頭和進口水流速度來給定;流動出口設置為尾水管的出口，采用壓力出口;體壁面均采用無滑移邊界條件.

1.4 求解方法

選取ANSYS中的CFX求解器進行水輪機內部全流道三維模擬計算.針對旋轉部件(轉輪)與靜止部件(活動導葉出口、尾水管進口)的交界面，坐標變換模型均采用瞬態動靜干涉(Transient Rotor-Stator)模式.對流項離散格式選用 High resolution格式，同時計算的每個時間步長中保證轉輪旋轉1°，計算迭代100次，收斂的殘差標準為1×10-5.考慮到流場中近壁區的局部低雷諾數性質，在LES求解中，近壁區的流動采用標準的壁面函數來處理.

2 計算結果及分析

數值模擬中，選取與Su等［22］的實驗測量相同的流動工況，即在偏離最優導葉開度a=18 mm的兩條活動導葉開度線a=10 mm和a=14 mm上，分別選取5個工況點進行LES(單相流和空化流動)計算，具體計算工況分布如表3所示.

表3 計算工況點參數

2.1 內流場和外特性驗證

為對比LES與RANS對混流式水輪機內部流場和空化性能預測的精度，Su等［23］對某一偏工況下的流場進行了三維全通道單相流LES，結果表明，LES比RANS能很好地再現流場中的空化渦.為了描述水輪機內部流場的特征，下面給出對活動導葉開度14 mm時單位轉速為65.4 r/min、單位流量為0.346 7 m3/s的計算結果，計算啟用空化模型，即為多相流LES計算，并將結果和未啟用的空化模型的單相LES進行比較.

圖3為兩種計算方法對無葉區附近流場的預測情況，選取沿導葉高度方向的中間平面，并以靠近固定導葉側的測量線端點為x軸起始點，圖中給出了基于平均速度場的流線及該測量線上量綱一的速度U/U0分布，其中U為當地絕對速度，U0為根據蝸殼進口流量計算得到的平均速度.

圖3 無葉區流場數值模擬與實驗結果對比(a=14 mm)

定量分析表明，活動導葉之后的尾跡流動在導葉壓力面一側(轉輪側)流速較高，而在吸力面一側(固定導葉側)流速較低，這是翼型擾流的顯著特征.通過對比文獻［22］中PIV實驗的流場分布，發現兩種方法獲得的活動導葉尾跡平均流場流線走勢與PIV實驗結果一致，LES空化計算獲得的測量線上速度值的大小與PIV實驗獲得的結果符合得更好.

圖4給出了不同活動導葉開度下LES計算結果與實驗結果的對比，包含了啟用和未啟用空化模型時的情況.可見，單位轉速與單位流量的關系與模型實驗符合很好，表明采用空化計算的LES能夠很好地預測水輪機模型的外特性，圖4(b)中活動導葉開度a=18 mm的計算結果作為對比.從圖4(a)和(b)可以看出，采用空化兩相LES計算的結果與實驗更加吻合，而單相LES的計算結果與實驗值有些許偏離.在上述工況下，同一開度下的水輪機單位流速與單位流量曲線接近線性變化關系，且3個開度下的特性曲線基本平行，即該范圍內，水輪機運行時的實際轉速與全流道內流量也應成線性關系.

圖4 不同活動導葉開度下的外特性對比

圖5給出了活動導葉開度a=10 mm時通過空化LES空化計算水輪機中真實流量-轉速(Q～N)與效率 -轉速(η-N)的關系曲線，其中a=18 mm的計算結果作為對比.由于在水頭不變的情況下，流動中n11與N成線性變化關系，圖5的結果也同時驗證了圖4的結論，即在上述工況變化范圍中，n11～Q11與Q～N均呈現線性關系.額定工況下，計算得混流式水輪機的水力效率接近95%，且水力效率隨著轉輪轉速的降低而減小.當轉速N＜800 r/min時，η迅速下降.在小開度a=10 mm工況，水力效率下降較小，但運行效率總體較低，約87%左右.

為了驗證空化LES對運行效率計算的準確程度，圖6給出了效率曲線與模型水輪機實驗獲得綜合特性曲線的對比.圖中數據點來自于圖5，且將實際轉輪轉速轉換成了單位轉速.

圖6中，計算得到的效率點(η，n11)由圓圈標示，圓圈中心標示了效率點的位置.可以看出，效率點均非常接近于模型實驗的開度線a=10 mm和a=18 mm，即空化計算的效率曲線與綜合特性曲線相符合.另外，對比模型綜合特性曲線可知，設計工況處效率曲線陡峭，而偏工況時，效率曲線比較平緩.可見，啟用空化模型的LES可以很好地再現并預測水輪機的外特性.

圖5 不同活動導葉開度下空化兩相流計算的轉速－流量與轉速－效率曲線

圖6 使用空化LES計算得到不同開度下效率曲線與綜合特性曲線對比圖

2.2 葉道渦預測及其演化

葉道渦是轉輪內部發生空化之后形成的渦結構，其潰滅時會對葉片表面產生嚴重的侵蝕，所以準確預測葉道渦的生成與發展對水輪機正常運轉具有重要意義.為了對葉道渦形態進行可視化，ANSYS軟件提供了許多渦結構特征提取方法，如螺旋度(Helicity)方法、渦流強度(Swirling strength)方法、渦量(Vorticity)方法、Q方法和λ2方法等.

用渦量準則顯示轉輪區的渦量分布.通過給定渦量的閾值，渦量準則可以正確判定出漩渦存在的區域，下面分別對啟用和不啟用空化模型時的工況進行分析.針對活動導葉開度a=10 mm、n11=59.9 r/min工況，先分別計算出空化LES和單相LES模擬結果的平均速度場，再采用0.5倍平均渦量 ＜|Ω|＞作為閾值進行葉道渦顯示，如圖7所示.可以看出，空化LES可以捕捉到更多的葉道渦，且更加符合模型實驗獲得的葉道渦形態［22］;但單相LES則不能很好地捕捉到沿流向的葉道渦運動.

圖7 平均流場的葉道渦可視化

針對葉道渦的演化過程，下面將進一步對渦結構的變化進行說明.如圖8所示，計算工況為活動導葉開度14 mm，單位轉速61.0 r/min，進行空化LES計算.

圖8 空化LES計算獲得的葉道渦演化過程

圖8中給出了水輪機轉輪內部葉道渦生成、發展和潰滅的一個典型演化過程，空化渦使用Q方法進行提取，并和相同工況下的模型實驗進行了比較.可見，葉道渦位于葉片通道之間，并初生于葉片進水邊附近，隨時間變化，轉輪順時針旋轉，空化渦體積逐漸增長(1～3)，達到最大之后，又逐漸減小直到消失(4～5).需要注意的是，從1～5，在葉片的出水邊(模型實驗未監測葉片出水邊)附近，有大量空化渦沿著展向分布，這是由于葉道渦在不斷生成積累過程中，隨主流流動而被帶到了出水邊.部分空化渦將被進一步帶到流動通道的下游，在尾水管中形成碎渦結構.

文獻［22］已通過實驗發現，轉輪在偏工況下運行時，其中葉道渦的變化呈現高頻特征，且其演化頻率為轉頻的整數倍，該倍數一般為轉輪葉片數.在上述計算中，設置時間步長為轉輪每轉動1°的時間，可以得到葉道渦變化的一個周期為轉輪轉動 24°，對應時間為 0.005 2 s，約 192.31 Hz.而該流動條件下，單位轉速為61.0 r/min，轉輪實際轉速 757.8 r/min，故轉頻為 12.63 Hz，可見葉道渦演化的周期正好約為轉頻的15倍，和葉片數相同.于是，上述數值模擬的結果和模型實驗的結果相符合，即使用空化LES計算能夠很好地預測水輪機偏工況運行時的葉道渦運動.

2.3 尾水管渦帶演化

尾水管渦帶是水輪機中另一種空化現象，從轉輪出口的中心開始，向肘管發展.由于轉輪出口水流的周期性運動，該渦帶也在不斷旋轉運動，隨不同工況的變化渦帶具有不同的形態，其運動會引發尾水管壁面的周期性壓力脈動.

針對活動葉開度10 mm、單位轉速64.9 r/min工況分別使用空化LES和單相LES計算，提取尾水管中空化渦帶結構，如圖9所示.將25℃下水的汽化壓力作為等值面提取渦形態，可見，空化LES計算得到的渦帶形態與模型實驗觀測到的結果接近，該工況下尾水管渦帶呈現螺旋狀.另外在該工況下單相LES出現了分岔的渦帶形態，說明單相流計算不能很好地處理氣液密度變化和壓力變化的關系，而空化兩相流則考慮了空泡密度與流場運動的關系，故能獲得較為滿意的空化形態.

為了獲取尾水渦帶的演化過程，對活動導葉開度14 mm、單位轉速61.0 r/min的工況進行空化LES計算，如圖10所示.圖中顯示為從尾水管開始的錐管段渦帶部分，事實上渦帶發展是從轉輪下方的泄水錐開始，為了顯示清楚渦帶的旋轉情況，特意截取尾水管中的一段空化渦帶.圖10表征了尾水管渦帶旋轉一周的演化過程，其所用時間約為轉輪旋轉2.5周的時間，可見，渦帶的演化是一個低頻過程，即約為轉頻的0.4倍，和文獻［22］的結論相符合.

圖9 尾水管渦帶形態對比

圖10 空化LES計算獲得的空化渦帶演化過程

由形態結構可看出，從泄水錐生成的空化渦被攜帶到下游時，會隨著轉輪的旋轉而變成螺旋狀結構，并不斷增長，獲得的空化渦結構和模型實驗比較符合.而中心渦帶在增長的同時，在接近尾水管壁面的部分有破碎渦生成，可以推測，部分破碎渦來源于轉輪出水邊集結的空化渦.這些破碎渦在渦帶旋轉一周的時間中，將從下游的肘管處瀉出，之后，渦帶形態將開始下一個周期變化.從機理上看，由于空化LES計算特別考慮了空化區域的壓力、密度變化和流場運動的關系，故能更好地處理流場中壓力梯度變化劇烈的區域.可見通過LES可以預測尾水管入口處壓力梯度較高的區域，和實際流動情況相符合.

綜上所述，采用空化模型的LES不僅可以準確地獲得混流式水輪機內部三維流場，更能很好地捕捉到空化區域的形態特征，尤其針對葉道渦和尾水管渦帶的形態演化，對水輪機的穩定運行具有重要的意義.

3 結論

1)LES可以預測尾水管入口存在壓力梯度較高的區域，且能捕捉到尾水管錐管段復雜的渦帶旋轉流動，而RANS則不能準確預測.

2)采用空化兩相流LES計算的外特性結果與實驗對比吻合度高，而單相流計算結果與實驗值偏差較大.

3)采用空化兩相流的LES可以很好地捕捉到轉輪葉片之間沿流向發展的葉道渦，與實驗結果相吻合，驗證了基于空化模型的LES方法在預測偏工況下混流式水輪機性能的可行性與準確性.

4)空化LES獲得的葉道渦演化過程為高頻運動，而尾水管渦帶為低頻運動，和模型實驗結果相符合.

［1］D?RFLER P，SICK M，COUTU A.Flow-induced pulsation and vibration in hydroelectric machinery［M］.London:Springer，2013.

［2］WEI X Z，SU W T，LI X B，et al.Effect of blade perforation on Francis hydro-turbine cavitation characteristics ［J］.Journal of Hydraulic Research，2014，52(3):412-420.

［3］ILIESCU M S，CIOCAN G D，AVELLAN F.Analysis of the cavitating draft tube vortex in a francis turbine using particle image velocimetry measurements in two-phase flow［J］.ASME Journal of Fluid Engineering，2008，130:21105.

［4］PALAFOX P，OIDFIELD M L G，LAGRAFF J E，et al.PIV maps of tip leakage and secondary flow fields on a low-speed turbine blade cascade with moving end wall［J］.Physics of Fluids，2008，130:11001.

［5］UZOL O，CHOW Y C，KATZ J，et al.Average passage flow field and deterministic stresses in the tip and hub regions of a multistage turbomachine［J］.Experiments in Fluids，2003，125(4):714-725.

［6］ UZOL O，KATZ J.Flow measurement techniques in turbomachinery［M］//Handbook of Experimental Fluids Mechanics.Heidelberg:Springer，2007:917-957.

［7］KATZ J，CHOW Y C，SORANNA F，et al.Experimental Characterization of Flow Structure and Turbulence in Turbomachines［C］//5th International Symposium on Fluid Machinery and Fluids Engineering.Jeju:［s.n.］，2012.

［8］李丹，陳次昌，季全凱，等.PIV測試水輪機尾水管錐管流場［J］.西華大學學報:自然科學版，2004(S1):228-230.

［9］王軍，付之躍，溫國珍，等.HL220改型裝置不同工況下尾水管進口流動狀態PIV試驗研究［J］.工程熱物理學報，2005，26(S1):93-96.

［10］陳次昌，李丹，季全凱，等.混流式水輪機尾水管內部流場的 PIV測試［J］.機械工程學報，2006，42(12):83-88.

［11］ZHANG L X，GUO Y K，WANG W Q.Large eddy simulation of turbulent flow in atrue 3D francis hydroturbine passage with dynamicalfluid-structure interaction［J］.International Journal for Numerical Methods in Fluids，2007，54(5):517-541.

［12］LIU X B，ZENG Y Z，CAO S Y.Numerical prediction of vortex flow in hydraulic turbine draft tube for LES［J］.Journal of Hydrodynamics，Series B，2005，17(4):448-454.

［13］黃劍峰，張立翔，王文全.混流式水輪機全流道三維非定常湍流場的動態大渦模擬［J］.中國農村水利水電，2009(1):101-108.

［14］WANG W Q，ZHANG L X，YAN Y.Large-eddy simulation of turbulent flow considering inflow wakes in a Francis turbine blade passage［J］.Journal of Hydrodynamics，Series B，2007，19(2):201-209.

［15］MA J M，WANG F J，TANG X L.Comparison of several subgrid-scale models for large-eddy simulation of turbulentflows in water turbine［C］//The 4thInternational Symposium on Fluid Machinery and Fluid Engineering.Beijing:［s.n.］，2008.

［16］NETO A D A，JESTER-ZUERKER J，JUNG A，et al.Evaluation of a Francis turbine draft tube flow at part load using hybrid RANS-LES turbulence modelling［J］.IOP Conf Series:Earth and Environmental Science，2012(15):062010.

［17］SENTYABOV A V，GAVRILOV A A，DEKTEREV A A，et al.Numerical investigation of the vortex core precession in a model hydro turbine with the aid of hybrid methods for computation of turbulent flows［J］.Thermophysics and Aeromechanics，2014， 21(6):707-718.

［18］王福軍.計算流體動力學分析:CFD軟件原理與應用［M］.北京:清華大學出版社，2004.

［19］ZWART P J，GERBER A G，BELAMRI T.A twophase flow model for predicting cavitation dynamics［C］//Fifth International Conference on Multiphase Flow.Yokohama:ICMF COP，2004:No.152.

［20］薛瑞，張淼，許戰軍，等.對不同空化模型的比較研究［J］.西北水電，2014(2):85-89.

［21］劉厚林，劉東喜，王勇，等.三種空化模型在離心泵空化流計算中的應用評價［J］.農業工程學報，2012，28(16):54-59.

［22］SU W T，LI X B，LI F C，et al.Experimental investigation on the characteristics of hydrodynamic stabilities in Francis hydro-turbine models［J］.Advances in Mechanical Engineering，2014，Article ID 486821.

［23］SU W T，LI F C，LI X B，et al.Assessment of LES performance in simulating complex 3D flows in turbomachines ［J］. Engineering Applications of ComputationalFluid Mechanics， 2012， 6(3):355-364.