載荷不確定條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化算法

趙軍鵬 王春潔

(北京航空航天大學(xué) 機(jī)械工程及自動化學(xué)院，北京100191)

結(jié)構(gòu)拓?fù)鋬?yōu)化通常是在給定的載荷工況下尋找材料的最佳分布.然而外界載荷一般存在不確定性，在確定性載荷工況下獲得的最優(yōu)結(jié)構(gòu)形式通常在抵抗載荷擾動時是比較脆弱的，因此有必要研究在載荷不確定條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化問題.現(xiàn)階段通常采用概率方法和非概率方法表示載荷的不確定性[1].概率方法通過指定載荷的概率分布來表示其不確定性，而非概率方法則通過指定載荷允許的變化范圍來表示其不確定性.在結(jié)構(gòu)拓?fù)鋬?yōu)化中，概率方法通常要求結(jié)構(gòu)在不確定載荷作用下柔度的均值與方差最小，而非概率方法則要求結(jié)構(gòu)的最大柔度最小.對于非概率載荷不確定性問題，一般利用橢球表示其不確定性，然后將最大最小問題轉(zhuǎn)化為半定規(guī)劃形式的優(yōu)化模型進(jìn)行求解[2－4];對于概率載荷不確定性問題，求解難度則集中在如何有效計算結(jié)構(gòu)柔度的均值與方差以及它們對設(shè)計變量的靈敏度[5－6].結(jié)構(gòu)柔度均值的計算已經(jīng)出現(xiàn)較為有效的方法[5，7]，然而其方差的計算問題還沒有得到很好的解決.對于桁架結(jié)構(gòu)，針對載荷各分量服從多元正態(tài)分布的情形已給出了結(jié)構(gòu)柔度均值與方差的解析表達(dá)式，并基于此發(fā)展了考慮柔度方差的拓?fù)鋬?yōu)化方法[6]，但是實際工程中載荷的不確定性通常采用載荷大小和作用方向描述，載荷各分量一般不服從多元正態(tài)分布.對于連續(xù)體結(jié)構(gòu)，針對載荷大小不確定的情形同樣發(fā)展了結(jié)構(gòu)拓?fù)鋬?yōu)化方法[8]，但由于未考慮載荷方向不確定性，因此在實際工程中可以應(yīng)用的情形有限.

本文針對載荷不確定條件下以結(jié)構(gòu)柔度均值與標(biāo)準(zhǔn)差的加權(quán)和最小為目標(biāo)的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題，首先利用線彈性結(jié)構(gòu)的位移疊加原理給出結(jié)構(gòu)柔度均值與方差的計算公式，其次在此基礎(chǔ)上給出結(jié)構(gòu)靈敏度分析方法，然后給出載荷不確定條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化算法，最后通過算例說明本文方法的有效性以及考慮載荷不確定性的結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果的穩(wěn)健性.本文的討論基于平面結(jié)構(gòu)問題，但很容易推廣到三維情形.

1 載荷不確定條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化問題

本文采用密度法進(jìn)行連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題的求解，并采用約束最小剛度的插值模型將單元的材料屬性Ee與單元密度ρe聯(lián)系起來[9]:

其中，E0為材料的彈性模量;Emin是為了防止剛度矩陣奇異而賦予無材料單元的很小的剛度值;p≥1為懲罰因子，可以迫使單元密度在優(yōu)化過程中取值趨向0或1，從而獲得構(gòu)型清晰的結(jié)構(gòu).

假設(shè)結(jié)構(gòu)被離散為N個單元，單元e的體積為ve，密度為ρe.則載荷不確定條件下的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題可以表示為

其中，μ(c)與σ(c)分別為結(jié)構(gòu)柔度的均值和標(biāo)準(zhǔn)差;非負(fù)常數(shù) β用于權(quán)衡σ(c)的重要性;Ku(ω)=f(ω)為結(jié)構(gòu)的有限元平衡方程;ω∈Θ表示載荷具有不確定性，ω為概率空間Θ上的任意實現(xiàn);Vmax為許用材料體積;ρ為由各單元密度組成的設(shè)計變量向量.

2 結(jié)構(gòu)柔度均值與方差的計算

