基于多分辨率-多邊形單元建模策略的多材料結構動剛度拓撲優化方法

江旭東，馬佳琪，熊 志，滕曉艷，王亞萍

(1.哈爾濱理工大學機械動力工程學院，哈爾濱 150080；2.哈爾濱工程大學機電工程學院，哈爾濱 150001)

在飛機、汽車制造與結構工程領域，多材料設計兼顧每一相基材料的有點，能夠保證結構具有綜合的力學性能，因而在飛機、汽車與橋梁結構中的應用日益增多[1-3]。拓撲優化方法能夠尋求滿足性能約束條件下材料的最優布局，極大地豐富了多材料結構的設計手段。與傳統單一材料的拓撲優化相比，多相材料的拓撲優化極大地擴大了設計空間，能夠提供綜合性能更優的最優解。隨著多材料增材制造技術的發展，依據最優拓撲的計算模型，使通過逐點和逐層方式打印多材料設計結構成為可能[4-6]。由此，多材料的拓撲優化問題引起了研究機構的廣泛關注。根據設計變量的定義方式，多材料拓撲優化方法通常歸納為2 類：基于單元的方法[7-18]和基于邊界的方法[19-23]。在基于單元的方法框架內，ZUO 等[7]、杜義賢等[8]、閆浩和吳曉明[9]通過構造材料性能的序列插值模型，實現了柔性機構和傳熱機構的多材料最優布局設計。俞燎宏等[10]和朱本亮等[11]采用交替主動相算法，將多相優化布局問題轉變為序列兩相材料優化子問題，然后通過并行計算策略逐級分層計算。目前，多體積約束下的多材料布局優化已拓展至穩態熱傳導[12]、熱力耦合[13-14]、頻域動力學[15]和空間桁架結構[16]的最優布局問題。最近，YANG和LI[17]、HUANG 和LI[18]則提出了單一質量約束下的多材料輕量化設計方法。在基于邊界方法的框架內，水平集法[19-21]和移動構件法(MMC)[22-23]因其能使邊界光滑清晰、便于提取設計構型等獨特優勢，也受到了一定的關注與研究。但是，現有的基于邊界的方法未能充分考慮新孔洞的創建，優化結果過度依賴于初始設計，難于獲得優化問題的全局最優解。

拓撲優化框架要求在消耗盡可能小的計算成本下獲得高分辨率的優化結果，而上述目標的實現取決于有限元求解器、自由度數量、優化建模策略、材料插值模型以及后處理等諸多因素。SUKUMAR[24]首次提出了多邊形有限元法，數值實驗表明它非常適于求解具有復雜幾何結構區域的連續介質力學問題。將多邊形單元應用于拓撲優化問題，能夠顯著地減小棋盤格和孤島效應等數值奇異性問題[25-28]。為了減少計算成本與提高求解精度，NGUYEN-XUAN[29]構建了多邊形有限元法與自適應網格重劃技術融合策略，CHAU 等[30]基于上述研究工作求解了多材料結構的靜態拓撲優化問題。另外，NGUYEN 等[31]提出了高分辨率拓撲優化方法 (multiresolution topology optimization，MTOP)，它采用多層級網格優化建模策略，即利用粗糙網格完成有限元分析，精細的重疊網格描述設計變量和密度變量空間，從而形成高分辨率的拓撲優化結果。最近，FILIPOV 等[32]將多邊形單元替代傳統單元作為分析單元，精確地估計復雜結構的動力學響應，高效地實現了特征值、受迫振動等結構頻域動力學問題的高分辨率拓撲優化。

在結構時域動響應方面，以提高結構動剛度[27,33]，降低結構振動幅度[34-35]和動應力幅值[28]等拓撲優化方法相繼提出。SHOBEIRI[33]在敏度分析中忽略了動能對單元刪除的影響，難于獲得慣性力影響顯著的動力學拓撲優化解。ZHAO 和WANG[34]、GIRALDO-LONDO?O 等[27-28]分別采用了先微分-后離散和先離散-后微分的伴隨敏度分析方法，但是，后者基于多邊形單元計算結構動響應，能夠處理任意曲線邊界結構的動力學優化問題。針對模型減縮技術在動力學拓撲優化過程中的有效性問題，ZHAO 和WANG[35]指出模態加速法在減少響應計算規模和預測精度上的平衡能力優于模態位移法，然而，如果高頻模態為主要影響模態時，模型減縮技術將失效。

