黃孟成 霍文杜棟宗 亮 劉,*暢郭 ,2旭 楊東,生 黃 佳

1 大連理工大學, 工業(yè)裝備結構分析國家重點實驗室, 工程力學系, 遼寧大連 116023

2 大連理工大學, 國際計算力學研究中心, 遼寧大連 116023

3 大連理工大學寧波研究院, 浙江寧波 315016

4 中國運載火箭技術研究院, 北京 100076

5 北京強度環(huán)境研究所，北京 100076

1 引 言

拓撲優(yōu)化可以通過在設計域內合理布置材料達到優(yōu)化結構性能的目的, 該方法已被廣泛應用于重大裝備和先進材料/結構的構型創(chuàng)新設計. 目前較為流行的拓撲優(yōu)化方法包括基于密度場的SIMP (solid isotropic material with penalization)方 法(Sigmund 2001)和ESO (evolutionary structural optimization)/BESO (bi-directional evolutionary structural optimization)方 法(Xie et al. 1993, Querin et al. 1998), 用水平集函數(shù)描述結構邊界的水平集方法(Wang et al. 2003)以及具有顯式幾何信息的MMC/MMV (moving morphable components/moving morphable void)方法(Guo et al. 2014, Zhang et al. 2017). 其中應用最為廣泛的SIMP方法將密度場用有限元網(wǎng)格離散,為了得到高分辨率設計和較為光滑的邊界需要充分加密網(wǎng)格, 其弊端是大大增加了有限元分析的計算量. 另外, 對于三維結構拓撲優(yōu)化問題, 細化網(wǎng)格帶來的“維度煩惱”是不可忽視的.

因此, 如何解決網(wǎng)格細化導致的有限元分析時間激增的問題, 是近年來拓撲優(yōu)化領域重點關注的方向之一: Borrvall和Petersson (2001)利用并行技術優(yōu)化設計了百萬單元的三維模型;Bruns和Tortorelli (2003)通過刪除低密度單元, 減少了分析自由度以提高計算效率; de Sturler等(2008)通過自適應技術減少低密度區(qū)域的網(wǎng)格數(shù), 節(jié)省了有限元計算時間. 此外, Nguyen等(2010)提出了多分辨率拓撲優(yōu)化方法(multi-resolution topology optimization, MTOP). 該方法通過將有限元分析網(wǎng)格與離散密度場解耦, 在粗網(wǎng)格(超單元)下進行有限元分析, 而應用細網(wǎng)格進行結構拓撲描述及優(yōu)化, 可以顯著減少有限元計算量, 因而被應用于涉及多尺度、幾何非線性等有 限 元 計 算 量 更 為 突 出 的 結 構 拓 撲 優(yōu) 化 問 題(Park et al. 2021, Ortigosa and Martínez-Frutos 2021, Yoo et al. 2021, Wang et al. 2020). 由于超單元的單元剛度陣根據(jù)其內部的平均密度逐個計算, 無法保證其有限元分析結果的精度, 導致相應的優(yōu)化結果對過濾半徑以及過濾方式十分敏感. 特別的, 在過濾半徑小于超單元尺寸時, 可能會出現(xiàn)棋盤格現(xiàn)象和QR (quick response)模式,即超單元內部出現(xiàn)了類似二維碼圖案的不連續(xù)密度分布(Gupta et al. 2018), 從而導致其單剛被過高估計. 為了克服多分辨率方法的以上缺陷, Groen等(2017)及Nguyen等(2017)建議借助高階有限體積法和高階有限元法以提高結構分析精度和提升結果的分辨率. 但它們都需要更多的密度分布參與單剛計算, 不可避免地增加了計算成本. 近來, 劉暢等(Liu et al. 2018)發(fā)展了基于MMC的多分辨率方法, 利用組件的幾何分布特點, 在一定程度上避免了QR模式的出現(xiàn). 但由于其超單元剛度計算仍采用類似MTOP方法的密度平均格式, 本質上仍存在剛度高估的問題.

