基于伴隨方法的單級低速壓氣機氣動設計優(yōu)化

羅佳奇，楊婧

1. 浙江大學 航空航天學院，杭州 310027 2. 德州農(nóng)工大學 機械工程系, College Station TX 77843

隨著高性能計算的快速發(fā)展，基于計算流體力學(CFD)的優(yōu)化設計技術(shù)在葉輪機械葉片氣動設計中已經(jīng)得到了應用?；贑FD的氣動設計優(yōu)化(ADO)不僅能夠縮短研制周期、降低成本，相對于依靠設計者經(jīng)驗的設計方法還能提高設計結(jié)果精度。對于壓氣機而言，通過優(yōu)化葉片氣動外形增加壓比、降低流動損失對于提高壓氣機氣動性能具有重要意義。

目前得到應用的優(yōu)化方法分為兩類：梯度方法和非梯度方法。進化算法[1-3]和代理模型[4]是典型的非梯度方法。Benini[2]采用進化算法對跨聲速壓氣機轉(zhuǎn)子葉片NASA Rotor 37的中弧線和厚度分布進行了優(yōu)化設計，Lian和Liou[3]采用進化算法和代理模型對跨聲速壓氣機轉(zhuǎn)子葉片NASA Rotor 67進行了多目標優(yōu)化設計，Samad等[4]采用多層次代理模型對轉(zhuǎn)子葉片的葉片掠、葉片傾斜及葉片彎進行了優(yōu)化設計。上述非梯度方法雖然在葉輪機械ADO中已經(jīng)得到應用，但是對于多設計參數(shù)的復雜ADO問題，尤其是需要通過求解Navier-Stokes方程計算流場時，將需要大量的計算資源，制約了這類方法的工程應用。

基于靈敏度分析的優(yōu)化方法在葉輪機械ADO中得到了應用，直接差分法(FDM)和伴隨方法是兩種常用的靈敏度計算方法。采用FDM計算靈敏度時，所需的流場計算次數(shù)和設計參數(shù)線性相關(guān)。在20世紀80年代，Jameson[5]成功地將伴隨方法應用于飛行器氣動外形優(yōu)化設計。采用伴隨方法，每個氣動函數(shù)的靈敏度計算只需要大約兩次流場計算，計算時間和設計參數(shù)個數(shù)基本無關(guān)。之后，基于伴隨方法的翼型、機翼以及翼身組合體ADO研究得到了迅速發(fā)展[5-9]。由于靈敏度計算的高效性和精確性，伴隨方法在葉輪機械ADO中的應用也越來越廣泛[10-19]。Luo等[17-18]開展了基于伴隨方法的跨聲速壓氣機葉片的多工況、多目標ADO研究。近年來，摻混面模型方法也應用于伴隨方程的數(shù)值求解，推動了基于伴隨方法的葉輪機械多排ADO研究的發(fā)展[13,14,16,20]。

本文將發(fā)展一種基于摻混面模型的葉輪機械多排伴隨方法，并應用于某型低速、低壓比壓氣機ADO。首先采用自主開發(fā)的初步設計方法設計4.5級壓氣機葉片的氣動外形，通過數(shù)值模擬發(fā)現(xiàn)：最后級靜子葉片的端壁和吸力面角區(qū)內(nèi)存在較大分離，嚴重影響了壓氣機效率。之后介紹適用于多排ADO的伴隨方法，并在近失速工況下對最后級靜子葉片的中弧線及安裝角進行優(yōu)化，優(yōu)化目標是降低出口熵增，同時約束出口流量。最后開展多工況ADO，在兩個不同轉(zhuǎn)速下進行最后級靜子葉片的外形優(yōu)化設計，分析外形優(yōu)化設計對設計工況、非設計工況氣動性能的影響。

1 壓氣機初步設計

該低速4.5級壓氣機由第1排進口導葉和后4級葉排構(gòu)成，給定設計壓比和等熵效率。在初步設計之前需要確定一些前提條件：單級多變效率和壓氣機相等，給定單級溫比，壓氣機進口到導葉出口的面積收縮率等。采用基于經(jīng)驗修正[21]的4.5級壓氣機初步設計的基本流程為：

