劉穎超，陳 榴，戴 韌，楊其國

(上海理工大學 能源與動力工程學院，上海 200093)

大功率凝汽式汽輪機排汽壓力低，低壓末級在高蒸汽比容和速度下工作。百萬千瓦級汽輪機的末級全速鋼制葉片長度已達1 300 mm，葉片頂部圓周速度約600 m/s。如果軸向絕對出口馬赫數(shù)接近1，則動葉頂部截面葉型相對出口馬赫數(shù)可達1.7左右[1]。末級動葉的氣動設計是個難題，末級的低輪轂比和超長葉片導致末級從葉根到葉頂?shù)牧鲃铀俣?、角度和葉型變化劇烈，尤其是葉片頂部截面葉型在滿足跨音速流動氣動設計時，還必須滿足嚴苛的許可應力要求。

跨音速葉柵流動中不僅有激波和膨脹波，還有波系之間的相互作用，激波與邊界層之間相互干涉，流動復雜而難以準確計算葉柵出流參數(shù)。早期計算方法是二元特征線法，其通過假設兩相鄰葉片喉部處音速線是一條直線來求解渦輪葉柵跨音速流動。Deich[2]通過分析得出特征線法在葉柵出口馬赫數(shù)低于1.4、低節(jié)弦比和薄直背葉型的葉柵流動中相當準確。而對于高節(jié)弦比、高反動度的葉柵頂部薄葉型，F(xiàn)orster[3]的實驗結果說明葉柵喉部的直線音速線假設已不成立，計算方法失效。

葉柵跨音速流動的實驗測量也十分困難。平面葉柵實驗中，葉片尾緣的激波在葉柵兩側自由流邊界反射，形成“反射激波”，影響葉柵出口流動參數(shù),在大安裝角和葉片間距大的葉柵中很容易發(fā)生這種現(xiàn)象。另一方面，在超音速氣流中，氣動探針難以精確測量當?shù)乜倝海捎眉y影法可以獲得激波形態(tài)，但是不能量化評定葉柵效率[4]。

基于時間推進法的CFD求解方法為分析跨音速葉柵流動提供了有力工具。分析激波相關流動，需要各種高精度格式，如總變差減小的差分格式(TVD)[5]、Roe迎風格式[6]、非自由參數(shù)耗散有限差分格式(NND)[7]等，目前對包含中等激波強度的跨音速葉柵流動計算基本沒有困難，大渦模擬(LES)[8]計算也是可行的。具有2階迎風格式的計算方法基本可以滿足葉柵跨音速流動分析的需要，結合適當?shù)膬?yōu)化算法，可以完成跨音速葉柵的優(yōu)化設計[9-11]。

系列實驗[12-13]證明，跨音速葉柵流動損失與葉型幾何形狀密切相關。葉型背弧的后半段曲率是重要的影響因素，大曲率出氣邊的葉型損失在超音速時是亞音速時的4倍，其中激波誘導邊界層分離是損失增加的主要原因，這時采用直線背弧是更好的選擇。前緣氣流加速主要由葉片形狀決定，并在某種程度上與壓比無關，而尾緣對流場的影響主要取決于壓比，尾緣形狀的影響較小。各種經典葉頂截面葉柵跨音速流動[14]的流動性能大致相當，主要損失是葉柵尾緣激波和吸力面上的反射波損失。

筆者應用CFD軟件分析某汽輪機末級動葉頂部截面葉型流動損失機理，結合葉柵氣動設計方法，通過參數(shù)優(yōu)化，調控葉柵流動中的激波強度和反射位置，降低葉柵內與激波有關的流動損失，提高葉柵流動效率。

1 葉柵流動分析的數(shù)值方法

1.1 驗證模型與邊界條件

流動分析選用CFX軟件，采用有限體積法求解雷諾平均N-S方程，守恒方程中的對流項和湍流輸運方程均選擇2階差分格式。選取類似末級動葉頂部截面的雙圓弧葉柵實驗[15]作為參照，驗證數(shù)值方法對跨音速葉柵流動特性計算的可靠性。該葉型相對柵距為1.5，安裝角為30°，最大相對厚度為0.1。

參考實驗工況，計算域進口和出口分別給定總壓p*和靜壓p2，壓比p2/p*=0.23，進口氣流角為35°，進口總溫T*=280 K。周向采用周期性邊界條件，葉片表面為絕熱邊界條件。

1.2 網格驗證

計算網格采用結構網格，計算域進口距葉型前緣距離為C(其中C為弦長)，出口距葉型尾緣距離為2C，邊界層伸長比小于1.1。如表1所示,采用3套不同質量網格進行流動計算，其中Cp為葉型表面壓力系數(shù),湍流模型選用SST模型，湍流度為5%，收斂后求解殘差小于10-6。

