基于自適應成長法的周期性加筋結構拓撲優化設計方法

董小虎 丁曉紅

上海理工大學機械工程學院，上海，200093

0 引言

周期性加筋結構廣泛應用于航天航空、船舶和汽車等領域，具有易加工、制造成本低、外形美觀等特點。為了得到力學性能良好的周期性結構，拓撲優化目前已成為上述領域中最常用的設計方法。劉書田等［1］研究了狹長結構的拓撲優化，提出了通過基本結構拼裝的狹長結構的拓撲優化方法，設計了狹長懸臂梁的拓撲形式。HUANG等［2］利用雙向漸進結構優化（bi-directional evolutionary structural optimization，BESO）方法對周期性結構優化設計進行了研究，由于數值計算不穩定，優化過程中易出現棋盤格、網格依賴性等問題，因此該方法引入了敏度過濾函數來抑制棋盤格和網格依賴性，從而獲得了較好的優化結果。榮見華等［3］提出了一種基于位移約束的類周期性連續體結構的拓撲優化設計方法，采用改進的對偶方法，得到了優化的類周期性結構拓撲形式。ZUO等［4］以圓環結構為對象，采用漸進結構優化（evolutionary structural optimization，ESO）方法對汽車輪轂進行拓撲優化設計。HE等［5］依據Shepard插值法，建立了單元與非結構設計點的關系，將BESO方法應用到不同的幾何尺寸和無規則有限元網格的周期性結構拓撲優化中。焦洪宇等［6-7］提出了一種基于固體各向同性材料懲罰（solid isotropic microstructure with penalization，SIMP）模型的周期性拓撲優化方法，設置了額外的約束條件，對平面矩形懸臂梁結構進行了研究，獲得了具有周期性的結構拓撲形式，并將該方法應用到橋式起重機主梁周期性拓撲優化中。孫士平等［8］針對狹長結構拓撲優化過程中存在的計算困難、數值不穩定等問題，利用拓撲優化方法和變量關聯技術開展了周期狹長結構拼裝模塊的構型優化。上述拓撲優化方法雖然可獲得較好的結構材料分布形式，但其拓撲形態往往過于復雜，若應用到加強筋分布設計中，無法直接獲得加強筋的實際分布形式，需通過后處理對加強筋分布進行識別，因此不易于直接用于加工制造。

仿生優化研究隨著計算機技術的普及而得到發展，在結構設計領域也得到了一定的應用，但人們直接針對周期性結構的研究還較少。丁曉紅等［9-11］通過研究自然界分枝系統的形成機理，提出了一種基于植物根系形態自適應成長規律的板殼結構加強筋設計的自適應成長法，并將該方法成功應用于以最小應變能和抗振設計為目標的板殼加筋優化設計中，通過改變加強筋高度，逐步實現加強筋的布局設計，從而獲得了離散的加強筋分布形式，且不需進行后處理。

本文通過在自適應成長法中引入周期性約束，提出了周期性結構的加筋布局設計方法，建立了以結構最小應變能為設計目標、以體積為約束的優化設計數學模型，并引入了敏度過濾函數解決數值不穩定等問題。同時以荷葉葉脈的分布為設計對象，得到與自然界荷葉葉脈相似的加筋分布形態，并進一步通過典型工程結構的設計案例，驗證了所提方法的有效性。

1 設計方法

1.1 周期性約束

為獲得具有周期性分布的結構最優拓撲形態，將整個優化區域劃分成P個子域，見圖1，其中圖1a為圓環結構，根據圓的周長等分成P=4個子域，圖1b為長條狀結構，沿著軸向等分成P=4個子域。對于周期性結構，子域與子域之間存在相同的子拓撲構型，但兩子域的拓撲尺寸大小可以不同，域設計數目P可由設計者根據需要自定。

圖1 優化區域分成P個子域（P=4）Fig.1 The optimized region is divided into P sub-domains（P=4）

1.2 優化數學模型

以整體結構應變能U最小為目標，以加強筋截面積A為設計變量，在給定的體積約束條件下，尋求結構力學性能最優。優化數學模型可表示為

式中，Ai為第i個活動筋板單元的截面積；n為筋板單元的總數；U（A）為結構應變能；V、V0分別為周期性加筋結構的總體積和初始體積；η為體積約束因子；Amin、Amax分別為設計變量Ai的取值下限和上限；t為某個子域編號；j為該子域單元編號；m為每一個子域內筋板單元的總數。

1.3 基結構

在初始設計模型中，要建立圖2所示的基結構，設計域采用板殼單元離散，用梁單元連接各個節點以模擬結構內的加強筋。初始截面高度h和寬度w均為極小值，即Amin為極小值（約為結構長度L或直徑D的萬分之一）。

圖2 基結構模型Fig.2 Model of the ground structure

