基于變密度法的ABAQUS-MATLAB 集成拓?fù)鋬?yōu)化方法

靳立涵JIN Li-han

（重慶交通大學(xué)機(jī)電與車輛工程學(xué)院，重慶 400074）

0 引言

結(jié)構(gòu)優(yōu)化共有三種優(yōu)化方式分別為：尺寸優(yōu)化、形狀優(yōu)化和拓?fù)鋬?yōu)化。其中拓?fù)鋬?yōu)化是一種在有限的空間中根據(jù)設(shè)計(jì)要求確定材料最優(yōu)布局的數(shù)學(xué)方法。近三十年來，隨著計(jì)算機(jī)技術(shù)的發(fā)展，拓?fù)鋬?yōu)化技術(shù)也取得了巨大的進(jìn)步。相關(guān)的從業(yè)人員和研究者提出了許多的拓?fù)鋬?yōu)化算法。

拓?fù)鋬?yōu)化算法中有三個(gè)主流方法分別是：變密度方法、雙向漸進(jìn)結(jié)構(gòu)法[1-2]和水平集方法[3]。對于變密度法，Sigmund[4]在MATLAB 平臺(tái)上介紹了基于固體各向同性材料懲罰（SIMP）的99 行二維拓?fù)鋬?yōu)化程序，簡潔易讀，很好地闡述了密度法。Andreassen[5]等在99 行代碼的基礎(chǔ)上提出了88 行拓?fù)鋬?yōu)化代碼，88 行算法在優(yōu)化相同算例時(shí)比99 行算法快100 倍。Liu 和Tovar[6]還在MATLAB 平臺(tái)上公開了一個(gè)169 行拓?fù)鋬?yōu)化程序，在保證代碼簡單和優(yōu)化效率高的同時(shí)實(shí)現(xiàn)了三維目標(biāo)的優(yōu)化。Ferrari 等人[7]提出了一種新的99 行和125 行拓?fù)鋬?yōu)化代碼，分別處理二維和三維拓?fù)鋬?yōu)化問題，是目前MATLAB 平臺(tái)上最高效的密度法拓?fù)鋬?yōu)化代碼。

本文將MATLAB 與ABAQUS 腳本接口進(jìn)行集成用以拓?fù)鋬?yōu)化，該集成方法的目的是充分利用MATLAB 平臺(tái)的矩陣計(jì)算能力和商業(yè)有限元軟件友好的可視化界面。該方法通過一段Python 代碼和一段MATLAB 代碼實(shí)現(xiàn)。這兩個(gè)代碼基于變密度法開發(fā)，通過MATLAB 程序調(diào)用有商業(yè)限元軟件實(shí)現(xiàn)了對目標(biāo)設(shè)計(jì)域在體積約束下的最小合規(guī)性的拓?fù)鋬?yōu)化。

1 SIMP 拓?fù)鋬?yōu)化方法理論

本文關(guān)于ABAQUS-MATLAB 拓?fù)鋬?yōu)化方法研究的目的，為了得到一個(gè)結(jié)合兩個(gè)軟件優(yōu)點(diǎn)并且能適用于基于FEM（有限元法）的多種拓?fù)鋬?yōu)化方法的集成框架。基于變密度方法（SIMP）的特點(diǎn)：發(fā)展時(shí)間長較為成熟、易拓展和灰度顯示。因此本文選擇變密度法作為集成框架的理論基礎(chǔ)。

SIMP 法拓?fù)鋬?yōu)化模型：

為了將離散型問題轉(zhuǎn)化為連續(xù)性優(yōu)化問題，變密度法提出了一種引入中間密度單元的方法進(jìn)行計(jì)算。但在實(shí)際制造過程中，一個(gè)單元內(nèi)只有填充材料和不填充材料兩種狀態(tài)，中間密度的存在使得模型無法制造。為此變密度法還提出了一個(gè)通過對中間密度單元進(jìn)行懲罰的方式來避免產(chǎn)生中間密度單元，使中間密度單元在迭代計(jì)算后趨向于實(shí)體單元或非實(shí)體單元。

在本文中討論的SIMP 表達(dá)式為：

其中xi是編號(hào)為i 的單元的相對密度。

SIMP 法優(yōu)化的數(shù)學(xué)模型為：

本方案是求解在體積或者質(zhì)量約束下的優(yōu)化模型的最小柔度，最終得到在約束條件下的最大剛度的優(yōu)化模型。其中X 是單元相對密度矢量為單元設(shè)計(jì)變量，C 是結(jié)構(gòu)柔度，F(xiàn) 和U 分別是載荷矢量和位移矢量，ki為單元?jiǎng)偠染仃嚕琸0為初始單元?jiǎng)偠染仃嚕瑅i為單元體積，f 是保留的體積分?jǐn)?shù)，V0是初始體積，xmin是設(shè)計(jì)變量的下限值，xmax是設(shè)計(jì)單元的上限值，m 為最大單元數(shù)。

