載荷敏度抑制下渦輪盤拓撲優化及結構演化

閆浩，吳曉明

廈門大學 航空航天學院，廈門 361102

航空航天飛行器結構承受慣性載荷，發動機渦輪工作時要承受極大的離心力，大型土木建筑中結構自重是不可忽略的重要載荷。自重和慣性力這類載荷的特點是結構與載荷相互耦合，結構的變化會引起載荷大小的變化，呈現“沒有結構則沒有載荷”的特征，稱為設計相關載荷。

針對設計相關載荷的結構拓撲優化方法得到了國內外學者的廣泛關注。Bruyneel和Duysinx采用基于梯度的移動漸近線法(Gradient-Based Method of Moving Asymptotes, GBMMA)求解，討論了固定載荷與自重的比例對拓撲結構的影響。高彤等采取對角二次逼近法(Method of Diagonal Quadratic Approximations，MDQA)進行優化求解，提出了一種可變參數的材料屬性有理近似模型(Rational Approximation of Material Properties，RAMP)，避免固體各向同性懲罰模型(Solid Isotropic Material with Penalization, SIMP)處理慣性載荷問題存在的“附屬效應”。張暉等研究了自重載荷作用下的連續體結構的拓撲優化問題，使用了移動漸近線(Method of Moving Asymptotes，MMA)方法求解，提出了RAMP材料插值模型和平均敏度過濾技術相結合的求解策略，以克服柔度隨單元密度變化的非單調行為。Yang等對漸進結構優化(Evolutionary Structural Optimization,ESO)和雙向漸進結構優化(Bi-directional Evolutionary Structural Optimization, BESO)拓撲優化方法中的固定載荷條件下的靈敏度分析和進化過程進行了改進，以適應設計相關載荷狀況。Huang和Xie提出了一種基于靈敏度數計算的BESO方法，算例驗證了該方法對含自重載荷結構拓撲優化的有效性。Lee和Martins對設計相關的壓力載荷作用下的拓撲優化問題，運用材料邊界識別方法進行載荷施加，并提出了一種求解載荷敏度的解析方法。易壘和文毅針對目標函數對設計變量的非單調性，提出在優化準則法的迭代式中對慣性載荷敏度進行修正。

自專利文獻[8]提出一種新型雙輻板渦輪盤以來，國內外眾多學者開展了針對航空發動機渦輪盤在設計相關載荷下的拓撲優化研究。董少靜等采用漸進結構優化算法對渦輪盤進行拓撲優化，得到了雙輻板渦輪盤的結構形式。張乘齊等先通過拓撲優化得到雙輻板渦輪盤基礎輪廓，再通過形狀優化對雙輻板渦輪盤形狀進行設計。Rindi等運用水平集方法 (Level Set Method,LSM)法對某型渦輪盤進行了優化；Xu等將導重法(Guide Weight, GW)引入連續體結構在慣性力作用下的拓撲優化，對一些典型的設計相關載荷拓撲優化算例進行了計算，并且與前人的研究結果進行了對比。陸山和趙磊基于ANSYS軟件，建立了雙輻板渦輪盤優化設計平臺，針對某高壓渦輪轉子，設計了2種雙輻板渦輪盤。

基于變密度結構拓撲優化模型，應用優化準則法推導慣性載荷下的求解迭代格式，針對設計相關載荷SIMP模型的材料附屬效應運用一種新的指數材料性能近似模型；提出一種新的密度過濾方法；針對設計相關載荷帶來的目標函數非單調問題提出一種載荷敏度抑制方法。同時討論了載荷敏度抑制因子與渦輪盤拓撲結構的關系。通過有限元建模，比較分析了所得到的單輻板和雙輻板渦輪盤結構。

1 結構拓撲優化模型

以結構柔度最小化為優化目標的變密度法結構拓撲優化數學模型如下：

(1)

式中：為設計域內有限元單元的相對密度向量；為單元總數；為整體結構柔度值；、與分別為結構整體剛度矩陣、位移向量與載荷向量；為單元體積；、分別為設計區域初始體積與體積分數；為設計變量的下限。

模型式(1)中目標函數對設計變量的靈敏度和基于K-T條件的設計變量迭代格式分別如式(2) 和式(3)所示：

(2)

(3)

2 材料插值模型與灰度抑制方法

2.1 材料插值模型

