限界法克服連續體結構拓撲優化載荷病態問題

俞燎宏，榮見華

（1.宜春學院物理科學與工程技術學院，江西 宜春 336000；2.長沙理工大學汽車與機械工程學院，湖南 長沙410076；3.工程車輛輕量化與可靠性技術湖南省高校重點實驗室，湖南 長沙 410076）

1 引言

自文獻[1]提出均勻化方法以后，連續體結構拓撲優化方法迅速發展，出現了實體各向同性材料懲罰法（SIMP）[2]、漸進結構優化法（ESO/BESO）[3-4]、水平集法（Level Set Method）[5-7]、獨立連續映射法（ICM）[8]和可移動變形組件法（MMC）[9]等。其中，SIMP方法目前應用最為廣泛。在實際工程中，很多結構通常受到多種工況載荷的作用。比如主要承載機械的結構上可能需要安裝傳感器或操作平臺等，故該結構除了強載荷作用外，還在偏離強載荷作用部位的某處存在小載荷作用。要想在優化過程中獲得性能良好的拓撲結構，必須同時考慮所有的載荷工況。

由于不同工況載荷的大小、作用方向和作用位置不同，對結構的作用效果也不盡相同。那么，在基于傳統優化方法的多工況載荷作用下結構拓撲優化過程中，可能出現沒有材料支承小載荷的優化結構，這就是“載荷病態”問題。為了解決載荷病態問題，文獻[10]提出了分層優化方法；文獻[11]提出了ICM 應力全局化方法；文獻[12]提出了分數模常法；文獻[13]對上述方法進行了綜合分析和評判，并提出了敏度分層過濾解決載荷病態的策略。該策略能夠得到較為清晰的結構，可以克服荷載病態問題，但是優化結果依賴于兩個關鍵系數。

多工況載荷作用下的結構柔順度拓撲優化問題本質上是一個多目標優化問題，關于多目標優化問題的研究綜述可見文獻[14]。針對以多工況載荷下結構柔順度最小為目標函數、結構體積為約束條件的優化問題，文獻[15]采用限界法提出了一種新的求解思路，能以較少的工作量尋求更優的結構拓撲。限界法[16-17]非常適合于不可微的最小最大優化問題。文獻[18]將限界法應用到振動連續體結構多特征頻率和頻率間隙最大化的拓撲優化設計。文獻[19]利用限界法開展了制造約束下復合材料框架最大基頻離散選材與結構拓撲優化。

經過對現有載荷病態處理方法分析和研究的基礎上，提出采用限界法來克服載荷病態問題。首先，給出以最大的結構柔順度最小化為目標，以體積為約束的拓撲優化模型。其次，引入限界變量將多個柔順度目標函數轉化為約束條件，結合變體積約束限技術，建立新的等效近似優化模型；針對同一工況內含有病態載荷的問題，在新建的優化模型中增加柔順度小量變化約束。然后，對目標函數和約束函數進行靈敏度計算，采用移動漸近線方法（MMA）優化求解。最后，給出了二維和三維算例，驗證了所提方法的可行性和有效性。

2 多工況載荷下的結構拓撲優化模型

2.1 材料插值與Heaviside過濾

通常，基于SIMP材料插值模型和Heaviside過濾技術，將結構的設計域離散成N個單元，每個單元賦予一個設計變量和一個相對密度變量。其中，第e號單元的相對密度變量x?e(e=1，2，…，N)與它周圍單元設計變量xi(i=1，2，…，N)之間存在如下關系：

當與Heaviside過濾技術一起使用時，常數加權可以更有效地減少了灰度，而線性加權則傾向于生成更平滑的拓撲構型。在文中采用常數加權形式。

為了獲得0/1分布清晰的結構拓撲，采用文獻[20]提出的光滑Heaviside函數，其表達式為：

式中：xˉe—單元物理變量；β—曲率參數。

