基于辛格插值模式的連續體結構拓撲優化研究*

李玉剛,楊 林

(1.貴州裝備制造職業學院,貴州 貴陽 550025;2.貴州大學 機械工程學院,貴州 貴陽 550025)

0 引 言

拓撲優化是指在設計領域內尋求最佳的材料分布,從而使所得到的結構達到最佳的或規定的性能。

自1988年BENDSOE M P等人[1]發表了關于均勻化方法的開創性研究論文以來,人們在此基礎上又發展了各種拓撲優化方法,例如,MLEJNEK H P等人[2]提出的變密度法、ALLAIRE G等人[3]提出的水平集法及XIE Yi-min等人[4]提出的漸進結構優化法。其中,變密度法以其數學理論嚴謹、物理意義明確等優點,因而得到了廣泛應用[5]。

HUANG Xiao-dong[6]指出了連續體結構拓撲優化問題本質上是0/1離散變量的組合優化。因此,采用變密度法將不可避免地產生在實際工程難以制造的中間密度單元[7]。BENDSOE M P等人[8]提出了采用懲罰函數法使中間變量對對應單元材料剛度的影響變小,從而使最終優化解趨于0或1。

目前,最流行的懲罰函數模式有:SIGMUND O[9]提出的固體各向同性懲罰微結構插值模式(SIMP),和STOLPE M等人[10]提出的材料屬性合理近似插值模式(rational approximation of material properties,RAMP)。但因兩種模式均具有單向懲罰特性,即總是以降低單元彈性模量和單元剛度來達到懲罰目的,不利于中間密度單元的消除,且易出現最終拓撲結構邊界不清晰、棋盤格等數值不穩定現象[11]。

針對傳統變密度法的不足,筆者提出一種基于辛格插值模式的連續體結構拓撲優化新方法,并通過算例來驗證該方法對二維及三維設計域拓撲優化的可行性,從而有效彌補傳統變密度法的缺陷。

1 辛格插值模式

針對傳統變密度法的單向懲罰性問題,理論上如果能構造一插值模式,其既可以提高“有用”單元彈性模量和單元剛度,又可以降低“無用”單元彈性模量和單元剛度,那么,此類具有雙向促進特性的插值模式比具有單向懲罰特性的插值模式更加高效、實用。

因此,筆者制定以下更加合理的中間單元處理策略:

首先,設定一參考值(對于本研究所涉及[0,1]區間內的設計變量,參考值取0.5),該參考值用以判斷中間密度單元彈性模量和單元剛度的促進方式。然后,提高密度值大于0.5的單元彈性模量和單元剛度,同時,降低密度值小于0.5的單元彈性模量和單元剛度;

然后,根據以上中間密度單元處理策略構造如下辛格插值模式:

利用辛格函數[12](即sinc函數)在[0,1]區間內具有較高斜率的單調特性(可使中間密度單元快速向實體單元進化或快速向空洞單元演化),構建如下材料屬性的辛格插值模式:

(1)

式中:xi—單元密度值;xmin—最小單元密度值;α(xi)—辛格插值函數;τ—雙向促進因子。

其中:τ既可促進xmin≤xi≤0.5的單元彈性模量和剛度提高,也可促進0.5

辛格插值模式雙向促進效果如圖1所示。

圖1 辛格插值模式雙向促進效果圖

SIMP插值模式單向懲罰效果如圖2所示。

圖2 SIMP插值模式單向懲罰效果圖

對比圖(1，2)可知:

隨著雙向促進因子τ或單向懲罰因子p的增大,兩種模式的懲罰或促進力度也隨之增強,但SIMP插值模式的增強力度顯然弱于辛格插值模式。

在基于SIMP插值模式的拓撲優化有限元分析過程中,通常利用基于單元密度設計變量的懲罰函數實現單元材料彈性模量的迭代,以達到修改單元剛度矩陣的目的。

相應地,辛格插值函數與單元剛度矩陣和彈性模量之間的關系可表示為:

(2)

式中:K0—實體材料剛度;Ki—優化后單元剛度;E0—實體材料彈性模量;Ei—優化后單元彈性模量。

在拓撲優化迭代過程中,因作為設計變量的單元密度是動態變化的,相應的單元材料彈性模量或單元剛度也是動態變化的,而事實上,懲罰函數的作用僅是懲罰具有中間密度單元材料的剛度或彈性模量,所以,更合理的策略應是對上一次迭代結果中具有中間密度的單元材料的剛度或彈性模量進行懲罰,而不是總以實體材料剛度和實體材料彈性模量為標準。

因此,式(2)可優化為:

(3)

相應地,在基于傳統SIMP懲罰函數拓撲優化有限元分析過程中,通常采用剛度矩陣對單元密度設計變量的偏導數表示單元剛度對設計變量的敏度。

其表達式如下:

(4)

