設(shè)計(jì)導(dǎo)向下的結(jié)構(gòu)拓?fù)鋬?yōu)化研究

周嘉侃，夏利娟，黃開(kāi)藝

（上海交通大學(xué)a.海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心；b.電力傳輸與功率變換控制教育部重點(diǎn)實(shí)驗(yàn)室，上海 200240）

0 引 言

結(jié)構(gòu)拓?fù)鋬?yōu)化是工程結(jié)構(gòu)設(shè)計(jì)的重要組成部分，具有效率高、靈活性高等特點(diǎn)。1988年，Bendsoe與Kikuchi 等[1]以引入微結(jié)構(gòu)單元的理念首次提出了均勻化拓?fù)鋬?yōu)化方法。在均勻化方法的基礎(chǔ)上，Mlejnek等[2]提出了變密度方法，將設(shè)計(jì)變量轉(zhuǎn)變?yōu)槿藶槎x的偽密度參數(shù)，并通過(guò)偽密度的取值決定單元的刪除或保留。變密度法常用的插值函數(shù)模型有SIMP 模型[3]與RAMP 模型[4]，其中SIMP 模型基于各向同性微結(jié)構(gòu)，以冪函數(shù)的形式建立偽密度與體積、剛度的關(guān)系。在拓?fù)鋬?yōu)化的研究過(guò)程中，一些學(xué)者針對(duì)變密度法的材料插值函數(shù)進(jìn)行改進(jìn)[5-8]，并對(duì)優(yōu)化過(guò)程中存在的棋盤(pán)格效應(yīng)、灰度單元等問(wèn)題給出了解決方案[9-10]。結(jié)構(gòu)拓?fù)鋬?yōu)化在船舶設(shè)計(jì)中具有廣泛的應(yīng)用，劉宏亮等[11]使用生長(zhǎng)進(jìn)化拓?fù)鋬?yōu)化算法進(jìn)行了船舶中剖面橫撐結(jié)構(gòu)的優(yōu)化；程遠(yuǎn)勝等[12]對(duì)肘板的拓?fù)湫问竭M(jìn)行了優(yōu)化，降低了應(yīng)力集中現(xiàn)象；張會(huì)新等[13]對(duì)船舶上層建筑板架進(jìn)行拓?fù)鋬?yōu)化，提高了板架的固有頻率。

在實(shí)際生產(chǎn)中，對(duì)于結(jié)構(gòu)的拓?fù)錁?gòu)型存在具體的設(shè)計(jì)需求，如船舶結(jié)構(gòu)中因維護(hù)需要而設(shè)立的人孔，因裝載與搬運(yùn)便利而設(shè)立的大開(kāi)孔，客機(jī)機(jī)艙地板梁預(yù)留鋪設(shè)電纜與管道的空洞，因阻隔需要而設(shè)置的外突結(jié)構(gòu)等。這些設(shè)計(jì)需求對(duì)結(jié)構(gòu)力學(xué)屬性的影響不可忽略，丁運(yùn)來(lái)等[14]對(duì)船體結(jié)構(gòu)中不同的開(kāi)孔結(jié)構(gòu)進(jìn)行了研究，計(jì)算了孔邊緣的應(yīng)力分布。而拓?fù)鋬?yōu)化的最終構(gòu)型因不可控制而難以應(yīng)用于實(shí)踐，目前關(guān)于改善拓?fù)錁?gòu)型的相關(guān)研究主要集中于解決工程制造性、棋盤(pán)格效應(yīng)等問(wèn)題，如李伯豪等[15]提出了一種基于參數(shù)化水平集的后處理方法，提高了拓?fù)錁?gòu)型的可制造性；Boutdin 與Borrvall等[16-17]提出了密度過(guò)濾法，通過(guò)修正偽密度改善拓?fù)錁?gòu)型，處理棋盤(pán)格效應(yīng)。盡管結(jié)構(gòu)的設(shè)計(jì)需求十分重要，但以既定設(shè)計(jì)需求構(gòu)型為導(dǎo)向的拓?fù)鋬?yōu)化方法仍有待進(jìn)一步研究。

