基于BESO方法的多孔微觀結(jié)構(gòu)材料優(yōu)化設(shè)計

閆曉磊，林逢春，花海燕，謝 露

（福建工程學(xué)院福建省汽車電子與電驅(qū)動重點(diǎn)實(shí)驗(yàn)室，福建 福州 350118）

1 引言

隨著我國工業(yè)的高速發(fā)展和科學(xué)技術(shù)的不斷進(jìn)步，各個領(lǐng)域?qū)Ω咝阅懿牧系男枨笕找嫫惹小Ｒ源藶楸尘埃嗫孜⒔Y(jié)構(gòu)材料因其具有相對密度小、強(qiáng)度高、重量輕等多種良好的物理性能和在生活中廣泛的應(yīng)用如蜂窩材料和泡沫金屬材料［1］，近年來成為關(guān)注的熱點(diǎn)［2］。由于多孔材料的宏觀性能取決于其內(nèi)部微觀結(jié)構(gòu)的形式。因此，可以通過設(shè)計其微結(jié)構(gòu)幾何構(gòu)型，來構(gòu)造特定的材料屬性。采用結(jié)構(gòu)拓?fù)鋬?yōu)化方法對材料微結(jié)構(gòu)進(jìn)行設(shè)計是實(shí)現(xiàn)特殊功能人工材料的一種有效途徑，也是近年來材料與工程領(lǐng)域的研究熱點(diǎn)［3］。一些學(xué)者已經(jīng)在這方面開展了十分有意義的工作，文獻(xiàn)［4］基于能量均勻化方法和水平集拓?fù)鋬?yōu)化方法設(shè)計具有最大體積模量和最大剪切模量的微結(jié)構(gòu)材料；文獻(xiàn)［5］基于獨(dú)立連續(xù)映射法（Independent Continuous Mapping，ICM）方法設(shè)計出具有良好導(dǎo)熱性能的微結(jié)構(gòu)材料；文獻(xiàn)［6］基于均勻化方法設(shè)計多孔材料，以剛度和負(fù)泊松比為優(yōu)化設(shè)計目標(biāo)，構(gòu)造多目標(biāo)拓?fù)鋬?yōu)化模型，利用改進(jìn)優(yōu)化準(zhǔn)則算法求解模型，得到兼具吸能和承載特性的微結(jié)構(gòu)材料。文獻(xiàn)［7］通過對具有負(fù)泊松比和負(fù)熱膨脹系數(shù)的微結(jié)構(gòu)進(jìn)行設(shè)計，實(shí)現(xiàn)了材料的多功能化；文獻(xiàn)［8］基于SIMP拓?fù)鋬?yōu)化方法設(shè)計具有負(fù)泊松比的微結(jié)構(gòu)材料；文獻(xiàn)［9］基于多相水平集模型的參數(shù)化水平集方法實(shí)現(xiàn)了微結(jié)構(gòu)形狀和拓?fù)涞难葑儯⑶以O(shè)計了具有極限熱彈性性能的微結(jié)構(gòu)材料；文獻(xiàn)［10］基于SIMP拓?fù)鋬?yōu)化方法設(shè)計了具有極限彈性性能的微結(jié)構(gòu)材料；陳文獻(xiàn)［11］針對多孔材料的高吸聲性能和具有可設(shè)計性特點(diǎn)，研究和設(shè)計了高聲控制材料與結(jié)構(gòu)，得到了一種具有較高聲能吸收率的多孔微結(jié)構(gòu)材料。

雙向漸進(jìn)結(jié)構(gòu)優(yōu)化（Bidirectional Evolutionary Structural Optimization，BESO）方法是一種常用的拓?fù)鋬?yōu)化方法，具有原理簡單、易于實(shí)現(xiàn)的特點(diǎn)［12］，但優(yōu)化結(jié)果存在鋸齒狀邊界問題［13］，無法直接用于精確幾何建模。針對此問題，這里對BESO方法進(jìn)行了改進(jìn)，然后基于改進(jìn)的BESO方法，對多孔周期微結(jié)構(gòu)材料進(jìn)行了優(yōu)化設(shè)計，得到了模量最大的多孔材料最優(yōu)微結(jié)構(gòu)拓?fù)錁?gòu)型。

2 均勻化方法及其有限元格式

