基于偶應力理論的雙向進化結構優化法*

喬赫廷，高 群，熊 展

（沈陽工業大學機械工程學院，沈陽 110870）

0 引言

在微/納尺度技術領域中，微機電系統（Micro-Electro-Mechanical System，MEMS）的結構設計逐漸成為研究的熱點問題。拓撲優化方法由于其能夠產生結構的創新式構型，可獲得比傳統的形狀和尺寸優化更好的優化結果而被廣泛應用于傳統宏觀機械結構概念設計中，成為機械領域中實現輕量化設計不可或缺的手段。但是在微觀機械結構領域中的發展卻較為遲緩。因此，為了在MEMS中實現結構的創新式構型設計，必須要研究在微觀尺度下的結構拓撲優化方法。

在宏觀機械領域，結構拓撲優化理論與方法逐步完善，拓撲優化方法被視為尋求創新結構設計的有力工具。S Mantovani等[1]針對方程式賽車轉向柱的設計問題，提出了一種以質量最小化為目標，同時遵循嚴格的結構約束的拓撲優化方法。Quhao Li等[2]對于薄壁結構承重構件的加強筋布局和截面輕量化設計問題，提出了一種并行拓撲優化方法，顯著提高了結構的性能，同時也降低了計算成本。畢政等[3]面向如何提升車輛底部防護組件的抗爆性能，以降低車身底板變形對車內乘員的威脅問題，提出了一種結合混合自動元胞法的拓撲優化設計方法，得到了防護組件中加強梁的最佳材料分布形式。同時，拓撲優化方法拓寬了其在微觀尺度上的應用，以設計滿足工作需求的微型結構。微機電系統中的結構設計也需要先進的方法來充分發揮其作用，因此進行微觀結構的設計顯得尤為重要。Abbas Homayouni-Amlashi等[4]為了設計一種從來自不同方向的面內力中獲取能量的壓電板能量采集器，采用拓撲優化方法尋找壓電板的最佳布局，以最大化電輸出并克服電荷抵消。Yi Gao等[5]考慮微觀結構的不確定性，提出了一種基于可靠性的拓撲優化框架，避免了不確定性對所得拓撲體積分數和局部特征的影響，實現精確的目標可靠性。此外，隨著材料微觀結構尺寸變小，非均質材料的有效特性取決于結構尺寸和材料微觀結構尺度的大小，這種現象稱為尺度效應。眾所周知，一些高階介質理論，例如偶應力理論可以提供合理的方法來解決這個問題，捕獲結構的力學性能和材料尺度參數，揭示優化設計中的尺度效應。Sergey P Pavlov等[6]在Timoshenko運動學假設和修正的偶應力理論的基礎上，提出了一種基于拓撲優化的納米梁微觀結構優化方法，獲得了納米梁在任意靜態和動態載荷以及不同邊界條件下的最佳拓撲，以增加其剛度。Wenzheng Su等[7]為了描述和檢驗周期性蜂窩狀固體微觀結構拓撲設計的尺寸相關性，以實現結構基頻最大化，在偶應力理論下建立優化模型，充分說明了周期性蜂窩狀固體微觀結構中優化結果的尺寸依賴性并改進了經典理論在尺度效應描述方面的弱點。

基于變密度（Solid isotropic material with penalization，SIMP）法的拓撲優化中，優化過程中會出現數值問題。如何解決數值問題已成為SIMP法優化求解中的重要組成部分。在這方面，許多研究者開展了非常有價值的工作。Quhao Li等[8]為了消除固有頻率拓撲優化過程中經常出現的數值問題，結合“有界公式”和“魯棒公式”建立了動態拓撲優化公式。Haitao Han等[9]提出了一種孔洞填充方法（HFM），用于控制最優結構中孔洞的存在，不僅限制了最優結構設計中的孔的數量，而且抑制了棋盤格和網格依賴現象。與SIMP法相比，雙向進化結構優化（Bi-directional evolutionary structural optimization，BESO）法允許添加和移除元素，收斂速度更快，效率更高，得到的優化結構更清晰，現如今已成為一種尋找最優解更穩健且更有效的算法。Mohsen Teimouri等[10]采用BESO法對新的混合實體網格結構進行拓撲優化設計，顯著提高了結構的機械性能。Yunkai Gao等[11]提出了一種改進的BESO法解決多約束拓撲優化問題，結果表明不同的約束組合產生不同的材料分布。Xiuyang Qian等[12]利用BESO法，對機翼截面進行拓撲優化設計，以去掉結構中的低效材料，證明了拓撲結構的最佳材料布局在降低應力集中和提高機翼壁板承載能力方面的表現更好。