第1步 給定效率初值，確定壓氣機進、出口流道幾何參數(shù)、速度三角形和氣動參數(shù)。

第2步 根據(jù)給定的溫升分布、靜子出口氣流角和流道參數(shù)等，進行逐級設計，迭代計算各級轉(zhuǎn)子、靜子出口截面的相關(guān)參數(shù)。

第3步 計算整機壓比和效率，判斷是否滿足設計條件，如果不滿足則用計算效率替代初始效率，重復設計。

第4步 采用簡單徑向平衡方程確定輪轂和機匣的相關(guān)參數(shù)。

第5步 確定壓氣機軸向尺寸，進行二維葉片造型和徑向堆積。

圖1為該4.5級壓氣機輪轂、葉中和機匣處的葉型，圖中：x表示 軸向，Rθ表示周向。值得注意的是：初步設計時，導葉和靜子葉片均采用直葉片設計。流場分析發(fā)現(xiàn)：直葉片設計使得在前三級的影響下，最后級靜子葉片的端壁和吸力面角區(qū)流動分離嚴重；因而，該初步設計的4.5級壓氣機存在較大優(yōu)化空間，可以通過葉片外形優(yōu)化設計抑制流動分離，提高壓氣機氣動性能。

圖1 3個不同徑向位置4.5級壓氣機葉片葉型Fig.1 Blade profiles of the 4.5-stage axial compressor at three different spans

2 氣動設計優(yōu)化

ADO的流場計算采用流場求解器Turbo90[22-24]，給定目標函數(shù)及優(yōu)化約束的數(shù)學建模以確定伴隨方程的邊界條件并實現(xiàn)數(shù)值求解。確定伴隨靈敏度后，采用簡單、易實現(xiàn)的最速下降法來進行優(yōu)化設計，氣動外形擾動步長根據(jù)經(jīng)驗給定。

2.1 流場數(shù)值模擬

壓氣機最后級的流場通過求解雷諾平均Navier-Stokes (RANS)方程及Spalart-Allmaras (SA) 湍流模型[25]確定。采用多重網(wǎng)格、隱式殘值光順、當?shù)貢r間步長等技術(shù)加速計算收斂。最后級進口流動參數(shù)如總壓、總溫、流動角由4.5級壓氣機流場數(shù)值模擬結(jié)果給定，出口處給定輪轂壓力并由徑向平衡方程及其他流動參數(shù)確定壓力徑向分布，轉(zhuǎn)-靜交接面上采用摻混面模型[26]實現(xiàn)流場連續(xù)求解，并采用無反射邊界條件[27]。本文所采用的Mixed-out摻混面模型能夠保證交接面上流量、動量和能量的通量守恒，其詳細數(shù)學形式參考附錄A。

計算采用H-型結(jié)構(gòu)網(wǎng)格，不考慮葉尖間隙影響，S1及S2流面網(wǎng)格如圖2所示。首先在3套不同的網(wǎng)格上進行流向和周向網(wǎng)格無關(guān)性研究，徑向網(wǎng)格數(shù)固定為73，由此確定收斂的流向網(wǎng)格數(shù)；之后在3套不同網(wǎng)格上進行周向和徑向網(wǎng)格無關(guān)性研究，葉面流向網(wǎng)格數(shù)固定，周向和徑向網(wǎng)格增長比為2。網(wǎng)格無關(guān)性研究的具體實現(xiàn)流程可參考本文作者前期工作[19]。表1為近失速工況下周向和徑向網(wǎng)格無關(guān)性研究結(jié)果，包括：歸一化流量、靜子效率、級效率、壓比。由表可知：相對于網(wǎng)格2與網(wǎng)格1的流場偏差，網(wǎng)格3與網(wǎng)格2的流場偏差明顯減小，流量、級效率的相對偏差分別為0.31%和0.46%，壓比偏差接近于零。后續(xù)研究將采用網(wǎng)格2，轉(zhuǎn)子和靜子的網(wǎng)格規(guī)模均為：葉面流向97、周向49、徑向57個網(wǎng)格點。

實際上，為了進一步證實網(wǎng)格無關(guān)性解的有效性，研究中對優(yōu)化后的最后級也采用了相同規(guī)模的密網(wǎng)格(網(wǎng)格3)進行流場數(shù)值模擬，其氣動性能也得到了提升。