通過式（1）和式（5）形成設計變量場x、過度變量場x?和物理變量場xˉ。結構有限元分析時，用物理變量xˉe替代相對密度變量x?e。

結構的單元體積和剛度矩陣通過公式（5）來獲得：

式中：Ve和Κe—第e號單元的體積和剛度矩陣；—第e號單元的固有剛度矩陣單元體積的懲罰函數文中選取剛度矩陣的懲罰函數其中，E0—實體材料的彈性模量；Emin—空洞材料的彈性模量，為了避免剛度矩陣奇異，這里選取Emin=10-9E0；p—懲罰因子，文中選取p=3。

2.2 多工況載荷下結構柔順度拓撲優化模型

在實際工程中，結構通常在復雜的工況載荷下工作，不同工況載荷對結構的要求也有差別。在多工況載荷作用下，以結構體積為約束條件，以最大的結構柔順度最小化為目標的拓撲優化模型可以表示為：

式中：L—工況載荷總數；Fl—第l號工況載荷向量；Cl和Ul—第l號工況載荷作用下結構的柔順度和位移向量；V0—整個結構的初始體積；f*—預先給定的體積比。

3 利用限界法克服載荷病態問題

3.1 基于限界法的拓撲優化模型

參考文獻[14-16]的限界法，引入一個限界變量z(z＞0)將優化模型式（6）中的多個目標函數轉化為約束條件，即將多目標優化問題轉換成單目標優化問題。為了協調Heaviside變換引起的結構體積沖突，同時確保每一迭代步的結構特性函數穩健變化，引入變體積限約束[21]。其轉換后的近似等效優化模型為：

式中：Cmax—所有工況載荷作用下的最大結構柔順度；VU(k)—設定的當前步結構體積限，其更新表達式為：

式中：V*—目標體積；V(k-1)—上一步得到的結構總體積；γ—體積限的經驗參數，其取值范圍為[0.005，0.025]；Q—Heaviside映射曲率參數β變化后的前三個迭代步的迭代步編號的集合。

3.2 柔順度小量變化約束策略

文獻[11]將載荷病態分為三類：第一類是多工況間有載荷病態，但工況內無載荷病態；第二類是僅在工況內部有載荷病態；第三類是多工況間有載荷病態，同時某工況內也有載荷病態。

為了便于描述，將引起載荷病態問題的工況載荷稱為“病態載荷”。通過數值算例表明，優化模型式（7）可以解決第一類載荷病態問題，但是不能解決同一工況內有載荷病態的問題，即不能解決第二類和第三類載荷病態問題。

為了解決同一工況內有載荷病態的問題，借鑒文獻[21]提出的方法：將柔順度小量變化要求作為約束條件，建立新的近似優化模型。文中將其稱為柔順度小量變化約束策略，其實現過程如下：

（1）將同一工況內的病態載荷單獨取出來作為新的虛載荷，并將其獨立加載在結構上。

（2）將虛載荷作用下的結構柔順度按照表達式（10）形成新的約束條件。

（3）將柔順度小量變化約束表達式（10）加入優化模型式（7），形成新的優化模型式（11）。

據董事長王茂莆介紹，公司創業之初，國內噴灌機產業面臨外資機構大規模并購的尷尬局面，但為了實現行業民族品牌的發展，華雨公司迎難而上，憑借“踏實、拼搏、責任”的企業精神和“誠信、共贏、開創”的經營理念，通過不懈努力，終于生產出具有我國自主知識產權的大型噴灌機設備，為贏得國產噴灌機市場定價權起到了積極的促進作用。王茂莆說，我國灌溉產業是個存量很大的市場，隨著十八屆三中全會的召開，“調整產業結構，轉變發展方式”已經被納入國家體制改革的宏觀計劃，未來新興灌溉產業將會有更大的發展空間，因此把握產業命脈的關鍵是打好基本功，夯實自身優勢。

