占金青，彭怡平，龍良明，劉 敏

(華東交通大學(xué) 機(jī)電與車輛工程學(xué)院，江西 南昌 330013)

0 引言

拓?fù)鋬?yōu)化是一種非常有效的結(jié)構(gòu)概念設(shè)計(jì)方法，通過在指定的載荷作用、邊界條件及約束條件下尋求宏觀結(jié)構(gòu)或材料微結(jié)構(gòu)的材料最佳布局，從而使其某種性能指標(biāo)達(dá)到最優(yōu)[1-2]?，F(xiàn)有研究大多集中在單一材料結(jié)構(gòu)拓?fù)鋬?yōu)化問題，但由于單一材料結(jié)構(gòu)通常難以實(shí)現(xiàn)特定的綜合性能要求，且實(shí)際工程結(jié)構(gòu)通常由多相材料組成，使得進(jìn)行多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)非常必要[3-4]。

Thomsen等[5]首先進(jìn)行了多相材料的拓?fù)鋬?yōu)化研究;Sigmund等[6]采用均勻化方法進(jìn)行負(fù)熱膨脹系數(shù)材料的微結(jié)構(gòu)多相材料拓?fù)鋬?yōu)化設(shè)計(jì);Yin等[7]利用峰值函數(shù)建立了多相材料拓?fù)鋬?yōu)化設(shè)模型;Sigmund[8]提出了基于固體各向同性材懲罰模型的多相材料拓?fù)鋬?yōu)化設(shè)計(jì)方法;Gao等[9]采用線性對(duì)等混合材料插值模型，進(jìn)行了多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì);張憲民等[10]提出了基于平行策略的多相材料拓?fù)鋬?yōu)化設(shè)計(jì)方法;Zuo等[11]采用序列冪函數(shù)插值模型進(jìn)行了多相材料連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化;杜義賢等[12]采用多重網(wǎng)格方法進(jìn)行多相材料拓?fù)鋬?yōu)化設(shè)計(jì)，其計(jì)算效率有所提高;龍凱等[13]采用獨(dú)立連續(xù)映射方法進(jìn)行多相材料結(jié)構(gòu)動(dòng)態(tài)拓?fù)鋬?yōu)化設(shè)計(jì);Taheri等[14]采用等幾何方法進(jìn)行多相材料的功能梯度結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì);Zhang等[15]提出了基于基結(jié)構(gòu)法的多相材材料結(jié)構(gòu)非線性拓?fù)鋬?yōu)化設(shè)計(jì)方法;Liu等[16]進(jìn)行了考慮尺寸約束的多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)研究。目前，多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化研究較多以各相材料體積或質(zhì)量份數(shù)為約束建立拓?fù)鋬?yōu)化模型;但僅以各相材料體積或質(zhì)量份數(shù)為約束的柔度優(yōu)化獲得多相材料拓?fù)浣Y(jié)構(gòu)難以滿足材料失效準(zhǔn)則。

目前，應(yīng)力約束拓?fù)鋬?yōu)化研究主要集中于單一材料結(jié)構(gòu)優(yōu)化設(shè)計(jì)。Le等[17]和Holmberg等[18]采用P范數(shù)對(duì)所有單元的應(yīng)力凝聚成一個(gè)全局應(yīng)力約束，進(jìn)行應(yīng)力約束結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì);王選等[19]采用改進(jìn)的雙向漸進(jìn)結(jié)構(gòu)優(yōu)化法進(jìn)行應(yīng)力約束拓?fù)鋬?yōu)化設(shè)計(jì)研究;占金青等[20]考慮熱載荷和機(jī)械載荷的共同作用，提出了基于全局應(yīng)力約束的熱彈性結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法;Amir等[21]提出了基于應(yīng)力約束的彈塑性材料結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)方法;Collet等[22]采用應(yīng)力約束進(jìn)行材料微結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì);Long等[23]采用二次規(guī)劃算法進(jìn)行考慮諧激勵(lì)作用下的應(yīng)力約束結(jié)構(gòu)拓?fù)鋬?yōu)化。最近，Guo等[24]通過引入應(yīng)力梯度函數(shù)來避免應(yīng)力引起的奇異性，采用水平集方法進(jìn)行多相材料結(jié)構(gòu)局部應(yīng)力約束拓?fù)鋬?yōu)化設(shè)計(jì)，并研究以基于應(yīng)力梯度放大的應(yīng)力優(yōu)化列式最小化為目標(biāo)的多相材料拓?fù)鋬?yōu)化設(shè)計(jì)方法;Chu等[25-26]采用應(yīng)力懲罰方法進(jìn)行多相材料結(jié)構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化設(shè)計(jì)，通過應(yīng)力懲罰項(xiàng)在優(yōu)化目標(biāo)中控制結(jié)構(gòu)局部應(yīng)力，該方法仍屬于局部約束，但是避免了大量的單元局部應(yīng)力約束。