圖2 最后級網(wǎng)格Fig.2 Grids for the last stage

表1 網(wǎng)格無關(guān)性解Table 1 Independent grid solutions

2.2 參數(shù)化建模

中弧線、厚度分布、安裝角是葉片的3個重要設計參數(shù)。本文的ADO研究中，只改變最后級靜子葉片的中弧線和安裝角徑向分布，其厚度分布不變。中弧線的變化通過疊加Hicks-Henne型函數(shù)[28]實現(xiàn)，安裝角變化在前期工作中已經(jīng)介紹[15]。選取10%、15%、20%、30%、70%、80%、85%、90%這8個徑向位置，每個徑向位置的安裝角都作為設計參數(shù)，葉型中弧線變化由8個型函數(shù)疊加確定；共有72個設計參數(shù)。

圖3為葉型中弧線變化和安裝角變化的幾何示意圖，其中：δ表示變化量，下標 camber 和 stagger 分別表示中弧線和安裝角。

葉型變化后，在原始網(wǎng)格基礎上采用插值網(wǎng)格擾動技術(shù)[7]重新生成計算網(wǎng)格：

xnew=xold+r(xW,new-xW,old)

(1)

式中：x表示網(wǎng)格坐標；下標 old和 new 分別表示原始和更新后的網(wǎng)格,W表示葉面網(wǎng)格點；r為歸一化的空間距離，在葉面網(wǎng)格點r=1。

圖3 設計參數(shù)變化Fig.3 Design parameter variations

2.3 目標函數(shù)

本文ADO的主要目標是降低流動損失，同時約束出口流量。葉輪機械的流動損失可以用熵增來衡量，因而目標函數(shù)定義為

I=Sgen+Λ|σ-1|

(2)

式中：σ表示設計流量與原始流量的比值；Λ為約束系數(shù)；Sgen表示平均熵增，定義為

(3)

其中：Δs表示每個單元上的熵增，其定義見文獻[15]。

基于伴隨方法的多工況ADO方法在前期工作中已經(jīng)介紹[17]，定義加權(quán)目標函數(shù)為

(4)

式中：I表示氣動參數(shù)；M為工況數(shù)；λi表示權(quán)重且滿足 ∑λi=1。相應地，靈敏度計算表達式為

(5)

式中：Gi表示不同工況下的氣動參數(shù)靈敏度。

本文研究采用最速下降法作為優(yōu)化方法，所需要的靈敏度由伴隨方法計算。接下來將重點介紹伴隨方法的基本原理。

3 伴隨方法

伴隨方法基本原理較為成熟，這里只作簡單介紹。定義I為氣動參數(shù)，一般與流場w和氣動外形相關(guān)，而氣動外形與設計參數(shù)V相關(guān)，即

I=I(w,V)

(6)

氣動參數(shù)關(guān)于設計參數(shù)的靈敏度為

(7)

式中：右端第1、第2項分別表示流場變分、氣動外形變分對靈敏度的影響，靈敏度計算的關(guān)鍵是流場變分的計算。如果采用FDM計算靈敏度，雖然原理簡單，但是計算效率偏低。值得注意的是，流場和設計參數(shù)滿足流體力學方程R，即

R(w,V)=0

(8)

類似于方程(7)，方程(8)可以寫成

(9)

伴隨方法的主要目標是通過將流體力學方程作為約束引入到靈敏度計算中，消除流場變分對靈敏度計算的影響。定義伴隨變量Ψ=[Ψ1,Ψ2,Ψ3,Ψ4,Ψ5]T，將伴隨變量與方程(9)相乘后再與方程(7)相減，有

(10)

若方程(10)右端第1項的系數(shù)為零，即

(11)

則有

(12)

方程(11)稱為方程(8)的伴隨方程，Ψ為伴隨變量，方程(12)為靈敏度計算公式且與流場變分無關(guān)。由此可知：只需求解一次流體力學方程和一次伴隨方程，即可由方程(12)計算氣動參數(shù)I的靈敏度。此外，伴隨方程是線性方程組，可以采用和RANS方程相同的數(shù)值格式進行求解[11,29]。

