引入水平集带的参数化水平集结构拓扑优化方法

魏鹏 蒋子润

(1.华南理工大学 土木与交通学院∥亚热带建筑科学国家重点实验室，广东 广州 510640；2.广州大学 广东省地震工程与应用技术重点实验室，广东 广州 510006)

至今为止，变密度法[13-14]仍是应用最为广泛的拓扑优化方法，常被用来求解柔顺度最小化问题[15]、柔顺机构优化问题[16]，可靠性优化问题[17]等。变密度法采用0-1之间连续分布的单元密度表示结构中的材料分布，0和1分别表示孔洞和实体材料，优化过程中中间密度逐渐趋近于0和1，使材料分布最终收敛于0-1解并得到最优结果。与水平集方法相比，由于中间密度的存在，变密度法的优化过程具有更高的连续性，对初始设计的依赖性更小，也更容易得到最优解。改进的参数化水平集方法[12]克服了传统水平集方法不能自动开孔的缺陷，但在优化过程中材料密度根据水平集曲面与零水平集的切割关系严格服从0-1分布，容易造成结构边界附近的拓扑突变，降低演化过程的连续性，导致优化结果难以收敛到最优解。当改变初始设计时(孔洞数量和分布)，常常得到不一样的优化结果，即优化结果对初始设计有依赖性。

为提高结构拓扑演化的连续性，Wei等[18]提出了水平集带方法，即在零水平集附近引入一个水平集带区域，区域内的材料密度通过水平集函数进行插值，区域外的材料密度则采用0-1分布。在优化过程中通过减小水平集带宽度逐渐去除区域内的中间材料密度，使目标函数与约束函数在优化结果收敛到0-1解前保持连续性，从而提高优化过程的稳定性，使优化结果更容易趋近于最优解。本文工作在文献[18]基础上，针对水平集方法中常见的初始设计依赖性问题进行了研究，结果表明水平集带的引入可以显著减轻水平集方法的初始设计依赖性。随后探讨了带宽减小速度对最优拓扑的影响，最后通过算例验证了该方法在采用不规则设计域和非结构化网格优化问题中的有效性。

1 引入水平集带的参数化水平集方法

本节将首先简单介绍参数化水平集方法，随后详细说明在参数化水平集方法中引入水平集带的拓扑优化方法。

1.1 参数化水平集方法

在水平集方法中，通过引入高一维的水平集函数的零等值面描述实体材料区域边界，即结构边界。采用水平集方法对实体材料区域内部、外部和边界的定义可以表示为：

(1)

式中，x∈D⊂{(x,y)|x,y∈}为设计域D中的任意一点，t是时间，φ(x,t)为水平集函数，Ω和∂Ω分别表示实体材料的内部区域和边界。

在乡村美食方面，里昂的联合检查组与旅游办公室合作认证的里昂传统餐厅的代名词“Bouchon Lyonnais”，可以证明餐厅不仅具有质量保证，还能烹制最典型和最正宗的里昂菜，对推动乡村饮食的规范化经营起到了重要作用。

在传统水平集方法中，需要通过求解以下Hamilton-Jacobi偏微分方程更新水平集函数：

(2)

(3)

式中，gk(x)为第k个径向基函数，αk(t)为相应的展开系数。将式(3)代入式(2)中，则更新水平集函数的Hamilton-Jacobi偏微分方程可表示为更新展开系数的常微分方程：

α(ti+1)=α(ti)+G-1ΔtVnδ(φ)

(4)

1.2 在参数化水平集方法中引入水平集带

由于水平集方法侧重于边界演化，有限元分析模型缺少单元中间密度，基于水平集方法的拓扑优化过程容易因为拓扑变化使目标函数与约束函数发生突变，导致优化过程不稳定，引起严重的局部最优问题。为解决此问题，Wei等[18]提出了水平集带方法。在传统水平集方法的零水平集附近引入一个包含一定宽度的带状区域，在这个区域内，材料参数由水平集函数值的插值函数确定。如图1所示，图1(a)描述了水平集函数曲面和水平集带，图1(b)显示了水平集带上下边界与水平集函数曲面的交线及材料密度分布。上下边界之间为水平集带区域，其中材料密度的分布与密度法类似，为0-1之间连续变化的数值，由插值函数与水平集函数值决定。水平集带的带宽是上边界和下边界之间的宽度，表达式如下：

图1 水平集带方法示意图Fig.1 Schematic of the level set band method

φb=|φu-φl|

(5)

式中，φu、φl分别代表水平集带的上边界和下边界。在引入水平集带的水平集方法中，单元密度通过密度插值函数得到，第i个单元的密度值为：

