許小奎 郭寶峰 金 淼

燕山大學(xué)先進(jìn)鍛壓成形技術(shù)與科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，秦皇島，066004

?

基于密度體積插值方法的結(jié)構(gòu)拓?fù)鋬?yōu)化

許小奎 郭寶峰 金 淼

燕山大學(xué)先進(jìn)鍛壓成形技術(shù)與科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，秦皇島，066004

針對變密度法結(jié)構(gòu)拓?fù)鋬?yōu)化中灰度單元的控制問題，提出了一種密度體積插值方法。該方法構(gòu)造的剛度與相對密度之間的線性關(guān)系保證了迭代中單元剛度變化的穩(wěn)定性；構(gòu)造體積與相對密度之間的懲罰關(guān)系以實(shí)現(xiàn)懲罰中間密度，同時(shí)更有利于在中間密度向兩端逼近的同時(shí)降低灰度單元的數(shù)量。應(yīng)用該插值方法對位移約束體積最小化問題進(jìn)行求解，所得結(jié)果顯示，增大懲罰程度，優(yōu)化過程浮動較小，算法相對穩(wěn)定。與應(yīng)用SIMP方法和RAMP方法的優(yōu)化結(jié)果相比較可知，在求解同一結(jié)構(gòu)拓?fù)鋬?yōu)化問題時(shí)，采用該方法且增大懲罰程度后的優(yōu)化結(jié)果中灰度單元減少明顯。

插值方法；變密度法；位移約束；拓?fù)鋬?yōu)化

0 引言

自1988年BENDSOE等[1]提出用于求解連續(xù)體拓?fù)鋬?yōu)化的均勻化方法以來，連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化一直是結(jié)構(gòu)優(yōu)化領(lǐng)域研究的熱點(diǎn)。變密度法[2-5]是在均勻化方法的基礎(chǔ)上發(fā)展形成的一種算法，該方法由于建模簡單、易于求解等特點(diǎn)而得到了廣泛應(yīng)用[6-7]。

變密度法中常用的密度插值方法是密度剛度插值方法，其原理是通過構(gòu)造相對密度與剛度之間的懲罰關(guān)系來達(dá)到消除中間密度的目的。典型的密度剛度插值方法有SIMP (solid isotropic microstructures with penalization)方法[2-3]和RAMP (rational approximation of material properties)方法[4]。由于優(yōu)化計(jì)算是建立在結(jié)構(gòu)有限元分析基礎(chǔ)之上的，密度剛度懲罰函數(shù)的作用主要體現(xiàn)在結(jié)構(gòu)有限元分析上，即優(yōu)化計(jì)算是基于懲罰剛度后的有限元分析結(jié)果的，優(yōu)化過程中相對密度變化時(shí)，單元剛度變化較大，導(dǎo)致算法穩(wěn)定性差且易使優(yōu)化落入局部最優(yōu)解，因此必須控制密度剛度函數(shù)的懲罰程度。對于SIMP方法，懲罰因子取值范圍一般為3～5[8]。然而限定懲罰程度帶來的后果是，中間密度不能夠得到更大程度的懲罰，導(dǎo)致優(yōu)化結(jié)果會存在相當(dāng)一部分灰度單元。

針對上述問題，本文提出了一種密度體積插值方法，并采用該方法求解了位移約束條件下的結(jié)構(gòu)拓?fù)鋬?yōu)化問題。

1 密度體積插值方法

密度體積插值方法中為實(shí)現(xiàn)對中間密度的懲罰，構(gòu)造單元體積與相對密度之間的關(guān)系為

Vi=Vi0ω(ρi)

(1)

其中，ρi為相對密度；Vi0和Vi分別為實(shí)體單元和插值后單元的體積；ω(*)為體積懲罰函數(shù)，對中間密度單元的體積向上進(jìn)行懲罰，映射到最小體積的優(yōu)化中，即為對目標(biāo)函數(shù)進(jìn)行了懲罰。

