輪盤概念設計中拓撲和形狀同時優化方法

范俊，尹澤勇，王建軍，米棟，3，閆成

(1.北京航空航天大學 能源與動力工程學院，北京100191;2.陸航研究所，北京101121;3.中航工業 航空動力機械研究所，株洲412002)

結構優化包括拓撲、形狀和尺寸優化3個階段，對于前2個優化階段，通常的優化順序是，先進行拓撲優化，在拓撲優化結果的基礎上，再進行形狀優化.在這種分步進行的優化過程中，需要人為設定一個大小一定的不可調的初始設計區域用于拓撲優化;并且在基于變密度法(即SIMP法)的拓撲優化過程結束后，還需要人為設定一個用于刪除中間密度單元的相對密度門檻值.然而，假如這兩種人為設定參數選擇不恰當，就有可能使得優化結果并不是最優解.

為了減小上述兩類人為選擇參數對優化結果的影響，本文提出了一種拓撲和形狀同時優化(STSO)的方法.國外學者 Deaton等［1］認為從拓撲優化結果到形狀優化初始設計模型的轉化，是未來連續體多學科拓撲優化研究重點之一.Ansola等［2］論述了一個形狀和拓撲組合優化的方法，在該方法中，交替進行形狀和拓撲優化步驟，而在每一個步驟中，先進行形狀優化再進行拓撲優化;這些步驟不斷重復，直到得到一個收斂的全局解.Hassani等［3］論述了一個板殼結構形狀和拓撲同時優化的方法，在該方法中，形狀和拓撲優化不是分開進行的，而是并行的.上述兩篇文獻除了優化流程不一致之外，在優化算法上也是不一樣的，Ansola采取的是單純形法，Hassani采取的是移動漸進線法.國內學者張衛紅等［4］為了解決壓力載荷作用下的結構輕量化設計問題，直接采用B樣條曲線描述壓力加載面，通過拓撲和形狀變量的聯合優化滿足了工程實際對結構輕量化與邊界的功能性與光滑性設計要求.雖然本文中的拓撲和形狀優化也是并行的，但與上述文獻只考慮應變能和應力優化響應不同，本文的優化模型中還考慮了以頻率為約束的三維實體結構(輪盤)動力學優化，采用的優化算法是序列二次規劃優化算法(SQP法).

輪盤是航空發動機中關鍵部件，已有文獻對輪盤進行了拓撲優化設計研究.文獻［5］以體積剛度比為目標，應力、破裂轉速、低循環疲勞壽命為約束建立拓撲優化數學模型，采用基于隨機抽樣敏度分析的雙向漸進結構拓撲優化方法，得到了航空發動機雙輻板形式的渦輪盤結構.文獻［6］以體積剛度比為目標、應力為約束建立了拓撲優化數學模型，采用基于自適應隨機抽樣策略的周期結構拓撲優化方法得到了多輻板風扇盤.文獻［7］以最小體積為目標、應力為約束建立了典型風扇盤子午面結構優化設計模型，得到了一個多輻板盤結構.從國外先進的發動機來看，許多高性能發動機寬弦風扇盤經過拓撲優化設計，選用了多輻板的盤轂混合式輪盤結構［6］，如CFM56-7，RB211-535E4等采用了雙輻板輪盤結構.但是，上述文獻一般只考慮了靜力學目標和約束，且拓撲和形狀優化都是分步進行的.為了設計出既輕且剛性適當又安全可靠的輪盤，有必要研究輪盤的模態特性，進行與實體輪盤頻率等有關的動力學優化設計.

本文首先建立了基于SIMP法的拓撲和形狀同時優化(STSO)數學模型，分析了相應的靈敏度，使用SQP法進行求解.然后，以板殼結構為例，進行了拓撲和形狀同時優化，將其結果與分步優化進行比對，探討了相對密度門檻值的選取和兩種優化變量之間的相互作用對最終優化結果的影響，說明拓撲和形狀同時優化的優點.最后，本文將同時優化的方法應用于輪盤概念設計中，對比分析了不同振型所對應的頻率約束時實體輪盤拓撲和形狀同時優化方法與單獨拓撲優化方法的不同結果.

1 STSO方法數學模型

