结构拓扑优化的一些基本概念和研究方法

结构优化的目的是使结构在满足一定使用功能要求下尽可能降低成本或取得最佳性能。这里的成本目前主要是材料, 实际上也可以包含建造、维护、运行, 甚至可以包含拆除和破坏损失等全生存周期成本。结构优化方法已经广泛应用于许多工程领域, 近几十年得到了迅速发展。目前许多大型有限元分析软件系统也都嵌入了结构优化设计模块。

结构优化分为尺寸优化, 形状优化和拓扑优化。结构尺寸优化和形状优化有明确的优化变量, 如横截面尺寸、板厚、结点坐标等, 可以容易地直接写出优化问题的数学表达式, 采用常规的数学优化方法求解。

结构拓扑优化研究连续体中开孔数量、形状和位置, 或杆系结构的结点数量、位置以及结点之间的杆件连接关系。结构拓扑优化问题一般很难直接找到明确的优化变量, 这是结构拓扑优化区别于尺寸优化和形状优化的最重要特征和困难所在。由于结构潜在拓扑形式可以有无限多, 难以用有限个变量描述, 采用数值优化方法求解时, 难以直接写出具体数学优化表达式。

结构拓扑优化早期研究工作主要是寻找拓扑优化结构的解析表达式, 代表性的工作是米歇尔理论, 将在第1 节介绍。但是, 由于其求解 和应用的困难, 后来的研究工作主要集中在数值求解方法, 将拓扑优化问题转化为0-1优化问题或尺寸和形状优化问题。如何选择优化 变量间接描述结构拓扑, 如何将结构拓扑优化问题转化为0-1优化问题或尺寸和形状优化问题是结构拓扑优化的最根本问题。各种结构 拓扑优化数值方法的本质区别就是这种转化方案的不同, 这些方法将在第2节介绍。前面所述的优化目标结构是连续体, 而第3节介绍 的是杆系结构。无论连续体还是杆系结构都会遇到一些数值优化不稳定等共同问题, 这些在第4节介绍。作为数值优化方法, 最后总 要将结构优化问题转化为数学优化问题。第5节介绍了与结构优化问题相关的数学优化求解方法。关于结构拓扑优化更详尽的内容可 以参考近几年的综述文献[1, 2, 3]。

Michell[4]最早研究了应力约束体积最小问题解析解。他首先证明了单工况应力约束体积最小结构具有的满应力性质, 后来被称为米歇尔准则。利用这个准则, 采用一种逆解法, 可以求出几个简单的拓扑优化结构[5, 6, 7, 8, 9, 10, 11, 12, 13, 14]。这些结构后来被统称为米歇尔桁架。米歇尔桁架不是一般意义上的桁架, 而是一种非均匀各向异性连续体。这种结构由杆件构成, 杆件只受轴力作用, 符合桁架的性质; 但是这些杆件之间无间隙地连续分布, 所以是连续体。这种类似桁架的连续体结构称为"类桁架连续体"。

与之对应的加肋板优化结构称作"类格栅"[15, 16, 17, 18, 19, 20]。如果加强肋承担了主要载荷, 忽略中面的内力, 则成为梁系优化问题。这两个问题的优化结果完全一样。另一类对应的空间问题称为"类拱格栅"[21, 22], 或"Prager结构"。他研究竖向载荷作用下的空间结构优化问题。竖向载荷作用在结构上, 作用点由水平面的投影位置给定。在结构优化过程中, 作用点的高度随结构一起变化。

米歇尔的研究工作具有里程碑的意义.

他发现了拓扑优化结构理论上一般是非均匀各向异性连续体, 而不是常规意义上的桁架。

例如, 图1(a)中的结构下边两端简支, 中间有集中力作用; 图1(b)给出的拓扑优化结构是由连续均匀分布的径向杆和沿圆周分布的集中杆件构成的类桁架连续体$^{[5]}$。图1(b)中的径向线只是示意性地画出连续分布杆件其中的一部分。采用后面将要介绍的各种数值优化方法可以得到如图1(c) -图1(f)所示离散桁架结构$^{[23]}$或拓扑等价的带孔板。图1中百分数是相对解析解的体积 增加量。可以看出, \结构拓扑越复杂, 结果越接近理论解。带有大量灰色区域的解析解是最优解。米歇尔理论的另一个重要意义在 于：他提供了对各种数值拓扑优化结果的评判标准, 解析解是拓扑优化结构材料体积目标函数下限。