多相材料時域動力學布局優化問題屬于強非自伴隨問題，需要考慮多材料結構的動態約束條件，致使敏度分析和優化問題求解具有挑戰性，因而目前的多相材料布局優化集中于靜態優化問題，而多相材料的時域動力學布局優化問題研究相對較少。因此，本文將多邊形網格與MTOP 的融合方法拓展至多材料結構的動態拓撲優化問題，建立融合框架內的多相材料高分辨率布局的時域動力學優化模型。為了避免先微分-后離散方法引起的靈敏度計算的相容性誤差，采用先離散-后微分的伴隨敏度分析方法。通過泰勒公式實現近似描述目標函數和體積約束函數，建立原問題的凸規劃子模型，采用基于敏度分離技術的ZPR方法求解子問題。最后，通過數值算例檢驗提出的方法在多相材料動力學布局優化的可行性。

1 多分辨率-多邊形網格建模策略

多邊形單元在動力學分析方面其求解精度優于傳統單元，因而將其作為動力學優化問題的響應計算單元[25-28]。同時，將多邊形單元劈分為精細的小單元形成設計變量與密度變量網格，并使兩者具有相同的位置坐標，繼而構建多分辨率-多邊形單元的優化建模策略，能夠實現粗糙位移網格條件下的高分辨率的拓撲優化構型設計(如圖1所示)。多分辨率-多邊形單元建模策略，利用粗糙的多邊形網格精確求解位移場，精細的設計變量與密度變量重疊網格用于構型優化設計，將有效地提升優化計算效率與優化結果的品質[26]。

在MTOP 框架內，為了精確計算多邊形位移單元的單元剛度矩陣，將形函數及其梯度的積分點設置在密度變量所在位置(也是設計變量網格的中心)。對于多材料問題，單元剛度和質量矩陣分別表示為：

式中：Nn為多邊形位移單元積分點的個數；y?,ij為密度變量在積分點處的估計值；B?為應變矩陣，D0為線彈性材料本構矩陣；A?,i為設計變量網格或密度變量網格的面積； ρ0為材料密度；m?E(y?,ij)和m?V(y?,i j)為m種材料的剛度與體積差值函數。

基于門檻投影函數[36]，體積插值函數表示為：

式中，p0>0為懲罰參數。然后，構造多材料剛度插值函數，表示為[38]：

式中，Ei為第i種材料的彈性模量。

由此，根據上述單元剛度與質量矩陣，多邊形位移單元的總體矩陣表示為：

為了抑制棋盤格與孤島現象，利用線性濾波方法獲得網格無關的優化結果，則有：

式中：Si為相應于密度變量單元i所占的子域；xn為與設計變量d?,nj對應的設計變量網格的中心坐標。線性權重函數定義為：

式中：rni為單元密度單元i和n的中心距；rmin為指定的過濾半徑。通過式(8)將設計變量線性加權處理而映射為密度變量，以此代入式(1)和式(2)中即可獲得單元剛度與質量矩陣。

此外，根據式(8)，密度變量對于設計變量的靈敏度為：

2 多材料結構動剛度拓撲優化

2.1 優化模型

以平均動柔度最小化為目標和多材料的體積占比為約束，建立多材料結構的動力學拓撲優化模型。假設設計域中包含m種材料，按彈性模量高低線性排列，在有限材料約束下，其動剛度優化模型為：

式中：d為設計變量矢量；Nt為動態激勵作用時間定分的區間數；Nc為多材料體積約束的數量；fi為t=ti時的動載荷向量；ui、u˙i、u¨i分別為相應的結構位移、速度、加速度響應；C=αrM+βrK為阻尼矩陣( αr、 βr為瑞利阻尼參數)； εj、 ηj、χj分別為單元指標集、設計變量指標集和多材料相數指標集，為第j種材料占設計域總體積的體積分數。

2.2 HHT-α 方法

HHT-α 方法作為廣義的Newmark-β 方法，它通過數值阻尼參數α 調控動力學控制方程的積分格式，使迭代過程無條件穩定，非常適于求解結構動力學問題。根據文獻[39]，HHT-α 方法將優化模型式(11)中的半離散形式的有限元方程修改為：

通過Newmark-β 有限差分關系，位移、速度場的更新格式為：

式中，β=(1+α)2/4 ，γ= (1+2α)/2為算法參數，合理選擇參數α 保證算法具有至少二階精度和無條件穩定。

將式(13)、式(14)代入式(12)，得到離散形式的控制方程的殘差，則有：

對于HHT-α 殘差方程式(15)，其初始時刻的表達式為：

2.3 敏度分析

根據優化模型式(11)，在MTOP 框架內目標函數對于設計變量的靈敏度為：

