一种高效的子单元拓扑优化方法：三维晶格的MATLAB实现

《Advances in Engineering Software》：An efficient sub-cell topology optimization method: MATLAB implementation for 3D lattices

　　赵力|王启涵|张明辉|高伟|罗振•开发用于三维晶格拓扑优化的MATLAB代码。•通过子单元设计变量提高效率。•为子单元公式推导灵敏度分析方法。•通过可重复的数值示例验证该方法。引言在轻量化结构设计[1]、[2]、[3]、[4]、能量吸收[5]、[6]以及生物医学植入物[7]领域中，设计具有特定性能的先进材料与结构至关重要。拓扑优化技术以及增材制造技术的发展为这类结构的設计提供了更高的灵活性和可制造性。近年来，人们提出了多种拓扑优化方法，包括带有惩罚项的各向同性固体材料法（SIMP）[8]、进化法[9]、水平集法[10]、[11]以及显式法[12]、[13]。在拓扑优化方法的发展过程中，配有开源教育代码的文献起到了重要作用，帮助读者更好地理解这些方法，并能将其快速应用于自身的设计问题中[14]。在基于SIMP的框架下，Sigmund[15]提出了一个包含99行代码的MATLAB程序，用于教学目的实现拓扑优化。后来该程序被改进为88行版本，显著提升了计算性能[16]。为了实现三维设计，Liu和Tovar[17]提出了一个包含169行代码的MATLAB程序，用于三维结构优化。尽管这些教育代码在结构优化方面非常有效，但在处理大规模、高分辨率的三维结构时，计算效率却成为一个重大挑战。尤其是，在解决合规性最小化问题时，求解线性方程组仍然是最耗时的步骤。为了解决这一问题，Amir等人[18]提出了一种集成多重网格预处理共轭梯度法（MGCG）的MATLAB实现方式，以减少计算时间。在此基础上，Aage等人[19]开发了一个完全并行的C++框架，进一步提升了计算效率。除了基于SIMP的方法外，还有其他几种拓扑优化框架也有了相应的教育代码。例如，双向进化结构优化法（BESO）[20]、[21]、水平集法[22]，以及移动变形组件法（MMC）[23]、[24]、[25]，都已有基于MATLAB或Python的实现代码，用以展示各自的方法论。除了在结构设计领域的应用外，拓扑优化技术还被拓展到晶格材料设计领域，尤其是用于设计具有所需物理特性的微观结构[26]、[27]。均匀化方法是计算周期性材料有效特性的有力工具。均匀化方法基本上分为两类：基于渐近展开的[28]、[29]以及基于能量的[30]、[31]均匀化方法。Andreassen和Andreasen[32]提出了一个用于二维问题数值均匀化计算的MATLAB代码，后来该代码又被扩展到三维问题[33]。随后，Ozdilek等人[34]在Dong等人[33]提出的MATLAB代码基础上开发了Python代码。Xia和Breitkopf[35]则基于能量均匀化方法开发了用于晶格材料设计的MATLAB代码。Gao等人[36]进一步将能量均匀化方法整合到并行多尺度拓扑优化框架中，用于设计复合结构。尽管现有的基于MATLAB的均匀化及材料优化代码已经提供了很有价值的教育工具，但其中许多代码主要侧重于有效性能评估，或是基于标准低阶离散化的二维材料设计。要将这类方法扩展到高分辨率的三维超材料设计仍然面临挑战，因为更大的设计空间通常需要大量的有限元，从而导致计算成本大幅上升。Gao等人[37]以及Nguyen等人[38]等人的先进研究已经证明了通过等几何分析进行三维正交各向异性超材料设计的可行性。然而，这些实现依赖于较为复杂的数值框架，初学者很难复现或修改。因此，目前仍缺乏一种基于更简单的基于密度法的易用教育框架，这样的框架能够在保持简洁且可重复的MATLAB实现的同时，提升三维晶格材料优化的设计精度。对于实际的三维晶格材料设计而言，解决大规模问题的效率至关重要，因为计算成本往往成为限制因素。在基于均匀化的晶格材料设计中，最耗时的步骤就是求解线性方程组。Kazakis和Lagaros[39]、[40]采用了模型降阶技术（MOR），以降低拓扑优化中全尺寸有限元分析的计算成本。另一种常见策略则是提升线性求解器的效率，比如使用带有多重网格或不完全乔列斯基预处理器的PCG方法[18]，或是针对逆均匀化问题使用基于GPU的多重网格求解器[41]。虽然这些技术有助于提升计算效率，但对于复杂的三维超材料结构来说，往往还需要进一步提升设计精度和建模灵活性。在这方面，高阶有限元为传统低阶有限元提供了一种替代方案，后者通常需要大幅细化网格才能达到类似的精度。在拓扑优化中，高阶有限元允许在单个元素内通过为每个积分点分配独立的材料参数，来实现材料属性的空间变化。为进一步提升设计精度，在多分辨率拓扑优化方法（MTOP）中，每个元素还可以被细分为多个子单元，每个子单元在多分辨率拓扑优化方法中都有独立的设計变量[42]、[43]、[44]。这种子单元策略也可以被引入基于有限单元法（FCM）[45]以及嵌入域法（EDM）[46]的拓扑优化框架中。除了常规有限元之外，FCM和EDM通常还会与p阶有限元方法[47]、[48]或等几何分析[49]相结合，以便处理具有不规则几何形状的复杂模型。这些研究表明，子单元离散化和高阶有限元分析是提升设计精度与分析成本之间平衡度的有效策略。然而，将这两种方法结合到基于MATLAB的紧凑型三维晶格材料拓扑优化框架中，目前仍较为有限。为填补这一空白，本研究开发了一种用于三维晶格材料设计的子单元拓扑优化框架（SCTOP）。所提出的框架将二次六面体有限元、子单元设计变量以及基于均匀化的拓扑优化方法整合到了基于密度的公式体系中。二次六面体有限元能够在相对粗糙的网格上提高均匀化分析的精度，而子单元设计变量则无需对有限元网格进行同等程度的细化，就能提升设计精度。此外，该方法还被实现了为简洁的教育用MATLAB代码，同时还提供了若干基准测试问题，包括设计最大剪模量、最大体积模量以及负泊松比的材料。本文的其余部分结构如下：第2节介绍了在拓扑优化背景下，将基于密度的方法与子单元设计变量以及数值均匀化相结合的理论基础。第3节和第4节分别描述了均匀化方法的MATLAB实现以及三维材料设计问题的处理方法。第5节给出了若干晶格材料设计的基准案例，以及相应的数值实现结果。最后，第6节总结了全文内容。方法部分首先介绍了提升计算效率的基本概念，随后阐述了所提框架中使用的均匀化流程以及过滤和投影方法。均匀化方法的MATLAB实现部分介绍了基于均匀化的方法对应的MATLAB实现代码（homo3d.m）。该方法用于计算材料属性及其相对于拓扑优化设计变量x的灵敏度。为确保计算的准确性与效率，代码中使用HEX20函数来生成单元设计域的三维网格。此外，还介绍了计算单元刚度矩阵和整体刚度矩阵的流程、为高阶有限元设置周期边界条件的方法，以及拓扑优化的相关内容。在本文中，我们定义了几类优化问题，包括设计最大剪模量、最大体积模量以及负泊松比的晶格材料。本节介绍了SCTOP.m函数的MATLAB实现方式。数值示例部分通过若干数值示例，验证了所提出的MATLAB代码的计算效率与实用性。本文中所使用的基体材料具有各向同性特性，其杨氏模量为E=1，泊松比为ν=0.3，单元尺寸为0.1×0.1×0.1，所有这些量均为无量纲值。对于晶格材料设计问题而言，设计变量的初始分布会对最优解产生很大影响，尤其是在设计超材料时。结论本文提出了一种高效且简洁的MATLAB实现方案，用于解决三维晶格材料设计中的计算效率问题。为同时兼顾精度与效率，该方案采用了高阶有限元，并将其细分为多个子单元，从而大幅降低了有限元模型中的自由度。在该框架下，我们实现了基于子单元的均匀化流程，随后进行了相应的灵敏度分析。CRediT作者贡献说明赵力：撰写——初稿、软件、方法论、研究实施、正式分析、数据整理、概念构建。王启涵：撰写——审阅与编辑、软件、方法论、概念构建。张明辉：撰写——审阅与编辑、方法论、概念构建。高伟：撰写——审阅与编辑、方法论、资金获取、概念构建。罗振：撰写——审阅与编辑、监督、资源协调、方法论、资金获取、概念构建。利益冲突声明作者声明自己不存在任何可能影响本文研究成果的已知财务利益或个人关系。致谢本文的研究工作得到了澳大利亚研究委员会项目DP210101353和DP240102559的支持。