基于變密度法的多材料與結(jié)構(gòu)一體拓?fù)鋬?yōu)化研究*

董 莉，吳曉明

(廈門(mén)大學(xué) 機(jī)電工程系，福建 廈門(mén) 361102)

0 引 言

結(jié)構(gòu)拓?fù)鋬?yōu)化是一種結(jié)構(gòu)概念設(shè)計(jì)方法，它可以在給定約束條件的情況下，根據(jù)優(yōu)化目標(biāo)，通過(guò)優(yōu)化算法計(jì)算得到結(jié)構(gòu)的材料最優(yōu)布局。拓?fù)鋬?yōu)化也是一種非常有效的結(jié)構(gòu)減重優(yōu)化方法，它可以得到比傳統(tǒng)的形狀優(yōu)化和尺寸優(yōu)化更輕的結(jié)構(gòu)，因此成為工程界研究熱點(diǎn)。

自從1988年BENDS?E等[1]提出了求解連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的均勻化方法以來(lái)，針對(duì)拓?fù)鋬?yōu)化的問(wèn)題得到了廣泛的研究。依據(jù)描述方式的不同, 目前常見(jiàn)的拓?fù)鋬?yōu)化方法包括均勻化方法、變密度法[2-3]、漸進(jìn)結(jié)構(gòu)優(yōu)化方法[4]，以及水平集方法[5]等。其中，變密度法是以單元密度作為設(shè)計(jì)變量，具有概念簡(jiǎn)單、設(shè)計(jì)變量少、收斂速度快、程序?qū)崿F(xiàn)簡(jiǎn)單的優(yōu)點(diǎn)。

較為成熟的研究大多是單一材料結(jié)構(gòu)拓?fù)鋬?yōu)化，然而單一材料結(jié)構(gòu)在工程實(shí)際中通常難以滿足一些需要特定綜合性能的工程結(jié)構(gòu)問(wèn)題，因此進(jìn)行多材料的結(jié)構(gòu)設(shè)計(jì)研究非常必要，梯度材料以及在梯度材料上的結(jié)構(gòu)設(shè)計(jì)優(yōu)化受到工程界廣泛關(guān)注。

多材料結(jié)構(gòu)拓?fù)鋬?yōu)化的研究始于THOMSEN[6]，之后一些學(xué)者用多種不同性能的材料和空洞進(jìn)行了多相材料拓?fù)鋬?yōu)化設(shè)計(jì)。SIGMUND和TORQUATO[7-8]針對(duì)兩相實(shí)體材料和空洞組成的柔性機(jī)構(gòu)拓?fù)?，以及具有極端熱膨脹特性的微結(jié)構(gòu)拓?fù)鋬?yōu)化進(jìn)行了研究；YIN等[9]用單個(gè)變量描述了多相材料條件下的材料屬性，提出了峰值函數(shù)多相材料插值模型；GAO等[10]采用線性對(duì)等混合材料插值模型，進(jìn)行了多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)；ZUO等[11]采用序列冪函數(shù)插值模型進(jìn)行了多相材料連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化。

另外有一些學(xué)者對(duì)功能梯度材料的拓?fù)鋬?yōu)化進(jìn)行了研究。陸丹和劉毅[12]采用雙向漸進(jìn)拓?fù)鋬?yōu)化(BESO)方法,引入了許用應(yīng)力,以實(shí)際應(yīng)力與許用應(yīng)力之比作為單元增刪的標(biāo)準(zhǔn)，結(jié)合梯度單元法實(shí)現(xiàn)了對(duì)功能梯度材料的拓?fù)鋬?yōu)化；XIA和WANG[13]基于水平集方法，以體分比和結(jié)構(gòu)邊界為設(shè)計(jì)變量,實(shí)現(xiàn)了功能梯度結(jié)構(gòu)的材料性能和拓?fù)洳季值牟⑿袃?yōu)化設(shè)計(jì)；邱克鵬和張衛(wèi)紅[14]采用了凸規(guī)劃求解策略以及周長(zhǎng)控制方法,實(shí)現(xiàn)了功能梯度MBB梁和功能梯度夾層結(jié)構(gòu)中夾芯的拓?fù)錁?gòu)型設(shè)計(jì)拓?fù)鋬?yōu)化。

這一類研究是在給定的功能梯度材料之上進(jìn)行的結(jié)構(gòu)設(shè)計(jì)，而結(jié)構(gòu)需要依據(jù)材料的狀況而定?，F(xiàn)代3D打印的多材料融合技術(shù)可以讓兩種甚至多種材料按照需求，以任意的配比進(jìn)行融合，且材料可以依據(jù)結(jié)構(gòu)的功能要求而變化，從而打破了結(jié)構(gòu)對(duì)材料的局限。因此，研究多材料與結(jié)構(gòu)一體設(shè)計(jì)優(yōu)化方法具有重要意義。