Giles和Pierce[29]指出了伴隨方程進、出口邊界條件確定方法：伴隨方程的特征線方向與流體力學方程相反。對于多排壓氣機葉片，轉(zhuǎn)子出口、靜子入口的伴隨方程邊界條件需要采用摻混面方法[13,16]確定。本文采用基于Mixed-out摻混面模型的伴隨摻混面方法，首先計算伴隨通量Fa1、Fa2、Fa3、Fa4和Fa5：

(13)

式中：r、p、H分別表示密度、壓力和總焓；ds為面積單元；u=[ux,uθ,ur] 表示速度矢量；n為單位法向矢量。之后由方程(13)可以計算平均伴隨變量Ψa，即

(14)

轉(zhuǎn)-靜交接面上仍采用Giles和Pierce提出的無反射邊界條件。具體實現(xiàn)流程為[10,13]：在亞聲速交接面下游(對應下游靜子入口)，伴隨特征變量Ψa1、Ψa2、Ψa3、Ψa4由內(nèi)場確定，Ψa5由交接面上游確定；在亞聲速交接面上游(對應上游轉(zhuǎn)子出口)，Ψa1、Ψa2、Ψa3、Ψa4由交接面下游確定，Ψa5由內(nèi)場外插確定。確定所有伴隨特征變量后，伴隨變量Ψa1、Ψa2、Ψa3、Ψa4、Ψa5可通過求解線性方程組確定。

4 結(jié)果分析

在基于伴隨方法的ADO研究中，首先需要對伴隨靈敏度進行精度驗證，將基于FDM的靈敏度收斂解作為精確解。圖4給出了30%葉高的靈敏度分布，前8個為中弧線設計參數(shù)靈敏度，第9個為安裝角變化靈敏度；AD表示伴隨方法。

整體上，從葉片前緣(LE)至葉片中部，伴隨靈敏度與FDM靈敏度較為接近；從葉片中部至葉片尾緣(TE)，伴隨靈敏度與FDM存在一定偏差，主要原因是沒有考慮湍流模型方程的伴隨方程求解，而是采用“定渦”假設(Frozen Eddy Viscosity)，即忽略湍流變量變分對靈敏度的影響。實際上，當流動分離較大時，忽略湍流影響會降低伴隨靈敏度精度[30]。圖4中，后3個中弧線設計參數(shù)(6,7,8)正處于靜子葉片分離區(qū)。雖然精度有所下降，但是伴隨靈敏度的分布與FDM基本一致，本文將采用基于“定渦”假設的伴隨方法進行ADO研究。

圖4 靈敏度對比Fig.4 Comparisons of gradients

4.1 近失速點單點優(yōu)化

該4.5級壓氣機最后級靜子角區(qū)存在流動分離，且在接近失速工況時分離區(qū)愈大。單點優(yōu)化設計在設計轉(zhuǎn)速近失速工況下進行，約束函數(shù)系數(shù)Λ= 0.4。圖5給出了等熵效率和流量的收斂曲線，流量最大相對偏差不超過0.1%，約束效果明顯。表2給出了優(yōu)化前后的熵增、流量及靜子和最后級的等熵效率，優(yōu)化后靜子和級效率分別增加了約0.29%、0.15%。

研究中采用相同規(guī)模的密網(wǎng)格(表1中的網(wǎng)格3)對優(yōu)化葉片流場進行數(shù)值模擬，表3給出了優(yōu)化前后密網(wǎng)格的計算結(jié)果。與原始葉片相比，優(yōu)化后靜子和級效率分別增加了0.28%、0.13%。這種較好的網(wǎng)格一致性結(jié)果進一步證實了本文網(wǎng)格無關(guān)性驗證的有效性及優(yōu)化設計結(jié)果的可靠性。

圖6為優(yōu)化前后靜子葉片的氣動外形。由圖6(a) 可知：優(yōu)化設計后自前緣至20%弦長，葉片中弧線曲率增加；自20%弦長至90%弦長，葉片中弧線曲率降低；葉片氣動外形的這種變化將使前緣附近流動加速、來流角與葉片前緣角偏差降低，有利于抑制葉背流動分離。

圖5 單點優(yōu)化熵增和流量收斂曲線Fig.5 Convergence history of entropy and mass flow rate under one point optimization