1.4 迭代公式

基于 KKT（Karush-Kuhn-Tucker）最優化條件推導出筋板成長公式，建立拉格朗日方程如下：

式中，λ為拉格朗日乘子。

在最優A*處，式（2）應滿足KKT必要條件。

當Amin＜ Ai＜ Amax時，根據式（3）得到

總應變能對設計變量的靈敏度為

整體結構的總體積V可表示為

式中，li為第i個筋板單元的長度。

將式（5）和式（7）代入式（4）可得

為了保證收斂性，引入步長因子α，取α=0.05，得到設計變量Ai的迭代公式如下：

式中，k為迭代步數。

拉格朗日乘子的表達式如下：

1.5 敏度過濾函數

拓撲優化過程中經常存在數值不穩定問題，如出現棋盤格、灰度單元等現象。SIGMUND［12］提出的敏度過濾法具有實施簡單、計算效率高等優點，被廣泛應用于拓撲優化設計中。該方法對所有在給定過濾半徑范圍內的各單元靈敏度進行加權平均，用來修正每個單元的靈敏度。目標函數對子域內設計變量Ai修正后的單元靈敏度為

式中，Hq為權重因子；Rmin為過濾半徑，取Rmin=2.5；ds（i，q）為過濾半徑內單元q至中心單元i的距離；N為過濾半徑內單元的總數目。

2 荷葉葉脈再生過程仿真

2.1 荷葉葉脈分支系統形態特征

圖3a中，荷葉葉片近似圓形，其中葉脈是莖分支結構的延伸，分布于整個葉片的網絡系統分支結構中，形成一個承載葉片的骨架，具有支撐整個葉片的功能，同時也最大化地展開了葉片，增加了葉片截獲太陽光的面積［13］。荷葉葉脈具有較強的韌性和硬度，能夠提高葉片對外界破壞（如小動物踐踏和暴雨冰雹等自然災害等）的抵抗能力。由此可見，葉脈分布直接影響荷葉對外界環境變化的適應能力。

葉脈經過逐級分支，形成網絡狀分叉分布結構見圖3b，葉脈從中央射出呈發散狀，每一個發散的分支葉脈在構成形態上幾乎相同，具有較強的自相似性。其中主脈的結構最完整且直徑最大，連通葉柄起主要支撐作用。葉脈直徑可以用來表征葉脈系統的機械支撐功能。不同級別的葉脈承擔不同程度的支撐作用。盡管末端葉脈也有一定的支撐力，但較粗的葉脈發揮主導作用。葉脈分布形態是植物適應生態系統的關鍵和決定性因素。

圖3 荷葉葉脈分支系統形態特征Fig.3 Branch system morphology of lotus leaf vein

2.2 模擬荷葉葉脈分布形態

將基于自適應成長法的周期性加筋結構拓撲優化設計技術運用到荷葉葉脈的再生長過程中，以整體結構應變能最小為目標進行拓撲優化設計。

圖4a中，荷葉由一根葉柄支撐挺出水面，在圖4b所示的設計模型的中間點施加全約束，用均布載荷表示荷葉承受暴雨冰雹的情況。

圖4 模擬荷葉葉脈分布形態設計模型Fig.4 Design model of imitated lotus leaf vein

圖5 所示為模擬“葉脈”自適應再生過程，依據自適應成長法的選種方式，選取中心點約束處為“種子”位置。“葉脈”首先從中心點“種子”處向外散發，沿著使整體結構應變能最小的方向成長與退化，逐步抵抗承載，成長至分歧截面積Amax后發生“分歧”；然后與各“分歧點”相連的“葉脈”繼續成長、退化或“消失”。該過程反復循環，直至最終整體結構應變能滿足如下收斂條件：

式中，ε為收斂容差，取ε=1×10-4；U為整體結構總應變能。

對比圖6a與圖6b可知，模擬“葉脈”的分布形態與真實的荷葉葉脈分布形態相似，“主葉脈”從中心點的“種子”出發，徑向生長，至邊緣處分歧出若干細枝，以抵抗受到的彎曲載荷，最大化展開并支撐“葉面”。最終結構的總應變能與初始結構的應變能之比U/U0為1.67×10-4。由荷葉葉脈的再生過程可知，荷葉葉脈分布能使荷葉在承受自然載荷作用下獲得整體剛度最大。圖6c為設計迭代歷程圖，可以看出，隨著迭代步數的增加，由于初始應變能U0和初始體積V0不變，故結構總應變能開始急劇減小，同時總體積迅速增大；隨后因筋板的截面積不斷地“生長”與“退化”，結構總應變能的減小幅度逐漸變緩并趨于穩定，而體積比值也逐漸趨于穩定，最終不再變化。

