基于改進(jìn)SIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法研究*

李啟宏,李海艷

(廣東工業(yè)大學(xué) 機(jī)電工程學(xué)院，廣東 廣州 510006)

0 引 言

變密度法[1]是連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化[2]的常用方法，是在工程領(lǐng)域中應(yīng)用最廣泛、商用最為成功的拓?fù)鋬?yōu)化方法。變密度法采用材料單元的密度作為設(shè)計(jì)變量，其設(shè)計(jì)變量的數(shù)量很少，容易通過編程實(shí)現(xiàn)，計(jì)算效率也高。

SIGMUND O等[3]對(duì)基于變密度法的建模方法進(jìn)行了卓有成效的研究，對(duì)不同材料中間密度懲罰方法進(jìn)行了分析與對(duì)比，最終提出了固體各向同性懲罰微結(jié)構(gòu)模型，即SIMP模型(solid isotropic microstructures with penali-zation , SIMP)[4];MARTINEZ J M[5]對(duì)SIMP法的理論收斂性進(jìn)行了研究;STOLPE M等[6]提出了材料屬性的有理近似模型(rational approximation of material properties , RAMP)，其使用有理式對(duì)材料中間密度進(jìn)行了懲罰;陳祥等[7]采用RAMP法建立了基于變密度理論的優(yōu)化模型。目前，SIMP法已經(jīng)在不同的結(jié)構(gòu)優(yōu)化問題中得到了成功的應(yīng)用，其引入一種假想的相對(duì)密度在0～1之間可變的材料，在一定的材料用量的條件下尋找具有最大剛度(結(jié)構(gòu)的柔順性最小)的結(jié)構(gòu)材料最佳分布形式。

根據(jù)文獻(xiàn)[8，9]的大量實(shí)驗(yàn)研究可發(fā)現(xiàn)，在使用SIMP法解決連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題時(shí)，其優(yōu)化結(jié)構(gòu)的剛度性能與力學(xué)響應(yīng)分析的計(jì)算精度成正比關(guān)系。在網(wǎng)格規(guī)模較小時(shí)，SIMP法的力學(xué)響應(yīng)分析的計(jì)算精度欠佳，導(dǎo)致基于SIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化得到的優(yōu)化結(jié)構(gòu)的剛度偏小。增大網(wǎng)格規(guī)?？商岣咂淞W(xué)響應(yīng)分析的計(jì)算精度，雖然優(yōu)化結(jié)構(gòu)的剛度得到了提高，但受限于網(wǎng)格依賴性等數(shù)值不穩(wěn)定問題[10]，其優(yōu)化結(jié)構(gòu)的細(xì)小分支增多，導(dǎo)致結(jié)構(gòu)的可制造性很差。

針對(duì)SIMP法的這些問題，筆者對(duì)SIMP法的數(shù)學(xué)模型進(jìn)行改進(jìn)，在力學(xué)響應(yīng)分析理論上引入新型力學(xué)數(shù)值求解方法：?jiǎn)卧⒎址?EDM)[11,12]以全體單元中心點(diǎn)應(yīng)力應(yīng)變積之和最小作為優(yōu)化目標(biāo)函數(shù)，再結(jié)合比例拓?fù)鋬?yōu)化法(PTO)[13]對(duì)改進(jìn)的數(shù)學(xué)模型進(jìn)行數(shù)值優(yōu)化求解；通過懸臂梁案例驗(yàn)證改進(jìn)的SIMP法，解決連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化存在的問題。

1 SIMP法數(shù)學(xué)模型的改進(jìn)

1.1 力學(xué)響應(yīng)分析方法

FEM是SIMP法力學(xué)響應(yīng)分析的理論基礎(chǔ)，具有求解簡(jiǎn)便、數(shù)值計(jì)算方法成熟等優(yōu)勢(shì)。同時(shí)FEM也具有一定的局限性，在網(wǎng)格規(guī)模較小時(shí)，其計(jì)算精度不高。本小節(jié)研究的EDM是GAO Xiao-wei等針對(duì)熱力耦合問題提出的一種新型力學(xué)數(shù)值求解方法，其采用離散有限單元的分析形式，且對(duì)節(jié)點(diǎn)受力平衡的二階偏微分方程[14]的求解具有更高的力學(xué)分析計(jì)算精度。本研究將EDM作為SIMP法數(shù)學(xué)模型的力學(xué)響應(yīng)分析方法，以提高計(jì)算精度，接下來推導(dǎo)其力學(xué)響應(yīng)分析公式。