表2 不同葉片的氣動性能參數(shù)Table 2 Aerodynamic performance parameters of different blades

表3 密網(wǎng)格氣動性能參數(shù)Table 3 Aerodynamic performance parameters of fine grid

為了比較分析優(yōu)化設計對最后級流場的影響，圖7給出了靜子吸力面附近的流線，圖8給出了最后級10%葉高S1流面的馬赫數(shù)云圖。由圖7可知：原始靜子輪轂、機匣角區(qū)均存在流動分離，輪轂角區(qū)尤為明顯；優(yōu)化設計使輪轂流動分離區(qū)朝尾緣和端壁移動，流動分離得到了較好抑制；機匣流動分離區(qū)也得到了較好抑制。由圖8可以更好地看出優(yōu)化設計對靜子葉片葉背流動分離的抑制作用。整體上，對于低速壓氣機，流動分離區(qū)的減小將是效率增加的主要原因。

圖6 優(yōu)化前后靜子葉片氣動外形Fig.6 Stator blade aerodynamic shape before and after optimization

圖7 靜子吸力面附近流線Fig.7 Streamlines near stator suction side

上述優(yōu)化設計能夠改善最后級在指定工況下的氣動性能，在壓氣機全工況范圍內(nèi)對氣動性能的影響更值得關(guān)注。圖9給出了最后級、靜子葉片的等熵效率η特征工作線。雖然在其他工況條件下氣動性能改善并不明顯，但是靜子、級效率在全工況范圍內(nèi)都得到一定的提升。若對最后級的轉(zhuǎn)子葉片也進行優(yōu)化設計，甚至對整個4.5級壓氣機進行優(yōu)化設計，通過調(diào)整上游葉排對下游流動的影響，將能進一步提升壓氣機氣動性能。

值得注意的是：單點優(yōu)化設計的工況點遠離失速邊界，葉背流動分離區(qū)較小，通過控制流動分離獲得效率提升有限。實際上，在前期工作中采用基于本征正交分解的混合模型方法對該型號壓氣機最后級也進行了優(yōu)化設計[31]：工況點距離失速邊界較近，流動分離較劇烈，通過最后級靜子安裝角優(yōu)化使等熵效率提高了1.47%。因此，雖然通過本文的優(yōu)化設計獲得的效率提升較為有限，但是由上述分析可知：優(yōu)化設計結(jié)果是可靠的，在工程中仍具有應用價值。

圖8 10%葉高S1流面馬赫數(shù)云圖Fig.8 Mach number contours at 10% span S1 stream surface

圖9 優(yōu)化前后效率工作線Fig.9 Adiabatic efficiency for both original and optimized blades

4.2 多工況優(yōu)化

多工況優(yōu)化設計研究中，將對最后級靜子葉片在100%設計轉(zhuǎn)速的近失速工況(P1)、140%設計轉(zhuǎn)速的近失速工況(P2)進行中弧線和安裝角的優(yōu)化設計。在每個轉(zhuǎn)速下分別計算其靈敏度，最后由方程(5)計算加權(quán)靈敏度，并由最速下降法實現(xiàn)優(yōu)化設計,權(quán)重λ1=λ2=0.5。

圖10給出了兩個不同轉(zhuǎn)速下(P1和P2)熵增和出口流量的收斂曲線。表4給出了優(yōu)化前后葉片氣動性能。值得注意的是：表4中 P2 轉(zhuǎn)速下原始和優(yōu)化葉片的熵增和流量都統(tǒng)一采用 P1 轉(zhuǎn)速下原始葉片的熵增和流量進行歸一化處理。多工況優(yōu)化設計使靜子效率、級效率在 P1 轉(zhuǎn)速下分別提升了0.29%和0.15%，在 P2 轉(zhuǎn)速下分別提升了0.24%和0.13%，而兩個轉(zhuǎn)速下流量偏差都低于0.15%，表明兩個轉(zhuǎn)速條件下流量都得了較好約束。此外，由表4可知高轉(zhuǎn)速條件下熵增增加，其主要原因是：高轉(zhuǎn)速條件下壓氣機負荷增大，角區(qū)分離加劇，流動損失增大。