基于變密度法，筆者提出一種雙材料與結(jié)構(gòu)一體拓?fù)鋬?yōu)化模型。

1 材料與結(jié)構(gòu)一體優(yōu)化模型

1.1 變密度法拓?fù)鋬?yōu)化模型

變密度法將材料密度值設(shè)置為在區(qū)間[0,1]之間連續(xù)變化的變量，其中，密度值為0時(shí)表示單元材料不存在，密度值為1時(shí)表示單元材料存在。

假設(shè)材料密度與材料彈性模量之間的關(guān)系如下：

E(xi)=Emin+xi(E0-Emin)xi∈[0,1]

(1)

式中：xi—第i個(gè)單元的密度；E(xi)—第i個(gè)單元的材料彈性模量；E0—材料存在時(shí)的彈性模量；Emin—材料不存在時(shí)的彈性模量，為防止出現(xiàn)剛度陣奇異而無(wú)法求解，可令Emin=E0/1 000。

由于在結(jié)構(gòu)中出現(xiàn)的中間密度單元是不存在的，為了解決中間密度值單元，筆者引入懲罰插值函數(shù)來(lái)解決中間單元問(wèn)題。SIMP法是一個(gè)冪函數(shù)作為懲罰函數(shù)，使得所有單元的密度值懲罰到0或者1附近。

SIMP法模型的表達(dá)式如下：

(2)

式中：p—懲罰因子；n—結(jié)構(gòu)單元總數(shù)。

引入SIMP插值模型后的材料彈性模量如下：

(3)

引入SIMP法插值模型后的結(jié)構(gòu)的剛度矩陣、柔度函數(shù)以及敏度函數(shù)如下所示：

(4)

式中：K—整體剛度矩陣；ki—單元?jiǎng)偠染仃?；C(x)—結(jié)構(gòu)柔度值；F—載荷向量；U—位移向量；C′(x)—結(jié)構(gòu)靈敏度。

以體積約束下柔度最小化(即剛度最大化)建立基于SIMP插值模型的變密度法拓?fù)鋬?yōu)化數(shù)學(xué)模型如下所示：

(5)

式中：Ω—結(jié)構(gòu)初始設(shè)計(jì)區(qū)域；n—離散后的單元總數(shù)目；xi—設(shè)計(jì)變量，表示離散后每個(gè)單元的密度值；ui—單元位移向量；vi—每個(gè)單元初始體積；V*—結(jié)構(gòu)總體體積約束。

為避免求解時(shí)出現(xiàn)總剛度矩陣奇異，取xmin=0.001，xi取[xmin，1]之間的連續(xù)值。

1.2 材料與結(jié)構(gòu)一體優(yōu)化模型

本文提出一種基于變密度法的雙材料與結(jié)構(gòu)一體拓?fù)鋬?yōu)化模型。結(jié)構(gòu)的每一個(gè)有限元單元可由A，B兩種材料以任意比例融合而成，材料A的彈性模量為Ea，材料B的彈性模量為Eb，設(shè)單元i中材料B的比例為αi，則單元i中材料A的比例為1-αi，則單元i的材料為Bαi+A(1-αi)。

基于SIMP插值模型的結(jié)構(gòu)與材料一體彈性模量如下：

(6)

剛度矩陣和柔度函數(shù)為：

(7)

(8)

結(jié)構(gòu)柔度最小化作為目標(biāo)函數(shù)，以結(jié)構(gòu)的體積和材料B在總體結(jié)構(gòu)中所占比例作為約束，構(gòu)建基于SIMP法插值模型的結(jié)構(gòu)與材料一體拓?fù)鋬?yōu)化模型如下式所示：

(9)

式中：Ω—結(jié)構(gòu)初始設(shè)計(jì)區(qū)域；xi—離散后每個(gè)單元的密度值；αi—單元中材料B的占有比例；n—結(jié)構(gòu)離散后的單元總數(shù)；C—結(jié)構(gòu)柔度；F—外力向量；U—位移；Eα—材料A的彈性模量；Eb—材料B的彈性模量；K—結(jié)構(gòu)剛度陣；V*—結(jié)構(gòu)目標(biāo)體積；μ—拓?fù)渫瓿山Y(jié)構(gòu)中材料B總的使用比例。

為了避免計(jì)算過(guò)程中剛度矩陣奇異，用Emin、xmin代替0；取Emin=0.001，xmin=0.001。

2 敏度分析與過(guò)濾

2.1 敏度分析

結(jié)構(gòu)的體積對(duì)單元密度的靈敏度為：

(10)

結(jié)構(gòu)的柔度對(duì)單元密度的靈敏度為：

