基于连续阶跃参考应力的SKO结构拓扑优化方法

侯剑云 丁晓红

上海理工大学，上海，200093

0 引言

结构拓扑优化是指寻求材料分布的最优拓扑形态，其目的是在设计域空间内寻求结构最佳的传力路线形式，以优化结构的某些性能或减轻结构的重量。结构拓扑优化已发展100多年，大体分为解析方法和数值方法［1］。解析方法的理论数学求解困难，不便于在工程实际中直接应用。近几十年来，计算机技术在结构分析中的普遍应用使结构优化数值方法得到了迅猛发展，它主要分为两大类：一类是可解决各类结构的尺寸、形状及拓扑优化，但计算效率和通用性并不理想的均匀化方法（homogenization method）［2］；另一类是以渐进结构优化方法（evolutionary structural optimization）［3－4］为代表的优化算法，它的特点是在宏观的角度上对结构材料分布进行设计，计算效率高、通用性好，具有启发式特征，能得到近似最优解。在工程实际中，近似最优解通常被广泛采纳，因此这一类算法具有良好的发展前景。

SKO（soft kill option）方法是一种启发式的算法，其基本原理是逐渐“软化”低应力的材料，“硬化”高应力的材料，使经过优化后的结构应力水平变得更均匀。它最初是德国Karlsruhe研究中心提出的［5－6］。目前，国内外对SKO方法的研究大体上分为两种，一是采用以体积率作为删除准则的方法［7］，二是以参考应力作为材料删除准则的方法［8］。前者通过多次删除一定的体积使优化结果收敛于目标体积，每次删除过程中将体积率对应的应力值作为删除标准逐渐软化和硬化材料。不同的目标体积可得到不同的拓扑结构，但体积率对应的应力值求取困难，运算耗时长。后者通过结构的应力分布选择参考应力逐渐软化和硬化材料，这种方法参考应力选取简单，运算耗时短，但不同的参考应力对优化结果影响较大。若参考应力取值过大会导致过删除，过小会导致优化的前后体积改变不大。丁晓红等［8］提出了以平均应力、过滤应力等为参考应力的删除准则。在迭代过程中，参考应力与结构的应力分布相关联，由于参考应力随着迭代次数的增加而增大，不同设计问题的变化规律是不同的，材料删除的规律具有不确定性，有可能存在大片区域的应力高于参考应力水平，材料无法继续删除，导致迭代过早收敛，仅达到局部最优。

局部最优是一种寻优过程中的早熟现象，需要避免。为此，本文提出一种基于连续阶跃参考应力的寻优迭代策略，可有效地改善早熟现象，得到最优结构。

1 连续阶跃参考应力SKO方法介绍

1.1 基本原理介绍

假定结构由不同的材料组成，将每个单元的弹性模量作为参数来改变结构的拓扑形态。材料的弹性模量被定义为温度的函数，即随着温度的升高，弹性模量变小，材料被“软化”，当材料“软”到一定程度，可认为材料被删除；同时随着温度的降低，材料的弹性模量变大，材料被“硬化”。温度T和弹性模量E间的关系可假设为线性关系，如图1所示，此处的温度没有物理意义，仅是单元弹性模量改变的控制器。在应用连续阶跃参考应力的SKO方法的过程中，会经过多阶的局部寻优过程，每阶寻优都经过数次迭代。温度的迭代公式为

式中，上标i、j表示第i阶寻优过程中第j次迭代；Tn（i，j）为第i阶寻优过程中第j次迭代时n节点的温度；σn（i，j－1）为第i阶寻优过程中第j－1次迭代时n节点的等效节点应力；σ（i，j）为第i阶局部寻优过程中第j次迭代的连续阶跃参考r应力；s为步长因子；Tmax、Tmin为迭代时的温度上下限，上限为100，下限为0。

图1 温度与弹性模量关系图

由式（1）可知，如果节点n的应力大于参考应力，则节点温度降低，材料的弹性模量增大，材料被“硬化”；否则节点温度升高，材料被“软化”。Emax为所选用材料的弹性模量，用以模拟固体区域 。 当 Tn（i，j－1）在［0，100］之间时，弹性模量和温度之间的关系满足如图1所示的线性关系。如果Tn（i ，j－1）≤0，则E＝Emax；当 Tn（i，j－1）≥100时，则E＝Emin。Emin非常小，接近结构空隙率的数值，经验证明，取Emin＝Emax／1000即可以确保数值稳定［7］。