但實際上,單元密度的敏度計算應滿足材料本構關系,即單元密度與單元剛度和彈性模量之間應是線性關系,因此,單元剛度對單元密度設計變量的敏度應按線性關系計算,然后再利用辛格插值模式對敏度進行處理。

所以,在拓撲優化過程中,應根據單元密度值大小不斷對迭代方向進行調整,這樣即使不采用過濾技術,也能得到清晰邊界的拓撲結構。

因此,基于辛格插值模式的敏度可表示為:

(5)

2 基于辛格插值模式的優化模型及算法

以連續體結構體積約束下,總體柔度最小化為優化目標,在給定邊界條件和載荷工況下,基于辛格插值模式的優化模型為:

(6)

式中:f—結構總體柔度函數;ki—第i單元優化后的剛度矩陣;ui—第i單元優化后的位移向量;V—優化后結構總體積;Vi—單元體積;V*—目標優化體積;F—作用力;K—總剛度矩陣;U—總位移矩陣;

其中:取xmin=0.001 5,以避免總剛度矩陣奇異。

結合式(5),剛度可采用兩種方式計算。其中,一種方式是只對敏度進行辛格插值,筆者稱該方法為第一類辛格插值法,如下式所示:

(7)

另一種方式是既對敏度進行辛格插值,又對剛度進行辛格插值,筆者稱該方法為第二類辛格插值法,如下式所示:

(8)

構造問題(6)的移動漸進算法(the method of moving asymptotes, MMA)[13，14]求解表達式為:

(9)

為方便表達,此處令:

(10)

(11)

并令t表示迭代次數,則展開模型表達式如下:

(12)

其中:

(13)

結合對偶理論,將移動近似凸子問題(13)轉換成易于計算求解的對偶問題,同時令λ表示對偶設計變量,則有:

(14)

其中:

(15)

其迭代終止條件如下:

(16)

其中:取δ=0.01。

3 算 例

3.1 算例一

基結構尺寸為80×50,厚度為1,材料彈性模量為E=1,泊松比為0.3,劃分為60×30個矩形四節點單元;固定點為左下角邊界點和右下角邊界點,且采用自由度全固定約束;集中載荷F=1,作用于上邊界中點;求體積約束分數為0.4時的結構最優拓撲。

為簡化計算,該例中所有設計參數均按無量綱處理,且不會影響優化有效性。

二維設計域模型及有限元模型如圖3所示。

圖3 二維設計域及有限元模型

數值不穩定現象[15]是影響連續體結構拓撲優化結果質量的關鍵因素之一。

在傳統變密度方法拓撲優化過程中,通常采用空間敏度過濾技術以消除優化過程中存在的數值不穩定現象,但該技術易使敏度在其過濾半徑內均勻化,即最終拓撲結構的拓撲邊界易出現模糊現象。

不采用過濾技術時,不同τ值的拓撲結果如圖4所示。

圖4 不采用過濾技術時,不同τ值的拓撲結果

根據圖4,不采用過濾技術時,辛格插值模式所得優化拓撲結果結構邊界清晰,即未出現網格依賴性、棋盤格等現象,且τ取不同值時的拓撲結構基本相同,所以雙向促進因子的取值對拓撲結果影響不大。

相應的,τ取不同值時的優化迭代過程如圖5所示。

圖5 τ取不同值時的優化迭代過程

從圖5可以看出:隨著迭代步數和τ取值的增加,目標函數值變化趨勢基本相同,均呈現出快速減小的規律,且最終收斂值趨于一致;

另外,τ=1～5時的收斂速度也基本相同,且均在第4步迭代后進入穩定收斂狀況。不同的是,隨著τ取值的增加,初始目標函數值快速減小(τ=1時,f1=97.543 7;τ=5時,f1=9.311 3)。

接下來,在辛格插值模式優化過程中,筆者同時采用空間敏度過濾技術(本例中,過濾半徑取1.2)。

同時采用敏度過濾技術時,不同τ值的拓撲結果如圖6所示。

圖6 同時采用敏度過濾技術時,不同τ值的拓撲結果

根據圖6,在辛格插值模式優化過程中,同時采用敏度過濾技術,雙向促進因子τ取不同值時,優化拓撲結果結構邊緣模糊、不清晰,與不用過濾方法時所得最終優化后拓撲結構基本一致。

3.2 算例二

二維設計域尺寸為40 mm×20 mm,厚度為1 mm,材料彈性模量E=200 GPa,劃分為40×20個矩形四節點單元;左邊界全部采用固定約束;集中載荷F=9 000 N作用于右邊界中點。

為了簡化計算過程,上述設計參數按無量綱處理,求具有最小體積和最大剛度的結構最優拓撲,設初始單元密度為0.9。

該算例二維設計域模型及有限元模型如圖7所示。

圖7 二維設計域及有限元模型

采用第一類辛格插值法所得最終拓撲結構如圖8所示。

