周健松，李海艷

（廣東工業(yè)大學(xué)機(jī)電工程學(xué)院，廣州 510006）

0 引言

任何具有一定使用功能的構(gòu)件都是由滿足要求的材料所制造的，在設(shè)計(jì)階段，就需要對(duì)該構(gòu)件在可能的外力作用下的內(nèi)部狀態(tài)進(jìn)行分析[1]。拓?fù)鋬?yōu)化作為結(jié)構(gòu)優(yōu)化的一種數(shù)學(xué)設(shè)計(jì)工具，是在給定的荷載條件約束和性能指標(biāo)下，對(duì)材料在給定區(qū)域內(nèi)的分布進(jìn)行優(yōu)化的一種方法。近年來，它在航空航天、機(jī)械、建筑、醫(yī)療和其他領(lǐng)域發(fā)展迅速，并取得了一定的研究成果[2-8]。目前，關(guān)于有限元的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化被廣泛應(yīng)用于各種行業(yè)軟件中，有限元的優(yōu)化問題，往往存在棋盤格和網(wǎng)格依賴性以及模糊區(qū)域的非0-1解的情況[9]，增大了實(shí)際結(jié)構(gòu)設(shè)計(jì)的難度。

為了防止棋盤格圖案和中間設(shè)計(jì)變量，國(guó)內(nèi)外有許多研究者提出各種可行方案。最早被提出的是用具有啟發(fā)式的平滑卷積濾波方法來解決去除棋盤格問題[10-11]，該過濾一開始只是應(yīng)用有限元網(wǎng)格目標(biāo)函數(shù)的靈敏度上，Sigmund 在此基礎(chǔ)上提出密度過濾方法，并將這兩種理論總結(jié)在MATLAB 上只有99 行和88 行代碼的拓?fù)鋬?yōu)化應(yīng)用。但這種簡(jiǎn)單的卷積過濾會(huì)產(chǎn)生大量灰色模糊區(qū)域，因此很多學(xué)者嘗試在基礎(chǔ)上采用一些后處理算法抑制中間值的設(shè)計(jì)變量，如Guest[12]采用的Heaviside 階躍函數(shù)和Sigmund利用圖像處理技術(shù)[13]，盡可能讓中間值的設(shè)計(jì)變量趨向0或1。為此需要不斷調(diào)整額外函數(shù)的參數(shù)確保收斂，增加迭代次數(shù)，這樣可能會(huì)增加拓?fù)鋬?yōu)化求解時(shí)間。

從過濾器的角度分析，拓?fù)鋬?yōu)化希望能在棋盤格區(qū)域內(nèi)進(jìn)行平滑處理，而在特征邊緣之間弱化平滑作用。本文基于圖像處理的方法，結(jié)合局部鄰域內(nèi)有限元之間的距離以及局部密度變化作為權(quán)重，將傳統(tǒng)的各向同性過濾轉(zhuǎn)為一種新的各向異性過濾密度方法。一般情況下，實(shí)體之間的過濾擴(kuò)散效果幾乎不變，但是特征邊界處的過濾效果會(huì)因?yàn)榭瞻讍卧c實(shí)體單元之間不連續(xù)，平滑擴(kuò)散作用大幅度降低，導(dǎo)致邊界單元只會(huì)傾向于{0，1}的設(shè)計(jì)變量解，有效減少灰色區(qū)域的可能性。本文介紹拓?fù)鋬?yōu)化的SIMP數(shù)學(xué)模型，以及基礎(chǔ)求解方法；提出各向異性的密度過濾方法比與經(jīng)典密度過濾方法進(jìn)行對(duì)比；通過常用的拓?fù)鋬?yōu)化案例進(jìn)行分析，驗(yàn)證該方法的可行性；并通過目標(biāo)柔度、灰度指標(biāo)以及求解時(shí)間展示其優(yōu)勢(shì)。

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