均勻化理論是由文獻(xiàn)［14-15］在20世紀(jì)70年代中期提出來的，其基本原理是基于攝動技術(shù)，把原問題轉(zhuǎn)化為一個微觀均勻化問題與一個宏觀均勻化問題［16］。均勻化理論可以細(xì)致地考慮材料的微觀結(jié)構(gòu)，在微結(jié)構(gòu)材料的性能研究中發(fā)揮著重要的作用，如圖1所示。宏觀尺度（x1-x2）下的尺寸遠(yuǎn)大于微觀尺度（y1-y2）下的尺寸，則從宏觀上，可以將微結(jié)構(gòu)材料（圖1（b））可以看成均質(zhì)材料（圖1（a））；從微觀上，該結(jié)構(gòu)材料又由周期分布的胞元（圖1（c））構(gòu)成。基于均勻化方法，對該周期胞元進(jìn)行分析，即可評價整個微結(jié)構(gòu)材料的宏觀性能。

圖1 周期微結(jié)構(gòu)材料均勻化過程Fig.1 Homogenization Process of the Periodic Microstructural Materials

在實(shí)際應(yīng)用中，均勻化方法常采用有限元技術(shù)來求解。在有限元框架下，周期胞元的宏觀等效彈性矩陣可計算為：

式中：Y—周期胞元的體積；M—區(qū)域Y離散的單元數(shù)；Ve—第e個單元的體積，Ke—第e個單元的體積剛度矩陣；χ0(ij)e—初始單位應(yīng)變場ε0(ij)對應(yīng)的單元節(jié)點(diǎn)位移向量；χ(ij)e—初始單位應(yīng)變場ε0(ij)引起的特征位移向量。對微結(jié)構(gòu)胞元施加周期邊界約束［17］，χ(ij)e可以通過求解均勻化平衡方程來確定。

式中：K—整體剛度矩陣，由單元剛度矩陣Ke集成得到；M—微結(jié)構(gòu)胞元內(nèi)劃分的單元數(shù)；Be—第e個單元的應(yīng)變位移矩陣；De—第e個單元的彈性矩陣；ce— 單元到整體的轉(zhuǎn)換矩陣；ε0—初始單位應(yīng)變場，對于二維情況ε0(11)=[ 1 0 0]T、ε0(22)=

3 BESO算法及改進(jìn)

3.1 優(yōu)化模型

以材料的體積模量K或剪切模量G最大化為目標(biāo)，以基體材料的體積為約束，多孔微結(jié)構(gòu)材料的拓?fù)鋬?yōu)化設(shè)計模型可描述為

式中：V*—基體材料的體積約束值；xe—設(shè)計變量代表第e個單元的密度；當(dāng)xe= 1—該單元為實(shí)單元；當(dāng)xe=xmin—該單元為空單元（為避免求解困難，xmin通常取接近于零的小值，這里xmin取10-3）；對于各向異性材料，這里定義體積模量K和剪切模量G分別為：

3.2 靈敏度分析及單元濾波

在BESO方法中，仍然采用SIMP材料插值模型：

式中：E0—基體材料的彈性模量；p—懲罰因子，這里取p= 3。

從優(yōu)化目標(biāo)函數(shù)式（3）可以看出，目標(biāo)函數(shù)K或G對單元設(shè)計變量xe的偏微分，實(shí)際上是宏觀等效彈性參數(shù)DH ijkl對設(shè)計變量xe的偏微分。利用均勻化平衡方程式（2）及材料插值模型式（6），對方程式（1）求關(guān)于設(shè)計變量xe的偏微分［18］，可得：

式中：K0—單元e填滿基體材料時的剛度矩陣，K0=∫VeBTe D0BedVe（D0—基體材料的彈性矩陣）。

確定了單元靈敏度之后，采用權(quán)重方法［19］，將單元靈敏度轉(zhuǎn)化為節(jié)點(diǎn)靈敏度，即：

式中：—濾波后的第i個節(jié)點(diǎn)靈敏度；αj—以節(jié)點(diǎn)i為圓點(diǎn)、以rmin為半徑的圓內(nèi)的第j個單元對應(yīng)的目標(biāo)函數(shù)靈敏度，αj=或—該濾波圓內(nèi)的單元總數(shù)；w(rij)—權(quán)重因子，定義為：

式中：rmin—過濾半徑；rij—節(jié)點(diǎn)i與單元j之間的中心距離。

3.3 邊界單元細(xì)化

