閆曉磊，陳佳文，張樹(shù)忠，張 勇，黃曉東

(1. 福建省汽車電子與電驅(qū)動(dòng)重點(diǎn)實(shí)驗(yàn)室(福建工程學(xué)院)，福建，福州 350118；2. 華僑大學(xué)機(jī)電及自動(dòng)化學(xué)院，福建，廈門 361021；3. 澳大利亞斯威本科技大學(xué)科學(xué)、工程與技術(shù)學(xué)院，維多利亞，Hawthorn 3122)

經(jīng)過(guò)近幾十年的發(fā)展，拓?fù)鋬?yōu)化已經(jīng)廣泛應(yīng)用于航天航空、機(jī)械工程、海洋工程等領(lǐng)域[1-3]。但拓?fù)鋬?yōu)化得到的結(jié)構(gòu)幾何形狀往往比較復(fù)雜，常用的制造工藝難以加工或成本太高[4-5]，因此，研究協(xié)同考量制造工藝約束的結(jié)構(gòu)設(shè)計(jì)方法具有明顯的工程實(shí)際意義。

在結(jié)構(gòu)設(shè)計(jì)中，極限尺寸是一類重要的制造約束。比如構(gòu)件的最小、最大尺寸以及孔洞的最小尺寸控制等。考慮尺寸約束是結(jié)構(gòu)拓?fù)鋬?yōu)化領(lǐng)域的一個(gè)難點(diǎn)，近年來(lái)，國(guó)內(nèi)外學(xué)者在這方面做了一些工作。

Poulsen[6]引入非負(fù)函數(shù)MOLE，對(duì)特定方向上的密度梯度施加約束，實(shí)現(xiàn)了最小尺寸控制。Guest 等[7]以節(jié)點(diǎn)密度為設(shè)計(jì)變量，通過(guò)調(diào)節(jié)節(jié)點(diǎn)到單元的映射半徑，實(shí)現(xiàn)拓?fù)浣Y(jié)構(gòu)的最小尺寸控制。類似地，白偉等[8]通過(guò)懲罰映射區(qū)域的單元平均密度，達(dá)到了控制結(jié)構(gòu)最大尺寸的目的。Zhou 等[9]在三重密度場(chǎng)(設(shè)計(jì)密度、濾波密度、映射密度)框架下，通過(guò)聯(lián)合調(diào)節(jié)濾波半徑、映射閾值及其范圍，提出了一種可調(diào)控最小尺寸的拓?fù)鋬?yōu)化方法。采用類似的處理方法，Dunning[10]將其引入基于參數(shù)化隱式函數(shù)的拓?fù)鋬?yōu)化方法，實(shí)現(xiàn)了最小尺寸控制。以上研究通過(guò)對(duì)設(shè)計(jì)變量施加約束，間接地隱式地實(shí)現(xiàn)了最小尺寸控制，但往往無(wú)法實(shí)現(xiàn)特征尺寸的精確控制[11]。

Zhang 等[12]基于SIMP(Solid Isotropic Material with Penalization)拓?fù)鋬?yōu)化方法，借助于拓?fù)浣Y(jié)構(gòu)骨架實(shí)現(xiàn)了最大、最小尺寸的精確控制；Guo 等[13]基于水平集方法，利用符號(hào)距離水平集函數(shù)提取了拓?fù)浣Y(jié)構(gòu)的骨架，并通過(guò)約束其極值，實(shí)現(xiàn)了拓?fù)浣Y(jié)構(gòu)的最大、最小尺寸控制；Niu 等[14]基于MMC(Moving Morphable Components)方法，通過(guò)在組件的寬度上施加相同的上限和下限，實(shí)現(xiàn)了結(jié)構(gòu)的等尺寸控制；Liu[15]基于水平集方法，通過(guò)對(duì)拓?fù)浣Y(jié)構(gòu)骨架的識(shí)別與分割，將拓?fù)浣Y(jié)構(gòu)分解為多個(gè)條狀組件，實(shí)現(xiàn)了結(jié)構(gòu)的動(dòng)態(tài)極限尺寸控制。以上研究，借助拓?fù)浣Y(jié)構(gòu)的骨架，顯式地實(shí)現(xiàn)了結(jié)構(gòu)特征尺寸的精確控制，但這些方法難以考慮骨架變化對(duì)設(shè)計(jì)變量靈敏度的影響[16]。

