基于TBR模型的高水頭混流式水輪機水力性能預測

孫龍剛，郭鵬程，麻 全，鄭小波，羅興锜

（西安理工大學水利水電學院，西安 710048）

0 引 言

動靜干涉（rotor-stator interaction，RSI）現象是葉輪機械的顯著特征，主要是由于轉子與定子的相對運動，前排葉片的尾跡區對后排葉片形成周期性擾動，引起葉輪內部速度及壓力場隨時間發生周期性變化，誘導幅值較高的壓力脈動[1-3]。Franke等澄清了RSI的產生機理，首次對 RSI作用下的壓力波模式進行了可視化研究，并對一個20個活動導葉的水泵水輪機進行了試驗驗證[4]。葉輪機械 RSI效應的捕捉，必須進行瞬態計算，而為了避免共振，葉輪機械轉子與定子葉片數往往互質，使得周期性邊界條件設置失效，此時必須對整圈葉片進行計算，保證轉子與定子之間的螺距比為1，從而確保上下游信息的準確傳遞，以提高計算精度[5]。江偉等[6]研究發現，隨著流量增加，離心泵半高導葉斷面間隙對葉輪和導葉之間的 RSI作用逐漸減弱，而葉輪與蝸殼之間作用逐漸增強。李偉等[7]采用PIV技術對混流泵內部流場進行了試驗研究，測試結果顯示受葉輪和導葉之間 RSI作用的影響，導葉進口處渦旋流動結構較為明顯，造成流道堵塞。麻全[8]數值研究了高水頭混流式水輪機穩定飛逸工況下水輪機不同位置的壓力脈動特性及 RSI效應。葉輪機械全流道數值計算相當耗費計算資源及時間，并且與計算域的離散、湍流模型，求解控制等均有關。水輪機部件相對較多，為更準確地捕捉內部流動特性，全流道計算的網格已達上千萬甚至更多[9]，而進行全流道瞬態計算，單個工況的計算時間通常要達數十天，若采用更高級的湍流模型及要求更高精度，時間會更長且消耗更多的CPU及內存[10-11]。

TBR（transient blade row）模型[12]提供了一種采用1個或者 2個流道計算葉輪機械瞬態效應的方法，并且保證一定的精度，從而解決了瞬態計算部分流道相鄰級之間角度比不為 1的問題，大大降低了葉輪機械瞬態分析的計算資源與時間消耗，使得瞬態分析能夠成為葉輪機械常規設計的有力工具。TBR模型提供了3種方法：PT（profile transform）法，TT（time transform）法和 FT（fourier transform）法[13]。PT法是最直接的變換法，是將分布在交界面上的變量按照角度比進行比例縮放從而保證通量守恒；TT法基于時間傾斜法，是一種顯示計算；而 FT方法采用傅里葉級數對周向邊界和轉靜邊界進行求解歷史的重構[14-15]。Zori 等[16]分別采用TT方法和PT方法對某一1.5級跨聲速壓縮機進行了數值模擬，對比分析不同方法的在氣體動力學參數及內部流動特征的捕捉能力。Connell等[17]分別以高壓動力渦輪級和低壓飛機渦輪級為研究案例，在計算資源的優化利用以及時間精度上對比分析不同參數對PT法、TT法以及FT法的影響，結果表明，TT法與FT法在時間精度上要高于PT法，而FT法在穩定性上要強于TT法。Robinson等[18]利用TT法對某一帶導葉擴壓器的離心壓縮機進行了數值計算，獲得了與試驗值比較一致的數值解。Cornelius等[19]評估了穩態計算方法和瞬態的 TT法在多級軸流壓縮機旋轉失速點捕捉和氣體動力穩定性預測的可靠性，與穩態計算比較，在接近失速區瞬態的 TT方法獲得了與試驗數據更匹配的速度線。此外，TT方法與FT方法還可應用于進口擾動、葉片顫振等。

目前TBR模型多用于渦輪機及壓縮機的性能預估及負荷預測上，且多用于可壓縮流動，而以水為工作介質的葉輪機械具有大密度特性，其動態液體力更大，RSI效應將更為顯著。2014年，挪威科技大學（Norwegian University of Science and Technology-NTNU）公布了名為Francis99（又名 Tokke模型）混流式模型水輪機的幾何模型，同時提供了詳細的水力效率、壓力脈動等試驗數據，便于對數值計算進行驗證、分析及比對[20-21]。Nicolle等[22]采用 TBR模型對 Francis99高水頭水輪機進行了性能預測，該文的計算結果表明壓力脈動幅值沿流線方向依次減小，與試驗結果差異較大。與Nicolle等的計算方法不同，本文采用不同的模擬域并重新生成高質量的網格，分別利用全通道方法、PT方法與FT方法對Francis99模型水輪機進行不可壓縮的RSI及葉片負荷數值預測，并與試驗值結果進行對比，以評估PT方法及FT方法在預測以水為工作介質的不可壓縮葉輪機械性能的可靠性。

