基于伴隨方法的單級低速壓氣機氣動設計優化

羅佳奇，楊婧

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

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

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

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

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

1 壓氣機初步設計

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

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

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

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

第4步 采用簡單徑向平衡方程確定輪轂和機匣的相關參數。

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

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

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

2 氣動設計優化

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

2.1 流場數值模擬

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

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

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

圖2 最后級網格Fig.2 Grids for the last stage

表1 網格無關性解Table 1 Independent grid solutions

2.2 參數化建模

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

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

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

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

(1)

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

圖3 設計參數變化Fig.3 Design parameter variations

2.3 目標函數

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

I=Sgen+Λ|σ-1|

(2)

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

(3)

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

基于伴隨方法的多工況ADO方法在前期工作中已經介紹[17]，定義加權目標函數為

(4)

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

(5)

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

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

3 伴隨方法

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

I=I(w,V)

(6)

氣動參數關于設計參數的靈敏度為

(7)

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

R(w,V)=0

(8)

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

(9)

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

(10)

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

(11)

則有

(12)

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

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

(13)

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

(14)

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

4 結果分析

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

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

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

4.1 近失速點單點優化

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

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

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

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

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

表3 密網格氣動性能參數Table 3 Aerodynamic performance parameters of fine grid

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

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

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

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

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

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

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

4.2 多工況優化

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

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

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

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

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

表4 不同轉速下的氣動性能參數Table 4 Aerodynamic performance parameters of different rotation speeds

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

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

5 結 論

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

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

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

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

附錄A

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

(A1)

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

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

(A2)

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

(A3)