張效溥,王仲林,田 杰,2,3,陳 勇,2,3,歐陽華,2,3
(1.上海交通大學機械與動力工程學院,上海200240;2.燃氣輪機與民用航空發動機教育部工程研究中心,上海200240;3.先進航空發動機協同創新中心,北京100191)
航空發動機適航認證條例CCAR 33.83 規定:每型發動機必須進行振動測試,以確定可能受機械或空氣動力導致激振的部件的振動特性在整個聲明的飛行包線范圍內是可接受的[1]。作為航空歷史上銷量最高的CFM56 系列發動機,雖然通過了適航認證且運行了上億小時,但仍在2016 年[2]和2018 年[3]分別發生了2 次嚴重的發動機事故。初步報告指出事故原因是由于風扇葉片的高周疲勞導致在靠近葉片根部發生失效斷裂,葉片飛出后造成非包容性損傷。可見風扇葉片的高周疲勞對發動機的安全性有很大影響。
國內外學者對航空發動機風扇葉片的振動特性和高周疲勞問題進行了大量研究。Cowles[4]總結了航空發動機高周疲勞產生的原因和研究方法;Ewins[5]歸納了旋轉葉輪機械的振動來源,重點關注共振問題的避免和控制技術;Ni 等[6]通過仿真和試驗方法研究了多模態激勵下風扇葉片的應力特性;馬艷紅等[7]分析了離心力產生的彎矩和扭矩與氣動力產生的彎矩復合作用對葉片靜應力的影響,結果表明可以通過調節葉片積疊線的彎掠構型,使得附加彎矩與扭矩載荷相互抵消,降低葉片最大靜應力;Deng 等[8]、Astrua 等[9]、張薊欣等[10]以靜應力和共振裕度為目標函數對葉片的積疊線進行優化,得到具有足夠安全系數的葉片設計;Li 等[11]、Wang 等[12]、Samad 等[13-14]、Myoren 等[15]學者利用CFD 對葉片的積疊線彎掠造型進行優化,得到具有更好氣動性能的葉片設計;Kou 等[16]和Simmons等[17]研究了氣動載荷對葉片靜頻和動頻的影響及共振下的應力分布,結果顯示氣動載荷對振動應力沒有明顯影響。上述基于葉片積疊線的優化研究主要考慮氣動性能參數,將靜應力和共振裕度作為限制條件,而對葉片振動應力及振動特性考慮較少。
本文通過航空發動機風扇葉片重心積疊線周向構型的優化設計,改善風扇的強度和振動特性。基于Kriging 代理模型和微種群遺傳算法,建立了航空發動機風扇“參數化建模-有限元仿真-強度和振動特性優化”一體化平臺。基于該平臺對某寬弦風扇葉片進行優化設計,并對比分析原型風扇和優化后風扇的振動特性。
本文采用Kriging 代理模型和微種群遺傳算法,通過Matlab 程序和批處理腳本,自動調用UG/NX 和ANSYS-APDL 等軟件,實現了優化流程的自動化運行,建立了航空發動機風扇“參數化建模—有限元仿真—強度和振動特性優化”的一體化平臺,如圖1 所示。
從圖中可見,首先利用均勻設計法構造設計變量矩陣,通過參數化建模、自動網格劃分、有限元仿真及后處理得到完整的樣本空間;其次建立代理模型,利用遺傳算法搜索并預測使目標函數最小的設計變量組合;再次通過種群迭代直到目標函數達到收斂準則;最終將滿足設計要求的最優解生成葉型,優化過程結束。
在計算資源有限的情況下,代理模型使用少量的數據點即可對整個樣本空間進行較高精度的預測[13]。本文采用Kriging 插值模型為代理模型,將微種群遺傳算法作為優化算法對風扇葉形進行優化設計,大大降低了優化所需要的計算資源和時間成本。
1.1.1 Kriging 模型的構建
Kriging 模型通過將樣本函數線性加權,得到響應值與設計量之間的近似插值函數,對最優設計變量進行預測,插值函數的基本形式為

式中:f(x)為未知的全局近似值;Z(x)為均方差為σ、均值為0 的高斯隨機函數。
空間內不同變量之間的相關性由協方差矩陣函數表示為

式中:R 為對角線對稱的相關矩陣;C 為矩陣中的相關函數;θ 為待訓練參數。本文相關系數選取各向同性的高斯指數模型,取光滑因子p=2,其表達式為

