單元微分法在拓撲優化中的運用

楊小輝

（廣東工業大學機電工程學院，廣州 510006）

0 引言

拓撲優化是一種根據給定的負載情況，約束和性能指標，優化給定區域中的材料分布的數學方法。拓撲優化的出現最初只是為了解決一般的機械設計問題。但是，隨著對拓撲優化的深入研究，其已被廣泛用于許多物理學科，包括固體力學、流體力學、熱傳導、電磁學等。此外，拓撲優化也已用于許多工程領域，例如運輸、建筑設計、復合材料等。

從Bends?e 和Kikuchi[1]提出均質化方法這一重要的拓撲優化方法到現在，拓撲優化領域已經出現了許多有效可行的方法，如帶懲罰的固體各向同性材料（SIMP）方法、演化方法、移動可變形成分（MMC）方法和漸進結構方法（MMA）[2]。

由于FEM幾乎可以用于描述所有復雜的工程問題，所以目前大多數關于拓撲優化和拓撲優化方法的研究都是基于有限元結構分析的，只有很少的研究集中在用其他有效的數值求解方法來代替FEM。因此，在這一領域中仍然需要大量的更深入的研究。

最近X.-W.Gao等[3-4]提出一種新的、強形式的數值求解方法，來求解二階偏微分方程的邊值問題。該方法分別借鑒了FEM[5-6]，有限塊體法（FBM）[7-8]和無網格法（MFM）[9-10]的部分思想，能快速有效得到穩定解。而且與FEM 相比，EDM具有兩個較為顯著的特點：（1）不需要任何數學原理或機械原理來構成方程組，因此EDM非常易于編寫；（2）無需任何積分。

本文將基于5 節點單元的EDM[11]，運用于基于SIMP 法的拓撲優化[12]，并用MATLAB 進行仿真計算，最后分析討論其拓撲結果，驗證其有效性和準確性。

1 理論基礎

1.1 單元微分法

EDM用于解決二階偏微分方程（PDEs）的邊值問題。其關鍵思想是用節點形函數的一階、二階偏導數來表示問題中的相關物理變量，并將問題計算域中的節點分為3類，不同的節點配置不同的平衡方程。具體步驟如下。

（1）用一系列5 節點等參單元劃分問題計算域，如圖1所示。

（2）推導節點形函數，并得出其一階、二階偏導數公式。

（3）用所推導的表達式表示問題中的相關物理變量，如節點坐標，節點位移等。

（4）進行節點分析，并將其分為內部節點、邊界節點、界面節點3類，如圖2所示。在不同的節點上配置不同的平衡方程（以彈性力學為例）：

①在內部節點配置應力平衡方程：

②在界面節點配置牽引力平衡方程（節點合力為0）：

tij＝ σijηj＝ 0

③在邊界節點配置牽引力平衡方程（節點受力平衡）：

tij＝ σijηj＝ ˉt

（5）統計所有節點，進行最終方程的組裝，其最終形式為： Ax ＝b ，并對其進行求解，得出結果。

圖1 2維5節點四邊形等參元素

圖2 節點的類型

1.2 拓撲優化

本文使用的拓撲優化理論是基于SIMP的結構拓撲優化理論，在SIMP 方法中，引入了相對密度為0～1 的偽可變材料。假設材料的宏觀彈性常數與其密度之間存在非線性關系，則介于0～1之間的元素會受到懲罰因子的約束。在一定數量的材料的條件下，找到具有最大剛度（結構的最小柔韌性）的結構材料的最佳分布形式。結構的順應性被視為目標函數，體積作為約束。則拓撲優化模型如下：

式中：X為單元密度；C( X )為結構順應性；F為載荷；U為節點位移；V ( X )為結構體積；V*為體積約束；K為整體剛獨矩陣。

結構的整體剛度矩陣由單元剛度矩陣組裝而成，即：

所以，對于結構順應性有：

1.3 基于EDM的拓撲優化

將EDM與基于SIMP的結構拓撲優化結合，其目標函數，即結構順應性可由以下等式表示：

式中：S為結構的受力面積；A為結構的整體系數矩陣；ae為單元的系數矩陣。

所以，拓撲優化目標函數，即結構的順應性可由如下等式表示：

2 數值案例

2.1 橋梁結構

圖3 所示為拱橋模型，長寬比設為4∶1，材料彈性模量E ＝1，泊松系數ν＝ 0.3。約束條件是結構的材料的體積不超過設計區域的體積的30%，以獲得滿足結構的最小順應性目的的最佳拓撲結構。計算域網格規模分別為80×20、240×60、320×80、400×100、480×120，并使用MATLAB 進行計算。對于EDM 的，使用上文中介紹的2 種計算方法分別計算結構的順應性，并與FEM 所得的結果進行比較。最終拓撲結構如圖4所示，具體的數值結果如表1 所示。從圖中可以看出，基于EDM 的拓撲優化得到的拓撲結構和基于FEM的拓撲優化得到的基本一致。而表中的數據顯示，在網格規模較小的時候，前者所需的迭代次數要多于后者，但是當網格規模較大時結果相反。

圖3 拱橋模型

圖4 拱橋最終拓撲結構圖

表1 拓撲優化數值結果

2.2 懸臂梁

如圖5 所示，用左側固定的懸臂梁模型簡單代表具有矩形立面的高層建筑物，該建筑物在左右兩邊高度方向上承受水平風荷載。該模型的長寬比為5∶1，材料彈性模量E ＝1，泊松系數ν＝0.3。約束條件是結構材料的體積不超過設計區域的體積的20%，以獲得滿足結構的最小順應性目的的最佳拓撲結構。計算網格分別為100×20、200×40、300×60、400×80 和500×100，由MATLAB 計算拓撲結構。在這個案例中下，僅使用EDM（3）進行計算，并將其結果與FEM 計算結果進行比較。最終拓撲結構如圖6所示，結果的具體數值如表2所示。

圖5 受對稱載荷的懸臂梁

圖6 懸臂梁最終拓撲結構圖

表2 拓撲優化數值結果

從上述結果可以看出，基于EDM的拓撲優化與基于FEM的拓撲優化的拓撲結構基本一致，但是前者所需的迭代次數要少于后者。

以上兩個拓撲優化案例證明了EDM的拓撲優化運用上的有效性和準確性，并且案例中基于EDM的拓撲優化所需的迭代次數要普遍少于基于FEM的拓撲優化所需的迭代次數。

3 結束語

本文將基于5 節點單元的EDM 運用于拓撲優化，并通過一系列數值案例進行數值仿真求解與分析。可以得到如下結論：

（1）案例中基于EDM的拓撲優化得到的拓撲結構和基于FEM的拓撲優化得到的基本一致；

（2）案例中基于EDM的拓撲優化所需的迭代次數普遍要少于基于FEM的拓撲優化；

（3）將EDM運用于拓撲優化是可行且有效的。

盡管當前的工作只是以少數的SIMP模型為例，但在將來基于EDM的拓撲優化工作會拓展到其他不同的模型和其他類型的拓撲優化方法中。