綜上所述，本文的研究目的是提出一種在微觀尺度下應用BESO法實施結構拓撲優化設計的方法。基于偶應力理論而不是經典連續介質理論構建了矩形8節點偶應力單元。在此基礎上考察了各向同性偶應力介質剛度最大化的優化結果，并與相同工況下的經典連續介質結果進行了對比。數值算例表明，目前基于偶應力的優化模型詳細說明了微型結構設計中優化結果的尺寸依賴性，并改進了經典優化模型在尺度效應描述方面的弱點。

1 偶應力理論

本文中，ui為線彈性偶應力模型中的位移矢量，κi為獨立旋轉矢量，χij=κj,i，為曲率張量。應變張量εij依賴于位移矢量和旋轉矢量。應力張量σij和偶應力張量mij通過本構方程與εij和χij存在聯系。本節闡述了在偶應力理論下的基本方程，其中包括幾何方程、平衡方程和本構方程，以及廣義有限元列式。與經典彈性理論相比，平面偶應力理論除了原有的σxx、σyy、σzz、τxy、τyx應力分量外，新增了mxz、myz兩個偶應力分量，采用符號表示如圖1所示。

圖1 二維偶應力模型

1.1 基本方程

幾何方程：

式中：εij為應變張量；uj,i為位移向量；ej i為置換符；κk為獨立的旋轉向量；χij為曲率張量；κj,i為旋轉向量。

平衡方程：

式中：σij,i為應力張量；ψj為體積力；mij,i為偶應力張量；ei jk為置換張量；σik為獨立的應力張量；φj為體力偶。

對于平面偶應力問題，剪應力是對稱的，分別用τs和τa表示，即（τxy≠τyx）：

對稱部分τs與剪切應變γx y相關：

式中：μ為剪切模量，μ=E/[ 2( 1+δ)]，E為楊氏模量，δ為泊松比。

反對稱部分τa與剛體局部轉角ω相關，并且它不是一個獨立的量，具體可以表示成如下形式：

式中：u和v分別為面內的位移和轉角。

與之對應的偶曲率為：

考慮各向同性材料，偶應力下的本構關系表示為：

式中：σi j、εk l分別為應力和應變；Dijkl為彈性本構矩陣。

邊界條件：

對于平面應變條件下的二維連續體，本構方程可以用矩陣形式改寫為：

式中：λ=D·δ/( 1+δ)( 1-2δ)，為拉梅常數；μg為類似于剪切模量的常數；lg為材料的特征長度。

1.2 有限元列式

本文研究了矩形8節點偶應力單元，其每個節點上分別有3個自由度，如圖2所示，對于偶應力模型，根據勢能原理，則單元總勢能泛函為：

圖2 平面矩形單元

位移場函數表述如下：

式中：ui為單元節點線位移；vi為節點角位移；ξ和η為局部坐標系中的自然坐標；Ni為形函數。

Ni在8節點單元中分別表示成：

將式（14）代入（13）取變分，并使其為0，得到的有限元列式可以簡寫成如下：

式中：Ke為單元剛度陣；Ue為節點位移列向量；Fe為節點載荷向量。

依據經典單元矩陣的表現形式，得出單元剛度矩陣，即：

2 雙向進化結構優化算法

2.1 優化模型

剛度是影響結構能否穩定工作的關鍵參數，通常情況下，結構優化的最終目的是實現剛度最大化，也就是說，在給定的材料約束條件下，需要考慮的問題是如何使結構柔度最小，以提高抵抗彈性變形的能力。其優化問題可以定義為如下形式：

式中：c(x)為平均柔順性；FT為力矢量的轉置；U為位移矢量；K為總剛度矩陣；V p為單個元素的體積；V*為規定的總體積；N為網格中的元素總數；xp為設計變量，為了防止剛度矩陣產生奇異問題，使用較小的xmin值來表示空洞元素。

2.2 求解策略

在優化過程中，BESO法的工作原理與SIMP法有所不同，它的材料去除方案允許添加或移除材料，這主要取決于元素靈敏度數值的大小，當靈敏度為正時，添加材料；當靈敏度為負時，移除材料，同時要通過靈敏度過濾方案來避免數值不穩定現象的出現。因此，使用下式進行濾波：

式中：Q為設計域中單元的節點總數；b ho為元素中節點h和o中心的距離；ιo為靈敏度數；?(b ho)為權重因子。

式中：R為濾波器半徑。

式（20）表示要將在該濾波半徑范圍內的元素靈敏度值作為新的靈敏度數值，以獲得更好的優化結果。

上述靈敏度濾波過程稱為空間平滑，但是存在收斂問題，為了保證優化過程的絕對收斂，引入時間平滑的概念，其主要思想是通過取靈敏度的平均值來防止優化過程中結構性能發生振蕩。它可以表示為：

其中t和t-1分別表示當前迭代次數和上一次迭代次數，在每次優化過程中，=φp將會被用到下一次迭代。

在BESO法中，為了確定下一次迭代的目標體積，采用如下方程計算：