表1 網格無關性比較Tab.1 Comparison of the grid independence

由表1可以看出，對出口馬赫數(shù)、激波前Cp峰值及相應位置的預測，3套網格之間略有差別，中等網格與細網格的預測效果基本一樣，表明中等網格的求解精度滿足性能分析要求。

1.3 湍流模型驗證

分別采用標準k-ε、k-ω和SST湍流模型進行定常計算，得到的葉型表面壓力分布Cp如圖1所示，其中p表示當?shù)貕毫?X表示葉型的橫坐標。與實驗對比發(fā)現(xiàn)，3種湍流模型預測葉型壓力面的壓力分布與實驗值均較吻合，差異不大，但對吸力面的預測效果存在偏差。標準k-ε模型對近壁區(qū)求解精度不高，在0.8C以后壓力一直降低，無法準確預測激波。k-ω和SST模型預測所得的吸力面激波前壓力高于實驗值，激波位置與實驗值吻合。在吸力面前緣附近，SST模型與實驗值更加接近，因此筆者選擇SST模型。

圖1 采用不同湍流模型計算所得葉型表面壓力分布Fig.1 Blade surface pressure distribution of different turbulent flow models

2 跨音速葉柵的流動分析

2.1 計算模型及邊界條件

以某汽輪機末級動葉跨音速頂部截面葉型為原葉型，其葉柵結構如圖2所示，葉片弦長為150 mm，相對柵距為0.926，最大厚度為6.09 mm，其他相關幾何參數(shù)如表2所示。

圖2 葉柵結構示意圖Fig.2 Schematic diagram of cascade structure

表2 葉型幾何參數(shù)Tab.2 Geometric parameters of the profile

計算域進口總壓p*為17 021 Pa，出口靜壓p2為4 845.4 Pa，進口總溫T*為308.1 K。不考慮末級流道在子午面的流線翹曲和子午加速率。湍流模型采用SST模型，湍流度為5%。計算網格仍采用上述結構網格。

2.2 葉柵流動的激波結構

圖3為原葉型表面和沿葉柵流道中線MN(見圖4中)的等熵馬赫數(shù)分布。圖4給出了葉柵流道內馬赫數(shù)分布云圖，由圖4能清晰看到膨脹波與激波結構。在吸力面前緣，氣流迅速達到超音速后再減速(AB1段)，產生了圖4中AB位置的激波。然后,超音速氣流在葉柵流道內流經相鄰尾緣壓力面附近，遇膨脹波K1加速，并產生反射膨脹波K2，加速至吸力面53%C處達到第2次速度峰值，即C點，等熵馬赫數(shù)接近1.8。之后遇到尾緣壓力面DC處斜激波減速(CD1段)，并在CE位置產生反射波(C2E1段減速)。由于出口馬赫數(shù)高，在尾緣吸力面FG處也形成一道激波。

圖3 葉柵內等熵馬赫數(shù)分布Fig.3 Isentropic Mach number distribution in the cascade channel

圖4 葉柵流道內馬赫數(shù)分布Fig.4 Mach number distribution in the cascade channel

尾緣處吸力面和壓力面激波通常稱為“燕尾波”。各種波的存在、傳播和相互作用，尾緣吸力面激波、反射波與尾跡(K3)干涉，以及激波誘導邊界層分離，使得葉柵流動更加復雜。

2.3 葉柵氣動損失分析

跨音速葉柵氣動損失來源包括激波損失、邊界層損失和尾緣損失(取決于尾緣厚度和基壓)。衡量和評價不同損失對提高葉柵氣動效率十分重要。

從熱力學分析，流場中不可逆熵增是損失來源，可以用熵增來衡量葉柵內各種損失?？缫羲偃~柵中復雜激波造成的氣動損失是熵增來源之一。利用激波前后的質量、動量和能量守恒，能計算激波引起的熵增Δsshock，表達式[16]如下：

(1)

式中：cv為比定容熱容；k為絕熱指數(shù)；Manormal為激波上游法向馬赫數(shù)，是由圖4中激波位置確定的平均變量。

以進口溫度T0與進口速度u0對激波熵增無量綱化,得到無量綱激波熵增Δsshock-normalized：

(2)

表3統(tǒng)計了葉柵內不同位置的激波損失。由表3可知,尾緣處激波引起的損失最大，占比58.5%，其次是反射波，這是因為渦輪的高壓比導致尾緣后基壓低，當流體在尾緣周圍膨脹到這種低壓，在吸力面和壓力面氣流相遇的地方被強激波重新壓縮[17]，使熵值增加。而反射波位于喉口后，是高馬赫數(shù)集中區(qū)域，損失增大。