ρi=H(φi)

(6)

文中采用光滑的Heaviside函数，表达式如下：

(7)

式中,h=φu=-φl，α和ε是趋近于0的正常数，本文算例均取0.001。

图2是密度插值方案的示意图。图2(a)表示当水平集带的带宽较大时，整个水平集函数曲面均处于水平集带的上下边界之内，则单元密度均处于ε和1之间；图2(b)表示当水平集带的带宽减小时，处于上边界以上的水平集曲面对应的单元密度等于1，处于下边界以下的水平集曲面对应的单元密度等于ε，上下边界之间的水平集曲面对应的单元密度处于ε和1之间。随着优化过程的进行，可以通过减小带宽，使灰度区域逐渐减小；图2(c)表示当水平集带的带宽下降为零时，即上下边界重合，整个结构的单元密度接近0-1分布，优化结果呈现出清晰、光滑的设计边界。

图2 基于水平集带方法的密度插值方法Fig.2 Density interpolation scheme based on the level set band method

2 拓扑优化模型

本文中仅考虑柔顺度最小化问题，目标函数和约束函数分别表示如下：

(8)

(9)

(10)

(11)

式中，b表示体力。驱动边界运动的法向速度场Vn由边界扩展到整个设计域，这由基于形状导数的灵敏度分析推导得到[4-5]，其表达式如下：

Vn=-ε(u):C:ε(u)-ε(u):C:ε(v)+

b·v+(·v)·n+κ(·v)-

(12)

式中，κ=·n为边界曲率，为控制体积约束的拉格朗日乘子，可以由拉格朗日乘子法或二分法计算得到。在本文的柔顺度最小化问题中，不考虑体力和表面分布力，所以法向速度场可以简化为下式：

(13)

通过自然速度扩展，可将上式计算得到的边界法向速度场扩展到整个设计域，代入式(4)中求出更新后的α值，再采用式(3)计算出水平集函数值。迭代过程中通过引入并逐渐改变水平集带带宽φb实现水平集带方法，得到不同区域材料的本构关系：

C=ρ(φ)C0

(14)

式中，C0是材料的原始本构关系。在每次更新带宽之后，需要通过式(7)重新计算结构单元密度分布，直至带宽减小为0，最终得到具有清晰、光滑边界的优化结果。优化过程的流程如图3所示。

图3 优化过程的流程图Fig.3 Flow chart of the optimization process

3 数值算例

本节采用数值算例验证引入水平集带的参数化水平集方法可以减轻初始设计对优化结果的影响，并将该方法拓展到不规则设计域优化问题。几个重要的参数如下：实体区域材料的杨氏模量E=1，孔洞的杨氏模量E=10-9，泊松比υ=0.3,当最终的目标函数值与前9个迭代步的目标函数值均相差小于0.1%时达到收敛，优化过程终止。算例中采用的径向基函数的类型为MultiQuadric(MQ)样条，其位置分布与有限元网格的节点重合。

3.1 算例1

如图4所示的简支梁设计域的长和宽分别是80和20，有限元网格划分为80×20，单元类型为四节点双线性正方形单元，体积分数设置为50%。径向基函数的节点数量为81×21，与有限元网格的节点重合。边界条件为两端下方简支，并在上端中点作用一个F=1N的集中力。本算例通过对5种不同的初始设计得到的优化结果进行比较，研究引入水平集带方法对初始设计依赖性的改进效果。5种不同的初始设计如表1所示，其中第5种初始设计表示随机分布的水平集函数初始值。水平集带的初始带宽为6，每一个迭代步减小0.1，最终带宽是0。

表1 不同初始设计对应的优化结果Table 1 Optimization results corresponding to different initial designs

图4 简支梁的设计域和边界条件Fig.4 Design domain and boundary condition of simply supported beam

表2给出了第1种初始设计对应的优化迭代过程。第2、3、4列分别对应零水平集、材料密度分布和水平集曲面。水平集函数值的初始范围是[-3,3]，初始水平集曲面全部处于水平集带的上下边界之内，因此初始的密度分布图是充满灰度的。优化过程中，水平集函数值及水平集带的带宽发生改变，部分水平集曲面逐渐处于水平集带的上下边界之外，这部分水平集曲面对应的单元密度值等于1或ε，密度分布逐渐趋近于黑白分明。当迭代至60步时，水平集带的带宽减小为0，当迭代至87步时，达到收敛条件，优化过程结束。

表2 第1种初始设计对应的优化过程Table 2 Optimization process corresponding to the first initial design

