非線性連續體拓撲優化方法綜述1)

文桂林 劉 杰 陳梓杰 魏 鵬 龍 凱 王洪鑫 榮見華 謝億民

* (燕山大學機械工程學院,河北秦皇島 066004)

? (廣州大學機械與電氣工程學院,廣州 510006)

** (華南理工大學土木與交通學院,廣州 510641)

?? (華北電力大學新能源學院,北京 102206)

***(長沙理工大學汽車與機械工程學院,長沙 410076)

??? (皇家墨爾本理工大學創新結構與材料研究中心,澳大利亞墨爾本 3001)

引言

連續體拓撲優化是指在給定邊界條件、載荷,以及各種約束條件下,根據性能指標要求,通過嚴格的數學理論得到結構最優材料布局的過程[1-3].它可以從力學上根本提升結構機械性能,其本質上屬于無窮維空間中變系數微分方程的0～1 控制問題,被公認為是最具挑戰性的結構優化問題.近年來,連續體拓撲優化在工程領域也展現出巨大發展潛力,已經在航空航天、海洋工程、機械工程等多個工程領域發揮了重要作用.圖1(a)～圖1(h)給出了連續體拓撲優化方法的典型工程應用.

圖1 連續體拓撲優化方法的典型工程應用: (a) 衛星支撐系統[4];(b) 海上風力發電機組導管架支撐結構[5];(c) 機翼結構[6];(d) 飛機擾流板[7];(e) 超長懸索橋[8];(f) 導彈氣動設計[9];(g) 飛機蒙皮拉伸成形模具[10];(h) 空客A320 機艙鉸鏈支架[11]Fig.1 Typical engineering application of continuum topology optimization method: (a) satellite support system[4];(b) offshore wind turbine jacket support structure[5];(c) wing structure[6];(d) aircraft spoiler[7];(e) super-long suspension bridge[8];(f) missile aerodynamic design[9];(g) drawing die for aircraft skin[10];(h) hinge bracket of Airbus A320 cabin[11]

均勻化方法開創了連續體拓撲優化的先河[12],其思想是將設計區域離散為大量單元,通過微孔結構幾何參數來確定單元的有無.由于具有設計變量多、優化最終構型為多孔結構導致難以制造等缺點,限制了其發展和應用,但是,其理論基礎嚴謹,且目前在材料設計領域還持續展現出活力[13-15].繼均勻化法之后,連續體拓撲優化方法得到了快速的發展,逐漸形成了變密度法,包括固體各向同性材料微結構/插值懲罰法(solid isotropic macrostructure/material with penalization method,SIMP)[16-17]和材料性能的合理近似法(rational approximation of material properties,RAMP)[18]、漸進結構優化方法(evolutionary structural optimization,ESO)[19-21]、水平集法(level set,LS)[22-24]、獨立-連續-映射法[25-26]、特征驅動法[27-28]、相場法[29-30]、拓撲導數法[31]以及移動可變形組件/空洞法[32-34]等.這些方法各有優缺點,共同促進了連續體拓撲優化領域的蓬勃發展.圖2給出了典型的連續體拓撲優化方法設計一個二維懸臂梁的最終構型對比圖,可見,各種方法均可產生清晰的拓撲構型.

圖2 懸臂梁設計問題及不同拓撲優化方法的優化型Fig.2 Cantilever design problem and optimization with different topology optimization methods

變密度法、漸進結構優化方法和獨立-連續-映射法由于采用像素點描述結構構型,會有結構邊界不光滑的問題.此外,變密度法會產生灰度單元,這些都增加了實際制造難度,不過可以通過細分單元、使用新的過濾策略等方式進行緩解.比如,文桂林團隊[35-36]提出了結構邊界單元自適應策略和積木博弈策略,有效地緩解了BESO (bi-directional evolutionary structural optimization,BESO)方法和SIMP 方法中結構邊界不光滑的問題,降低了后處理時邊界經驗處理帶來的誤差.變密度法發展較早,相對于其他方法在理論層面上較為成熟;而謝億民團隊[19-21]提出的漸進結構優化方法具有算法理念簡單、收斂速度快和邊界清晰等優點,且成功開發了商用軟件Ameba,已在建筑工程、機械工程和航空航天等領域的創新設計方面展現出來強大的實際工程解決能力[37].水平集法、特征驅動法、相場法、拓撲導數法和移動可變形組件/空洞法屬于基于邊界演化發展的拓撲優化方法,這類方法的優勢在于能夠得到光滑和清晰的優化結構邊界,且幾何信息可以通過函數表示,易于提取最終構型的幾何信息;然而,由于計算通常基于固定的有限元網格,動響應和靈敏度的計算精度方面有一定的提升空間.張衛紅團隊[27-28]提出的特征驅動拓撲優化方法將結構特征貫穿于結構建模、分析和優化設計的整個流程,即,將CAD 特征模型嵌入到結構拓撲模型,從源頭上保證了結構的特征屬性.郭旭團隊[32-34]提出的移動可變形組件/空洞法給出了一種顯式的水平集演化模式,大幅度地減少了設計變量的數目,提高了計算效率,且優化結果可與CAD 軟件進行直接連接,易于制造.此外,作為LS 法的一種,參數化水平集方法可以改善傳統水平集方法的初始設計依賴性缺陷,然而,由于材料分布不存在中間密度,材料密度突變易造成不連續的拓撲變化,致使優化過程不穩定.基于此,Wei 等[38]結合變密度法的優勢,提出了水平集帶方法.該方法的思路是,利用具有一定寬度的水平集帶范圍替代傳統零水平集面,水平集帶區域內的材料密度通過插值可以連續地分布于[0,1]范圍內.