SIMP法是拓撲優化方法的一種，由于其數學論證嚴密、便于實施等優點而得到廣泛應用［8-9］.本文的STSO方法將基于SIMP法來展開研究.

SIMP法建立在結構有限元模型的基礎上，設計變量是各單元的相對密度，是在［0，1］區間的連續變量，第 i個單元的相對密度用 ρi來表示(i=1，2，…，N)，N指單元總數.設單元的原始密度是 ρ，優化后密度是 ρe，則存在關系式 ρi= ρe/ρ;SIMP法中，單元材料屬性是單元相對密度的指數函數，設第i個單元初始和優化后的彈性模量分別是E0，Ei(ρi)，第i個單元初始和優化后的剛度矩陣分別是K0i，Ki(ρi)，則存在關系式:

式中P為懲罰因子，其作用是對中間相對密度進行懲罰，使結構單元相對密度盡可能趨于0或者1.用于航空發動機輪盤拓撲和形狀優化的數學模型是以體積、應力和頻率為約束，以最小化結構應變能為目標.

式中，C為應變能;F為載荷;U為位移;σeq為當量應力;σ0.1為 0.1% 屈服強度;V0i為第i個單元的初始體積;V為優化之后結構體積;V0為結構原始體積;V0iρi為優化后單元的等效體積;g(f)為頻率f的函數;K為全局剛度陣;M為全局質量陣;λ為特征值;φ為對應的特征向量;Di為第i個節點的位移矢量;D為一個位移常量;j和為形狀優化變量yj的上下限.

2 靈敏度分析

靈敏度分析包括拓撲優化變量和形狀優化變量靈敏度分析兩個部分，本文著重分析目標和約束對拓撲優化變量的靈敏度.

2.1 頻率靈敏度

頻率約束的靈敏度，是通過多自由度無阻尼系統自由振動的運動方程，利用特征值法求得.對式(2)中的運動方程，通過對特征向量用質量矩陣進行規格化［10-11］，則特征值對設計變量的導數可表示為

式中，λj為第j階特征值;φj為第j階特征值對應的特征向量.由于

式中，M0i為第i個單元的初始質量矩陣;Mi為第i個單元優化后的質量矩陣，所以

將式(5)代入式(3)，可得

由 λj=f2j，可得

2.2 應變能、體積和位移靈敏度

由于式(2)中外力F是一個常量，則有

由于結構剛度矩陣的對稱性，并將式(8)代入可得應變能的靈敏度為

容易得到，體積約束的靈敏度為

由式(8)可得位移靈敏度為

2.3 應力靈敏度

現階段應力相關的優化問題面臨著諸多的難題［12］，如優化過程中由于應力約束不連續導致的應力奇異現象、局部應力約束過多導致的超大計算規模等計算難題.隋允康等［13］等利用Mises強度理論，提出了應力約束全局化策略，將局部的應力約束問題轉化為結構整體的應變能約束問題.榮見華等［14］在優化迭代循環的每一輪子循環迭代求解開始時，通過形成和引進新的位移和應力約束限，自動構建設計變量移動限，將結構應力約束歸并為幾個最可能的有效應力約束，從而減少應力靈敏度的分析量.París等［15］總結了連續體拓撲優化中應力約束靈敏度分析的一般方法.

式(2)中第i號單元的VonMises應力近似地可以表示為

式中，Bhi為第i號單元的第h個高斯積分點的幾何矩陣;ni為第i號單元的高斯積分點數;Ui為僅與第i號單元所有自由度相應的U中元素構成的單元位移矢量.則應力靈敏度為

由式(11)可以求得?Ui/?ρi為?U/?ρi中與第i號單元所有自由度相應元素構成的矢量，將其代入式(13)，則可以求得應力靈敏度.

對于形狀優化變量的靈敏度，并不能十分明確地列出［3］.在本文中，采用一種半解析的辦法，這樣，形狀優化變量的靈敏度一部分由解析解推出來，一部分由有限差分近似推出來，具體的推導見文獻［16-18］.

在求得目標和約束函數的靈敏度基礎上，運用SQP法進行求解，SQP算法具體原理可參考文獻［19-20］.

3 分步優化和同時優化對比分析

以圖1中的板殼結構為例來進行對比分析，待優化的板長度為100mm，寬度為10mm，厚度為1 mm，劃分成1000個1×1的單元.

