面向雙向漸進(jìn)結(jié)構(gòu)拓?fù)鋬?yōu)化法的高效靈敏度過濾算法研究

尉洋，王春彭，陳興，曹月昊，姚松

(1.軌道交通安全教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長沙410075；2.軌道交通安全關(guān)鍵技術(shù)國際合作聯(lián)合實(shí)驗(yàn)室，湖南 長沙410075；3.軌道交通列車安全保障技術(shù)國家地方聯(lián)合工程研究中心，湖南 長沙410075)

自1988年Bends?e和Kikuchi提出均勻化拓?fù)鋬?yōu)化以來，拓?fù)鋬?yōu)化理論得到了巨大的發(fā)展[1]，各種基于梯度的算法也相繼提出，如變密度法(SIMP)[2-3]、水平集法(Level-set)[4]、雙向漸進(jìn)優(yōu)化法(BESO)[5-6]以及近年提出的浮動投影方法(FPTO)[7-8]等。隨著計(jì)算機(jī)性能的不斷提升以及增材制造技術(shù)的逐步完善，拓?fù)鋬?yōu)化已然成為結(jié)構(gòu)優(yōu)化設(shè)計(jì)領(lǐng)域中的熱點(diǎn)問題[9]。拓?fù)鋬?yōu)化的本質(zhì)是凸優(yōu)化問題，是在仿真分析的基礎(chǔ)上，找到材料在設(shè)計(jì)空間當(dāng)中的最優(yōu)分布。拓?fù)鋬?yōu)化技術(shù)應(yīng)用在結(jié)構(gòu)的概念設(shè)計(jì)階段，能夠突破傳統(tǒng)優(yōu)化手段的設(shè)計(jì)極限，最大限度地拓展設(shè)計(jì)空間[10]。在連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化過程中，經(jīng)常存在著數(shù)值不穩(wěn)定現(xiàn)象[11-12]，如棋盤格和網(wǎng)格依賴性等問題。前者會導(dǎo)致拓?fù)浣Y(jié)果可制造性降低，后者會使得計(jì)算結(jié)果不可靠。為了避免出現(xiàn)這些現(xiàn)象，可用的手段有：高階單元法、周長約束法、局部梯度約束法以及網(wǎng)格過濾技術(shù)[13-14]。其中，過濾技術(shù)不需要在優(yōu)化過程中加入其他約束，易于實(shí)現(xiàn)的同時(shí)計(jì)算穩(wěn)定，被認(rèn)為是最有前景的方法，得到了廣泛的應(yīng)用。網(wǎng)格過濾技術(shù)包括靈敏度過濾、密度過濾以及偏微分過濾，前2種是目前的主流過濾方法。靈敏度過濾是基于鄰域單元的靈敏度大小對中心單元靈敏度進(jìn)行修正，一般選用外接方形(二維)或外接立方體(三維)作為單元鄰域，然后進(jìn)行單元距離與權(quán)重的計(jì)算，如圖1所示，關(guān)鍵一步是確定待過濾單元的鄰域單元及其對應(yīng)的權(quán)重系數(shù)。在規(guī)則網(wǎng)格的背景下，F(xiàn)ERRARI發(fā)布的教育代碼中給出了依據(jù)單元規(guī)則編號確定鄰域單元及權(quán)重的方法[15]，具有較好的效率，但其適用對象僅限于單元和節(jié)點(diǎn)編碼按照預(yù)設(shè)規(guī)則設(shè)置的平面矩形結(jié)構(gòu)。對于復(fù)雜結(jié)構(gòu)，單元節(jié)點(diǎn)編碼不再有序，其不再適用。另一種方法是利用枚舉[16]，結(jié)合單元形心間距通過循環(huán)來確定單元鄰域，實(shí)現(xiàn)簡單方便。然而實(shí)際上絕大部分單元不在過濾半徑之內(nèi)，浪費(fèi)了大量的計(jì)算資源，其計(jì)算復(fù)雜度為O(n2)，對于具有大量網(wǎng)格的模型而言，枚舉的效率難以被接受。另外，當(dāng)結(jié)構(gòu)內(nèi)部中存在有細(xì)小非設(shè)計(jì)域孔洞時(shí)，僅僅通過單元形心距離確定中心單元鄰域的做法，無法準(zhǔn)確判斷單元之間的連通性以及識別孔洞邊界，以至于得到不準(zhǔn)確的單元鄰域。為了準(zhǔn)確快速找到過濾單元的鄰域，減少前處理階段中過濾算法的耗時(shí)，探尋解決細(xì)小非設(shè)計(jì)域孔洞問題的方法，針對BESO法，本文提出一種快速搜索策略，能夠適應(yīng)任意非規(guī)則網(wǎng)格且易拓展至三維。相對于枚舉策略而言，能夠顯著提高搜索效率。通過計(jì)算表明，在設(shè)計(jì)空間中存在微小非設(shè)計(jì)域孔洞的情況下，使用這種過濾策略的BESO法能得到性能更優(yōu)的拓?fù)錁?gòu)型。

