一种用于防止几何非线性拓扑优化中病态条件的特征值稳定化方法

《Energy Research & Social Science》：An eigenvalue stabilization method for preventing Ill-conditioning in geometrically nonlinear topology optimization

　　罗云峰|袁伟|胡静宇|李宗卓|李曲豪|刘树天
中国山东大学机械工程学院先进装备与金属成形技术国家重点实验室，济南，250061

**摘要**
由于网格变形引起的数值奇异性在几何非线性拓扑优化中仍然是一个具有挑战性的问题。大多数传统方法试图通过防止低刚度区域的网格变形来提高收敛性。然而，即使在显著应力作用下的高刚度固体区域也是变形的高风险区，这是这些方法应用于更大变形情况的关键限制。为了解决这个问题，本文提出了一种基于直接数学识别易变形元素的几何非线性拓扑优化方法。具体来说，提出了一种改进的特征值稳定技术，首先对元素刚度矩阵进行特征值分解，然后使用经验公式修正其中的小负特征值，从而防止全局刚度矩阵变得病态。为了不改变物理问题，在载荷侧添加相应的修正项以平衡修正后的刚度影响。该方法无缝集成到拓扑优化框架中，仅在有限元分析期间需要局部稳定，而不会改变敏感性分析或优化迭代。其主要优点包括提高收敛性、保持问题的物理完整性、易于实现以及具有很高的灵活性。数值示例证明了所提出方法在处理具有严重网格变形的几何非线性拓扑优化问题中的有效性。

**引言**
几何非线性拓扑优化将大变形引起的几何非线性效应整合到拓扑优化中，以设计能够适应显著形状变化的高性能结构，已成为过去几十年的一个热门研究方向。自从Yuge等人[1]关于几何非线性条件下的结构布局的工作以来，该领域发展迅速，目前应用于许多先进领域，如柔性机构[2,3]、非线性执行器[4,5,6]、具有预定和可编程机械响应的非线性结构[7,8]、大幅变形和可编程的超材料[9,10,11,12]、可拉伸电连接器[13]等。
在大变形拓扑优化中的一个关键挑战是在迭代过程中经常发生在低刚度或中等密度区域的严重网格变形[14,15]。这种变形会导致有限元分析的牛顿-拉夫森迭代中的数值不稳定[16,17]。尽管现代计算力学工具（如成熟的商业软件包Abaqus）能够稳健地处理大变形结构分析，但其可靠的收敛性仍然基本上依赖于简单且固定的结构配置。然而，拓扑优化过程会产生数百个中间设计，这些设计通常具有复杂的几何形状和刚度分布的突变。确保所有中间设计的稳定非线性收敛对于成功的优化过程至关重要，但这一要求严重限制了优化所能处理的最大变形范围。因此，即使经过几十年的发展，大变形拓扑优化在涉及极端变形的实际应用（如柔性电子[18]和软体机器人[19,20]）中仍然很少得到应用。
先前的研究广泛认识到低刚度区域容易发生网格变形，并基于这一观察开发了各种方法来提高收敛稳定性。例如，Buhl等人[21]和Pedersen等人[22]提出在优化过程中忽略低密度区域内节点的内力以抑制变形。Bruns和Tortorelli[23]直接移除了低密度元素以完全消除其影响。Yoon和Kim[24]引入了连接所有元素的虚拟弹簧，赋予它们高刚度以防止过度变形。此外，使用超弹性材料模型可以在一定程度上提高低密度元素的变形能力[25,26,27,28]。例如，Luo等人[29]通过叠加一层超弹性元素增强了低密度区域的剪切阻力。同样，Scherz等人[30]最近引入了一种与变形梯度的条件数相关的超弹性本构材料模型用于数值稳定。值得注意的是，Wang等人[15]提出了一种能量插值方案，在低密度区域应用线性弹性模型，在高密度区域应用超弹性模型。通过抑制低刚度区域的几何非线性效应，这种方法显著提高了计算收敛性，并因其简单性而成为主流方法。在此基础上，Song等人[31]进一步引入了一种基于Neo-Yeoh超弹性的插值方案，选择性地将Yeoh模型应用于低刚度区域以增强其剪切变形能力。显然，这些防止低刚度区域变形的方法需要预定义的刚度阈值来区分“低”和“高”密度。然而，如果阈值设置得太低，则无法完全缓解元素变形的风险；如果设置得太高，则会改变实际的物理问题。这种固有的矛盾从根本上限制了基于低刚度元素易变形假设的方法在推动可管理变形边界方面的潜力。因此，进一步提高大变形拓扑优化的变形能力仍然是该领域的一个关键挑战。
实际上，网格变形是由载荷条件和局部刚度的共同作用引起的。即使是在相对较高刚度的区域，在大的局部载荷下（例如，在载荷作用点或应力集中点）也可能发生变形。相反，低刚度区域在小的载荷下可能不会变形。因此，尽管防止低刚度元素的变形可以在一定程度上提高收敛性，但由于高刚度固体区域的变形，这种方法仍然难以应用于更大变形的情况。一种更有前景的方法可能是直接从数学或机械角度识别真正有变形风险的元素，并施加有针对性的局部稳定。特征值稳定技术[32]就是这样一种实现直接识别和处理易变形元素的方法，其中识别出切线刚度矩阵的小特征值元素并进行局部修改以进行稳定。Loehnert等人[32]首次在扩展的有限元方法中提出了所谓的特征值稳定技术。他们的方法通过使用元素刚度矩阵的特征值分解来解决病态问题，根据特征值的大小选择性地稳定问题模式，从而在保持原始求解框架的同时提高收敛性。随后，Garhuom等人[33]成功将其扩展到有限单元方法（FCM），有效缓解了矩阵病态并提高了FCM中大变形分析的稳健性。然而，这些工作主要关注小特征值对收敛性的影响，而没有完全解决在几何非线性问题中频繁出现的负特征值，这些负特征值也会严重损害数值稳定性。
为此，本文提出了一种改进的特征值稳定方法，同时考虑了小特征值和负特征值的影响，然后将该改进方法成功集成到大变形拓扑优化过程中。它具有几个关键优势：首先，它能够直接且精确地识别高变形风险的元素，从而可能处理涉及更大变形范围的设计问题；其次，其实现非常直接，仅在有限元分析期间需要引入局部稳定，而不会改变敏感性分析或优化迭代过程；第三，它具有很高的通用性和灵活性，因为它不依赖于非线性FEM求解器和优化策略，因此可以很容易地与各种现有的分析和优化方法结合使用。通过几个基准数值示例证明了所提出方法的有效性。