漸進(jìn)結(jié)構(gòu)優(yōu)化方法(Bi-directional Evolutionary Structural Optimization, BESO)是一種啟發(fā)式的拓?fù)鋬?yōu)化方法，具有原理簡(jiǎn)單、效率高、易實(shí)現(xiàn)的特點(diǎn)[17-18]。因此，本文嘗試以BESO 方法為基礎(chǔ)，通過(guò)細(xì)化邊界單元及提取拓?fù)浣Y(jié)構(gòu)骨架，提出一種高效的拓?fù)鋬?yōu)化方法，避免繁瑣的公式推導(dǎo)的同時(shí)，實(shí)現(xiàn)最小尺寸的精確控制。

1 BESO 算法改進(jìn)

BESO 拓?fù)鋬?yōu)化方法，雖然具有效率高、原理簡(jiǎn)單的優(yōu)點(diǎn)，但優(yōu)化結(jié)果存在鋸齒狀邊界問(wèn)題[19]，無(wú)法直接用于CAD 幾何建模，也不利于特征尺寸的控制。對(duì)此，本文首先對(duì)傳統(tǒng)BESO 方法進(jìn)行改進(jìn)。

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

以最小化結(jié)構(gòu)平均柔度為目標(biāo)為例，拓?fù)鋬?yōu)化模型可定義為：

式中：xi為第i個(gè)單元的拓?fù)湓O(shè)計(jì)變量，當(dāng)xi=1 表示該單元為實(shí)單元，當(dāng)xi=xmin表示該單元為空單元(為避免求解困難，xmin通常取接近于零的小值，本文xmin取10-3)；C為結(jié)構(gòu)的平均柔度；U為節(jié)點(diǎn)位移向量；K為整體剛度矩陣；Vf*為目標(biāo)體積分?jǐn)?shù)；Vi為單元i的面積或體積；M為設(shè)計(jì)域單元總數(shù)。

1.2 靈敏度分析及單元濾波

在BESO 算法中，仍然采用SIMP 插值模型，即第i個(gè)單元的彈性模量可表示為：

式中：E0為基體材料的彈性模量；p為懲罰因子(本文取3)。

根據(jù)目標(biāo)函數(shù)式(1)，再結(jié)合插值模型式(2)，可求得目標(biāo)平均柔度C對(duì)拓?fù)湓O(shè)計(jì)變量xi的偏微分為：

有了靈敏度分析式(3)，單元敏度數(shù)可以定義為：

為了抑制棋盤格模式和網(wǎng)格依賴性，同時(shí)將單元敏度數(shù)轉(zhuǎn)化為節(jié)點(diǎn)敏度數(shù)，這里采用權(quán)重濾波方法[20]：

式中：rmin為過(guò)濾半徑；rji為節(jié)點(diǎn)j和單元i中心的距離。

1.3 邊界單元細(xì)化

為了得到光滑的邊界，通常的做法是加密有限元網(wǎng)格，但網(wǎng)格增加勢(shì)必會(huì)造成計(jì)算量的急劇增加。本文采用形函數(shù)插值方法，對(duì)結(jié)構(gòu)拓?fù)溥吔缟系膯卧M(jìn)行虛擬細(xì)化，在不增加有限元網(wǎng)格的情況下，達(dá)到光滑結(jié)構(gòu)拓?fù)溥吔绲哪康?。下面給出具體細(xì)化過(guò)程。

首先，將結(jié)構(gòu)邊界上的某個(gè)單元?jiǎng)澐譃镹e個(gè)小單元(該小單元并不參與有限元計(jì)算)，則第k個(gè)小單元的靈敏度可由該邊界單元的節(jié)點(diǎn)靈敏度計(jì)算得到：

在優(yōu)化過(guò)程中，并不是所有單元都進(jìn)行上述細(xì)化處理，而僅細(xì)化結(jié)構(gòu)拓?fù)溥吔缟系膯卧?。判斷單元是否是邊界單元，要根?jù)節(jié)點(diǎn)靈敏度的大小排序確定，即邊界單元的節(jié)點(diǎn)靈敏度滿足：

2 拓?fù)浣Y(jié)構(gòu)骨架提取及最小尺寸控制

2.1 骨架提取算法

圖形骨架具有中心性、同倫、單位性且不會(huì)改變?cè)瓐D形的拓?fù)涮匦?。因此，采用骨架?duì)尺寸精確定義以及控制具有明顯優(yōu)勢(shì)。本文所采用拓?fù)浼?xì)化算法[22]來(lái)提取拓?fù)浣Y(jié)構(gòu)骨架。