假設(shè)結(jié)構(gòu)承受n個不確定載荷，并且每個載荷的大小hi和作用方向與x方向的夾角θi的概率分布均已知.設(shè)f2i－1和f2i分別為單獨在第i個不確定載荷作用點沿x方向和沿y方向施加單位載荷時的整體載荷向量，并且在f2i－1和f2i作用下結(jié)構(gòu)的位移向量分別為u2i－1和u2i，則根據(jù)疊加原理，對于任意的載荷:

相應(yīng)的結(jié)構(gòu)位移為

其中

從而結(jié)構(gòu)的柔度為

為了書寫簡潔，在下面的討論中記m=2n，并引入ξi的二階和四階中心矩:

ξij與ξijkl的解析表達(dá)式通常很難給出，但可以通過蒙特卡羅方法對它們的值進(jìn)行精確估計.

結(jié)構(gòu)柔度的均值為結(jié)構(gòu)在各種可能載荷工況下柔度的期望:

結(jié)構(gòu)柔度的方差可以表示為結(jié)構(gòu)在各種可能載荷工況下柔度平方的期望與結(jié)構(gòu)柔度均值的平方之差:

由此可知，在ξij與ξijkl已知的情況下，只要獲得結(jié)構(gòu)在載荷工況f1，f2，…，fm作用下的位移向量u1，u2，…，um，就可以求得結(jié)構(gòu)柔度的均值與方差，而并不需要對所有可能工況進(jìn)行結(jié)構(gòu)分析.注意到cij為向量fi與uj的內(nèi)積，因此相對于結(jié)構(gòu)有限元分析，其計算量可以忽略不計.

3 結(jié)構(gòu)靈敏度分析

根據(jù)式(9)和式(10)可知，結(jié)構(gòu)柔度均值與標(biāo)準(zhǔn)差對密度ρe的偏導(dǎo)數(shù)分別為

引入變量:

則目標(biāo)函數(shù)對ρe的偏導(dǎo)數(shù)為

利用拉格朗日乘子法容易證明:

體積約束對于設(shè)計變量的靈敏度為

4 載荷不確定條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化算法

綜合以上內(nèi)容可以給出載荷不確定條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化算法，其流程如圖1所示.首先利用蒙特卡羅方法計算ξij與ξijkl;然后離散設(shè)計域并初始化設(shè)計變量;在此基礎(chǔ)上對結(jié)構(gòu)進(jìn)行多載荷工況有限元分析;根據(jù)有限元分析結(jié)果計算結(jié)構(gòu)柔度均值與方差并進(jìn)行靈敏度分析;為了消除結(jié)構(gòu)拓?fù)鋬?yōu)化過程中經(jīng)常出現(xiàn)的棋盤格等現(xiàn)象，采用靈敏度過濾器[10]對靈敏度信息進(jìn)行過濾，然后利用移動漸近線[11]方法更新設(shè)計變量并判斷優(yōu)化過程是否收斂.若收斂，則迭代結(jié)束;否則重新進(jìn)行有限元分析等直至優(yōu)化過程收斂.

從圖1中可以看出，在載荷不確定條件下進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化，如果承受n個不確定載荷，則每次迭代需要對2n個工況進(jìn)行有限元分析，而不考慮載荷不確定性時則僅需要進(jìn)行單個工況有限元分析.因此載荷不確定條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化算法計算量要大于確定性條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化算法.

圖1 載荷不確定條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化流程

5 數(shù)值算例

本節(jié)以兩個算例來說明載荷不確定條件下的最優(yōu)結(jié)構(gòu)與確定條件下的最優(yōu)結(jié)構(gòu)不同，并且前者要比后者能夠更好地抵抗載荷的擾動.兩個算例均采用單位雙線性正方形單元對設(shè)計域進(jìn)行離散.蒙特卡羅方法計算 ξij與 ξijkl時樣本數(shù)目為1000000.實體與空單元的彈性模量分別為E0=1，Emin=10－9，材料的泊松比 μ 均為 0.3，懲罰因子p=3.優(yōu)化結(jié)束的準(zhǔn)則為每個單元的偽密度變化量均不超出0.01.

5.1 簡單立柱結(jié)構(gòu)的拓?fù)鋬?yōu)化

