基于联合算法的连续体结构拓扑优化新方法

齐 峻，武 星

（长安大学工程机械学院，陕西 西安 710064）

0 引言

拓扑优化算法自20世纪被提出以来，一直在不断更新。拓扑优化的目标是在满足各种设计约束的情况下，找到材料的最优分布[1]。双向渐进结构优化法（BESO）是一种将渐进结构优化法（ESO）和材料添加策略算法（AESO）相结合的方法[2]。最初的ESO主要是去除最低应力单元，经过发展，BESO在ESO原有的基础上可实现增加单元。原始的硬杀伤ESO/BESO方法消除了设计领域的元素，通常，作为设计变量的元件不能直接从设计领域中消除，除非软元件完全等同于空心元件。换言之，对于一个新的拓扑优化问题，开发一种软杀伤BESO方法，然后探索硬杀伤元件的可能性可能会更可靠[3-5]。

近年来，遗传算法（GAs）[6]也不断地发展，并逐渐应用于拓扑优化问题[7]。由于GAs的随机性，其直接应用于拓扑优化问题时会产生棋盘格现象。Liu等[8]很好地意识到了GAs和ESO结合的问题，但由于在设计过程中无法实现材料的增加，得到的结构可能无法满足使用要求。Zuo等[9]提出了GAs和BESO相结合的拓扑优化方法，实现了迭代过程中材料的增减。NSGA-II是在GAs的基础上发展起来的，通过快速的非支配排序和拥挤计算、竞争选择和精英策略来实现快速优化过程。在全局优化方面，NSGA-II显现出比GAs更大的优势。

本文主要部分说明如下：第1节，阐述了拓扑优化的基本问题；第2节，根据传统的灵敏度分析，推导出了固体和软单元的灵敏度数值；第3节，简要说明了本文程序的原理和数值实现；第4节，通过具体的数值算例证明所提出的新方法的稳定性和有效性。

1 问题陈述

在拓扑优化问题中，往往追求给定体积分数下的最大刚度。刚度通常用应变能数值表示，对于具有n个离散单元的拓扑优化问题可以表示为：

式中，C为平均柔度，在设计中，通过计算平均柔度来达到最大刚度。u为全局位移向量，K为全局刚度矩阵。vi是单元i的体积，V*是目标体积。

2 灵敏度数值计算和过滤方案

2.1 灵敏度数值计算

设计变量通过基于梯度的方法更新。对于单相固体材料-空洞优化问题，通过推导目标函数可以得到单元灵敏度值：

式中，p表示惩罚指数，xi表示单相固体材料的单元密度值，当某单元固相材料处于激发状态时，xi=1，否则，关闭此激活状态，单元最终转化为空洞状态，此时xi=xmin。ui表示第i个单元的位移矢量，ki表示第i个单元的刚度矩阵。

灵敏度值可用式（5）表示：

2.2 过滤方案

采用灵敏度滤波策略对灵敏度进行修正。这不仅有效减少了棋盘格现象，还通过过滤半径内物理单元的灵敏度值推导出空单元的灵敏度，实现单元的增删。本文的过滤策略是将与节点关联的所有元素的灵敏度值均质化为节点的元素灵敏度值，如式（6）所示：

式中，H为设计域的节点总数，w(rik)为权重因子，rik表示元素i的中心到节点k的距离。权重因子可以表示为：

式中，rmin表示过滤器半径。为了提高方法优化的收敛性，采用考虑历史迭代信息的收敛策略来修改灵敏度值。具体操作是计算本次迭代和上次迭代的灵敏度的平均值，将其作为灵敏度的更新值。

式中， ,ˆigα 和 ,1ˆig-α 分别为单元i第g次和g-1次迭代的灵敏度值。最后，根据该灵敏度值测量各单元对结构性能的影响，然后将体积作为拓扑优化设计的约束条件。

3 算法的提出与实现

3.1 算法提出

在本文提出的算法中，每个元素都被视为一个独立的个体进行分析，在迭代过程中保持种群维数稳定，允许新种群代替原种群，但不允许合并。每个单元都由一定长度的二进制字符串表示。初始字符串在本文中定义为全“1”的数组，表示所有单元均为固相材料。一个全为“0”的数组意味着该单元将在迭代中被删除。数字“1”代表材料属性的强度。交叉和变异算子只有“0”和“1”之间的转换。

通过初始化种群，实现设计域内所有个体的编码，单位编码的初始值为长度为6的矩阵。对模型进行预处理后，按式（5）计算各元素的灵敏度值。此时，可以根据每个迭代步骤的目标体积分数设置灵敏度阈值，并根据元素灵敏度值进行排序。对于具有相同灵敏度的单元，计算它们的拥挤程度，使同级个体之间存在清晰的层级关系。如果种群中目标函数的最大值等于目标函数的最小值，则将种群的拥挤程度设置为无穷大。一般来说，拥挤的表达式如下：

