比較基于化整交融應力拓撲優化諸解法的效果1)

彭細榮 隋允康 葉紅玲 鐵 軍

* (湖南城市學院土木工程學院,湖南益陽 413000)

? (北京工業大學材料與制造學部,北京 100124)

** (天津財經大學理工學院,天津 300222)

引言

連續體結構應力拓撲優化由于面對局部性能約束,較之全局性能約束問題,具有更多求解困難.對此,諸多學者例如Le 等[1]總結了主要困難,一是應力奇異性,二是應力局部特性導致敏度計算量大,約束數目多.

本文在贊同這一看法的同時,認為首先應當區分困難的主次,且不可漏掉相關的困難.一是主要困難:約束數目多,應力的局部性特點,要求每個單元的約束都不可違背,于是導致約束總數劇增,以靜力問題為例,約束數目等于單元總數與載荷工況數目的乘積,從而敏度計算量難以接受,造成建立優化模型的極大困難;二是現有的有些解法存在個別單元應力超限的現象,這緣于應力較強非線性特征,在形狀與拓撲突變區域出現應力集中,也必須解決;三是應力奇異性的困難,根本原因在于,在迭代中傾向于被刪除的子區域或單元集,其中的拓撲變量和應力分別趨于0 的速度決定了是否出現應力奇異性現象.旨在解決應力奇異性問題,Cheng 和Guo[2]提出了ε-松馳法,Bruggi[3]提出了qp-松馳法等.

為了克服約束數目多的困難,以往的研究提出了分部、化整和集成3 種解法.

分部方法[4-8]即滿應力方法,或應力的零階近似處理.隋允康等[4-5]提出了基于ICM 方法的應力零階近似處理的方法.程耿東等[6-7]提出了修改的滿應力方法求解二維連續結構的應力拓撲優化問題.Guan 等[8]用雙向進化結構優化方法(BESO)處理應力約束問題,依據單元應力值進行刪除或加入.應力零階近似的方法避免了求應力敏度,約束則轉化為設計變量動態下限.但Zhou 和Sigmund[9]指出,滿應力方法僅對設計區域無幾何限制的拓撲優化問題有效,對有幾何限制的問題,無法得到合理拓撲.

化整方法[10-12]將應力約束集合成為全局化的約束.隋允康等[10-11]將局部應力約束轉化為結構應變能約束,或按結構畸變比能進行轉化處理[12].但化整解法存在個別單元應力超限的現象,一些迭代的算法被提出[10-12],但是解決方案屬于準則方法,不能令人滿意.

集成方法[1,13-25]采用集成函數如模函數(或稱為P 范數)或 Kreisselmeier-Steinhauser (K-S)函數等將諸多的局部應力約束集成為一個或幾個全局化的函數.K-S 函數最初應用于控制系統的設計[13],后來被應用于工程結構優化中,將多目標或大量約束集成為單一目標或約束.Yang 和Chen[14]最先提出了用K-S 函數將眾多應力集成為一個應力函數,進行結構拓撲優化設計.隋允康等[15]為解決多工況下應力約束的連續體結構拓撲優化問題,引入了K-S 函數對應力約束進行集成化處理,并進而提出拋物型凝聚函數對應力約束進行集成化處理[16].Le 等[1]及Luo 等[17]基于SIMP 法用集成處理方式求解應力拓撲優化問題.Picelli 等[18]采用集成函數及水平集方法解決應力拓撲優化問題.Xia 等[19]提出了一種基于集成函數的BESO解法.王選等[20]在BESO 方法中采用K-S 函數集成應力,用拉格朗日乘子法求解,高云凱等[21]使用K-S 函數集成結構畸變比能.Wares和Hideyuki[22]對最大應力極小化的結構拓撲優化問題,應用K-S 函數對應力集成,并給出了二階近似,應用牛頓梯度法進行求解.但集成的方法亦難以保證應力處處滿足,增大集成函數的參數可以更好逼近應力最大值,但又易引起優化求解中的數值困難.于是一些解決此問題的算法被提出.Andrew 等[23]以翼盒結構質量極小化應力約束優化問題為例,比較研究了K-S 函數及模函數等不同集成函數及其參數對優化結果的影響,算例結果表明采用誘導指數集成函數并分為多個子域凝集的效果最好.Zhang等[24]對機翼結構使用代理模型進行全局優化,對比了K-S 函數對應力約束的不同凝集處理策略,比較了固定K-S 函數參數及自適應參數、全部結構應力約束集成及劃分多個子域應力約束凝集,結果表明采用較大的固定K-S 函數參數及多個子域集成的策略效果最好.Graeme[25]對應力拓撲優化問題提出了基于K-S 函數的自適應方法,可自動調整K-S 函數的參數及應力約束集成的子域劃分.但是以上所提出的諸多解決方案均不能令人非常滿意.