板的左端受全約束，板的右端上頂點處施加100 N的集中力.材料的彈性模量是 2.1×105MPa，泊松比是0.3.優化的設計區域是整個黃色區域.為了對比分步優化和同時優化結果，相比于式(2)，減少了頻率和應力約束，將優化目標改為體積最小，約束是板的右端下頂點的位移小于3 mm，建立數學模型見式(14).

式中d101為第101個節點的垂直位移，該節點位于板的右端下頂點.

圖1 結構優化的初始有限元模型Fig.1 Original finite element model for structural optimization

3.1 先拓撲后形狀分步優化

圖2(a)描述了通常一般優化順序:先進行拓撲優化，在拓撲優化結果的基礎上，再進行形狀優化.其中，在基于SIMP法的拓撲優化過程中，并不刪除單元，而是通過改變單元的相對密度值來滿足約束條件.在拓撲優化過程結束后，需要人為設定相對密度門檻值，刪除低于該值的單元后，得到最終拓撲優化結果，用于形狀優化.

圖2 結構優化流程Fig.2 Flow of structure optimization

按照這個優化順序，先對圖1中的板殼結構進行拓撲優化，迭代25步收斂，取相對密度門檻值為0.16，刪除相對密度小于該值的單元，得到的拓撲優化結果如圖3(a)所示.然后，定義該結果中的形狀優化變量.本文出于演示論證的目的，只定義了一個形狀優化變量.設形狀變量為圖形上邊界中點處節點的位置，該節點可上下移動.同樣采用SQP算法，迭代3步后收斂，得到的結果如圖3(b)所示.從結果看出，在拓撲結構一定的情況下，進行形狀優化，其拓撲構型并沒有改變，只是形狀的邊界尺寸發生了變化，優化的效果不明顯.

圖3 拓撲優化結果及基于拓撲優化的形狀優化結果Fig.3 Topology optimization result and shape optimization result based on topology optimization

因此，拓撲優化結果對于后續的優化有著重要的作用，而人為選定的相對密度門檻值對拓撲優化結果有著重要的影響.為了減小這種主觀選擇參數對結果的影響，有必要拓撲和形狀同時優化.

3.2 先形狀后拓撲分步優化

不同的初始設計模型會產生不同的拓撲優化結果，為了得到較好的拓撲優化結果，文獻［2］認為，應該先通過形狀優化來確定初始設計模型，再進行拓撲優化.本文將以算例表明，先形狀后拓撲的優化流程也并不一定能得到最優解.

對于圖1所示模型，先進行形狀優化，以確定用于拓撲優化的初始設計模型.定義兩個形狀優化變量:一是模型上邊界右端頂點處節點的位置，該節點可以上下移動;二是模型上邊界中間點處節點的位置，該節點可以上下移動.優化在迭代5步后收斂，得到的形狀優化結果如圖4(a)所示.將該形狀優化結果作為初始設計模型，再進行拓撲優化，迭代5步后收斂，得到的拓撲優化結果如圖4(b)所示.從圖中可以看出，在形狀優化結果基礎上再拓撲優化，發生的改變并不多.這也表明，經過形狀優化確定的初始設計模型，不一定是最合適的用于拓撲優化的初始設計模型;同時從圖3(a)和圖4(b)的比對看出，不同的初始設計模型進行拓撲優化會產生不同的優化解.

圖4 形狀優化結果及基于形狀優化的拓撲優化結果Fig.4 Shape optimization result and topology optimization result based on shape optimization

3.3 拓撲和形狀同時優化

在分步優化過程中，人為設定的初始設計區域和相對密度門檻值，對最終優化結果有著重要的影響，為了減小這種人為設定參數的影響，本節使用拓撲和形狀同時優化的方法對該板殼結構進行優化.

對圖1中的初始設計模型，進行拓撲和形狀同時優化，其中，形狀優化變量與3.2節保持一致.優化的流程如圖2(b)所示，迭代59步后得到STSO優化結果如圖5(a)所示，優化結果類似于桁架結構.該結果的重量已經較小，可直接進入詳細設計階段對應的尺寸優化.

圖5 形狀和拓撲同時優化結果Fig.5 Simultaneous shape and topology optimization result