根據插值函數無偏性的特點,推導出未知點處預估函數值為

式中:r 為空間內未知點與樣本點之間的相關向量,β 及β 的標準差為
1.1.2 Kriging 模型訓練
Kriging 插值模型優劣的關鍵在于參數的選取,本文采取最大似然估計法作為訓練的依據,最終將問題轉化為使下式最大

由于無法直接求解θ,因此需要采用合適的優化算法。比較常見的算法為模式搜索算法,但是該算法對初始設計變量的依賴性很高[18],會影響Kriging 模型的訓練結果。本文采用微種群遺傳算法,相比傳統的遺傳算法,微種群遺傳算法的種群較小,遺傳操作采取錦標賽選擇和均勻交叉算子,并且略去了變異操作,因此擁有更高的收斂速度和全局尋優能力[19-20]。
1.1.3 最優解的求取
本文采取EI 函數來求得最優解,其定義為

式中:Φ 和Ψ 分別為標準正態分布的分布函數和密度函數;和分別為x 處的代理模型估計值及標準差;ymin為當前樣本庫中的最優函數值。
EI 函數可以有效地表征目標函數的改善量,因此在約束范圍內求解出使EI 函數達到最大值的x 即為預測最優解。與代理模型的訓練相似,本文采用微種群遺傳算法求解x,有效克服了傳統懲罰函數法中收斂性較差、懲罰因子難以選取的問題[21],大幅度減小了計算量。
為了在優化過程中實現自動建模和劃分網格,需要將葉片模型參數化處理。首先將葉片以流道線為界分成葉型和葉根2 部分,僅對葉型部分進行優化。為了提取原型葉片的截面和積疊線,以風扇旋轉中心為軸,利用圓柱面和圓錐面分別截取葉型,得到16 組葉型截面,通過積分法求得重心積疊線的坐標。其次,選取第7、11 和15 截面的重心作為控制點[10,22],如圖2所示。這樣既保證葉型與伸根段的光順連接,又避免對葉尖造型造成較大影響。優化過程中改變控制點的坐標,通過樣條插值得到新的積疊線,通過曲面網格法生成新葉型。最后,將新葉型與葉根進行光滑性拼接后得到新葉片模型。在ANSYS-APDL 中通過掃掠生成結構化網格,包括18976 個單元和25602 個節點,得到網格質量良好的有限元模型,如圖3 所示。在有限元計算中,對榫頭上表面節點施加法向約束,對榫頭前端面施加軸向約束。

圖2 葉片型線的提取和積疊線的生成

圖3 葉片有限元模型及網格劃分
1.3.1 設計參數
文獻[7]中的結果表明:葉片積疊線的彎掠構型使得附加彎矩與扭矩載荷相互抵消,從而降低葉片最大靜應力。雖然風扇葉片的氣動載荷和靜態變形存在耦合關系,但相比于葉片的掠型,彎型的小偏移量對氣動參數影響較小,氣動參數的微小變化也不會對葉片的振動應力產生明顯影響。文獻[23]對比了單獨離心載荷與離心和氣動載荷同時作用下風扇葉片的固有頻率、共振裕度及振動應力,結果顯示:氣動載荷對葉片的振動應力等特性影響較小。為此本文僅對葉片的周向積疊線進行小范圍優化,且不考慮氣動載荷的影響。
在葉片參數化建模中截取16 個截面,如果將所有截面的周向坐標作為設計參數,很可能導致葉片表面呈現波紋狀不光順的造型。因此設計參數選取積疊線上3 個截面重心作為控制點,每個截面的周向坐標作為設計參數,通過3 次樣條生成新的積疊線。為了保證與葉根的光順連接且葉尖位置造型上與原始葉片沒有太大的區別,靠近葉根的1~5 號截面重心設為固定點,16 號截面和15 號截面保持相同的位移。這樣的設計參數可保證盡可能減小對氣動性能的影響。
1.3.2 目標參數
大涵道比渦扇發動機風扇葉片由于尺寸較大且轉速較高,在工作中會產生很大的離心應力,這對葉片的結構強度提出很高要求。現有涉及葉片重心積疊線優化的文獻中均將最大靜應力作為約束條件或目標參數,但由于風扇葉片在工作中受到來自發動機進口的多種激勵,如流場畸變、地面吸渦等,還受到流道內渦流和湍流的激勵。這些來自氣動載荷的激勵通常具有一定的頻率特征,使風扇葉片產生劇烈的振動甚至是顫振,對風扇的安全性和可靠性有很大影響。因此,為了使優化后葉片的振動特性滿足要求,除最大靜應力Smax之外,還引入模態振動應變能密度指數作為優化目標。
應變能密度可表征應變狀態。應變能密度準則常被用來判斷彈塑性材料的失效斷裂及裂紋擴展[24-26]。本文選擇的是應變能密度指數κ,其定義為葉片表面最大應變能密度ρmax,SEDe與平均應變能密度的比值