本文提出一種基于可分離應(yīng)力插值模型的多相材料結(jié)構(gòu)全局應(yīng)力約束拓?fù)鋬?yōu)化設(shè)計(jì)方法。采用可分離應(yīng)力插值模型計(jì)算多相材料結(jié)構(gòu)的剛度和應(yīng)力，以結(jié)構(gòu)體積最小化作為目標(biāo)函數(shù);采用應(yīng)力松弛方法來消除應(yīng)力奇異解現(xiàn)象;利用改進(jìn)的P范數(shù)求解各相材料結(jié)構(gòu)的最大應(yīng)力，并以各相材料結(jié)構(gòu)的最大應(yīng)力作為約束，建立多相材料結(jié)構(gòu)全局應(yīng)力約束拓?fù)鋬?yōu)化模型，進(jìn)一步采用移動(dòng)漸近線優(yōu)化算法進(jìn)行優(yōu)化問題求解。

1 多相材料結(jié)構(gòu)全局應(yīng)力約束拓?fù)鋬?yōu)化模型

1.1 可分離應(yīng)力插值模型

傳統(tǒng)的固體各項(xiàng)同性材料懲罰(Solid Isotropic Material with Penalization, SIMP)模型中，每個(gè)單元只能對(duì)應(yīng)一個(gè)單元?jiǎng)偠群鸵粋€(gè)單元應(yīng)力，不能對(duì)單獨(dú)各相材料結(jié)構(gòu)的應(yīng)力進(jìn)行計(jì)算。采用基于單元堆積方法的可分離應(yīng)力插值模型[26-27]，在單元應(yīng)力求解過程中，每個(gè)單元能夠?qū)?yīng)于每相材料都有一個(gè)應(yīng)力，應(yīng)力的個(gè)數(shù)與材料的數(shù)目有關(guān)，如圖1所示。多相材結(jié)構(gòu)單元的剛度矩陣可以表示為

(1)

(2)

(3)

式中Ei和υi分別為第i相材料的彈性模量和泊松比。

1.2 設(shè)計(jì)變量過濾方法

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化容易出現(xiàn)網(wǎng)格依賴、棋盤格等現(xiàn)象，采用設(shè)計(jì)變量過濾能夠避免這些數(shù)值的不穩(wěn)定性，獲得清晰的拓?fù)浣Y(jié)構(gòu)[28]。設(shè)計(jì)變量xe,i通過相鄰設(shè)計(jì)變量xj,i的加權(quán)平均，獲得對(duì)應(yīng)于第i相材料的過濾單元密度ρe,i。過濾單元密度ρe,i也稱為物質(zhì)密度，用來計(jì)算單元的剛度和應(yīng)力，

(4)

式中：Ωe表示與單元e質(zhì)心距離小于過濾半徑rmin的所有單元集合，如圖2所示；wj為權(quán)重系數(shù)，

(5)

式中：rj為單元j和單元e的質(zhì)心之間的距離;對(duì)于不在子區(qū)域Ωe的權(quán)重系數(shù)wj=0。因此，將單元權(quán)重系數(shù)wj由原子區(qū)域單元擴(kuò)展至所有單元，設(shè)計(jì)變量過濾可改寫為

(6)

式中:NE為有限元單元的數(shù)目;Wej為任一單元的權(quán)重系數(shù)。

