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

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

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

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

1 BESO法

以實(shí)際常見(jiàn)的受到體積約束的剛度最大問(wèn)題為例，建立基于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ì)變量，表示單元的相對(duì)密度，實(shí)體單元取值為1，為了保證整體剛度矩陣的非奇異性，空單元取較小值xmin，一般為0.001。

引入材料插值方案：

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

其中：K0i為實(shí)體狀態(tài)下的單元?jiǎng)偠染仃嚒?/p>

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

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

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

靈敏度過(guò)濾公式為：

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

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

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

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

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

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

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

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

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

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

3 新策略靈敏度過(guò)濾流程

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

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

4 過(guò)濾策略的比較

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

4.1 網(wǎng)格適用性

依據(jù)規(guī)則編號(hào)確定單元鄰域具有非常大的局限性，即僅適用于規(guī)則網(wǎng)格，在非規(guī)則網(wǎng)格或任意網(wǎng)格的背景下公式法無(wú)法適用。枚舉的方法能夠適應(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次過(guò)濾耗時(shí)的平均值做比較(不含導(dǎo)出單元?單元信息時(shí)間)，計(jì)算環(huán)境為L(zhǎng)enovo Y7000 2020，Intel 10750處理器，16 GB內(nèi)存，在win10系統(tǒng)下，Python版本為3.7，耗時(shí)如表1所示。分析可知，僅在過(guò)濾半徑內(nèi)的單元才會(huì)對(duì)中心單元的靈敏度產(chǎn)生影響，處于過(guò)濾半徑之外的單元對(duì)其靈敏度影響為0，而鄰域單元相比設(shè)計(jì)空間而言數(shù)量很少，由此可見(jiàn)，利用枚舉直接計(jì)算方法會(huì)極大浪費(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所示，為便于行人和車(chē)輛通過(guò)，留有一層厚度為0.6 m的非設(shè)計(jì)域。由于對(duì)稱性僅計(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%，過(guò)濾半徑為0.6 m，收斂容差為0.01，有限元過(guò)程采用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)型相似，柔度一致。但對(duì)比過(guò)濾及有限元用時(shí)可知，雖只在優(yōu)化開(kāi)始時(shí)組成過(guò)濾權(quán)重矩陣，但當(dāng)問(wèn)題擴(kuò)大到一定規(guī)模之后，這僅僅一次的鄰域搜索會(huì)消耗大量的計(jì)算時(shí)間，遠(yuǎn)超有限元迭代用時(shí)，相比而言本文策略具有顯著優(yōu)勢(shì)。

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

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

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

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

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

圖7顯示了運(yùn)用2種不同的策略對(duì)圓盤(pán)進(jìn)行優(yōu)化得到的最終拓?fù)浣Y(jié)構(gòu)，新過(guò)濾策略能夠很好地避免數(shù)值不穩(wěn)定的現(xiàn)象，得到理想的拓?fù)錁?gòu)型，如圖7(a)。枚舉過(guò)濾策略得到的結(jié)構(gòu)大體相同，如圖7(b)，但是在縫隙周?chē)霈F(xiàn)了“孤島”的現(xiàn)象，這是由于縫隙外層單元在載荷附近，位移較大，因此單元靈敏度相對(duì)較大，枚舉時(shí)外側(cè)單元的影響會(huì)越過(guò)縫隙，使得縫隙內(nèi)側(cè)單元加權(quán)后的靈敏度也相對(duì)較大，以至在優(yōu)化過(guò)程中不會(huì)被刪掉。而運(yùn)用新策略搜索鄰域時(shí)，縫隙外邊的單元不在內(nèi)側(cè)單元鄰域范圍內(nèi)，因此不會(huì)出現(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)系快速確定單元鄰域的策略，通過(guò)建立中心單元——從屬單元——鄰域單元的逐層搜尋模式迅速確定單元鄰域，減少了搜索范圍，其搜索速度與枚舉策略相比有較大優(yōu)勢(shì)。

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

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