考虑承载面最大位移约束的结构拓扑优化方法

龙凯，陈卓，谷春璐，王选

1. 华北电力大学 新能源电力系统国家重点实验室，北京 102206 2. 合肥工业大学 土木与水利工程学院，合肥 230009

连续体结构拓扑优化设计指在给定的区域内，寻求结构的某种布局(如结构有无孔洞、孔洞位置、数量和结构连接方式等)，使其能够在满足指定约束条件下的设计目标最优[1-2]。拓扑优化在早期研究中以全局性性能如结构柔顺度、振动频率为设计指标，发展到考虑局部应力、疲劳损伤等性能，从解决工程实际问题的角度出发，以改善结构综合性能[3]。

位移或变形量是结构产品刚度的重要衡量指标。Rong等[4-5]在拓扑优化设计中引入多点位移约束，实现了对结构变形的控制。Huang和Xie[6]提出采用柔顺度目标和节点位移约束以实现对全局和局部刚度的控制。Zuo和Xie[7]提出一个全局化位移指标，采用双向进化式拓扑优化方法实现了结构刚度设计。Qiao和Liu[8]提出了几何平均位移指标用于衡量拓扑优化结构的刚度。Li和Kim[9]从工程实际出发，提出了多相材料结构位移约束下的轻量化设计列式和求解算法。Long等[10]提出了多点位移约束拓扑优化问题的二次规划列式与求解算法。针对局部翘曲变形，朱继宏等[11]提出了多点保形的概念，用来衡量多点间相对位移变化程度和变形协调。基于保形指标的拓扑优化设计方法在考虑变形方向性、谐响应振动和几何非线性问题中得到进一步拓展[12-14]。

在工程问题中，有一类工程结构承受分布力载荷作用，例如受风载荷的建筑、承受气动力的飞行翼面、舵体等[3]。此类结构承载表面的节点变形量需要满足刚度设计要求。由于承载面节点数目庞大，并且在优化迭代过程中，最大位移发生的位置会发生变化，若采用多点位移约束的处理方式，不仅敏度求解中的伴随工况数量庞大，有限元分析计算量大，而且会造成约束方程数目增多而难以优化求解。

针对上述问题，提出考虑承载面最大位移约束的拓扑优化模型。遵循独立连续映射(Independent Continuous Mapping, ICM)建模方法[15-16]，引入倒变量中间变量，将原有拓扑优化问题转换为一系列严格正定的二次规划问题，并采用高效、稳健的序列二次规划算法求解。将提出模型获得的优化结果与传统的体积比约束下柔顺度最小化列式结果、多点位移约束下体积最小化结果做对比，说明提出列式与求解方法的可行性和优越性。

1 体积比约束下的拓扑优化列式

在现有的拓扑优化方法中，以体积比为约束、柔顺度最小化为目标的拓扑优化列式最为典型，数学表达式为[1]

find:ρ

minimize:c=FTu

(1)

式中:ρ为密度矢量;分量ρe为单元e的相对密度值，用于表征单元e的有无；c为系统柔顺度值；K、u和F分别为总刚度阵、位移列阵和载荷列阵；V和V0分别为优化后体积和初始结构体积；f为设定体积比；密度下限ρmin用于避免结构分析中的数值奇异性，这里取值ρmin=10-3；NE为设计域内单元总数。

SIMP插值模型通过建立弹性模量(或刚度阵)与密度变量的假设关系，驱使在优化过程中密度变量向0或1靠近，其表达式为[1]

(2)

柔顺度c对密度变量ρe的一阶导数为[1]

(3)

式中:ue为单元e的位移列阵。拓扑优化列式(1)在得到结构响应量敏度的前提下，将敏度数值代入到连续型优化算法，即可实现优化求解。在拓扑优化研究领域中，常见的优化求解算法有准则法(Optimality Criteria, OC)[17]、凸线性化(Convex Linearization, CONLIN)[18]、移动渐进线法(Method of Moving Asymptotes, MMA)[19-20]等。

2 考虑承载面最大位移约束的拓扑优化列式

2.1 拓扑优化列式

尽管列式(1)在拓扑优化研究领域得到了广泛的应用，但存在着以下缺点[1]：① 柔顺度值缺乏工程意义；② 体积比是人为给定的数值，但实际工程应用中，拓扑优化结果仅用于结构传力路径的确定，工程可用的结构需要进一步通过尺寸优化、形状优化等后续工作完成。

对于承受分布载荷的工程结构，承载面的变形量需要进行严格控制。假设承载面上每一个节点位移满足约束限值，以结构轻量化设计为目标的拓扑优化列式为

find:ρ

minimize:V/V0

(4)