1.3 拓?fù)鋬?yōu)化模型

采用可分離應(yīng)力插值模型求解結(jié)構(gòu)的單元?jiǎng)偠燃皯?yīng)力，以結(jié)構(gòu)的體積最小化作為目標(biāo)函數(shù)，以各相材料的最大應(yīng)力和其中NM-1相材料的體積作為約束，建立基于應(yīng)力約束的多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化模型，其拓?fù)鋬?yōu)化模型為

s.t.

KU=F;

(7)

1.4 全局應(yīng)力約束

單元e對(duì)應(yīng)于第i相材料的Von Mises等效應(yīng)力可以表示為

=[(σe,i)TVσe,i]1/2。

(8)

(9)

(10)

對(duì)于多相材料結(jié)構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化問題，結(jié)構(gòu)設(shè)計(jì)域被離散為NE個(gè)單元，每個(gè)單元對(duì)應(yīng)多個(gè)應(yīng)力約束，這也意味著約束數(shù)目非常巨大，且約束靈敏度分析計(jì)算量非常龐大。采用改進(jìn)的P范數(shù)將所有的單元Von Mises應(yīng)力凝聚化為一個(gè)近似等于結(jié)構(gòu)最大Von Mises應(yīng)力的全局應(yīng)力，從而大大提高計(jì)算效率。第i相材料結(jié)構(gòu)的P范數(shù)應(yīng)力可表示為

(11)

式中P為范數(shù)參數(shù)。當(dāng)P趨近無窮大，P范數(shù)應(yīng)力等于最大應(yīng)力值。但是P不是越大越好，當(dāng)P越大，非線性越明顯，甚至將導(dǎo)致優(yōu)化問題求解困難，因此P的取值不能太大。文獻(xiàn)[17]采用P范數(shù)方法，是以單元的體積作為比例因子，具有內(nèi)置縮放的作用，有利于應(yīng)力約束優(yōu)化問題的收斂；但是采用國際單位制時(shí)，單元的體積特別小，比例因子內(nèi)置縮放作用消失。本文以單元密度ρe,i作為比例因子的改進(jìn)的P范數(shù)，采用國際單位制時(shí)也能夠具有很好的內(nèi)置縮放作用。

為了使得P范數(shù)應(yīng)力接近最大應(yīng)力，采用自適應(yīng)約束縮放方法[20]修改P范數(shù)應(yīng)力：

(12)

式中Ci為自適應(yīng)約束縮放系數(shù)，當(dāng)?shù)綌?shù)k≥1時(shí)

(13)

2 靈敏度分析與優(yōu)化算法

2.1 靈敏度分析

采用基于梯度的優(yōu)化算法求解應(yīng)力約束拓?fù)鋬?yōu)化問題，并進(jìn)行目標(biāo)函數(shù)和約束的靈敏度分析是必要的。

目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的靈敏度可直接求導(dǎo)得到：

(14)

對(duì)式(12)求導(dǎo)，可得含有第i相材料結(jié)構(gòu)的最大應(yīng)力約束對(duì)單元設(shè)計(jì)變量xe,i的靈敏度：

(15)

分別對(duì)式(15)各項(xiàng)進(jìn)行求解，則有

(16)

(17)

對(duì)式(8)求導(dǎo)有

(18)

式(18)可改寫為

(19)

(20)

由于作用載荷恒定，對(duì)有限元平衡方程兩邊同時(shí)對(duì)單元密度ρe,i求導(dǎo)，則有

(21)

由式(21)可得

(22)

將式(22)代入式(20)，則有

(23)

將式(23)和式(19)代入式(17)有

(24)

令

(25)

將式(24)和式(16)代入式(15)，則有

(26)

令有下列伴隨方程

(27)

由式(27)求出λe，式(26)可改寫為

(28)

其中

(29)

2.2 優(yōu)化算法

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化問題通常主要采用優(yōu)化準(zhǔn)則算法(Optimality Criteria method, OC)和數(shù)學(xué)規(guī)劃算法求解。OC算法[29]簡單明了，但是難以適用于多約束的多相材料結(jié)構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化問題。移動(dòng)漸近優(yōu)化算法(Method of Moving Asymptotes，MMA)[30]適用于復(fù)雜的多約束優(yōu)化問題求解，具有很好的魯棒性。本文采用MMA算法進(jìn)行多相材料結(jié)構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化問題求解。