目前連續體拓撲優化在解決線性問題中發展較為成熟,已經成功解決了諸如傳熱、增材制造、振動抑制、穩健性和可靠性設計等問題[39-50].然而,實際工程中存在大量的非線性問題,若近似為線性問題來解決,會導致較大的誤差,甚至是錯誤的結果.實際結構設計中的非線性主要體現為三種: 材料非線性、幾何非線性和邊界非線性(接觸非線性).圖3 分別為三種非線性問題的示意圖.圖3(a)為材料非線性問題,其中材料的應力應變之間不再是線性關系,常見的如彈塑性、超彈性等;圖3(b)為幾何非線性問題,當結構發生幾何大變形、大撓度或大轉動時,力和位移曲線不再是線性關系,即剛度矩陣不再是常數;圖3(c)為邊界非線性問題,結構邊界在分析過程發生變化.基于此,本文將分別對材料非線性、幾何非線性和邊界非線性三種連續體拓撲優化方法進行綜述.需要指明的是,本文重點在宏觀力學范疇內對非線性連續體拓撲優化進行綜述,對于考慮微觀結構的非線性拓撲優化問題可以參考文獻[51].

圖3 三種結構非線性Fig.3 Three kinds of structural nonlinearity

1 材料非線性連續體拓撲優化

材料非線性主要體現為材料應力-應變關系的非線性,在實際工程中比較常見的非線性材料為彈塑性材料和超彈性材料.當彈塑性材料受力變形導致內部應力超過比例極限時,其應力-應變將呈現非線性關系,隨著應力水平的進一步提升,當應力超過屈服應力時,將產生不可恢復的塑性變形.而超彈性材料具有更為復雜的非線性應力-應變關系,在大變形狀態下卸載可恢復初始狀態.如何將不同類型非線性材料的應力-應變關系引入到拓撲優化列式,以及求解由此為計算帶來的系列難題,是材料非線性連續體拓撲優化區別于線性問題的關鍵.本部分將以彈塑性材料和超彈性材料兩個類別對材料非線性連續體拓撲優化相關研究進行綜述.

1.1 彈塑性材料連續體拓撲優化

Maute 等[52-53]較早地在SIMP 法框架下,利用Von Mises 應力屈服準則推導了彈塑性材料塑性部分的非線性應力-應變關系,其中需要分別建立彈塑性材料的楊氏模量、硬化模量和屈服應力的插值模型,然后將塑性部分與彈性部分相結合,得到彈塑性材料的應力-應變關系,該方法需要權衡三種插值模型的懲罰系數,見圖4(a).在此基礎上,Wallin 等[54]發展了彈塑性結構能量吸收最大化的拓撲優化方法,見圖4(b);Kato 等[55]研究了體積約束下彈塑性復合材料能量吸收能力最大化的拓撲優化方法,見圖4(c);Li 等[56-57]提出了抗斷裂的吸能彈塑性結構拓撲優化設計和循環荷載作用下的耗能彈塑性結構拓撲優化設計方法,見圖4(d);Blachowski 等[58]研究了應力約束下彈塑性結構拓撲優化新方法.彈塑性材料模型可簡化為由楊氏模量、硬化模量和屈服應力確定的雙線性應力-應變形式,楊氏模量和硬化模量分別為第一線性段和第二線性段的斜率,屈服應力作為兩個線性段的分段應力值.但是,同樣需要分別建立三個相關材料屬性的插值模型.以此為基礎,Mayer 等[59]和Patel 等[60]研究了彈塑性材料拓撲優化設計在吸能方面的應用;Guimar?es 等[61]將拓撲優化應用在采用血管成形術的支架細胞設計中;Lee 等[62]采用等效靜載荷法高效地解決了非線性拓撲優化設計問題;Yoon 等[63]發展了一種單元連接參數化方法,為彈塑性材料拓撲優化設計提供了一種新的有效思路,但是該方法也會增加額外的剛度矩陣計算,一定程度上犧牲了計算效率,見圖4(e);豆麟龍等[64]發展了一種基于漸進結構優化方法的彈塑性材料拓撲優化設計方法;Luo 等[65]提出了一種移動等值面閾值法,解決了彈塑性材料的拓撲優化設計問題,該方法允許刪除的孔洞單元參與設計變量的更新,并在隨后的迭代中重新出現,因此,具有相對較高的計算精度;Kongwat 等[66]發展了基于比例法的彈塑性材料拓撲優化方法.此外,在均勻化法框架下,Yuge 等[67]建立了一個材料張量的二維數據庫,以解決彈塑性材料拓撲優化問題.Bogomolny等[68]采用Drucker-Prager 屈服準則建立彈塑性模型,提出了一種針對鋼筋混凝土結構的概念設計新方法.以最大損傷為約束,Li 等[69]提出了一種結構最優能量吸收的拓撲優化方法;隨后,基于非局部彈塑性損傷模型,他們發展了結構抗破壞拓撲優化方法,所采用的非局部彈塑性損傷模型通過Von Mises 應力屈服準則與一個非局部損傷函數耦合實現[70],見圖4(f).Zhao 等[71]通過構造一個虛構的非線性彈性模型,并利用漸近分析的方法考慮滿足Von Mises 應力屈服準則的非線性材料響應,并以此建立了兩類優化列式,即應變能最大化和載荷系數最大化,見圖4(g).此外,徐斌團隊近期在BESO 框架下,發展了考慮應力約束的材料非線性拓撲優化方法[72-73].需要指出,除了具有非線性應力-應變關系,加載路徑相關性是彈塑性材料連續體拓撲優化的另一個特點.該特點導致求解優化目標函數和約束函數對設計變量的靈敏度信息變得困難,其敏度分析過程需采用考慮路徑相關的伴隨求解策略.可見,彈塑性材料連續體拓撲優化相關研究中,通常采用Von Mises 應力屈服準則推導彈塑性材料塑性部分的非線性應力-應變關系,也有一些研究使用簡化的雙線性模型.對于一些更加復雜的問題,可對Von Mises 應力屈服準則塑性材料模型進行進一步修正.