拓?fù)浼?xì)化算法的核心就是逐層刪去結(jié)構(gòu)邊界上的簡(jiǎn)單點(diǎn)(Simple point)，所謂簡(jiǎn)單點(diǎn)就是刪除之后不改變圖形拓?fù)涮卣鞯狞c(diǎn)。在有限元框架下，若某邊界實(shí)單元v滿足下面的公式，則視該單元為簡(jiǎn)單點(diǎn)[23-24]：

式中：N3(v)為單元v的26 個(gè)鄰域單元，即3×3×3基體，如圖1 所示；δ 為變量的改變；G為該單元基體N3(v)的歐拉示性數(shù)；O為基體N3(v)的連接體數(shù)；S為實(shí)體單元。若第一式成立表示該單元滿足歐拉特性；若第二式成立表示刪除該邊界單元后連通性不改變。

圖 1 N2(v)和N3(v)基體Fig. 1 Frameworks for N2(v) and N3(v)

在實(shí)際操作中，第一步是選出設(shè)計(jì)域的邊界，若N3(v)的26 鄰域單元中只有一個(gè)實(shí)單元，則單元v為邊界點(diǎn)中的端點(diǎn)，直接保留，繼續(xù)第二步。

第二步驗(yàn)證歐拉特性。將3×3×3 基體N3(v)劃分為8 個(gè)2×2×2 基體N2(v)，則單元v的歐拉特性可計(jì)算為[25]：

式中，G6為單元v的6 鄰域歐拉特性，在第i個(gè)N2(v)基體中，單元v的6 鄰域歐拉示性數(shù)的變化量可計(jì)算為：

式中，ei、fi、octi分別為單元i所在N2(v)基體中實(shí)單元所形成的邊、面和基體的數(shù)目。

第三步驗(yàn)證連通性。可采用八叉樹(shù)算法[26]，計(jì)算邊界實(shí)單元v的連接體數(shù)，若計(jì)算得到的連接體數(shù)為1，同時(shí)邊界單元v刪除時(shí)式(11)為0，則該單元為簡(jiǎn)單點(diǎn)，直接刪除，否則，保留。

2.2 骨架毛刺去除

利用第2.1 節(jié)過(guò)程，可提取到每一步優(yōu)化后拓?fù)浣Y(jié)構(gòu)的骨架，但在拓?fù)浣Y(jié)構(gòu)邊緣處的骨架可能會(huì)生成毛刺[27]，使骨架形狀失真。因此，在得到拓?fù)浣Y(jié)構(gòu)骨架后，還需要進(jìn)行去毛刺處理。去除毛刺的原理比較簡(jiǎn)單，若骨架端點(diǎn)到連結(jié)點(diǎn)的長(zhǎng)度小于閾值Lth，即認(rèn)為該骨架分支為毛刺，進(jìn)行刪除。需要注意的是，在去除毛刺時(shí)，可能會(huì)將原連結(jié)點(diǎn)轉(zhuǎn)變?yōu)樾碌亩它c(diǎn)，從而產(chǎn)生新的毛刺，因此，毛刺的去除須分級(jí)去除，直到無(wú)毛刺可去除。本文毛刺閾值Lth取10 個(gè)單元的長(zhǎng)度。

2.3 最小尺寸控制

在迭代過(guò)程中，得到拓?fù)浣Y(jié)構(gòu)的骨架單元后，以最小控制尺寸dmin為直徑，以骨架單元的形心為球心，掃掠整個(gè)骨架，即可確定最小尺寸對(duì)應(yīng)的控制域B(dmin)。在優(yōu)化過(guò)程中，始終保證形心在控制域B(dmin)內(nèi)的單元全部為實(shí)單元，從而實(shí)現(xiàn)拓?fù)浣Y(jié)構(gòu)最小尺寸的精確控制。

此外，需要說(shuō)明的是，拓?fù)鋬?yōu)化過(guò)程中，拓?fù)錁?gòu)型的演化是通過(guò)空洞的生成與合并形成的，如果直接施加最小尺寸約束，可能會(huì)造成優(yōu)化結(jié)果早熟[12]，得不到最優(yōu)解。為避免這一現(xiàn)象，本文采用逐漸強(qiáng)化約束的方法，即構(gòu)造最小尺寸約束為：

