張效溥，王仲林，田 杰，2，3，陳 勇，2，3，歐陽華，2，3

（1.上海交通大學機械與動力工程學院，上海200240；2.燃氣輪機與民用航空發(fā)動機教育部工程研究中心，上海200240；3.先進航空發(fā)動機協(xié)同創(chuàng)新中心，北京100191）

0 引言

航空發(fā)動機適航認證條例CCAR 33.83 規(guī)定：每型發(fā)動機必須進行振動測試，以確定可能受機械或空氣動力導致激振的部件的振動特性在整個聲明的飛行包線范圍內(nèi)是可接受的[1]。作為航空歷史上銷量最高的CFM56 系列發(fā)動機，雖然通過了適航認證且運行了上億小時，但仍在2016 年[2]和2018 年[3]分別發(fā)生了2 次嚴重的發(fā)動機事故。初步報告指出事故原因是由于風扇葉片的高周疲勞導致在靠近葉片根部發(fā)生失效斷裂，葉片飛出后造成非包容性損傷?？梢婏L扇葉片的高周疲勞對發(fā)動機的安全性有很大影響。

國內(nèi)外學者對航空發(fā)動機風扇葉片的振動特性和高周疲勞問題進行了大量研究。Cowles[4]總結(jié)了航空發(fā)動機高周疲勞產(chǎn)生的原因和研究方法；Ewins[5]歸納了旋轉(zhuǎn)葉輪機械的振動來源，重點關(guān)注共振問題的避免和控制技術(shù)；Ni 等[6]通過仿真和試驗方法研究了多模態(tài)激勵下風扇葉片的應(yīng)力特性；馬艷紅等[7]分析了離心力產(chǎn)生的彎矩和扭矩與氣動力產(chǎn)生的彎矩復合作用對葉片靜應(yīng)力的影響，結(jié)果表明可以通過調(diào)節(jié)葉片積疊線的彎掠構(gòu)型，使得附加彎矩與扭矩載荷相互抵消，降低葉片最大靜應(yīng)力；Deng 等[8]、Astrua 等[9]、張薊欣等[10]以靜應(yīng)力和共振裕度為目標函數(shù)對葉片的積疊線進行優(yōu)化，得到具有足夠安全系數(shù)的葉片設(shè)計；Li 等[11]、Wang 等[12]、Samad 等[13-14]、Myoren 等[15]學者利用CFD 對葉片的積疊線彎掠造型進行優(yōu)化，得到具有更好氣動性能的葉片設(shè)計；Kou 等[16]和Simmons等[17]研究了氣動載荷對葉片靜頻和動頻的影響及共振下的應(yīng)力分布，結(jié)果顯示氣動載荷對振動應(yīng)力沒有明顯影響。上述基于葉片積疊線的優(yōu)化研究主要考慮氣動性能參數(shù)，將靜應(yīng)力和共振裕度作為限制條件，而對葉片振動應(yīng)力及振動特性考慮較少。

本文通過航空發(fā)動機風扇葉片重心積疊線周向構(gòu)型的優(yōu)化設(shè)計，改善風扇的強度和振動特性?；贙riging 代理模型和微種群遺傳算法，建立了航空發(fā)動機風扇“參數(shù)化建模-有限元仿真-強度和振動特性優(yōu)化”一體化平臺?；谠撈脚_對某寬弦風扇葉片進行優(yōu)化設(shè)計，并對比分析原型風扇和優(yōu)化后風扇的振動特性。

1 優(yōu)化流程及建模

本文采用Kriging 代理模型和微種群遺傳算法，通過Matlab 程序和批處理腳本，自動調(diào)用UG/NX 和ANSYS-APDL 等軟件，實現(xiàn)了優(yōu)化流程的自動化運行，建立了航空發(fā)動機風扇“參數(shù)化建?！邢拊抡妗獜姸群驼駝犹匦詢?yōu)化”的一體化平臺，如圖1 所示。