如圖2a所示，問題的設(shè)計域為100×100的矩形區(qū)域，底部邊界固定，結(jié)構(gòu)頂端中心承受不確定載荷F.該載荷的大小為10，作用方向與水平方向的夾角服從均值為－π/2(豎直向下)，標(biāo)準(zhǔn)差為π/48的正態(tài)分布.材料的許用體積為設(shè)計域的20%.

結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果如圖2b～圖2d所示.載荷確定條件下材料全部聚集在載荷作用線附近，形成立柱狀結(jié)構(gòu)，這與工程直覺相符，但該結(jié)構(gòu)的水平方向承載能力比較弱.考慮載荷不確定性得到的結(jié)構(gòu)均類似于兩個頂端相交桿件組成的結(jié)構(gòu)，這樣的結(jié)構(gòu)雖然豎直方向承載能力有所下降，但是水平方向承載能力卻得到很大提高，從而增強(qiáng)了結(jié)構(gòu)的穩(wěn)健性.

圖2 拓?fù)鋬?yōu)化設(shè)計區(qū)域及不同條件下的優(yōu)化結(jié)果

結(jié)構(gòu)拓?fù)鋬?yōu)化過程收斂曲線如圖3所示，優(yōu)化結(jié)果的目標(biāo)函數(shù)值如表1所示.

圖3 簡單立柱結(jié)構(gòu)拓?fù)鋬?yōu)化收斂曲線

表1 簡單立柱結(jié)構(gòu)優(yōu)化結(jié)果的目標(biāo)函數(shù)值

可以看出，確定性載荷條件下得到的優(yōu)化結(jié)果在考慮載荷不確定性時的目標(biāo)函數(shù)值要高于不確定載荷條件下得到的優(yōu)化結(jié)果的目標(biāo)函數(shù)值，這證明了本文算法的有效性.

載荷不確定性條件下的最優(yōu)結(jié)構(gòu)與優(yōu)化目標(biāo)中柔度標(biāo)準(zhǔn)差的權(quán)重有關(guān)，權(quán)重增大，兩個桿件之間的夾角也增大，結(jié)構(gòu)柔度由于載荷不確定性而引起的波動會減小，但是柔度均值可能會加大.鑒于篇幅的限制，本文不對β如何合理取值的問題進(jìn)行深入討論.

5.2 米歇爾結(jié)構(gòu)的拓?fù)鋬?yōu)化

如圖4a所示，問題的設(shè)計域為240×120的矩形區(qū)域，底部兩端均簡支，并且承受3個不確定載荷 P1，P2，P3.3 個載荷的大小 hi和作用方向與水平方向夾角θi均服從正態(tài)分布，它們的均值μhi，μθi和標(biāo)準(zhǔn)差σhi，σθi如表 2 所示.材料的許用體積為設(shè)計域的30%.

表2 載荷大小和方向的均值與標(biāo)準(zhǔn)差

取目標(biāo)函數(shù)中的權(quán)重β=1，結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果如圖4所示.由于結(jié)構(gòu)的設(shè)計域、邊界條件以及作用載荷均左右對稱，因此優(yōu)化后的結(jié)構(gòu)也均左右對稱.其中只考慮載荷大小不確定性時獲得的結(jié)構(gòu)與載荷確定條件下獲得的結(jié)構(gòu)構(gòu)型基本一致，同時考慮載荷大小與方向不確定性時獲得的結(jié)構(gòu)與只考慮載荷方向不確定性時獲得的結(jié)構(gòu)構(gòu)型基本一致.對本例而言，載荷方向不確定性比大小不確定性對結(jié)構(gòu)的影響更大.考慮載荷方向不確定性得到的結(jié)構(gòu)中間兩根較粗的桿件直徑有所減小，兩根較細(xì)的桿件與水平方向的夾角也變小，而其下方呈現(xiàn)出連接左右兩個端點以及3個載荷作用點的桿件.這樣的結(jié)構(gòu)雖然豎直方向承載能力有所下降，但是卻更加穩(wěn)健.這是由于載荷作用方向的不確定性，其可能含有水平方向分量，而水平方向桿件具有較好的承受水平方向載荷的能力，其存在使得結(jié)構(gòu)可以更好地抵抗載荷的不確定性，從而使得結(jié)構(gòu)更加穩(wěn)健.