表1展示了不同的初始设计对应的优化结果，从左到右分别为初始设计、引入水平集带方法的优化结果密度分布、引入水平集带方法的优化结果零水平集、未引入水平集带方法的优化结果零水平集。由第2、3列可知，水平集带宽度为0时，零水平集和密度分布一致。从第3列可知，即使初始设计不同，引入水平集带后也可以得到形状和结构形式相似的优化结果，而从第4列可知，未引入水平集带时，得到的优化结果相差较大，因此可以看出引入水平集带方法可以有效地减轻初始设计对优化结果的影响。表3给出了引入和未引入水平集带方法时优化结果的目标函数值。由表3可知，引入水平集带方法时，不同初始设计对应的目标函数值更加接近，而且大部分情况下小于未引入水平集带方法时对应的目标函数值，说明引入水平集带方法减小了参数化水平集方法对初始设计的依赖性并容易得到更优的优化结果。

表3 引入和未引入水平集带方法的优化目标函数值Table 3 Objective function values with and without the level set band

3.2 算例2

在3.1节中，算例1的优化结果说明引入水平集带的参数化水平集方法可以有效地减轻初始设计的影响，本算例将研究不同的水平集带宽度减小速度对优化结果的影响。仍然采用算例1中的简支梁优化问题，边界条件、荷载、单元类型和数量与算例1相同，采用表2中的前四种初始设计，在每一种初始设计中设置不同的水平集带宽度减小速度，初始带宽均为6，考虑3种迭代步长每步分别减小0.1、0.2、0.4，最终带宽均为0。

表4展示了不同初始设计在不同的水平集带宽度减小速度下得到的最优拓扑。可以看出，除了虚线框内的两个最优拓扑不一致，其他结果均具有相似的形状和拓扑形式，说明最优结果对水平集带宽度减小的速度并不敏感。结合算例1中的结论可知，水平集带的存在有效地减轻了初始设计对优化结果的影响。从表5可知，当水平集带宽度减小速度不一致时，不同初始设计对应的目标函数值变化幅度很小，大部分结果小于未引入水平集带时的目标函数值，进一步论证了优化结果的合理性。

表4 不同的水平集带宽度减小速度对应的最优拓扑Table 4 Optimal topologies corresponding to the different speed reduction of the level set band

表5 不同水平集带带宽减小速度对应的目标函数值Table 5 Objective function values corresponding to the different speed reduction of the level set band

3.3 弯曲悬臂梁

参数化水平集方法的优势之一是可以方便地求解不规则设计域的优化问题，因此，本节将把水平集带方法应用到采用不规则设计域的优化问题中。图5给出了弯曲悬臂梁的设计域和边界条件，弯曲悬臂梁的左端固定，在右下角作用一个F=1N 的集中力。由于设计域是不规则的，所以采用线性三角形单元，网格数量为5 970，节点数量为3 113，体积分数为55%。径向基函数的节点数量为3 113，与有限元网格的节点重合。水平集函数初始值为1，水平集带的初始带宽为6，每迭代步减小0.2，最终带宽是0。

图5 弯曲悬臂梁的设计域和边界条件Fig.5 Design domain and boundary condition of the serpentine beam

图6是弯曲悬臂梁算例的非结构化有限元网格和拓扑优化结果。边界光滑的优化结果类似于由若干粗细不一的杆件组成的桁架结构。图7是弯曲悬臂梁算例的目标函数和体积分数的迭代曲线，展示了良好的收敛性能。该结果与文献[20]和文献[21]中的结果相似，表明本文的优化结果是合理的，验证了水平集带方法可以直接拓展到非结构化网格中用于求解实际工程问题。

图6 弯曲悬臂梁的最优拓扑Fig.6 Optimal topology of the serpentine beam

图7 弯曲悬臂梁优化问题的迭代曲线Fig.7 Iterative curves of the optimization problem for the serpentine beam

4 结语

文中分析了水平集方法中拓扑变化的不连续性问题，并提出通过引入水平集带方法可以有效地减轻初始设计对优化结果的影响，并能一定程度上得到更优的目标函数值。数值算例验证了水平集带方法可以有效地减轻水平集方法对初始设计的依赖性,并显示优化结果对水平集带宽度的减小速度不敏感。最后采用弯曲悬臂梁算例说明该方法可以直接拓展到不规则设计域优化问题中用于求解复杂的实际工程问题，并展示出良好的收敛性能。

文中的数值算例没有考虑表面力的作用，由于在水平集带的带宽减小到0之前，结构没有清晰的边界，因此表面力的存在会对带宽的大小有一定的限制。此外，本文方法还可用于处理带有中间密度的材料如孔隙材料的优化问题，这些也是后续将要开展的研究工作。