因為3.1節和3.2節所述分步優化過程中，第2步優化的效果相比于第1步要小得多，圖5(b)中描述了兩種分步優化方法的第1步優化結果(體積)和同時優化結果隨優化步數變化趨勢.從圖中可以看出，相比于分步優化，同時優化后的結構體積要小得多，這說明同時優化方法的效果是最好的.

此外，從上述的優化結果對比分析看，形狀和拓撲優化是一個整體，是相互作用的兩個過程，割裂地進行優化有可能得不到全局最優解.

4 STSO方法在輪盤概念設計階段的應用

4.1 初始設計模型

使用STSO方法優化的輪盤是一個全三維模型，該模型所參考的航空發動機渦輪盤子午面如圖6(a)所示［4］.其中，中心孔半徑 Rb=65 mm，輪緣半徑Rr=278.6 mm，輪緣寬度Wr=44 mm，輪緣高度Hr=11.4mm，輪轂寬度Wb=91mm;輪緣均布載荷Pr=169.39 MPa，工作轉速 N=13333 r/min;輪盤所用的材料為 GH4169，文中計算采用了600°C時該材料的性能參數.

圖6(b)是用于結構優化的初始設計模型，其中藍色非設計區域與參考盤的輪緣部分是一致的;綠色設計區域寬度參考輪轂寬度，將其暫定為Wb=160 mm，綠色設計區域高度為Rr－Rb－Hr=202.2 mm.

STSO方法在輪盤應用中的優化數學模型除了不考慮位移約束外，與式(2)所示相同.需要說明的是，本文目的主要是闡述STSO方法和概念，如圖2(c)所示，該方法是應用于輪盤概念設計階段，其優化結果距實際應用還有較長距離.因此，雖然輪盤不同的地方應有不同的應力約束條件，但是為了突出動力學優化以及對比STSO方法和單獨拓撲優化方法的優劣，本文將各個地方的靜力學要求簡化為0.1%屈服極限的6/10.

另外，在STSO方法中，形狀優化設計變量為圖6(b)中4個角點(1，2，3，4 點)水平方向的位置，為了獲取最佳的初始設計區域，將其變化范圍設為［－80，80］;拓撲優化設計變量為綠色區域單元的相對密度值.

4.2 STSO方法在輪盤概念設計階段的應用

圖6中初始設計模型的第4階和第5階頻率是679 Hz和1060 Hz，其振型如圖7所示，分別是節圓和節徑振動.從第3節可知，分步優化中第2步優化的效果相比于第1步要小得多，并且同時優化方法所得的優化結果如有必要也可以用形狀和尺寸優化方法進行進一步的細節設計，因此，本節將運用STSO方法和單獨拓撲優化方法對該輪盤進行結構優化，對比分析兩種不同方法下不同振型對應的頻率約束對輪盤優化結果結構形式的影響.

圖6 初始設計模型Fig.6 Initial design model

圖7 初始設計模型的第4和第5階振型Fig.7 Fourth and fifth mode of vibration of initial design model

4.2.1 降低節徑振動頻率

將式(2)中的頻率約束g(f)≤0確定為降低節徑振動頻率約束為f5≤900.單獨進行拓撲優化，得到的優化結果見圖8(a)所示;拓撲和形狀同時優化的結果如圖8(b)所示.

圖8 降低節徑振動頻率的優化結果對比Fig.8 Comparison of optimization results with lower frequency of nodal diameter vibration

從結果看，兩種方法得到優化結果拓撲形式是不一樣的，同時優化結果拓撲形式較單獨拓撲優化結果更為簡潔.

圖9分別列出了拓撲和形狀同時優化，以及單獨拓撲優化的應變能、頻率、體積分數(優化后的體積與初始體積的比值)隨迭代步數變化情況.從圖中可以看出，兩種優化的頻率和體積約束在最后迭代步均分別達到相應上限值，而優化的目標——應變能則不同，拓撲和形狀同時優化結果的應變能較單獨拓撲優化的應變能更小，減小了約2.4%，并且其收斂速度較快.

圖9 降低節徑振動頻率的優化結果隨迭代步數變化趨勢對比Fig.9 Variation tendency comparison of optimization results with lower frequency of nodal diameter vibration in different iteration steps

