王 曦，趙重年，王文強

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

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

1 拓撲優(yōu)化基本理論

1.1 變密度法

結(jié)構(gòu)優(yōu)化設計是用系統(tǒng)的、目標定向的過程與方法代替?zhèn)鹘y(tǒng)設計，其目的在于尋求既經(jīng)濟又適用的結(jié)構(gòu)形式，以最少材料、最低造價實現(xiàn)結(jié)構(gòu)的最佳性能。一般地，優(yōu)化問題的數(shù)學模型可表示為：

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

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

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

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

1.2 插值模型

變密度法定義一個經(jīng)驗公式（式）來表達相對密度與彈性模量E的非線性關系，通過這種關系使設計變量與目標函數(shù)產(chǎn)生聯(lián)系，則結(jié)構(gòu)拓撲優(yōu)化問題轉(zhuǎn)化為材料的最優(yōu)分布問題。

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

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

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

2 COMSOL中的拓撲優(yōu)化

2.1 優(yōu)化目標與約束條件的設置

在考慮結(jié)構(gòu)優(yōu)化的問題時，需要特別注意，對復雜工程結(jié)構(gòu)的拓撲優(yōu)化是一個計算量相當大的問題，因此，在軟件中實現(xiàn)優(yōu)化時要進行模型簡化和等效條件轉(zhuǎn)化。

（1）目標函數(shù)

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

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

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

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

優(yōu)化過程中，定義當前結(jié)構(gòu)應變能與初始應變能的比值Ws/Ws0為目標函數(shù)，則比值最小表示結(jié)構(gòu)剛度最大。

（2）約束條件

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

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

2.2 拓撲優(yōu)化中的密度模型及不穩(wěn)定現(xiàn)象的處理

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

此外，在拓撲優(yōu)化中，會出現(xiàn)一些與模型本身無關的數(shù)值不穩(wěn)現(xiàn)象，主要有中間密度單元，棋盤格式現(xiàn)象，網(wǎng)格依賴性等。

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

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

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

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

2.3 參數(shù)對優(yōu)化結(jié)果的影響

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

圖2 MBB梁幾何模型

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

表格1 懲罰因子對MBB梁拓撲優(yōu)化結(jié)果的影響

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

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

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

圖3 β取不同值時對MBB梁拓撲優(yōu)化結(jié)果的影響

圖4 β取不同值時對MBB梁拓撲優(yōu)化的求解迭代次數(shù)

3 聯(lián)合仿真

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

3.1 對一般結(jié)構(gòu)優(yōu)化的實現(xiàn)

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

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

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

3.2 幾種結(jié)構(gòu)的測試

為了檢驗程序的可行性和通用性，通過幾種典型的三維結(jié)構(gòu)測試程序的優(yōu)化結(jié)果。

3.2.1 L型負載梁

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

圖5 L型負載梁結(jié)構(gòu)及拓撲優(yōu)化結(jié)果

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

3.2.2 簡單懸臂梁

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

圖6 懸臂梁結(jié)構(gòu)及拓撲優(yōu)化結(jié)果

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

由上兩例優(yōu)化結(jié)果可看出，本程序拓撲優(yōu)化的結(jié)果為整體結(jié)構(gòu)的概念形狀，對于邊界的處理較為粗糙，雖然不適于實際的生產(chǎn)加工，但其作為工程設計的初始結(jié)構(gòu)是可以接受的?？梢砸源藘?yōu)化結(jié)果為基礎，建立參數(shù)化模型，根據(jù)具體的載荷大小和實際工藝，對結(jié)構(gòu)進行形狀優(yōu)化和尺寸優(yōu)化。

4 結(jié)束語

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