3 數(shù)值算例

本章中，多相材料結(jié)構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化設(shè)計(jì)方法在MATLAB軟件中編程實(shí)現(xiàn)，通過兩個(gè)數(shù)值算例來驗(yàn)證其有效性。所有算例中，實(shí)體材料1和實(shí)體材料2的彈性模量分別為190 GPa和210 GPa，泊松比μ均為0.3，兩相材料的屈服應(yīng)力極限值分別為σ1=241 MPa和σ2=358 MPa。設(shè)計(jì)變量過濾的半徑rmin=2.5，范數(shù)參數(shù)P=8。

3.1 L型梁結(jié)構(gòu)

L型梁結(jié)構(gòu)的設(shè)計(jì)域、作用載荷及邊界條件如圖3所示。L型梁結(jié)構(gòu)尺寸為100 mm×100 mm，厚度為1 mm，固定梁結(jié)構(gòu)的左上端，施加作用荷載F大小為250 N，均勻分布在結(jié)構(gòu)右上端的6節(jié)點(diǎn)上，如圖3的局部放大所顯示，以避免載荷作用區(qū)域應(yīng)力集中。材料2的體積份數(shù)約束β2=0.15，初始結(jié)構(gòu)域離散成6 400個(gè)四節(jié)點(diǎn)平面應(yīng)力單元。

為驗(yàn)證所提方法的有效性，對(duì)L型梁進(jìn)行有、無應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)。無應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化采用體積約束的結(jié)構(gòu)柔度最小化模型，各相材料體積份數(shù)約束的大小根據(jù)應(yīng)力約束拓?fù)鋬?yōu)化結(jié)果得到，其他設(shè)計(jì)參數(shù)與應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化相同。有無應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化的L型梁結(jié)構(gòu)的拓?fù)鋱D及應(yīng)力分布圖分別如圖4和圖5所示。拓?fù)浣Y(jié)構(gòu)中紅色區(qū)域代表實(shí)體材料1(彈性模量小)，藍(lán)色區(qū)域代表實(shí)體材料2(彈性模量大)，白色區(qū)域代表空洞。

由圖4和表1可知，無應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化獲得的L型梁拓?fù)浣Y(jié)構(gòu)位于內(nèi)凹角拐角出現(xiàn)明顯應(yīng)力集現(xiàn)象，材料1和材料2結(jié)構(gòu)的最大應(yīng)力分別為305.719 MPa和428.738 MPa，均超過材料1和材料2的屈服極限值。然而，應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化獲得的拓?fù)浣Y(jié)構(gòu)在拐角處出現(xiàn)了過渡弧形邊界(如圖5)，并且結(jié)構(gòu)的應(yīng)力分布更加均勻，有效地抑制了應(yīng)力集中現(xiàn)象。由圖4b和圖5b可見，有些不同相材料連接區(qū)域存在應(yīng)力不連續(xù)現(xiàn)象，這主要是位于不同相材料連接區(qū)域出現(xiàn)兩種相材料堆積的單元，且存在兩材料堆積的單元在應(yīng)力分布圖中只顯示應(yīng)力較大的相材料應(yīng)力分布，從而導(dǎo)致有些不同材料連接區(qū)域存在應(yīng)力不連續(xù)現(xiàn)象。由表1可知，在應(yīng)力約束拓?fù)鋬?yōu)化獲得拓?fù)浣Y(jié)構(gòu)中，材料1結(jié)構(gòu)和材料2結(jié)構(gòu)的最大應(yīng)力分別240.669 MPa和358.008 MPa，均能很好地滿足應(yīng)力約束極限值，滿足結(jié)構(gòu)的靜強(qiáng)度要求；與無應(yīng)力約束拓?fù)鋬?yōu)化相比，應(yīng)力約束多相材料拓?fù)鋬?yōu)化獲得L型梁結(jié)構(gòu)的柔度較大，但是其拓?fù)浣Y(jié)構(gòu)的強(qiáng)度更好，結(jié)構(gòu)應(yīng)力分布更加均勻。

