基于模板的子结构多分辨率拓扑优化

黄孟成 霍文杜栋宗 亮 刘,*畅郭 ,2旭 杨东,生 黄 佳

1 大连理工大学, 工业装备结构分析国家重点实验室, 工程力学系, 辽宁大连 116023

2 大连理工大学, 国际计算力学研究中心, 辽宁大连 116023

3 大连理工大学宁波研究院, 浙江宁波 315016

4 中国运载火箭技术研究院, 北京 100076

5 北京强度环境研究所，北京 100076

1 引 言

拓扑优化可以通过在设计域内合理布置材料达到优化结构性能的目的, 该方法已被广泛应用于重大装备和先进材料/结构的构型创新设计. 目前较为流行的拓扑优化方法包括基于密度场的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), 用水平集函数描述结构边界的水平集方法(Wang et al. 2003)以及具有显式几何信息的MMC/MMV (moving morphable components/moving morphable void)方法(Guo et al. 2014, Zhang et al. 2017). 其中应用最为广泛的SIMP方法将密度场用有限元网格离散,为了得到高分辨率设计和较为光滑的边界需要充分加密网格, 其弊端是大大增加了有限元分析的计算量. 另外, 对于三维结构拓扑优化问题, 细化网格带来的“维度烦恼”是不可忽视的.

因此, 如何解决网格细化导致的有限元分析时间激增的问题, 是近年来拓扑优化领域重点关注的方向之一: Borrvall和Petersson (2001)利用并行技术优化设计了百万单元的三维模型;Bruns和Tortorelli (2003)通过删除低密度单元, 减少了分析自由度以提高计算效率; de Sturler等(2008)通过自适应技术减少低密度区域的网格数, 节省了有限元计算时间. 此外, Nguyen等(2010)提出了多分辨率拓扑优化方法(multi-resolution topology optimization, MTOP). 该方法通过将有限元分析网格与离散密度场解耦, 在粗网格(超单元)下进行有限元分析, 而应用细网格进行结构拓扑描述及优化, 可以显著减少有限元计算量, 因而被应用于涉及多尺度、几何非线性等有 限 元 计 算 量 更 为 突 出 的 结 构 拓 扑 优 化 问 题(Park et al. 2021, Ortigosa and Martínez-Frutos 2021, Yoo et al. 2021, Wang et al. 2020). 由于超单元的单元刚度阵根据其内部的平均密度逐个计算, 无法保证其有限元分析结果的精度, 导致相应的优化结果对过滤半径以及过滤方式十分敏感. 特别的, 在过滤半径小于超单元尺寸时, 可能会出现棋盘格现象和QR (quick response)模式,即超单元内部出现了类似二维码图案的不连续密度分布(Gupta et al. 2018), 从而导致其单刚被过高估计. 为了克服多分辨率方法的以上缺陷, Groen等(2017)及Nguyen等(2017)建议借助高阶有限体积法和高阶有限元法以提高结构分析精度和提升结果的分辨率. 但它们都需要更多的密度分布参与单刚计算, 不可避免地增加了计算成本. 近来, 刘畅等(Liu et al. 2018)发展了基于MMC的多分辨率方法, 利用组件的几何分布特点, 在一定程度上避免了QR模式的出现. 但由于其超单元刚度计算仍采用类似MTOP方法的密度平均格式, 本质上仍存在刚度高估的问题.

为改善基于SIMP的多分辨拓扑优化方法在分析精度、棋盘格现象和QR模式方面的缺陷,本文将超单元视为子结构, 通过静态凝聚得到其单元刚度阵以保证有限元分析精度, 并进一步根据拓扑优化过程中子结构的密度分布特征构建子结构的模板库, 从而省去了子结构单刚的重复计算, 进一步提高了计算效率.

2 超单元刚度矩阵的静态凝聚

图1为子结构的静态凝聚示意图(Guyan 1965), 蓝色为被凝聚的节点, 用下标1表示; 红色节点为保留的节点, 用下标2表示.

图1

子结构的平衡方程可分块表示为

由式(1)可得

若被凝聚节点上无外力, 即f1=0, 式(2)代入式(1)可得

定义 子结构的平衡方程可以缩聚为保留节点的平衡方程

3 超单元模板库

本文通过调用超单元的模板库来替代子结构凝聚过程. 若直接将超单元内的密度值二值化(非0即1), 其对应模板库数目会过于庞大(以图2(a)所示含有5 × 5密度网格的超单元为例, 共计 225种模板). 为缩减模板库的规模, 结合数值经验, 将图2(b)所示的子结构4个角点单元保持原有密度、区域2, 4, 5, 6, 8内的单元密度分别取平均, 而后将这9个独立密度根据阈值二值化,从而可减缩至512种模板.

为更高效地检索到凝聚刚度阵, 可以将二值化的9个密度值按图2(b)中的顺序形成二进制数 , 并将其对应的十进制数作为序号存储在模板库以方便提取.

图2

4 灵敏度分析

在使用模板库进行有限元计算时, 对超单元密度场按照图2(b)所示进行平均和二值化处理,可能会导致参与有限元计算的密度分布与真实的密度场分布存在差异. 因此须对目标函数的灵敏度做如下调整

