带有稳定性约束的几何非线性框架结构拓扑优化

《Computers & Structures》：Topology optimization of geometrically nonlinear frame structures with stability constraints

　　

摘要
本研究提出了一种具有稳定性约束的几何非线性框架优化计算策略，旨在解决两个关键挑战：非线性仿真中由严重结构变形引起的数值不稳定，以及降低临界载荷因子估计精度的虚假模态。这些问题通常出现在低刚度区域，并在拓扑发生剧烈变化时显现。以往的研究在离散设置中应对这些挑战的能力有限，特别是在截面消失累积的情况下。为此，本研究的优化采用双重设计变量，包括定义梁布局的拓扑变量和确定尺寸的截面变量。因此，连续体拓扑优化中采用的方法可以为本研究提供灵感。首先，提出了一种基于能量插值的自适应迭代方案以获得平衡解，并采用基于特征值分析的近似方法评估临界载荷因子。其次，拓扑变量的Heaviside投影实现了刚度贡献单元和非刚度贡献单元的快速聚类，并提出了一种基于模态应变能的统一策略来抑制虚假模态。此外，开发了一种双层优化程序以生成具有标准构件的优化设计，并给出了若干算例以验证所提方法的有效性。

引言
在过去的几十年中，拓扑优化因其能在满足结构性能要求的同时实现显著的材料节约而广受欢迎，然而大多数现有研究集中于连续体自由形态设计，其非标准构件难以施工[1]，[2]，[3]。因此，一些研究建议进行桁架/框架优化以获得优化设计[4]，[5]，[6]。这激发了对桁架/框架优化的研究，其中地面结构方法被广泛采用以识别杆/梁的最优布局[4]，[5]，[6]。关键在于，这些工作大多在弹性假设下进行，排除了几何非线性行为及相关的二阶效应，这种简化可能导致设计缺陷。主要障碍在于优化过程中低刚度单元引起的数值困难（例如严重变形和虚假模态）。为弥补这一不足，本文开发了一种几何非线性框架拓扑优化计算策略，明确纳入整体稳定性约束。

地面结构方法因其能提供扩展的初始设计域而受到广泛关注，在该框架下已完成许多优秀作品以寻求满足特定要求的经济设计[4]，[5]，[6]。除早期的刚度相关设计外，许多研究指出，优化过程中出现的细长构件可能导致结构失稳，并主张对每根杆件施加欧拉屈曲约束以强制局部稳定性[7]，[8]。然而，该方法在优化公式中引入了大量局部约束，显著增加了计算负担。此外，最优桁架结构的整体稳定性无法得到保证，因为某些共线连接的桁架单元可能合并为更长的单元，这种现象被称为"铰链消除"。随后，一些研究者建议在桁架优化中包含几何缺陷或重叠杆件以实现结构屈曲控制[9]，[10]，[11]。然而，由于未考虑整体稳定性要求，结构仍面临失稳破坏的潜在风险。

整体结构稳定性可通过刚度矩阵的正定性进行评估。在弹性条件下，通过假设几何刚度与外力成线性比例关系，可以建立近似特征值问题[11]，[12]：

KL + λi Kg φi = 0

其中KL和Kg分别为整体弹性刚度和几何刚度矩阵。特征对λi, φi按升序排列，最小特征值λ1通常用于近似临界载荷因子。对λ1施加约束已被广泛用于增强桁架拓扑优化的稳定性。然而，由于桁架单元仅承受轴力，忽略弯曲和扭转，因此局部失稳可能无法被捕获。为此，一些研究者主张使用框架单元以提高最优设计的实用性[13]，[14]，[15]。

大多数现有研究是在线性弹性假设下进行的，但屈曲本质上属于几何非线性问题，线性近似往往无法准确预测结构承载能力。因此，近期研究主张将几何非线性分析纳入拓扑优化框架，以获得更准确的临界载荷因子[16]，[17]，[18]。相应地，开发了几种非线性分析求解器，使拓扑优化能够考虑大变形和初始几何缺陷，并已被证明能产生可靠的高性能优化设计[19]，[20]。几何非线性拓扑优化的一个主要挑战是低刚度单元的过度变形，这会导致数值不稳定，甚至使Newton-Raphson（NR）迭代程序无法收敛。这种变形在连续体拓扑优化中常被报告（图1a），已有大量研究致力于解决此问题[21]，[22]，[23]，[24]。在这些解决方案中，能量插值方案将空体单元建模为弹性单元，已被提出并广泛采用于连续体拓扑优化中以减轻此类过度变形（图1b）[21]。此外，随着拓扑演化，由于NR程序中使用的均匀载荷增量往往无法追踪平衡路径，非线性响应变得不可靠，从而损害结构可靠性[25]，[26]。然而，这些问题在框架拓扑优化中更为严重，因为仅使用截面变量时，弯曲刚度的退化速度远快于轴向刚度。据作者所知，很少有研究致力于解决框架拓扑优化中的这一问题。