圖1 靈敏度過濾技術(shù)Fig.1 Sensitivity filtering method

1 BESO法

以實(shí)際常見的受到體積約束的剛度最大問題為例，建立基于BESO法的求解模型如下：其中：C為結(jié)構(gòu)柔度；F為外載荷向量；U為節(jié)點(diǎn)位移向量；V*為結(jié)構(gòu)體積約束上限；N為單元個(gè)數(shù)；Vi為單元體積；K為結(jié)構(gòu)整體剛度矩陣。xi為設(shè)計(jì)變量，表示單元的相對密度，實(shí)體單元取值為1，為了保證整體剛度矩陣的非奇異性，空單元取較小值xmin，一般為0.001。

引入材料插值方案：

E0表示實(shí)體材料的彈性模量；p為懲罰系數(shù)，一般取值為3。通過伴隨法和鏈?zhǔn)角髮?dǎo)法則求得目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的靈敏度為：

其中：K0i為實(shí)體狀態(tài)下的單元剛度矩陣。

優(yōu)化迭代停止條件為：在過去10次迭代中，結(jié)構(gòu)柔度的變化足夠小。

2 利用節(jié)點(diǎn)單元從屬關(guān)系確定單元鄰域

式(3)所求得的單元靈敏度在直接使用時(shí)會存在著顯著的數(shù)值不穩(wěn)定現(xiàn)象，使得結(jié)構(gòu)中出現(xiàn)棋盤格或網(wǎng)格依賴性等問題。為了避免這些現(xiàn)象出現(xiàn)，需要進(jìn)行靈敏度過濾。

靈敏度過濾公式為：

其中：Hij為權(quán)重因子，計(jì)算公式為：

d ist(i,j)表示單元i的中心到單元j中心的距離；rmin是過濾半徑。

對于拓?fù)鋬?yōu)化過濾技術(shù)而言，最為關(guān)鍵且耗時(shí)最長的部分在于確定單元的搜索鄰域N，本文提出一種利用從中心單元——從屬單元——鄰域單元的逐層搜尋模式快速準(zhǔn)確地搜尋求解。

當(dāng)結(jié)構(gòu)建模、網(wǎng)格劃分完成之后，有限元模型中會儲存有表示節(jié)點(diǎn)位置的節(jié)點(diǎn)-坐標(biāo)信息、單元所有節(jié)點(diǎn)的單元-節(jié)點(diǎn)信息E={n1,…,nNe},通過E可以方便地推導(dǎo)出節(jié)點(diǎn)所在單元的節(jié)點(diǎn)-單元信息N={e1,…,eNn}。其中，ni表示單元i中節(jié)點(diǎn)的編號(i=1,…,Ne，Ne為結(jié)構(gòu)中單元數(shù))；ej表示節(jié)點(diǎn)j所屬單元的編號(j=1,…,Nn，Nn為結(jié)構(gòu)中節(jié)點(diǎn)數(shù))。定義f為e到E的一個(gè)映射，g為n到N的一個(gè)映射，即：

其中：unique為去除集合中重復(fù)單元。

以如圖2的二維四邊形單元為例。首先確定中心單元ei，通過單元-節(jié)點(diǎn)信息找到單元ei包含的節(jié)點(diǎn)(ni1,ni2,ni3,ni4)，通過節(jié)點(diǎn)-單元信息找到擴(kuò)充后的單元編號，如式(6)所示。循環(huán)所有設(shè)計(jì)域單元組合后得到每個(gè)單元的相鄰單元信息，至此構(gòu)建出中心單元與從屬單元之間的聯(lián)系，如式(7)所示。當(dāng)過濾半徑給定時(shí)，反復(fù)索引調(diào)用式(7)使得搜索范圍向外輻射，通過自適應(yīng)收斂條件能夠針對不同網(wǎng)格密度區(qū)域設(shè)置合理的迭代層數(shù)，便能迅速確定單元鄰域，最終構(gòu)建出從屬單元到鄰域單元之間的聯(lián)系，如式(8)和式(9)所示。收斂得到的鄰域單元去除重復(fù)之后即為該單元的最終搜索鄰域。