BESO方法是一種離散的拓?fù)鋬?yōu)化方法，得到的結(jié)構(gòu)存在明顯的鋸齒狀邊界。為了得到光滑的邊界，通常的做法是加密有限元網(wǎng)格，但網(wǎng)格加密勢必會造成計算量的急劇增加。這里采用形函數(shù)插值方法，對結(jié)構(gòu)拓?fù)溥吔缟系膯卧M(jìn)行虛擬細(xì)化，從而在不降低計算效率的情況下，達(dá)到光滑結(jié)構(gòu)拓?fù)溥吔绲哪康摹?/p>

單元具體細(xì)化過程為：首先將結(jié)構(gòu)邊界上的某個單元劃分為Ne個小單元（該小單元并不參與有限元計算），則第k個小單元的靈敏度可由該邊界單元的節(jié)點(diǎn)靈敏度計算得到：

式中：Nj(xk，yk)—邊界單元在小單元k中心坐標(biāo)(xk，yk)處的形函數(shù)—邊界單元第j個節(jié)點(diǎn)的靈敏度。

在優(yōu)化過程中，并不是所有單元都進(jìn)行上述細(xì)化處理，而僅細(xì)化結(jié)構(gòu)拓?fù)溥吔缟系膯卧Ｅ袛鄦卧欠袷沁吔鐔卧鶕?jù)節(jié)點(diǎn)靈敏度的大小排序確定，即邊界單元的節(jié)點(diǎn)靈敏度滿足

式中：—第e個單元的節(jié)點(diǎn)靈敏度向量；αth—節(jié)點(diǎn)靈敏度閥值，在滿足當(dāng)前迭代步的材料體積約束條件下，可以通過二分法求得［20］。

另外，在優(yōu)化過程中，邊界單元的體積為：

式中：N′—邊界單元e中小單元靈敏度αˉk高于閥值αth的單元總數(shù)。

3.4 改進(jìn)BESO算法步驟

根據(jù)前面分析，改進(jìn)BESO算法的主要計算步驟可概括為：

（1）設(shè)計域網(wǎng)格離散，初始化優(yōu)化參數(shù)：收斂精度ε、濾波半徑rmin、目標(biāo)體積約束V*、進(jìn)化率ER、單元細(xì)化個數(shù)Ne等；

（2）均勻化有限元分析，根據(jù)式（1）計算等效彈性模量，根據(jù)式（7）分析目標(biāo)函數(shù)的單元靈敏度；

（3）根據(jù)式（8）計算節(jié)點(diǎn)靈敏度；

（4）確定當(dāng)前迭代的目標(biāo)體積分?jǐn)?shù)：

（5）節(jié)點(diǎn)靈敏度排序，根據(jù)當(dāng)前材料體積約束，采用二分法確定靈敏度閥值αth，細(xì)化邊界單元；

（6）單元（包括小單元）刪、添，更新結(jié)構(gòu)；

（7）返回（2），直到滿足體積約束和目標(biāo)函數(shù)收斂精度。

4 周期微結(jié)構(gòu)材料設(shè)計

4.1 二維情況

對于長為100mm、寬為100mm的正方形設(shè)計域，將其離散為（50×50）個4節(jié)點(diǎn)矩形平面應(yīng)力單元。 其中，基體材料的楊氏模量E=1，泊松比m=0.3；同時以（10×10）網(wǎng)格對邊界單元進(jìn)行細(xì)化；其它優(yōu)化參數(shù)為ER=2%、濾波半徑rmin= 3、體積分?jǐn)?shù)約束為總設(shè)計域的30%，即=30%。

以最大化材料的等效剪切模量為優(yōu)化目標(biāo)，得到的材料微結(jié)構(gòu)胞元拓?fù)浣Y(jié)構(gòu)及其等效彈性矩陣，如圖2所示。從圖中可以看出，為了使剪切模量最大化，微結(jié)構(gòu)胞元中的基體材料幾乎嚴(yán)格按照45?角分布。目標(biāo)函數(shù)迭代收斂曲線及材料微結(jié)構(gòu)拓?fù)涞淖兓^程，如圖3所示。從圖中可以看出，目標(biāo)剪切模量是隨著材料的減少而逐漸遞減的。

圖2 具有最大剪切模量的材料微結(jié)構(gòu)胞元拓?fù)浼捌涞刃椥跃仃嘍H（Vf*=30%）Fig.2 Topology of the Material Microstructural Unit Cell with Maximum Shear Modulus and its Effective Elasticity Matrix DH（Vf*=30%）