这里的解析方法是一种逆解法, 不具有一般性, 而且仅限于单工况应力约束体积最小问题。非均匀各向异性材料也不便于工程应用。 从工程应用角度考虑, 人们更倾向于得到各向同性均匀连续体。所以, 近几十年拓扑优化领域主要研究工作集中在带孔各向同性均匀 连续体的数值拓扑优化方法。

带孔的各向同性均匀连续体拓扑优化方法的核心问题是如何描述孔。因为拓扑优化问题中孔的数量也需要优化, 所以在结构拓扑优化中很难事先布置一定数量的孔, 然后去优化孔的形状和位置。

由于结构优化一般都采用有限元方法分析结构, 所以用有限元直接表示孔最方便。首先将设计域划分为有限单元, 通过在任意位 置删除部分有限单元 可以近似表示任意数量和形状的孔。

第1种方法最直接, 就是将单元的"有"和"无"作为一个[0, 1]取值的离散优化变量。Xie等[24, 25]提出的进化结构优化方法, 每次迭代删除一定比例(称为放弃比)效率较低的单元形成空洞。迭代过程中放弃比逐步改变一个量(称为进化比)。后来又发展出的双向进化结构优化方法[26]还可以增加单元, 由于其方法简单有效而被广泛接受。

第2种采用放松的方法, 将[0, 1]取值的离散变量变为[0-1]的连续变量。Bendsoe等[27]提出的均匀化方法在每个单元内嵌入一个矩形孔, 每个孔的长和宽以及角度作为优化设计变量。当孔充满单元时则单元为空; 反之, 当孔大小为零时保留单元。带孔单元的刚度矩阵需要数值方法计算, 其值实际是单元内部均匀化的结果, 这正是这个名称的来历。对于不同尺寸的孔需要事先计算出各种单元刚度矩阵。后来提出的多阶材料模型假设材料沿几个不同方向分布不同密度的微结构。多阶材料模型可以是二维或三维微结构, 不仅可以描述更复杂的材料, 而且避免了单元刚度矩阵的复杂计算[27, 28]。

放松的方法会得到大量中间密度的灰色区域。虽然带有灰色区域的结构更容易接近理论解, 理论上可能更优, 但往往不是工程上期望的结果。为此, 通过罚中间密度的方法抑制中间密度。优化方法采用这样一种人工材料[29]

$$E({\pmb x}) = \rho ({\pmb x})^pE_0 (1)$$

式中$\rho ({\pmb x})$是单元材料密度,${\pmb x}$表示位置, $E_{0}$和$E$是材料的弹性模量和计算时实际采用的弹性模量值, 指数$p> 1$是罚参数。罚参数的作用是在密度变小时单元刚度变小的速度更快, 使得中间密度会趋于零, 有效抑制了中间密度。这就形成了采用罚函数的各向同性实体或微结构优化方法(solid isotropic material or microstructure with penalization, SIMP)。有1个99行的小程序可以帮助大家理解这个方\linebreak 法[30]。

隋允康提出独立连续拓扑变量及映射变换方法, 针对式(1)的关系研究了几个不同方案。针对几何尺寸、允许应力、刚度等几个不同性质的物理或几何量采用不同的独立的映射关系[31, 32]。该方法可以用于连续体和刚架结构的优化。

水平集方法 定义一个称为水平集函数的标量函数$\varPhi({\pmb x})$。结构拓扑用这个函数表示