圖2 快速尋找單元鄰域Fig.2 Fast search elements neighborhoods

網(wǎng)格過濾算法來源于數(shù)字圖像處理中的降噪技術(shù)，應(yīng)當(dāng)考慮結(jié)構(gòu)或單元之間的連通性。當(dāng)結(jié)構(gòu)中出現(xiàn)有細(xì)小非設(shè)計(jì)域孔洞時(shí)，傳統(tǒng)過濾策略若要準(zhǔn)確識別，需要在搜索過程中反復(fù)調(diào)用和讀取結(jié)構(gòu)的實(shí)際邊界信息，代碼實(shí)現(xiàn)繁瑣，顯著增加編程難度。一般的做法是按單元形心間距離進(jìn)行搜索，不考慮孔洞所帶來的影響，如圖3(a)所示。

圖3 含微小非設(shè)計(jì)域孔洞時(shí)傳統(tǒng)策略與新策略的比較Fig.3 Comparison of traditional and new strategies when there are small non-design domains holes

對于這種情況，新策略有天然的優(yōu)勢：有限元模型中的單元-節(jié)點(diǎn)信息和節(jié)點(diǎn)-單元信息已經(jīng)包含了結(jié)構(gòu)實(shí)際邊界信息，當(dāng)搜索到幾何邊界時(shí)，若要繼續(xù)進(jìn)行，由于孔洞沒有單元，因此搜索區(qū)域不會越過細(xì)小孔洞，從而得到相對正確的單元鄰域，如圖3(b)所示。

3 新策略靈敏度過濾流程

圖4展示了運(yùn)用新策略進(jìn)行靈敏度過濾的詳細(xì)流程圖。

圖4 新策略靈敏度過濾流程圖Fig.4 New strategy sensitivity filtering flow chart

4 過濾策略的比較

本節(jié)從網(wǎng)格適用性、過濾用時(shí)以及帶初始設(shè)計(jì)孔洞3個(gè)方面比較幾種過濾策略的優(yōu)劣，均選擇使用Python語言在ABAQUS平臺上進(jìn)行二次開發(fā)，ABAQUS版本為2016，拓?fù)鋬?yōu)化算法為BE‐SO法。

4.1 網(wǎng)格適用性

依據(jù)規(guī)則編號確定單元鄰域具有非常大的局限性，即僅適用于規(guī)則網(wǎng)格，在非規(guī)則網(wǎng)格或任意網(wǎng)格的背景下公式法無法適用。枚舉的方法能夠適應(yīng)任意網(wǎng)格，并已經(jīng)得到了廣泛的應(yīng)用。利用新策略尋找鄰域的策略也能適應(yīng)任意網(wǎng)格，可直接應(yīng)用至三維結(jié)構(gòu)。

4.2 計(jì)算用時(shí)

選擇5種不同規(guī)模的立方體網(wǎng)格數(shù)量，取10次過濾耗時(shí)的平均值做比較(不含導(dǎo)出單元?單元信息時(shí)間)，計(jì)算環(huán)境為Lenovo Y7000 2020，Intel 10750處理器，16 GB內(nèi)存，在win10系統(tǒng)下，Python版本為3.7，耗時(shí)如表1所示。分析可知，僅在過濾半徑內(nèi)的單元才會對中心單元的靈敏度產(chǎn)生影響，處于過濾半徑之外的單元對其靈敏度影響為0，而鄰域單元相比設(shè)計(jì)空間而言數(shù)量很少，由此可見，利用枚舉直接計(jì)算方法會極大浪費(fèi)計(jì)算資源。

表1 2種策略在不同網(wǎng)格數(shù)量下的計(jì)算耗時(shí)Table 1 Calculation time of the two strategies under different number of elements s

以一橋梁結(jié)構(gòu)優(yōu)化設(shè)計(jì)[17]為例，進(jìn)一步比較2種策略的計(jì)算用時(shí)：結(jié)構(gòu)初始設(shè)計(jì)如圖5所示，為便于行人和車輛通過，留有一層厚度為0.6 m的非設(shè)計(jì)域。由于對稱性僅計(jì)算1/4模型，單元尺寸為0.2 m，劃分為147 500個(gè)單元，彈性模量為20 GPa，泊松比取0.2。體積分?jǐn)?shù)為設(shè)計(jì)域的10%，目標(biāo)為結(jié)構(gòu)剛度最大，刪除率取5%，添加率最大為5%，過濾半徑為0.6 m，收斂容差為0.01，有限元過程采用6核并行計(jì)算。