表3 葉柵不同位置激波損失Tab.3 Shock loss at different positions of the cascade

圖5給出了葉型邊界層厚度分布，在經過相鄰尾緣壓力面激波后，吸力面位移厚度突增形成波峰，表明激波干涉邊界層，造成局部邊界層分離，邊界層增厚，氣動損失增加。葉片邊界層內黏性效應導致能量耗散，是不可逆過程，邊界層損失用熵增Δsbl計算，定義[16]如下：

圖5 葉型邊界層厚度分布Fig.5 Boundary layer thickness distribution of the blade

(3)

由進口溫度與進口速度對邊界層引起的無量綱熵增為：

(4)

通過上式計算，邊界層無量綱熵增損失為0.157 8，其中還有壓力側激波與吸力面邊界層干涉而附加的損失，而無量綱激波熵增損失為0.255 4，能量損失更大，故削弱激波強度、減小激波相關損失對提升葉柵氣動性能是首要的。通過調整優(yōu)化葉型幾何形狀，改變曲率分布，控制尾緣激波及其反射波，以降低葉柵內激波損失。

3 葉型氣動優(yōu)化方法

3.1 葉型的參數(shù)化

采用中弧線疊加厚度分布生成葉型[18]。中弧線采用3階Bezier曲線表示，如圖6所示。n階Bezier曲線表達式如下：

圖6 采用3階Bezier曲線表示的葉型中弧線Fig.6 Camber line of the blade represented by the third-order Bezier curve

(5)

式中：p(t)表示構成Bezier曲線上的點;Bi,n(t)為Bernstein函數(shù)；Pi為直角坐標系下控制點的位置矢量。

葉型中弧線首、末兩點由葉型前緣點、尾緣點和安裝角確定。前后緣的斜率設定為進、出氣流角的約束，控制點P2、P3只能分別沿直線P0P2和P3P4移動，故中弧線形狀由2個參數(shù)XP2、XP3控制。

葉型厚度分布如圖7所示，以最大厚度控制點T3為界，左側用2階Bezier曲線表示，右側用3階Bezier曲線表示。通過葉型前緣厚度和尾緣厚度確定首、末2點T1、T6?？刂泣cT2、T3、T4的厚度均為最大厚度，葉型厚度分布由控制點T3、T4和T5確定，故葉型厚度分布由5個參數(shù)XT3、YT3、XT4、XT5、YT5控制。

圖7 Bezier曲線表示葉型厚度分布Fig.7 Thickness distribution of the blade represented by Bezier curve

3.2 基于Kriging代理模型的優(yōu)化

采用基于克里金法(Kriging)的代理模型，對動葉頂部截面葉型進行氣動優(yōu)化。Kriging代理模型[19]具有精確的全局近似設計空間，能準確地預測線性和非線性函數(shù)，其表達式如下：

(6)

試驗設計(DOE)是構建Kriging代理模型的第一步，建立設計空間生成樣本。采用最優(yōu)拉丁超立方設計在設計變量范圍(見表4)內生成30個樣本，作為Kriging代理模型的訓練樣本集。

表4 設計變量取值范圍Tab.4 Value range of design variables

優(yōu)化的目標函數(shù)選擇能評價葉柵氣動損失的出口總壓損失系數(shù)Cp,t，其具體定義為：

(7)

式中：p1為進口靜壓；p02為出口總壓。

葉型前緣和尾緣直徑分別為2.5 mm和1.4 mm。設定的約束條件是葉型截面面積的變化相對于原葉型不超過±10%。在Isight軟件上，基于DOE生成樣本，選擇Kriging代理模型建立目標函數(shù)的近似模型，并在設計變量范圍內，采用多島遺傳算法尋優(yōu)，得到最低Cp,t所對應的設計變量數(shù)值,然后用CFX校核最優(yōu)葉型的性能。

3.3 優(yōu)化結果分析

經過優(yōu)化計算，最終獲得優(yōu)化后葉型，圖8給出了優(yōu)化前后葉型的幾何構型對比,圖中采用葉型弦長對葉型軸向和周向坐標進行無量綱化。優(yōu)化后葉型截面面積比原葉型小9.6%，最大厚度為5.78 mm，位置基本不變。在葉型80%C～90%C處，優(yōu)化葉型壓力面曲率明顯減小，型線為直線，有利于減小壓力梯度，削弱尾緣壓力面激波強度。通過Kriging模型預測得到Cp,t為0.215 6，CFX實際計算所得Cp,t為0.216 7，代理模型不確定度為0.51%。

圖8 優(yōu)化前后葉型對比Fig.8 Comparison between optimized and baseline profiles