圖3 目標(biāo)剪切模量收斂曲線及材料微結(jié)構(gòu)拓?fù)渥兓^程Fig.3 Evolution History of the Objective Shear Modulus and the Corresponding Material Microstructural Topologies

為了進(jìn)一步說明改進(jìn)BESO方法的有效性，這里在優(yōu)化參數(shù)相同的情況下，將改進(jìn)BESO方法與傳統(tǒng)BESO方法和SIMP方法進(jìn)行優(yōu)化對比，對比結(jié)果，如表1所示。從表中結(jié)果可以看出，與SIMP 方法相比，改進(jìn)BESO 方法和傳統(tǒng)BESO 方法一樣，都可以得到清晰的拓?fù)溥吔纾珎鹘y(tǒng)BESO結(jié)果中存在明顯的鋸齒狀邊界，而前者結(jié)果中邊界非常光滑；此外，從目標(biāo)函數(shù)上看，改進(jìn)BESO的結(jié)果也是最好的，明顯優(yōu)于SIMP的結(jié)果。由此證明了改進(jìn)BESO方法的有效性。

表1 不同優(yōu)化方法的二維微結(jié)構(gòu)材料優(yōu)化結(jié)果對比Tab.1 Comparison of Optimized Results for the Two-Dimensional Microstructure Material with Different Topology Optimization Methods

4.2 三維情況

對于三維微結(jié)構(gòu)材料，設(shè)定其周期胞元的設(shè)計域?yàn)檫呴L30mm的正方體。該設(shè)計域被離散為（30×30×30）個8節(jié)點(diǎn)正方體單元；邊界單元細(xì)化成（10×10×10）個網(wǎng)格的小單元；基體材料參數(shù)及其他優(yōu)化參數(shù)同二維情況。

基于改進(jìn)BESO方法，以微結(jié)構(gòu)材料等效體積模量最大化為設(shè)計目標(biāo)，表2給出了不同材料體積約束下的三維優(yōu)化結(jié)果。從中可以看出，材料微結(jié)構(gòu)拓?fù)渚哂星逦⒐饣倪吔纾欣谔岣咄負(fù)浣Y(jié)構(gòu)的CAD模型重構(gòu)精度；此外還可以看到，最優(yōu)微結(jié)構(gòu)拓?fù)涑省翱諝ぁ睜睿瑸榱俗畲蠡挚贵w積壓縮，基體材料總體上趨向軸向分布。

表2 基于改進(jìn)BESO方法的最大化三維微結(jié)構(gòu)材料體積模量優(yōu)化結(jié)果Tab.2 Optimized Results of the 3D Microstructural Material for Maximizing the Bulk Modulus Based in the Improved BESO Method

以微結(jié)構(gòu)材料等效剪切模量最大化為設(shè)計目標(biāo)的材料微結(jié)構(gòu)優(yōu)化結(jié)果，如表3所示。從中可以看出，材料微結(jié)構(gòu)拓?fù)錁?gòu)型依然呈“空殼”狀，且具有清晰、光滑的邊界，但基體材料分布情況明顯與最大化體積模量的優(yōu)化結(jié)果不同，此時，基體材料整體沿空間45?平面方向分布，以抵抗最大的剪切力。

表3 基于改進(jìn)BESO方法的最大化三維微結(jié)構(gòu)材料剪切模量優(yōu)化結(jié)果Tab.3 Optimized Results of the 3D Microstructural Material for Maximizing the Shear Modulus Based on the Improved BESO Method

5 結(jié)論

針對BESO方法邊界不光滑的問題，通過采用形函數(shù)插值細(xì)化邊界單元的方法，對BESO方法進(jìn)行了改進(jìn)；在此基礎(chǔ)上，運(yùn)用改進(jìn)的BESO方法對周期結(jié)構(gòu)材料進(jìn)行了優(yōu)化設(shè)計研究。研究結(jié)果表明，在不降低計算效率的情況下，改進(jìn)的BESO方法可以獲得更優(yōu)的優(yōu)化目標(biāo)和邊界更光滑的結(jié)構(gòu)拓?fù)洌换诟倪M(jìn)BESO拓?fù)鋬?yōu)化方法可以實(shí)現(xiàn)周期微結(jié)構(gòu)材料最大體積模量和最大剪切模量設(shè)計，并獲得相應(yīng)的微結(jié)構(gòu)拓?fù)錁?gòu)型。這里研究為實(shí)現(xiàn)人工微結(jié)構(gòu)材料功能化設(shè)計提供了一種有效途徑。