如果工況內有多個病態載荷，則添加對應的多個柔順度小量變化約束表達式。

4 靈敏度計算及優化求解

4.1 目標函數和約束函數的靈敏度計算

在優化模型式（11）中，目標函數fo(x)對物理變量的導數為：

結構柔順度Cl(x)對物理變量xˉe的導數為：

式中：Ul，e—僅與結構第e號單元自由度相關的Ul的分量構成的矢量。

柔順度小量變化約束函數?s(x)對物理變量的導數為：式中：—第s個病態載荷作用下的位移向量僅與結構第i號單元自由度相關的的分量構成的矢量。

結構體積V(x)對物理變量的導數為：

由前文可知，結構的任意性能函數f(x)對設計變量xi的導數為：

4.2 MMA優化求解

優化模型式（7）和式（11）采用移動漸近線方法（MMA）求解，具體求解過程參見文獻[22]。對于載荷病態優化問題，由于不同工況載荷的作用效果相差較大，優化過程中結構柔順度會發生劇烈振蕩。為了抑制振蕩行為，同時協調Heaviside變換引起的急劇變化，參考文獻[23]，對MMA算法進行適當的修改。（1）修改前兩次迭代上、下移動漸近線表達式（19）中的移動漸近系數s0，修改后的s0如式（20）；（2）將式（21）引入到MMA算法中，通過和xml來控制設計變量xi的變化。

式中：曲率參數β的初始值β(0)=1，每迭代50 步按β(m)=1.5*β(m-1)進行更新。

5 優化算例

5.1 算例1

文獻[11]采用的基結構，尺寸為2000mm×400mm，厚度為9mm，彈性模量E=68.89GPa，泊松比v=0.3，左右兩端由鉸鏈支承，如圖1所示。受到兩組工況載荷作用，工況1：兩個集中載荷=200N分別作用于上邊界距離梁端100mm處的A點和C點；工況2：一個集中載荷F2=20000N 作用于梁上邊界中點B。結構劃分為100×20個矩形網格單元，體積比f*=0.5。

圖1 基結構Fig.1 Basic Structure

本算例工況間有載荷病態但工況內無載荷病態，屬于第一類載荷病態問題。本算例過濾半徑rmin=2.0Δmi（nΔmin是網格單元邊長的最小尺寸），收斂經驗參數ε=0.002，采用所提方法獲得的優化歷程和優化結果，如圖2所示。文獻[11]方法獲得的優化結果，如圖3所示。

圖2 優化歷程和優化結果Fig.2 Optimization Process and Results

圖3 文獻[11]方法獲得的優化結果Fig.3 Optimization Results Obtained by Reference[11]

其中圖3（a）為反演前的拓撲值分布圖（紅色單元的拓撲值為1，其它顏色單元拓撲值介于0與1之間），圖3（b）為反演后的最優拓撲圖形。從圖2和圖3可見，優化結果整體構型基本相似，但支撐病態載荷傳力路徑的構件有所不同。采用所提方法不需要進行反演就能夠獲得0/1分布清晰的拓撲結構，如圖2（d）所示。離散度指標md=2.63%，病態載荷的傳力路徑清晰。所提方法和文獻[11]方法優化得到最大工況（工況2）的最小柔順度（最小應變能）分別為20.04N·m和22.37N·m。因此，所提方法可獲得柔順度更小的優化結果，更符合柔順度最小化的優化目標。該算例表明，對于工況間有載荷病態同時工況內無載荷病態的拓撲優化問題，采用所提方法是可行且有效的。

5.2 算例2