另一个关键问题在于临界载荷因子的评估。虽然基于弧长技术的直接平衡路径追踪方法能准确识别实际临界载荷，但其在多次优化迭代中的极高计算成本限制了其在拓扑优化中的应用，特别是对于大规模三维工程结构[27]，[28]，[29]。另一种方法是Lindgaard等人提出的路径追踪方法，通过特征值分析预测即将到来的分叉点，但该方法仅作为临界载荷最大化的目标函数，而非约束条件[30]，[31]。此外，基于刚度降阶的隐式方法已被用于估计优化填充结构的临界载荷[32]。然而，这些方法仅在优化后分析中实施，并未集成到优化过程中。为在精度和计算成本之间取得平衡，近期研究广泛主张在非线性平衡点（也称为几何非线性特征值分析）执行特征值分析以近似临界载荷因子。其基本思想是，当切线刚度矩阵失去正定性时，结构发生失稳，一种标准公式如下[33]：

KT - λi I φi = 0

其中KT为整体切线刚度矩阵，I为单位矩阵。然而，该公式在协调网格下执行，因此不适用于采用非协调网格的地面结构方法。另一种替代方案是将收敛平衡步的切线刚度信息外推到临界点。该方法假设几何刚度矩阵Kg（或与位移刚度Ku一起）随载荷因子线性变化。临界载荷因子可使用以下特征值公式评估[34]，[35]，[36]：

KL + Ku + λi Kg φi = 0 或 KL + λi Kg + Ku φi = 0 （3）

虽然直接方法能提供更准确的临界载荷因子预测，但其在拓扑优化中的应用鲁棒性较差，因为模态跳跃可能导致无法识别真正的临界载荷。因此，在对拓扑优化中不同屈曲分析方法进行全面比较后，推荐采用近似特征值方法（公式3）[18]。

对λ1施加约束的一个关键挑战在于，最小算子本质上不可微，并且随着结构拓扑演化可能发生所谓的"模态切换"现象[37]，[38]。为此，多阶屈曲载荷因子约束结合聚合函数方法已被用于形成等价的微分约束。此外，频繁提到的另一个问题是低刚度单元引起的"虚假屈曲模态"，该问题类似于动力学优化中的"局域化模态"[39]，[40]。在连续体拓扑优化中，一种常见的补救方法是使用材料惩罚模型来修正低密度区域的单元刚度和/或质量矩阵，以适当匹配并缓解低阶虚假模态。此外，一些研究者认识到，这些低刚度单元在屈曲模态中的参与度应控制在较低水平[41]，[42]。基于这一思路，开发了一种基于低刚度单元模态能量比的虚假模态识别方法[13]，[23]，[39]。然而，该方法需要若干额外参数来识别低刚度单元并定义能量比阈值，这限制了其适用性。

此外，经典地面结构方法通常使用截面参数作为设计变量，当梁的截面趋近下界时该梁被消除。然而，这种公式化触发了虚假模态问题，原因在于上述离散设置中截面消失累积时轴向刚度和弯曲刚度的差异性退化。一些研究者认识到这一问题并提出引入额外的二元变量来表示梁的存在与否[23]，[43]，[44]。然而，这种补救方法不能完全解决框架拓扑优化中的虚假模态问题，因为底层截面参数仍可趋近于零，导致残余低刚度构件。

在本文中，优化框架采用拓扑变量表示梁的存在，并将截面变量投影到非消失截面的标准离散构件尺寸上。因此，连续体结构的计算策略可被借鉴以在将稳定性要求纳入几何非线性框架拓扑优化时解决网格变形和虚假模态问题。具体而言，开发了一种基于载荷控制方法和能量插值的自适应迭代方案以减轻低刚度单元的过度变形。随着Heaviside投影的尖锐度增加，可以实现刚度贡献单元和非贡献单元的快速聚类，并提出了一种基于非贡献单元模态应变能的统一策略来抑制虚假屈曲模态。最后，提出了一种基于MMA算法的双层嵌套优化程序，以在几何非线性场景下求解优化问题。给出了若干数值算例以展示所提方法的优势。

本文其余部分结构如下：第2节介绍几何非线性框架分析和临界载荷因子评估的近似特征值分析。第3节建立拓扑变量和截面变量的投影和插值模型，并建立稳定性约束下的综合优化模型。第4节开发自适应迭代方案和一致的虚假屈曲模态识别方法，以解决低刚度构件引起的数值困难。第5节推导目标和约束函数对拓扑和截面设计变量的灵敏度，并详细介绍嵌套优化求解程序。第6节给出若干数值算例以验证所提方法的可行性。第7节总结全文。