克服約束數目多的困難其他處理策略還有Guo等[26]提出的主動集策略,Zhang 等[27-28]提出了結構邊界曲率信息及應力梯度信息加權的應力集成策略,Xia 等[29]提出懲罰加權的應力集成策略,Fernando等[30]基于增加拉格朗日法的應力懲罰處理策略等.

為了徹底解決應力超限現象和同時提高求解效率.對應力、局部穩定及疲勞壽命等局部性能結構優化問題,隋允康等[31-33]提出了互逆規劃理論,闡述了單目標多約束模型(s 方模型)及多目標單約束模型(m 方模型)之間存在互逆關系,并應用互逆規劃理論解決結構拓撲優化建模中的合理性問題,進而基于互逆規劃的理論提出了交融優化的3 種解法——分部-集成、化整-集成、集成-集成解法.

本文的工作針對應力拓撲優化的交融優化解法,由于分部-集成解法的前半部分屬于準則方法,不必再討論,集成-集成解法的研究將另文給出,本文的研究聚焦在化整交融(化整-集成)方法的研究上,并且同化整方法進行比較.

首先,敘述集成方法即s 方模型的化整解法,然后,對m 方模型集成解法,提出了乘子法、一階近似及二階近似等不同解法,最后,敘述了基于化整處理的的s 方模型和K-S 函數集成處理的m 方模型交替交融的迭代解法,亦即化整交融解法.用數值算例對比了不同m 方模型集成解法的效率,選擇了求解效率最佳的二階近似法對m 方模型進行求解,在此基礎上,對比了化整解法與化整交融解法的求解效果.

1 s 方模型化整解法[34-35]

應力約束下最小化結構重量(或體積)化整解法的優化模型,屬于s 方模型,即

式中,t=[t1,t2,···,tN] 為拓撲設計向量,W(t) 為結構重量,wi(ti) 為單元重量,ψil(t) 為i 單元在l 號載荷工況下的局部性能如單元應力,N 為載荷工況總數,L為載荷工況總數.t 為防止結構分析時產生奇異而設置的拓撲變量下限值.

當應力約束化整為全局化應力約束時,有

其中,eil為i 單元在l 載荷工況下的應變能,Vs為結構體積,E,μ 和 σ 分別為材料彈性模量、泊松比和許用應力.

式(2)的左端為結構應變能,其中單元應變能為

其中 σil,εil與 Pil為i 單元在l 載荷工況下的應力、應變與節點力向量,Ki與 Ωi分別為i 單元的剛度矩陣及子域.

將拓撲變量 ti代入單元剛度矩陣中,得到

再代入式(3)得

在第 ν 次迭代中,由靜定化假設,該輪中的內力不變,于是得應變能函數為

將式(6)代入式(7)得

將式(8)代入式(2)得

其中Ve為單元體積,N(ν)為第 ν 次迭代的主動拓撲變量總數,取 N(ν-1)近似替代,為結構許用應變能.

由此,s 方優化模型式(1)中的局部應力約束轉化為L 個全局化的應變能約束,得到應力全局化的s 方優化模型

2 m 方集成模型的不同逼近解法

化整交融解法的第二步是多目標的m 方模型求解算法.m 方模型是以結構重量(或體積)為約束,單元最大局部性能如應力極小化為目標的優化模型,即

式中符號意義與式(1)相同.

對式(12)中的目標函數,采用K-S 函數進行集成化處理,即為m 方模型的集成解法[33],得到集成轉化后的優化模型

以下給出m 方集成化處理模型式(13)的3 種解法,為后文敘述的方便,每種解法在括號內給出了其簡稱,分別如下:

(1) Lagrange 乘子解法(乘子法);

(2) Taylor 一階近似逼近的解法,采用移動漸近線法[36]求解(MMA 法);

(3) Taylor 二階近似逼近的解法,采用序列二次規劃法求解(SQP 法).

2.1 Lagrange 乘子解法

對于應力約束,式(13)中

為了求出規劃(13)的解,先列出它的Lagrange 函數

由式(15)得鞍點條件

從式(16)得

從而有

定義

按式(20)代入式(17)得

由式(22)得

由式(23)代入式(20)得

其中 βi見式(21).注意,式(24)的 ti取因而式(24)左端 ti取,于是式(24)就構成了一個迭代式.迭代求解中當 t(v+1)與 t(v)充分接近后收斂,每次小循環迭代中都用 0 ≤ti≤1 修改迭代解.