葉片表面連續(xù)曲率影響葉片性能[20-21]。圖9比較了優(yōu)化前后葉型吸力面曲率的分布，優(yōu)化后葉型吸力面曲率整體更加連續(xù)平滑，在20%C～75%C處曲率均低于原葉型，75%C后優(yōu)化后葉型曲率略大于原葉型，從而決定葉型表面等熵馬赫數(shù)的分布形狀(見圖10)。由圖10可知，優(yōu)化后葉型前緣吸力面的速度減小，吸力面上的第2次速度峰值降低并向后移動，延遲了壓力峰值點，尾緣壓力面激波入射位置后移。從圖11優(yōu)化前后沿葉柵流道中線MN的等熵馬赫數(shù)分布也可以看到，反射波位置后移，激波強度削弱。同時，優(yōu)化后葉型接近尾緣的壓力面維持了較高壓力，起到類似后加載的作用，有利于降低葉型氣動損失。圖12給出了優(yōu)化后葉柵流道內馬赫數(shù)分布。由圖12可以看出,優(yōu)化后葉柵喉口后高馬赫數(shù)集中區(qū)域縮小，吸力面反射波區(qū)域縮小，強度減弱。

圖9 優(yōu)化前后葉型吸力面曲率分布對比Fig.9 Comparison of curvature distribution of blade suction surface before and after optimization

圖10 優(yōu)化前后葉型表面等熵馬赫數(shù)對比Fig.10 Comparison of isentropic Mach number distribution of the blade surface before and after optimization

圖11 優(yōu)化前后沿葉柵流道中線MN的等熵馬赫數(shù)分布對比Fig.11 Comparison of isentropic Mach number distribution along the centerline MN of the cascade channel before and after optimization

圖12 優(yōu)化后葉柵流道內馬赫數(shù)分布Fig.12 Contour distribution of Mach number distribution in cascade channel of optimized blade

圖13比較了優(yōu)化前后葉型不同位置的激波損失，優(yōu)化后葉型的各處與激波相關的損失均下降，前緣AB處激波基本消失，反射波損失減少26%，尾緣處激波損失減少34%，總的激波損失相比原葉型下降約32%。同時，激波強度減弱能減少對邊界層的干擾，邊界層損失降低。優(yōu)化后葉型的邊界層損失Δsbl-normalized為0.130 8，相比原葉型降低了17%。

圖13 優(yōu)化前后葉型不同位置激波損失對比Fig.13 Comparison of shock loss at different positions of the blade profile before and after optimization

圖14和圖15給出了優(yōu)化前后葉柵出口1.2倍軸向弦長Cax處尾跡速度和尾跡氣流角沿相對柵距的分布,圖中采用葉柵出口截面平均速度對尾跡速度進行無量綱化。由圖14和圖15可以看出，優(yōu)化葉型低速幅值降低,尾跡寬度有所減小，尾跡氣流角幅值有所降低，有利于減小尾跡摻混損失。優(yōu)化后，吸力面?zhèn)任槽E氣流角的變化更加緩慢，壓力面?zhèn)任槽E速度降低，尾跡氣流角的變化相對平緩，表明尾緣燕尾波減弱，尾流干涉減弱。

圖14 優(yōu)化前后葉型尾跡速度分布Fig.14 Circumferential velocity distribution of the blade before and after optimization

圖15 優(yōu)化前后葉型尾跡氣流角分布Fig.15 Flow angle circumferential distribution of the blade before and after optimization

圖16比較了優(yōu)化前后葉型沿流向方向的總壓損失系數(shù)。由圖16可以看出，優(yōu)化后葉型的總壓損失降低，特別是在中間吸力面后，相鄰尾緣壓力面激波干擾減弱，降幅增大，表明優(yōu)化后葉型的氣動性能提升。優(yōu)化后葉型出口總壓損失比原葉型降低約24%，質量流量增加4.8%。

圖16 優(yōu)化前后葉型流向方向的總壓損失系數(shù)對比Fig.16 Comparison of total pressure loss in the blade flow direction before and after optimization

4 結 論

(1)末級動葉跨音速頂部截面的葉柵流動損失主要來自激波，激波帶來的熵增損失比邊界層產生的熵增損失多61%，則跨音速葉柵主要通過削弱激波來減少氣動損失。

(2)優(yōu)化葉型曲率分布是降低激波損失的關鍵。葉型尾緣的厚度決定了尾緣燕尾波的強度，吸力面的曲率決定了反射波的位置和強度。優(yōu)化后葉型激波產生的熵增損失相比原葉型下降約32%，葉型出口總壓損失系數(shù)比原葉型減少約24%，優(yōu)化后葉型的氣動性能明顯提升。