文獻[13]采用的初始結構，為300mm×200mm的長方形平面區域，左右兩邊受固定支撐，如圖4所示。單工況內受兩個載荷作用，集中載荷F1=1000N作用于上邊界中點A，方向向下；集中載荷F2=1N作用于下邊界中點B，方向向上。結構劃分為60×40個矩形網格單元，體積比f*=0.5。本算例僅在工況內有載荷病態，屬于第二類載荷病態問題。對于這一類問題，將病態載荷F2作為虛載荷，并獨立加載到結構中的B點。本算例中柔順度小量變化經驗參數?*=0.002，采用所提方法獲得的優化結果，如圖5所示。如果初始結構僅受到集中載荷F1的作用，該問題不屬于載荷病態問題，優化獲得的結果，如圖6所示。

圖4 初始結構Fig.4 Initial Structure

圖5 優化結果Fig.5 Optimization Result

圖6 僅載荷F1作用下的優化結果Fig.6 Optimization Result Only Under F1

由圖5和圖6可知，優化結果中集中載荷F1的傳力路徑基本相似，圖5僅比圖6多了病態載荷F2的傳力路徑。可見，對于僅在工況內有載荷病態的問題，按提出的柔順度小量變化約束策略不會改變主要載荷的傳力路徑，而且得到了清晰的病態載荷傳力路徑，表明該方法是可行的。

采用文獻[13]方法，不同參數a和b對應的優化結果，如圖7所示。對應的最小柔順度分別為4975.556N·m、4993.970N·m、4940.627N·m和4966.611N·m。由圖5和圖7可知，載荷F1的傳力路徑基本相似，載荷F2的傳力路徑略有差異；圖5結構0/1分布更加清晰，離散度指標md=0.15%。采用所提方法計算得到載荷F1作用下的最小柔順度為4739.244N·m，比文獻[13]方法計算得到的柔順度值更小。該算例表明，對于同一工況內有載荷病態的拓撲優化問題，采用所提方法是可行且有效的。

圖7 文獻[13]方法優化結果Fig.7 Optimization Results Obtained by Reference[13]

5.3 算例3

初始結構為600mm×400mm×40mm 的長方體三維區域，左右兩端面受固定支撐，如圖8所示。

圖8 初始結構Fig.8 Initial Structure

受到2組工況載荷作用，工況1：集中載荷F1=10000N豎直向下作用于結構頂面中心A點，F′1=10N豎直向上作用于結構底面中心B點；工況2：集中載荷F2′=F2′′=10N豎直向下分別作用于結構頂面C點和D點。劃分60×40×4個矩形網格單元，體積比f*=0.4。

本算例工況間有載荷病態同時某工況內也有載荷病態，屬于第三類載荷病態問題。對于這一類問題，將工況1中的病態載荷F1′作為虛載荷，并獨立加載到結構中的B點，柔順度小量變化經驗參數?*=0.005。采用所提方法獲得的優化結果，如圖9 所示。其中，X-Z平面視圖，如圖9（a）所示；30°三維視圖，如圖9（b）所示。由圖9 可知，優化結果單元0/1 分布清晰，離散度指標md=3.7%，病態載荷F′1、F2′和F2′′的傳力路徑被保留并清晰展現。可見，按照提出的方法能很好地克服工況間有載荷病態同時某工況內也有載荷病態的拓撲優化問題，即能解決第三類載荷病態問題。該算例也表明，所提方法在三維空間結構優化中也是可行且有效的。

圖9 三維結構優化結果Fig.9 Optimization Results of 3D Structure

6 結論

針對連續體結構拓撲優化中的載荷病態問題，通過限界法結合柔順度小量變化約束策略和變體積約束限技術，基于移動漸近線方法（MMA），提出了一種克服載荷病態問題的新的優化求解方法。數值算例結果表明：

（1）所提方法可以清晰的展現病態載荷的傳力路徑，能夠有效克服三種不同類型的載荷病態問題；（2）在相同優化條件下，與文獻[11-13]方法相比，所提方法可高效地獲得0/1分布更清晰、性能更優的結構拓撲；（3）該方法不僅適用于二維平面結構，同時也適用于三維空間結構，為載荷病態問題的優化求解提供了新的思路。