为了解决上述困难，这里提出的解决方案是采用承载面上最大节点位移约束取代多点位移约束方程，即

find:ρ

minimize:V/V0

(5)

式(5)中的最大值函数具有不可导性，从而无法利用现有的连续型优化算法求解。受最大应力约束拓扑优化问题的启发[21-22]，此问题可采用包络函数包括KS函数、p范数(p-norm)和p平均(p-mean)函数对最大值函数进行连续光滑处理，这里选择采用KS函数，最大值约束的代理模型为

(6)

式中:μ为KS包络函数参数，并满足:

(7)

根据数值经验，当μ取值较大时，容易引起优化求解中的数值振荡；当μ取有限值时，包络函数值与真实最大值之间具有一定的差异性。为了消除这种差异性，引入修正系数以调整两者间的关系，即[23-24]

(8)

式中:cKS为调整系数，数学表达式为

(9)

2.2 敏度分析

对式(6)求偏导可得

(10)

由链式求导法则得

(11)

自由度j对应的位移表达为位移向量的函数，即

(12)

将式(12)代入到式(11)可得

(13)

对平衡方程Ku=F两边求导可得

(14)

假设载荷F与密度变量无关，则

(15)

将式(15)代入到式(13)中可得

(16)

建立伴随方程

(17)

式中：λ为伴随向量，则敏度表达式为

(18)

由式(18)可得，当采用提出的最大位移指标时，伴随方程对应一次有限元分析，相比较多点位移约束问题，可以大大减少伴随方程求解计算量。

2.3 独立连续映射法显式化过程

与常见拓扑优化方法不同，ICM以各类结构响应量为约束，重量最小化为目标。通常情况下，柔顺度、节点位移和应力等各类结构响应量难以求其二阶敏度，在获取一阶敏度信息的情况下，可获取其一阶近似表达。在引入倒变量情形下，目标体积函数可以解析获取对设计变量的二阶敏度值。综合目标函数和约束函数的近似表达，可构建出二次规划问题的模型，配合稳健高效的序列二次规划(Sequential Quadratic Programming, SQP)算法求解。本节将遵循ICM建模方法，实现目标函数和约束函数的显式化表达。

定义设计变量为密度变量ρe的倒数函数，即

(19)

由此推导体积函数对设计变量的一和二阶敏度表达式为

(20a)

(20b)

式中:ve为单元e的体积。

由式(20b)组成的海塞矩阵(Hessian)的表达式为

H=

(21)

与优化列式(4)对比可知，若直接采用ρe作为设计变量，体积函数对密度变量ρe的一阶敏度为常数，二阶敏度为0。当采用倒变量形式的设计变量xe时，由式(21)易得，海塞矩阵非对角线元素为0且对角线元素恒定大于0，即海塞矩阵严格正定且具有可分离性质。

利用链式求导法则可得

(22)

式中:λe为单元伴随向量。

基于上述敏度信息，位移约束函数采用一阶泰勒展开得到显式表达式：

(23)

式中:上标l代表第l轮优化迭代。

体积目标函数采用二阶泰勒展开得到显式表达式，并忽略其中的常数项可得

find:x

(24)

式中:矩阵B与体积函数的一阶敏度相关。优化列式(24)是标准的二次规划模型，为了克服拓扑优化中设计变量庞大的困难，通常转化为对偶模型并采用SQP求解[15-16, 25]，优化迭代至与体积比相对变化率满足收敛条件，即

|V(l)-V(l-1)|/V(l)≤ε

(25)

收敛率ε通常为1‰。

借鉴数字图像处理中的过滤技术消除拓扑优化中常见的棋盘格现象，具体过程可参考文献[15-16，24]。值得注意的是，由于过滤平均的影响，将导致结构边界存在着一定数量的灰度单元。

3 数值算例与讨论

本节采用数值算例来验证提出模型的可行性和有效性。在各算例中，材料的弹性模量和泊松比值分别为1 MPa和0.3，平面结构厚度为1 mm，单元尺寸为1 mm。初始结构的单元密度值取1，即设计区域全部为实体材料。

算例1如图1所示的长悬臂梁结构，结构尺寸为240 mm×84 mm，顶部4层单元为非设计区域，即单元密度值恒定为1。左端全部固定，顶部受到总载荷值为100 kN的力作用，载荷均匀分布在各个节点上。假设顶层y方向最大位移许用值为8 mm。KS函数参数μ=20。采用提出的拓扑优化列式和求解算法，所得到的拓扑优化结果如图2所示。

图1 长悬臂梁结构示意图Fig.1 Illustration of long cantilever structure