$$\varPhi ({\pmb x})\left\{ \begin{array}{ll} { \geqslant 0}\, ,\ \ & \hbox{实体}\\ { = 0}\, , & \hbox{边界}\\ { < 0} \, , & \hbox{空} \end{array} \right. (2)$$

式中${\pmb x}$表示优化设计变量。通过改变这个函数来改变域的内外边界实现结构拓扑优化[33, 34]。产生新孔是水平集方法的一个比较困难的问题。泡泡法 通过一个特征函数决定在适当的位置结构中插入孔[35]。几何描述方法中改变边界的能力是十分重要的性能。

Guo等[36]和Zhang等[37]提出可移动变形组件或孔洞法。他们事先在设计域内布置几个弹性组件或孔洞, 通过优化弹性组件或孔洞的形状和尺寸实现结构拓扑优化。相对基于单元描述或基于结点的描述方法显著减少了设计变量, 提高了优化效率, 减少了后处理工作。

Bendsoe 等[38]提出的自由材料优化方法 将材料的弹性张量的所有分量都作为优化设计变量, 通过优化弹性张量的分布场函数实现结构拓扑优化。这种材料具有最大的优化空间, 但与真实材料不对应, 难以实现, 一般需要进一步的后处理[39, 40, 41, 42, 43]。

从第1章介绍可知, 理论上, 拓扑优化结构是非均匀各向异性类桁架连续体, 而本章前2节数值优化方法旨在得到工程上经常需要的均匀各向同性带孔连续体, 这实际上是在最优解附近寻找近优解。这类方法试图将理论上最优和工程上适用两个相互矛盾的目标放在一起考虑, 引起一系列困难：为了得到等厚板限制了中间密度, 限制中间密度又导致了单元铰接, 为此又要控制梯度或周长等。解决这些矛盾的另一个思路是将这两个相互矛盾的目标放在两个独立的过程中实现：先实现理论上的最优, 再根据工程上的需要离散化, 这就是基于类桁架材料模型的优化方法[44, 45, 46]。

基于类桁架材料模型的优化方法根据拓扑优化的类桁架性质建立一种类桁架材料模型。采用有限元分析方法和数学规划方法可以得到十分接近理论解的拓扑优化类桁架结构, 这种类桁架材料结构与杆件具有明确的对应关系。可以保留部分离散杆件成为桁架结构[46]; 也可以取舍部分区域的连续分布杆件成为带孔连续体。这个过程称为离散化过程。

上一章是研究连续体, 本章讨论桁架和刚架结构, 统称杆系结构。杆系结构拓扑优化方法目前主要是基结构方法：在设计域内布置足够多的结点, 所有结点之间都连接杆件。通过优化这些杆件的横截面面积, 删除横截面面积过小的杆件实现拓扑优化。基结构方法完全回避了拓扑变量的选择问题。假设按照$n$行$n$列规则布置共$n^{2}$个结点, 所有结点之间都连接杆件, 在平面问题中共$n^{2}(n^{2} -1)/2$个杆件, 在空间问题中共$n^{3}(n^{3}-1)/2$个杆件。杆件数量会随着结点网格数量增加很快。如果结点在各个方向都按照$n =100$布置, 平面问题是近5$\times$10$^{7}$个杆, 空间问题是5$\times$10$^{11}$个杆。这个规模的优化问题太大了, 不适于解决实际问题, 目前主要用来检验其他数值优化结果。即使基结构中布置了很多杆件, 仍有可能漏掉最优的杆件。结合结构形状优化方法, 甚至建立杆件与结点的增加策略是一个可行的方法。如果能够根据对初始结构的初步分析结果合理布置杆件也可以显著减少杆件数量。

杆系结构和连续体结构优化研究工作在一定程度上具有一致性。当连续体材料占据的空间相对于设计域空间的比例趋于零时, 拓扑优化的连续体趋于桁架。同时, 他们也存在一些共性问题。

在应力约束优化问题中, 当拉、压允许应力不相等时, 可能会出现奇异性问题。因为当杆件截面趋近于零时, 计算应力并 不趋近于零, 存在所谓极限应力, 表现为可行域中的各部分的维度不一致, 数学优化过程的强非线性、非凸性。这些问题使得寻找全局最优解变得非常困难。

为了便于理解, 这里举一个非常简单的例子, 如图2所示一个两端固定的杆件, 中间受到一个轴向集中力。为了简单起见, 假设杆长为2, 拉、压允许应力分别为1和2, 集中力的大小为1。两段杆件的横截面积分别计作$x_{1}$和$x_{2}$, 可行域由$D_1 = \{(0, x_2 )\vert x_2 \geqslant 1 / 2\}$和$D_2 = \{(x_1 , x_2)\vert x_1 + x_2 \geqslant1\}$共同构成。前者$D_{1}$是1维空间, 后者$D_{2}$是2维空间。优化点(0, 1/2)在1维空间的端点, 远离$D_{2}$空间。如果初始点选在$D_{2}$, 则常规的数值迭代方法很难通过$D_{1}$达到优化点。

Cheng等[47, 48, 49, 50]将常用的应力约束改写成内力约束形式, 并做适当放松, 使退化的可行子空间被扩充, 问题得到很大程度的解决, 取得重要进展。奇异性问题不仅存在于桁架, 也存在于连续体; 不仅存在于应力约束, 也存在稳定等约束问题[51]。

但是, 一般地, 单工况下允许应力为常数时不存在奇异性问题。

许多结构拓扑优化方法研究柔度(外力功)最小化问题。这种优化问题的约束条件只有一个, 而且对局部参数不敏感, 称为全局约束条件。但是, 像强度这样最普遍的应力约束条件, 要求结构在所有位置都满足, 对局部参数敏感。离散化的数值优化方法需要满足大量的约束条件, 过多的约束条件会增加优化求解的难度。为了减少约束条件, 可以将这些大量局部约束写成一个全局整体约束形式[51, 52]

$$G_{ks} = \dfrac{1}{p}\ln \sum_i {\rm e}^{pR_i } \leqslant \sigma _{ s} (3)$$

或

$$G_{kk} =\Big [\dfrac{1}{p}\sum_i (R_i )^p \Big ]^{1 / p} \leqslant \sigma _{ s} (4)$$

式中

$$R_i = \dfrac{f_i (\sigma )}{f_{\max } (\sigma )}\, f_{\max } (\sigma ) = \mathop {\max }\limits_i f_i (\sigma ) (5)$$

$f_i (\sigma )$是$i$点等效应力, $\sigma _{ s}$是允许应力, $p$是1个较大的数, 累加号表示对所有约束位置累加。较大的$p$可以增加最大应力的权重, 减少局部应力超限的机会, 但也会更容易引起迭代震荡。

以SIMP为代表的拓扑优化方法为了得到清晰的(0-1)结构, 避免大面积的灰色区域使用式(1)的罚函数方法抑制中间密度, 抑制中间密度导致棋盘格、单元铰接以及局部极小值等数值不稳定问题。

为了帮助理解这个问题, 讨论图3(a)所示矩形域(假设无限宽, 忽略边界影响)受均布力作用。假设设计域内允许应力为1, 均布外力集度为0.5, 上边为刚性不可设计域。应力约束优化结构是厚度为0.5的板。但是, 限制厚度为0或1时, 如果如图3(b)和图3(c)所示只有一半材料是同样的优化结构, 优化结果不唯一。实际优化计算中具体得到哪一个结构很大程度上取决于有限单元数量。Cheng等[53]最早在研究加肋板优化问题时发现了这种单元网格依赖性问题。

如果采用4结点矩形单元, 由于位移模式的双线性假设, 图3(d)的结构也 同样最优, 这就是著名的棋盘格现象。同理, 图1(b)的连续均匀径向分布杆在抑制中间密度的情况下也会出现图3(b)和图3(c)的单元依赖性和图3(d)的棋盘格现象。提高有限元的精度可以减轻棋盘格现象, 但是, 在优化后期一些区域单元数量逐渐减少。由于铰接区域的单元数过少导致计算误差过大, 需要借助周长控制、局部梯度控制、过滤等其他方法[54, 55]。

工程结构实际工作环境经常具有大量不确定性。在确定条件下的最优结构在不确定性环境下可能就不再是最优了, 因此结构拓扑优化有必要考虑这些不确定性因素。不确定性因素可以与结构有关或无关, 前者的分析难度更大。可以采用鲁棒性设计优化和基于可靠性的设计优化[56, 57], 可以基于概率的和非概率的方法[58]。

考虑不确定性的优化方法远比确定性优化方法困难, 而且优化结果也不一样。如图4所示, 集中力作用点$O$点距上边固定边距离给定的情况下, 如果确定载荷的方向在$ - 45^\circ \leqslant \theta \leqslant 45^\circ$灰色区域内, 优化结构的杆件与竖直方向的夹角$\alpha = 45^\circ $; 如果是在同样范围内的不确定载荷, 则$\alpha =35.3^\circ $[59]。这是由于在不确定载荷下, 确定的结构需要在所有不确定载荷下都安全; 而确定性载荷下的结构只需满足一个确定载荷条件。因此, 不能用确定性载荷下的拓扑优化结果直接评判不确定载荷的拓扑优化结果。

结构拓扑优化方法不仅可以实现结构优化, 还可以进行材料微结构优化形成轻质高强材料。为了能够进一步提高结构整体性能, 同时实现宏观结构和材料微结构的多尺度优化则是一个具有吸引力的研究方向。当然, 这种优化需要两个层次的嵌套, 运算量巨大, 可制造性也是一个问题。采用周期性微结构可以增加实用性、降低计算成本[60]。微结构的材料还在散热性能方面表现出一些优良性能[61]。

前面介绍的各种结构拓扑优化数值方法最终都要求解数学优化问题。求解数学优化问题主要有准则法, 序列规划方法和启发式算法。

应用最早的是满应力准则的应力比公式

$$\rho _i^{k + 1} = \rho _i^k \sigma _i^k / \sigma _{ s} (6)$$

式中$\sigma _i^k $是工作应力, $\rho _i^k$是作为优化设计变量的材料参数, 上角标是迭代指标, 下角标表示位置。满应力准则用于应力约束优化问题具有很高的计算效率, 对优化设计变量不敏感, 特别适合大规模拓扑优化问题。在单工况应力约束下体积最小结构符合满应力准则, 他同时也是体积约束柔度最小结构。静定结构甚至可以一次迭代达到最优解。对于多工况或多约束的超静定结构, 拓扑优化结构并不满足满应力准则。通过构造拉格朗日函数利用极值条件也可以导出一些特定问题的优化准则, 用于结构拓扑优化分析。

优化问题的目标函数和约束函数如果都是线性函数, 则属于线性规划问题, 可以采用成熟的单纯形法求解, 通过有限次迭代得到最优解。对于非线性规划问题, 可以利用泰勒级数将目标函数和约束函数在当前点近似展开成线性规划或二次规划子问题, 用序列线性规划或序列二次规划 求解。如果将函数展开为变量的倒数形式可以改善计算性能[62]。

Fleury[63]和Zhang等[64]将目标函数和约束函数用凸线性化方法展开。Svanberg[65]提出移动渐近线方法, 移动渐近线方法后来发展出了广义移动渐近线方法, 基于梯度移动渐近线方法 和全局收敛移动渐近线方法 等多个版本$^{[66]}$。这些方法的展开式保证了优化子问题的凸性, 有较好的收敛性, 目前在结构优化领域取得良好的效果。

系列序列规划方法以及一些准则法都需要计算目标函数和约束函数的导数(敏度)。基于敏度的方法一般具有较快的收敛速度和较高的优化效率, 适合较大规模的结构拓扑优化, 特别是利用二阶导数会进一步提高优化收敛速度。但是, 导数计算本身的计算量往往很大, 特别是二阶导数计算量更大。因此, 一般只使用1阶导数, 或再加上2阶导数海森矩阵的对角元。

敏度分析有3个基本方法[67]：虚载荷法, 状态空间法和设计空间法, 他们针对的表达式不尽相同。优化问题本身如果是非凸问题, 基于导数的方法有可能会陷入局部最优解。

遗传算法 借鉴遗传学的研究成果, 基于概率的搜索方法, 通过复制、交叉和变异几种操作逐步寻找优化点[68, 69, 70, 71]。 模拟退化算法[72]等都是属于启发式算法, 这些方法都不需要导数信息。虽然在理论上, 这些方法经过无限次地迭代就可以接近1的概率找到全局最优解, 但实际运算时也可能出现"早熟"等问题, 并不一定收敛到最优解。而且这些算法的收敛速度都很慢, 计算效率不高, 不适于大规模的寻优, 不推荐直接用于拓扑优化问题的求解。

本文对结构拓扑优化领域一些基本概念进行了简单评述, 包括力学模型, 数学优化方法以及其中的一些重要问题。本文主要从结构力学角度进行了总结。实际上, 结构拓扑优化方法目前已广泛应用于许多工程领域, 大量研究工作集中在了多学科优化问题, 包括航空、材料、微机电、热力学、流体力学、断裂力学等诸多领域, 本文都没有涉及。但本文讨论的内容是这些研究的重要基础。

周克民，华侨大学土木工程学院教授，博士生导师，1985年毕业于天津大学应用力学专业，后获得天津大学固体力学硕士和结构工程博士学位。著有《类桁架材料模型及结构拓扑优化》《有限元法与程序设计》等。发表结构优化方面的论文数十篇，多项软件著作权。主持国家自然科学基金3 项。担任福建省力学学会副理事长，《力学与实践》编委。

The authors have declared that no competing interests exist.