式中:上標m 表示第m 階模態。
應變能密度定義為

式中:εe、σe和Ve分別為單元的應變、應力和體積。
指數κmax的相對值表征了不同葉片設計的應變能密度的集中程度。葉片的κmax值越小,其應變能分布越均勻,最大應變值相對較低,葉片的疲勞壽命更長。由于優化過程中的種群最優解是通過目標參數值對比得到的,所以應變能密度指數適用于優化過程中表征葉片的振動特性。
在發動機工作中,風扇葉片的第1 階彎曲模態會受到側風、進氣畸變、吸渦等激勵,而且顫振也經常發生在第1 階模態。第1 階模態的振幅和能量都明顯高于其他階次,所以在工程中更關注第1 階模態的振動。本文在優化中將第1 階模態的應變能密度指數κ1max作為代表振動特性的目標參數。考慮到靜應力和應變能密度指數的數量級和量綱不同,分別通過原始模型的參數值對目標參數進行歸一化處理,再加權得到最終的目標參數f,本文選取權重系數c1=c2=1。

基于上述優化流程,將某寬弦風扇葉片作為優化對象,進行周向積疊線優化,其設計及材料參數見表1。在有限元靜力分析中施加轉速載荷,模態分析和諧響應分析分別采用Block Lanczos 法和完全法,均考慮了離心力產生的預應力。

表1 原型葉片的材料及設計參數
7、11 和15 號截面重心的周向偏移量分別記為x1、x2和x3。文獻[27]表明:葉片彎向構型會促使附面層低能流體向葉中移動,對喘振裕度、氣動噪聲等產生一定影響。因此,為了保證葉型的光滑性并使其氣動性能影響盡可能小,設置其變化范圍分別為-2~2 mm、-3.5~3.5 mm 和-5~5 mm。則該優化問題可表示為

使用代理模型優化算法,先分別以靜應力(參數A)和應變能密度指數(參數B)為單目標函數進行優化,再對二者加權的復合目標函數(C)進行優化。經過50 次迭代后,優化結束。原型與C 參數優化后的葉型積疊線對比如圖4 所示。從圖中可見,3 個控制點都向吸力面偏移,中間的控制點偏移比兩端的大,葉片整體呈“C”型。樣本庫內最優值隨著迭代次數變化如圖5 所示。

圖4 優化前、后的重心積疊線

圖5 樣本庫最優值與迭代次數的變化
優化后的設計參數見表2。從表中可見,單目標參數優化的結果在改善了目標參數的同時,導致另一參數的惡化,而雙目標參數加權的優化結果使得2 個參數都得到了改善。最大靜應力相比原型葉片減小了5.45%,第1 階應變能密度指數減小了5.94%。

表2 優化葉片與原型葉片參數對比
優化前、后的等效靜應力如圖6 所示。從圖中可見,原型葉片的2 個高應力區域分別位于葉根前緣和葉型中部。靜應力的最大值位于吸力面的葉根前緣,葉片中部的靜應力具有寬弦風扇葉片典型的“靶心”式分布。相比之下,優化葉片的最大應力也出現在葉根前緣,但是數值減小了5.45%,葉片安全系數則由1.29 提高到1.37。同時,葉片中部的高應力區向葉尖方向移動;葉根前緣處較高應力區的應力水平也得到了較明顯的改善。

圖6 優化前、后葉片Von Mises 靜應力
優化前、后的第1 階模態應變能密度分布如圖7所示。從圖中可見,原形在葉片中心有1 個較高幅值區域。而在優化葉型中,該區域幅值有所減小。同時,原型葉片的應變能密度集中于前緣,而優化葉片高應變能密度區域分成了2 個應變能密度相對較低的區域,最大值點位置沒有明顯變化,且峰值降低了5.94%。說明葉片的應變能密度分布得到了改善。

