基于改進雙向漸進法的輕量化設計方法

王穎　張爭艷　雷焱陽　唐瑀　孫立新

摘要 雙向漸進法是近年來結構拓撲優化領域的一個熱點，針對現有結構拓撲優化方法容易出現的網格依賴性和局部極值等數值不穩定現象，以及計算速率慢的問題，提出一種雙重過濾增刪的雙向漸進法，簡稱為DBESO方法。該方法通過對典型的Michell結構、短懸臂結構和MBB梁結構進行不同過濾增刪次數的計算，得出二次增刪過濾為最佳的過濾增刪次數，即在執行靈敏度過濾和單元的增刪的基礎上，將結果存入中間變量再對靈敏度過濾及單元的增刪進行計算，并且設定一個判定結構性能的指標進行收斂。在MATLAB編程環境下，通過二維短懸臂結構和三維簡支梁結構的拓撲優化算例，證明DBESO方法在保證力學性能不變的情況下，有效地提高結構輕量化的收斂速率，更好地抑制雙向漸進方法的網格依賴性和局部極值等數值不穩定現象。

關 鍵 詞 雙向漸進法;輕量化;中間變量;靈敏度;MATLAB

中圖分類號 TB21? ? ?文獻標志碼 A

Lightweight design method based on improved BESO method

WANG Ying， ZHANG Zhengyan， LEI Yanyang， TANG Yu， SUN Lixin

Abstract The Bi-directional evolutionary structural optimization method is a hot topic in the field of structural topology optimization. To solve the problem of slow iteration rate and numerical instability such as mesh-dependency and local minima in the existing structural topology optimization method， the paper proposes a Bi-directional evolutionary structural optimization method of double sensitivity filter and double addition and deletion （ DBESO）. In this method， through the different addition and deletion filtering of typical Michell structure， short cantilever structure and MBB beam structure were calculated. It is concluded that the second addition and deletion filtering is the best filtering addition and deletion times. The sensitivity filter and the element addition and deletion are calculated again through an intermediate variable after the implementation of the sensitivity filter and the element addition and deletion. And set an index to determine the performance of the structure of the convergence condition. By the topology optimization example of two-dimensional short cantilever beam and three-dimensional simply supported beam structure in MATLAB programming environment，verify under the condition of ensuring the same mechanical performance， the convergence rate of the lightweight structure is effectively improved and the numerical instability of the Bi-directional evolutionary structural optimization method is preferably eliminated.

Key words BESO; lightweight; intermediate variable; sensitivity; MATLAB

結構輕量化是目前研究輕量化方法的主要方向之一，其最大的分支之一是連續體結構輕量化。由于連續體結構廣泛應用于工程中，因此對該理論的研究也較為廣泛[1]。連續體結構輕量化主要包括尺寸優化、形狀優化和拓撲優化。在尺寸和形狀優化中，主要對截面特性和設計域邊界進行優化。拓撲優化除了對設計域邊界進行優化外，還規定了連續體結構中空腔的數量、尺寸和位置，并在工業領域，例如航空航天和汽車領域都存在有巨大的潛力[2-4]。近20年來連續體結構拓撲優化獲得快速發展，成為工程創新設計的強有力工具。各種拓撲優化方法，例如密度法/各向同性實體材料懲罰法（solid isotropic material with penalization， SIMP）[5]、漸進結構優化法（evolutionary structural optimization， ESO）[6]、獨立連續映射法 （independent continuous mapping method， ICM）[7]，基于幾何邊界描述的水平集方法（level set method， LSM）[8]，以及最近發展的移動可變形組件法（moving morphable components， MMC）[9]，都用來解決各種結構拓撲優化設計問題[10]。