懲罰函數(shù)ω可為指數(shù)型懲罰函數(shù)：

(2)

亦可為有理型懲罰函數(shù)：

(3)

其中，p和q分別為兩種形式懲罰函數(shù)的懲罰因子。如圖1所示，ω1中懲罰因子p取值越小,懲罰程度越大，ω2中懲罰因子q取值越大,懲罰程度越大。懲罰程度越大，越有利于中間密度向兩端逼近及降低灰度單元的數(shù)量。

(a)ω1

(b)ω2圖1 ω函數(shù)懲罰曲線圖Fig.1 The curves of penalty function ω

密度體積插值方法中，構(gòu)造彈性模量與相對密度之間為線性函數(shù)：

Ei=ρiE0

(4)

式中，E0和Ei分別為實(shí)體單元和插值后單元的彈性模量。

這種線性關(guān)系即意味著單元剛度與相對密度之間為線性關(guān)系，反映在迭代過程中，便是單元剛度隨相對密度穩(wěn)定變化。

如果將密度體積插值方法中的中間密度材料看成是一種由空洞和實(shí)體材料組成的兩相復(fù)合材料，則當(dāng)實(shí)體材料的占比為ω時(shí)，由兩相材料的Hashin-Shtrikman上下限[9]可知，體積模量κ和剪切模量μ需要滿足：

(5)

(6)

(7)

(8)

式中，κ0和μ0分別為實(shí)體材料的體積模量和剪切模量；ν0為實(shí)體材料的泊松比。

將式(7)、式(8)代入到式(5)、式(6)，并且設(shè)定中間密度材料的泊松比與實(shí)體材料的一致，整理可得

(9)

(10)

將ω的兩種形式(式(2)和式(3))代入到式(9)和式(10)，可以得到中間密度材料滿足Hashin-Shtrikman上下限時(shí)，指數(shù)型和有理型懲罰函數(shù)中懲罰因子的取值范圍:

(11)

(12)

當(dāng)泊松比ν0取1/3時(shí)，可以得到

(13)

密度體積插值方法不同于密度剛度插值方法中通過弱化中間密度材料的剛度以懲罰中間密度，而是采用對體積進(jìn)行懲罰，映射到體積最小化問題上，即通過對中間密度所對應(yīng)的目標(biāo)函數(shù)進(jìn)行懲罰，以達(dá)到懲罰中間密度的目的。對比密度剛度插值方法和密度體積插值方法的優(yōu)化迭代過程，對于相同的步長，即相同的相對密度變化，密度剛度插值方法因其剛度與相對密度呈懲罰關(guān)系而導(dǎo)致單元剛度變化劇烈且不穩(wěn)定；密度體積插值方法中剛度與相對密度呈線性關(guān)系，因而單元剛度變化穩(wěn)定。單元剛度變化的穩(wěn)定關(guān)系著計(jì)算過程的穩(wěn)定，所以密度體積插值方法具有更穩(wěn)定的優(yōu)化過程。

2 優(yōu)化模型與求解策略

位移約束條件下求結(jié)構(gòu)最小體積的優(yōu)化問題可以表示為

(14)

由密度體積插值方法，可得體積敏度為

(15)

通過虛載荷法可以求解位移敏度為

(16)

式中，vUj和vuji分別為虛載荷法求解節(jié)點(diǎn)位移敏度時(shí)的虛載荷總位移矩陣和虛載荷單元位移矩陣；ki0為實(shí)體單元剛度矩陣；ui為第i個(gè)單元的位移向量；U為全局位移矩陣。

根據(jù)Kuhn-Tucher最優(yōu)化條件，在最優(yōu)解處滿足：

(17)

式中，λj為拉格朗日乘子。

根據(jù)式(17)及設(shè)計(jì)變量上下限構(gòu)造迭代公式如下：

(18)