在節(jié)點(diǎn)方程配置上，傳統(tǒng)的配置方法(traditional collocation method ,TCM)[15]通常是在所有節(jié)點(diǎn)中配置平衡方程，并在邊界節(jié)點(diǎn)處添加邊界條件。而EDM在微分單元的內(nèi)部節(jié)點(diǎn)中配置微分平衡方程，在交界節(jié)點(diǎn)和邊界節(jié)點(diǎn)中配置牽引平衡方程。

對(duì)于微分單元的內(nèi)部節(jié)點(diǎn)，其微分平衡控制方程如下：

(1)

式中：bi—體力；Ω—設(shè)計(jì)域；D—單元彈性矩陣；N—微分單元的形函數(shù)矩陣；u—單元節(jié)點(diǎn)位移；ξ—單元內(nèi)部的等參坐標(biāo)，對(duì)于平面二維問題有ξ=(ξ，η)。

對(duì)于由多個(gè)單元共有的交界節(jié)點(diǎn)，其包含該節(jié)點(diǎn)的所有表面，滿足牽引平衡條件：

(2)

交界節(jié)點(diǎn)的牽引平衡方程如下式：

(3)

對(duì)于位于結(jié)構(gòu)外部的邊界節(jié)點(diǎn)，牽引力作用于該節(jié)點(diǎn)，也滿足牽引平衡條件。

邊界節(jié)點(diǎn)牽引力平衡方程為：

(4)

將結(jié)構(gòu)設(shè)計(jì)離散為有限單元，對(duì)整體節(jié)點(diǎn)進(jìn)行編號(hào)，將單元的節(jié)點(diǎn)配置方程按節(jié)點(diǎn)編號(hào)疊加到相應(yīng)的位置，形成全局節(jié)點(diǎn)配置方程組；最終方程組中的未知數(shù)在全局節(jié)點(diǎn)序列中進(jìn)行編號(hào)，并將邊界條件的指定位移和牽引力代入節(jié)點(diǎn)配置方程中。

最終的全局節(jié)點(diǎn)配置方程組形式如下：

Au=b

(5)

式中：u—整體節(jié)點(diǎn)的位移；b—整體節(jié)點(diǎn)的牽引力；A—由每個(gè)節(jié)點(diǎn)方程的系數(shù)構(gòu)成的系數(shù)矩陣。

將式(5)作為SIMP法數(shù)學(xué)模型的力學(xué)響應(yīng)分析式。

1.2 目標(biāo)函數(shù)優(yōu)化

SIMP法是以柔順性最小為優(yōu)化目標(biāo)進(jìn)行拓?fù)鋬?yōu)化，其柔順性越小，優(yōu)化結(jié)構(gòu)的剛度越大。提高SIMP法的力學(xué)響應(yīng)分析計(jì)算精度可得到柔順性更小的優(yōu)化結(jié)構(gòu)，而理論上更小的柔順性常常與結(jié)構(gòu)的合理性相矛盾。當(dāng)網(wǎng)格規(guī)模很大時(shí)，雖然SIMP法的力學(xué)分析計(jì)算精度很高，但同時(shí)優(yōu)化結(jié)構(gòu)會(huì)出現(xiàn)很多不同的局部最優(yōu)結(jié)構(gòu)，表現(xiàn)為優(yōu)化結(jié)構(gòu)出現(xiàn)很多細(xì)小的分支。而過多的細(xì)小分支不符合實(shí)際的生產(chǎn)理念，其可制造性很差。

SIMP法數(shù)學(xué)模型是連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化最為常用的密度插值模型[16]，以柔順性最小為優(yōu)化目標(biāo)，即追求優(yōu)化結(jié)構(gòu)的綜合形變最小。受該思想原理的啟發(fā)，以追求優(yōu)化結(jié)構(gòu)的最佳綜合應(yīng)力應(yīng)變表現(xiàn)為出發(fā)點(diǎn)，再結(jié)合EDM含內(nèi)部節(jié)點(diǎn)的特點(diǎn)，筆者提出以全體單元中心點(diǎn)應(yīng)力應(yīng)變積之和最小作為SIMP法數(shù)學(xué)模型的優(yōu)化目標(biāo)。

將單元中心點(diǎn)應(yīng)力應(yīng)變積與單元的插值密度x和材料彈性模量E相關(guān)聯(lián)，則優(yōu)化目標(biāo)函數(shù)構(gòu)造如下：

(6)

式中：xi—單元密度設(shè)計(jì)變量；N—單元數(shù)量；p—懲罰因子；σi—單元中心點(diǎn)的應(yīng)變；εi—單元中心點(diǎn)的應(yīng)變；Bi—單元應(yīng)變矩陣；E0—實(shí)體材料彈性模量；Emin—空洞材料彈性模量。