其中，雙向漸進結構優化方法（BESO）是在ESO方法和AESO方法綜合的基礎上提出的，不僅可以刪除靈敏度值較低的單元，而且在結構的靈敏度值較高區域周圍可以適當的添加單元，完善了ESO方法移除多余的單元后不能添加的缺陷。在近些年的BESO算法研究中，許多學者以目標函數靈敏度為基礎，從不同的角度對單元靈敏度進行了改進[11]。Huang等[12]采用了考慮靈敏度歷史迭代信息的方式，提高了靈敏度計算的精度。Huang等[13]利用插值策略提出了帶有懲罰參數的靈敏度修正方法。Huang等[14-15]和Nguyen等[16]分別提出了基于增加位移約束和基于非線性材料的靈敏度計算公式。Ghabraie等[16]通過提高泰勒級數更準確地計算了單元靈敏度。孫艷想[17]研究了一種平滑處理技術的雙向漸進結構拓撲優化。然而對于BESO方法容易出現的網格依賴性、局部極值等數值不穩定現象以及迭代速率慢的問題，目前還沒有得到良好的解決。

綜上所述，針對BESO優化算法計算過程中會導致數值不穩定現象，增刪過程中易出現誤差以及迭代速率較慢的問題，基于平滑處理技術的雙向漸進法（SBESO），提出雙重過濾增刪的雙向漸進方法（DBESO）。該方法是在靈敏度過濾和增刪后重復進行靈敏度過濾和增刪，并且設定一個測試結構性能的指標來進行收斂迭代，從而提高算法迭代的運算速率，消除數值不穩定現象，保障優化結果的穩定性和正確性。

1 優化模型的建立及靈敏度分析

1.1 數學模型的建立

改進的DBESO方法是在體積約束條件下，以最大靜剛度為目標函數，優化問題的數學模型描述為

式中：[f]是結構的載荷矢量;[u]是結構的位移矢量;[C]是結構的平均靜柔順度;[V?]是結構的約束體積;[Vi]是每個單元的體積;[N]是結構離散后總的單元個數;[xi]是設計變量;[ηj（TGS）]是權重函數，選用的是線性函數[ηj（TGS）=αj+β]，將[TGS]區域的設計變量分劃為0-1之間連續的數值，實現了算法計算過程中的平滑性和穩定性。

1.2 靈敏度分析

基于材料插值模型，使目標函數兩邊同時對[xi]進行求導，推導得出靈敏度計算公式為

式中，p表示為懲罰因子。當xi=1時，靈敏度數值的計算與p沒有關系;當[xi=1×ηj（TGS）]時，p越大，單元的靈敏度數值也越接近于0;當xi=0時，靈敏度值為0。

2 DBESO方法的改進策略

2.1 靈敏度過濾方法

由于靈敏度過濾方法易實現且具有高效率性，所以本文的濾波技術采用靈敏度過濾方法。首先計算單元節點靈敏度的數值，節點靈敏度可以表示為

式中：[M]是與[j]節點有關的單元總數目;[ωi]表示[i]單元的權重因子，可以用下式表示：

式中，[rij]表示[i]單元中心到[j]節點的距離。當單元[i]越接近節點[j]時，[ωi]的數值越大，說明[i]單元的靈敏度對[j]節點的靈敏度影響越大。

然后，將節點的靈敏度值轉換為單元的靈敏度值，設定某一單元為[i]單元，其靈敏度是通過以[i]單元中心為圓心，[rmin]為過濾半徑的區域Ω內節點靈敏度的計算而獲得的。Ω區域不隨網格的大小的變化而變化，并且Ω區域應該包含至少一個單元。Ω區域的每個節點靈敏度都會影響[i]單元的靈敏度，因此，[i]單元的靈敏度可以表示為

式中：[K]表示Ω區域中的節點數;[ωij]代表權重因子。

2.2 單元的增加和刪除策略