為了避免狀態變量對于設計變量的導數計算，使用伴隨變量法完成目標函數的敏度分析。類似于HHT-α 殘差方程式(15)，將Newmark-β 有限差分關系式(13)、式(14)變換為殘差形式：

根據式(19)、式(20)，?Pi/?d?,n j=0 和?Qi/?d?,nj=0以及初始條件?u0/?d?,nj=0 和?u˙0/?d?,nj=0，同時，滿足伴隨方程式(23)和式(24)使?u¨0/?d?,nj、?u¨i/?d?,nj、?u˙i/?d?,nj和?ui/?d?,n j等項消失，則式(21)化簡為：

將HHT-α 殘差方程式(15)、初始條件式(17)和Newmark-β 有限差分關系式(13)、式(14)代入伴隨方程式(23)、式(24)，則有：

由于優化模型式(11)中使用了體積插值函數V=與剛度插值函數E=，那么通過鏈式求導法則，拉格朗日函數與約束函數的靈敏度為：

我們這一代人一般都是三口之家。三居室的分配是這樣的：夫妻一間，孩子一間，書一間。書跟人平起平坐，甚至因書太多，書房得足夠寬敞，孩子只能屈尊居于小房間。

另外，根據優化模型式(11)，目標函數和約束函數對于V=和E=的敏度表示為：

結合式(1)、式(2)和式(7)，HHT-α 殘差方程對于V=和E=的偏導數表示為：

由此，根據式(25)～式(28)確定伴隨變量，隨后將式(31)～式(37)代入式(29)和式(30)，則可確定拉格朗日函數與約束函數的靈敏度。

2.4 ZPR 設計變量更新方案

利用凸近似方法，優化模型式(11)的近似子問題：

式(42)表明，每一個設計變量僅與一個約束函數相關，因而拉格朗日函數L可采用分離變量的方法求極值。鑒于Li的可分離特性，其一階最優條件為：

將式(45)代入式(42)，獲得子問題的對偶目標函數，根據Li(d(λi),λi)的相互獨立性，尋求拉格朗日乘子 λi使其達到極大值，則Li(d(λi),λi)的駐點條件表示為：

式(46)為關于 λi的非線性方程，可通過區間收縮法或迭代法求解。

3 典型算例

本節通過懸臂梁、簡支梁、L 型結構與Michell型懸臂梁4 個典型數值算例，對比4 種工況：①靜態載荷；② 載荷持續時間為0.05 s；③ 載荷持續時間為0.03 s；④ 2 種～10 種不同彈性模量(Ei=E0(m-i+1)/m，i=1,2,···,m，m=2,3,···,10，如圖2 所示)多相材料的優化結果。分析該優化方法的可行性和動態載荷作用時間對優化結果的影響機制。

圖2 多相材料的彈性模量Fig.2 Elastic modulus of multiphase materials

對于任一種多材料優化工況，每一相材料體積約束的上限vˉj=0.45/m(j=1,2,···,m)。另外，為了確保HHT-α 方法至少具有二階精度和無條件穩定性，所用算例均使用α=0.05 ，β=(1+α)2/4，γ= (1+2α)/2。

3.1 懸臂梁

如圖3 所示，懸臂梁右側自由邊中點處承受正弦半波載荷f(t)=f0sin(πt/tf)，載荷幅值f0=1 kN，初始設計域的結構尺寸為L×H×h=8 m×4 m×0.01 m。多相材料的基準彈性模量E0=200 GPa，任一材料相的泊松比νi=0.3，密度ρi=7800 kg/m3，瑞利阻尼參數αr=10 ，βr=1×10-5。

圖3 懸臂梁Fig.3 Cantilever beam

圖4 為懸臂梁雙材料設計的迭代歷史(載荷作用時間tf=0.05 s)。目標函數漸進收斂于0.0294 N·m，軟、硬兩相材料均漸進收斂于體積約束的上限：

圖5 為不同載荷持續時間工況下雙材料懸臂梁問題的優化設計結果。靜力學的優化拓撲更接近于動力學tf=0.05 s的結果，但是，動力學tf=0.03 s的優化拓撲顯著不同于前兩者的結果。與靜力學解相比，動力學解分布硬的材料相或三角形結構用于強化載荷作用點處的剛度，以抵制動載荷的快速波動。如圖6 所示，靜力學工況的位移與動柔度響應均高于動載荷工況，從而驗證了前面的假設。另外，tf=0.05 s工況下的設計將硬材料相置于載荷作用端，使其產生較小的變形，因而具有更高的動剛度。

圖5 雙材料懸臂梁問題動靜態載荷優化結果對比Fig.5 Comparison of dynamic and static load optimization results of bimaterial cantilever beam