圖8 第一類辛格插值法所得最終拓撲結構

采用第二類辛格插值法所得最終拓撲結構如圖9所示。

圖9 第二類辛格插值法所得最終拓撲結構

采用SIMP法所得最終拓撲結構如圖10所示。

圖10 SIMP法所得最終拓撲結構

對比圖(8～10)可知:第二類辛格插值方法所得最終拓撲結構與SIMP方法相似,而第一類辛格插值方法拓撲結構更加清晰、簡潔。

第一類辛格插值法優化目標函數值、作用點位移、體積及柔度變化曲線,如圖11所示。

圖11 第一類辛格插值法優化目標函數值、作用點位移、體積及柔度變化曲線

由圖11可知:第一類辛格插值方法優化結果體積為初始體積的0.264倍,在優化過程中,開始時不進行辛格插值計算迭代,直到目標函數值趨于穩定(相對誤差小于0.000 01),得到數學上最優解,迭代步30步為轉換步,之后開始采用辛格插值第一類辛格插值法優化策略,以消除中間密度單元;自30～32步之間迭代有波動,這是由于對中間密度單元進行雙向調整所致,之后計算趨于穩定。

因進行雙向調整,更合理地分布了材料,使得在體積相對一致的情況下,目標函數值快速下降,結構柔度值快速減小,增大了剛度,同時在相同外力作用下,作用點位移快速從-4.406 mm變化到-1.327 mm。

第二類辛格插值法和SIMP方法優化目標函數值、作用點位移、體積及柔度變化曲線,如圖12所示。

圖12 第二類辛格插值法和SIMP方法優化目標函數值、作用點位移、體積及柔度變化曲線

由圖12可知:第二類辛格插值方法和SIMP方法優化拓撲近似,迭代收斂結果大小近似,收斂步數亦相近,但SIMP方法計算數值波動比第二類辛格插值方法大。

綜上所述,與第一類辛格插值法方法優化所得目標函數值(2.526e6)相比,第二類辛格插值方法、SIMP方法所得的目標函數值(9.451e6,1.295e7)要大數倍。所以,第一類辛格插值法方法更加敏捷,可得到更小的目標函數值。

3.3 算例三

三維設計域的尺寸為60×30×10,材料彈性模量為E=1,泊松比為0.33,劃分為60×30×10個立方體八節點單元;一端完全約束,且在另一端自由邊緣向下施加垂直的分布載荷F;求體積約束分數為0.4時的結構最優拓撲。

為簡化計算,該算例中所有設計參數均按無量綱處理,且不會影響優化的有效性。

該算例三維設計域模型及有限元模型分別如圖13所示。

圖13 三維設計域及有限元模型

為了證明所提方法在處理三維問題中的有效性,筆者采用基于辛格插值模式的MMA算法進行三維設計域的拓撲優化,研究雙向促進因子τ取不同值時的拓撲優化結果及優化計算的迭代收斂情況。

τ取不同值時的最終拓撲結構如圖14所示。

圖14 τ取不同值時的最終拓撲結構

由圖14可知:采用基于辛格插值模式的MMA算法進行三維設計域拓撲優化,并經isosurface函數后處理,所得拓撲結構可靠、清晰且不失真,可直接利用3D打印技術制造;且雙向促進因子τ=3和τ=4時的最終結構拓撲較τ=1、τ=2及τ=5時更簡潔。

相應的,τ取不同值時的優化迭代過程如圖15所示。

圖15 τ取不同值時的優化迭代過程

由圖15可知:隨著τ取值的增大,目標函數初始值不斷增大,迭代次數也隨之增多,但目標函數值卻呈現出先減后增的趨勢,且當τ=3時,可以取得最小目標函數值456.449 6;該趨勢與辛格插值模式的雙向促進特性相吻合。

4 結束語

針對傳統基于變密度法的求解方法的單向懲罰特性缺陷,筆者提出了一種求解連續體結構拓撲優化問題的新方法,并結合算例驗證了該方法對二維及三維設計域拓撲優化的可行性。

研究結論如下:

(1)與以往采用過濾技術以消除拓撲優化過程中存在的數值不穩定現象的技術相比,在不采用過濾技術的情況下,該新方法既可消除數值不穩定現象,又可得到具有清晰邊界的拓撲結構;

(2)在二維設計域拓撲優化迭代過程中,與傳統變密度拓撲優化求解方法相比,該新方法所得的最終結構拓撲邊界更加清晰,初始目標函數值、迭代步數和最終目標函數值更小;且與第二類辛格插值法和SIMP插值法相比,第一類辛格插值法更加敏捷,可得到更小的目標函數值;在三維設計域拓撲優化迭代過程中,該新方法所得拓撲結構可靠、清晰且不失真。

在后續的研究工作中,筆者將進一步開展該新方法與3D打印技術有效結合的關鍵技術研究,并將其應用于工程實踐。