圖4 彈塑性非線性拓撲優化: (a) 彈塑性和彈性拓撲優化設計簡支梁結構對比[52];(b) 有限塑性變形最大塑性應變能最大化設計[54];(c) 彈塑性吸能最大簡支梁結構優化設計[55];(d) 抗斷裂吸能L 型結構設計[56];(e) 單元連接參數化策略示意圖[63];(f) 抗材料損傷能量吸收結構優化設計[70];(g) 考慮應變能最大和載荷系數最大的材料非線性拓撲優化設計[71]Fig.4 Elastoplastic nonlinear topology optimization: (a) comparison between elastoplastic and elastic topology optimization design of a simply supported beam structure[52];(b) maximum plastic strain energy design with finite plastic deformation[54];(c) optimal design of a simply supported beam structure with maximum elastic-plastic energy absorption[55];(d) design of a fracture resistant and energy absorbing L-shaped structure[56];(e) schematic diagram of element connection parameterization strategy[63];(f) optimal design of energy absorbing structure resistant to material damage[70];(g) nonlinear topology optimization design of materials considering maximum strain energy and maximum load coefficient[71]

1.2 超彈性材料連續體拓撲優化

相對于彈塑性材料模型,超彈性材料具有更為復雜的應力-應變關系,通常采用應變能函數描述其本構關系[74],比較經典的包括Mooney-Rivlin 模型、Neo-Hookean 模型、Kirchhoff-St.Venant 模型和Ogden 模型等.本部分根據本構模型分類綜述超彈性材料連續體拓撲優化方法.需要指出的是,超彈性行為往往包含了幾何非線性,關于幾何非線性拓撲優化相關問題和研究進展見第2 部分.

基于Mooney-Rivlin 模型,Lee 等[75]對考慮靜動態特性的橡膠隔振器進行拓撲優化設計,見圖5(a);文獻[76-77]在水平集方法框架下研究了考慮幾何非線性的超彈結構拓撲優化方法,見圖5(b);Deng 等[78]認為在設計大變形超彈性材料時,Von Mises 應力破壞準則不準確,進而提出一種基于畸變能的連續體拓撲優化方法(見圖5(c)),其中在低密度區域,利用了Wang 等[79]提出的應變能插值方法;Ortigosa 等[80-81]解決了SIMP 和LS 框架下超彈性材料大應變拓撲優化設計中算法不穩定問題;薛日野等[82]推廣了移動可變形空洞法,實現了超彈結構拓撲優化設計,見圖5(d).利用Neo-Hookean 模型表征超彈本構模型,基于多凸本構模型與松弛函數法,Lahuerta 等[83]實現了超彈結構的幾何非線性拓撲優化,該方法在優化過程中無需去除孔洞單元,具有較高計算精度;James 等[84]將連續體拓撲優化方法應用到超彈性雙穩態心血管支架的優化設計,見圖5(e);充分利用無網格法解決大變形問題,Zhang 等[85]提出了一種有限元法與無網格伽遼金法結合的超彈結構拓撲優化方法,見圖5(f);Wallin 等[86]研究了在幾何非線性拓撲優化中使用割線剛度和切線剛度對最終優化構型的影響,研究發現,在低荷載水平下,兩種結構構型差別不大,而在大變形情況下,則存在顯著差異;受附加超彈材料技術的啟發[87],張憲民團隊[88]發展了與商業化軟件ANSYS 結合的超彈結構幾何非線性拓撲優化方法,并提供了MATLAB 程序,良好地提升了幾何非線性拓撲優化方法應用實際工程結構設計的能力,見圖5(g);在此基礎上,Ye 等[89]在ICM 框架下,解決了應力約束的非線性拓撲優化設計.然而,該方法在解決三維非線性問題時會存在不收斂的問題,基于此,受文獻[87-88]啟發,文桂林團隊[90]近期發展了一種基于多分辨率框架的超彈材料拓撲優化設計方法,用單核PC 機成功解決了千萬分辨率的三維非線性結構優化設計問題,該方法具有高效和保精度的優點,并可以很好地解決三維非線性問題,見圖5(h);Capasso 等[91]提出了應力約束下的柔性機構非線性拓撲優化方法;Zhang 等[92]改進了BESO方法,解決了發生大變形的超彈材料設計問題,見圖5(i).基于Kirchhoff-St.Venant 模型,Bruns 等[93-94]提出了一種適定的考慮幾何大變形的連續體拓撲優化數學模型,解決了結構非線性拓撲優化設計問題,并應用到柔性機構設計;Klarbring 等[95]研究了包含非零位移約束下的超彈結構拓撲優化問題,對比了多種超彈性材料模型對最終拓撲材料分布的影響,發現基于Kirchhoff-St.Venant 模型所得的優化構型與其他包含物理壓縮行為的模型有所差異,且力學性能較差;基于變密度法、水平集法和擴展有限元法,Geiss 等[96]發展了一種可以優化設計大變形4D 打印結構的有效方法.然而,對于有些實際具體工程問題,需要對經典超彈性材料本構模型進行修正.Luo 等[97]研究了在無摩擦接觸支承下超彈結構的拓撲優化問題,其中超彈性材料的應變能函數通過人工懲罰模型表示.為了考慮超彈性材料的非線性,Chi 等[98]在SIMP 法框架中提出一種改進的Ogden通用本構模型,該模型可同時考慮平面應變和平面應力兩種情況,并通過改變材料參數能夠使材料在平面應變和平面應力情況下獲得一定范圍的非線性材料行為.綜上所述,超彈性材料具有十分復雜的應力-應變關系,對于基于有限單元法的連續體拓撲優化的優化列式建立、靈敏度推導以及算法穩定性都是一個挑戰.另外,在實際工程問題中,只通過單軸拉伸試驗往往無法得到足夠的樣本點來表征超彈材料行為,與經典本構模型有較大誤差,需要同時開展單軸試驗、雙軸試驗、平面試驗、體積試驗等.