3.2 算法实现

本文算法步骤如下：

1）通过有限元思想进行预处理，主要是网格划分、施加载荷和设置边界条件[10]。

2）初始化种群，个体由一系列“1”码组成，根据灵敏度计算公式加上惩罚因子和过滤方案的实现来计算灵敏度。指定惩罚指数（p）和体积进化率（ER）等参数。

3）根据灵敏度值迭代进化种群。主要进行快速的非支配排序和交叉变异操作。

4）删除设计域中单元基因串中仅包含数字“0”的元素。如果上次迭代中删除的元素在当前生成中具有高灵敏度值，则将该元素添加到设计域中。

5）根据灵敏度级别和主导关系将实体材料属性分配给设计域中的其余元素。

6）检查设计域的剩余体积和收敛标准。

7）重复步骤2至步骤6，直到达到目标体积和收敛标准。

对于收敛准则，当满足体积约束时，考虑了两个连续优化迭代之间最大柔顺度的绝对变化，两次连续优化迭代的收敛性之差小于0.005。本文中S=5，即连续优化迭代次数为5。

4 数值算例

本节通过数值算例来验证本文提出的算法的有效性、稳定性和优越性。考虑了一种虚拟材料，材料的物理性质如表1所示。所涉及的算法均在Win10 Intel®CoreTM i5-9300 CPU、2.40 GHz、16 G运行内存的桌面上执行。

表1 虚拟材料的物理性质

如图1所示，考虑了尺寸为50 mm×30 mm的二维悬臂梁结构。悬臂梁结构左侧完全固定，右侧中间施加1 N的集中荷载。将体积约束设置为初始体积的35%。网格大小设置为1，共1 500个单元，惩罚指数p=3，过滤器半径rmin=1.5 mm，体积进化率ER=0.02。

图1 悬臂梁模型

4组悬臂梁结构最终迭代结果如图2所示。从图2可以看出，模型会产生良好的布局效果，边界更清晰。由于NSGA-II具有随机性，每次迭代的结果可能会有所不同，但比较4次运行的结果，结构布局显示出相似之处。悬臂梁模型在设置相同参数条件下通过BESO方法所迭代优化出的最终结果如图3所示。4组悬臂梁最终平均柔度结果如表2所示。通过比较表2中4个运行结果的目标函数，相对于平均值的误差和最大误差均不超过1%，4个运行结果的最终体积分数为35%，满足约束条件。通过数值算例分析，证明该算法具有良好的稳定性和收敛性。

图2 4组悬臂梁结构最终迭代结果

图3 BESO方法下悬臂梁最终结果

表2 4组悬臂梁最终平均柔度结果

与BESO方法相比，本文提出的新方法具有更强的进化能力。传统的BESO方法在每次迭代过程中直接增删满足条件的单元，但不同的是本文提出的新方法引入了NSGA-II的演化机制。在迭代过程中，不会直接删除敏感度阈值以下的单元格或增加单元格。相反，将“0”和“1”的数量用作添加和删除的标准。不断地通过“0”和“1”之间的演化，当一个单元出现全“0”的字符串时，就意味着一次删除操作。该算法避免了每次迭代过程中的多次增删改查，有利于算法逐步优化直至收敛。

本文提出的新算法与原始BESO方法的对比如表3所示，设置相同的约束、材料和参数。由于本文提出的新方法具有很强的全局搜索能力，保证了优化的效果。本文算法与原始BESO方法相比，整体在效果上和原始BESO方法近乎相同。考虑到遗传算法具有随机性，其在灵敏度计算、快速非支配排序这两个方面每次的运行时间会较长，但是通过对比本文算法与原始BESO方法的总运行时间与总迭代次数，发现二者几乎无差距，证明了本算法的效率和可靠性。

表3 悬臂梁最终平均柔度比较结果

5 结论

研究小组提出了一种新的结构轻量化设计拓扑优化方法。通过将NSGA-II算法融入BESO方法中，保证了算法的全局优化效果，实现了结构的刚度最大化，获得了干净的拓扑结构。通过数值算例分析，表明本文提出的新方法也能得到类似于BESO的结构布局，证明了算法的稳定性和有效性。本文算法通过结合遗传算法达到了全局最优搜索的效果，在相同约束条件下与传统BESO算法比较，发现本文算法可以提高结构刚度，证明了本文算法的优越性。显然，研究小组提出的方法可以更好地满足复杂工况下的轻量化结构设计，可以应用于各个领域的工程问题，提高结构的性能，对工程结构的优化具有指导意义。