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

1 大連理工大學(xué), 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室, 工程力學(xué)系, 遼寧大連 116023

2 大連理工大學(xué), 國際計(jì)算力學(xué)研究中心, 遼寧大連 116023

3 大連理工大學(xué)寧波研究院, 浙江寧波 315016

4 中國運(yùn)載火箭技術(shù)研究院, 北京 100076

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

1 引 言

拓?fù)鋬?yōu)化可以通過在設(shè)計(jì)域內(nèi)合理布置材料達(dá)到優(yōu)化結(jié)構(gòu)性能的目的, 該方法已被廣泛應(yīng)用于重大裝備和先進(jìn)材料/結(jié)構(gòu)的構(gòu)型創(chuàng)新設(shè)計(jì). 目前較為流行的拓?fù)鋬?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ù)描述結(jié)構(gòu)邊界的水平集方法(Wang et al. 2003)以及具有顯式幾何信息的MMC/MMV (moving morphable components/moving morphable void)方法(Guo et al. 2014, Zhang et al. 2017). 其中應(yīng)用最為廣泛的SIMP方法將密度場用有限元網(wǎng)格離散,為了得到高分辨率設(shè)計(jì)和較為光滑的邊界需要充分加密網(wǎng)格, 其弊端是大大增加了有限元分析的計(jì)算量. 另外, 對于三維結(jié)構(gòu)拓?fù)鋬?yōu)化問題, 細(xì)化網(wǎng)格帶來的“維度煩惱”是不可忽視的.

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

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

2 超單元?jiǎng)偠染仃嚨撵o態(tài)凝聚

圖1為子結(jié)構(gòu)的靜態(tài)凝聚示意圖(Guyan 1965), 藍(lán)色為被凝聚的節(jié)點(diǎn), 用下標(biāo)1表示; 紅色節(jié)點(diǎn)為保留的節(jié)點(diǎn), 用下標(biāo)2表示.

圖1

子結(jié)構(gòu)的平衡方程可分塊表示為

由式(1)可得

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

定義 子結(jié)構(gòu)的平衡方程可以縮聚為保留節(jié)點(diǎn)的平衡方程

3 超單元模板庫

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

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

圖2

4 靈敏度分析

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

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

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

5 數(shù)值算例

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

圖3

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

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

表 1 結(jié)果對比

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

表 2 密度閾值對比

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

圖4

6 結(jié) 論

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

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