式中C为柔顺度,p为惩罚因子分别为实体单元的刚度阵以及单元节点对应的位移列向量.由于仅作为中间变量参与有限元分析, 体积分数的计算以及变量更新均对应于因而无需修改体积灵敏度及收敛准则.

进一步, 为了使结构在迭代过程中更加贴合模板库中的模板(即减少中间密度单元)和让收敛过程更加平稳, 采用逐步减小体积的策略, 即

5 数值算例

为验证本文所提方法的有效性, 这里采用SIMP方法和本文方法分别求解了图3所示悬臂梁的最小柔顺度优化问题. 不失一般性, 所有几何、材料以及载荷参数均采用无量纲形式, 其中设计域尺寸为 6×3, 上表面受均布力q=1, 体积分数上限材料弹性模量E=1,Emin=0.001, 泊松比ν=0.3 ; 惩罚因子p=3, 过滤半径为1.5倍细网格宽度; 模板的密度阈值为0.45, 体积

图3

θ=0.007缩 减 比 率 . 计 算 环 境 为Dell-7730工 作 站, Intel(R) Xeon(R) E-2186M 2.90 GHz CPU,32 GB RAM, Windows 10 OS, MATLAB 2018a.

采用四节点双线性单元将结构离散为 600×300个 超单元, 每个超单元包含 5×5的细网格, 并将上表面的两层密度单元固定为实体. 表1分别给出了本文方法和 3000×1500网格对应的SIMP拓扑优化结果. 值得注意的是, 模板法所得最优结构利用模板库和精确分析的目标函数分别为10.69和10.71(相对误差仅为0.187%), 相比SIMP法优化结果(目标函数为10.54)的目标函数只高出了1.61%; 值得说明的是, 虽然模板法与SIMP方法得到优化设计的结构柔度值较为接近, 但结构构型存在明显差异, 这本质上是由于结构拓扑优化问题的非凸性而存在众多分布差异较大但性能相近的局部最优解. 此外本算例采用的收敛条件为两种算法都在最大步数跳出了循环, 因此对比了前20步平均单次循环时间, 模板法(98.95 s)相比SIMP方法(115.67 s)单次迭代时间缩短了14.45%, 这说明了模板法具有较高的有限元分析精度和效率.

表 1 结果对比

为了说明密度阈值对优化结果的影响, 还分别计算了阈值为0.1, 0.45和0.8的优化结果. 如表2所示, 不同的密度阈值会对优化结果产生一定的影响, 密度阈值较大时结构细节较多, 但结构轮廓和目标函数较为接近.

表 2 密度阈值对比

进一步, 按照模板法的相同参数设置, 应用MTOP方法和MMA算法(Svanberg 1987)求解了本问题, 所得结果如图4. 图4(a)所示优化结果利用超单元和精确分析的目标函数分别为10.71和82.17, 相对误差高达86.9%; 并且由于过滤半径小于超单元的尺寸, 出现了如图4(b)所示的棋盘格现象和QR模式. 这是由于MTOP方法在计算超单元单刚时将其内部密度场进行平均而无法识别超单元内密度场的具体分布导致的; 虽然将过滤半径增大至超单元尺寸可以有效避免这些问题, 但这也将导致无法获得小于超单元特征尺寸的结构细节, 从而浪费了高分辨率的密度场离散, MTOP方法中过滤半径对QR模式影响的详细讨论请参见Gupta 等 (2018)的研究.而模板法的超单元刚度阵是利用子结构静态凝聚形成的, 考虑超单元内密度场的分布, 从而有效解决了MTOP方法中超单元刚度被高估的问题, 进而从根本上解决了QR模式、打破了多分辨率优化方法过滤半径的相关限制.

图4

6 结 论

为节省结构拓扑优化中的有限元分析成本并保证计算精度, 本文在MTOP框架的基础上,基于子结构静态凝聚得到超单元的刚度阵, 并根据拓扑优化过程中子结构的密度分布特征构建了子结构的模板库. 可以实现在整体刚度阵的维数减少的同时省去超单元单刚的重复计算、提升了有限元分析效率, 并有效抑制了MTOP方法中的QR模式. 通过求解受均布载荷的悬臂梁最小柔顺度问题并分别与SIMP法和MTOP法的结果进行对比, 验证了基于模板的子结构多分辨率拓扑优化的有效性. 另一方面, 由于算例在于展示模板法在采用超单元进行高效结构分析的同时仍具备获得结构细节的能力, 因而并未考虑过多细节对结构稳定性的影响. 事实上, 将模板法推广至涉及特征值求解的结构屈曲性能优化设计问题也具有重要意义, 这将是后续工作的重要方向之一. 然而因为SIMP方法优化过程中不可避免地存在中间密度单元, 本方法对密度场进行的二值化处理带来了误差, 为改善这一点可进一步结合0−1离散变量拓扑优化方法(Liang et al.2019).

致 谢 本工作得到了国家自然科学基金(11821202, 11732004, 12002073, 12002077), 国家重点研发计划(2020YFB1709401, 2016YFB0201601), 大连理工大学科研启动项目(DUT20RC(3)020)和博士后科学基金(2020T130078, 2020M680944)的支持.