1 計算理論及模型

1.1 流動控制方程及k-ω SST模型

非穩態的 RANS（Reynolds Average Navier-Stokes equations）方程如下[23]

式中t為時間，s; ρ為密度，kg/m3；p為壓力，Pa；v為渦黏度；Ui與Uj分別為i和j方向上的時均速度，m/s;ui′和u′j分別為脈動速度，為Reynolds應力張量。

本文采用k-ω SST湍流模型來封閉方程組（1）。

式中vt表示湍動渦黏度，為待求量，Sij為變形率張量，δijk為Kronecker符號。

k-ω SST湍流模型中湍動能k和比耗散率ω的輸運方程如式（3）所示[24]

式中F1與F2為混合函數，Pk為湍動生成項，α，α1，β，β*，σk，σω，σω2等均為方程組閉合系數。

k-ω SST湍流模型在邊界層使用k-ω湍流模型，在其余區域應用k-ε湍流模型，可較好地捕捉葉輪機械的流動分離現象。

1.2 TBR模型

由于TT方法僅僅適用于可壓縮流動，故本文只進行了PT法與FT法的數值計算。如圖1所示為PT法與FT法原理示意圖。

圖1 PT法與FT法原理示意圖Fig.1 Principle diagram for profile transform (PT) and Fourier transform (FT) methods

如圖1a所示，PT方法是最簡單最直接的方法，其原理是按照式（4）對變量進行縮放，從而達到上下游通量守恒的目的，但是其對角度比有嚴格的限制，要求接近1。

式中1pφ與2pφ分別為交界面上游及下游通量，l為角度比。

FT方法對角度比沒有嚴格限制，使用一種采樣方法對求解歷史數據進行重構，是將前幾次計算生成的數據作為一個采樣點用來進行下一次計算，然后對采樣點進行傅里葉級數展開，獲得所關心的物理量在某一位置上的函數關系式。FT方法采用相變邊界條件，其圓周方向在不同時刻彼此成周期。按照式（5）轉子與定子在中間面上進行流動變量的采樣，從而計算出傅里葉系數[22]，如圖1b所示。進一步，FT方法利用基于周向角度差的相角滯后及轉輪速度，對周期面上的流動變量用傅里葉系數進行重新構建，如式（6）和（7），因此不用存儲完整的計算信號，提高了計算速度。

式中sφ，1pφ，2pφ分別為采樣面，2個周期面上的流動變量，t與Δt分別為時間和時間步長，Am為傅里葉系數，m為采樣數，ω為角速度，PR和PS分別為一個流道內轉子與定子的周向角度，R為半徑。

1.3 計算模型及網格生成

Franciss99模型主要部件包括一個內嵌14個固定導葉的蝸殼、28個活動導葉，帶有15個長葉片及15個短葉片的轉輪及一個彎肘形尾水管。模型機與原型轉輪出口直徑分別為0.349和1.778 m，兩者間的縮放比為1:5.1。最優工況下，模型機組活動導葉開度α為9.84°，試驗水頭H為11.91 m，轉速n為35.12 rad/s，按照IEC 60193，流量系數QED及轉速系數nED分別為0.15和0.18，對應的原型機組水頭及出力分別為377 m和110 MW[25-27]。

計算域離散采用精度較高的塊結構化六面體網格對計算域進行離散，對活動導葉及轉輪葉片，整體采用“H”形塊，然后在翼型壁面使用一個“O”形塊包裹葉片以便提供更好的網格質量，并且對所有的壁面處進行局部加密處理以期望獲得更好的壁面求解。采用 5種不同密度的網格進行網格無關性驗證，網格總數分別為306.7萬、651.2萬、960.4萬、1 227.9萬和1 523.4萬，圖2所示為水輪機效率及扭矩隨網格數的變化規律。由無關性測試結果可知，當網格數大于 960萬時，水輪機水力效率與轉輪扭矩隨網格數的增加不再變化，因此考慮到計算資源與計算精度的折衷，數值計算網格數最終選取為960.4萬。如圖 3所示為數值計算采用的轉輪及活動導葉網格示意圖。