本文采用的增加刪除策略是基于SBESO方法進行不斷地改進而得到的。首先，對所有單元的靈敏度進行排序，劃分成[n]個小組。然后根據靈敏度的大小設置閥值[αthdel]和[αthadd]，按照閥值把結構所有單元劃分在3個區域[T]，[TGS]和[TLS]。如果滿足[αi≤αthdel]，該單元歸為[TLS]區域，如果滿足[αi>αthadd]，該單元歸為[T]區域，如果滿足[αthdel<αi≤αthadd]，該單元歸為[TGS]區域。最后，對[T]區域中的單元進行單元的添加，[TLS]區域中的單元完全移除，[TGS]區域的單元按照權重函數[ηj（TGS）]進行修改再判定是否增刪。即將單元靈敏度劃分成[n]個小組，定義[n]組中的p%單元完全被移除和增加，（1-p%）單元的移除和重新返回通過一個線性函數來實現如圖1所示，單元密度的高低用一個權重函數[η（j）]來表示，[0≤η（j）≤1]。該函數用以判斷[TGS]區域所有單元靈敏度的重要程度，當單元的靈敏度越大，那么對應的權重函數也就越大。這就表明，在下一次迭代中這些單元有可能會重新恢復到結構中，而當單元靈敏度越小，對應的權重函數就越小，下一次迭代中這些單元可能完全從結構中移除。

2.3 DBESO方法的改進策略

靈敏度過濾方法是當前處理棋盤格現象、局部極值和網格依賴性等數值不穩定現象的最有效手段，但依舊無法根除數值不穩定現象，只能在抑制該現象方法的基礎上進行不斷改進，更好地抑制該現象的出現。而單元的增刪是對單元靈敏度進行排序，進而進行單元的增刪。由于需要不斷地添加和刪除，所以出現誤差的可能性較高。

本文針對上述問題，提出在一次迭代中進行多次過濾增刪的方法，用以加強迭代過程中對數值不穩定現象的抑制作用，減少單元增刪的誤差。

由于進行多次過濾增刪過多次會引起迭代速率減慢，局部極值變多的現象，需要通過研究設定最優過濾增刪次數。由此設計1次、2次、3次等多次過濾增刪，對其進行基于SBESO方法的優化迭代得到的結果進行對比，選用典型的Michell結構、短懸臂結構和MBB梁結構，幾何尺寸分別為[L×D=96? mm×48? mm]、[L×D=80 mm×50 mm]和[L×D=120 mm×60 mm]，[F=100 N]，具體模型和作用力位置和方向如圖2所示，其他初始條件不變，對3個數值算例的3種優化結果總結如表1所示。

通過表1可以看出，對Michell結構、短懸臂梁結構和MBB梁結構進行優化，最優結果隨著靈敏度過濾增刪次數的不斷增加，網格依賴性的抑制作用明顯得到了加強。通過3種數值算例的平均柔順度變化曲線可以看出，在2次靈敏度過濾增刪的情況下局部極值問題較好的得以解決。對于平均柔順度值來說，隨著靈敏度過濾增刪次數的增多而越來越大，2次過濾增刪的情況下柔順度值還較為穩定，與1次靈敏度過濾增刪相似，但3次靈敏度過濾增刪的情況下，平均柔順度值差距較大。并且，通過計算運算時間，可以得出單次運行時間在2次過濾增刪時比1次過濾增刪時增大了1%～6%，在3次過濾增刪時比1次過濾增刪時增大了6%～10%，可以看出，2次過濾增刪對計算速率影響較小，所以選用2次過濾增刪為最優的過濾增刪次數。

綜上所述，選用雙重過濾增刪的策略，設定中間變量[x?]，用來儲存第1次過濾增刪后的變量情況;然后將中間變量[x?]代入靈敏度計算公式中，按步驟進行第2次過濾增刪;最后繼續進行優化迭代的其他步驟。

2.4 收斂準則

雖然雙重過濾增刪會使其在迭代中數值不穩定現象得到很好的抑制，但在收斂速度方面并沒有提高，所以需定義一個結構性能指標，來監控性能是否滿足要求并且可以提高迭代速率。對于優化算法來說，由于體積和柔順度是結構性能好壞的展現，所以，定義了一個和體積與柔順度2個參數有關的性能指標來對結構性能進行監控。結構性能參數I的表達式如下：