基于此，本文通過(guò)構(gòu)造基于單元靈敏度的修正函數(shù)，改進(jìn)SIMP變密度法，建立設(shè)計(jì)導(dǎo)向下的結(jié)構(gòu)拓?fù)鋬?yōu)化模型。之后通過(guò)Python語(yǔ)言編程加以實(shí)現(xiàn)，使用短懸臂梁經(jīng)典算例對(duì)模型進(jìn)行可行性分析，最后對(duì)某油船貨油艙板架結(jié)構(gòu)進(jìn)行設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化。

1 SIMP法拓?fù)鋬?yōu)化模型

1.1 優(yōu)化算法思想

本文主要選取以體積為約束條件，結(jié)構(gòu)剛度最大化的拓?fù)鋬?yōu)化問(wèn)題。首先使用有限元軟件對(duì)連續(xù)體進(jìn)行單元離散，賦予單元偽密度。選取最優(yōu)準(zhǔn)則法作為優(yōu)化算法，通過(guò)庫(kù)恩塔克條件建立優(yōu)化準(zhǔn)則并生成迭代格式，對(duì)單元偽密度進(jìn)行迭代。迭代完成后，保留偽密度較大的單元，剔除偽密度較小的單元。

1.2 靈敏度計(jì)算與迭代公式

本文選取的SIMP拓?fù)鋬?yōu)化模型可表述為

式中，V*為給定的體積約束，U為位移向量，K與Ki分別為總體剛度矩陣與單元?jiǎng)偠染仃?，Ki0為初始單元?jiǎng)偠染仃嚕瑃為設(shè)計(jì)變量向量，p為權(quán)系數(shù)。根據(jù)文獻(xiàn)[18]，構(gòu)造拉格朗日函數(shù)L(t)：

式中，β、λj、χk為拉格朗日乘子。由于迭代作用的偽密度介于上下界之間，可得λj=χk= 0，根據(jù)庫(kù)恩塔克條件，

式中目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的靈敏度可通過(guò)對(duì)有限元方程K(t)U(t)=F與式（1）中目標(biāo)函數(shù)求導(dǎo)，并聯(lián)立化簡(jiǎn)得到

結(jié)合式（3）庫(kù)恩塔克條件，可構(gòu)造迭代因子：

式中拉格朗日乘子β可根據(jù)文獻(xiàn)[19]的二分法計(jì)算，迭代公式如下所示：

式中，tmin為避免剛度矩陣奇異而設(shè)置的最小偽密度。為保持迭代過(guò)程的穩(wěn)定，引入移動(dòng)極限常數(shù)δ和阻尼因子η，SIMP法迭代公式為

式中迭代參數(shù)的選取如表1所示：

表1 SIMP法迭代參數(shù)Tab.1 Iterative parameters of SIMP

1.3 靈敏度過(guò)濾與收斂準(zhǔn)則

采用靈敏度過(guò)濾法處理棋盤(pán)格效應(yīng)，對(duì)式（4）的靈敏度進(jìn)行過(guò)濾：

式中：?c(t)*/?ti為修正后單元i的靈敏度；?c(t)/?tf為修正前單元f的靈敏度；Hf=rmin-d(i,f)，rmin為過(guò)濾半徑，取為1.5倍的單元最小直徑，d(i,f)表示單元i中心至單元f中心的距離；集合B表示所有與單元i距離小于rmin的單元集合，即{f|1 ≤f≤n,Hf＞0} 。

采用兩次迭代設(shè)計(jì)變量最大分量的相對(duì)誤差作為收斂準(zhǔn)則，同時(shí)加入最大迭代次數(shù)以保證優(yōu)化迭代的終止：

式中，ε取0.01，最大迭代次數(shù)取50，當(dāng)?shù)K止時(shí)，以0.5作為單元增刪的臨界密度。