2 ABAQUS-MATLAB 集成方法及實(shí)現(xiàn)

圖1 是ABAQUS-MATLAB 集成框架示意圖，其中MATLAB 程序完成拓?fù)鋬?yōu)化模型求解，包括靈敏度過濾和設(shè)計(jì)變量更新；Python 程序由MATLAB 調(diào)用后會(huì)執(zhí)行ABAQUS 有限元計(jì)算內(nèi)核，實(shí)現(xiàn)每個(gè)迭代步驟中有限元求解。在運(yùn)行之前，CAE 文件、MATLAB 代碼和Python 代碼需要放在同一個(gè)文件夾中，其余文件由程序生成。

圖1 ABAQUS-MATLAB 集成方案

根據(jù)結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)要求，在ABAQUS/CAE GUI中建立有限元模型，模型包括網(wǎng)格單元、材料屬性、邊界條件和載荷應(yīng)用，其中材料屬性可以被更新。MATLAB 代碼完成求解初始化，包括有限元模型的單元總數(shù)和迭代的相關(guān)參數(shù)。隨后，在循環(huán)迭代開始時(shí)進(jìn)行拓?fù)鋬?yōu)化，Python 程序讀取MATLAB 代碼生成的各單元設(shè)計(jì)變量，根據(jù)（1）式的材料插值模型創(chuàng)建各單元對應(yīng)的楊氏模量，并更新有限元模型的材料屬性。

Python 程序調(diào)用ABAQUS 內(nèi)核，求解有限元模型，并從結(jié)果文件中提取相關(guān)參數(shù)。MATLAB 代碼根據(jù)優(yōu)化相關(guān)參數(shù)完成靈敏度過濾和設(shè)計(jì)變量更新，在獲取新的靈敏度后根據(jù)收斂準(zhǔn)則判斷是否收斂，如果不收斂則進(jìn)行下一次循環(huán)迭代，直到收斂。*.XLS 文件是MATLAB 和ABAQUS之間傳輸數(shù)據(jù)的介質(zhì)，該文件是通過代碼調(diào)用Microsoft Excel 生成的。

2.1 MATLAB 程序的實(shí)現(xiàn)