從圖中可見，首先利用均勻設(shè)計法構(gòu)造設(shè)計變量矩陣，通過參數(shù)化建模、自動網(wǎng)格劃分、有限元仿真及后處理得到完整的樣本空間；其次建立代理模型，利用遺傳算法搜索并預(yù)測使目標函數(shù)最小的設(shè)計變量組合；再次通過種群迭代直到目標函數(shù)達到收斂準則；最終將滿足設(shè)計要求的最優(yōu)解生成葉型，優(yōu)化過程結(jié)束。

1.1 基于Kriging 模型的優(yōu)化算法

在計算資源有限的情況下，代理模型使用少量的數(shù)據(jù)點即可對整個樣本空間進行較高精度的預(yù)測[13]。本文采用Kriging 插值模型為代理模型，將微種群遺傳算法作為優(yōu)化算法對風扇葉形進行優(yōu)化設(shè)計，大大降低了優(yōu)化所需要的計算資源和時間成本。

1.1.1 Kriging 模型的構(gòu)建

Kriging 模型通過將樣本函數(shù)線性加權(quán)，得到響應(yīng)值與設(shè)計量之間的近似插值函數(shù)，對最優(yōu)設(shè)計變量進行預(yù)測，插值函數(shù)的基本形式為

式中：f(x)為未知的全局近似值；Z(x)為均方差為σ、均值為0 的高斯隨機函數(shù)。

空間內(nèi)不同變量之間的相關(guān)性由協(xié)方差矩陣函數(shù)表示為

式中：R 為對角線對稱的相關(guān)矩陣；C 為矩陣中的相關(guān)函數(shù)；θ 為待訓練參數(shù)。本文相關(guān)系數(shù)選取各向同性的高斯指數(shù)模型，取光滑因子p=2，其表達式為

根據(jù)插值函數(shù)無偏性的特點，推導出未知點處預(yù)估函數(shù)值為

式中：r 為空間內(nèi)未知點與樣本點之間的相關(guān)向量，β 及β 的標準差為

1.1.2 Kriging 模型訓練

Kriging 插值模型優(yōu)劣的關(guān)鍵在于參數(shù)的選取，本文采取最大似然估計法作為訓練的依據(jù)，最終將問題轉(zhuǎn)化為使下式最大

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

1.1.3 最優(yōu)解的求取

本文采取EI 函數(shù)來求得最優(yōu)解，其定義為

式中:Φ 和Ψ 分別為標準正態(tài)分布的分布函數(shù)和密度函數(shù)；和分別為x 處的代理模型估計值及標準差；ymin為當前樣本庫中的最優(yōu)函數(shù)值。

EI 函數(shù)可以有效地表征目標函數(shù)的改善量，因此在約束范圍內(nèi)求解出使EI 函數(shù)達到最大值的x 即為預(yù)測最優(yōu)解。與代理模型的訓練相似，本文采用微種群遺傳算法求解x，有效克服了傳統(tǒng)懲罰函數(shù)法中收斂性較差、懲罰因子難以選取的問題[21]，大幅度減小了計算量。

1.2 風扇葉片參數(shù)化建模

為了在優(yōu)化過程中實現(xiàn)自動建模和劃分網(wǎng)格，需要將葉片模型參數(shù)化處理。首先將葉片以流道線為界分成葉型和葉根2 部分，僅對葉型部分進行優(yōu)化。為了提取原型葉片的截面和積疊線，以風扇旋轉(zhuǎn)中心為軸，利用圓柱面和圓錐面分別截取葉型，得到16 組葉型截面，通過積分法求得重心積疊線的坐標。其次，選取第7、11 和15 截面的重心作為控制點[10，22]，如圖2所示。這樣既保證葉型與伸根段的光順連接，又避免對葉尖造型造成較大影響。優(yōu)化過程中改變控制點的坐標，通過樣條插值得到新的積疊線，通過曲面網(wǎng)格法生成新葉型。最后，將新葉型與葉根進行光滑性拼接后得到新葉片模型。在ANSYS-APDL 中通過掃掠生成結(jié)構(gòu)化網(wǎng)格，包括18976 個單元和25602 個節(jié)點，得到網(wǎng)格質(zhì)量良好的有限元模型，如圖3 所示。在有限元計算中，對榫頭上表面節(jié)點施加法向約束，對榫頭前端面施加軸向約束。

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