其中，k為迭代次數(shù)，α為松弛因子，用來控制設(shè)計(jì)變量變化幅度，保證迭代穩(wěn)定收斂，本文中取α為0.5。Mi的表達(dá)式為

(19)

當(dāng)ω取式(2)中ω1時(shí)，可得到

(20)

當(dāng)ω取式(3)中ω2時(shí)，可得到

(21)

迭代公式中的拉格朗日乘子λj可以通過對偶方法[10-11]進(jìn)行求解。

為解決優(yōu)化過程中常出現(xiàn)的棋盤格[12]問題，采用基于敏度的過濾技術(shù)[13]對位移敏度進(jìn)行過濾處理。對于任意單元b，以其中心為圓心，以rmin為半徑作圓，圓內(nèi)單元參與單元b的過濾計(jì)算，過濾后敏度為

(22)

權(quán)重因子Hi為

Hi=rmin-dist(i,b) {i∈N|dist(i,b)≤rmin)}

(23)

其中，dist(i,b)為過濾范圍內(nèi)單元i和中心單元b的距離，N為過濾范圍內(nèi)單元的數(shù)量。過濾處理起到均化作用，可有效消除數(shù)值奇異現(xiàn)象，過濾半徑rmin一般取單元尺寸的1～3倍。

3 數(shù)值算例及分析

3.1 算例1

如圖2所示，設(shè)計(jì)域大小為80 mm×40 mm，厚度為1 mm，下端面左右兩端固定，中間位置作用豎直向下載荷F，大小為1000 N。位移約束為載荷作用點(diǎn)豎直向下方向上的位移，大小為0.1 mm。材料的彈性模量E0=210 GPa，泊松比ν0=0.3，過濾半徑rmin=2.5 mm，初始設(shè)計(jì)變量取1。圖3a、圖3c、圖3e是基于指數(shù)型懲罰函數(shù)ω1的優(yōu)化結(jié)果，懲罰因子p分別為0.3，0.2，0.1。圖3b、圖3d、圖3f是基于有理型懲罰函數(shù)ω2的優(yōu)化結(jié)果，懲罰因子q分別為0.7，0.8，0.9。圖4所示為采用不同懲罰因子結(jié)構(gòu)體積的變化過程。表1為算例1優(yōu)化結(jié)果數(shù)據(jù)。

圖2 算例1設(shè)計(jì)域和邊界條件Fig.2 Design domain and boundary conditions for example 1

(a)p=0.3 (b)q=0.7

(c)p=0.2 (d)q=0.8

(e)p=0.1 (f)q=0.9圖3 算例1優(yōu)化結(jié)果Fig.3 Topology results for example 1

(a)指數(shù)型懲罰函數(shù)

(b)有理型懲罰函數(shù)圖4 算例1體積變化過程Fig.4 Volume change process of example 1

表1 算例1優(yōu)化結(jié)果數(shù)據(jù)對比

從圖3中的優(yōu)化結(jié)果可以看出，基于指數(shù)型和有理型懲罰函數(shù)得到的優(yōu)化結(jié)果拓?fù)湫问揭恢?。從圖4可以看出，加大懲罰程度之后，迭代過程沒有出現(xiàn)較大浮動，算法較穩(wěn)定。從表1中數(shù)據(jù)可以看出，基于指數(shù)型密度體積插值方法的優(yōu)化，隨懲罰因子p的減小(即懲罰程度加大)，優(yōu)化結(jié)果的體積小幅度增大，灰度單元減少，迭代次數(shù)減少?；谟欣硇兔芏润w積插值方法的優(yōu)化，隨懲罰因子q的增大(即懲罰程度加大)，優(yōu)化結(jié)果的體積小幅度增大，灰度單元減少，迭代次數(shù)減少。對比兩種懲罰形式的迭代次數(shù)，可以發(fā)現(xiàn)有理型懲罰的迭代次數(shù)相對較少，在優(yōu)化效率上占有優(yōu)勢。p=0.1和q=0.7時(shí)，優(yōu)化結(jié)果的體積和灰度占比基本相同。

