吳貝尼，夏利娟

（1.上海交通大學海洋工程國家重點實驗室，上海200240；2.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240）

0 引 言

結構拓撲優(yōu)化是在給定的設計空間里確定滿足給定約束的最優(yōu)材料分布以求獲得最佳結構性能。在過去幾十年中，各種優(yōu)化方法已被成功運用到各類工程結構拓撲優(yōu)化設計中，包括均勻化方法[1]、變密度法[2-4]、雙向漸進結構優(yōu)化方法（BESO）[5-7]、獨立連續(xù)映射方法[8]和水平集方法[9]等。其中，BESO方法原理簡單，收斂性好，通用性強，易于編程實現(xiàn)，因此備受青睞。

BESO 方法分為硬刪除（hard-kill）方法[10]和軟刪除（soft-kill）方法[11-12]，前者直接刪除低效單元，因此計算效率高，而后者引入帶懲罰項人工材料插值模型，并非完全移除單元，而是賦予其非常低的材料屬性。兩種BESO 方法（本文中稱為標準BESO 方法）在優(yōu)化求解的過程中同時考慮刪去低效單元和增添高效單元，有很大幾率獲得最優(yōu)結構，但是二者增刪單元是依據某種優(yōu)化準則的，每次迭代優(yōu)化受控于進化率，從而可能導致誤刪最優(yōu)拓撲結構中的高效單元而無法恢復，最終得到的只是局部最優(yōu)解。

遺傳算法基于適者生存思想，有很強的全局尋優(yōu)能力，而且遺傳算法中優(yōu)勝劣汰的思想十分契合結構拓撲優(yōu)化中高效單元保留、低效單元淘汰的原則。因此，將遺傳算法結合到BESO方法中，能有效克服BESO方法的缺陷，進一步提高優(yōu)化求解質量。

基于此，本文提出了一種基于改進遺傳算法的軟刪除雙向漸進結構優(yōu)化法（G-BESO），采用一種考慮單元密度歷史信息的材料插值模型以增強恢復誤刪高效單元的能力，利用遺傳算法中的交叉和變異操作啟發(fā)式更新結構狀態(tài)以提高全局尋優(yōu)能力。結合有限元平衡方程建立了連續(xù)體結構優(yōu)化模型并給出具體求解流程，最后用經典算例驗證了方法的有效性和穩(wěn)健性，同時對比標準BESO 方法表明了其優(yōu)越性。

1 基于改進遺傳算法的雙向漸進結構優(yōu)化模型構建

1.1 算法思想

首先利用有限單元將連續(xù)體結構離散化，設計區(qū)域為種群，離散單元為個體，以離散的單元密度值作為設計變量（0 或1），其中0 代表孔洞單元，1 表示實體單元；對于個體（單元）使用二進制編碼，編碼中1的數目代表著其重要程度即適應度。適應度決定個體編碼的更新方式，個體之間通過交叉、變異更新其編碼，編碼決定個體保留與淘汰。通過這種方式，實現(xiàn)種群的進化即結構的優(yōu)化。

標準BESO方法中每一步迭代優(yōu)化通常需要刪除由進化率決定的固定體積的單元，完全根據靈敏度大小機械地增添單元。本文提出的G-BESO 方法基于優(yōu)勝劣汰原則，根據由適應度決定的單元編碼中1的數量啟發(fā)式增添不定體積的單元。具體而言，若單元為實單元，且其編碼值均為0，則從結構中刪除；若為空單元，且其編碼中1的數量達到一定值，則添加到結構中。

1.2 數學模型

本文研究連續(xù)體結構在給定體積約束下的最大剛度設計問題，拓撲優(yōu)化模型可以描述為

確定設計變量：

式中：F、u分別是載荷向量和位移向量；K是剛度矩陣；V*是給定的體積約束；xi是設計變量，其值為0或1，表征空/實單元。

本文提出的G-BESO 方法是一種軟刪除方法，采用了一種新的考慮單元密度歷史信息的材料插值模型。引入文獻[13]中單元權重系數的概念，結合SIMP 材料插值模型，建立了單元權重系數I與單元密度X、單元材料屬性E的關系式：

1.3 靈敏度計算與收斂準則

單元靈敏度值代表其對整個結構剛度貢獻的大小，關系著單元的刪除與保留，決定最終結構拓撲形狀。對于剛度最大優(yōu)化問題，基于所提出的材料插值模型的單元應變能靈敏度為

式中，ui，Ki，Ii分別是單元i的位移向量、剛度矩陣和權重系數。

為解決棋盤格問題和網格依賴性，本文采用文獻[10]中的靈敏度過濾方法：