圖5 超彈性材料拓撲優化: (a) 橡膠隔震器拓撲優化[75];(b) 水平集法線性和非線性優化設計懸臂梁結果對比[77];(c) 基于畸變能的拓撲優化設計[78];(d) 基于移動可變形孔洞方法設計懸臂梁結構[82];(e) 雙穩態心血管支架優化設計[84];(f) 有限元和無網格耦合拓撲優化方法設計懸臂梁結構[85];(g) 柔性機構幾何非線性設計[88];(h) 高分辨率三維非線性拓撲優化設計類橋梁結構[90];(i) 改進的BESO方法優化設計懸臂梁結構[92]Fig.5 Topology optimization of hyperelastic materials: (a) topology optimization of rubber isolator[75];(b) comparison of linear and nonlinear optimization design of a cantilever by level set method[77];(c) topology optimization design based on distortion energy[78];(d) design of a cantilever structure based on MMV method[82];(e) optimal design of a bistable cardiovascular stent[84];(f) design of a cantilever structure by finite element and meshless coupling topology optimization method[85];(g) geometric nonlinear design of a flexible mechanism[88];(h) high resolution 3D nonlinear topology optimization design of bridge structure[90];(i) optimization design of a cantilever structure by improved BESO method[92]

2 幾何非線性連續體拓撲優化

對使用過程中會發生大位移或大轉動的工程結構進行連續體拓撲優化設計時,需考慮幾何非線性的影響.幾何非線性連續體拓撲優化問題中,有限元響應求解的剛度矩陣具有位移相關性,常采用Newton-Raphson 迭代求解其非線性平衡方程.基于變密度法的幾何非線性連續體拓撲優化中,中低密度單元區域剛度較低,當結構發生大變形時,這些材料區域易發生過度畸形,致使Newton-Raphson 迭代中的切線剛度矩陣失去其正定性質,無法求解非線性平衡方程,進而導致非線性拓撲優化無法收斂.而其他經典拓撲優化方法與變密度法類似,空洞區域的低剛度同樣會帶來數值求解困難.若將結構中除實體單元外的其他單元統稱為“偽材料”單元,幾何非線性連續體拓撲優化的難題主要在于如何解決單元過度畸形導致的數值不穩定.幾何非線性包括大位移小應變和大位移大應變兩種情況,前者稱為弱幾何非線性,而后者稱為強幾何非線性.目前關于幾何非線性連續體拓撲優化的相關研究大部分針對弱幾何非線性情況,若無特別說明,本文默認指弱幾何非線性情況.本部分分別基于“偽材料”單元、單元畸形和切線剛度矩陣非正定這三個特征對幾何非線性連續體拓撲優化相關研究進行綜述.

2.1 基于排除“偽材料”單元思路的幾何非線性連續體拓撲優化

在拓撲優化迭代循環過程中人為地將“偽材料”單元排除在外,可以有效避免非線性有限元計算失敗.基于這一思想,Buhl 等[99]認為低密度單元節點位移的持續振蕩是導致Newton-Raphson 迭代無法收斂的原因,他們在SIMP 法框架中提出“收斂準則松弛”法,將被低密度單元包圍的節點排除在Newton-Raphson 迭代過程的收斂準則之外,見圖6(a).隨后將該方法應用于大位移柔性機構拓撲優化設計[100],研究發現,與線性拓撲優化構型相比,考慮幾何非線性拓撲優化構型的輸出性能增益可達2.5 倍.基于該思路,李兆坤等[101]實現了多輸入多輸出的柔性機構幾何非線性拓撲優化設計,見圖6(b);在參數化水平集框架下,Luo 等[102]進一步將該思想推廣到求解大位移柔性機構拓撲優化問題;為了解決幾何缺陷對結構穩定性的影響,Jansen 等[103]推廣了文獻[99] 的“收斂準則松弛”法,發展了一種穩健性拓撲優化方法.基于SIMP 法,Bruns 等[104]提出一種“單元移除法”,該方法可將刪除的低密度單元重新引入的策略,根據迭代歷史中單元密度是否低于所規定密度這一規則對單元進行計數,若該單元在此次迭代中密度小于所規定密度,則次數減1,反之加1,并通過單元次數與所設定閾值的比較判斷是否對單元狀態進行切換,見圖6(c).借鑒“單元移除法”,Luo 等[65]利用移動等值面設置閾值,低密度單元在有限元計算中被消除,而在計算結構響應時被重新引入.Cho 等[105]認為外力的加載方式帶來較多“偽材料”單元導致優化過程不容易收斂,因此,他們采用了位移加載方式來緩解該問題.基于無網格法,Cho 等[106]通過重構核形狀函數對位移場和體密度場進行離散,使形狀控制點能夠位于設計域任何部分,這種做法可將低密度域排除在有限元響應分析之外,解決了幾何非線性拓撲優化數值求解數值不穩定問題,見圖6(d);在SIMP 法框架下,Du等[107]利用無網格伽遼金法將單元設計域轉化為節點設計域,發展了一種新的幾何非線性拓撲優化方法,并應用于大位移柔性機構優化設計,該方法由于不需要單元進行插值,可避免中低密度區域單元過度畸變問題;在水平集法框架下,Yamada 等[108]提出了一種利用粒子交互模型結合移動粒子半隱式思路,成功解決了幾何非線性拓撲優化數值困難;基于無網格伽遼金法,亢戰團隊[109]提出了一種獨立密度場插值的幾何非線性拓撲優化方法,該方法采用了能量收斂準則,有效解決了由于低密度局部數值不穩定帶來的收斂困難.此外,在BESO 框架下,Huang等[110-112]系統地研究了位移荷載和力載荷作用下幾何非線性拓撲優化問題,見圖6(e).值得提出的是,徐斌團隊近年來在BESO 框架下,進一步發展了幾何非線性拓撲優化方法,成功解決了應力約束等問題[113-114].此外,張衛紅團隊[115]提出了一種考慮連接點載荷約束和幾何非線性的裝配體結構拓撲優化設計方法,他們采用超單元壓縮方法對低剛度區域自由度進行壓縮,以避免算法不收斂問題.