若變密度拓撲優化方法中插值模型函數為()，則彈性模量、自重或離心載荷等設計相關載荷與單元相對密度之間的關系為

(4)

式中：為實體材料的彈性模量；0為實體材料第單元的慣性力。因此，載荷與彈性模量之間的比值如式(5)所示：

(5)

在SIMP插值模型中()=(>1)，為SIMP插值模型中的材料懲罰因子,由式(5)可知，在低密度單元中，設計相關載荷與彈性模量的比值會非常大，如圖1所示。這種由于插值模型引起的弱材料單元無法承載自身慣性載荷現象，被稱為“材料附屬效應”。這種附屬效應會導致計算無法進行或者使優化結果不穩定，也可能會使優化的拓撲結構存在不合理的分支結構。

圖1 SIMP模型下載荷與彈性模量比值Fig.1 Ratio of load to elastic modulus in SIMP model

為了克服這種材料附屬效應，文獻[1]應用了一種改進的分段SIMP插值模型方法，使得在低密度單元兩者的比例不再趨于無限大。文獻[2]采用RAMP模型，其材料插值函數如式(6)所示:

(6)

為RAMP插值模型中的材料懲罰因子,載荷與彈性模量的比值如式(7)所示：

(7)

由式(7)繪制RAMP模型對應的載荷與彈性模量之間的比值示意圖，如圖2所示。由圖2可知在RAMP模型中，低密度單元上兩者比值始終在有限值范圍內。

圖2 RAMP模型下載荷與彈性模量比值Fig.2 Ratio of load to elastic modulus in RAMP model

文獻[15]提出了一種材料性能指數近似模型(Exponential Approximation of Material Properties，EAMP)，其材料插值函數如式(8)所示:

(8)

式中:為EAMP插值模型中的材料懲罰因子。對應的函數圖像如圖3所示。

圖3 EAMP材料插值模型函數Fig.3 EAMP material interpolation model function

其載荷與彈性模量的比值如式(9)所示:

(9)

根據式(9)繪制EAMP模型對應的載荷與彈性模量之間的比值示意圖，如圖4所示。

由式(9)及圖4可知，EAMP模型中，在低密度區間內，載荷與彈性模量的比值也始終保持在有限值范圍內，同時算例表明EAMP材料插值模型能夠有效提高優化迭代的收斂速度。

圖4 EAMP模型下載荷與彈性模量比值Fig.4 Ratio of load to elastic modulus in EAMP model

2.2 新的灰度抑制方法

結構邊界不清晰，灰度單元多是拓撲優化結果中普遍存在的問題，基于EAMP材料插值模型，提出了一種灰度抑制方法。基于此方法的新的迭代格式如式(10)所示：

(10)

為了驗證基于EAMP模型的灰度抑制方法的可行性，對經典懸臂梁結構分別進行不使用灰度抑制和使用EAMP灰度抑制進行拓撲優化。

如圖5所示，設計域為長、寬0.5的矩形平板結構，將設計域離散為80×40個有限單元，結構左側固定位移全約束，右側受到向下的固定載荷。

圖5 懸臂梁結構示意圖Fig.5 Schematic diagram of cantilever structure

為了比較度量拓撲優化結構中的單元灰度，運用文獻[16]中的表征結構灰度值如式(11)所示，越小則結構的整體單元灰度越小。

(11)

圖6為不使用灰度抑制和使用EAMP灰度抑制方法得到的拓撲優化結構，各自的值、優化步數、目標函數如表1所示。

由圖6可知，使用EAMP灰度抑制之后的拓撲結構明顯更加清晰，單元灰度更小。由表1可知，在目標函數接近的情況下，使用EAMP灰度抑制方法比未進行灰度抑制的結構灰度值降低了84.5%，迭代收斂的步數也降低了70%。

圖6 優化結構圖Fig.6 Optimization structure chart

表1 灰度抑制結果Table 1 Gray suppression results

算例表明基于EAMP材料插值模型的灰度抑制方法的有效性，此方法可以有效地減少拓撲結構中的灰度單元并且可以更快達到收斂。

3 載荷敏度抑制方法

由式(2)可知，當載荷為設計相關載荷時，目標函數對單元密度的敏度式不再恒為負，即目標函數不再單調。為此提出一種載荷敏度抑制算法，在式(2)載荷敏度項前添加一個抑制因子1/(+1)。添加抑制因子之后，目標函數對設計變量的靈敏度為