有限元拓?fù)鋬?yōu)化需要將材料域離散化成有限元網(wǎng)格，通過有限元分析尋找最優(yōu)材料分布。目前最常用的有限元網(wǎng)格拓?fù)鋬?yōu)化數(shù)學(xué)模型是基于SIMP（Solid Isotropic Microstructure with Penalization） 模型基礎(chǔ)上進(jìn)行的[14]，本文利用最小化柔度作為例進(jìn)行拓?fù)鋬?yōu)化。假設(shè)單元中設(shè)計(jì)變量是處于0～1之間，對(duì)應(yīng)于設(shè)計(jì)域內(nèi)材料分布為空隙或者實(shí)體的程度，材料的物理屬性與相對(duì)密度之間存在非線性的關(guān)系。利用懲罰系數(shù)p，使得0～1中間相對(duì)密度剛度降低，從而抑制設(shè)計(jì)變量趨向中間密度值。在一定的材料數(shù)量約束下，在設(shè)計(jì)域中尋求特定邊界約束和負(fù)載下，最小化柔度（最大剛度）的材料分配情況如式（1）所示。

式中：ρ為有限單元網(wǎng)格的相對(duì)密度；C（ρ）為目標(biāo)結(jié)構(gòu)柔度；U為整體變形位移；K為整體剛度矩陣；F為負(fù)載；V為單元體積；Vvol為體積約束值。由于整體剛度矩陣由各單元組裝而成，目標(biāo)函數(shù)也可以寫成：

式中：E為材料彈性變形的彈性模量；k0為單元?jiǎng)偠染仃?；ui為第i個(gè)單元變形位移；N為設(shè)計(jì)域離散劃分的單元數(shù)量。由于SIMP懲罰系數(shù)的作用，通常設(shè)置非常小的彈性模量Emin= 10-9E0，防止當(dāng)單元相對(duì)密度ρi=0，導(dǎo)致剛度矩陣發(fā)生奇異。一般情況下，懲罰系數(shù)p=3。

在單一的設(shè)計(jì)變量的情況下，會(huì)選擇優(yōu)化準(zhǔn)則算法作為求解算法[15]。

式中：λ為拉格朗日約束非負(fù)系數(shù)，控制結(jié)構(gòu)體積不大于預(yù)設(shè)值。

2 密度過濾

一般低階單元的拓?fù)鋬?yōu)化存在棋盤格問題，如圖1所示，該案例是經(jīng)典的MBB梁結(jié)構(gòu)，長(zhǎng)寬比3∶1，左側(cè)邊界對(duì)橫向約束，右下角垂直向上固定，受左上角垂直向下的單位力。若采用四節(jié)點(diǎn)單元，其有限元求解得到只是節(jié)點(diǎn)處位移變形，根據(jù)節(jié)點(diǎn)之間位移與剛度的關(guān)系，其連續(xù)性只適用于節(jié)點(diǎn)處，無法保證單元之間邊界具有連續(xù)性。

圖1 無過濾的MBB拓?fù)鋬?yōu)化結(jié)果

2.1 經(jīng)典各向同性過濾器

Sigmund 提出的過濾器的概念，在原有設(shè)計(jì)變量ρ的基礎(chǔ)上額外增加了實(shí)際密度變量ρ，并提出密度過濾方法[16]。利用圖像去噪的卷積算法，對(duì)其平滑化處理。

式中：Ω為單元的過濾范圍所有單元集合，{Ω:≤rmin}；H為關(guān)于各單元之間距離的函數(shù)。

與密度相關(guān)材料屬性也會(huì)隨之變化。

利用Sigmund提出88行密度過濾的拓?fù)鋬?yōu)化代碼[16]，得到過濾后的拓?fù)鋬?yōu)化結(jié)果，往往存在大量灰色區(qū)域（設(shè)計(jì)變量不為0 或1，可參考圖3（a））。在實(shí)際應(yīng)用中，這種拓?fù)鋬?yōu)化結(jié)果難以判斷其結(jié)構(gòu)邊界。

2.2 新各向異性過濾器

經(jīng)典卷積過濾算法的結(jié)果存在灰色邊界最大的原因在于，在邊界處空白單元（相對(duì)密度為0）和黑色單元（相對(duì)密度為1）具有相似的距離，兩者權(quán)重相當(dāng)，因此經(jīng)過密度過濾后產(chǎn)生大量灰色區(qū)域。若過濾器本身具有識(shí)別邊界特征能夠有效減少灰色區(qū)域的分布，雙邊濾波平滑器在去噪過程中，通過非線性組合附近的圖像值，根據(jù)相似度修正權(quán)重因子，從而保留邊緣[17]。本文結(jié)合圖像去噪的雙邊去噪技術(shù)，在將原有距離作為權(quán)重的基礎(chǔ)上，添加不同單元之間密度差作為權(quán)重，形成各向異性的過濾器。