圖3 葉片有限元模型及網(wǎng)格劃分

1.3 設(shè)計參數(shù)與目標參數(shù)

1.3.1 設(shè)計參數(shù)

文獻[7]中的結(jié)果表明：葉片積疊線的彎掠構(gòu)型使得附加彎矩與扭矩載荷相互抵消，從而降低葉片最大靜應(yīng)力。雖然風扇葉片的氣動載荷和靜態(tài)變形存在耦合關(guān)系，但相比于葉片的掠型，彎型的小偏移量對氣動參數(shù)影響較小，氣動參數(shù)的微小變化也不會對葉片的振動應(yīng)力產(chǎn)生明顯影響。文獻[23]對比了單獨離心載荷與離心和氣動載荷同時作用下風扇葉片的固有頻率、共振裕度及振動應(yīng)力，結(jié)果顯示：氣動載荷對葉片的振動應(yīng)力等特性影響較小。為此本文僅對葉片的周向積疊線進行小范圍優(yōu)化，且不考慮氣動載荷的影響。

在葉片參數(shù)化建模中截取16 個截面，如果將所有截面的周向坐標作為設(shè)計參數(shù)，很可能導致葉片表面呈現(xiàn)波紋狀不光順的造型。因此設(shè)計參數(shù)選取積疊線上3 個截面重心作為控制點，每個截面的周向坐標作為設(shè)計參數(shù)，通過3 次樣條生成新的積疊線。為了保證與葉根的光順連接且葉尖位置造型上與原始葉片沒有太大的區(qū)別，靠近葉根的1～5 號截面重心設(shè)為固定點，16 號截面和15 號截面保持相同的位移。這樣的設(shè)計參數(shù)可保證盡可能減小對氣動性能的影響。

1.3.2 目標參數(shù)

大涵道比渦扇發(fā)動機風扇葉片由于尺寸較大且轉(zhuǎn)速較高，在工作中會產(chǎn)生很大的離心應(yīng)力，這對葉片的結(jié)構(gòu)強度提出很高要求。現(xiàn)有涉及葉片重心積疊線優(yōu)化的文獻中均將最大靜應(yīng)力作為約束條件或目標參數(shù)，但由于風扇葉片在工作中受到來自發(fā)動機進口的多種激勵，如流場畸變、地面吸渦等，還受到流道內(nèi)渦流和湍流的激勵。這些來自氣動載荷的激勵通常具有一定的頻率特征，使風扇葉片產(chǎn)生劇烈的振動甚至是顫振，對風扇的安全性和可靠性有很大影響。因此，為了使優(yōu)化后葉片的振動特性滿足要求，除最大靜應(yīng)力Smax之外，還引入模態(tài)振動應(yīng)變能密度指數(shù)作為優(yōu)化目標。

應(yīng)變能密度可表征應(yīng)變狀態(tài)。應(yīng)變能密度準則常被用來判斷彈塑性材料的失效斷裂及裂紋擴展[24-26]。本文選擇的是應(yīng)變能密度指數(shù)κ，其定義為葉片表面最大應(yīng)變能密度ρmax,SEDe與平均應(yīng)變能密度的比值

式中：上標m 表示第m 階模態(tài)。

應(yīng)變能密度定義為

式中：εe、σe和Ve分別為單元的應(yīng)變、應(yīng)力和體積。

指數(shù)κmax的相對值表征了不同葉片設(shè)計的應(yīng)變能密度的集中程度。葉片的κmax值越小，其應(yīng)變能分布越均勻，最大應(yīng)變值相對較低，葉片的疲勞壽命更長。由于優(yōu)化過程中的種群最優(yōu)解是通過目標參數(shù)值對比得到的，所以應(yīng)變能密度指數(shù)適用于優(yōu)化過程中表征葉片的振動特性。

