基于變密度法的ABAQUS-MATLAB 集成拓撲優(yōu)化方法

靳立涵JIN Li-han

（重慶交通大學機電與車輛工程學院，重慶 400074）

0 引言

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

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

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

1 SIMP 拓撲優(yōu)化方法理論

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

SIMP 法拓撲優(yōu)化模型：

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

在本文中討論的SIMP 表達式為：

其中xi是編號為i 的單元的相對密度。

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

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

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

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

圖1 ABAQUS-MATLAB 集成方案

根據(jù)結(jié)構(gòu)拓撲優(yōu)化設(shè)計要求，在ABAQUS/CAE GUI中建立有限元模型，模型包括網(wǎng)格單元、材料屬性、邊界條件和載荷應(yīng)用，其中材料屬性可以被更新。MATLAB 代碼完成求解初始化，包括有限元模型的單元總數(shù)和迭代的相關(guān)參數(shù)。隨后，在循環(huán)迭代開始時進行拓撲優(yōu)化，Python 程序讀取MATLAB 代碼生成的各單元設(shè)計變量，根據(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ù)收斂準則判斷是否收斂，如果不收斂則進行下一次循環(huán)迭代，直到收斂。*.XLS 文件是MATLAB 和ABAQUS之間傳輸數(shù)據(jù)的介質(zhì)，該文件是通過代碼調(diào)用Microsoft Excel 生成的。

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

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

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

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

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

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

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

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

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

A-M 平臺與單ABAQUS 平臺運算對比：

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

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

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

表1 兩個平臺各項用時對比

4 輪轂的拓撲優(yōu)化

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

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

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

5 結(jié)論

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