3.2 算例2

如圖5所示，設(shè)計(jì)域大小為80 mm×40 mm，厚度為1 mm，下端面左端固定，右端豎直方向約束，A、B、C位置作用豎直向下的載荷F，大小為1000 N。位移約束為A、B、C三點(diǎn)豎直向下方向上的位移，大小都為0.15 mm。材料的彈性模量E0=210 GPa，泊松比ν0=0.3，過濾半徑rmin=2.5 mm。圖6a、圖6b分別是基于本文密度體積插值方法的兩種懲罰形式的優(yōu)化結(jié)果。圖6c、圖6d分別是基于SIMP和RAMP插值方法的優(yōu)化結(jié)果。表2所示為基于不同插值方法的優(yōu)化結(jié)果數(shù)據(jù)對比。

圖5 算例2設(shè)計(jì)域和邊界條件Fig.5 Design domain and boundary conditions for example 2

(a)指數(shù)型密度體積 (b)有理型密度體積 插值方法 插值方法

(c)SIMP (d)RAMP圖6 算例2優(yōu)化結(jié)果Fig.6 Topology results for example 2

由圖6可以看出，基于本文密度體積插值方法的優(yōu)化結(jié)果與SIMP和RAMP方法的優(yōu)化結(jié)果拓?fù)湫问揭恢?，而?yōu)化結(jié)果明顯較清晰。從表2可以看出，本文密度體積插值方法的優(yōu)化結(jié)果與SIMP和RAMP方法的優(yōu)化結(jié)果相比，結(jié)構(gòu)體積明顯較小，且灰度單元較少。這是由于基于密度體積插值方法的優(yōu)化中可以取得對中間密度更大程度的懲罰。

表2 算例2優(yōu)化結(jié)果數(shù)據(jù)對比

4 結(jié)語

密度體積插值方法是通過構(gòu)造單元體積與相對密度之間的懲罰關(guān)系實(shí)現(xiàn)懲罰中間密度的。該方法在位移約束求最小體積優(yōu)化問題上的應(yīng)用結(jié)果顯示，隨密度體積懲罰函數(shù)懲罰程度的增大，優(yōu)化結(jié)果較SIMP和RAMP方法的灰度單元更少，且沒有出現(xiàn)懲罰程度增大，優(yōu)化不穩(wěn)定甚至落入局部最優(yōu)解等現(xiàn)象。

[1] BENDSOE M P, KIKUCHI N. Generating Optimal Topologies in Structural Design Using a Homogenization Method[J]. Computer Methods in Applied Mechanics and Engineering,1988,71(2):197-224.

[2] BENDSOE M P. Optimal Shape Design as a Material Distribution Problem[J]. Structural Optimization,1989,1(4):193-202.

[3] BENDSOE M P, Sigmund O. Material Interpolation Schemes in Topology Optimization[J]. Archive of Applied Mechanics,1999,69(9/10):635-654.

[4] STOLPE M, SVANBERG K. An Alternative Interpolation Scheme for Minimum Compliance Topology Optimization[J]. Structural and Multidisciplinary Optimization,2001,22(2):116-124.

[5] 羅震, 陳立平, 黃玉盈, 等. 連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)[J]. 力學(xué)進(jìn)展,2004,34(4):463-476. LUO Zhen, CHEN Liping, HUANG Yuying, et al. Topological Optimization Design for Continuum Structures[J]. Advances in Mechanics,2004,34(4):463-476.

[6] 楊陳, 郝志勇, 劉保林, 等. 變密度法在氣缸體輕量化設(shè)計(jì)中的應(yīng)用[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版),2009,43(5):916-919. YANG Chen, HAO Zhiyong, LIU Baolin, et al. Application of Variable Density Method to Light-weight Design of Cylinder Block[J]. Journal of Zhejiang University (Engineering Science),2009,43(5):916-919.