在發(fā)動機工作中，風扇葉片的第1 階彎曲模態(tài)會受到側(cè)風、進氣畸變、吸渦等激勵，而且顫振也經(jīng)常發(fā)生在第1 階模態(tài)。第1 階模態(tài)的振幅和能量都明顯高于其他階次，所以在工程中更關(guān)注第1 階模態(tài)的振動。本文在優(yōu)化中將第1 階模態(tài)的應(yīng)變能密度指數(shù)κ1max作為代表振動特性的目標參數(shù)。考慮到靜應(yīng)力和應(yīng)變能密度指數(shù)的數(shù)量級和量綱不同，分別通過原始模型的參數(shù)值對目標參數(shù)進行歸一化處理，再加權(quán)得到最終的目標參數(shù)f，本文選取權(quán)重系數(shù)c1=c2=1。

2 優(yōu)化設(shè)計算例

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

表1 原型葉片的材料及設(shè)計參數(shù)

7、11 和15 號截面重心的周向偏移量分別記為x1、x2和x3。文獻[27]表明：葉片彎向構(gòu)型會促使附面層低能流體向葉中移動，對喘振裕度、氣動噪聲等產(chǎn)生一定影響。因此，為了保證葉型的光滑性并使其氣動性能影響盡可能小，設(shè)置其變化范圍分別為-2～2 mm、-3.5～3.5 mm 和-5～5 mm。則該優(yōu)化問題可表示為

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

圖4 優(yōu)化前、后的重心積疊線

圖5 樣本庫最優(yōu)值與迭代次數(shù)的變化

優(yōu)化后的設(shè)計參數(shù)見表2。從表中可見，單目標參數(shù)優(yōu)化的結(jié)果在改善了目標參數(shù)的同時，導致另一參數(shù)的惡化，而雙目標參數(shù)加權(quán)的優(yōu)化結(jié)果使得2 個參數(shù)都得到了改善。最大靜應(yīng)力相比原型葉片減小了5.45%，第1 階應(yīng)變能密度指數(shù)減小了5.94%。

表2 優(yōu)化葉片與原型葉片參數(shù)對比

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

圖6 優(yōu)化前、后葉片Von Mises 靜應(yīng)力

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

圖7 優(yōu)化前、后1 階應(yīng)變能密度指數(shù)

圖8 優(yōu)化前、后前5 階振型

除了作為目標參數(shù)的靜應(yīng)力和應(yīng)變能指數(shù)以外，風扇葉片的固有頻率及振型、共振裕度和weaklink 也是衡量振動特性的重要指標。葉片前5 階振型如圖8所示。圖中從左至右依次為第1 階周向彎曲模態(tài)、第1 階軸向彎曲模態(tài)、第1 階扭轉(zhuǎn)模態(tài)、第2 階彎曲模態(tài)及第2 階彎扭組合模態(tài)。結(jié)果顯示優(yōu)化前、后的前4 階振型基本不變，第5 階振型在后緣處有一定的差別，可見高階復雜振型對積疊線的變化更敏感。

原形葉片和優(yōu)化后葉片在不同轉(zhuǎn)速下的固有頻率如圖9 所示。從圖中可見，二者在前6 階的固有頻率基本一致。由坎貝爾圖可以得到葉片在各轉(zhuǎn)速下的共振裕度。對于工程上較關(guān)注的第1 階模態(tài)，從轉(zhuǎn)速為60%～105%內(nèi)，優(yōu)化前、后葉片的共振裕度都保持在20%以上，相互差別極小，對比見表3。因此優(yōu)化設(shè)計沒有對共振裕度產(chǎn)生明顯影響。