表1 不同優(yōu)化模型獲得的L梁拓?fù)鋬?yōu)化結(jié)果

3.2 懸臂梁結(jié)構(gòu)

懸臂梁結(jié)構(gòu)的設(shè)計(jì)域、作用載荷及邊界條件如圖6所示。懸臂梁結(jié)構(gòu)尺寸為200 mm×100 mm，厚度為1 mm；結(jié)構(gòu)的左端固定，均布載荷施加結(jié)構(gòu)的右端中部，載荷的大小為800 N，均勻分布5單元節(jié)點(diǎn)上。材料2的體積份數(shù)約束β2=0.25。初始結(jié)構(gòu)域離散成20 000個(gè)四節(jié)點(diǎn)平面應(yīng)力單元。

有無應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化的懸臂型梁結(jié)構(gòu)拓?fù)鋱D及應(yīng)力分布圖分別如圖7和圖8所示。同樣，由圖7和表2可知，無應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化獲得的懸臂梁拓?fù)浣Y(jié)構(gòu)中，材料1和材料2結(jié)構(gòu)的最大應(yīng)力分別為306.312 MPa和490.421 MPa，也均超過材料1和材料2的屈服極限值。然而，應(yīng)力約束多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化獲得的懸臂梁拓?fù)浣Y(jié)構(gòu)有很大不同，邊界形狀突變得到抑制(如圖8)，且結(jié)構(gòu)的應(yīng)力分布更加均勻，有效地抑制了應(yīng)力集中現(xiàn)象。由表2可知，在應(yīng)力約束拓?fù)鋬?yōu)化獲得優(yōu)化結(jié)構(gòu)中，材料1結(jié)構(gòu)和材料2結(jié)構(gòu)的最大應(yīng)力分別為240.719 MPa和345.430 MPa，均沒有超過材料屈服極限值，滿足結(jié)構(gòu)的靜強(qiáng)度要求。同樣，與無應(yīng)力約束拓?fù)鋬?yōu)化相比，應(yīng)力約束多相材料拓?fù)鋬?yōu)化獲得優(yōu)懸臂梁結(jié)構(gòu)的柔度較大，但是其拓?fù)浣Y(jié)構(gòu)的強(qiáng)度更好，結(jié)構(gòu)應(yīng)力分布更加均勻。

表2 不同優(yōu)化模型獲得的懸臂梁拓?fù)鋬?yōu)化結(jié)果

4 結(jié)束語

本文基于可分離應(yīng)力插值模型提出了多相材料結(jié)構(gòu)全局應(yīng)力約束拓?fù)鋬?yōu)化設(shè)計(jì)方法。以結(jié)構(gòu)體積最小化作為目標(biāo)函數(shù)，將可分離應(yīng)力插值模型用于計(jì)算多相材料結(jié)構(gòu)的剛度和應(yīng)力，實(shí)現(xiàn)計(jì)算每相材料結(jié)構(gòu)的應(yīng)力，采用改進(jìn)的P范數(shù)求解結(jié)構(gòu)的最大應(yīng)力，并以各相材料結(jié)構(gòu)的最大應(yīng)力作為約束，建立了多相材料結(jié)構(gòu)全局應(yīng)力約束拓?fù)鋬?yōu)化模型，進(jìn)行優(yōu)化問題求解，通過數(shù)值算例驗(yàn)證了其可行性。

算例結(jié)果表明，與無應(yīng)力約束比較，應(yīng)力約束的多相材料結(jié)構(gòu)拓?fù)鋬?yōu)化獲得的拓?fù)浣Y(jié)構(gòu)有所不同，能很好地滿足各相材料應(yīng)力約束，滿足強(qiáng)度失效準(zhǔn)則；能夠抑制應(yīng)力集中現(xiàn)象，使得多相材料結(jié)構(gòu)應(yīng)力分布更加均勻。本文是在確定性載荷條件下進(jìn)行的研究工作，未來將考慮載荷不確定性因素對(duì)多相材料結(jié)構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化結(jié)果的影響。