圖5 “葉脈”的自適應再生過程Fig.5 Regrowth process of the lotus leaf

圖6 設計結果對比與分析Fig.6 Comparison and analysis of the design result

3 典型算例與分析

3.1 不同子域數目的圓環結構周期性筋板分布形態設計

首先考慮子域數目P=4時，承受指向圓心的均布面內載荷作用的圓環結構，見圖7a。設計模型尺寸為D∶d=4∶1，其中D為外圓直徑，d為內圓直徑。圖7b中，在4個約束點處選取“種子”，加強筋從“種子”點處開始發散成長，逐步連接至圓環的邊界，再依據自身的設計敏度過濾函數，形成由疏到密的加強筋布局形態，從而直接改善整個圓環結構力學性能。最終加強筋的設計分布結果見圖7c，此時結構的整體應變能與初始結構的應變能之比U/U0為1.27×10-2。圖7d為圓環結構的體積和應變能比值的迭代歷程曲線，可以看出，隨著迭代步數的增加，結構總應變能減小，總體積增大，2條曲線最終均趨于穩定。

圖7 P=4時周期性圓環結構的設計結果Fig.7 Design results for a periodic ring structure when P=4

圖8 a所示為子域數目P=5時承受指向圓心的均布面內載荷作用的圓環結構，其他設計參數均與P=4時的設計模型相同，“種子”位置見圖8b。圖8c所示為筋板最終的分布形態。圖8d為圓環結構的體積和應變能比值的迭代歷程曲線，可以看出，隨著迭代步數的增加，圓環結構的總體積增大，結構總應變能減小，在U/U0為1.19×10-2時收斂。

圖8 P=5時周期性圓環結構的設計結果Fig.8 Design results for a periodic ring structure when P=5

對比上述2個不同子域數目的加強筋板布局的設計結果可以看出，兩者每個子域的設計形態大致形同，起初加強筋板均從“種子”處開始發散成長，連接至結構邊緣處時，筋板數目逐漸增多，而圖7c的分支數量多于圖8c的分支數量，這是由于后者圓環的子域數目（P=5）比前者（P=4）多了一根主支。本例結果表明，自適應成長法可根據工程實際需要設置周期子域數目，以得到合理的加強筋分布形式。

3.2 不同幾何形狀的懸臂條狀結構周期性筋板分布形態設計

圖9a所示為平面矩形懸臂條狀結構設計模型，其中左端固定，右端作用有垂直載荷F。結構子域數目P=4，設計模型尺寸為L∶W=4∶1，其中L為矩形長，W為矩形寬。“種子”選取位置見圖9b。最終結構的總應變能與初始應變能比值U/U0為8×10-3。圖9c為圓環結構的體積和應變能比值的迭代歷程圖，可以看出，隨著迭代步數的增加，圓環結構的總體積增大，結構總應變能減小，2條曲線最終均趨于平緩。

圖9 平面矩形懸臂條狀結構的設計結果Fig.9 Design results for a planar rectangular cantilever bar structure

圖10 a所示為平面梯形懸臂條狀結構設計模型，結構子域數目P=4，設計模型尺寸為L∶W1∶W2=8∶2∶1，其中L為梯形的長，W1為梯形左端最大寬度，W2為梯形右端最小寬度。圖10b為“種子”選取位置與加強筋最終分布形態圖，其應變能比值U/U0收斂于1.49×10-2。圖10c為圓環結構的體積和應變能比值的迭代歷程圖，同理可以看出，隨著迭代步數的增加，結構總應變能減小，總體積增大，2條曲線最終均收斂。

圖10 平面梯形懸臂條狀結構的設計結果Fig.10 Design results for a planar trapezoid cantilever bar structure

上述2個例子考慮的是不同幾何形狀的懸臂結構加筋設計，周期子域數目P相同。雖設計模型的幾何形狀不同，但每個子域的筋板數目相等，均呈交叉型周期性分布。但隨著條狀結構沿軸向方向的寬度變化，加強筋交叉的角度會發生適應性變化。

4 結論

（1）通過引入周期性約束，將自適應成長法應用到周期性加筋結構的拓撲優化設計中。構建了優化設計模型，推導了基于KKT最優化條件的尋優迭代公式，并引入敏度過濾函數。

（2）通過對荷葉葉脈再生過程模擬表明，自然界中真實的荷葉葉脈確實能使荷葉獲得最大的整體剛度。

（3）針對周向受壓的圓環結構，分析了不同周期子域數目對設計結果的影響，并設計了不同形狀的懸臂條狀結構的周期性加強筋分布。

（4）設計結果表明，所提方法的迭代公式簡單，設計參數少，求解方便，筋板布局形態清晰，不需進行后處理即可得到筋板的實際布局，且可用于不同載荷作用下的各種結構，適用范圍廣。