3 拓?fù)鋬?yōu)化最小尺寸控制算法步驟

根據(jù)上面分析，基于BESO 的最小尺寸控制拓?fù)鋬?yōu)化算法可概括為以下步驟：

步驟1：設(shè)計(jì)域網(wǎng)格離散，初始化設(shè)計(jì)參數(shù)：收斂精度ε、目標(biāo)體積分?jǐn)?shù)Vf*、進(jìn)化率ER、過(guò)濾半徑rmin及最小尺寸dmin等；

步驟2：對(duì)結(jié)構(gòu)進(jìn)行有限元方法分析，根據(jù)式(3)分析目標(biāo)函數(shù)的單元靈敏度；

步驟3：根據(jù)式(5)計(jì)算節(jié)點(diǎn)靈敏度；

步驟4：確定當(dāng)前迭代步材料目標(biāo)體積分?jǐn)?shù)：

步驟5：節(jié)點(diǎn)靈敏度排序，細(xì)化邊界單元；

步驟6：?jiǎn)卧?包括小單元)刪、添；

步驟7：提取拓?fù)浣Y(jié)構(gòu)的骨架；

步驟8：根據(jù)式(13)，確定控制域B(dmin)，施加最小尺寸約束；

步驟9：更新結(jié)構(gòu)；

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

4 實(shí)例分析

下面通過(guò)實(shí)例來(lái)驗(yàn)證所提出的最小尺寸控制拓?fù)鋬?yōu)化算法。為簡(jiǎn)化起見(jiàn)，以下所有實(shí)例中，材料的彈性模量和泊松比都分別取為E=1 GPa、ν=0.3；進(jìn)化率ER=0.01；另外，二維單元細(xì)化規(guī)模為10×10 個(gè)小單元、三維單元細(xì)化規(guī)模為10×10×10 個(gè)小單元。

4.1 算例1. L 型梁結(jié)構(gòu)

首先，為驗(yàn)證算法的有效性，對(duì)于如圖2 所示L 型梁結(jié)構(gòu)，右端中點(diǎn)受大小為1 kN 的集中載荷作用，上邊界固支；初始設(shè)計(jì)域被離散為5120 個(gè)矩形有限元網(wǎng)格；過(guò)濾半徑rmin=2 mm，最小控制尺寸為dmin=2 mm。

圖 2 L 型梁結(jié)構(gòu)設(shè)計(jì)域Fig. 2 Design domain of a L-shape beam

表 1 L 型梁不同體積分?jǐn)?shù)下的最小尺寸控制優(yōu)化結(jié)果Table 1 The optimized results with minimum length scale control under different volume fractions for the L-shape beam

圖3 給出了表1 對(duì)應(yīng)的最小尺寸約束下不同目標(biāo)體積分?jǐn)?shù)約束時(shí)的目標(biāo)函數(shù)收斂曲線。優(yōu)化過(guò)程體積分?jǐn)?shù)逐漸減少到目標(biāo)體積分?jǐn)?shù)，目標(biāo)柔度逐漸增大直至收斂，收斂過(guò)程十分平穩(wěn)。

4.2 算例2. 懸臂梁結(jié)構(gòu)

先來(lái)研究最小尺寸約束施加過(guò)程，即參數(shù)β 的影響，此時(shí)固定最小尺寸約束為dmin=2 mm。圖5 對(duì)比了β=5、β=10、β=20 情況及無(wú)尺寸約束的優(yōu)化結(jié)果。從圖5 中可以看出，最小尺寸約束的施加對(duì)優(yōu)化結(jié)果具有顯著的影響。β 較小時(shí)，由于整個(gè)優(yōu)化過(guò)程最小尺寸約束都在起作用，此時(shí)優(yōu)化得到的拓?fù)浣Y(jié)構(gòu)中存在很多桿件，并且目標(biāo)柔度變大(相對(duì)于無(wú)尺寸約束情況)；隨著β 的增加，最小尺寸約束相當(dāng)于施加在優(yōu)化過(guò)程的后期，此時(shí)拓?fù)浣Y(jié)構(gòu)中桿件逐漸變少，結(jié)構(gòu)趨于簡(jiǎn)化，目標(biāo)柔度也逐漸接近無(wú)尺寸約束情況。事實(shí)上，在拓?fù)錁?gòu)型演化中，空洞的生成與合并過(guò)程也就是最小尺寸“桿件”的產(chǎn)生與斷裂過(guò)程，最小尺寸約束的施加勢(shì)必會(huì)影響這一過(guò)程，造成優(yōu)化過(guò)早收斂[10,12]，而通過(guò)松弛施加最小尺寸約束可以有效避免這一問(wèn)題。