(11)

結(jié)構(gòu)的柔度對(duì)單元中材料B的比例靈敏度為：

(12)

2.2 過(guò)濾技術(shù)及優(yōu)化算法

為避免結(jié)果出現(xiàn)棋盤(pán)格現(xiàn)象和網(wǎng)格依賴性問(wèn)題，筆者采用敏度過(guò)濾技術(shù)，即將目標(biāo)單元和鄰接單元進(jìn)行加權(quán)平均，將得到的敏度值代替原來(lái)的值，從而實(shí)現(xiàn)過(guò)濾。

結(jié)構(gòu)的柔度對(duì)單元密度的靈敏度過(guò)濾如下：

(13)

(14)

式中：rmin—過(guò)濾半徑；dist(i,f)—目標(biāo)單元和鄰接單元的中心距。

筆者采用OC優(yōu)化算法求解式(9)優(yōu)化問(wèn)題。OC優(yōu)化算法推導(dǎo)過(guò)程簡(jiǎn)單、收斂速度快，是常用的優(yōu)化算法。通過(guò)引入拉格朗日乘數(shù)因子，構(gòu)造拉格朗日函數(shù)。根據(jù)當(dāng)設(shè)計(jì)變量取極值時(shí)拉格朗日函數(shù)應(yīng)滿足K-T必要條件，以及拉格朗日函數(shù)關(guān)于設(shè)計(jì)變量和乘數(shù)因子的駐值條件，可得到設(shè)計(jì)變量迭代格式(中間推導(dǎo)過(guò)程不再贅述)：

(15)

(16)

式中：m—移動(dòng)限制；η—阻尼系數(shù)，可以在保證收斂穩(wěn)定性的同時(shí)提升迭代效率。

Bi，Di如下式：

(17)

(18)

式中：λ1，λ2—拉格朗日乘子，可由二分法計(jì)算得到。

迭代終止條件為兩次迭代之間設(shè)計(jì)變量差值的絕對(duì)值小于等于一個(gè)值change，即定義終止條件為max(|xk+1-xk|)≤change，其中：k—迭代次數(shù)。

基于OC法的多材料與結(jié)構(gòu)一體拓?fù)鋬?yōu)化方法流程如圖1所示。

圖1 基于OC法的多材料與結(jié)構(gòu)一體拓?fù)鋬?yōu)化方法流程圖

3 算例及結(jié)果分析

筆者選取兩個(gè)算例進(jìn)行數(shù)值運(yùn)算，驗(yàn)證提出設(shè)計(jì)方法的有效性。兩個(gè)算例均使用筆者所提出新的模型計(jì)算，進(jìn)行結(jié)構(gòu)與材料一體拓?fù)鋬?yōu)化。

將結(jié)構(gòu)設(shè)計(jì)區(qū)域離散為四節(jié)點(diǎn)平面應(yīng)力單元，每個(gè)單元的相對(duì)密度值為xmin至1之間的連續(xù)值，每個(gè)單元由A、B兩種材料構(gòu)成；單元中B材料的比例為αi，A材料的比例則為1-αi，αi—0到1之間的連續(xù)變量。以每個(gè)單元的相對(duì)密度xi和單元中材料B的比例αi為設(shè)計(jì)變量，以體積和拓?fù)渫瓿山Y(jié)構(gòu)中材料B總的使用比例μ為約束，運(yùn)用OC優(yōu)化算法進(jìn)行求解。

3.1 懸臂梁

懸臂梁結(jié)構(gòu)設(shè)計(jì)域及邊界條件如圖2所示。

圖2 懸臂梁結(jié)構(gòu)設(shè)計(jì)域及邊界條件

圖2中，設(shè)計(jì)域尺寸為40×20，懸臂梁結(jié)構(gòu)左端固定，受力作用點(diǎn)集中在右側(cè)中間的一個(gè)單元節(jié)點(diǎn)上，初始結(jié)構(gòu)設(shè)計(jì)域離散成40×20個(gè)四節(jié)點(diǎn)平面應(yīng)力單元。材料A的彈性模量Ea=0.5，材料B的彈性模量Eb=1。體積分?jǐn)?shù)為50%，過(guò)濾半徑為1.5，懲罰因子為3，μ為0.4，移動(dòng)限制為0.2；阻尼因子為1/2，設(shè)計(jì)變量初始值x為0.5，α為0.4。

根據(jù)筆者提出結(jié)構(gòu)與材料拓?fù)鋬?yōu)化算法，得出的懸臂梁優(yōu)化結(jié)果如圖3所示。

圖3 懸臂梁優(yōu)化結(jié)果