圖6 最優懸臂梁載荷作用點處垂向位移和動柔度時間歷程Fig.6 Time history of vertical displacement and dynamic flexibility at the loading point of optimal cantilever structure

圖7 為tf=0.05 s時懸臂梁多材料拓撲優化結果。外部拓撲結構基本一致，內部拓撲結構依靠交叉肋結構強化，隨著材料相數目的增加，低彈性模量材料向固定端分布，特別是對于4 種材料和10 種材料工況，固定端增加了強化結構平衡慣性力的影響。

圖7 懸臂梁的多材料拓撲優化Fig.7 Multi material topology optimization of cantilever beam

圖8 為懸臂梁多材料最優拓撲的動柔度時間歷程，多相材料設計的動柔度變化基本一致，平均動柔度幾乎相同(如圖7 所示)。由于7 相多材料設計的最軟相材料位于內部加強肋的交叉點處，使結構的平均動柔度高于其他多相材料設計。

圖8 懸臂梁多材料最優拓撲的動柔度時間歷程Fig.8 Time history of dynamic flexibility for optimized multi-material topology of cantilever beam

3.2 簡支梁

如圖9 所示，簡支梁上端自由邊中點處承受余弦半波載荷f(t)=f0cos(πt/tf)，載荷幅值f0=1 kN，初始設計域的結構尺寸為：L=12 m，H=2 m，h=0.01 m。多相材料的基準彈性模量E0=200 GPa，任一材料相的泊松比νi=0.3，密度ρi=7800 kg/m3，瑞利阻尼參數αr=10 ，βr=1×10-5。

圖9 簡支梁Fig.9 Simply supported beam

圖10 為不同載荷持續時間工況下雙材料簡支梁問題的優化設計結果。3 種工況的優化拓撲皆是由兩側支撐桿和中間環型結構組成，但是兩相材料分布迥然不同。靜力學工況的優化拓撲在中間環型結構處集中了高彈性模量材料，因而它的動剛度大于動力學工況優化拓撲的動剛度；但是，動力學工況的優化拓撲在兩側支撐桿集中了高彈性模量材料，能夠進一步抑制慣性載荷的影響，致使其動剛度大于靜力學工況優化拓撲的動剛度。如圖11 所示，靜力學工況的位移與動柔度響應介于兩種動載荷工況優化拓撲的相應動響應之間，驗證了前面的結論。另外，隨著動力學載荷作用時間的減少，阻尼對于結構的能量耗散作用有所減弱。

圖10 雙材料簡支梁問題動靜態載荷優化結果對比Fig.10 Comparison of dynamic and static load optimization results of bimaterial simply supported beam

圖11 最優簡支梁載荷作用處垂向位移和動柔度時間歷程Fig.11 Time history of vertical displacement and dynamic flexibility at the loading point of the optimal simply supported beam structure

圖12 為tf=0.05 s 時簡支梁的多材料拓撲優化結果，多相材料設計的平均動柔度基本一致。多相材料結構的整體拓撲基本一致，中間環型強化結構有所差異，特別是4 種材料、6 種材料和10 種材料的優化拓撲差異顯著。隨著材料相數目的增加，中間環型強化結構的下端桁架部分布置更多的低彈性模量材料，兩側支撐桿結構布置更多的高彈性模量材料，使兩者能夠達到剛度的最優匹配。

3.3 L 型結構

如圖13 所示，L 型結構右側自由邊中點處承受正弦半波載荷f(t)=f0sin(πt/tf)，載荷幅值f0=1 kN，初始設計域的結構尺寸為：L1=H2=12 m，L2=H1=18 m，厚度h=0.01 m。相材料的基準彈性模量E0=200 GPa，任一材料相的泊松比νi=0.3，密度ρi=7800 kg/m3，瑞利阻尼參數αr=50 ，βr=3×10-5。

圖13 L 型結構Fig.13 L-shaped structure

圖14 為不同載荷持續時間工況下雙材料L 型結構問題的優化設計結果。動力學工況的優化拓撲在內側直角應力集中處分布了更多的高彈性模量材料，而靜力學工況的優化拓撲在直角拐角處增加了加強筋強化結構。與靜力學優化拓撲相比，動力學最優拓撲在載荷作用點附近分布了更多的材料，特別是外側的三角形構型具有更大的截面尺寸，因而具有更大的動剛度。如圖15 所示，靜力學工況的位移與動柔度響應均高于動載荷工況，其動剛度偏低，從而驗證了前面的假設。另外，動力學tf=0.05 s工況下的優化拓撲在內側直角處堆積了更多的高彈性模量材料，使得其動剛度大于tf=0.03 s工況下優化拓撲的動剛度。