式中，N為過濾范圍內節(jié)點總數，rmin為過濾半徑，rij為單元i中心與節(jié)點j之間的距離。

同時，為穩(wěn)定優(yōu)化過程，當前單元靈敏度需要根據上一步所得值進行修正[11]：

結構拓撲優(yōu)化是一個不斷更新單元狀態(tài)的迭代過程，直到滿足體積約束和收斂準則，迭代終止。本文采用文獻[11]中推薦的收斂終止公式：

式中：τ為收斂因子即迭代精度，本文取0.001；k為迭代步；N為常數，一般取為5。表征結構的應變能在連續(xù)10次以上的迭代中變化量足夠小即認為得到穩(wěn)定結構。

2 基于改進遺傳算法的雙向漸進結構優(yōu)化求解流程

2.1 技術手段

為提高G-BESO 方法的實用性，使之能用于優(yōu)化復雜的工程結構，本文以Python 為編程手段，通過調用HyperMesh軟件更新有限元模型，調用Nastran求解器進行有限元計算，可實現(xiàn)對實際工程結構的拓撲優(yōu)化設計。

2.2 求解步驟

基于上一章構建的拓撲優(yōu)化模型，給出G-BESO方法的具體求解步驟如下：

（1）使用HyperMesh對連續(xù)體結構進行網格劃分，定義設計區(qū)域，施加邊界條件和載荷。

（2）利用Nastran進行結構有限元分析，計算過濾后的單元靈敏度。

（3）根據靈敏度將單元劃分為三個集合：高效單元集、過渡單元集和低效單元集。若為第一次進化，對設計區(qū)域單元進行編碼：對所有單元賦予少量0、較多1 的編碼。計算進度指標Prg，交叉概率Pc，變異概率Pm。

式中，Vcur，V*分別為結構當前體積分數和約束體積分數。

為解決遺傳算法隨機性引發(fā)的跳躍解、非可行解（結構不協(xié)調）等問題，參照文獻[14]，引入懲罰因子如下，使得Pc與Pm隨結構優(yōu)化逐漸增大：

式中，Penc，Penm分別是預先選定的交叉和變異懲罰因子。

（4）個體（單元）之間實施交叉、變異以更新單元編碼。

類似于精英保留策略，高效單元不作為父代（但作為母代）參與交叉直接進入下一代，以避免高效單元意外喪失導致優(yōu)化陷入局部最優(yōu)解。其他單元作為父代隨機選擇單元（母代）進行一次配對交叉，產生的子代之一進入下一代，方式為兩點交叉。選擇同一集合中的單元配對交叉的概率為Pc，選擇其他集合中單元的概率為( 1 -Pc)/2。

變異操作只作用于高效單元集和低效單元集，直接改變單元二進制編碼中的數位值[14]。對于每個數位來說，其變異的概率設置為Pm，如果產生的隨機數值小于Pm，高效單元編碼中0 變?yōu)?，低效單元中的1變?yōu)?。這樣，變異可以使高效單元編碼中含有更多的1，而低效單元編碼中含有更多的0，結構逐步優(yōu)化。

（5）根據單元編碼中1的數量確定單元狀態(tài)，更新設計變量，調用HyperMesh修改結構。

盡管遺傳算法能以較大概率找到全局最優(yōu)解，由于其固有的隨機性，仍可能限于局部最優(yōu)。為控制每一次進化結果，每一步交叉變異更新單元編碼的過程進行若干次（10次），根據應力離散度[15]進行評估，從中預選出最優(yōu)解，更新設計變量。需要指出的是，此處不需要進行有限元計算，只需進行的簡單代數運算并不耗時。

應力離散度表征結構受力的均勻程度。隨著優(yōu)化的進行，結構中低效單元逐漸刪除，高效單元得以保留，應力離散度值越來越低。因此每次迭代選擇應力離散度最小的種群為當前種群參與后續(xù)迭代優(yōu)化。式中，N表示實單元數量，σVMi為實單元i的相當應力，σVMmax為結構最大相當應力。

（6）判斷體積約束和收斂條件是否滿足，若滿足則終止優(yōu)化，否則返回至第（2）步。

G-BESO方法流程如圖1所示。

圖1 G-BESO方法流程圖Fig.1 G-BESO flow chart

3 算例分析

3.1 L型支架優(yōu)化

圖2所示為經典L 型支架優(yōu)化問題。支架的上端固支，右側頂端受豎直向下的集中力F=4 N。材料彈性模量E=206 GPa，泊松比ν=0.3。體積分數約束為0.4。將結構劃分為1 600個4×4 mm 四邊形單元。過濾半徑rmin=12 mm，交叉概率Pc?min=0.6，Pc?max=1.0，變異概率Pm?min=0.6，Pm?max=1.0，權重因子δ=0.5。為便于與標準BESO 法進行對比，相同參數設為一致，BESO法的其他參數設置為：進化率取1%～5%，最大添加率6%。為驗證G-BESO方法的穩(wěn)健性，連續(xù)計算5次。