圖3(a)為計(jì)算所得的拓?fù)浣Y(jié)構(gòu)形狀，顏色為黑色的單元為保留單元，白色單元為去掉單元，灰色單元為密度值處于中間的單元，該類單元可根據(jù)要求保留或刪去；

圖3(b)為計(jì)算所得材料B在拓?fù)湫螤罨A(chǔ)上每個(gè)單元中的比例值αi的示意圖，顏色為黑色的單元值αi為1，代表該類單元由B材料構(gòu)成；顏色為灰色的單元由A、B兩種材料按照單元中B材料的比例為αi，A材料的比例為1-αi混合構(gòu)成。

在圖3(a)中保留單元基礎(chǔ)上的圖3(b)白色單元值為0，代表此類單元由A材料構(gòu)成。

迭代終止條件為change=0.001時(shí)，計(jì)算歷時(shí)20.687 825 s，迭代339步，柔度值為82.817 3。

由優(yōu)化結(jié)果可以看出，筆者提出的優(yōu)化方法所得拓?fù)浣Y(jié)構(gòu)形狀邊界清晰、分布合理，整體優(yōu)化穩(wěn)定。在材料受力單元附近與結(jié)構(gòu)上下邊界受力變形大的部分，都采用了彈性模量較大的B材料，表明得到的材料比例分布也是十分合理的。

3.2 MBB梁

MBB梁結(jié)構(gòu)設(shè)計(jì)域及邊界條件如圖4所示。

圖4 MBB梁結(jié)構(gòu)設(shè)計(jì)域及邊界條件

圖4中，設(shè)計(jì)域尺寸為40×20，懸臂梁結(jié)構(gòu)左端固定，受力作用點(diǎn)集中在一個(gè)單元節(jié)點(diǎn)上，初始結(jié)構(gòu)設(shè)計(jì)域離散成40×20個(gè)四節(jié)點(diǎn)平面應(yīng)力單元。

材料A的彈性模Ea=0.5，材料B的彈性模量Eb=1；體積分?jǐn)?shù)為50%，過(guò)濾半徑為1.5，懲罰函數(shù)為3，μ為0.4，移動(dòng)限制為0.2，阻尼因子為1/2，設(shè)計(jì)變量初始值x為0.5，α為0.4。

根據(jù)筆者提出的結(jié)構(gòu)與材料拓?fù)鋬?yōu)化算法，得出的MBB梁優(yōu)化結(jié)果如圖5所示。

圖5 MBB梁優(yōu)化結(jié)果

圖5(a)為計(jì)算所得的拓?fù)浣Y(jié)構(gòu)形狀，黑色的單元為保留單元，白色單元為去掉單元，灰色單元為密度值處于中間的單元，此類單元可根據(jù)要求保留或刪去；

圖5(b)為計(jì)算所得在拓?fù)浣Y(jié)構(gòu)形狀的基礎(chǔ)上B材料比例值，黑色單元由B材料構(gòu)成，灰色單元由A、B兩種材料按照單元中B材料的比例為αi，A材料的比例為1-αi混合構(gòu)成，在圖5(a)中保留單元基礎(chǔ)上的圖5(b)白色單元由A材料構(gòu)成。

迭代終止條件為change=0.001時(shí)，計(jì)算歷時(shí)16.914 929 s，迭代266步，柔度值為93.509 0。

由優(yōu)化結(jié)果可以看出，筆者提出的優(yōu)化方法所得拓?fù)浣Y(jié)構(gòu)形狀邊界清晰、分布合理，整體優(yōu)化穩(wěn)定；在材料受力單元附近與結(jié)構(gòu)下邊界受力變形大的部分，都采用了彈性模量較大的B材料，表明得到的材料比例分布也是十分合理的。

4 結(jié)束語(yǔ)

基于SIMP法插值模型，筆者建立了單元中兩種材料不同比例與單元彈性模量之間的關(guān)系，在變密度法中，以單元密度作為設(shè)計(jì)變量的基礎(chǔ)上，加入材料在單元構(gòu)成中所占的比例為設(shè)計(jì)變量，以結(jié)構(gòu)柔度最小化為目標(biāo)，結(jié)構(gòu)體積和其中一類材料在結(jié)構(gòu)中所占比值為約束，建立了基于變密度法的結(jié)構(gòu)與材料一體拓?fù)鋬?yōu)化方法模型，并采用OC優(yōu)化算法進(jìn)行了求解。

相關(guān)算例及研究結(jié)果表明，該方法可以計(jì)算出結(jié)構(gòu)的拓?fù)錁?gòu)型，同時(shí)可得到兩種材料比例數(shù)值結(jié)果；且拓?fù)錁?gòu)型邊界清晰，可以實(shí)現(xiàn)結(jié)構(gòu)與材料一體的協(xié)同優(yōu)化。