基于密度體積插值方法的結構拓撲優化

許小奎 郭寶峰 金 淼

燕山大學先進鍛壓成形技術與科學教育部重點實驗室，秦皇島，066004

?

基于密度體積插值方法的結構拓撲優化

許小奎 郭寶峰 金 淼

燕山大學先進鍛壓成形技術與科學教育部重點實驗室，秦皇島，066004

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

插值方法；變密度法；位移約束；拓撲優化

0 引言

自1988年BENDSOE等[1]提出用于求解連續體拓撲優化的均勻化方法以來，連續體結構拓撲優化一直是結構優化領域研究的熱點。變密度法[2-5]是在均勻化方法的基礎上發展形成的一種算法，該方法由于建模簡單、易于求解等特點而得到了廣泛應用[6-7]。

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

針對上述問題，本文提出了一種密度體積插值方法，并采用該方法求解了位移約束條件下的結構拓撲優化問題。

1 密度體積插值方法

密度體積插值方法中為實現對中間密度的懲罰，構造單元體積與相對密度之間的關系為

Vi=Vi0ω(ρi)

(1)

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

懲罰函數ω可為指數型懲罰函數：

(2)

亦可為有理型懲罰函數：

(3)

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

(a)ω1

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

密度體積插值方法中，構造彈性模量與相對密度之間為線性函數：

Ei=ρiE0

(4)

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

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

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

(5)

(6)

(7)

(8)

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

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

(9)

(10)

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

(11)

(12)

當泊松比ν0取1/3時，可以得到

(13)

密度體積插值方法不同于密度剛度插值方法中通過弱化中間密度材料的剛度以懲罰中間密度，而是采用對體積進行懲罰，映射到體積最小化問題上，即通過對中間密度所對應的目標函數進行懲罰，以達到懲罰中間密度的目的。對比密度剛度插值方法和密度體積插值方法的優化迭代過程，對于相同的步長，即相同的相對密度變化，密度剛度插值方法因其剛度與相對密度呈懲罰關系而導致單元剛度變化劇烈且不穩定；密度體積插值方法中剛度與相對密度呈線性關系，因而單元剛度變化穩定。單元剛度變化的穩定關系著計算過程的穩定，所以密度體積插值方法具有更穩定的優化過程。

2 優化模型與求解策略

位移約束條件下求結構最小體積的優化問題可以表示為

(14)

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

(15)

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

(16)

式中，vUj和vuji分別為虛載荷法求解節點位移敏度時的虛載荷總位移矩陣和虛載荷單元位移矩陣；ki0為實體單元剛度矩陣；ui為第i個單元的位移向量；U為全局位移矩陣。

根據Kuhn-Tucher最優化條件，在最優解處滿足：

(17)

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

根據式(17)及設計變量上下限構造迭代公式如下：

(18)

其中，k為迭代次數，α為松弛因子，用來控制設計變量變化幅度，保證迭代穩定收斂，本文中取α為0.5。Mi的表達式為

(19)

當ω取式(2)中ω1時，可得到

(20)

當ω取式(3)中ω2時，可得到

(21)

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

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

(22)

權重因子Hi為

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

(23)

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

3 數值算例及分析

3.1 算例1

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

圖2 算例1設計域和邊界條件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優化結果Fig.3 Topology results for example 1

(a)指數型懲罰函數

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

表1 算例1優化結果數據對比

從圖3中的優化結果可以看出，基于指數型和有理型懲罰函數得到的優化結果拓撲形式一致。從圖4可以看出，加大懲罰程度之后，迭代過程沒有出現較大浮動，算法較穩定。從表1中數據可以看出，基于指數型密度體積插值方法的優化，隨懲罰因子p的減小(即懲罰程度加大)，優化結果的體積小幅度增大，灰度單元減少，迭代次數減少。基于有理型密度體積插值方法的優化，隨懲罰因子q的增大(即懲罰程度加大)，優化結果的體積小幅度增大，灰度單元減少，迭代次數減少。對比兩種懲罰形式的迭代次數，可以發現有理型懲罰的迭代次數相對較少，在優化效率上占有優勢。p=0.1和q=0.7時，優化結果的體積和灰度占比基本相同。

3.2 算例2

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

圖5 算例2設計域和邊界條件Fig.5 Design domain and boundary conditions for example 2

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

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

由圖6可以看出，基于本文密度體積插值方法的優化結果與SIMP和RAMP方法的優化結果拓撲形式一致，而優化結果明顯較清晰。從表2可以看出，本文密度體積插值方法的優化結果與SIMP和RAMP方法的優化結果相比，結構體積明顯較小，且灰度單元較少。這是由于基于密度體積插值方法的優化中可以取得對中間密度更大程度的懲罰。

表2 算例2優化結果數據對比

4 結語

密度體積插值方法是通過構造單元體積與相對密度之間的懲罰關系實現懲罰中間密度的。該方法在位移約束求最小體積優化問題上的應用結果顯示，隨密度體積懲罰函數懲罰程度的增大，優化結果較SIMP和RAMP方法的灰度單元更少，且沒有出現懲罰程度增大，優化不穩定甚至落入局部最優解等現象。

[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] 羅震, 陳立平, 黃玉盈, 等. 連續體結構的拓撲優化設計[J]. 力學進展,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] 楊陳, 郝志勇, 劉保林, 等. 變密度法在氣缸體輕量化設計中的應用[J]. 浙江大學學報(工學版),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] 朱劍峰, 林逸, 陳瀟凱, 等. 汽車變速箱殼體結構拓撲優化設計[J]. 吉林大學學報(工學版),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] 羅震, 陳立平, 鐘毅芳. 基于凸規劃方法的計算機輔助結構拓撲優化設計[J]. 計算機輔助設計與圖形學學報,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

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

TH122;O343

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

許小奎，男，1990年生。燕山大學機械工程學院博士研究生。主要研究方向為成形設備結構分析與優化設計。E-mail:xxkysu@126.com。郭寶峰，男，1958年生。燕山大學機械工程學院教授、博士研究生導師。金 淼，男，1968年生。燕山大學機械工程學院教授、博士研究生導師。