式中：[Vi]是第[i]個迭代步驟結構的體積;[Vk+1]是第[i+1]個迭代步驟結構的體積。可以看出，[I]值越大，結構的性能越好。所以，用[I]來對該算法進行收斂。得出收斂準則為

式中：[Ii]是當前迭代步的性能指標;[Ii-1]是上一迭代步的性能指標;error是收斂因子。在結構滿足式（8）的收斂準則時，結構獲得最優拓撲結構，停止迭代。

3 改進DBESO方法優化流程

圖3為雙重過濾增刪BESO（DBESO）方法的迭代流程圖。由圖3可以得出改進的雙重增刪過濾BESO（DBESO）優化方法的迭代流程為：

1）設定初始參數，建立有限元模型，設定負載和邊界約束條件;

2）采用靈敏度公式計算每一個單元的靈敏度，進行靈敏度過濾和更新;

3）按照靈敏度的大小設定增刪閥值，然后將單元依據靈敏度大小劃分在3個區域T，TLS，TGS;對T區域中的單元進行添加，對TLS區域中的單元進行移除，按照權重函數[ηj（TGS）]對TGS區域的設計變量進行增刪判定;

4）將增刪結果設為中間變量[x?];

5）用[x?]進行靈敏度再過濾和更新，然后進行單元的再增刪;

6）計算下一步迭代體積并更新單元設計變量;

7）計算性能指標，判定結構性能。

8）重復步驟2-7，在結構達到目標體積和滿足收斂準則時，終止迭代。

4 改進DBESO方法數值算例驗證

4.1 二維數值算例短懸臂梁結構

定義優化設計模型為二維短懸臂梁結構，其幾何參數為[L×D=80? mm×50? mm]，材料的特性為楊氏模量為[2×105? MPa]，泊松比為0.3。邊界約束條件和載荷如圖2b）所示，左端固定，右端中點承受向下恒定載荷[F=100? N]。對結構采用80×50四節點有限元平面應力單元來離散結構，采用的權重函數為線性函數[ηj（TGS）=αj+β]。初始參數為：[V?=0.5]，[er=0.02]，[rmin=5 mm]，[penal=3]，[p=0.8]。

本文將該短懸臂梁進行120×76以及160×100的網格劃分，結果如圖4所示，對比圖中結果，對于改進后DBESO方法網格數為80×50、120×76以及160×100的計算出的優化結果的網格是一樣的，所以改進后的方法明顯地減少了網格依賴性。對于改進前SBESO方法網格數為120×76以及160×100計算出的結果明顯比網格數80×50計算出的結果網格數要多，網格依賴性嚴重，并且隨著網格劃分數的增大，優化結果的網格數越來越多。由此可以看出，改進后DBESO方法對于抑制網格依賴性等數值不穩定現象有了明顯的提高。

為了驗證改進后DBESO算法的優點，本文將基于改進前平滑的SBESO方法和改進后雙重增刪過濾的DBESO方法對該算例進行計算。在材料參數和約束條件不變的情況下，2種方法最終優化結果如圖5和圖6所示。并將2種方法結果以及不同網格劃分下的結果進行總結對比，如表2所示。

觀察圖4、圖5、圖6和表2，可以得出以下結論。

1）對于平均柔順度而言，在相同網格劃分的情況下，改進前后2種方法的最終柔順度值相近，改進前的平均柔順度曲線有明顯的數值波動，數值穩定性較差。改進后的平均柔順度曲線比改進前的更為平滑，無明顯的數值波動，較為平穩。

2）對于計算效率來說，在保持平均柔順度相近的情況下，對于80×50網格劃分的情況，改進后DBESO方法迭代速率提高了29%;對于120×76網格劃分情況，改進后DBESO方法迭代速率提高14.4%;對于160×100網格劃分情況，改進后DBESO方法迭代速率提高14.7%。