圖14 雙材料L 型結構問題動靜態載荷優化結果對比Fig.14 Comparison of dynamic and static load optimization results of bimaterial L-shaped structure

圖15 最優L 型結構載荷作用點處的垂向位移和動柔度時間歷程Fig.15 Time history of vertical displacement and dynamic flexibility at the loading point of optimal L-shaped structure

圖16 為tf=0.05 sL 型結構的多材料拓撲優化結果。多相材料結構的整體拓撲構型基本一致，載荷作用點附近的三角形構型的強化方式為：增加內側棱邊的截面寬度和內部添加加強筋。隨著材料相數目的增加，高彈性模量材料分布于固定端垂向平行桿附近，低彈性模量材料分布于載荷作用點附近的三角形構型。但是，由于上述兩種強化作用，致使載荷作用點附近的三角形構型與固定端垂向平行桿構型能夠達到剛度匹配，從而保證整體結構的動剛度最優。對于6 相和10 相多材料設計，最軟相材料位于載荷作用點附近的三角形構型的內部棱邊，致使結構的平均動柔度均高于其他多相材料設計。

圖16 L 型結構的多材料拓撲優化Fig.16 Multi-material topology optimization of L-shaped structure

3.4 Michell 型懸臂梁

如圖17 所示，Michell 型懸臂梁右側自由邊中點處承受正弦半波載荷f(t)=f0sin(πt/tf)，載荷幅值f0=1 kN，初始設計域的結構尺寸為：L=5 m，H=4 m，h=0.01 m，R=1 m。多相材料的基準彈性模量E0=200 GPa，任一材料相的泊松比νi=0.3，密度ρi=7800 kg/m3，瑞利阻尼參數αr=2.5 ，βr=4.5×10-4。

圖17 Michell 型懸臂梁Fig.17 Michell cantilever beam

圖18 為不同載荷持續時間工況下雙材料Michell 型懸臂梁問題的優化設計。與靜力學工況的優化設計相比，動力學工況的優化拓撲內部加強肋的具有更大的截面尺寸，載荷作用點附近聚集了更多的高彈性模量材料，因而后者具有更高的動剛度。此外，動力學tf=0.05 s工況下的優化拓撲內部加強肋的寬度尺寸較大，其方位更接近于垂向，有利于抑制動載荷作用下的變形，相對于tf=0.03 s工況下優化拓撲具有更高的動剛度。如圖19 所示，靜力學工況的位移與動柔度響應均高于動載荷工況，其動剛度偏低，而動力學tf=0.05 s工況的上述動響應最低，其動剛度最高，與優化拓撲的分析結果一致。

圖18 雙材料Michell 懸臂梁問題動靜態載荷優化結果對比Fig.18 Comparison of dynamic and static load optimization results of bimaterial Michell cantilever beam

圖19 最優Michell 型懸臂梁結構載荷作用點處的垂向位移和動柔度時間歷程Fig.19 Time history of vertical displacement and dynamic flexibility at the loading point of optimal Michell cantilever beam structure

圖20 為Michell 型懸臂梁的多材料拓撲優化結果。外環拓撲極為相似，內部強化構型差異顯著。高比例的硬材料分布于外部拓撲，而軟材料在內部聯接于內部拓撲。另外，載荷作用點附近也集中了高彈性模量材料，抵消慣性載荷的動態影響。對于9 相和10 項多材料設計，最軟相或次最軟相材料位于載荷作用點附近的三角形構型的內側棱邊，致使結構的平均動柔度均高于其他多相材料設計。

圖20 Michell 型懸臂梁的多材料拓撲優化Fig.20 Multi-material topology optimization of Michell cantilever beam

4 結論

本文基于多分辨率-多邊形單元建模策略，提出了求解多材料結構動響應問題的拓撲優化方法，研究結論如下：

(1) 通過懸臂梁、簡支梁、L 型結構和Michell型懸臂梁等典型數值算例，驗證了所提多材料結構時域動剛度拓撲優化方法的可行性。

(2) 動態載荷的作用時間對于優化結果具有重要的影響。隨著載荷作用時間的減少，結構內慣性力的影響更加顯著，致使多材料結構的優化拓撲存在明顯的差異。

(3) 將本文所提出的多材料結構動力學優化方法延伸至超材料結構的動力學設計問題，可以實現高分辨率胞元-宏觀結構構型的動力學多尺度協同優化設計。

基于提出的動力學多材料拓撲優化框架，作者將在后面工作中考慮動力學載荷作用下多材料之間的界面損傷與解耦，提出多材料結構界面失效問題的動力學拓撲優化方法。