數值計算域各過流部件網格信息如下，蝸殼（包括活動導葉）、活動導葉、轉輪及尾水管網格的最小角度分別為21.4°，24.7°，28.8°和35.4°，對應的網格縱橫比為77、96、242和83，最小網格質量分別為0.32、0.41、0.45和0.62。第一層網格距壁面的距離分別為 0.3，0.1，0.15與0.3 mm，在最優工況計算得到的壁面最大y+值（ y+為第1層網格距離壁面的無量綱距離）分別約為 65、88、67和87，平均y+分別為28、32、16和22。

圖2 網格無關性測試Fig.2 Mesh independence test

圖3 過流部件網格示意圖Fig.3 Grid views of different flow passage components

2 計算工況及設置

在最優工況（best efficiency point—BEP）下分別采用全通道（Full），單流道的 PT方法以及雙通道的 FT方法進行水力性能預測，3種不同計算策略模擬域如圖4所示。由于Francis99模型轉輪為長短葉片形式，轉輪一個流道必須包含一個長葉片與一個短葉片，故將 2個活動導葉與一對長短葉片的組合作為PT方法計算域，同理，FT方法計算域為4個活動導葉與2對長短葉片的組合。

圖4 三種模擬策略計算域示意圖Fig.4 Views of fluid domain for three simulation strategies

首先進行全通道穩態計算，獲得一個穩定的初場。數值計算進口給定質量流量，出口指定靜壓，壁面均采用光滑、無滑移條件。瞬態計算動靜交界面為“Transient Rotor Stator”，對流采用高階求解格式，瞬態模型則采用二階向后歐拉模式，收斂標準設為最大殘差小于0.001。為消除進出口條件不同帶來的誤差，將全通道提取出的PC（profile boundary conditions）條件作為PT法及FT法的進出口條件，進口為速度矢量，出口為靜壓，數值計算中僅僅與計算域重合的PC邊界條件是有效的。計算過程中，通過設置一系列監測點來動態監測積分變量及壓力、速度等值的變化，一方面可以對動量方程提供反饋，另一方面起到診斷計算是否正確的目的。圖 5所示為試驗裝置壓力傳感器安裝位置[23]，數值計算在相同位置布置壓力測點以進行驗證和對比。

圖5 試驗傳感器與數值監測點位置示意圖Fig.5 Locations of experimental pressure probe and monitoring points of numerical simulation

如圖5所示，試驗測量與數值計算分別進行了6個不同位置的壓力測量，其中VL01位于轉輪與活動導葉之間的無葉區，DT11和DT21安裝在尾水管內，安裝位置Z=-0.305 8 m，剩余3個壓力傳感器P42、P71和S51安裝在轉輪葉片上，隨轉輪做旋轉運動。其中 P42和 P71在葉片壓力面，P42位于葉片中心位置，P71在近葉片出水邊且靠近輪轂一側，S51位于吸力面約葉展2/3、靠近輪轂處。

3 結果與分析

為綜合對比不同計算方案在數值計算精度上的預測能力，本文主要從水動力學特性、轉輪內部壓力負荷分布、時變壓力信息及其頻譜特性等方面進行分析。在配置相同工作站（戴爾T7910，至強雙處理器E5-2680 V3，24核心，內存80G）且計算步數、收斂準則相同的條件下，全通道計算、PT方法以及FT方法計算時間比約為1:0.375:0.23，表明PT法與FT法對計算資源的需求大大減小。

針對Francis99模型，國外進行了一系列的數值及試驗研究。Trivedi等[25]對Francis99模型進行了詳細的模型試驗及數值研究，Lucien等[28]進行了Francis99水輪機部分流道的穩態及全流道的瞬態計算，Buron等[29]采用非線性諧波法模擬了Francis99模型的RSI效應。上述研究獲得比較一致的結果為：數值模擬的水輪機水力效率均高于試驗值；不同位置測點上的平均壓力在轉輪與活動導葉之間的無葉區及葉片表面均高于試驗值，而尾水管內則低于試驗值；測點壓力脈動頻譜特性方面，葉片正面中心位置測點幅值最高、無葉區次之。本文主要通過與試驗數據以及上述研究的計算結果進行對比分析，以評估PT方法以及FT方法在水輪機水力性能預測的可靠性。

3.1 水力效率