圖4 米歇爾結(jié)構(gòu)拓?fù)鋬?yōu)化問題及優(yōu)化結(jié)果

結(jié)構(gòu)拓?fù)鋬?yōu)化過程收斂曲線如圖5所示，優(yōu)化結(jié)果的目標(biāo)函數(shù)值如表3所示.

表3 米歇爾結(jié)構(gòu)優(yōu)化結(jié)果的目標(biāo)函數(shù)值

可以看出，確定性載荷條件下得到的優(yōu)化結(jié)果在考慮載荷不確定性時的目標(biāo)函數(shù)值一般高于不確定載荷條件下得到的優(yōu)化結(jié)果的目標(biāo)函數(shù)值，這證明了本文算法的有效性.注意到載荷大小不確定條件下優(yōu)化結(jié)果的目標(biāo)函數(shù)值與確定條件下優(yōu)化結(jié)果的目標(biāo)函數(shù)值相差不超過0.5%，二者的優(yōu)化結(jié)果也幾乎相同.

圖5 米歇爾結(jié)構(gòu)拓?fù)鋬?yōu)化收斂曲線

6 結(jié)論

1)與確定性載荷作用下的結(jié)構(gòu)拓?fù)鋬?yōu)化方法相比，本文方法得到的優(yōu)化結(jié)果在載荷不確定條件下具有更小的柔度均值和方差，因此更加穩(wěn)健;

2)優(yōu)化目標(biāo)中柔度均值與標(biāo)準(zhǔn)差的相對權(quán)重對結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果有較大影響;

3)對于本文中的算例，相對于載荷大小不確定性，方向不確定性對優(yōu)化結(jié)果有更顯著的影響.

致謝 衷心感謝瑞典皇家理工學(xué)院的Krister Svanberg教授提供移動漸近線法的代碼.

References)

[1]Calafiore G C，Dabbene F.Optimization under uncertainty with applications to design of truss structures[J].Structural and Multidisciplinary Optimization，2008，35(3):189 －200

[2]Ben-Tal A，Nemirovski A.Robust truss topology design via semidefinite programming[J].SIAM Journal on Optimization，1997，7(4):991－1016

[3]李東澤，于登云，馬興瑞.不確定載荷下的桁架結(jié)構(gòu)拓?fù)鋬?yōu)化[J].北京航空航天大學(xué)學(xué)報，2009，35(10):1170 －1173

Li Dongze，Yu Dengyun，Ma Xingrui.Truss topology optimization with uncertain loading scenarios[J].Journal of Beijing University of Aeronautics and Astronautics，2009，35(10):1170 － 1173(in Chinese)

[4]De Gournay F，Allaire G，Jouve F.Shape and topology optimization of the robust compliance via the level set method [J].Esaim-Control Optimisation and Calculus of Variations，2008，24(1):43－70

[5]Alvarez F，Carrasco M.Minimization of the expected compliance as an alternative approach to multiload truss optimization[J].Structural and Multidisciplinary Optimization，2005，29(6):470－476

[6]Carrasco M，Ivorra B，Ramos A.A variance-expected compliance model for structural optimization[J].Journal of Optimization Theory and Applications，2012，152(1):136 －151

[7]Dunning P D，Kim H A，Mullineux G.Introducing loading uncertainty in topology optimization[J].AIAA Journal，2011，49(4):760－768

[8]Dunning P D，Kim H A.Robust topology optimization:minimization of expected and variance of compliance[J].AIAA Journal，2013，50(11):2656 －2664

[9]Andreassen E，Clausen A，Schevenels M，et al.Efficient topology optimization in MATLAB using 88 lines of code[J].Structural and Multidisciplinary Optimization，2011，43(1):1 －16

[10]Sigmund O.On the design of compliant mechanisms using topology optimization[J].Mechanics of Structures and Machines:An International Journal，1997，25(4):493 － 524

[11]Svanberg K.The method of moving asymptotes:a method for structural optimization[J].International Journal for Numerical Methods in Engineering，1987，24(2):359 －373