本文所編寫的MATLAB 程序結(jié)構(gòu)如圖2 所示。當(dāng)開始優(yōu)化時(shí)，初始化設(shè)計(jì)模塊對設(shè)計(jì)域內(nèi)的設(shè)計(jì)變量和優(yōu)化參數(shù)進(jìn)行初始化。初始化后的設(shè)計(jì)變量經(jīng)設(shè)計(jì)變量輸出模塊以電子表格的形式輸出到工作目錄。使用系統(tǒng)命令調(diào)用Python 程序執(zhí)行有限元分析。待有限元分析結(jié)束后使用有限元結(jié)果處理模塊將有限元結(jié)果組裝成靈敏度信息。使用網(wǎng)格過濾子程序?qū)`敏度信息進(jìn)行過濾。過濾后的靈敏度信息經(jīng)過設(shè)計(jì)變量更新子程序的處理得到新的設(shè)計(jì)變量。將新舊設(shè)計(jì)變量通過收斂準(zhǔn)則計(jì)算是否收斂，如果收斂則輸出優(yōu)化結(jié)果否則程序就回到設(shè)計(jì)變量輸出模塊繼續(xù)執(zhí)行。

圖2 MATLAB 程序結(jié)構(gòu)

2.2 Python 程序的實(shí)現(xiàn)

Python 程序的結(jié)構(gòu)如圖3 所示。該程序可以分為有限元前處理、有限元分析和有限元后處理三部分。當(dāng)接到MATLAB 程序發(fā)出的系統(tǒng)命令時(shí)，該程序開始執(zhí)行。在前處理階段程序的操作對象為.CAE 文件，使用MATLAB 輸出的設(shè)計(jì)變量對模型中所有單元的材料屬性進(jìn)行更新。使用有限元分析模塊將更新后的模型提交有限元分析并等待分析結(jié)束。分析結(jié)束后使用有限元后處理模塊，將有限元分析結(jié)果中的彈性應(yīng)變能以電子表格的形式輸出到根目錄。如果處于第一次迭代中Python 程序還會(huì)輸出模型中所有單元的中心坐標(biāo)。

圖3 Python 程序結(jié)構(gòu)

3 優(yōu)化算例驗(yàn)證

根據(jù)圖4（a）中模型示意圖，創(chuàng)建一個(gè)長度為60 寬為20 的二維懸臂梁結(jié)構(gòu)。梁的左側(cè)施加對稱約束而梁的右下角施加完全固定約束。在梁的右上角施加一個(gè)方向向下大小為1 的力。在材料方面，創(chuàng)建的材料屬性中將彈性屬性中的楊氏模量值設(shè)定為1 泊松比為0.33。在劃分網(wǎng)格時(shí)，設(shè)定每個(gè)網(wǎng)格為網(wǎng)格尺寸的邊長為1 形狀為四邊形的四節(jié)點(diǎn)網(wǎng)格，因此在梁中共有1200 個(gè)網(wǎng)格單元如圖4（b）中所示。完成上述步驟后，將該模型的名稱、工作目錄、模型數(shù)據(jù)庫名稱、網(wǎng)格單元數(shù)等信息輸入拓?fù)鋬?yōu)化集成程序中開始進(jìn)行拓?fù)鋬?yōu)化，進(jìn)過68 次迭代后最終的結(jié)果如圖4（c）所示。

圖4 二維算例的設(shè)計(jì)域與優(yōu)化結(jié)果

A-M 平臺(tái)與單ABAQUS 平臺(tái)運(yùn)算對比：

在經(jīng)過修改后，上述ABAQUS-MATLAB 拓?fù)鋬?yōu)化平臺(tái)（簡稱A-M 平臺(tái)）可以實(shí)現(xiàn)雙向漸進(jìn)結(jié)構(gòu)法進(jìn)行拓?fù)鋬?yōu)化。有學(xué)者曾使用Python 程序在ABAQUS 平臺(tái)中實(shí)現(xiàn)了雙向漸進(jìn)結(jié)構(gòu)優(yōu)化方法[8]。在同一硬件條件下，使用同一優(yōu)化對象與相同的優(yōu)化參數(shù)對兩個(gè)平臺(tái)的運(yùn)算時(shí)間進(jìn)行了對比，所得到優(yōu)化結(jié)果如圖5 所示，其中（a）為A-M 平臺(tái)的優(yōu)化結(jié)果（b）則是ABAQUS 平臺(tái)的優(yōu)化結(jié)果。

圖5 A-M 平臺(tái)與ABAQUS 平臺(tái)懸臂梁優(yōu)化結(jié)果

從表1 中可以直觀地看到兩個(gè)平臺(tái)在同一條件下，對同一對象優(yōu)化所耗費(fèi)的時(shí)間。在總用時(shí)方面A-M 平臺(tái)遠(yuǎn)小于ABAQUS 平臺(tái)，反映在單次迭代的平均用時(shí)上也是A-M 平臺(tái)耗時(shí)最少約為ABAQUS 平臺(tái)的三分之一左右。造成這種結(jié)果的原因是因?yàn)锳BAQUS 平臺(tái)的計(jì)算是基于Python 語言，由于該語言在處理矩陣計(jì)算時(shí)需要完全遍歷的特點(diǎn)導(dǎo)致在網(wǎng)格過濾階段ABAQUS 平臺(tái)的耗時(shí)要遠(yuǎn)遠(yuǎn)大于A-M 平臺(tái)。

表1 兩個(gè)平臺(tái)各項(xiàng)用時(shí)對比

4 輪轂的拓?fù)鋬?yōu)化

使用A-M 集成優(yōu)化方法對輪轂進(jìn)行拓?fù)鋬?yōu)化。如圖6（a）所示優(yōu)化前的輪轂結(jié)構(gòu)由白色的內(nèi)外環(huán)非設(shè)計(jì)區(qū)域和中間青色的設(shè)計(jì)區(qū)域組成，其中設(shè)計(jì)區(qū)域使用鋁合金材料非設(shè)計(jì)區(qū)域使用剛體建模，各部件的尺寸在圖中標(biāo)注。

圖6 輪轂結(jié)構(gòu)優(yōu)化示意圖

圍繞外環(huán)施加四個(gè)等大的切向力，內(nèi)環(huán)內(nèi)側(cè)則完全固定，三個(gè)部件之間使用Tie 約束。設(shè)計(jì)區(qū)域的網(wǎng)格單元共7360 個(gè)如圖6（b）所示。優(yōu)化目標(biāo)體積分?jǐn)?shù)為20%，過濾半徑為2.0。經(jīng)過159 次迭代之后得到了圖6（c）所示的優(yōu)化結(jié)果。該優(yōu)化結(jié)果滿足預(yù)設(shè)的體積約束，材料分布在集中力施加點(diǎn)與內(nèi)環(huán)周圍且以中軸線對稱。

5 結(jié)論

提出和實(shí)現(xiàn)了一種基于變密度方法的ABAQUSMATLAB 集成平臺(tái)。通過經(jīng)典算例驗(yàn)證了該集成平臺(tái)的可行性。經(jīng)過修改后該平臺(tái)可以用雙向漸進(jìn)結(jié)構(gòu)法進(jìn)行拓?fù)鋬?yōu)化。與另一種拓?fù)鋬?yōu)化平臺(tái)對比得到了本集成平臺(tái)耗時(shí)更少的結(jié)果，證明了A-M 集成方法的高效性。使用該集成方法以輪轂為對象進(jìn)行拓?fù)鋬?yōu)化，得到了結(jié)構(gòu)合理造型美觀的優(yōu)化結(jié)果。