由图2可知，拓扑优化结果顶部的变形尽可能保持水平且上下起伏，这表明可能多处节点的位移达到约束限值。最终结构体积比为0.309，最大节点位移值为7.999 mm。

为了验证提出方法的有效性，采用拓扑优化列式(5)和MMA算法求解，拓扑优化结果如图3所示。

由图2和图3对比可知，2种不同方法下的拓扑优化结果及变形特点类似。MMA算法求解下的最优结构体积比为0.310，最大节点位移值为8.001 mm。上述结果也表明，本文提出的优化列式具有可行性。

图2 所提方法的拓扑优化结果Fig.2 Optimized topology obtained from the proposed method

图3 采用拓扑优化列式(5)和 MMA算法的拓扑结果Fig.3 Optimized topology obtained from Eq.(5) and MMA algorithm

为了对比与传统拓扑优化结果的异同，以体积比0.309为约束条件，采用刚度最大化列式(1)，得到的拓扑优化结果如图4所示，优化结构顶部受载区域最大位移为10.167 mm，比设定的许用位移值8 mm大27.1%。顶面结构变形保持为直线，朝着自由端倾斜，这与图2和图3中最优结构的变形特点明显不同。上述结果表明，传统的刚度最大化列式(1)不具有控制局部最大变形的能力。

图4 拓扑优化列式(1)下的拓扑优化结果Fig.4 Optimized topology using Eq.(1)

算例2如图5所示的平面结构，结构尺寸为180 mm×120 mm，底部4层单元为非设计区域，即单元密度值恒定为1,下端两角点全约束。下端面受到合力F=100 kN作用，载荷均匀分配到下端面每个节点上，假设底层y方向最大位移许用值为8 mm，最优拓扑与结构变形如图6所示。

图5 平面结构示意图Fig.5 Illustration of plane structure

图6 所提列式下的拓扑优化结果与变形Fig.6 Optimized topology obtained from the proposed equation

为了对比与多点位移约束下拓扑优化结果的异同，设置底部距离左端点1/6，2/6，3/6，4/6和5/6长度位置点位移约束条件，形成多点位移约束拓扑优化列式，同样采用ICM方法求解[15-16]，最优拓扑与结构变形如图7所示。位移约束设置位置对应的节点位移均达到许用值8 mm，但底部最大节点位移为9.433 mm，这说明即使采用多点位移约束，也无法严格控制承载面内每一点的变形。值得注意的是，多点位移约束的敏度分析包含5个伴随工况，即对应5次结构有限元计算，而本文提出的优化列式将多点位移约束凝聚为单约束函数，仅包含1个伴随工况。上述结果证明了提出方法的有效性和优越性。

图7 多点位移约束下的拓扑优化结果Fig.7 Optimized topology obtained from multiple nodal displacement constraints

算例3三维结构如图8所示，结构尺寸为90 mm×40 mm×28 mm，底部四角点全约束，顶部受到总载荷F=26 390 kN的力作用，载荷均匀分配到上端面每个节点上。顶部2层单元为非设计区域，即单元密度值恒定为1。设顶部区域y向最大许用位移为18 mm。采用提出的方法得到的拓扑优化结果如图9(a)所示，最优结构体积比为0.125，承载面最大节点位移为17.996 mm；设定0.125为体积比约束,采用柔顺度最小优化列式(1),拓扑优化结果如图9(b)所示，承载面最大节点位移为19.056 mm。

图8 三维结构示意图Fig.8 Illustration of 3D structure

由图9可知，不同的拓扑优化列式下对应的拓扑优化结果有所不同。类似于二维结构，本文提出方法对对该三维结构具有控制承载面最大位移的效果。即在相同的材料用量下，采用柔顺度最小化目标，优化结构的节点最大位移值超过许用位移值约5.89%。算例结果证明了提出列式和优化求解算法在三维结构优化问题中的适用性。

图9 不同列式下的拓扑优化结果Fig.9 Optimized topology obtained from different equations

4 结 论

针对工程中受分布载荷的连续体结构，提出承载面最大位移约束下的拓扑优化设计模型。该模型通过KS包络函数实现了多点位移最大值的凝聚，推导了相应的敏度表达式。基于ICM方法，形成了具有正定可分离海塞矩阵的二次规划列式，采用序列二次规划法优化求解，得到以下结论：

1) 与传统的体积比约束下柔顺度最小化列式对比，提出方法具有控制局部区域最大位移的效果。优化结果证明了提出方法的有效性。

2) 与多点位移约束对比，提出方法更精确、有效控制承载面的最大变形量。且将多点位移约束凝聚为单一约束，不仅减少了伴随工况计算工作量，且有利于优化求解。

3) 提出的优化列式和相应求解算法有望在几何非线性位移约束问题中进一步拓展。