4.2.2 提高節徑振動頻率

明確式(2)中的頻率約束，提高節徑振動頻率約束為f5≥1100，優化結果如圖10所示.圖中給出的是優化后相對密度的分布情況(刪除了較小相對密度單元的優化結果)，紅色表示優化后相對密度較高的單元顏色，藍色表示優化后相對密度較低的單元顏色.

在圖10(a)和圖10(b)中，當使用單獨拓撲優化方法時，優化結果出現了開孔現象.而在圖10(c)和圖10(d)中，當使用STSO方法時，優化結果與參考盤類似，是一個單幅板的輪盤，并未出現開孔現象.

相較于單獨拓撲優化方法的結果，使用STSO方法得到的結果，更符合現役發動機少有開孔輪盤的現狀.

圖11也分別列出了拓撲和形狀同時優化，以及單獨拓撲優化的應變能、頻率、體積隨迭代步數變化情況.從圖中可以看出，兩種優化的體積約束在最后迭代步均達到相應上限值，而優化的目標——應變能則不同，拓撲和形狀同時優化結果的應變能較單獨拓撲優化的應變能更小，減小了約8.5%，并且其收斂速度較快.

圖10 提高節徑振動頻率的優化結果對比Fig.10 Comparison of optimization results with higher frequency of nodal diameter vibration

圖11 提高節徑振動頻率的優化結果隨迭代步數變化趨勢對比Fig.11 Variation tendency comparison of optimization results with higher frequency of nodal diameter vibration in different iteration steps

4.2.3 提高節圓振動頻率

明確式(2)中的頻率約束，提高節圓振動頻率約束為f4≥700，優化結果如圖12所示.從圖中可以看出，與單獨拓撲優化結果類似，拓撲和形狀同時優化得到了類似于雙幅板形式的結構.二者不同的是，同時優化結果的幅板之間的間隔更小.

圖12 提高節圓振動頻率的優化結果對比Fig.12 Comparison of optimization results with higher frequency of umbrella vibration

圖12(d)所示的雙幅板渦輪盤拓撲形式，經過完善約束條件，以及在該優化結果基礎上進一步的形狀和尺寸優化，是能夠應用于工程實踐的.美國空軍 CRRT計劃［5］中采用了類似于圖12(d)(如圖12)的雙幅板渦輪盤而取得了明顯效果.

圖13分別列出了拓撲和形狀同時拓撲優化，以及單獨拓撲優化的應變能、頻率、體積隨迭代步數變化情況.從圖中可以看出，兩種優化的體積約束在最后迭代步均達到相應上限值，而優化的目標——應變能則不同，拓撲和形狀同時優化結果的應變能較單獨拓撲優化的應變能更小，減小了約28.3%，并且其收斂速度較快.

4.2.4 降低節圓振動頻率

明確式(2)中的頻率約束，降低節圓振動頻率約束為f4≤500，優化結果如圖14所示.

圖13 提高節圓振動頻率的優化結果隨迭代步數變化趨勢對比Fig.13 Variation tendency comparison of optimization results with higher frequency of umbrella vibration in different iteration steps

圖14 降低節圓振動頻率的拓撲和形狀同時拓撲優化結果Fig.14 Results of simultaneous topology and shape optimization with lower frequency of umbrella vibration

形狀優化變量使得初始設計模型中輪盤中心孔處的厚度不斷變化，但是如果厚度過小，如圖14(b)所示，會使該處單元的邊長比過大，而導致在優化過程中由于網格畸變過大，使得優化過程提前終止，這樣的優化結果是無效的.圖15描述了此時中心孔處單元邊長比達到117，超過了上限值(該值等于100).

圖15 邊長比過大的單元Fig.15 Elements with overlarge aspect ratio

為了得到合適的優化結果，需修改形狀優化設計變量的取值范圍，由前述的［－80，80］改為［－25，25］，使用同時優化方法得到的結果如圖16(c)、圖16(d)所示.與單獨拓撲優化方法所得結果類似，均為單幅板形式的拓撲.

圖16 降低節圓振動頻率的優化結果對比Fig.16 Comparison of optimization results with lower frequency of umbrella vibration

使用同時優化方法所得結構的應變能較單獨拓撲優化方法所得結構，減小3.25%，如圖17所示.