圖2 L型支架初始結構Fig.2 L-bracket initial design

由于標準BESO 方法中的兩種方法（硬刪除方法和軟刪除方法）的計算結果類似，為簡潔計算，以下均僅列出前者的結果。為便于描述，標準BESO 方法中的硬刪除BESO 方法簡稱為H-BESO，采用本文提出的材料插值模型的軟刪除BESO 方法簡稱為S-BESO，本文提出的基于改進遺傳算法的雙向漸進結構優(yōu)化方法簡稱為G-BESO。

三種方法所得優(yōu)化構型和結果數據如表1和表2所示。

表1 不同優(yōu)化方法對L型支架優(yōu)化所得拓撲構型Tab.1 Topological layouts of L-bracket based on different methods

表2 不同優(yōu)化方法對L型支架優(yōu)化結果數據（柔順度單位：×10-3 N?mm）Tab.2 Optimization data of L-bracket based on different methods(Compliance unit:×10-3 N?mm)

對比S-BESO和H-BESO方法，隨著進化率由1%增大到5%，前者均能得到最優(yōu)拓撲構型，而后者在進化率增大到5%時優(yōu)化失?。ū?），這體現(xiàn)出本文提出的考慮材料歷史信息的插值模型的優(yōu)勢。該插值模型基于單元密度-單元權重系數-材料彈性模量插值，單元的密度是逐步減小而非（在標準BESO 方法中）從1 直接減到0 或某極小值（0.001），且部分誤刪的高效單元也能通過密度正向生長而重新被添加到結構中，因此單元誤刪率更低。

此外，本文提出的基于改進遺傳算法的G-BESO 方法不僅在5 次運行中均收斂到最優(yōu)結構，而且平均柔順度8.395 9×10-3N?mm 也優(yōu)于其他兩種方法所得值，證明了其穩(wěn)健性和有效性；平均迭代步數44.4明顯少于另外兩種方法所需步數，表明其優(yōu)化效率也更高。

3.2 短懸臂梁優(yōu)化

圖3所示為經典短懸臂梁算例。梁的左端固支、右端中點承受豎直向下集中力F=100 N。梁的尺寸為長L=120 mm，寬H=60 mm，厚度t=1 mm；材料彈性模量E=206 GPa，泊松比ν=0.3。體積分數約束為0.5。將結構劃分為60×30 個四邊形單元。過濾半徑rmin=6 mm，其他參數均與上例一致。

從最終優(yōu)化得到的拓撲構型（表3）上看，H-BESO 方法受進化率的影響較大，5次優(yōu)化下只有2次（3%、4%）得到最優(yōu)構型，S-BESO 方法受進化率的影響相對較小，3次（2%、4%、5%）得到最優(yōu)構型，再次表明本文建立的材料插值模型的優(yōu)勢；受益于遺傳算法的全局尋優(yōu)能力，G-BESO 方法5 次運行均得到最優(yōu)構型，表明G-BESO 方法的穩(wěn)健性和有效性；從計算效率和優(yōu)化目標值（表4）上看，G-BESO方法平均迭代步數38.6、平均柔順度1.5516 N?mm均優(yōu)于其他兩種方法。

圖3 短懸臂梁初始結構Fig.3 Short cantilever initial design

表3 不同優(yōu)化方法對短懸臂梁優(yōu)化所得拓撲構型Tab.3 Topological layouts of short cantilever based on different methods

續(xù)表3

表4 不同優(yōu)化方法對短懸臂梁優(yōu)化結果數據（柔順度單位：N?mm）Tab.4 Optimization data of short cantilever based on different methods(Compliance unit:N?mm)

4 結 論

本文提出的G-BESO 方法綜合了漸進結構優(yōu)化法，能得到簡單清晰的拓撲構型以及遺傳算法全局尋優(yōu)能力強的優(yōu)勢，計算結果表明：

（1）建立的基于單元密度-單元權重系數-材料彈性模量插值的材料插值模型，考慮了單元密度的歷史信息，相較標準BESO方法單元誤刪率更低，優(yōu)化結果更好。

（2）對基本遺傳算法所做的改進及與BESO 方法的結合，賦予G-BESO 方法更強的全局尋優(yōu)能力、更高的計算效率。

（3）G-BESO方法能穩(wěn)定得到最優(yōu)拓撲構型而且計算效率高，具有應用到實際工程結構優(yōu)化設計之中的潛力和價值。