圖6 基于排除“偽材料”單元思路的幾何非線性拓撲優化: (a) “收斂準則松弛”思路(白色單元四周的節點被排除)[99];(b) 多輸入多輸出柔順機構優化設計[101];(c) 使用“單元移除法”前后設計柔性結構構型對比[104];(d) 基于無網格法進行懸臂梁設計[106];(e) 基于BESO 的幾何非線性拓撲優化設計[111]Fig.6 Geometrically nonlinear topology optimization based on eliminating "pseudo material" elements: (a) "convergence criterion relaxation" idea (nodes around white elements are excluded)[99];(b) optimal design of a multi input and multi output compliant mechanism[101];(c) comparison of flexible structure configuration before and after using "element removal method"[104];(d) design of a cantilever based on meshless method[106];(e) geometrically nonlinear topology optimization based on BESO[111]

綜上所述,“收斂準則松弛”法發展較早,且思想簡便,應用較為廣泛,但其為了實現數值穩定性,直接忽略了低密度單元的有限元響應,導致最終優化結果精度不高.“單元移除法”直接將單元進行移除,同樣影響優化結果精度,而低密度單元在有限元計算中被消除,在計算結構響應時被重新引入的策略可以一定程度提高計算精度.考慮到低密度單元數值不穩定出現在非線性有限元求解環節,而并非拓撲優化環節,無網格分析技術提供了另一種有效的手段,具有較高精度,但存在計算量大、效率較低等缺點.

2.2 基于改變單元畸形思路的幾何非線性連續體拓撲優化

解決幾何非線性連續體拓撲優化數值不穩定性的另一個思路是將單元的過度變形轉化為其他方式體現或者對單元變形進行縮放,以抑制單元過度畸形.基于這一思想,Yoon 等[116-117]利用單元連接參數化法分別在SIMP 和水平集法框架下進行了幾何非線性拓撲優化設計,見圖7(a)和圖7(b).單元連接參數化法使用零長度一維彈性連桿連接相鄰單元,利用彈性連桿變形代替平面單元變形,由于平面單元保留了原始材料特性,可在結構發生較大變形時避免過度畸變.彈性連桿的存在表示相鄰單元的連接,最終保留優化迭代結束時相互連接的單元,這種方法稱為外部單元連接參數化方法,與之對應,他們后來又發展了內部單元連接參數化方法,并擴展到基頻優化問題[118-119].在此思想基礎上,Moon 等[120]發展了單元連接參數化qp 松弛法,以解決Von Mises 應力約束下的體積最小化幾何非線性拓撲優化設計問題.在水平集法框架下,利用Delaunay三角剖分技術,Ha 等[76]提出了一種網格重劃方法,并采用超彈性材料,一定程度上避免了由于大應變引起的單元內部過度變形.此外,網格局部細化策略[35]也可以一定程度上緩解幾何非線性拓撲優化過程中的單元過度扭曲問題.van Dijk 等[121]基于SIMP 法提出了單元變形縮放技術,將單元的節點位移分離成保留的旋轉剛體部分和允許縮放的變形部分,根據單元密度大小對其內部單元位移進行插值實現單元變形縮放,從而避免低密度單元過度畸形,見圖7(c).

圖7 基于改變單元畸形思路的幾何非線性拓撲優化: (a) 單元連接參數化法在SIMP 框架下懸臂梁設計[116];(b) 單元連接參數化法在水平集框架下懸臂梁設計[117];(c) 利用單元變形縮放技術進行懸臂梁設計[121]Fig.7 Geometrically nonlinear topology optimization based on the idea of changing element deformity: (a) designing a cantilever under SIMP framework by element connection parameterization method[116];(b)design of a cantilever under level set frame by element connection parameterization method[117];(c) cantilever design using the element deformation scaling technology[121]

由此可見,基于改變單元畸形思路的幾何非線性連續體拓撲優化方法可以一定程度上緩解優化過程中的數值不穩定難題,但在計算效率和計算精度上仍需要進一步的提升.單元連接參數化方法在幾何非線性拓撲優化問題求解中展現出了良好的潛能,但亦存在額外的剛度矩陣計算導致其計算效率較低等問題.網格重劃分和局部網格細化策略可以一定程度緩解數值不穩定難題,但是較難從根本上解決.單元變形縮放技術中單元內力對外部位移場比較敏感,致使其處理幾何非線性拓撲優化問題時有一定的局限性.