圖17 降低節圓振動頻率優化結果的應變能隨迭代步數變化趨勢對比Fig.17 Variation tendency comparison of strain energy of optimization results with lower frequency of umbrella vibration in different iteration steps

總的來看，STSO方法的優點是系統、準確、收斂快.它的局限性除了易引起網格畸變過大外，還在于將STSO方法應用于輪盤優化中時，并不能直接將優化結果應用于工程實踐中.這是由于STSO方法只能包含初始設計模型中的形狀優化變量，對于STSO方法的優化結果中新出現的拓撲構型，依然需要設置相應的形狀優化變量，進行進一步的形狀和尺寸優化等細化設計.因此，從這個方面看，STSO方法只是對應于結構的概念設計階段，如圖2(c)所示.

從概念設計階段的結果到工程實際應用的輪盤，需要再考慮的因素較多，比如滿足其他靜力學約束，降低某階次頻率同時、提高另外階次頻率，熱應力約束等.本文僅僅是運用STSO方法考慮了振型對應的頻率約束對同時優化方法所得結果的影響，得到了幾種概念設計方案，在本文所用STSO方法得到的優化結果的基礎上，還需要進一步的形狀和尺寸優化.

5 結論

本文首先建立了STSO方法的數學模型，并推導了目標和約束函數的靈敏度，運用SQP算法進行求解.用板殼結構優化的例子對比分析了同時優化和分步優化結果.將同時優化的方法用于輪盤動力學優化，分析了不同振型時頻率約束對同時優化方法所得結果的影響，并將之與單獨拓撲優化結果進行比較，結果表明:

1)拓撲和形狀同時優化方法較分步優化減少了人為參與的環節，從優化過程看，更為客觀科學.

2)拓撲和形狀同時優化考慮了形狀優化變量和拓撲優化變量之間的相互作用，回避了文獻中所提出的拓撲和形狀優化先后順序問題，從整體看，更為系統化.

3)將拓撲和形狀同時優化的方法應用于板殼結構優化和實體輪盤結構動力學優化時，結果表明，相較于分步拓撲優化和單獨拓撲優化，同時優化方法的結果較分步優化效果更好，且收斂更快.

4)相較于輪盤單獨拓撲優化方法結果，使用STSO方法得到的結果開孔較少甚至不開孔，更加符合現役發動機少有開孔輪盤的現狀.

5)不同振型對應的頻率約束不同時，得到了不同的優化結果.降低節徑振動頻率，均得到了輻條式拓撲;提高節圓振動頻率，均得到了雙幅板式拓撲;降低節圓振動頻率，均得到了單幅板式拓撲.

6)由于形狀優化變量取值范圍選取不當，可能導致網格畸變過大而得到無效優化解，此時需要重新選擇形狀優化變量的取值范圍.

References)

［1］ Deaton J D，Grandhi R V.A survey of structural and multidisciplinary continuum topology optimization:post 2000［J］.Struct Multidisc Optim，2014，49(1):1-38.

［2］ Ansola R，Canales J，Tarrago J A，et al.An integrated approach for shape and topology optimization of shell structures［J］.Computer ＆ Structures.2002，80(5):449-458.

［3］ Hassani B，Tavakkoli S M，Ghasemnejad H.Simultaneous shape and topology optimizationof shell structures［J］.Struct Multidisc Optim.2013，48(1):221-233.

［4］ 張衛紅，楊軍剛，朱繼宏.壓力載荷下的結構拓撲-形狀協同優化［J］.航空學報，2009，30(12):2335-2341.Zhang W H，Yang J G，Zhu J H.Simultaneous topology and shape optimization of pressure loaded structures［J］.Acta Aeronautica et Astronautica Sinica，2009，30(12):2335-2341(in Chinese).

［5］ 劉超.渦輪盤結構拓撲與形狀優化方法研究［D］.南京:南京航空航天大學，2010.Liu C.Topology and shape optimization method for turbine disk［D］.Nanjing:Nanjing University of Aeronautics and Astronautics，2010(in Chinese).

［6］ 宋健，溫衛東，崔海濤，等.航空發動機多輻板風扇盤拓撲優化與形狀優化設計技術［J］.推進技術，2013，34(9):1188-1196.Song J，Wen W D，Cui H T，et al.Topology and shape optimization method for multi-web fan disk in aero-engine［J］.Journal of Propulsion Technology，2013，34(9):1188-1196(in Chinese).