式中：σr與過濾距離的參數(shù)相關(guān)，σr=rmin∕3；σm與單元相對(duì)密度值相關(guān)，一般可取1∕3。比較式（6）和（9），此時(shí)，權(quán)重包括兩個(gè)成分——空間和密度。第一個(gè)權(quán)值是關(guān)于中心單元[x，y]與局部區(qū)域單元[x+i，y+j]之間的幾何距離。由于σr與過濾距離的參數(shù)相關(guān)，可視作空間距離的一個(gè)標(biāo)準(zhǔn)，該部分與原本的卷積過濾算法區(qū)別不大，一般情況取值為3～8，也可根據(jù)實(shí)際制造加工的最小物理參數(shù)進(jìn)行適當(dāng)調(diào)整。因此，過濾半徑可作為結(jié)構(gòu)最小特征值，隨著網(wǎng)格數(shù)量的增減而變動(dòng)，解決拓?fù)鋬?yōu)化的網(wǎng)格依賴性問題。第二個(gè)權(quán)值是基于中心樣本[x，y]與[x+i，y+j]的設(shè)計(jì)變量差異，不同單元之間密度差并不固定，因此過濾區(qū)域內(nèi)所有單元權(quán)重并不相同，其目的是改善邊界處的過濾效果，使得最終結(jié)果更傾向黑白解。不難看出σm的取值會(huì)影響拓?fù)鋬?yōu)化灰度的情況，σm越小，灰色模糊區(qū)域就越少；相反，當(dāng)σm越大趨近于無窮時(shí)，則該方法的結(jié)果接近與傳統(tǒng)卷積過濾的結(jié)果。

以拓?fù)鋬?yōu)化結(jié)果邊界為例，假設(shè)邊界處單元設(shè)計(jì)變量密度值接近1，邊界附近代表實(shí)體的單元權(quán)重更大。反之，若邊界處的單元設(shè)計(jì)變量密度值接近0，對(duì)于孔隙的單元權(quán)重更大。其平滑擴(kuò)散因子會(huì)隨著迭代過程不斷改進(jìn)，最終在邊界方向上得到收斂，使得設(shè)計(jì)變量更容易趨向于0或1，有效避免灰色區(qū)域的邊界。而存在棋盤格問題區(qū)域是交錯(cuò)分布的黑白格，仍然能發(fā)揮過濾器的平滑作用，在迭代過程逐漸消除灰色區(qū)域，最終形成保留無棋盤格的黑白解結(jié)果。

3 案例分析

以上述經(jīng)典的MBB 梁作為案例，將其拓展測(cè)試該方法的可行性，采用不同數(shù)量的網(wǎng)格進(jìn)行對(duì)比分析，同時(shí)分析該方法網(wǎng)格依賴性問題。除了目標(biāo)柔度之外，Sigmund還提出一個(gè)測(cè)量拓?fù)鋬?yōu)化結(jié)果灰色區(qū)域情況的指標(biāo)Mnd[13]：

顯然，當(dāng)ρ越接近0.5時(shí)，Mnd越接近1。相反地，當(dāng)ρ越接近0 或者1 時(shí)，Mnd越接近0。因此該指標(biāo)能夠代表整體拓?fù)鋬?yōu)化結(jié)果的灰色區(qū)域情況。在下面案例比較中，會(huì)把目標(biāo)柔度、灰度指標(biāo)以及求解速度展示出來。

3.1 MBB梁以及網(wǎng)格依賴性

以經(jīng)典MBB 梁作為案例，按比例選擇以下4 種不同網(wǎng)格數(shù)量為例子，進(jìn)行拓?fù)鋬?yōu)化。