2 設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化模型

2.1 設(shè)計(jì)導(dǎo)向算法思想

對(duì)拓?fù)錁?gòu)型的需求可轉(zhuǎn)化為控制設(shè)計(jì)區(qū)域內(nèi)某些子域中單元的存在性，可通過(guò)修改迭代過(guò)程中對(duì)應(yīng)單元的偽密度實(shí)現(xiàn)。如需子域A中的單元在拓?fù)鋬?yōu)化的結(jié)果中存在，則在每一次迭代中，不斷增加A中單元的偽密度，提升單元被保留的概率，反之亦然，最終達(dá)成拓?fù)鋬?yōu)化迭代與設(shè)計(jì)導(dǎo)向之間的平衡點(diǎn)。

2.2 設(shè)計(jì)導(dǎo)向修正函數(shù)

本文使用設(shè)計(jì)導(dǎo)向修正函數(shù)Di控制迭代中單元密度的調(diào)整，使拓?fù)涞蛟O(shè)計(jì)需求構(gòu)型轉(zhuǎn)變，修正方式如下：

式中，t為調(diào)整后的單元密度，tD為修正幅值，用于控制導(dǎo)向幅度的大小。設(shè)計(jì)需求導(dǎo)向修正函數(shù)Di需滿足式（11），保證調(diào)整前后結(jié)構(gòu)總體積的恒定：

當(dāng)單元被均勻離散時(shí)，不同單元之間的體積差異可被忽略，可近似為= 0。

設(shè)需要調(diào)整的子域?yàn)镸，作如下定義：若設(shè)計(jì)構(gòu)型要求子域M中的單元存在，則稱子域M為存在域，反之則稱其為剔除域。表2列出了本文構(gòu)造的設(shè)計(jì)導(dǎo)向修正函數(shù)的表達(dá)式。

表2 設(shè)計(jì)導(dǎo)向修正函數(shù)Tab.2 Design-oriented correction functions

表2 中，常數(shù)P用于滿足式（11）；Si為單元i的靈敏度；Sinc_max、Sinc_min、Sexc_max、Sexc_min分別為子域M內(nèi)外的單元靈敏度最大與最小值，用于靈敏度歸一化處理；正數(shù)q為形狀因子，用于控制修正函數(shù)的形狀，q越大表示修正函數(shù)對(duì)靈敏度的變化越敏感。

修正函數(shù)根據(jù)單元的靈敏度確定單次迭代中對(duì)單元偽密度修正的量。以剔除域M為例，為降低子域M中單元的偽密度，同時(shí)應(yīng)盡可能減少因改變密度而導(dǎo)致的目標(biāo)函數(shù)增加，修正函數(shù)采用如下策略進(jìn)行單元偽密度調(diào)整：

（1）對(duì)于剔除域M內(nèi)的單元，降低處于存在狀態(tài)的單元密度，且單元靈敏度越低，調(diào)整幅度越大，反之亦然。

（2）對(duì)于剔除域M外的單元，增加處于刪除狀態(tài)的單元密度，且單元靈敏度越高，調(diào)整幅度越大，反之亦然。

2.3 設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化流程及程序?qū)崿F(xiàn)

通過(guò)上節(jié)構(gòu)造的設(shè)計(jì)導(dǎo)向修正函數(shù)，在每一次迭代中，單元的偽密度將會(huì)被修正。具體的優(yōu)化流程如圖1所示，優(yōu)化求解步驟如下：

圖1 設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化流程Fig.1 Design-oriented topology optimization process

（1）建立有限元模型，定義設(shè)計(jì)區(qū)域、設(shè)計(jì)導(dǎo)向的子域等，施加邊界條件與載荷。

（2）進(jìn)行有限元分析計(jì)算，提取結(jié)果文件中各單元的剛度矩陣與位移向量。