［7］ 李倫未，陸山.基于ANSYS的多輻板風扇盤結構優化設計技術［J］.航空動力學報，2011，26(10):2245-2250.Li L W，Lu S.Structure optimum design techniques for multi-web fan disk based on ANSYS［J］.Journal of Aerospace Power，2011，26(10):2245-2250(in Chinese).

［8］ 榮見華，唐國軍，羅銀燕，等.考慮位移要求的大型三維連續體結構拓撲優化數值技術研究［J］.工程力學，2007，24(3):20-27.Rong J H，Tang G J，Luo Y Y，et al.A research on the numerical topology optimization technology of large three-dimensional continuum structures considering displacement requirements［J］.Engineering Mechanics，2007，24(3):20-27(in Chinese).

［9］ 左孔天，陳立平，鐘毅芳，等.基于人工材料密度的新型拓撲優化理論和算法研究［J］.機械工程學報，2004，40(12):31-37.Zuo K T，Chen L P，Zhong Y F，et al.New theory and algorithm research about topology optimization based on artificial material density［J］.Chinese Journal of Mechanical Engineering，2004，40(12):31-37(in Chinese).

［10］ 鄒春江，左孔天，向宇，等.基于SIMP方法微電容加速度計結構固有頻率拓撲優化［J］.科學技術與工程，2011，11(29):7086-7091.Zhou C J，Zuo K T，Xiang Y，et al.Natural frequency topology optimization for structure of micro-capacitive accelerometer based on SIMP method［J］.Science Technology and Engineering，2011，11(29):7086-7091(in Chinese).

［11］ 葉紅玲，沈靜嫻，隋允康.頻率約束的三維連續體結構動力拓撲優化設計［J］.力學學報，2012，44(6):1037-1045.Ye H L，Shen J X，Sui Y K.Dynamic topological optimal design of three-dimensional continuum structures with frequencies constraints［J］.Chinese Journal of Theoretical and Applied Mechanics，2012，44(6):1037-1045(in Chinese).

［12］ 張維聲.基于水平集方法的應力相關拓撲優化問題研究［D］.大連:大連理工大學，2013.Zhang W S.Studies on stress-related topology optimization problem via level set approach［D］.Dalian:Dalian University of Technology，2013(in Chinese).

［13］ 隋允康，葉紅玲，彭細榮.應力約束全局化策略下的連續體結構拓撲優化［J］.力學學報，2006，38(3):364-370.Sui Y K，Ye H L，Peng X R.Topological optimization of continuum structure under the strategy of globalization of stress constraints［J］.Chinese Journal of Theoretical and Applied Mechanics，2006，38(3):364-370(in Chinese).

［14］ 榮見華，葛森，鄧果，等.基于位移和應力靈敏度的結構拓撲優化設計［J］.力學學報，2009，41(4):518-529.Rong J H，Ge S，Deng G，et al.Structure topological optimization based on displacement and stress sensitivity analysis［J］.Chinese Journal of Theoretical and Applied Mechanics，2009，41(4):518-529(in Chinese).

［15］ Paris J，Navarria F，Colominas I，et al.Stress constraints sensitivity analysis in structural topology optimization［J］.Computer Methodsin Applied Mechanicsand Engineering，2010，199(133):2110-2122.

［16］ Kai-Uwe B，Matthias F，Fernass D.Approximation of derivatives in semi-analytical structural optimization［J］.Computers ＆Structures，2008，86(13):1404-1416.

［17］ De Boer H，Van Keulen F.Refined semi-analytical design sensitivities［J］.International Journal of Solids and Structures，2000，37(46):6961-6980.

［18］ Van Keulen F，Haftka R T，Kim N H.Review of options for structural design sensitivity analysis.Part 1 linear systems［J］.Computer Methods in Applied Mechanics and Engineering，2005，194(30):3213-3243.

［19］ Fletche R.Practical methods of optimization［M］.New York:John Wiley and Sons，1987:304-319.

［20］ Nocedal J，Wright.S J.Numerical optimization［M］.New York:Springer Series in Operations Research，1999:529-561.