基于辛格插值模式的連續(xù)體結(jié)構(gòu)拓撲優(yōu)化研究*

李玉剛,楊 林

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

0 引 言

拓撲優(yōu)化是指在設(shè)計領(lǐng)域內(nèi)尋求最佳的材料分布,從而使所得到的結(jié)構(gòu)達到最佳的或規(guī)定的性能。

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

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

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

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

1 辛格插值模式

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

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

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

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

利用辛格函數(shù)[12](即sinc函數(shù))在[0,1]區(qū)間內(nèi)具有較高斜率的單調(diào)特性(可使中間密度單元快速向?qū)嶓w單元進化或快速向空洞單元演化),構(gòu)建如下材料屬性的辛格插值模式:

(1)

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

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

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

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

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

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

對比圖(1，2)可知:

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

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

相應(yīng)地,辛格插值函數(shù)與單元剛度矩陣和彈性模量之間的關(guān)系可表示為:

(2)

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

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

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

(3)

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

其表達式如下:

(4)

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

所以,在拓撲優(yōu)化過程中,應(yīng)根據(jù)單元密度值大小不斷對迭代方向進行調(diào)整,這樣即使不采用過濾技術(shù),也能得到清晰邊界的拓撲結(jié)構(gòu)。

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

(5)

2 基于辛格插值模式的優(yōu)化模型及算法

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

(6)

式中:f—結(jié)構(gòu)總體柔度函數(shù);ki—第i單元優(yōu)化后的剛度矩陣;ui—第i單元優(yōu)化后的位移向量;V—優(yōu)化后結(jié)構(gòu)總體積;Vi—單元體積;V*—目標優(yōu)化體積;F—作用力;K—總剛度矩陣;U—總位移矩陣;

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

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

(7)

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

(8)

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

(9)

為方便表達,此處令:

(10)

(11)

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

(12)

其中:

(13)

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

(14)

其中:

(15)

其迭代終止條件如下:

(16)

其中:取δ=0.01。

3 算 例

3.1 算例一

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

為簡化計算,該例中所有設(shè)計參數(shù)均按無量綱處理,且不會影響優(yōu)化有效性。

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

圖3 二維設(shè)計域及有限元模型

數(shù)值不穩(wěn)定現(xiàn)象[15]是影響連續(xù)體結(jié)構(gòu)拓撲優(yōu)化結(jié)果質(zhì)量的關(guān)鍵因素之一。

在傳統(tǒng)變密度方法拓撲優(yōu)化過程中,通常采用空間敏度過濾技術(shù)以消除優(yōu)化過程中存在的數(shù)值不穩(wěn)定現(xiàn)象,但該技術(shù)易使敏度在其過濾半徑內(nèi)均勻化,即最終拓撲結(jié)構(gòu)的拓撲邊界易出現(xiàn)模糊現(xiàn)象。

不采用過濾技術(shù)時,不同τ值的拓撲結(jié)果如圖4所示。

圖4 不采用過濾技術(shù)時,不同τ值的拓撲結(jié)果

根據(jù)圖4,不采用過濾技術(shù)時,辛格插值模式所得優(yōu)化拓撲結(jié)果結(jié)構(gòu)邊界清晰,即未出現(xiàn)網(wǎng)格依賴性、棋盤格等現(xiàn)象,且τ取不同值時的拓撲結(jié)構(gòu)基本相同,所以雙向促進因子的取值對拓撲結(jié)果影響不大。

相應(yīng)的,τ取不同值時的優(yōu)化迭代過程如圖5所示。

圖5 τ取不同值時的優(yōu)化迭代過程

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

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

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

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

圖6 同時采用敏度過濾技術(shù)時,不同τ值的拓撲結(jié)果

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

3.2 算例二

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

為了簡化計算過程,上述設(shè)計參數(shù)按無量綱處理,求具有最小體積和最大剛度的結(jié)構(gòu)最優(yōu)拓撲,設(shè)初始單元密度為0.9。

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

圖7 二維設(shè)計域及有限元模型

采用第一類辛格插值法所得最終拓撲結(jié)構(gòu)如圖8所示。

圖8 第一類辛格插值法所得最終拓撲結(jié)構(gòu)

采用第二類辛格插值法所得最終拓撲結(jié)構(gòu)如圖9所示。

圖9 第二類辛格插值法所得最終拓撲結(jié)構(gòu)

采用SIMP法所得最終拓撲結(jié)構(gòu)如圖10所示。

圖10 SIMP法所得最終拓撲結(jié)構(gòu)

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

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

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

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

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

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

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

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

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

3.3 算例三

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

為簡化計算,該算例中所有設(shè)計參數(shù)均按無量綱處理,且不會影響優(yōu)化的有效性。

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

圖13 三維設(shè)計域及有限元模型

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

τ取不同值時的最終拓撲結(jié)構(gòu)如圖14所示。

圖14 τ取不同值時的最終拓撲結(jié)構(gòu)

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

相應(yīng)的,τ取不同值時的優(yōu)化迭代過程如圖15所示。

圖15 τ取不同值時的優(yōu)化迭代過程

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

4 結(jié)束語

針對傳統(tǒng)基于變密度法的求解方法的單向懲罰特性缺陷,筆者提出了一種求解連續(xù)體結(jié)構(gòu)拓撲優(yōu)化問題的新方法,并結(jié)合算例驗證了該方法對二維及三維設(shè)計域拓撲優(yōu)化的可行性。

研究結(jié)論如下:

(1)與以往采用過濾技術(shù)以消除拓撲優(yōu)化過程中存在的數(shù)值不穩(wěn)定現(xiàn)象的技術(shù)相比,在不采用過濾技術(shù)的情況下,該新方法既可消除數(shù)值不穩(wěn)定現(xiàn)象,又可得到具有清晰邊界的拓撲結(jié)構(gòu);

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

在后續(xù)的研究工作中,筆者將進一步開展該新方法與3D打印技術(shù)有效結(jié)合的關(guān)鍵技術(shù)研究,并將其應(yīng)用于工程實踐。