水力效率是水輪機整機性能的表征參數，表 1統計了 3種不同計算方法獲得的水力效率與試驗值，數值效率計算采用進出口總壓差法。由效率統計表可知，3種數值方法計算的水力效率均與試驗值吻合得比較好，數值解均高于試驗值，這是由于數值計算對模型進行了一定的簡化，如未考慮轉輪上冠和下環的迷宮環密封，未計容積損失、圓盤摩擦損失等。試驗效率值為93.4%，全通道方法為94.6%，而PT和FT方法均為95.8%，3種數值計算效率偏差分別為 1.28%，2.57%，2.57%。PT與 FT方法的結果相對于試驗測試偏差高于全通道計算，分析認為主要由兩種原因所致。一是PT與FT方法僅僅對活動導葉與轉輪的部分流道進行模擬，忽略了蝸殼、活動導葉及固定導葉區域的損失；二是PT與FT方法進出口條件均為由全通道計算結果提取出來的變量邊界條件（profile boundary），計算條件的設置可以認為與全通道相同，而這 2種方法均為非全通道計算，活動導葉及轉輪區域固壁面的水力損失均小于全通道方法，故其計算效率偏差略大。

表1 試驗與3種不同數值方法的水力效率統計結果Table 1 Statistics results of test hydraulic efficiency and numerical hydraulic efficiency with three different simulation methods

3.2 壓力及速度分布

為進一步驗證數值模擬的可靠性，本文對比了由活動導葉與轉輪之間無葉區至尾水管內各監測點上的平均壓力的計算值和試驗值，如圖6所示。BEP工況下，不同位置測點平均壓力數值解和實測值的差異較小，在無葉區和轉輪葉片上，平均壓力的數值解均高于實測值。對全通道方法而言，其預測的尾水管中平均壓力的計算值小于實測值。全通道計算方法、FT及PT方法捕捉的壓力均值僅僅在無葉區有較小區別，而轉輪葉片上的數值差異幾乎可以完全忽略，表明3種計算方法在預測計算域內平均壓力上效果相當，且與實測值以及Trivedi等[22,25,28-29]的結果比較一致。

圖6所示分別為采用PT方法、FT方法以及全通道方法計算的 0.5倍葉片高度活動導葉與轉輪流道內局部靜壓分布（圖 7左側）、轉輪區域內速度云圖及速度矢量分布（圖 7右側）。由靜壓分布可以清楚地看到，由活動導葉進口至轉輪出口，壓力值均勻減小，未出現突變，表明做功良好。PT法與全通道計算在靜止域與旋轉域內的壓力顯示幾乎一致，而FT法與其余兩種方法在活動導葉進口的靜壓分布上有區別，主要表現在兩方面，一是活動導葉頭部對壓力分布的影響，二是總體壓力值偏低。3種計算方法速度矢量分布比較一致，轉輪流道內流線順暢，未出現漩渦及回流等現象，表明能量轉率效率較高。相對于PT方法及全通道計算，FT方法提取的速度矢量在長葉片進口側壓力面上輕微沖擊壁面，分析認為FT方法出現不同于其余2種方法是由于其所采用的相變條件向上游的反饋所致。

圖6 試驗與3種計算方法平均壓力比較Fig.6 Average pressure comparison of experimental test and numerical simulation with three simulation methods

圖7 三種不同模擬方法流道壓力及速度分布Fig.7 Pressure and velocity distribution at mid-span plane by three simulation methods

為定量分析流道內壓力分布情況，圖 8顯示了葉片上壓力負荷沿弦長的分布。

圖8 轉輪葉片負荷分布Fig.8 Blade load distribution for runner

由圖8可知，單通道的PT方法以及雙通道的FT方法模擬結果相當，且均獲得了與全通道結果比較一致的負荷分布，在相對弦長小于0.8時，PT與FT方法在葉片壓力面上的負荷大于全通道方法，其余位置則吻合的較好。結果分析表明，PT與FT方法均可較準確捕捉葉片表面負荷分布。

3.3 葉片合力及扭矩分布

水力效率及不同位置監測點壓力均值，均為采用時均法獲得，其對外特性及時均流動特性的捕捉比較準確，但水輪機內部流動特征具有強烈的非穩態性，因此應進行時變量的深入分析。圖 9所示為一組長短葉片所受合力及扭矩隨時間變化規律。由時變葉片受力及扭矩可知，FT與PT法預測結果均小于全通道計算，然而FT方法與全通道比較接近，在幅值的一側與全通道計算有差別，而另一側則吻合地較好，且波峰及波谷出現時間間隔相同，僅僅存在一定的相位差，PT方法僅僅捕捉到較小波動且其最大最小值時間間隔與全通道有差別。可以斷定，利用采樣方法及相變邊界條件的FT法能較好地還原時變信號，而PT方法誤差較大。