2.3 直接影響切線剛度矩陣非正定的幾何非線性連續體拓撲優化

直接影響切線剛度矩陣非正定特性是另一種解決幾何非線性連續體拓撲優化數值不穩定的思路.Bruns 等[94]建議額外添加超彈性材料防止單元過度畸形導致的數值收斂困難.在多材料多輸出大位移柔性機構拓撲優化設計中,Saxena[122]利用遺傳算法隱式地規避計算非線性變形時的不收斂解問題,見圖8(a).Cho 等[123]在水平集法框架下,采用均勻材料屬性和隱式邊界轉換為顯式邊界的策略實現了幾何非線性求解的數值穩定性.Kawamoto[124]利用Levenberg-Marquardt 法替代Newton-Raphson法對非線性方程進行迭代求解,避免了優化迭代過程中的不收斂問題,見圖8(b).基于等效靜載荷思路,Lee 等[62]將非線性靜態響應近似等效為某一些時刻靜載荷的線性靜態響應,有效解決了幾何非線性連續體拓撲優化設計問題.為了保證切線剛度矩陣的正定性,Lahuerta 等[83,125]采用多凸性本構模型結合松弛函數實現了過度變形單元的穩定.Yoo 等[126]使用改進的蟻群優化算法實現了幾何非線性拓撲優化.Gomes 等[127]將序列分段線性規劃算法引入到幾何大變形結構拓撲優化問題求解.Wang 等[79]基于SIMP 法提出一種新的能量插值策略,在低剛度單元區域和實體單元區域分別采用小變形理論和大變形理論對彈性能密度進行插值,保證了切線剛度矩陣正定性.后來,Lazarov 等[128]將該方法應用到微機電系統的大位移柔性機構的穩健性拓撲優化設計.Wallin 等[86]比較了使用割線剛度和切線剛度對最終優化構型的影響,他們發現,當結構發生小變形時,最終構型基本一致,而當發生大變形時,則會出現明顯的區別.Zhu 等[129]發展了移動可變形組件法,成功解決了幾何非線性拓撲優化問題,見圖8(c);Silva 等[130]提出考慮應力約束、制造不確定性和幾何非線性的柔性機構拓撲優化方法;Zhu 等[131]使用開源軟件平臺FreeFEM 編寫了幾何非線性拓撲優化MATLAB 程序.在SIMP 法框架下,Luo 等[87]提出附加超彈材料技術,在中低密度區域中加入一種與原始單元共用節點的超彈性單元,通過所附加超彈性材料的高體積模量和剪切能力可以有效地緩解中低密度單元畸形導致的數值失穩問題.基于附加超彈材料技術,Chen 等[88]公布了幾何非線性結構拓撲優化MATLAB 程序,將幾何非線性拓撲優化算法與商業分析軟件ANSYS 結合,一定程度上推進了拓撲優化算法的工程實用性.在此基礎上,龍凱等[132]提出了一種控制結構最大位移的幾何非線性拓撲優化方法,具有良好的工程實用性,見圖8(d).為了進一步提高優化結果的精度,Liu 等[133]提出了一種改進的附加超彈性技術,并實現了大位移柔性機構的幾何非線性拓撲優化設計,見圖8(e).此外,Dunning[134]提出一種基于共旋轉法的幾何非線性拓撲優化方法,這種方法的優點在于切線剛度矩陣具有固有正定性,見圖8(f).需要指出的是,柔性機構一般利用幾何大變形來實現大的位移或力傳遞,其拓撲優化需要考慮幾何非線性,然而也應考慮諸如多輸入多輸出需求、類鉸鏈、灰度單元等問題[2].比如,榮見華團隊[135]近期提出了一種考慮多輸入多輸出的無鉸鏈全解耦柔性機構拓撲優化的新方法以解決上述問題,可將該方法進一步推廣到考慮幾何非線性拓撲優化設計問題.

圖8 基于直接影響切線剛度矩陣非正定思路的幾何非線性拓撲優化: (a) 利用遺傳算法優化大位移柔性機構[122];(b) 基于Levenberg-Marquardt 法優化直線運動機構[124];(c) 基于MMC 的懸臂梁優化[129];(d) 控制最大位移的懸臂梁設計[132];(e) 改進的附加超彈性技術優化設計柔性夾持機構[133];(f) 共旋轉優化設計不同的位移反向器[134]Fig.8 Geometrically nonlinear topology optimization based on the idea of directly affecting the non positive definite tangent stiffness matrix:(a) optimization of a large displacement flexible mechanism using genetic algorithm[122];(b) optimization of a linear motion mechanism based on Levenberg-Marquardt method[124];(c) optimization of a cantilever based on MMC method[129];(d) design of a cantilever to control maximum displacement[132];(e) optimization design of a flexible clamping mechanism based on improved additional hyperelasticity technology[133];(f) optimal design of different inverter mechanisms with co-rotational method[134]

綜上所述,上述方法可以有效影響切線剛度矩陣非正定特性,進而一定程度上緩解幾何非線性拓撲優化數值不穩定問題,但是,也存在進一步改進的空間.比如,能量插值策略在低密度單元區域的線性有限元分析會一定程度上影響最終優化結果的精度;附加超彈材料技術額外地增加超彈性單元同樣會影響優化結果的精度,且關鍵參數的經驗選取對最終結果有一定的影響;Levenberg-Marquardt 法的實用性取決于迭代過程中啟發式更新方案的選擇;序列分段線性規劃算法需要較多的計算成本;共旋轉法由于需要進行額外的矩陣轉換計算,計算效率較低.需要指出的是,基于排除“偽材料”單元和基于改變單元畸形思路來解決幾何非線性拓撲優化數值不穩定問題的核心還是為了解決非線性有限元計算過程中切線剛度矩陣非正定問題,所以,上述三種策略在本質上是一樣的.