振動應(yīng)力裕度基于Goodman 曲線，在考慮靜應(yīng)力的基礎(chǔ)上可評估振動應(yīng)力到許用應(yīng)力的裕量，是評價葉片高周疲勞特性的重要指標。Goodman 曲線如圖10 所示，直線上方的區(qū)域代表高周疲勞的危險區(qū)域?？紤]葉片材料屬性和加工的不確定度，需要在設(shè)計中留有足夠應(yīng)力裕度。振動應(yīng)力裕度M 定義為考慮靜應(yīng)力下，實際振動應(yīng)力σv與許用振動應(yīng)力σA的相對差值

圖9 優(yōu)化前、后葉片的故有頻率

表3 優(yōu)化前、后轉(zhuǎn)速為60%~110%時第1 階頻率裕度

為了對比優(yōu)化前、后的應(yīng)力裕度，在諧響應(yīng)分析中對葉片施加周期性全局加速度載荷，載荷頻率為葉片的第1 階固有頻率，計算得到振動應(yīng)力從而求得第1 階模態(tài)的應(yīng)力裕度。結(jié)果顯示優(yōu)化后的第1 階振動應(yīng)力裕度從66.81%提高到70.46%。

在應(yīng)力裕度的基礎(chǔ)上，對葉片高周疲勞的失效位置，即weaklink 進行評估。通過仿真計算得到每個節(jié)點的靜態(tài)應(yīng)力和模態(tài)應(yīng)力，再通過對模態(tài)應(yīng)力的等比例縮放，使得剛好有1 個節(jié)點在Goodman 曲線上。這1 點就代表了某個固有模態(tài)階次振動下最先發(fā)生高周疲勞失效的位置，即weaklink 點。優(yōu)化前、后葉片的第1 階weaklink 結(jié)果如圖11 所示。從圖中可見，原型葉片有3 個節(jié)點（A、B 和C）同時接近Goodman 曲線，均位于前緣根部。優(yōu)化后葉片的weaklink 點仍然為這3 個點。由于優(yōu)化葉片靜應(yīng)力水平較低，3 個點在坐標圖上明顯向左偏移，且對應(yīng)的第1 階振動應(yīng)力裕度均有所增大，說明葉片的高周疲勞特性得到了改善。

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

圖11 優(yōu)化前、后葉片的Weaklink

3 結(jié)論

本文對風扇葉片重心積疊線彎曲構(gòu)型進行了參數(shù)化建模，使用Kriging 代理模型和微種群遺傳算法搭建了風扇葉片強度和振動特性的一體化平臺，對某寬弦風扇葉片的積疊線進行優(yōu)化設(shè)計，得到了如下結(jié)論：

（1）基于重心積疊線周向構(gòu)型對風扇葉片進行參數(shù)化建模和自動化網(wǎng)格劃分。結(jié)合Kriging 代理模型和微種群遺傳算法，綜合考慮靜應(yīng)力和應(yīng)變能密度指數(shù)為優(yōu)化目標函數(shù)，建立了航空發(fā)動機風扇“參數(shù)化建?！邢拊抡妗獜姸群驼駝犹匦詢?yōu)化”一體化平臺。

（2）基于該平臺對某型寬弦風扇葉片改型后的直葉片進行了優(yōu)化，得到優(yōu)化設(shè)計的各項目標參數(shù)均有改善。優(yōu)化后的葉片靜應(yīng)力在葉片中心高應(yīng)力區(qū)向葉尖偏移，葉根前緣的應(yīng)力集中得到改善，且最大靜應(yīng)力比原型葉片的減小了5.45%，第1 階應(yīng)變能密度指數(shù)減小了5.94%，相比于單目標優(yōu)化結(jié)果更加全面有效。

（3）對比優(yōu)化前、后葉片的振動特性發(fā)現(xiàn)，葉片的振型、固有頻率和weaklink 點沒有明顯變化，第1 階振動應(yīng)力裕度從66.81%提高到了70.46%。說明優(yōu)化流程有效地改善了葉片的強度和高周疲勞特性，且對其他振動特性沒有負面影響。