基于COMSOL的拓撲優化及聯合仿真

王 曦，趙重年，王文強

（1．陸軍軍事交通學院研究生隊，天津300161；2．陸軍軍事交通學院軍事交通運輸研究所，天津300161；3．陸軍軍事交通學院，天津 300161）

拓撲優化是根據負載情況、約束條件（如應力、位移、頻率和重量等）和性能指標（剛度、重量等），利用有限元分析和優化方法，使設計域達到最優材料布局的一種結構優化方法[1]。目前，連續體拓撲優化的研究已經較為成熟，現有的商業有限元軟件基本都可以進行拓撲優化，但每款軟件的優化模塊各有區別，關于參數的設置也多種多樣，通用性較差，同時，對用戶的使用要求較高，不適于拓撲優化的利用和推廣，一定程度上遲滯了結構優化的發展，因此，提高拓撲優化的在仿真中的通用性和工程應用中的應用效率對結構優化具有重要意義。

1 拓撲優化基本理論

1.1 變密度法

結構優化設計是用系統的、目標定向的過程與方法代替傳統設計，其目的在于尋求既經濟又適用的結構形式，以最少材料、最低造價實現結構的最佳性能。一般地，優化問題的數學模型可表示為：

式中：X為設計變量，f（X）為目標函數，gi（x），hj（x）為約束條件。

假設設計域由有限個單元組成，變密度法的基本思想就是定義一種密度可變的材料單元填充設計域，取值定義在（0，1）區間內，取0表示孔洞，取1表示實體，則優化問題轉變為0-1問題。

在拓撲優化中，希望優化后的結構在保持最大應力承受極限的情況下質量減少，對于大多數工程上的零件，均采用統一結構制造，則質量最小可以等價為體積最小。優化的目標一般是結構的剛度，為了在數學上方便凸優化計算，將剛度最大等效為柔度最小，其模型表達為：

式中：C（x）為柔度，V 為設計域初始體積，f為體積約束因子。

1.2 插值模型

變密度法定義一個經驗公式（式）來表達相對密度與彈性模量E的非線性關系，通過這種關系使設計變量與目標函數產生聯系，則結構拓撲優化問題轉化為材料的最優分布問題。

式中，f（xi）又稱為材料插值模型，設計變量xi控制著材料的取舍，為了更好的獲得拓撲結構的形狀和邊界，插值模型應使xi的取值向區間兩端分化。目前，應用最為廣泛的兩種插值模型為SIMP模型[2]（式）和RAMP[3]模型（式）：

兩種模型都是通過引入懲罰系數，使得變量更快的逼近端點值，前者采用指數函數的形式，引入懲罰因子p，后者采用有理數函數的形式，引入權重因子q。兩者性能由圖1可知。

圖1 兩種材料插值模型懲罰效果圖

2 COMSOL中的拓撲優化

2.1 優化目標與約束條件的設置

在考慮結構優化的問題時，需要特別注意，對復雜工程結構的拓撲優化是一個計算量相當大的問題，因此，在軟件中實現優化時要進行模型簡化和等效條件轉化。

（1）目標函數

拓撲優化的目標函數一般為剛度最大，即柔度最小。用結構在極限條件下的總變形能來表示結構剛度最大時的狀態。應變能可以精確的平衡載荷所做的功，因此，應變能最小可以使施加載荷所引起的位移最小，有效地得到柔順性最小的目標。

通過有限元知識可知，結構的總變形能為各個單元變形能量的總和。對由微小體元dΩ=dxdydz組成的結構體，在受到各個方向的正應力和切應力時的應變能為：

COMSOL中對于應變能的描述體現為指標形式，如公式：

式中：solid.Sl為應力，solid.eel為應變。

優化過程中，定義當前結構應變能與初始應變能的比值Ws/Ws0為目標函數，則比值最小表示結構剛度最大。

（2）約束條件

從概念上講，拓撲優化只是最小化所用材料的數量，以初始設計域的體積分數來定義約束條件為：

式中：f為體積分數，V為初始設計域的體積。

2.2 拓撲優化中的密度模型及不穩定現象的處理

在COMSOL中的拓撲優化模塊中，建立了變密度法為基礎的密度模型，用戶可以根據需要選擇SIMP或RAMP材料插值模型。設計變量相對密度xi用材料體積因子θ表示，即每個單元格實體材料的體積占比。

此外，在拓撲優化中，會出現一些與模型本身無關的數值不穩現象，主要有中間密度單元，棋盤格式現象，網格依賴性等。

對于棋盤格現象和網格依賴性，較為常用的方法是敏度過濾法，通過提取相鄰單元信息，對當前單元的取值產生影響，距離越近，影響越大。但是對相鄰單元提取信息的計算量相當龐大，COMSOL中采用Sigmund提出的亥姆霍茲過濾器作為敏度過濾的改進方案[4]，僅需要提取所需單元的網格信息，用最小過濾半徑與體積因子的拉普拉斯算子的乘積對材料體積因子進行過濾，提高了計算效率。

式中：θf表示過濾后的材料體積因子；θc表示系統控制的材料體積因子；Rmin表示最小過濾半徑，一般取最小網格的1/2或1/4。