2.2 一階近似的MMA 解法

約束函數為線性函數,其函數及一階導數為

很多學者都在研究高管薪酬激勵與技術創新之間的關系，其理論各不相同，但核心基本相似：為獲超額利潤與競爭優勢，企業傾向于高管做技術創新，而高管進行技術創新的程度與企業的薪酬激勵相關，因此，提高高管薪酬是企業技術創新的基本手段。技術創新是企業發展的根本，技術創新需要資金等條件的保障。高管是企業發展的關鍵因素，讓高管利益與企業利益進行掛鉤，是企業進行技術創新突破的關鍵。通過薪酬激勵，促使高管引領企業技術創新，是保障企業科學發展的重要手段[4]。

移動漸近線法(MMA)是基于一階導數求解的算法[36],有了目標函數及約束函數的原函數及其一階導數如式(26)～ 式(29)所示,利用MMA 求解器可直接對模型式(25)進行求解.

2.3 二階近似的序列二次規劃解法

優化模型仍為式(13),轉換到x 空間里的式(25).對式(25)中目標函數進行二階Taylor 近似,記原函數和一階偏導數分別為式(26) 及式(27) 所示.后面將采用可分離規劃解法迭代逼近式(25)的準確解,為此將式(27)所含求和項中的xi暫時取為于是原函數的Hessian 矩陣中只有對角線上的二階導數為非零值,其為

其中xi(k)的初始值取為下標(k)為小循環標志,代表一次分析之后,通過二階近似迭代逼近求解;上標(v)為大循環標志,代表結構分析的第v 次,亦即一次大迭代尋優過程.

對原函數式(26)的二階近似式略去對優化沒有影響的常數項,為

為了求出規劃式(33) 的解,列出其Lagrange函數

由式(34)得鞍點條件

由式(35)可求得主動變量為

解得

將式(40)代入式(38)得主動變量為

把式(41)由主動集延伸到非主動集,則有下式

式(42)的計算伴隨著主動集的更新,直到主動集不變為止.

3 化整交融解法的計算流程

化整交融模型的解法是通過交替求解s 方模型及m 方模型,以達到更好收斂速度的目標.s 方為化整處理法建立的優化模型,m 方為集成處理法建立的優化模型.一次結構分析的大循環內,先進行s 方化整模型式(11)的求解,再進行m 方模型公式(13)的求解,流程如圖1 所示.

圖1 化整與集成解法流程圖Fig.1 The flow chart of globalization-aggregation method

4 數值算例

4.1 算例1:m 方模型K-S 函數集成應力拓撲優化

如圖2 所示,基結構平面體的尺寸為50 mm×100 mm×10 mm,采用 50×100 網格.材料彈性模量為E=210 GPa,泊松比 ν=0.3 .左邊界固定,右邊界中點作用集中載荷 F=5000 N 作用于上邊界中點,為避免應力集中,分散作用在鄰近的3 個節點上.體積比約束為0.2%,收斂精度取0.01.

圖2 兩工況算例分析及優化的模型Fig.2 Model for analysis and optimization of the example with two load cases

單獨的m 方模型存在優化模型提法的不合理性[34-35,37].此算例的設計僅為了比較求解m 方集成模型式(25)的3 種不同求解方法——乘子法、MMA法及SQP 法的效果.

3 種求解方法得到的最優拓撲如圖3 所示,對應的Mises 應力分布如圖4 所示.迭代次數、體積比約束值、目標函數應力值及收斂精度對比如表1中所示.目標及約束函數的迭代歷史曲線如圖5 所示.從結果對比可以看出,在0.01 的收斂精度下,乘子法及SQP 法均能收斂到清晰的拓撲圖形,兩種方法下的最大Mises 應力接近(見表1),拓撲圖形也相似(見圖3(a)及圖3(c)),但SQP 的收斂速度更快.MMA 方法在0.01 的收斂精度下得不到清晰的拓撲圖形(見圖3(b)),相比于乘子法及SQP 法,其收斂時的Mises 應力也非常大(見表1 及圖5(b)).當將MMA法的收斂精度提高到0.001 時,經過217 次迭代收斂,其拓撲圖形及Mises 應力分布如圖6 所示,此時的最優拓撲與乘子法及SQP 法得到的最優拓撲是相似的(對比圖6(a) 與圖3(a) 及圖3(c)),Mises應力為191.192 MPa,與乘子法及SQP 法的對應值也是接近的(見表1).可見,MMA 法在提高一個量級收斂精度的情況下,才能得到與乘子法及SQP 法相同的結果.