為改善基于SIMP的多分辨拓撲優(yōu)化方法在分析精度、棋盤格現(xiàn)象和QR模式方面的缺陷,本文將超單元視為子結構, 通過靜態(tài)凝聚得到其單元剛度陣以保證有限元分析精度, 并進一步根據(jù)拓撲優(yōu)化過程中子結構的密度分布特征構建子結構的模板庫, 從而省去了子結構單剛的重復計算, 進一步提高了計算效率.

2 超單元剛度矩陣的靜態(tài)凝聚

圖1為子結構的靜態(tài)凝聚示意圖(Guyan 1965), 藍色為被凝聚的節(jié)點, 用下標1表示; 紅色節(jié)點為保留的節(jié)點, 用下標2表示.

圖1

子結構的平衡方程可分塊表示為

由式(1)可得

若被凝聚節(jié)點上無外力, 即f1=0, 式(2)代入式(1)可得

定義 子結構的平衡方程可以縮聚為保留節(jié)點的平衡方程

3 超單元模板庫

本文通過調用超單元的模板庫來替代子結構凝聚過程. 若直接將超單元內的密度值二值化(非0即1), 其對應模板庫數(shù)目會過于龐大(以圖2(a)所示含有5 × 5密度網(wǎng)格的超單元為例, 共計 225種模板). 為縮減模板庫的規(guī)模, 結合數(shù)值經驗, 將圖2(b)所示的子結構4個角點單元保持原有密度、區(qū)域2, 4, 5, 6, 8內的單元密度分別取平均, 而后將這9個獨立密度根據(jù)閾值二值化,從而可減縮至512種模板.

為更高效地檢索到凝聚剛度陣, 可以將二值化的9個密度值按圖2(b)中的順序形成二進制數(shù) , 并將其對應的十進制數(shù)作為序號存儲在模板庫以方便提取.

圖2

4 靈敏度分析

在使用模板庫進行有限元計算時, 對超單元密度場按照圖2(b)所示進行平均和二值化處理,可能會導致參與有限元計算的密度分布與真實的密度場分布存在差異. 因此須對目標函數(shù)的靈敏度做如下調整

式中C為柔順度,p為懲罰因子分別為實體單元的剛度陣以及單元節(jié)點對應的位移列向量.由于僅作為中間變量參與有限元分析, 體積分數(shù)的計算以及變量更新均對應于因而無需修改體積靈敏度及收斂準則.

進一步, 為了使結構在迭代過程中更加貼合模板庫中的模板(即減少中間密度單元)和讓收斂過程更加平穩(wěn), 采用逐步減小體積的策略, 即

5 數(shù)值算例

為驗證本文所提方法的有效性, 這里采用SIMP方法和本文方法分別求解了圖3所示懸臂梁的最小柔順度優(yōu)化問題. 不失一般性, 所有幾何、材料以及載荷參數(shù)均采用無量綱形式, 其中設計域尺寸為 6×3, 上表面受均布力q=1, 體積分數(shù)上限材料彈性模量E=1,Emin=0.001, 泊松比ν=0.3 ; 懲罰因子p=3, 過濾半徑為1.5倍細網(wǎng)格寬度; 模板的密度閾值為0.45, 體積

圖3

θ=0.007縮 減 比 率 . 計 算 環(huán) 境 為Dell-7730工 作 站, Intel(R) Xeon(R) E-2186M 2.90 GHz CPU,32 GB RAM, Windows 10 OS, MATLAB 2018a.

