基于COMSOL的拓?fù)鋬?yōu)化及聯(lián)合仿真

王 曦，趙重年，王文強(qiáng)

（1．陸軍軍事交通學(xué)院研究生隊(duì)，天津300161；2．陸軍軍事交通學(xué)院軍事交通運(yùn)輸研究所，天津300161；3．陸軍軍事交通學(xué)院，天津 300161）

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

1 拓?fù)鋬?yōu)化基本理論

1.1 變密度法

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

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

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

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

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

1.2 插值模型

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

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

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

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

2 COMSOL中的拓?fù)鋬?yōu)化

2.1 優(yōu)化目標(biāo)與約束條件的設(shè)置

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

（1）目標(biāo)函數(shù)

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

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

COMSOL中對(duì)于應(yīng)變能的描述體現(xiàn)為指標(biāo)形式，如公式：

式中：solid.Sl為應(yīng)力，solid.eel為應(yīng)變。

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

（2）約束條件

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

式中：f為體積分?jǐn)?shù)，V為初始設(shè)計(jì)域的體積。

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

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

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

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

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

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

式中：β表示投影斜率，θβ表示投影點(diǎn)。

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

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

圖2 MBB梁幾何模型

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

表格1 懲罰因子對(duì)MBB梁拓?fù)鋬?yōu)化結(jié)果的影響

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

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

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

圖3 β取不同值時(shí)對(duì)MBB梁拓?fù)鋬?yōu)化結(jié)果的影響

圖4 β取不同值時(shí)對(duì)MBB梁拓?fù)鋬?yōu)化的求解迭代次數(shù)

3 聯(lián)合仿真

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

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

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

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

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

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

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

3.2.1 L型負(fù)載梁

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

圖5 L型負(fù)載梁結(jié)構(gòu)及拓?fù)鋬?yōu)化結(jié)果

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

3.2.2 簡(jiǎn)單懸臂梁

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

圖6 懸臂梁結(jié)構(gòu)及拓?fù)鋬?yōu)化結(jié)果

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

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

4 結(jié)束語

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