1.2 局部寻优过程和连续阶跃参考应力

在SKO方法的寻优过程中，有可能出现局部最优结构，它存在着一些材料堆积区域，此时若提高参考应力，则寻优过程又可继续进行。本文提出一种连续阶跃参考应力的概念，使结构能在寻优过程中自动调节参考应力的取值，从而逐步使结构趋优。其计算公式为

其中，ζ（i，j）为第i阶寻优过程中第j次迭代的权重系数；连续阶跃参考应力σr（i，j）是以第i阶寻优过程中初次迭代的平均应力σa（i）为基数乘以一个权重系数ζ（i，j）得 到的。 权 重 系 数ζ（i，j）计 算公 式如下：

式 中 ，ζ0（i）为第i阶寻优权重系数初值；α为常数，取0.1；Δζ为权重系数增量。

第i阶权重系数初值ζ0（i）为第i－1阶寻优过程中最后一次迭代的权重系数，即

由以上公式可知，连续阶跃参考应力的变化趋势与权重系数基本一致。权重系数ζ（i，j）随着寻优阶数的增加而逐阶增大。以ζ0（0）＝0.85、α＝0.1、Δζ＝0.05为参数，分别经过4阶局部寻优过程，来说明权重系数的变化，由于每阶迭代次数不确定，为便于说明，假定每阶寻优的迭代次数均为10。如图2所示，曲线划分为5个区域。j在进入下一阶局部寻优过程时重新计数，i在原来的基础上增加一次，每阶寻优之初，Δζαj与Δζ刚好相等，式（4）、式（5）式变为ζ（i，0）＝ζ0（i），ζ0（i）＝，故局部寻优过程的第一步迭代与上一阶局部寻优过程的最后一步迭代的权重系数相等，函数不发生突变。而随着迭代次数的增加，Δζαj很快衰减为零，最终权重系数ζ（i，j）等于第i阶权重系数初值ζ（i）与权重系数增量Δζ之和。区域1、3、4、5的权

0重系数均满足以上所述的提升方式，而区域2经过的两次迭代都达到局部最优，说明权重系数增大得不够，需再次提升，才能进入下一阶寻优过程，这是权重系数的自动调节过程。

图2 权重系数ζ变化图

1.3 局部最优判断标准和收敛条件

在寻优过程中，当结构前后两次迭代的体积变化小于给定的局部体积容差ε1时，认为结构达到一个局部最优状态；当小于给定的总体体积容差ε时，认为得到最终的优化结果。显然，ε1要大于ε。实验证明，ε1取100～1000倍的ε时可确保数值稳定。这里的体积指的是量纲一化的名义体积。SKO方法中结构的设计变量是弹性模量E，定义第i阶寻优过程中第j次迭代的名义体积为

式 中 ，En（i，j）为第i阶寻优过程中第j次迭代时n单元的弹性模量；hn和An分别为n单元的厚度和面积。

将名义体积量纲一化得到第i阶寻优过程中第j次迭代时量纲一体积指标，其公式为

2 算法流程

本文的算法流程如图3所示。具体步骤如下：

（1）建立初始物理模型，并确定设计空间，实际结构中某些具有特殊要求的边界或面需保留而不作为设计区域。

（2）赋予结构具有如图1所示的材料属性，并用有限元离散网格施加约束和载荷，图3中的T0为环境 温 度，和 T（i，j）n一样，T0没有实际的物理意义。

（3）进行结构线性静力有限元分析，提取节点应力，通过式（3）～式（5）计算连续阶跃参考应力，权重系数ζ（i，j）不得超过规定的权重系数上限ζmax。

图3 算法流程

（4）按式（1）、式（2）实现“软硬化”单元的操作。

（5）未达到局部最优时，更新j，重复步骤（3），当达到局部最优结构后，更新i、j，再重复步骤（3），直到满足上文所述的迭代终止条件，终止迭代。

3 设计参数选择与算例