圖3 算例1 最優拓撲Fig.3 Optimized topology of the Example 1

圖4 算例1 最優拓撲的Mises 應力分布Fig.4 Mises stress distribution of of the optimized topology for the Example 1

表1 m 方模型不同求解方法對比Table 1 Comparison of different methods to solve m-model

圖5 算例1 目標及約束迭代歷史Fig.5 The iteration history of the objectives and constraints for the Example 1

圖6 算例1 收斂精度0.001 時MMA 法得到的最優拓撲及Mises 應力分布Fig.6 Optimized topology and Mises stress distribution for the Example 1 by MMA method with the convergence accuracy 0.001

從此算例對比中可知:3 種方法中SQP 法的收斂速度最快,乘子法與之接近.而MMA 法的收斂速度遠低于乘子法及SQP 法,其原因是MMA 方法是一階近似的求解方法.由此,在本文在化整-集成方法中,對m 方集成模型的求解選用SQP 解法.

4.2 算例2:s 方模型化整交融應力拓撲優化

此算例的條件與算例1 相同.但是求解的優化模型為s 方模型,即應力約束下極小化結構體積的優化模型.結構允許應力為200 MPa.分別用單獨的化整解法及化整-集成解法兩種方法進行求解.

單獨的化整解法經過19 次收斂,最大Mises 應力為201.743 MPa,略微超出了允許應力200 MPa.化整-集成解法經過17 次收斂,最大Mises 應力為199.704 MPa.兩種方法得到的最優拓撲如圖7 所示,對應的Mises 應力如圖8 所示,目標及約束迭代歷史曲線如圖9 所示.兩種方法的最優拓撲及Mises應力分布均是接近的,化整-集成解法比單獨的化整解法收斂速度略快.

圖7 算例2 最優拓撲Fig.7 Optimized topology of the Example 2

圖8 算例2 最優拓撲的Mises 應力分布Fig.8 Mises stress distribution of the optimized topology for the Example 2

圖9 算例2 目標及約束迭代歷史.Fig.9 The iteration history of the objectives and constraints for the Example 2

4.3 算例3:多工況s 方模型化整交融應力拓撲優化

如圖10 所示,基結構平面體的尺寸為200 mm×100 mm×10 mm,采用 200×100 網格.材料彈性模量為 E=210 GPa,泊松比 ν=0.3 .左右邊界固定,工況1:集中載荷 F1=2600 N 作用于上邊界中點;工況2:集中載荷 F2=2600 N 作用于下邊界中點.劃分為200×100個矩形單元.應力約束為200 MPa,收斂精度取0.001.

圖10 兩工況算例分析及優化的模型Fig.10 Model for analysis and optimization of the example with two load cases

此兩工況應力約束拓撲優化問題在文獻[8-9]中采用單獨的化整解法進行過求解.在此用它來比較多工況情況下化整解法與化整-集成解法的效果.

化整解法經過41 次迭代收斂,最優體積比為35.895%,最大Mises 應力為200.771 MPa.化整-集成解法經過40 次迭代收斂,最優體積比為24.055%,最大Mises 應力為200.171 MPa.兩種解法得到的最優拓撲如圖11 所示,對應的兩個工況下的Mises 應力分布如圖12～ 圖13 所示.經對比可以看到:兩種解法得到的最優拓撲構型相差較遠.雖然兩種解法收斂速度差不多,但在同樣滿足應力約束的條件下,化整交融解法得到的最優拓撲總體積更小,表現出更好的尋優能力.

圖11 算例3 最優拓撲Fig.11 Optimized topology of the Example 3

圖12 算例3 化整方法最優拓撲的Mises 應力分布Fig.12 Mises stress distribution of of the optimized topology for the Example 3 by globalization approach

圖13 算例3 化整-集成方法方法最優拓撲的Mises 應力分布Fig.13 Mises stress distribution of the optimized topology for the Example 3 by globalization-integration approach

5 結論

對連續體結構應力約束拓撲優化問題,本文提出了求解m 方模型的乘子法及SQP 法,并與MMA法進行了對比,結果表明乘子法及SQP 法求解效率遠高于MMA 法.在此基礎上,構建了應力約束拓撲優化問題的化整交融解法,其中求解m 方集成模型選用效果最好的SQP 解法.與文獻[8-9]中原所提出的化整解法進行了對比.單工況及多工況算例均表明:化整交融解法求解效率與化整解法基本相當,但尋優能力更好,得到的最優拓撲結構在滿足應力約束條件下更輕.

順便指出:本文的方法對于其他類型的局部性能例如疲勞或局部穩定約束下的結構拓撲優化都是適用的,相關研究將另文發表.