3）對于網格依賴性來說，相對于改進前優化方法，改進后方法在抑制網格依賴性方面有了明顯的提高，減小了對零件的性能影響。

對于二維結構來說，改進的DBESO方法在提高計算效率方面以及抑制局部極值和網格依賴性等數值不穩定現象均有明顯的優化。

4.2 三維數值算例簡支梁結構

定義優化簡支梁的尺寸為[50? mm×20? ?mm×4? ?mm]。模型如圖7所示，其底面4個頂點為簡支約束，頂面正中心施加垂直向下均布載荷[F=10 N]，密度[ρ=1 g/mm3]。將設計域劃分為50×20×4的網格，材料特性有：[E=1? GPa]，泊松比[μ=3.0]，所需的參數為：[V?=0.5]，[er=0.02]，[rmin=3 mm]，[penal=3]，[p=0.9]。

本文將該短懸臂梁進行66×26×6以及60×24×4的網格劃分，結果如圖8所示，對比圖中結果，對于改進后DBESO方法網格數為50×20×4、66×26×6以及60×24×4的計算出的優化結果的網格是一樣的，所以改進后的方法明顯地減少了網格依賴性。對于改進前SBESO方法網格數為66×26×6以及60×24×4計算出的結果明顯比網格數50×20×4計算出的結果網格數要多，網格依賴性嚴重，并且隨著網格劃分數的增大，優化結果的網格數增多。由此可以看出，改進后DBESO方法對于抑制網格依賴性等數值不穩定現象有了明顯的提高。

在MATLAB編寫好的改進優化算法中，首先編寫好該算例的有限元模型、邊界條件、受力情況和位移約束情況，然后對改進前的SBESO方法和改進后的DBESO方法執行優化計算，得到改進前后最優結果圖如圖8所示，改進前后的平均柔順度與迭代次數的曲線圖如圖9所示，改進前后的性能指標與迭代次數的曲線圖如圖10所示，對在不同網格劃分下的兩種方法結果總結如表3所示。

觀察圖8、圖9、圖10和表3，可以得出以下結論。

1）對于平均柔順度而言，在相同網格劃分的情況下，改進前后方法的最終柔順度值相近，改進前的平均柔順度曲線有明顯的數值波動，數值穩定性較差，改進后的平均柔順度曲線比改進前的更為平滑，無明顯的數值波動，較為平穩。

2）對于計算效率來說，在保持平均柔順度相近的情況下，對于50×20×4網格劃分的情況，改進后DBESO方法迭代速率提高了1.8%;對于60×24×4網格劃分情況，改進后DBESO方法迭代速率提高0.6%;對于66×26×6網格劃分情況，改進后DBESO方法迭代速率提高9.5%。

3）對于性能指標來說，改進后的曲線比改進前的更為平滑，無明顯的數值波動，較為平穩，改進前的曲線有明顯的數值波動，數值穩定性較差，且迭代前后性能指標最終值相近。

4）對于網格依賴性來說，相對于改進前SBESO方法，改進后DBESO方法在抑制網格依賴性方面有了明顯的提高。

對于三維結構來說，改進的DBESO方法在抑制局部極值和網格依賴性等數值不穩定現象有了明顯的優化，且計算速率也有了小的提高。

5 結論

本文在現有SBESO方法的基礎上，提出了雙重過濾增刪的DBESO方法。該方法是在原有過濾增刪的基礎上經過計算驗證后，將結果存入中間變量，再進行雙重過濾增刪，并且在收斂時設定1個對結構進行性能判定的指標來進行收斂，用于提高收斂速率。通過二維短懸臂和三維簡支梁結構作為數值算例進行計算，得出結果與現有的SBESO方法進行對比，驗證了該DBESO方法的優越性和可行性。

1）對于網格依賴性問題，相較于SBESO方法明顯的網格依賴性問題，可以看出改進后的DBESO方法更好的抑制了網格依賴性的問題。