在追求綜合形變最小時(shí)，柔順性目標(biāo)函數(shù)因忽略實(shí)際生產(chǎn)理念，過多的細(xì)小分支導(dǎo)致優(yōu)化結(jié)構(gòu)的可制造性很差。筆者提出的目標(biāo)函數(shù)更加注重優(yōu)化結(jié)構(gòu)每個(gè)單元的應(yīng)力和應(yīng)變的狀態(tài)，為結(jié)合新型數(shù)值求解方法以減少細(xì)小分支數(shù)量的研究打下基礎(chǔ)。

2 改進(jìn)的SIMP法

2.1 數(shù)學(xué)模型的構(gòu)建

由于前面筆者已分別用式(5，6)取代SIMP法數(shù)學(xué)模型中，原有的FEM力學(xué)響應(yīng)分析理論和最小柔順性目標(biāo)函數(shù)，至此，已完成對(duì)SIMP法數(shù)學(xué)模型的改進(jìn)，并提出一種改進(jìn)的SIMP法(ISIMP法)，其數(shù)學(xué)模型構(gòu)建如下：

(7)

式中：x—單元密度矩陣；xmin—單元密度最小值，取0.001；C(X)—優(yōu)化目標(biāo)函數(shù)；V—優(yōu)化后的結(jié)構(gòu)體積；V0—結(jié)構(gòu)的初始體積；f—體積約束參數(shù)。

2.2 數(shù)值優(yōu)化求解方法

PTO法是由Biyikli于2015年提出的一種新型拓?fù)鋬?yōu)化數(shù)值求解算法。在應(yīng)力約束問題和最小柔順性問題中，根據(jù)每個(gè)單元所占總體應(yīng)力和總體柔順性的應(yīng)力比和柔順性比，PTO法將材料分配給相應(yīng)單元。而全體單元中心點(diǎn)應(yīng)力應(yīng)變積之和目標(biāo)函數(shù)涉及到最能反應(yīng)受力結(jié)構(gòu)力學(xué)狀態(tài)的單元應(yīng)力與應(yīng)變，所以筆者結(jié)合PTO法對(duì)ISIMP法數(shù)學(xué)模型進(jìn)行數(shù)值優(yōu)化求解，在每次迭代優(yōu)化中，將材料根據(jù)每個(gè)單元中心點(diǎn)應(yīng)力應(yīng)變積占全體單元應(yīng)力應(yīng)變積之和的比例大小，分配給相應(yīng)單元。其基本步驟如下：

(1)在計(jì)算單元中心點(diǎn)，應(yīng)力應(yīng)變積占總體單元應(yīng)力應(yīng)變積之和的比例為：

(8)

式中：Ci—第個(gè)單元的應(yīng)力應(yīng)變積；Cpro—單元貢獻(xiàn)比例；n—單元總數(shù)。

(2)在內(nèi)循環(huán)中，將剩余體積按單元的應(yīng)力應(yīng)變積占總體單元應(yīng)力應(yīng)變積和的比例大小分配給每個(gè)單元，其分配公式為：

(9)

(3)當(dāng)剩余體積少于0.001時(shí)，退出內(nèi)循環(huán)，完成一次PTO法迭代更新。此處定義歷史平衡系數(shù)α，本文取α=0.5，用于平衡當(dāng)前更新得出的密度與上一次更新得出的密度，其更新方案可表示為：

(10)

2.3 應(yīng)用與算法設(shè)計(jì)

為解決SIMP法在求解連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題時(shí)存在的缺陷，筆者將ISIMP法應(yīng)用在連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化上?；贗SIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法的主要步驟為：

(1)將設(shè)計(jì)域離散為有限單元網(wǎng)格,并對(duì)網(wǎng)格節(jié)點(diǎn)進(jìn)行編號(hào)；定義約束條件及載荷等邊界條件；

(2)按照不同的節(jié)點(diǎn)類型配置節(jié)點(diǎn)方程，根據(jù)節(jié)點(diǎn)編號(hào)組裝全局節(jié)點(diǎn)的EDM方程組；

(3)求解EDM方程組得到節(jié)點(diǎn)位移矢量，計(jì)算全體單元中心應(yīng)力應(yīng)變積之和目標(biāo)函數(shù)；

(4)并結(jié)合PTO法對(duì)ISIMP數(shù)學(xué)模型進(jìn)行數(shù)值優(yōu)化求解，以更新單元密度；