圖5 設(shè)計(jì)域與非設(shè)計(jì)域、幾何尺寸及邊界條件Fig.5 Design and non-design domains,geometric dimensions and boundary conditions

優(yōu)化結(jié)果如表2所示，運(yùn)用2種策略進(jìn)行優(yōu)化迭代次數(shù)相同，有限元用時(shí)及最終構(gòu)型相似，柔度一致。但對比過濾及有限元用時(shí)可知，雖只在優(yōu)化開始時(shí)組成過濾權(quán)重矩陣，但當(dāng)問題擴(kuò)大到一定規(guī)模之后，這僅僅一次的鄰域搜索會消耗大量的計(jì)算時(shí)間，遠(yuǎn)超有限元迭代用時(shí)，相比而言本文策略具有顯著優(yōu)勢。

表2 2種策略優(yōu)化結(jié)果對比Table 2 Comparison of the optimization results of the two strategies

4.3 設(shè)計(jì)域存在微小孔洞的柔度優(yōu)化問題

以一個(gè)三維算例來對比說明新策略在初始設(shè)計(jì)中有微小非設(shè)計(jì)域孔洞時(shí)的適用性，優(yōu)化背景為考慮在一定體積約束下柔度最小問題。幾何結(jié)構(gòu)如圖6所示。

圖6 優(yōu)化初始構(gòu)型Fig.6 Optimized initial structure

圓盤厚度為10 mm，半徑為55 mm，中間挖有R5 mm的通孔，距離受力位置1 mm處各有一圓心角為60°，厚度為1 mm的細(xì)小縫隙，貫穿圓盤。結(jié)構(gòu)約束及載荷如圖6紅線所示，灰色部分為非設(shè)計(jì)域，中間圓柱面固定，周圍分布有4個(gè)線均布力，力的大小為110 N。有限元離散采用的單元尺寸為1 mm，材料彈性模量為210 GPa，泊松比為0.3?？紤]目標(biāo)體積約束為20%，刪除率取1%，添加率最大為1%，過濾半徑為3 mm，收斂容差為0.001。

圖7顯示了運(yùn)用2種不同的策略對圓盤進(jìn)行優(yōu)化得到的最終拓?fù)浣Y(jié)構(gòu)，新過濾策略能夠很好地避免數(shù)值不穩(wěn)定的現(xiàn)象，得到理想的拓?fù)錁?gòu)型，如圖7(a)。枚舉過濾策略得到的結(jié)構(gòu)大體相同，如圖7(b)，但是在縫隙周圍出現(xiàn)了“孤島”的現(xiàn)象，這是由于縫隙外層單元在載荷附近，位移較大，因此單元靈敏度相對較大，枚舉時(shí)外側(cè)單元的影響會越過縫隙，使得縫隙內(nèi)側(cè)單元加權(quán)后的靈敏度也相對較大，以至在優(yōu)化過程中不會被刪掉。而運(yùn)用新策略搜索鄰域時(shí)，縫隙外邊的單元不在內(nèi)側(cè)單元鄰域范圍內(nèi)，因此不會出現(xiàn)上述情況。

圖7 2種策略最終拓?fù)錁?gòu)型Fig.7 Final topology of the two strategies

圖8給出了2種策略的優(yōu)化迭代曲線，顯然，由于消除了“孤島”現(xiàn)象的存在，新策略得到的最終拓?fù)錁?gòu)型柔度要低于枚舉策略所得到的結(jié)構(gòu)，性能差別達(dá)到4%。

5 結(jié)論

1)提出一種利用節(jié)點(diǎn)單元從屬關(guān)系快速確定單元鄰域的策略，通過建立中心單元——從屬單元——鄰域單元的逐層搜尋模式迅速確定單元鄰域，減少了搜索范圍，其搜索速度與枚舉策略相比有較大優(yōu)勢。

2)新策略的實(shí)現(xiàn)過程簡單方便，能夠適應(yīng)任意網(wǎng)格并且容易拓展至三維結(jié)構(gòu)，結(jié)合此策略的過濾算法能夠有效抑制棋盤格和網(wǎng)格依賴性的現(xiàn)象。

3)對于初始結(jié)構(gòu)設(shè)計(jì)中含有微小非設(shè)計(jì)域孔洞的問題，基于BESO法，使用新策略確定的單元鄰域進(jìn)行后續(xù)優(yōu)化時(shí)的結(jié)果能夠有效消除結(jié)構(gòu)中的“孤島”現(xiàn)象，相對于傳統(tǒng)策略，新策略能夠得到性能更優(yōu)的構(gòu)型。