除了圖2（a），其他的測(cè)試過濾過濾半徑與網(wǎng)格數(shù)量成比例。這是因?yàn)閳D2（a）的網(wǎng)格數(shù)量太少，等比例選擇的過濾半徑太小會(huì)導(dǎo)致過濾效果不明顯（當(dāng)rmin取值太小，過濾器失去本身的價(jià)值），因此適當(dāng)增大過濾半徑，但整體拓?fù)鋬?yōu)化結(jié)構(gòu)仍然具有高度相似性。因此可知，該方法具有一定的網(wǎng)格獨(dú)立性，其拓?fù)鋬?yōu)化結(jié)構(gòu)不會(huì)隨有限元網(wǎng)格數(shù)量影響而發(fā)生明顯差異。不同網(wǎng)格數(shù)量MBB 優(yōu)化結(jié)果數(shù)據(jù)如表1 所示。由表可知，網(wǎng)格數(shù)量較小的情況下能夠更快速求解，這是因?yàn)檫^濾半徑和網(wǎng)格數(shù)量增大使得有限元矩陣行列數(shù)增加，在過濾器運(yùn)行時(shí)間也相應(yīng)增加。但是由于該方法得到的拓?fù)鋬?yōu)化解很少出現(xiàn)非0-1 解，因此作為緩沖的灰色網(wǎng)格區(qū)域較小，適當(dāng)增加有限元網(wǎng)格數(shù)量有利于得到視覺上較平滑的拓?fù)鋬?yōu)化結(jié)構(gòu)，同時(shí)防止過濾半徑較小導(dǎo)致棋盤格現(xiàn)象。

圖2 不同數(shù)量網(wǎng)格下MBB梁的拓?fù)鋬?yōu)化結(jié)果

表1 不同網(wǎng)格數(shù)量MBB優(yōu)化結(jié)果數(shù)據(jù)

3.2 懸臂梁

第二個(gè)案例選擇懸臂梁，假設(shè)懸臂梁長(zhǎng)寬比4∶1，左側(cè)邊界單元的全部固定約束，同時(shí)只有右下角有一個(gè)垂直向下的單位負(fù)載力。在本文中，選用200×50的網(wǎng)格數(shù)量進(jìn)行拓?fù)鋬?yōu)化。

將新過濾方法與最常見的密度過濾方法[16]和 帶Heaviside階躍函數(shù)過濾[12]的結(jié)果作比較，結(jié)果如圖3 所示。由圖可知，相比于其他兩種密度過濾方法，新的各向異性過濾方法明顯不存在灰色模糊區(qū)域，得到一個(gè)清晰邊界結(jié)果，其結(jié)構(gòu)更加細(xì)致、明確。不同過濾方法的MBB 結(jié)果數(shù)據(jù)如表2 所示。由表可知，在相同的情況下，新的各向異性過濾能在很短時(shí)間和迭代次數(shù)內(nèi)得到滿足收斂條件的結(jié)果，同時(shí)灰度測(cè)量值Mnd遠(yuǎn)小于其他兩個(gè)方法的結(jié)果，而目標(biāo)柔度與帶Heaviside 階躍函數(shù)密度過濾結(jié)果相差不大。與傳統(tǒng)密度過濾方法相比，只用15%的時(shí)間得到目標(biāo)柔度降低10%的拓?fù)鋬?yōu)化結(jié)果。

圖3 不同過濾器的懸臂梁的拓?fù)鋬?yōu)化結(jié)果

表2 不同過濾方法的MBB結(jié)果數(shù)據(jù)

4 結(jié)束語

本文以有限元拓?fù)鋬?yōu)化作為研究對(duì)象，考慮過濾器對(duì)優(yōu)化結(jié)果的影響，提出一個(gè)新的各向異性的過濾器。針對(duì)拓?fù)浣Y(jié)構(gòu)在有限元網(wǎng)格之間存在一定的不連續(xù)性，所提方法可修正過濾區(qū)域之間的權(quán)重。基于各向異性的特點(diǎn)，新方法能夠穩(wěn)定快速求解，使得邊界平滑過濾得到修正，有效解決拓?fù)鋬?yōu)化中低階網(wǎng)格的棋盤格和網(wǎng)格依賴性問題。同時(shí)，在拓?fù)鋬?yōu)化過程中收斂成一個(gè)穩(wěn)定的黑白解，減少模糊灰度區(qū)域。與傳統(tǒng)的密度過濾方法相比，本文所提方法求解時(shí)間、目標(biāo)柔度以及灰度指標(biāo)都得到大幅度的改善。與帶階躍函數(shù)的密度過濾器相比，在得到差距不大的目標(biāo)函數(shù)的拓?fù)浣Y(jié)果情況下，本文所提方法的求解效率明顯更高。

為了進(jìn)一步研究該方法對(duì)實(shí)際拓?fù)鋬?yōu)化的影響，接下來將建立三維空間的拓?fù)鋬?yōu)化模型，并且優(yōu)化過濾器權(quán)重計(jì)算方式，進(jìn)一步提高計(jì)算效率。