[7] 朱劍峰, 林逸, 陳瀟凱, 等. 汽車變速箱殼體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)[J]. 吉林大學(xué)學(xué)報(bào)(工學(xué)版),2013,43(3):584-589. ZHU Jianfeng, LIN Yi, CHEN Xiaokai, et al. Structural Topology Optimization Based Design of Automotive Transmission Housing Structure[J]. Journal of Jilin University (Engineering and Technology Edition),2013,43(3):584-589.

[8] 羅震, 陳立平, 鐘毅芳. 基于凸規(guī)劃方法的計(jì)算機(jī)輔助結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2005,17(4):774-781. LUO Zhen, CHEN Liping, ZHONG Yifang. Computer-aided Topology Optimization Using Convex Programming Method[J]. Journal of Computer-aided Design and Computer Graphics,2005,17(4):774-781.

[9] HASHIN Z, SHTRIKMAN S. A Variational Approach to the Theory of the Elastic Behaviour of Multiphase Materials[J]. Journal of the Mechanics and Physics of Solids,1963,11(2):127-140.

[10] BECKERS M. Topology Optimization Using a Dual Method with Discrete Variables[J]. Structural Optimization,1999,17(1):14-24.

[11] JOG C S. A Dual Algorithm for the Topology Optimization of Non-linear Elastic Structures[J]. International Journal for Numerical Methods in Engineering,2009,77(4):502-517.

[12] SIGMUND O, PETERSSON J. Numerical Instabilities in Topology Optimization: a Survey on Procedures Dealing with Checkerboards, Mesh-dependencies and Local Minima[J]. Structural Optimization,1998,16(1):68-75.

[13] SIGMUND O. Morphology-based Black and White Filters for Topology Optimization[J]. Structural and Multidisciplinary Optimization,2007,33(4):401-424.

(編輯 王艷麗)

Structural Topology Optimization Based on Density-volume Interpolation Scheme

XU Xiaokui GUO Baofeng JIN Miao

Key Laboratory of Advanced Forging & Stamping Technology and Science, Yanshan University, Qinhuangdao,Hebei,066004

A density-volume interpolation scheme applied in the variable density method was proposed, which was used for controlling grayscale elements in the topology optimization of continuum structures. In the new interpolation scheme, a linear relationship was established between the element stiffness and the relative density to ensure that the element stiffness changed smoothly with the changes of relative density in each iteration. And a penalty function was established between the element volume and the relative density in order to penalize the intermediate densities, which is better for reducing the number of grayscale elements acompained with the intermediate densities approaching both ends.The new interpolation scheme was introduced into the minimum volume optimization problems subjected to displacement constraints. The examples show that the optimization process changes little when increasing the penalty and it indicates that this algorithm is relatively stable. When solving the topology optimization problems with the same structures, it shows that there are less grayscale elements in the optimization results by using the new interpolation scheme than the results optimized by the SIMP and RAMP method. It is because a large penalty may be implemented on intermediate density in the new interpolation scheme.

interpolation scheme; variable density method; displacement constraint; topology optimization

2016-10-18

國家自然科學(xué)基金資助項(xiàng)目(51575474)；河北省自然科學(xué)基金資助項(xiàng)目(E2015203223)

TH122;O343

10.3969/j.issn.1004-132X.2017.11.003

許小奎，男，1990年生。燕山大學(xué)機(jī)械工程學(xué)院博士研究生。主要研究方向?yàn)槌尚卧O(shè)備結(jié)構(gòu)分析與優(yōu)化設(shè)計(jì)。E-mail:xxkysu@126.com。郭寶峰，男，1958年生。燕山大學(xué)機(jī)械工程學(xué)院教授、博士研究生導(dǎo)師。金 淼，男，1968年生。燕山大學(xué)機(jī)械工程學(xué)院教授、博士研究生導(dǎo)師。