(12)

(13)

抑制因子作用于優化過程的每一次迭代中，其作用如圖7所示，其中實線為剛度敏度，點劃線為抑制前的載荷敏度，虛線為抑制后的載荷敏度。經過載荷敏度抑制使目標函數保持單調。

圖7 載荷敏度及其抑制算法示意圖(局部放大)Fig.7 Schematic diagram of load sensitivity and its suppression algorithm (Local amplification)

對于抑制后的敏度式(12)，采用Sigmund敏度過濾方法，得到迭代格式(10)中的迭代系數為

(14)

優化求解的流程圖如圖8所示。

圖8 優化求解流程圖Fig.8 Flow chart for optimization solution

為了驗證此方法的有效性，對照文獻[3，5，7]中的自重載荷算例，在相同計算條件下使用所提出的載荷敏度抑制算法進行計算，對拓撲優化結果進行比較。表2為文獻中算例的計算參數設置。

表2 文獻算例參數Table 2 Reference example parameters

表3為文獻與載荷敏度抑制算法所得到的拓撲結構。由表3可知，在計算參數相同的情況下，載荷敏度抑制算法與文獻算例的結構相似。文獻[5]給出了算例的目標函數值3.82和收斂步數80，本文載荷敏度抑制算法得到的目標函數值3.79，收斂迭代步數30，目標函數值略低且收斂速度更快，說明載荷敏度抑制算法在自重載荷拓撲優化中的有效性。

表3 不同算法得到的拓撲優化結構Table 3 Topology optimization structure obtained by different algorithms

4 發動機渦輪盤拓撲優化

航空發動機渦輪盤工作時要承受極大的離心載荷，是一種典型的設計相關載荷的結構優化問題。在美國 “綜合高性能渦輪發動機技術”(Integrated High Performance Turbine Technology，IHPTET)計劃的第3階段中，驗證了雙輻板渦輪盤在結構傳力路徑等方面的優勢，且其較單輻板渦輪盤質量減輕、轉速提高。1999年Cairo申請了雙輻板渦輪盤結構的專利，2001年Burge發明了應用于壓氣機的雙輻板輪盤結構，2007年Harding和Curtiss發明一種盤轂尺寸較大的雙輻板渦輪盤結構。

相關研究已經表明，在相同的載荷和設計邊界條件下，雙輻板渦輪盤較單輻板渦輪盤具有降低最大應力、變形小等優勢。但是從拓撲優化算法角度，單輻板與雙輻板結構之間如何演化，以及這樣的演化受到哪些因素的影響未見文獻報道。以下根據本文提出的載荷敏度抑制算法進行討論。

4.1 抑制因子對渦輪盤拓撲結構的影響

抑制因子1/(+1)中的由式(13)定義，與單元的載荷敏度和剛度敏度相關，由式(12)和式(14) 可以看出，抑制因子的大小決定了載荷敏度在求解迭代格式(10)中的影響。有學者為避免設計相關載荷目標函數的非單調問題，在優化迭代計算時直接忽略載荷敏度項。為討論載荷敏度對最終拓撲結構的影響，將式(12)中的抑制因子設置為1/(+)，稱為抑制系數，當從小到大變化時，載荷敏度在迭代格式(10)中的影響逐步變小，討論此時拓撲結構的演化。

將渦輪盤的徑向截面設置為一個高寬為×2的矩形平面，離散為50×100個四邊形有限單元。矩形平面繞左側距離為的軸旋轉，轉速為，如圖9所示。

圖9 結構示意圖Fig.9 Structure diagram

運用本文所提出的算法，其中EAMP插值模型和灰度抑制的懲罰因子分別為、，敏度過濾半徑為，體積分數為，計算參數的具體數值如表4所示。

表4 算例3參數Table 4 Parameter of Example 3

在抑制系數取0.1、1、10、100、1 000、10 000、100 000和無窮大(即不考慮迭代式中的載荷敏度項)的情況下，對應的渦輪盤拓撲結構演化規律如圖10所示，其中橫坐標為抑制系數的常用對數。