在迭代过程中，不同的参考应力对优化结果的影响较大。本文的连续阶跃参考应力因权重系数而阶跃式地由小到大变化。影响权重系数的参数有初始权重系数ζ（0）0，权重系数增量Δζ和权重系数上限ζmax。初始权重系数ζ（0）0取得过大，在开始时便使材料删除得过多，由于结构的后续分析基于删除材料后的结构，会造成某些需“硬化”的材料无法恢复，从而最终不能达到最优拓扑结构。权重系数增量Δζ和权重系数上限ζmax，它们分别控制参考应力的增长幅值和上限。本文通过受两面内力悬臂梁的拓扑优化例子探讨不同的权重系数参数对优化结果的影响。

3.1 受两面内力悬臂梁的拓扑优化

如图4a所示的初始物理模型，悬臂梁的长宽比为2，板的一端固定，受两垂直向下的力P。从图4b所示的初始模型应力分布图可以看出，固定端的两角和受力点存在应力集中。在不考虑应力集中区域的情况下，结构的应力值集中在［0.02，51.2］MPa之间。

图4 长方形悬臂梁受两面内力模型

权重系数参数分别按表1所示的4种方案进行选取，局部体积容差ε1和总体体积容差ε如表1所示。

表1 参数取值表

方案1中初始权重系数ζ0（0）与权重系数上限ζmax相等，即采用无连续阶跃参考应力的SKO方法，始终以平均应力的0.8倍作为参考应力进行拓扑优化。优化后材料分布如图5a所示。结构的中间左上区域材料堆积，对比应力分布图（图5b），该区域的应力接近白色，说明应力水平不高，但却大于参考应力而未被软化。若不考虑应力集中区域，存在材料的区域应力值集中在［0，51.6］MPa之间，相比于初始模型的应力水平均匀程度无明显改善，优化结果仅得到局部最优。图5c为量纲一体积指标V、权重系数ζ和量纲一化后的参考应力指标σr的迭代进程图。权重系数和参考应力不变，结构的体积随迭代次数的增加而降低，直到达到总体体积容差而终止，此次历经了221次迭代，最终体积为初始体积的0.61。

图5 方案1的结果图

方案2采用了连续阶跃参考应力的方案，初始权 重 系 数 为ζ0（0）为0.5，权重系数增量Δζ为0.05，权重系数上限ζmax为1.25。优化后材料分布如图6a所示，与图5a相比有较多的细部结构，其应力云图6b在不考虑应力集中区域的情况下，存在材料的区域应力值集中在［17.2，51.5］MPa之间，相比于方案1的应力水平均匀程度有所提高，得到较为理想的近似最优解。图6c为性能指标的迭代进程图。权重系数连续阶跃增大最终达到权重系数上限ζmax。参考应力变化与权重系数基本一致。体积随着迭代次数的增加一直减小，最终趋于总体体积容差而终止。与方案1相比，此次历经了325次迭代，结构的最终体积变为原来的0.41，比方案1多删除一部分。

方案3与方案2相比，权重系数增量Δζ从0.05提升为0.12。优化后的材料分布如图7a所示。方案3所得到的优化结果与方案2的优化结果差别不大，其应力分布图（图7b）与方案2的应力分布差别不大，存在材料区域的应力值集中在［17.2，51.4］MPa之间，相比于方案2的应力水平均匀程度相当，同样得到较为理想的近似最优解。由图7c可知，权重系数的逐阶增幅相应比方案2高，直到趋于权重系数上限ζmax。参考应力增大趋势与权重系数基本一致。体积指标随着迭代次数的增加逐渐趋于水平，最终达到总体体积容差而终止。迭代次数由方案2的325次减少为方案3的229次。迭代效率明显所提高，最终体积为初始体积的0.40。

图6 方案2的结果图

图7 方案3的结果图

方案4与方案2相比，初始权重系数ζ0（0）从0.5增大为1.4，权重系数增量Δζ不变，ζmax取2.0。图8a为优化后的材料分布图，材料分布与方案2和方案3相比有很大不同，它的最终体积虽然减少了很多，为原来的0.23，但由应力分布图（图8b）可知，存在材料的区域应力值集中在［17.1，102.6］MPa之间，应力水平的均匀程度不如方案2和方案3，得到的结果收敛于次优解。因为当初始权重系数ζ0（0）取得过大，迭代初始删除材料过多，造成了后续需“硬化”的材料无法恢复。从图8c可知，权重系数逐阶增大直至趋于权重系数上限ζmax。参考应力在第一阶的增幅很大，后续的增幅与权重系数一致。体积指标随着迭代次数的增加逐渐趋于水平，最终达到总体体积容差而终止。

