李啟宏,李海艷

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

0 引 言

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

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

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

針對SIMP法的這些問題，筆者對SIMP法的數學模型進行改進，在力學響應分析理論上引入新型力學數值求解方法：單元微分法(EDM)[11,12]以全體單元中心點應力應變積之和最小作為優(yōu)化目標函數，再結合比例拓撲優(yōu)化法(PTO)[13]對改進的數學模型進行數值優(yōu)化求解；通過懸臂梁案例驗證改進的SIMP法，解決連續(xù)體結構拓撲優(yōu)化存在的問題。

1 SIMP法數學模型的改進

1.1 力學響應分析方法

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

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

對于微分單元的內部節(jié)點，其微分平衡控制方程如下：

(1)

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

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

(2)

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

(3)

對于位于結構外部的邊界節(jié)點，牽引力作用于該節(jié)點，也滿足牽引平衡條件。

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

(4)

將結構設計離散為有限單元，對整體節(jié)點進行編號，將單元的節(jié)點配置方程按節(jié)點編號疊加到相應的位置，形成全局節(jié)點配置方程組；最終方程組中的未知數在全局節(jié)點序列中進行編號，并將邊界條件的指定位移和牽引力代入節(jié)點配置方程中。

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

Au=b

(5)

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

將式(5)作為SIMP法數學模型的力學響應分析式。

1.2 目標函數優(yōu)化

SIMP法是以柔順性最小為優(yōu)化目標進行拓撲優(yōu)化，其柔順性越小，優(yōu)化結構的剛度越大。提高SIMP法的力學響應分析計算精度可得到柔順性更小的優(yōu)化結構，而理論上更小的柔順性常常與結構的合理性相矛盾。當網格規(guī)模很大時，雖然SIMP法的力學分析計算精度很高，但同時優(yōu)化結構會出現(xiàn)很多不同的局部最優(yōu)結構，表現(xiàn)為優(yōu)化結構出現(xiàn)很多細小的分支。而過多的細小分支不符合實際的生產理念，其可制造性很差。

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

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

(6)

式中：xi—單元密度設計變量；N—單元數量；p—懲罰因子；σi—單元中心點的應變；εi—單元中心點的應變；Bi—單元應變矩陣；E0—實體材料彈性模量；Emin—空洞材料彈性模量。

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

2 改進的SIMP法

2.1 數學模型的構建

由于前面筆者已分別用式(5，6)取代SIMP法數學模型中，原有的FEM力學響應分析理論和最小柔順性目標函數，至此，已完成對SIMP法數學模型的改進，并提出一種改進的SIMP法(ISIMP法)，其數學模型構建如下：

(7)

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

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

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

(1)在計算單元中心點，應力應變積占總體單元應力應變積之和的比例為：

(8)

式中：Ci—第個單元的應力應變積；Cpro—單元貢獻比例；n—單元總數。

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

(9)

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

(10)

2.3 應用與算法設計

為解決SIMP法在求解連續(xù)體結構拓撲優(yōu)化問題時存在的缺陷，筆者將ISIMP法應用在連續(xù)體結構的拓撲優(yōu)化上?；贗SIMP法的連續(xù)體結構拓撲優(yōu)化方法的主要步驟為：

(1)將設計域離散為有限單元網格,并對網格節(jié)點進行編號；定義約束條件及載荷等邊界條件；

(2)按照不同的節(jié)點類型配置節(jié)點方程，根據節(jié)點編號組裝全局節(jié)點的EDM方程組；

(3)求解EDM方程組得到節(jié)點位移矢量，計算全體單元中心應力應變積之和目標函數；

(4)并結合PTO法對ISIMP數學模型進行數值優(yōu)化求解，以更新單元密度；

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

基于ISIMP法的連續(xù)體結構拓撲優(yōu)化流程圖如圖1所示。

圖1 基于ISIMP法的連續(xù)體結構拓撲優(yōu)化流程圖

3 懸臂梁實驗

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

懸臂梁拓撲優(yōu)化實驗的基本參數設置為：結構彈性模量E0=210，材料的泊松比μ=0.3，懲罰因子p=3，過濾半徑r=1.5，材料的體積約束為0.4。

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

圖2 懸臂梁模型圖

3.1 實驗基本設定

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

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

3.2 懸臂梁拓撲優(yōu)化

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

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

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

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

接著，在不同網格規(guī)模下，筆者分別進行了基于SIMP法與ISIMP法的懸臂梁拓撲優(yōu)化實驗。

基于ISIMP法，在不同網格規(guī)模下的懸臂梁拓撲優(yōu)化的優(yōu)化結構柔順性的變化曲線圖，如圖3所示。

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

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

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

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

最后，在小規(guī)模網格和大規(guī)模網格下，筆者對懸臂梁的實驗結果進行分析。

懸臂梁的拓撲優(yōu)化實驗結果圖如圖4所示。

圖4 懸臂梁拓撲優(yōu)化實驗結果圖

由表1可知，在小規(guī)模網格實驗中，相比FEM，EDM的力學響應分析的計算精度優(yōu)勢更大。

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

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

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

4 結束語

本研究對SIMP法的數學模型進行了改進，提出了一種基于改進SIMP法的連續(xù)體結構拓撲優(yōu)化方法；通過懸臂梁案例對改進的SIMP法進行了驗證，解決了連續(xù)體結構拓撲優(yōu)化存在的問題。

與基于SIMP的方法相比，基于ISIMP法的連續(xù)體結構拓撲優(yōu)化方法的優(yōu)勢有：

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

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