（3）計(jì)算設(shè)計(jì)變量靈敏度，進(jìn)行靈敏度過(guò)濾，計(jì)算迭代因子、拉格朗日乘子，代入迭代公式計(jì)算得到新的設(shè)計(jì)變量值。

2.5.2.2 農(nóng)業(yè)防治首先進(jìn)行合理輪作，不能和瓜類(lèi)、向日葵、茄科等作物輪作，可以和禾本科、大豆、甜菜等輪作。其次是進(jìn)行深耕秋翻，將種子翻入20cm以下土壤深處，可以減輕食葵列當(dāng)種子的萌發(fā)。最后就是及時(shí)鏟除列當(dāng)，在列當(dāng)開(kāi)花結(jié)籽前，結(jié)合除草將列當(dāng)鏟除。

（4）統(tǒng)計(jì)子域內(nèi)外的設(shè)計(jì)變量靈敏度最值，計(jì)算設(shè)計(jì)導(dǎo)向修正函數(shù)，再次更新設(shè)計(jì)變量。

（5）判別是否滿足收斂條件，如未滿足，則修改模型，返回步驟（2）；如滿足，則迭代結(jié)束，根據(jù)單元?jiǎng)h除的臨界密度獲得對(duì)應(yīng)拓?fù)錁?gòu)型。

本文使用Hypermesh建立有限元模型，并輸入Nastran進(jìn)行有限元計(jì)算，從計(jì)算結(jié)果文件中提取各個(gè)單元的剛度矩陣與位移向量，通過(guò)Python編程實(shí)現(xiàn)設(shè)計(jì)變量的更新，再根據(jù)所得的結(jié)果修改有限元模型，實(shí)現(xiàn)設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化迭代過(guò)程。

3 算例及分析

3.1 短懸臂梁算例

如圖2所示，懸臂梁左端固支約束，右邊緣中點(diǎn)受豎直向下方向的集中載荷作用，梁的長(zhǎng)度為200 mm，寬度為100 mm，集中載荷大小為1000 N。劃分網(wǎng)格為40×20，材料彈性模量為206 000 MPa，泊松比為0.3，優(yōu)化體積分?jǐn)?shù)為0.5。

圖2 短懸臂梁Fig.2 Short cantilever beam

圖3 為無(wú)設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化結(jié)果，優(yōu)化后柔順度為333.745 N·m，將優(yōu)化構(gòu)型的桁架根部設(shè)置為剔除域。選取“導(dǎo)向比”評(píng)價(jià)設(shè)計(jì)導(dǎo)向的效果，其定義為設(shè)計(jì)子域內(nèi)達(dá)成導(dǎo)向的單元數(shù)與子域單元總數(shù)之比，如圖3的導(dǎo)向比為3/16；選取“柔順比”評(píng)價(jià)優(yōu)化后的力學(xué)性能，定義為設(shè)計(jì)導(dǎo)向介入與不介入的柔順度之比。導(dǎo)向比越高，柔順比越低，表示設(shè)計(jì)導(dǎo)向效果越好。選取不同參數(shù)對(duì)短懸臂梁進(jìn)行設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化，得到表3與表4所示的結(jié)果。

圖3 無(wú)設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化結(jié)果Fig.3 Topology optimization result without design guidance

表3 不同設(shè)計(jì)導(dǎo)向參數(shù)下短懸臂梁優(yōu)化結(jié)果構(gòu)型Tab.3 Optimal configurations of short cantilever beam with different parameters

表4 不同設(shè)計(jì)導(dǎo)向參數(shù)下短懸臂梁優(yōu)化數(shù)據(jù)Tab.4 Optimization data of short cantilever beam with different parameters

3.2 油船板架算例