(5)若連續(xù)兩次迭代的單元密度最大變化小于0.01，或迭代次數(shù)達(dá)到100次，則終止迭代優(yōu)化，并輸出單元密度分布結(jié)果；否則重復(fù)步驟(3，4)。

基于ISIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化流程圖如圖1所示。

圖1 基于ISIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化流程圖

3 懸臂梁實(shí)驗(yàn)

懸臂梁(實(shí)驗(yàn)無量綱)長(zhǎng)寬比為3∶1，梁左側(cè)固定約束，梁的上邊界施加大小為1 000的均勻分布載荷(將均勻分布載荷等效為集中力作用到節(jié)點(diǎn)上)。

懸臂梁拓?fù)鋬?yōu)化實(shí)驗(yàn)的基本參數(shù)設(shè)置為：結(jié)構(gòu)彈性模量E0=210，材料的泊松比μ=0.3，懲罰因子p=3，過濾半徑r=1.5，材料的體積約束為0.4。

懸臂梁模型圖如圖2所示。

圖2 懸臂梁模型圖

3.1 實(shí)驗(yàn)基本設(shè)定

為保證懸臂梁始終受均布載荷，實(shí)驗(yàn)中將懸臂梁上邊厚度為寬度乘0.2的區(qū)域設(shè)置為密度不更新區(qū)域，即使該區(qū)域單元密度保持為1。為保持?jǐn)?shù)值求解方法的一致性，筆者基于SIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化同樣結(jié)合PTO法進(jìn)行數(shù)值優(yōu)化求解。為了對(duì)基于SIMP法與ISIMP法的懸臂梁優(yōu)化結(jié)構(gòu)的性能進(jìn)行對(duì)比，該實(shí)驗(yàn)以優(yōu)化結(jié)構(gòu)的柔順性作為對(duì)比指標(biāo)。優(yōu)化結(jié)構(gòu)的柔順性越小，其剛度則越大。

對(duì)基于ISIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的單元密度分布結(jié)果，筆者同樣使用FEM進(jìn)行力學(xué)分析，計(jì)算優(yōu)化結(jié)構(gòu)的節(jié)點(diǎn)位移，從而計(jì)算優(yōu)化結(jié)構(gòu)的柔順性大小。

3.2 懸臂梁拓?fù)鋬?yōu)化

首先，在不同網(wǎng)格規(guī)模下，筆者分別使用SIMP法與ISIMP法的力學(xué)響應(yīng)分析理論基礎(chǔ)，對(duì)懸臂梁進(jìn)行力學(xué)響應(yīng)分析，通過對(duì)比懸臂梁右下角節(jié)點(diǎn)的Y方向位移uy(該實(shí)驗(yàn)以3 000×1 000的超大規(guī)模網(wǎng)格進(jìn)行FEM力學(xué)響應(yīng)分析的結(jié)果作為假定理論值，uy=-215.298 0)，來比較SIMP法的FEM與ISIMP法的EDM在力學(xué)響應(yīng)分析上的計(jì)算精度。

懸臂梁右下角節(jié)點(diǎn)的Y方向位移及誤差如表1所示。

表1 懸臂梁右下角節(jié)點(diǎn)的Y方向位移

由表1可知，EDM的力學(xué)響應(yīng)分析的計(jì)算誤差隨網(wǎng)格規(guī)模的增大而減小，具有良好的收斂性。與FEM相比，EDM在同一網(wǎng)格規(guī)模的力學(xué)響應(yīng)分析中具有更高的計(jì)算精度；在網(wǎng)格規(guī)模較小的力學(xué)響應(yīng)分析中，精度提升更為明顯。

接著，在不同網(wǎng)格規(guī)模下，筆者分別進(jìn)行了基于SIMP法與ISIMP法的懸臂梁拓?fù)鋬?yōu)化實(shí)驗(yàn)。

基于ISIMP法，在不同網(wǎng)格規(guī)模下的懸臂梁拓?fù)鋬?yōu)化的優(yōu)化結(jié)構(gòu)柔順性的變化曲線圖，如圖3所示。

由圖3可知，在迭代優(yōu)化過程中，基于ISIMP法的懸臂梁拓?fù)鋬?yōu)化的柔順性逐步減小，具有良好的收斂性，由此可以驗(yàn)證ISIMP法在連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化上的可行性。

兩種方法的懸臂梁優(yōu)化結(jié)構(gòu)的柔順性，以及基于ISIMP法的優(yōu)化結(jié)構(gòu)柔順性，相對(duì)于基于SIMP法的優(yōu)化結(jié)構(gòu)柔順性的降低率，如表2所示。