由圖10可見，渦輪盤徑向截面的拓撲結構隨抑制系數的增大(載荷敏度在迭代式中的影響變小)，由雙輻板結構逐漸演化為單輻板結構。抑制系數取1 000和10 000之間存在一個臨界點，在此臨界點附近拓撲結構從雙輻板向單輻板突變。應用二分法計算的臨界點為2 301和2 302，如圖11所示。

圖10 抑制系數變化時拓撲結構演化圖Fig.10 Topological structure evolution diagram with variation of suppression coefficient

圖11 臨界點附近的拓撲結構演化圖Fig.11 Topological structure evolution graph near critical point

當抑制系數取2 301時，結構為雙輻板結構，取2 302時，結構為單輻板結構。為了進一步探究抑制系數對結構演化的作用，分別計算抑制系數取2 301和2 302時，優化迭代過程中拓撲結構的演變過程。

從圖12可以看到，2個抑制系數下拓撲結構在迭代65步之前得到的結構是相似的，在65～72步之間結構發生不同演變，=2 301時結構的上下輻板得到加強，逐步演化為雙輻板結構。=2 302時結構中間的輻板得到加強，逐步演化為單輻板結構。

圖12 迭代過程中拓撲結構的演變圖Fig.12 Evolution graph of topological structure in iterative process

通過典型大慣性載荷算例3中抑制系數大小對拓撲結構演化的討論可以看出，一方面設計相關載荷的載荷敏度對最終的拓撲優化結構有關鍵性作用。當抑制系數較小，即載荷敏度在迭代式中保留較大時，渦輪盤徑向截面優化拓撲為雙輻板結構。抑制系數較大，即載荷敏度在迭代式中保留較小時，拓撲結構為單輻板結構。另一方面，從圖10可以看出，渦輪盤從雙輻板結構向單幅板結構演化過程中，目標函數逐步增加，即結構柔度增大，剛度減小。但圖11顯示，在抑制系數的臨界點附近，雙輻板結構可能與單輻板結構剛度相當甚至更差。

將算例3得到的2種渦輪盤徑向截面拓撲結構設計成渦輪盤模型進行有限元分析，對其強度和剛度進行比較分析。

4.2 幾何建模與有限元分析

算例3中當取值非常小時渦輪盤的拓撲結構為典型的雙輻板結構，當取值無窮大時渦輪盤的拓撲結構為單輻板結構，見圖13。對2種拓撲結構分別建立渦輪盤的幾何模型，如圖14所示，對2個幾何模型分別施加相同的約束條件，進行有限元分析，得到應力、應變能結果以及模型的質量如表5所示。

圖13 渦輪盤徑向截面拓撲結構Fig.13 Topological structure of radial section of turbine disk

圖14 渦輪盤結構幾何模型Fig.14 Geometric model of turbine disk structure

2種渦輪盤結構在相同約束和載荷條件下的應力和應變能云圖如圖15和圖16所示。由表5可知，在徑向截面面積相同的情況下，雙輻板渦輪盤整體的質量比單輻板減輕16%。圖15顯示了單、雙輻板的應力分布狀況，最大應力均在輪心位置，這與文獻[24]得到的結果相同。相比于單輻板，雙輻板的最大應力下降了23%。由圖16 可以看出，雙輻板的應變能分布均勻，整體應變能低于單輻板。

表5 有限元分析結果Table 5 Results of finite element analysis

圖15 渦輪盤應力云圖Fig.15 Stress distribution of turbine disk

圖16 渦輪盤應變能云圖Fig.16 Strain energy distribution of turbine disk

5 結 論

1) 針對設計相關載荷拓撲優化中存在的“材料附屬效應”問題，運用一種指數型EAMP材料插值模型，并且基于此模型構建了一種新的灰度抑制方法。算例表明，該方法具有單元灰度少、結構清晰、收斂速度快的特點。

2) 對設計相關載荷變密度拓撲優化中存在的目標函數非單調問題，提出一種載荷敏度抑制方法。算例表明了該方法對自重和大慣性載荷下結構拓撲優化的有效性。

3) 應用本文提出的算法對發動機渦輪盤的徑向截面進行拓撲優化。通過抑制系數的討論，說明了設計相關載荷的載荷敏度對渦輪盤的徑向截面拓撲結構的影響，揭示了傳統單輻板渦輪盤結構與近年廣泛研究的雙輻板渦輪盤結構之間的演化規律和內在聯系。在相同約束和離心載荷條件下對兩種結構渦輪盤進行了有限元分析。