圖4（a）為一艘320 000 DWT 油船中部貨油艙有限元艙段計(jì)算模型。如圖4（b）所示，優(yōu)化區(qū)域位于橫艙壁相鄰肋位，包含船底肋板與舭肘板。由于肋板靠近橫艙壁，因此僅考慮圖5 所示的4 種橫向極限載荷工況，各工況的權(quán)重系數(shù)根據(jù)單工況下結(jié)構(gòu)應(yīng)變能大小確定。圖6（a）為無(wú)設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化構(gòu)型，在橫向載荷工況下，船底肋板由于承載來(lái)自海水與貨物的壓力，應(yīng)力分布集中，單元傾向于被保留。在船底肋板處設(shè)置剔除域，使用設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化算法進(jìn)行優(yōu)化。根據(jù)短懸臂梁算例結(jié)果，形狀因子選取ln6，修正幅值選取0.4，得到優(yōu)化后的構(gòu)型如圖6（b）所示，優(yōu)化數(shù)據(jù)如表5所示。

圖4 油船有限元艙段模型Fig.4 Finite element model of oil tanker cabin

圖5 油船艙段載荷工況Fig.5 Load conditions of oil tanker cabin

表5 油船板架優(yōu)化數(shù)據(jù)Tab.5 Optimization data of oil tanker frames

圖6（a）中，無(wú)設(shè)計(jì)導(dǎo)向的拓?fù)鋬?yōu)化構(gòu)型中船底肋板采用全封閉式，未考慮因檢修或維護(hù)需要而設(shè)置的人孔。圖6（b）中，隨著設(shè)計(jì)導(dǎo)向的介入，船底肋板出現(xiàn)開(kāi)孔，而舭肘板處開(kāi)孔減小，可說(shuō)明在船底肋板開(kāi)孔的同時(shí)需要加固舭肘板。算例導(dǎo)向后的柔順比小于1.2，表明設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化模型能夠結(jié)合開(kāi)孔的設(shè)計(jì)需求，得到更加合理的拓?fù)錁?gòu)型以指導(dǎo)設(shè)計(jì)，進(jìn)一步體現(xiàn)模型的可行性與穩(wěn)定性。通過(guò)對(duì)比不同拓?fù)鋬?yōu)化結(jié)果，綜合考慮柔順度、結(jié)構(gòu)重量、導(dǎo)向比等指標(biāo)，進(jìn)行重新設(shè)計(jì)，得到新的結(jié)構(gòu)拓?fù)湫问饺鐖D6（c）所示。

圖6 油船板架優(yōu)化結(jié)果構(gòu)型Fig.6 Optimal configurations of oil tanker frames

4 結(jié) 論

本文通過(guò)構(gòu)造修正函數(shù)，將其引入SIMP 變密度拓?fù)鋬?yōu)化算法，建立了設(shè)計(jì)導(dǎo)向下的結(jié)構(gòu)拓?fù)鋬?yōu)化模型。采用以單元靈敏度為依據(jù)的修正函數(shù)，盡可能保留結(jié)構(gòu)的力學(xué)性能，實(shí)現(xiàn)了拓?fù)錁?gòu)型向設(shè)計(jì)需求方向引導(dǎo)，并通過(guò)算例驗(yàn)證了模型的可行性和穩(wěn)定性。數(shù)值計(jì)算可以得到如下結(jié)論：

（1）在拓?fù)鋬?yōu)化程序中引入設(shè)計(jì)導(dǎo)向流程能夠在維持一定柔順比的情況下，獲得較高導(dǎo)向比的拓?fù)錁?gòu)型，因此采用單元偽密度修正方法是可行的。

（2）模型參數(shù)的選取十分重要，過(guò)大的形狀因子會(huì)導(dǎo)致柔順比上升。對(duì)于本文算例，修正幅值選取介于0.3至0.4，形狀因子選取介于ln2至ln6可獲得穩(wěn)定、導(dǎo)向合理的構(gòu)型。

（3）設(shè)計(jì)導(dǎo)向下的拓?fù)鋬?yōu)化模型在求解過(guò)程中能夠自發(fā)地進(jìn)行材料優(yōu)化配置，可為實(shí)際工程結(jié)構(gòu)的優(yōu)化設(shè)計(jì)提供參考方案。