基于BESO 算法的渦輪盤拓撲優化1)

孟子皓 * 丁 超 * 任毅如 * 米 棟 錢正明 張立章

*(湖南大學機械與運載工程學院，長沙 410082)

?(中國航發湖南動力機械研究所，湖南株洲 412002)

渦輪盤作為航空發動機的核心部件之一，長期在高溫高壓的惡劣條件下運轉，還需要承受高轉速產生的離心載荷，其結構和性能對于航空發動機的可靠性和性能有決定性影響，對渦輪盤進行結構優化和輕量化設計能有效地提高航空發動機的推質比和燃油效率。

拓撲優化是在給定負載、約束和性能指標的條件下，在設計區域內尋找材料最優分布的結構優化方法，被廣泛應用于航空航天等領域[1-3]。Xie等[4]首次提出的漸進結構優化方法（evolutionary structural optimization，ESO 方法），已廣泛運用于各個工程領域中[5-7]。Querin 等[8]基于ESO 算法，提出了雙向漸進結構優化策略（bidirectional evolutionary structural optimization， BESO），這種改進的算法解決了ESO 算法中被刪除的單元無法恢復的缺陷，實現了優化過程中對設計域內低效單元的刪減的同時又增加了高效單元。在此基礎上，為了解決結構慣性載荷作用下結構拓撲優化設計問題，Yang 等[9]利用BESO 算法研究了重力載荷作用下平面矩形板的結構拓撲優化，并得到了“拱橋”結構的優化結果。為了提高結果的收斂效率， Huang 等[10]在BESO 算法的基礎上，引入材料屬性有理近似模型(rational approximation of material properties，RAMP)，解決了結構自重載荷作用下結構拓撲優化設計問題。

國內外對渦輪盤的結構優化設計進行了一定的研究，Bhavikatti 等[11]基于旋轉圓盤的研究，使用五次多項式定義了圓盤橫截面的形狀，并通過有限元方法對圓盤進行了應力分析，采用非線性規劃的方法，即改進的順序線性規劃移動劃分優化技術，研究了圓盤的形狀優化設計。Rindi 等[12]提出了使用水平集方法（level set method，LSM）來解決由于引入偽材料密度而引起的問題，并將其優化方法應用到渦輪盤構件。Shen 等[13]利用ESO 方法優化結構拓撲，提出了雙輻板渦輪盤形狀，并通過旋轉三維光彈性試驗進行了驗證。閆浩等[14]提出了一種載荷敏度抑制算法，解決了設計相關載荷的靈敏度計算問題和變密度拓撲優化中存在的目標函數非單調問題，并從優化算法的角度揭示了單輻板與雙輻板渦輪盤結構之間的聯系及演化規律。

渦輪盤結構的傳統設計方法大部分缺乏理論基礎，主要依靠大量的設計經驗和反復的實驗基礎得到，研發周期長、成本高。拓撲優化設計有設計變量多樣性、優化構型可靠性、研制成本低以及設計效率高等優勢，引入拓撲優化技術進行渦輪盤結構設計，可以有效避免傳統設計方法的缺陷。國內外目前對渦輪盤拓撲優化設計研究不充分，相較而言，當前針對渦輪盤的研究所呈現的迭代步數多，優化效率較低，優化過程和優化結果無法滿足飛行器結構設計要求的高速發展。本文主要圍繞航空發動機渦輪盤結構的拓撲優化設計，針對渦輪盤的工作狀態和受載情況，將 BESO算法應用于受離心載荷作用下渦輪盤的結構優化，提高了設計效率。采用ANSYS 參數化設計語言（ANSYS parametric design language，APDL）進行二次開發，無需在有限元分析軟件和計算軟件之間反復輸入和輸出數據，在有限元軟件上實現了渦輪盤的BESO 拓撲優化設計，開發的拓撲優化APDL 代碼具有拓展到其他大型設計相關載荷的工程問題的前景。