圖9 葉片合力及扭矩隨時間變化規律Fig.9 Temporal variation of force and torque on blade

3.4 壓力脈動頻譜特性

本文進一步對數值計算及模型試驗各監測點瞬時壓力數據進行快速傅里葉變換（FFT），獲得各個測點的頻譜特性。

模型試驗獲得的時變壓力信息中混摻著電網頻率及其諧波，往往在頻率中占據主導作用，對原始信號有較強的干擾作用。因此，在進行試驗壓力脈動及其頻譜分析時，必須進行50 Hz及其諧波的濾波處理以消除電網頻率的影響。本文采用Matlab對試驗原始信號進行濾波處理。圖10為試驗條件下、BEP工況濾波前后傳感器P42頻譜圖。由圖10可知，傳感器采集的原始信號的幅值在100，200及300 Hz明顯高于其他頻率，采用濾波器對信號進行去噪處理，有效去除了電網頻率的干擾。

圖10 測點P42濾波前后頻譜分析對比Fig.10 Spectral analysis of P42 before and after filtering

圖11 所示為BEP工況下全通道計算、PT方法、FT方法以及試驗結果計算的不同監測位置的壓力脈動頻譜圖，為時變壓力數據經過快速傅里葉變換（FTT）而來。圖11中，橫軸坐標為經過歸一化處理的無量綱量，為實際頻率與對應工況下轉頻的比值，即轉頻的倍數，縱坐標為壓力脈動幅值。

圖11 壓力脈動頻譜結果顯示，FT方法幾乎再現了全通道方法的所有特征，2種不同的數值方案均準確捕捉到了靜止域與旋轉域測點的主頻，即VL01主頻為30fn（fn為轉頻），P42、S51、P71主頻為28fn，而PT方法預測的所有主頻均為30fn，表明FT方法較準確地還原了原始信號的求解歷史，而PT方法由于采用在交界面上對變量進行比例縮放的策略，從而引起頻率變化。此外，FT方法與全通道算的明顯區別是各個測點上均出現了諧波分量，這是由于信號重構采用了傅里葉分解的緣故。幅值方面，FT法與全通道法幅值大小變化趨勢一致，即P42壓力幅值最大，VL01次之，然后為S51以及P71。PT法同樣為P42最大，而S51次之，然后依次為VL01與P71。相對于試驗值，FT方法與全通道在VL01上幅值高于試驗值，其余測點則與試驗值比較一致，而PT方法VL01上與試驗值比較吻合，而其余測點均小于試驗值。通過與Nicolle等[22]的壓力脈動頻譜特性結果對比發現，本文所采用的FT方法捕捉到的壓力脈動特性在主頻及幅值上與試驗值及全通道計算均有較好的一致性，應該與計算域的選取以及網格離散有關。綜合對比分析表明，FT方法可以以較少的計算資源捕捉到與全通道相當的RSI效應。

4 結 論

本文以Francis99高水頭混流式模型水輪機為研究對象，采用TBR模型中單通道PT（profile transform）法和雙通道FT（fourier transform）法對水輪機水力性能以及負荷進行了預測，并與試驗值及全通道結果進行對比，獲得了以下結論：

1）不同的計算方法在水輪機外特性效率的計算上均有較高的精度，水力效率均高于試驗結果，而部分通道計算由于忽略了蝸殼、固定導葉及尾水管的水力損失，故相對于試驗值偏差略大。

2）計算域內不同位置平均壓力僅僅在導葉與轉輪之間的無葉區有差別，其余位置幾乎相同。計算域內由活動導葉進口至轉輪出口壓力梯度均勻變化，在葉片負荷分布上PT法、FT法以及全通道法差別不大，PT法與FT法與全通道計算的區別主要出現在相對弦長小于 0.8葉片壓力面上。

3）FT法預測的葉片受力和扭矩脈動量及其頻率與全通道更接近，而PT方法偏差較大。

4）FT方法在頻譜特性的預測上與全通道結果相當一致，對不同測點主頻及其幅值均預測的較準確，而PT方法由于交界面上進行了變量的比例縮放，產生了頻率變化，且預測幅值小與試驗值。

綜合評估分析可知，FT方法在水輪機外特性、內部流動特征以及時變信號的捕捉上與全通道計算相當，其計算時間僅為全通道3/8，該方法在水輪機的瞬態流動數值計算中具有一定的潛力和優勢。