將過濾后的設計變量進行投影，可以減少中間密度單元[5]。選擇基于雙曲正切函數的投影可以得到較為清晰的拓撲圖像（如公式10）。

式中：β表示投影斜率，θβ表示投影點。

2.3 參數對優化結果的影響

本文通過二維MBB梁結構在的拓撲優化討論懲罰因子p和投影斜率β對優化結構的影響。如圖2為MBB梁結構，矩形梁長6 m，高1 m，兩端支撐，中間施加垂直向下300 kN的力，結構材料為普通鋼鐵，設置體積因子上限為0.5。

圖2 MBB梁幾何模型

相同條件下，表格1為懲罰因子p取2～5時的優化結果對比。

表格1 懲罰因子對MBB梁拓撲優化結果的影響

由表中圖像可知，懲罰因子的增大可以控制灰度單元，得到更加清晰的拓撲邊界。也就是說，p值可以取到足夠大，以獲得足夠清晰的邊界。但同時迭代次數也相應增多，計算時間變長，收斂速度減慢，兩者相互矛盾。工程中，一般取p=3，使得優化結果較為清晰的同時也能保證計算效率。

β作為雙曲正切投影的關鍵參數，控制著結果整體的清晰度，在相同條件下，分別取β為1～12時的結果進行評估，圖3選取了5個取值的優化結果，圖4為迭代次數變化趨勢。

由上圖表可知，受雙曲正切函數性質的影響，隨著投影斜率的增大，收斂速度加快，圖像也更加清晰，但當β大于8時，迭代次數趨于平穩，不再減少，且結構出現奇異，不能得到正確的拓撲結構，因此，β一般取8最為合適。

圖3 β取不同值時對MBB梁拓撲優化結果的影響

圖4 β取不同值時對MBB梁拓撲優化的求解迭代次數

3 聯合仿真

由上文可知，對于COMSOL中的拓撲優化大致分為兩個步驟，第一步是對原結構的分析，得到初始結構的應變能參考值，第二步才是拓撲優化的相關設置。繁瑣的二次設置勢必會影響工作效率，因此，如果將拓撲優化的關鍵步驟集成在一個程序內，那么對于一般的結構，只要完成其相應的前處理，就可以方便地實現優化。

3.1 對一般結構優化的實現

COMSOL源自于MATLAB的工具箱（原FEMLAB），通過其自身的APL語法，將從建模開始的所有步驟編入MATLAB中，再與控制語句相結合，實現模型的各種后處理。

首先，在COMSOL界面中完成初始結構模型的前處理，包括幾何模型搭建，材料、邊界條件和網格的設置，保存為.m文件，通過MATLAB控制語句調節程序流程，對初始結構進行有限元分析。然后，將第一步分析得到的結構應變能提取出來，作為參數放入第二步優化中。最后，在MATLAB中調試好優化求解器進行求解。本程序設置體積約束、懲罰因子、投影點、投影斜率、優化域及非優化域等6個參數，也可以根據需要進行添加或減少。

對于優化的結果，可以在MATLAB的圖片查看器中分析，也可以將結果連接到COMSOL界面中進行處理。

3.2 幾種結構的測試

為了檢驗程序的可行性和通用性，通過幾種典型的三維結構測試程序的優化結果。

3.2.1 L型負載梁

如圖5，三維結構L型梁長寬比為0.9，厚度為0.1 m。梁的上端固定，右側施加一個垂直向下的載荷，材料為鋼鐵，體積約束為0.3。其優化結果如圖5所示。

圖5 L型負載梁結構及拓撲優化結果

由優化結果可看出，在L型梁的上半部分只保留了兩側的受力懸臂，下半部分得到3個大小不等孔洞，其內部基本掏空，只保留了部分連接兩個表面的結構。拓撲結構較為清晰，邊界明顯，孔洞雖然較為規律，仍有優化空間，整體來看，是比較成功的優化結果。

3.2.2 簡單懸臂梁

如圖6所示，三維結構懸臂梁，長0.8 m，高0.5 m，厚0.1 m。懸臂梁一端固定，另一端施加一個垂直向下的變載荷，材料為鋼鐵，體積約束為0.5。其優化結果如圖6所示。

圖6 懸臂梁結構及拓撲優化結果

由優化結果可看出，梁的右上角由于受力較小，全部挖空，中間掏出一個三角形的孔洞，整體呈現桁架結構，符合實際物理情況，同樣是比較成功的優化結果。

由上兩例優化結果可看出，本程序拓撲優化的結果為整體結構的概念形狀，對于邊界的處理較為粗糙，雖然不適于實際的生產加工，但其作為工程設計的初始結構是可以接受的。可以以此優化結果為基礎，建立參數化模型，根據具體的載荷大小和實際工藝，對結構進行形狀優化和尺寸優化。

4 結束語

本文研究了基于變密度法的拓撲優化理論及在COMSOL軟件中的應用，討論研究了幾種重要參數對優化結構的影響，得出了最佳的設置方案，并利用MATLAB將COMSOL中繁瑣的操作集成到一個程序中，開發了可以對一般結構實現拓撲優化的通用程序，從而為工程設計人員提供可靠的、指導性的結構設計方案，避免了復雜繁瑣設計步驟，提高了程序的普適性，縮短了設計周期，對工程設計有一定促進意義。