其中M為進化體積比，當V＞V*時，則根據上式計算并更新得到下一次迭代的目標體積，當V＜V*時，則結構的體積將在之后迭代中保持不變。

圖3所示為BESO法的循環流程，包括確定目標函數和目標體積約束的詳細步驟。

3 數值算例

在這一節中，分別給出了幾個例子來驗證方法的可行性。所有算例均采用的是基于偶應力下的8節點矩形單元求解，并且假設所有固體材料的楊氏模量為1，泊松比為0.3，其進化體積比設為0.04。

優化過程的收斂準則表示如下：

其中ρc(x)設為0.5%，?V設為0.4%，Z為整數，取值為5。此外，如果迭代次數達到200次，則優化將被終止。

3.1 算例一

本例模型為長L=90 mm，寬H=30 mm的Michell桁架，底部兩端固定約束，并在中間施加豎直向下的載荷，如圖4所示。將該模型離散化為90×30個8節點偶應力單元，體積約束為0.5。首先從圖5所示的變形示意圖上可知，結構自身各部分的變形程度分配較為平均。符合Michell桁架的工作機理以及微型結構的工作機制。其次，研究了不同尺度下模型的優化結果，將偶應力材料的特征長度lg與梁高度H的比值稱為尺度參數，記為ζ。如表1所示，依次考察了ζ為0.1，0.2，0.4，0.6，0.9以及1時的優化結果。隨著ζ的增加，優化結果會有明顯的差異。當ζ在0.1～0.4之間，拓撲結構幾乎相同，具體地說，結構內部的孔洞變得越來越圓滑。當ζ=0.6時拓撲結構發生了顯著的改變，孔洞由2個變為1個。當ζ=0.9時，孔洞徹底消失，最終得到穩定的結構。應該指出，當ζ增加到與域的維數相當時，拓撲結構會發生很大的變化。繪制桁架在ζ=0.6時的平均柔度、體積分數和拓撲優化的演變歷史，以研究BESO法的質量，如圖6所示。通過比較，隨著體積分數的減小，平均柔度先增大，然后在達到目標體積后收斂到一個幾乎恒定的值。經過48次迭代，最終收斂到穩定的拓撲結構。

圖4 底部兩端固支Michell桁架結構

圖5 Michell桁架結構優化結果變形示意圖

表1 Michell桁架的最佳拓撲結構

圖6 Michell桁架的平均柔度和體積分數的演變歷史

3.2 算例二

在這個例子中，采用長為L=80 mm、高為H=30 mm的懸臂梁作為研究對象，左端固定約束，上方承受均布載荷，如圖7所示。宏觀設計域離散化為160×60個8節點偶應力單元，材料的體積約束為0.4。最優構型的變形如圖8所示，最大位移發生在結構上端，由此可以得出在微觀尺度下懸臂梁的結構設計符合最大剛度要求。同時為了比較本文采用的BESO法和SIMP法的優化設計，在不同尺度下，應用BESO法求解MBB梁問題。從表2中可以看出，隨著尺度參數的增大，MBB梁的柔度逐漸變小，迭代時間和次數減少，說明尺度的變化會對優化過程造成影響，并且宏觀尺度和微觀尺度下結構的性能有明顯的差異。

圖7 承受均布線載荷的懸臂梁

圖8 均布線載荷懸臂梁優化結果變形示意圖

表2 不同尺度下承受均布線載荷的懸臂梁的結果比較

表3所示為ζ取0.1、0.6和1尺度下的迭代歷史圖，其中2列和4列為算例模型同一尺度下不同算法的迭代過程，可以看出SIMP法在迭代過程中會出現灰度單元，這一點充分說明了本文采用的BESO法比SIMP法在數值求解中穩定性更好。并且該表的3～5列反映了材料微觀尺度下結構的特性，因此偶應力介質拓撲優化具有明顯的尺度效應。隨著精密加工工藝的出現，往往需要設計微型結構的尺寸或與之配合的構件，而基于偶應力理論的拓撲優化方法在這一應用領域表現得更出色。

表3 迭代歷史圖

4 結束語

本文通過將偶應力理論與BESO法結合，針對微型結構實行最小柔度拓撲優化設計，事實證明了偶應力介質下結構存在顯著的尺度效應，解決了經典介質理論無法得到微觀結構的信息，提供了一種微觀結構設計的新方法。最終的拓撲結果取決于尺度參數的大小，且不同尺度下結構的表現形式是不同的，當尺度參數相差較大時，宏觀尺度下產生的是桁架結構，而微觀尺度下會出現不規則的橢圓形孔。并且偶應力介質具有捕獲微觀結構信息的能力，因此尺度效應明顯，這一點是經典連續介質無法做到的。此外，與SIMP法相比，本文所用的BESO法不會產生灰度單元，數值穩定性更好。