3 邊界非線性連續體拓撲優化

實際工程結構系統中,各結構之間往往發生接觸,接觸體的變形和接觸邊界的摩擦使得部分邊界條件伴隨加載過程而發生非線性變化,諸如此類結構拓撲優化設計需要考慮邊界非線性.邊界非線性一般發生在具有接觸邊界條件的結構變形問題中,如汽車碰撞、齒輪的嚙合、電機組合轉子的組裝和軸承接觸等.本部分對邊界非線性連續體拓撲優化相關研究進行綜述.

通過將次梯度優化算法引入到拓撲優化中,Petersson 等[136]優化了單邊接觸條件下板結構,使其剛度最大,該方法具有較好的收斂性和較高的計算效率,見圖9(a).Mankame 等[137-138]借助正則化接觸模型對輔助接觸式柔性機構進行了拓撲優化設計,正則化接觸模型利用平滑逼近的單邊位移約束處理接觸相互作用;后來,針對具有給定非光滑路徑的輔助接觸式柔性機構,他們提出了一種考慮大變形的正則化法向接觸模型順序優化方法,以克服該機構的非光滑響應求解難題,見圖9(b).在接觸邊界條件下,Fancello[139]研究了考慮局部失效的質量最小多工況拓撲優化問題,在每個迭代步和加載工況下采用接觸求解器以獲得平衡變形構型.為了研究考慮結構件相互接觸時的拓撲優化設計問題,Peng 等[140]在給出接觸模型求解算法的基礎上,建立含有接觸約束的拓撲優化列式.Desmorat[141]利用均勻熱力學勢的概念,引入接觸模型,建立了無摩擦單邊接觸柔度最小化的拓撲優化列式,并通過固定網格和局部靈敏度計算策略進行有效求解,見圖9(c).為了求得具有給定摩擦的單邊接觸彈性體拓撲優化問題的數值解,My?liński[142-143]利用控制位移場的二階橢圓型變分不等式描述考慮摩擦的接觸問題,并通過拓撲導數和水平集法進行求解,見圖9(d).在SIMP 框架下,Str?mberg 等[144]開展了無摩擦單邊剛性支座接觸的連續體拓撲優化問題,其中,利用增廣拉格朗日法和光滑近似方法對Signorini 接觸條件進行改進,以描述無摩擦接觸問題.隨后,為了提高該方法的計算效率,Str?mberg[145]提出了一種以最大勢能為目標的優化列式.他還進一步將該思想擴展到解決具有制造約束和單邊接觸約束的結構連續體拓撲優化問題[146],見圖9(e).基于線性等效載荷法思路,Yi 等[147]發展了考慮接觸非線性的連續體拓撲優化方法.基于顯式水平集法和擴展有限元法,在無限小應變和線彈性假設下,Lawry 等[148]提出了一種解決彈性體之間存在滑移接觸的連續體拓撲優化方法,見圖9(f);在此基礎上,他們又將該方法擴展到解決有限應變問題[149],見圖9(g).利用材料掩膜覆蓋策略,Kumar等[150-151]對考慮自接觸和相互接觸兩種工況的柔性機構進行了拓撲優化設計.基于他們前面提出的附加超彈材料技術[87],并通過建立非線性彈簧模型模擬無摩擦接觸支承,Luo 等[97]進一步研究了無摩擦接觸支承下超彈結構的拓撲優化設計問題,研究發現,超彈結構發生大變形將導致結構的接觸邊界條件改變.針對無摩擦單邊接觸條件下結構剛度最大拓撲優化問題,Kanno[152]從拉格朗日對偶理論角度出發,建立了一種新的優化列式.基于三場SIMP 模型,通過引入接觸壓力方差約束,Niu 等[153]發展了一種無摩擦接觸彈性結構的拓撲優化方法,該方法可提高接觸壓力均勻性,見圖9(h).在假設接觸區域已知的前提下,Kristiansen 等[154]利用p范數的目標函數控制接觸壓力分布和變化,同時,為了控制接觸壓力,他們在耦合牛頓解中采用了基于拉格朗日乘子的接觸公式.Ma 等[155]提出了一種解決接觸問題的等效靜力位移法,該方法通過在靜力分析中施加一組位移模擬接觸力,所施加位移組的反力等于非線性分析中的接觸力.但是,該方法局限于尺寸優化問題.針對多個三維可變形物體互相接觸的拓撲優化問題,Fernandez 等[156]使用砂漿逐段接觸算法離散接觸面,該方法在模擬大變形固-固接觸時具有較好的魯棒性,見圖11(i).基于庫侖摩擦定律描述摩擦行為,Niu 等[157]發展了一種考慮摩擦接觸的彈性結構剛度最大化設計拓撲優化方法.最近,張衛紅團隊[158]研究了具有最大接觸壓力約束的彈性接觸的拓撲優化設計問題;他們采用了KS 凝聚函數表征一定接觸區域的最大接觸壓力,并將之作為約束引入到拓撲優化列式.