章节摘要

非线性屈曲分析
虽然考虑整体稳定性的线性弹性框架拓扑优化已被广泛研究，但此类方法本质上忽略了轴力-弯矩耦合效应。这一疏漏导致优化后的细长受压构件在轴向受压时出现严重的弯曲刚度退化，进而导致整体结构刚度丧失甚至结构失稳。本节简要描述基于Green-Lagrangian应变的几何非线性框架模型，然后……

优化问题描述
本文采用经典地面结构方法生成初始网格化框架结构，候选构件由连接节点对形成。地面结构的复杂度由连接等级lv决定，定义为两个连接节点之间允许的最大跨度。连接等级对4×4节点网格地面结构的影响如图2所示。为在设计域复杂度和计算成本之间取得平衡，采用……

网格变形和虚假模态问题的缓解方法
显然，非线性屈曲特征值分析的准确性严重依赖于NR迭代的收敛性，而这在优化过程中容易受到低刚度构件引起的单元变形破坏。此外，低刚度构件诱导的虚假屈曲模态是损害临界载荷因子评估精度的另一个关键挑战。

为解决这些问题，首先开发了一种基于荷载控制方法的自适应迭代方案。

**灵敏度分析与求解过程**

鉴于上述优化模型中存在大量的拓扑和尺寸设计变量，基于梯度的优化器，尤其是广泛使用的移动渐近线法（MMA）算法，更适合求解所建立的优化问题。因此，本节推导了目标函数和约束函数的解析灵敏度，并开发了一种基于MMA算法的嵌套优化流程以高效求解所提出的问题。

**(a) 末端柔度的灵敏度分析**

为了避免直接计算位移场关于设计变量的导数，采用伴随法计算末端柔度的灵敏度。为此，通过用伴随拉格朗日乘子向量左乘控制平衡方程，然后将所得乘积附加到原始末端柔度上，构造辅助拉格朗日函数：

C?(ρ,t) = C(ρ,t) ? v^T (F_int(ρ,t) ? F_ext)

**(b) 稳定性约束函数的灵敏度分析**

对特征值方程求导得到：

?K_u/?ρ?_i + K_g/?ρ?_i + μ_j ?K_L/?ρ?_i + ?μ_j/?ρ?_i K_L φ_j = 0

考虑到切线刚度矩阵的对称性，用特征向量左乘方程两边，第j个特征值关于设计变量的导数可表示为：

?μ_j/?ρ?_i = ?1/(φ_j^T K_L φ_j) · φ_j^T (?K_u/?ρ?_i + ?K_g/?ρ?_i + Σ ?K_u/?u_l · ?u_l/?ρ?_i + ?K_g/?u_l · ?u_l/?ρ?_i + μ_j ?K_L/?ρ?_i) φ_j

同时，聚合稳定性约束函数关于设计变量的导数给出如下：

**数值算例**

本节首先展示弹性条件下的优化设计作为基准设计。随后，研究考虑几何非线性效应的优化配置，并检验其在不同稳定性要求下的性能，以表明所提出方法的可行性。

对于所有数值算例，地面结构的候选构件按照中国《钢结构设计标准》进行选择。

**结论**

本文提出了一种具有稳定性要求的几何非线性框架拓扑优化的计算策略。通过引入拓扑和尺寸变量，提出了一种自适应迭代求解过程和一致的非真模态识别方案，以处理优化过程中与低刚度单元相关的数值困难。主要结论总结如下：

（1）耦合的拓扑-尺寸设计公式产生了...

**CRediT作者贡献声明**

赵雷：原创草稿撰写、方法学、资金获取、形式分析、概念化。
袁航宇：可视化、验证、软件。
薛金丹：可视化、验证、调查研究。
罗维峰：审稿与编辑、可视化、软件。
张海良：验证、软件、调查研究。
容轩培：验证、软件、形式分析。
韩燕：审稿与编辑、软件、方法学。
容建华：审稿与编辑、方法学。

**利益冲突声明**

作者声明不存在任何已知的竞争性财务利益或个人关系可能影响本文所报告的工作。

**致谢**

本研究由中国国家自然科学基金（批准号：12102066）、教育部桥梁工程安全控制重点实验室开放基金（长沙理工大学）（21KB08）资助。非常感谢审稿人对本文的评论。

赵雷 | 袁航宇 | 薛金丹 | 罗维峰 | 张海良 | 容轩培 | 韩燕 | 容建华