采用四節(jié)點雙線性單元將結構離散為 600×300個 超單元, 每個超單元包含 5×5的細網(wǎng)格, 并將上表面的兩層密度單元固定為實體. 表1分別給出了本文方法和 3000×1500網(wǎng)格對應的SIMP拓撲優(yōu)化結果. 值得注意的是, 模板法所得最優(yōu)結構利用模板庫和精確分析的目標函數(shù)分別為10.69和10.71(相對誤差僅為0.187%), 相比SIMP法優(yōu)化結果(目標函數(shù)為10.54)的目標函數(shù)只高出了1.61%; 值得說明的是, 雖然模板法與SIMP方法得到優(yōu)化設計的結構柔度值較為接近, 但結構構型存在明顯差異, 這本質上是由于結構拓撲優(yōu)化問題的非凸性而存在眾多分布差異較大但性能相近的局部最優(yōu)解. 此外本算例采用的收斂條件為兩種算法都在最大步數(shù)跳出了循環(huán), 因此對比了前20步平均單次循環(huán)時間, 模板法(98.95 s)相比SIMP方法(115.67 s)單次迭代時間縮短了14.45%, 這說明了模板法具有較高的有限元分析精度和效率.

表 1 結果對比

為了說明密度閾值對優(yōu)化結果的影響, 還分別計算了閾值為0.1, 0.45和0.8的優(yōu)化結果. 如表2所示, 不同的密度閾值會對優(yōu)化結果產生一定的影響, 密度閾值較大時結構細節(jié)較多, 但結構輪廓和目標函數(shù)較為接近.

表 2 密度閾值對比

進一步, 按照模板法的相同參數(shù)設置, 應用MTOP方法和MMA算法(Svanberg 1987)求解了本問題, 所得結果如圖4. 圖4(a)所示優(yōu)化結果利用超單元和精確分析的目標函數(shù)分別為10.71和82.17, 相對誤差高達86.9%; 并且由于過濾半徑小于超單元的尺寸, 出現(xiàn)了如圖4(b)所示的棋盤格現(xiàn)象和QR模式. 這是由于MTOP方法在計算超單元單剛時將其內部密度場進行平均而無法識別超單元內密度場的具體分布導致的; 雖然將過濾半徑增大至超單元尺寸可以有效避免這些問題, 但這也將導致無法獲得小于超單元特征尺寸的結構細節(jié), 從而浪費了高分辨率的密度場離散, MTOP方法中過濾半徑對QR模式影響的詳細討論請參見Gupta 等 (2018)的研究.而模板法的超單元剛度陣是利用子結構靜態(tài)凝聚形成的, 考慮超單元內密度場的分布, 從而有效解決了MTOP方法中超單元剛度被高估的問題, 進而從根本上解決了QR模式、打破了多分辨率優(yōu)化方法過濾半徑的相關限制.

圖4

6 結 論

為節(jié)省結構拓撲優(yōu)化中的有限元分析成本并保證計算精度, 本文在MTOP框架的基礎上,基于子結構靜態(tài)凝聚得到超單元的剛度陣, 并根據(jù)拓撲優(yōu)化過程中子結構的密度分布特征構建了子結構的模板庫. 可以實現(xiàn)在整體剛度陣的維數(shù)減少的同時省去超單元單剛的重復計算、提升了有限元分析效率, 并有效抑制了MTOP方法中的QR模式. 通過求解受均布載荷的懸臂梁最小柔順度問題并分別與SIMP法和MTOP法的結果進行對比, 驗證了基于模板的子結構多分辨率拓撲優(yōu)化的有效性. 另一方面, 由于算例在于展示模板法在采用超單元進行高效結構分析的同時仍具備獲得結構細節(jié)的能力, 因而并未考慮過多細節(jié)對結構穩(wěn)定性的影響. 事實上, 將模板法推廣至涉及特征值求解的結構屈曲性能優(yōu)化設計問題也具有重要意義, 這將是后續(xù)工作的重要方向之一. 然而因為SIMP方法優(yōu)化過程中不可避免地存在中間密度單元, 本方法對密度場進行的二值化處理帶來了誤差, 為改善這一點可進一步結合0?1離散變量拓撲優(yōu)化方法(Liang et al.2019).

致 謝 本工作得到了國家自然科學基金(11821202, 11732004, 12002073, 12002077), 國家重點研發(fā)計劃(2020YFB1709401, 2016YFB0201601), 大連理工大學科研啟動項目(DUT20RC(3)020)和博士后科學基金(2020T130078, 2020M680944)的支持.