圖9 邊界非線性連續體拓撲優化: (a) 利用次梯度優化算法設計單邊接觸條件下板結構[136];(b) 借助正則化接觸模型優化設計輔助接觸式柔性機構[137];(c) 單側無摩擦接觸結構剛度優化結構[141];(d) 分段常數水平集法解決單側接觸問題[143];(e) 考慮制造約束和單邊接觸約束的優化設計[146];(f) 考慮滑動接觸界面的拓撲優化設計[148];(g) 有限應變下雙邊接觸拓撲優化結果[149];(h) 可提高接觸壓均勻性的拓撲優化設計[153];(i) 多個三維可變形物體互相接觸拓撲優化設計[156]Fig.9 Topology optimization with boundary nonlinearity: (a) utilizing sub-gradient optimization algorithm to design the plate structure under the condition of unilateral contact[136];(b) optimization design of an auxiliary contact flexible mechanism with regularized contact model[137];(c) stiffness optimization of one side frictionless contact structure[141];(d) solving one side contact problem by the piecewise constant level set method[143];(e) optimal design considering manufacturing constraints and unilateral contact constraints[146];(f) topology optimization design considering sliding contact interface[148];(g) topology optimization with bilateral contact under finite strain[149];(h) topology optimization design to improve the uniformity of contact pressure[153];(i) topology optimization design of a multiple 3D deformable object in contact with each other[156]

由此可見,接觸帶來的強非線性和不連續性對優化建模和求解算法提出了巨大挑戰.目前研究已經初步解決了基本的結構優化問題,但是,還存在一定的局限性,比如,主要局限于剛度最大設計問題,需要進一步向其他優化設計目標問題擴展,以適應實際工程需求;通常假設接觸區域已知,這與實際工程問題往往有較大出入,亟需進一步的研究.

值得提出的是,近期羅陽軍等[159-160]提出了一種基于材料場級數展開的結構非梯度拓撲優化方法,該方法通過模型降階的思想可為非線性拓撲優化問題(包括材料非線性、幾何非線性和接觸非線性)提供一個全新的求解思路,已成功解決了類如軟體聲子晶體、柔性電子和薄膜等非線性結構優化設計,以及動力學優化設計等連續體拓撲優化領域充滿挑戰性的問題[161-164].此外,引入其他力學約束或者性能需求的非線性拓撲優化問題將變得更加復雜,但也是解決結構實際工程結構優化設計問題必須要考慮的.比如,工程結構設計過程中存在固有的客觀不確定性和主觀不確定性,進行結構優化設計時需要慎重考慮,邱志平團隊近期在考慮不確定性線性結構拓撲優化領域做了大量工作[165-167].另外,文桂林團隊[44-49]近年來在穩健性拓撲優化設計方面取得了一些進展,成功將云模型引入到拓撲優化列式,得到了綜合主觀和客觀不確定性的創新性設計;將區間優化問題轉化為等效的確定性多工況優化設計問題,大大提升了計算效率,如圖10(a)～10(c)所示.而考慮不確定性的非線性拓撲優化設計是一個更具挑戰和值得探索的領域.

圖10 穩健性拓撲優化相關研究: (a) 外載力方向不確定性穩健性優化設計[44];(b) 區間不確定性穩健性優化設計[46];(c) 基于虛擬隨機載荷的穩健性拓撲優化設計[48];(d) 外載力位置不確定性穩健性優化設計[49]Fig.10 Research on robust topology optimization: (a) robust topology optimization (RTO) with uncertain direction of external load[46] and compared with that obtained from deterministic topology optimization(DTO);(b) interval uncertainty robust topology optimization[46];(c) robust topology optimization based on pseudo random load[48];(d) robust topology optimization with uncertainty ofexternal load position[49]

4 結論和展望

航空航天、航海工程、機器人等裝備領域結構不管從自身材料屬性還是使用要求來看,都存在大量的非線性問題,即,材料非線性、幾何非線性和邊界非線性.而連續體拓撲優化方法為諸如這些結構的設計提供了一種優于傳統經驗方法的創新手段,因此,研究非線性拓撲優化設計方法具有重要意義.本文詳細介紹了國內外在非線性連續體拓撲優化方面的研究進展,可以發現,雖然非線性拓撲優化方法近年來取得了長足的進展,但還尚未成熟,亟需進一步的發展.

(1)材料非線性連續體拓撲優化和幾何非線性連續體拓撲優化相對于邊界非線性連續體拓撲優化研究較多,原因在于接觸問題的強非線性和不連續性對優化模型的建立以及優化算法提出了更高的挑戰,特別對于高速碰撞結構和材料的拓撲優化設計問題,需要進一步的研究.材料非線性拓撲優化研究大部分局限于理想彈塑性假設,需要發展適用于更加復雜應力-應變行為的方法,特別是超彈材料結構設計研究.幾何非線性拓撲優化大多局限于大變形小應變問題,對于大變形大應變,甚至高應變率問題需要進一步探索;此外,目前方法大多數在拉格朗日法描述框架下進行,非線性拓撲優化循環中有限元響應計算存在精度差的問題,因此,基于歐拉法的幾何非線性拓撲優化方法值得進一步的探索.

(2)目前非線性連續體拓撲優化還局限于靜力學或準靜態系統,而實際工程結構,如衛星、火箭等在發射或者使用過程中,常常表現為復雜的非線性動力學行為[168-169],考慮復雜非線性動力學行為的連續體拓撲優化方法還鮮有人研究,是將來研究的重點方向.此外,非線性連續體拓撲優化由于需要求解大量的非線性方程組,計算效率往往很低,限制了其工程應用,基于多分辨率的思路可以實現高效率和保精度的拓撲優化設計[90,170-173],但該方法仍存在諸如在過濾半徑較小的情況容易出現棋盤格現象和QR 模式等問題[174],需要進一步解決.機器學習、并行高性能計算和多重網格法也是克服非線性連續體拓撲優化計算量大的潛在研究方向[6,8,175-178].