圖 3 最小尺寸控制下的目標(biāo)柔度及體積分?jǐn)?shù)迭代曲線 (dmin=2 mm)Fig. 3 Iterative history curves of the objective compliance and the volume fraction with the minimum length scale control (dmin=2 mm)

圖 4 懸臂梁設(shè)計(jì)域Fig. 4 Design domain of the cantilever beam

圖 5 不同β 值的拓?fù)鋬?yōu)化(dmin=2 mm)Fig. 5 Topology optimization results with different index β (dmin=2 mm)

下面研究最小控制尺寸dmin對(duì)結(jié)果的影響。圖6 對(duì)比了最小尺寸分別為dmin=2 mm、dmin=2.5 mm、dmin=3 mm 及無(wú)尺寸約束情況下的拓?fù)鋬?yōu)化結(jié)果，此時(shí)，β 設(shè)定為20。從圖6 中可以看出，隨著最小控制尺寸dmin的增大，結(jié)構(gòu)中的最細(xì)桿件逐漸變粗，并且都滿足了最小尺寸控制要求；但隨dmin的增大，結(jié)構(gòu)拓?fù)渥兊脧?fù)雜，目標(biāo)柔度也逐漸增大；可見(jiàn)，最小控制尺寸dmin的大小對(duì)優(yōu)化結(jié)果具有顯著的影響，設(shè)置最小控制尺寸可以避免細(xì)小支桿的產(chǎn)生，提高結(jié)構(gòu)的穩(wěn)定性，但過(guò)大的最小控制尺寸可能使設(shè)計(jì)結(jié)果性能變差且結(jié)構(gòu)復(fù)雜。

4.3 算例3. 三維支座

圖 6 不同最小尺寸約束的拓?fù)鋬?yōu)化結(jié)果Fig. 6 Topology optimization results under different minimum length scale controls

圖 7 三維支座設(shè)計(jì)域Fig. 7 Design domain of a 3D support

圖8 對(duì)比了最終的拓?fù)鋬?yōu)化結(jié)果，其中圖8(a)為無(wú)尺寸控制的拓?fù)鋬?yōu)化結(jié)果；圖8(b)為施加了最小尺寸控制的拓?fù)鋬?yōu)化結(jié)果；圖8(c)為對(duì)應(yīng)的拓?fù)涔羌芙Y(jié)構(gòu)；對(duì)比圖8(a)與圖8(b)，可以看到，支座中最細(xì)桿的尺寸可以被精確控制，但施加最小尺寸控制后，結(jié)構(gòu)的整體剛度略有下降(目標(biāo)柔度略有上升)；圖8(d)為體積分?jǐn)?shù)及目標(biāo)函數(shù)收斂曲線，從圖中可以看出，施加最小尺寸控制后，迭代收斂過(guò)程依然很平穩(wěn)。

圖 8 三維支座施加最小尺寸控制的拓?fù)鋬?yōu)化結(jié)果Fig. 8 Optimized results with minimum length scale control for the 3D support

5 結(jié)論

本文在改進(jìn)BESO 算法的基礎(chǔ)上，通過(guò)提取拓?fù)浣Y(jié)構(gòu)的骨架，提出了一種可控制最小特征尺寸的拓?fù)鋬?yōu)化方法。二維和三維算例表明：

(1) 該方法可以顯式地實(shí)現(xiàn)結(jié)構(gòu)最小尺寸精確控制；

(2) 合理地設(shè)置最小控制尺寸可以避免結(jié)構(gòu)產(chǎn)生細(xì)小支桿，提高結(jié)構(gòu)的穩(wěn)定性，但過(guò)大的最小控制尺寸會(huì)使結(jié)構(gòu)變得復(fù)雜且性能變差；

(3) 最小尺寸約束的施加對(duì)拓?fù)錁?gòu)型演化過(guò)程影響顯著，本文通過(guò)最小尺寸約束松弛施加的方法，有效避免了優(yōu)化過(guò)程早熟的問(wèn)題。

在后面工作中，作者將進(jìn)一步對(duì)拓?fù)浣Y(jié)構(gòu)孔洞尺寸及其他制造工藝約束的協(xié)同控制進(jìn)行研究。