**本文的其余部分组织如下：**
第2节详细介绍了特征值稳定增强拓扑优化方法，然后在第3节提供数值示例，在第4节给出结论。

**数值示例**
**非线性拓扑优化公式**
本工作的主要创新在于将特征值稳定技术引入大变形结构的拓扑优化中，从而解决了由网格变形引起的数值奇异性问题。应当注意的是，本文采用的特征值稳定技术是一种适用于非线性控制方程的通用解决方案技术，理论上可以应用于任何拓扑优化方法。然而，由于篇幅限制，这里不再详细讨论。

**数值示例**
在本节中，提供了几个数值示例来证明所提出方法的有效性。所有示例中的参数保持一致。为了简化，所有设计域都使用1×1的均匀网格和元素尺寸。对于柔顺性优化问题，材料属性固定如下：杨氏模量E=3GPa，泊松比μ=0.4；而对于柔顺性逆变器设计问题，这些参数分别设置为3 MPa和0.3。

**结论**
本文通过引入一种基于与现有方法根本不同视角的新方法，解决了大变形拓扑优化中严重网格变形的关键挑战。我们超越了变形与低刚度区域相关的传统假设，提出了一种改进的特征值稳定技术，通过数学诊断直接识别易变形元素。

**数据可用性**
数据可根据合理请求提供。

**作者贡献声明**
罗云峰：撰写——原始草稿、方法论、研究。
袁伟：软件、研究。
胡静宇：软件、方法论。
李宗卓：撰写——审阅与编辑、软件。
李曲豪：撰写——审阅与编辑、方法论。
刘树天：撰写——审阅与编辑、监督。