表2 懸臂梁優(yōu)化結(jié)構(gòu)的柔順性

結(jié)合表(1，2)數(shù)據(jù)可知：SIMP法的FEM與ISIMP法的EDM的力學(xué)響應(yīng)分析的計(jì)算精度，都隨著網(wǎng)格規(guī)模的增大而提高，懸臂梁優(yōu)化結(jié)構(gòu)的柔順性都隨之減小，即剛度隨之增大；在同一網(wǎng)格規(guī)模下，EDM的力學(xué)響應(yīng)分析的計(jì)算精度比FEM更高；基于ISIMP法的懸臂梁拓?fù)鋬?yōu)化得到了柔順性更小的優(yōu)化結(jié)構(gòu)，即懸臂梁優(yōu)化結(jié)構(gòu)的剛度性能得到了提升。

最后，在小規(guī)模網(wǎng)格和大規(guī)模網(wǎng)格下，筆者對(duì)懸臂梁的實(shí)驗(yàn)結(jié)果進(jìn)行分析。

懸臂梁的拓?fù)鋬?yōu)化實(shí)驗(yàn)結(jié)果圖如圖4所示。

圖4 懸臂梁拓?fù)鋬?yōu)化實(shí)驗(yàn)結(jié)果圖

由表1可知，在小規(guī)模網(wǎng)格實(shí)驗(yàn)中，相比FEM，EDM的力學(xué)響應(yīng)分析的計(jì)算精度優(yōu)勢(shì)更大。

由圖4(a)可知：在網(wǎng)格規(guī)模為60×20實(shí)驗(yàn)中，兩種方法的優(yōu)化結(jié)構(gòu)大致相同，但基于ISIMP法的優(yōu)化結(jié)構(gòu)柔順性相對(duì)于SIMP法的減小了5.41%，剛度明顯提升；

圖4(c～d)記錄了網(wǎng)格規(guī)模為300×100的懸臂梁拓?fù)鋬?yōu)化的迭代次數(shù)為20次、40次和60次時(shí)，單元中心點(diǎn)應(yīng)力應(yīng)變積與單元插值密度的分布。由此可知，基于ISIMP法的懸臂梁拓?fù)鋬?yōu)化在前40次迭代中，逐步形成穩(wěn)定的傳力路徑區(qū)域，且在傳力路徑區(qū)域上的單元中心點(diǎn)應(yīng)力應(yīng)變積明顯高于其他區(qū)域；再結(jié)合PTO法在數(shù)值優(yōu)化求解中，將材料根據(jù)單元中心點(diǎn)應(yīng)力應(yīng)變積占全體單元應(yīng)力應(yīng)變積之和的比例大小，分配給相應(yīng)單元的密度插值方法，那么可知在傳力路徑區(qū)域上，絕大多數(shù)單元將在后續(xù)的迭代優(yōu)化中繼續(xù)保持實(shí)體材料密度。

對(duì)比圖4(b)中的網(wǎng)格規(guī)模為300×100的懸臂梁最終優(yōu)化結(jié)構(gòu)可知，基于ISIMP法的懸臂梁優(yōu)化結(jié)構(gòu)有效地減少了細(xì)小分支的數(shù)量，提高了優(yōu)化結(jié)構(gòu)的可制造性。

4 結(jié)束語

本研究對(duì)SIMP法的數(shù)學(xué)模型進(jìn)行了改進(jìn)，提出了一種基于改進(jìn)SIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法；通過懸臂梁案例對(duì)改進(jìn)的SIMP法進(jìn)行了驗(yàn)證，解決了連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化存在的問題。

與基于SIMP的方法相比，基于ISIMP法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法的優(yōu)勢(shì)有：

(1)在同一網(wǎng)格規(guī)模下，得益于EDM更高精度的力學(xué)響應(yīng)分析，該方法可得到柔順性更小的優(yōu)化結(jié)構(gòu)，即優(yōu)化結(jié)構(gòu)的剛度性能得到了提升；并且在網(wǎng)格規(guī)模較小時(shí)，剛度提升更加明顯；

(2)全體單元應(yīng)力應(yīng)變積之和目標(biāo)函數(shù)更加注重優(yōu)化結(jié)構(gòu)每個(gè)單元的應(yīng)力和應(yīng)變的狀態(tài)，再結(jié)合PTO法進(jìn)行數(shù)值優(yōu)化求解的方法，在大規(guī)模網(wǎng)格下，可有效減少細(xì)小分支的數(shù)量，提高優(yōu)化結(jié)構(gòu)的可制造性。