图8 方案4的结果图

由以上4种方案的优化结果可知，初始权重系数ζ0（0）如果取值过低，迭代过程会被延长，但每次迭代的应力变化会更平稳，这样可以避免应力变化过大引起误删除，故一般在［0.5，1］之间取值。Δζ作为权重系数增量，它决定了每阶寻优过程中参考应力增大幅度，为提高效率可适当提高权重系数增量Δζ，但提高过大也会导致结构的过删除，一般在［0.05，0.15］之间取值。显然权重系数上限ζmax要大于初始权重系数ζ0（0），但取得过大会导致过删除，一般在［1.0，2.5］之间取值。

综上所述，在结构优化过程中，利用连续阶跃参考应力的SKO方法可以更快更好地得到近似最优解，并且得到的优化结果边界清晰，结构合理，应力水平更均匀。

3.2 拱桥的拓扑优化

如图9所示，初始模型为拱桥的平面模型，假设桥底到桥面顶部的距离为H，桥面的曲率半径为R，R／H ＝5，桥底左右两端固支，桥上端面受均布载荷P。

图9 拱桥平面模型

图10a为拱桥优化后材料分布图，图10b为赵州桥实景图，图10c为迭代过程中性能参数的变化图，初始权重系数ζ（0）0为0.5，权重系数增量Δζ为0.12，经过195次迭代结束，权重系逐阶增大，达到ζmax后不再增加，体积指标随着迭代次数增加而逐渐趋于平稳，最后达到总体体积容差而终止。最终体积为原来的0.41，得到的结果与实际拱桥结构相似。

图10 模型的优化结果图

3.3 受水平力的直立穿孔板的拓扑优化

如图11所示，初始模型为一直立的平板，中间有两长方形孔，H／L ＝2，H1／L1＝2，H／H1＝3，平板下端固定，左侧受两水平集中载荷P。

图11 受水平力的直立穿孔板的初始模型

图12a为受水平力的直立穿孔板的优化后材料分布图。如图12b所示，初始权重系数ζ（0）0为0.5，权重系数增量Δζ为0.12，第2、第3阶局部寻优过程中权重系数增大幅度超过权重系数增量，这是因为迭代过程中出现相邻几次迭代均处于局部最优，权重系数自动调节至更高，以便能进入下阶寻优。随着迭代次数的增加，体积变化逐渐趋于平稳，最后达到总体体积容差而终止。最终体积为原来的0.34，结构得到近似最优解。

图12 模型的优化结果图

4 结束语

本文提出的基于连续阶跃参考应力的SKO方法，能在寻优过程中自动调节参考应力值，有效地解决寻优过程中的材料堆积问题，使结构趋于最优。算例表明：基于连续阶跃参考应力的SKO优化方法实现对平面应力的线弹性结构的拓扑优化设计，方法易于实现，收敛性好，得到的结果结构合理，边界清晰，优化后的结构应力水平变得更均匀。该方法也同样适用于三维实体模型的拓扑优化。

［1]周克民，李俊峰，李霞.结构拓扑优化研究方法综述［J］.力学进展，2005，35（1）：69－76.

［2]Bendsoe M P，Kikuchi N.Generating Optimal Tpologies in Structural Design Using Homogenization Method［J］.Computer Methods in Applied Mechanics and Engineering，1988，71（2）：197－224.

［3]Xie Y M，Steven G P.A Simple Evolutionary Procedure for Structural Optimization［J］.Computer ＆Structure，1993，49：885－896.

［4]荣见华，谢忆民，姜节胜，等.渐进结构优化设计的现状与进展［J］.长沙交通学院学报，2001，17（3）：17－23.

［5]Mattheck C.Design in Nature：Learning from Trees［M］.Berlin：Springer－Verlag，1998.

［6]Harzheim L，Graf G.A Review of Optimization of Cast Parts Using Topology Optimization I－Topology Optimization without Manufacturing Constraints［J］.Structural and Multidisciplinary Optimization，2005，30：491－497.

［7]Harzheim L，Graf G.Topshape：An Attempt to Create Design Proposals including Manufacturing Constraints［J］.International Journal of Vehicle Design，2002，28（4）：389－409.

［8]丁晓红，程莉.基于SKO方法的满应力结构拓扑优化设计［J］.中国机械工程，2009，20（8）：1765－1770.