1 BESO 算法

1.1 結構拓撲優化模型

以結構應變能最小化為優化目標的BESO 法的優化問題數學模型可以表達為

1.2 材料插值方案

常用的材料插補方案有固體各向同性懲罰模型(solid isotropic material with penalization，SIMP)和RAMP模型等，分別如式（2）和式（3）所示。Huang 等[15]將SIMP 模型應用于求解離散拓撲優化問題的BESO 方法。后來，有研究者研究了RAMP 模型在求解慣性載荷拓撲優化問題中的應用，并說明RAMP 的優勢[16-18]更適用于求解慣性載荷拓撲優化問題。本文采用的是RAMP 模型。

式中，E(xi) 為第i個單元的楊氏模量，E0為初始設置的楊氏模量，p和q為懲罰因子。

1.3 靈敏度計算

當載荷為離心載荷、重力載荷等設計相關載荷時，單元增加或刪除而引起的結構柔順性的變化可以表示為

定義單元靈敏度αi為

式中，ui為單元i的節點位移矢量，Ki0為初始單元i的剛度矩陣，ω為轉速，ri為單元i到旋轉軸的距離，vi為單元i的體積，ρ0為單元初始密度，fˉ 為離心載荷的方向矩陣。為了保證結構收斂，應該選擇較大的懲罰因子，但如果懲罰因子無限大，當x(i)=xmin時，靈敏度為0，這就等同于“硬殺傷”BESO 方法。

在漸進結構優化設計中，優化結果往往會出現類似于圖1(a)所示的棋盤格現象，優化結果圖中孔洞材料單元與實體材料單元交錯分布的“棋盤格現象”。為消除棋盤格現象導致的產品制造困難的缺陷，采用基于獨立網格濾波技術的敏度過濾方法[15]對結構優化過程中單元敏度進行過濾處理。計算原始單元靈敏度αi平均分配得到節點靈敏度αj；通過比較第i個單元中心點位置到節點的距離rij和過濾半徑rmin之間的大小得到加權函數ω(rij) 的值。之后將節點靈敏度αj與權函數ω(rij) 代入式（6），得到新的單元敏度值αi′，該過濾技術會使得空白單元有機會轉化為實體單元。

圖1 基于敏度過濾的棋盤格抑制結構Fig.1 Checkerboard suppression structure based on sensitivity filtering

使用上述敏度過濾方法后優化結果如圖1(b)所示，對比圖1(a)，該拓撲優化結構中的棋盤格現象基本已經消除。

1.4 離心載荷下算例計算

算例為如圖2 所示的渦輪盤初始結構模型，其幾何模型的長（L）和寬（H）分別為100 mm和50 mm，厚度為1 mm，彈性模量E=200 GPa，泊松比ν=0.3，密度ρ=7.8 g/cm3?，F將懸臂梁整體繞旋轉軸Y作勻速轉動，轉速ω0=1000 rad/s，r0=20 mm。其中黑色部分為非設計優化區域，目標體積分數為40%?，F以結構整體剛度最大為優化目標，結構的總體積作為體積約束進行離心載荷作用下的拓撲優化設計。

圖2 渦輪盤初始結構模型Fig.2 Initial structure model of turbine disk

圖3 給出了拓撲結果的演化圖，通過離心載荷下渦輪盤拓撲優化設計的文獻資料可以看出，渦輪盤拓撲結構的演化進程是逐步向孔洞狀的雙輻板結構進化。經過34 次迭代即可獲得收斂結構。結構總柔度下降，滿足總體剛度提升的優化目標。

圖3 拓撲結果演化圖Fig.3 Evolution diagram of topological results

文獻結果如圖4 所示。對比該算例的結果，二者結果相似度較高，說明該算法是可行的，相較而言，該算法僅需34 步迭代即可獲得收斂的優化結果，優化效率明顯提高。

圖4 渦輪盤離心載荷下的拓撲優化結果[14]Fig.4 Topology optimization results of turbine disk under centrifugal load[14]