圖7 優化前、后1 階應變能密度指數

圖8 優化前、后前5 階振型
除了作為目標參數的靜應力和應變能指數以外,風扇葉片的固有頻率及振型、共振裕度和weaklink 也是衡量振動特性的重要指標。葉片前5 階振型如圖8所示。圖中從左至右依次為第1 階周向彎曲模態、第1 階軸向彎曲模態、第1 階扭轉模態、第2 階彎曲模態及第2 階彎扭組合模態。結果顯示優化前、后的前4 階振型基本不變,第5 階振型在后緣處有一定的差別,可見高階復雜振型對積疊線的變化更敏感。
原形葉片和優化后葉片在不同轉速下的固有頻率如圖9 所示。從圖中可見,二者在前6 階的固有頻率基本一致。由坎貝爾圖可以得到葉片在各轉速下的共振裕度。對于工程上較關注的第1 階模態,從轉速為60%~105%內,優化前、后葉片的共振裕度都保持在20%以上,相互差別極小,對比見表3。因此優化設計沒有對共振裕度產生明顯影響。
振動應力裕度基于Goodman 曲線,在考慮靜應力的基礎上可評估振動應力到許用應力的裕量,是評價葉片高周疲勞特性的重要指標。Goodman 曲線如圖10 所示,直線上方的區域代表高周疲勞的危險區域。考慮葉片材料屬性和加工的不確定度,需要在設計中留有足夠應力裕度。振動應力裕度M 定義為考慮靜應力下,實際振動應力σv與許用振動應力σA的相對差值

圖9 優化前、后葉片的故有頻率

表3 優化前、后轉速為60%~110%時第1 階頻率裕度

為了對比優化前、后的應力裕度,在諧響應分析中對葉片施加周期性全局加速度載荷,載荷頻率為葉片的第1 階固有頻率,計算得到振動應力從而求得第1 階模態的應力裕度。結果顯示優化后的第1 階振動應力裕度從66.81%提高到70.46%。
在應力裕度的基礎上,對葉片高周疲勞的失效位置,即weaklink 進行評估。通過仿真計算得到每個節點的靜態應力和模態應力,再通過對模態應力的等比例縮放,使得剛好有1 個節點在Goodman 曲線上。這1 點就代表了某個固有模態階次振動下最先發生高周疲勞失效的位置,即weaklink 點。優化前、后葉片的第1 階weaklink 結果如圖11 所示。從圖中可見,原型葉片有3 個節點(A、B 和C)同時接近Goodman 曲線,均位于前緣根部。優化后葉片的weaklink 點仍然為這3 個點。由于優化葉片靜應力水平較低,3 個點在坐標圖上明顯向左偏移,且對應的第1 階振動應力裕度均有所增大,說明葉片的高周疲勞特性得到了改善。

圖10 修正的Goodman 曲線[4,28]

圖11 優化前、后葉片的Weaklink
本文對風扇葉片重心積疊線彎曲構型進行了參數化建模,使用Kriging 代理模型和微種群遺傳算法搭建了風扇葉片強度和振動特性的一體化平臺,對某寬弦風扇葉片的積疊線進行優化設計,得到了如下結論:
(1)基于重心積疊線周向構型對風扇葉片進行參數化建模和自動化網格劃分。結合Kriging 代理模型和微種群遺傳算法,綜合考慮靜應力和應變能密度指數為優化目標函數,建立了航空發動機風扇“參數化建模—有限元仿真—強度和振動特性優化”一體化平臺。
(2)基于該平臺對某型寬弦風扇葉片改型后的直葉片進行了優化,得到優化設計的各項目標參數均有改善。優化后的葉片靜應力在葉片中心高應力區向葉尖偏移,葉根前緣的應力集中得到改善,且最大靜應力比原型葉片的減小了5.45%,第1 階應變能密度指數減小了5.94%,相比于單目標優化結果更加全面有效。
(3)對比優化前、后葉片的振動特性發現,葉片的振型、固有頻率和weaklink 點沒有明顯變化,第1 階振動應力裕度從66.81%提高到了70.46%。說明優化流程有效地改善了葉片的強度和高周疲勞特性,且對其他振動特性沒有負面影響。