圖10 多工況優(yōu)化熵增和流量收斂曲線Fig.10 Convergence history of entropy and mass flow rate under multi-points optimization

圖11給出了P2轉(zhuǎn)速下最后級10%葉高S1流面的馬赫數(shù)云圖。值得注意的是：圖11的最大相對馬赫數(shù)為0.45，遠高于圖8的0.375。對比圖11和圖8發(fā)現(xiàn)：隨著轉(zhuǎn)速的提高，原始葉片葉背分離更嚴重，如圖8(a)和圖11(a)所示；高轉(zhuǎn)速下，通過優(yōu)化設計靜子吸力面的分離區(qū)也得到了較好抑制；此外，相對于圖8，高轉(zhuǎn)速下分離區(qū)的控制效果更加明顯。

上述結(jié)果表明：通過單點、多工況氣動外形優(yōu)化設計能夠提升近失速工況最后級的等熵效率。非設計轉(zhuǎn)速下氣動外形設計優(yōu)化對設計轉(zhuǎn)速下氣動性能的影響需要進一步深入分析。圖12給出了單點、多工況優(yōu)化后靜子、最后級在非設計轉(zhuǎn)速全工況范圍內(nèi)的效率變化。整體上，多工況優(yōu)化設計能使靜子、最后級的等熵效率在非設計轉(zhuǎn)速全工況范圍內(nèi)得到提升。主要原因是：靜子端壁和吸力面角區(qū)的流動分離在不同轉(zhuǎn)速條件下都存在，通過氣動外形設計優(yōu)化抑制流動分離的基本原理是一致的：減小來流角和葉片前緣角的偏差、流動加速等。此外，由于ADO在近失速工況下進行，當工作點遠離失速邊界時，效率提升不斷降低；當工作點向阻塞邊界靠近時，效率提升不斷增加，這一點有待深入研究。

表4 不同轉(zhuǎn)速下的氣動性能參數(shù)Table 4 Aerodynamic performance parameters of different rotation speeds

圖11 140%設計轉(zhuǎn)速下10%葉高S1流面馬赫數(shù)云圖Fig.11 Mach number contours at 10% span S1 stream surface with 140% design rotation speed

圖12 單工況、多工況優(yōu)化效率提升Fig.12 Improvements of adiabatic efficiency for single and multi-points optimization

5 結(jié) 論

1) 采用伴隨方法計算梯度優(yōu)化方法所必需的靈敏度信息，即使對于多設計參數(shù)的ADO問題，優(yōu)化效率仍非常高；對于強湍流問題，忽略湍流影響的伴隨靈敏度精度會稍微下降，但仍可以用于氣動外形設計優(yōu)化工作。

2) 通過單/多工況優(yōu)化設計，最后級靜子葉片中弧線曲率在前緣附近增加，降低了相對來流角并且使葉片吸力面流動加速，對葉片端壁和吸力面角區(qū)的流動分離具有較好的抑制效果。

3) 最后級及靜子葉片效率的提升主要來源于流動分離的抑制，在全工況范圍內(nèi)葉片等熵效率都得到了提升；此外，多工況設計優(yōu)化能夠在不同轉(zhuǎn)速條件下實現(xiàn)葉片全工況范圍內(nèi)氣動性能的提升。

本文的研究工作初步實現(xiàn)了基于黏性伴隨方法的葉輪機械多葉排ADO，對于后續(xù)開展考慮湍流影響的伴隨方法研究、基于伴隨方法的多級ADO研究提供了理論和技術(shù)基礎。

附錄A

Mixed-out是葉輪機械多排流場數(shù)值模擬的一種主要的摻混面方法,該方法的基本表達式為

(A1)

式中：H表示總焓；p、ρ分別表示壓力、密度；u、n、ds分別表示速度矢量、單位外法線矢量、單元面積矢量。

方程(A1)經(jīng)過整理后，可以簡化成一個關(guān)于平均壓力的二次方程。確定平均壓力之后，再由方程(A1)即可計算其他平均流動參數(shù)：

(A2)

式中：下標x、θ、r分別表示軸向、周向和徑向；平均壓力表達式中根號前面的“+”對應亞聲速流動解，根號前面的“-”對應超聲速流動解；參數(shù)a、b、c分別表示為

(A3)