2 渦輪盤拓撲優化設計

2.1 渦輪盤幾何模型和有限元模型的構建

某型發動機渦輪盤的3D 幾何模型如圖5 所示，若直接進行網格劃分會出現網格質量差的問題，從而影響計算效率和求解精度?，F對渦輪盤幾何模型進行簡化：（1）除去非設計優化域的葉片結構，在渦輪盤邊緣施加等效節點力近似替代葉片對渦輪盤的作用，假設葉片載荷均勻分布于渦輪盤外緣表面[19]，可通過式（8）進行計算；（2）渦輪盤為旋轉對稱結構，故選取截面進行優化可以提升優化效率。簡化渦輪盤模型用APDL創建如圖5 所示。

圖5 簡化的渦輪盤Fig.5 Simplified turbine disk

式中，σrb為均勻分布于渦輪盤外緣表面的葉片載荷，F為葉片和榫頭凸塊自身的離心載荷，該載荷為固定載荷，在拓撲過程中保持不變，R為榫槽底半徑，H為渦輪盤輪緣的軸向寬度。

選取簡化后渦輪盤的截面進行優化設計區域的劃分如圖6 所示，其中黑色部分為非設計區域，網格單元類型選用平面4 節點的PLANE42 的單元，按照平面結構四邊形網格自動劃分。

圖6 優化設計區域的劃分Fig.6 Division of optimal design domain

僅考慮盤身和葉片的離心載荷，對渦輪盤的內徑接觸面施加固定約束，對整體結構施加繞X軸的轉速ω0。該模型在簡化過程中刪除了葉片結構，通過計算得到葉片的離心力以等效節點載荷的形式作用在輪盤外徑上，如圖7 所示。

圖7 渦輪盤的離心載荷、葉片等效載荷和邊界約束條件Fig.7 Centrifugal load, blade equivalent load and boundary constraint condition of turbine disk

2.2 渦輪盤在離心載荷下的拓撲優化結果分析

在完成優化前處理操作步驟后，接下來針對優化主循環的迭代結果進行分析，如圖8 所示。迭代過程如圖9 所示。

圖8 渦輪盤在離心載荷下拓撲優化的演化結果Fig.8 Evolution results of topology optimization of turbine disk under centrifugal load

圖9 拓撲結果演化圖Fig.9 Evolution diagram of topological results

通過優化結果的應力分布圖和應變能分布圖可以看出，優化后的結構內部單元von Mises 應力和應變能分布更均勻，結構總柔度減小，剛度有所提高。其次，通過對比優化前后結構最大應力值與最大單元應變能云圖可以發現，優化后的渦輪盤結構的最大應變能密度較優化前要低48%，且應變能分布更均勻。最大von Mises 等效應力比優化前降低了25%。最后，從優化目標來看，渦輪盤結構在減重26%的條件下，結構總應變能降低了48%，滿足提高結構整體剛度的目的。

3 結論

通過ANSYS APDL 有限元軟件建立了離心載荷下BESO 拓撲優化算法模型，對渦輪盤進行拓撲優化，得到以下結論。使用BESO 算法對設計相關載荷的優化問題求解，優化過程中迭代步數較少，同時避免了在計算軟件和有限元軟件之間切換，提升了拓撲優化效率。將渦輪盤模型進行簡化，避免了因葉片結構復雜導致參數化建模困難和效率低下的問題，提升了優化設計效率，同時通過施加葉片等效載荷的方法，考慮葉片在離心載荷下對渦輪盤的影響。在相同的外載荷條件下，優化后的渦輪盤比傳統渦輪盤的質量降低了26%，最大von Mises 等效應力比優化前降低了25%，總應變能降低了48%，滿足提高結構整體剛度的目的，滿足了某型號發動機的輕量化設計對渦輪盤的要求。進一步驗證了雙輻板渦輪盤更適合于高推質比的發動機。