2）對于局部極值現象，相較于SBESO方法的大的局部極值波動，改進后DBESO方法的迭代曲線更為平順，幾乎不存在局部極值問題。

3）對于計算結果與效率問題，相較于SBESO方法，改進后的DBESO方法的計算效率有了明顯提高。

4）對于應用的廣泛性來說，在二維和三維結構的優化中，改進的DBESO方法比SBESO方法在提高迭代速率和抑制網格依賴性、局部極值等數值不穩定問題方面有了明顯的優化。

結果證明了改進的DBESO方法在二維結構和三維結構都明顯優于原算法，更好地解決了數值不穩定導致的網格依賴性和局部極值等問題，且在迭代效率方面也有了提高。

參考文獻：

[1]? ? 宋武林. 三維連續體雙向漸進結構拓撲優化[D]. 哈爾濱：哈爾濱工程大學，2017.

[2]? ? 徐文鵬，苗龍濤，劉利剛. 面向3D打印的結構優化研究進展[J]. 計算機輔助設計與圖形學學報，2017，29（7）：1155-1168.

[3]? ? LEDERMANN C，ERMANNI P，KELM R. Dynamic CAD objects for structural optimization in preliminary aircraft design[J]. Aerospace Science and Technology，2006，10（7）：601-610.

[4]? ? SHOBEIRI V. Topology optimization using bi-directional evolutionary structural optimization based on the element-free Galerkin method[J]. Engineering Optimization，2016，48（3）：380-396.

[5]? ? BENDS?E M P，SIGMUND O. Topology Optimization：Theory，Methods，and Applications[M]. Berlin，Heidelberg：Springer，2004.

[6]? ? XIE Y M，STEVEN G P. A simple evolutionary procedure for structural optimization[J]. Computers & Structures，1993，49（5）：885-896.

[7]? ? 隋允康，葉紅玲，彭細榮，等. 連續體結構拓撲優化應力約束凝聚化的ICM方法[J]. 力學學報，2007，39（4）：554-563.

[8]? ? WANG M Y，WANG X M，GUO D M. A level set method for structural topology optimization[J]. Computer Methods in Applied Mechanics and Engineering，2003，192（1/2）：227-246.

[9]? ? GUO X，ZHANG W S，ZHONG W L. Doing topology optimization explicitly and geometrically—A new moving morphable components based framework[J]. Journal of Applied Mechanics，2014，81（8）：081009.

[10]? SIGMUND O，MAUTE K. Topology optimization approaches[J]. Structural and Multidisciplinary Optimization，2013，48（6）：1031-1055.

[11]? 劉娟. SIMP算法和BESO算法的關鍵技術研究[D]. 桂林：桂林電子科技大學，2016.

[12]? HUANG X，XIE Y M. Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[J]. Finite Elements in Analysis and Design，2007，43（14）：1039-1049.

[13]? HUANG X，XIE Y M. Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials[J]. Computational Mechanics，2008，43（3）：393-401.

[14]? HUANG X，XIE Y M. Evolutionary Topology Optimization of Continuum Structures[M]. Chichester，UK：John Wiley & Sons，Ltd，2010.

[15]? HUANG X D，XIE Y M. A further review of ESO type methods for topology optimization[J]. Structural and Multidisciplinary Optimization，2010，41（5）：671-683.

[16]? NGUYEN T，GHABRAIE K，TRAN-CONG T. Applying bi-directional evolutionary structural optimisation method for tunnel reinforcement design considering nonlinear material behaviour[J]. Computers and Geotechnics，2014，55：57-66.

[17]? 孫艷想. 基于BESO方法的結構動力學拓撲優化研究[D]. 哈爾濱：哈爾濱工程大學，2017.

收稿日期：2019-03-23

基金項目：國家自然科學基金（61802108）;河北省自然科學基金（E2017202296）

第一作者：王穎（1993—），女，碩士研究生。通信作者：孫立新（1964—），男，教授，Sunlixin100@126.com。