一种自由度缩减和收敛加速的高效等几何拓扑优化方法

杨雨豪 郑 伟 王英俊，2

1.华南理工大学聚合物新型成型装备国家工程研究中心，广州，510641

2.华中科技大学数字制造装备与技术国家重点实验室，武汉，430074

0 引言

拓扑优化是一种在给定载荷与约束的条件下寻找最优材料分布使结构目标性能最佳的结构优化方法[1]。拓扑优化设计可在设计域内自由改变拓扑结构，具有高设计自由度，目前已成为结构设计、增材制造领域的研究热点，并广泛应用于船舶机械、航空航天、军工及建筑土木等领域[2-3]。常见的传统拓扑优化方法有变密度法[4]、水平集方法(level set method，LSM)[5]、渐进优化法(evolutionary structural optimization，ESO)[6]、移动变形组件方法(moving morphable component，MMC)[7]和等几何拓扑优化[8]等。其中，等几何拓扑优化是一种基于等几何分析(isogeometric analysis，IGA)[9]的新型结构拓扑优化方法。与其他拓扑优化方法相比，等几何拓扑优化将CAD、CAE与结构拓扑优化无缝融合，为实现结构的设计、分析、优化的一体化奠定了理论基础。

近年来，等几何拓扑优化得到了广泛的研究与发展。WANG等[8]和GAO等[10]分别对等几何拓扑优化进行了较为全面的综述。GAO等[11]提出了一种具有密度分布增强函数的等几何拓扑优化方法来显示结构的拓扑构型，该方法可用于解决一些数值优化问题(如多材料结构和机械材料的合理设计等)。在通用代码方面，GAO等[12]公开了一套高效、紧凑的等几何拓扑优化MATLAB代码，为学习者提供了便利。

虽然等几何拓扑优化可避免网格离散，且具有高精度、高效率等优点，但优化中的迭代过程仍需要消耗大量时间。随着工业的发展，实际工程应用中的优化问题变得日益复杂，导致优化问题的计算规模也相应增大，因此，提高计算效率成为相关研究人员重点关注的内容。目前主要从以下两个方面提高拓扑优化的计算效率：一是从硬件角度，引入CPU/GPU并行计算[13-14]以提高运算速度；二是从高效算法角度。从算法层面提高优化计算效率可从根本上减小拓扑优化的时间复杂度[15]，因此，本文主要从算法角度出发，开展等几何拓扑优化高效算法的相关研究。

拓扑优化的求解过程需要反复迭代，而每次迭代都需完成一次结构性能分析，因此，提高拓扑优化的效率可从减小计算量和加快优化收敛速度两方面开展研究。在拓扑优化过程中，有限元分析是耗时占比较大的一个步骤，如果能够缩减有限元方程的自由度，可有效减小计算量，且自由度规模的减小也有利于优化迭代收敛。近年来，国内外学者提出了一系列用于缩减拓扑优化中有限元方程自由度的方法。WANG等[16]实现了基于多层网格的等几何拓扑优化有限元求解，有效减小了有限元方程的规模；杜义贤等[17]将序列插值模型和多重网格方法应用于拓扑优化，缩减了有限元求解的自由度；GOGU[18]提出了一种动态构造简化基的方法，有效缩短有限元求解步骤的耗时。

在拓扑优化过程中，迭代收敛时所需的次数和时间也是影响整体优化时间的重要因素。为了加速拓扑优化问题的收敛，可在优化迭代过程中缩减部分设计变量以加速迭代过程。KIM等[19]提出一种设计变量可缩减方法用于减少拓扑优化迭代次数，有效缩短迭代收敛所需时间。LIAO等[20]提出了一种局部更新策略，只更新实体区域及其邻近边界区域，在减小更新设计变量所需计算量的同时减少了优化迭代次数。除了减少设计变量以加快收敛速度外，还有一些方法可以加速迭代收敛，如廉睿超等[21]提出一种灰度单元双重惩罚方法，加速了中间密度单元的二值化；许小奎等[22]提出了两种形式的对比度增强算子，驱动中间密度向两级进行分化。此外，拓扑优化也可以和机器学习等热门方法相结合来提高计算效率[23-24]。

本文从性能分析方程的自由度缩减和优化收敛加速两个方面提高等几何拓扑优化的计算效率。

1 理论基础

1.1 NURBS基函数

等几何分析将NURBS引入CAE描述分析物理场，实现了CAD建模与CAE分析的统一表达。在拓扑优化领域，用NURBS基函数插值优化变量场，可进一步将CAD建模、CAE分析与结构拓扑优化纳入统一框架，有望实现结构设计、分析与优化的一体化。

通常使用递归方式来定义B样条的基函数。Ξ=(ξ1，ξ2，…，ξn+p+1)是一个参数空间中的非递减实数数列，被称作节点向量，其中，n是控制点和基函数的数量，p是样条的次数。由Cox-de Boor递推公式[25]定义p次B样条基函数：

(1)

其中，Bi，p(ξ)为第i个p次B样条基函数。向每个B样条基函数引入一个权重ωi，便可得到单变量的NURBS基函数：

(2)

根据张量积性质，可以构造二维或多维的NURBS基函数Ni，p(ξ)如下：

(3)

p=(p1，p2，…，pdp)

式中，ξ为参数坐标；i为张量积结构中的索引位置；p为所有参数方向的多项式次数；Nim，pm为沿参数方向m的B样条基函数，可由式(2)计算得到；dp为参数空间的维度；Ξ(m)为该参数空间内的节点向量，表示参数坐标；nm为对应的函数个数；pm为沿参数方向m的多项式次数。

1.2 基于SIMP方法的拓扑优化

固体各向同性材料惩罚(SIMP)法是密度法拓扑优化中最为常用的材料插值方法。在SIMP方法中，每个单元被赋予一个介于0、1之间的单元密度，单元密度xe和弹性模量Ee之间的关系如下：

(4)

式中，Emin为弹性模量最小值(Emin>0)，给空单元赋予该最小值可防止刚度矩阵发生奇异；E0为单元密度为1的实心单元的弹性模量；β为惩罚系数，用于实现中间密度的二值化，通常取β= 3。

在经典最小柔度目标优化问题中，优化问题的数学模型如下：

(5)

式中，C(x)为目标函数，此处表示结构柔度(是实际柔度值的两倍)；x为设计变量；n为设计域中设计变量的个数；U为总体位移向量；F为总体受力(载荷)向量；K为总体刚度矩阵；ue为单元e的位移向量；ke为单元e的单元刚度矩阵；V(x)为优化过程中的总体积；V0为设计域的总体积；VF为体积比约束(VF∈[0，1])。

最优化准则方法(OC)是求解式(5)优化问题的一种高效方法。在求解过程中，设计变量x可以不断迭代，并在迭代过程中逐步逼近最优解，其迭代过程依赖于下式所示的算法：

(6)

(7)

式中，M为最大步长即设计变量每次迭代的最大变化量(M>0)；η为数值衰减系数(取η=1/2)；λ为拉格朗日算子，用于在迭代过程中保证体积约束，通常可由二分法计算获得[26]。

目标函数C和体积V对xe的导数分别为

(8)

(9)

其中，式(9)建立的前提是每个单元都有单位体积。

1.3 基于控制点的等几何拓扑优化

与传统的SIMP方法不同，等几何拓扑优化的数值计算是在控制点上进行的，因此，基于等参的思想，在等几何拓扑优化中参数坐标为ξ的设计变量x可由附近的控制点求出：

(10)

式中，Ni为影响该参数坐标的第i个控制点的NURBS基函数；xi为控制点密度。

在基于密度法的等几何拓扑优化中，单元的刚度矩阵ke可由下式计算：

(11)

在等几何拓扑优化中，控制点密度x(ic)可由下式计算：

(12)

式中，xi为设计域中编号为i的控制点的密度值；ic对应编号为i的单元的中心；ci为影响单元i的所有控制点；cij为影响单元i的第j个控制点；Ncij(ic)为控制点cij的NURBS基函数。

单元i及对应的控制点关系如图1所示，在次数为2的二维控制点网格中，单元i受到周围9个控制点(图中红色控制点)的影响。

图1 单元i及其控制点的对应关系示意图

因此，传统SIMP法中求解灵敏度的式(8)可以改写为

(13)

式中，ei为被控制点i影响的所有单元；eij为ei中的第j个单元。

一个关于控制点i及所控制的单元的位置关系示意图见图2，图中以二维控制点网络为例，红色控制点i的影响范围是它周围的9个单元。

图2 控制点i及其所控制的单元示意图

2 高效等几何拓扑优化方法

2.1 自由度缩减方法

拓扑优化中的关键一步是求解式(5)中的方程KU=F，并得到位移向量U。本节介绍两种自由度缩减方法以加速等几何拓扑优化中的控制点位移求解，分别是基于位移变化的自由度缩减和基于空单元的自由度缩减。

2.1.1基于位移变化的自由度缩减

以单元规模为64×32的悬臂梁等几何拓扑优化为例，观察优化过程中位移向量的变化情况，如图 3所示，发现大部分自由度对应的位移在迭代后期几乎保持不变。由此可知，一旦某自由度的位移变化趋于稳定后，它在后续的迭代过程中几乎保持不变。因此，将保持位移不变的自由度进行缩减，能减小性能分析方程的求解规模，从而加快优化进程。

图3 随机选取的自由度位移随迭代次数变化曲线

基于位移变化的自由度缩减可分为两个步骤：①判断各个自由度的位移是否已经保持稳定；②对位移已稳定的自由度进行缩减。

首先，本文提出了一种多次迭代内自由度位移相对变化的量化表示：

(14)

式中，D为在过去的几次迭代中位移的相对变化量；l为当前迭代次数；U(l)为当前迭代的位移向量；L=5，即式(14)可以衡量过去10次迭代位移变化的波动情况。

然后，通过对比D与设置的阈值ε的大小，即可判断某个自由度对应的位移是否已经保持稳定，即当D(i)<ε时，可认为在i处的位移已经保持稳定，可进行自由度缩减，反之不能。通常情况下，随着迭代的进行，可缩减的自由度会逐渐增加。

一旦判断某自由度处的位移在数次迭代处于稳定状态，即可在方程求解过程中移除该自由度。下面以小规模(10×10)方程为例，阐明迭代过程中基于位移变化缩减自由度的过程：

(15)

式中，Ki，j、Ui、Fi(i，j=1，2，…，10)分别为总体刚度矩阵中第i行第j列的元素、自由度i的位移和自由度i的外载荷。

假设通过式(14)判断U9在最近的数次迭代中基本保持不变，对上面的分析方程进行自由度缩减，可得到缩减规模(9×9)的方程：

(16)

基于位移变化的自由度缩减流程中，部分自由度位移被直接代入下一次迭代计算中，可能会引入一定误差。随着迭代的进行，累计误差可能会大到足以影响等几何拓扑优化结果，因此本文引入一种误差的量化表示方法：

(17)

式中，er为等几何分析求解误差；K、F为自由度缩减之前的全局刚度矩阵和力向量；U*为规模缩减后等式的求解结果位移向量。

当er大于设定的误差阈值时，下次迭代将不再进行基于位移变化的自由度缩减，以保证等几何分析求解的误差在设定阈值内。随机选取的自由度位移随迭代次数的变化曲线见图3。

2.1.2基于空单元的自由度缩减

本方法的首要问题是判断哪些自由度在等几何分析中是可以被缩减的。以二维的等几何拓扑优化为例，图4所示为一个由12个矩形单元、30个控制点和60个自由度描述(前30个为各控制点横向自由度，后30个为各控制点纵向自由度)的等几何拓扑优化设计域。图4中，每个单元的编号ei写在单元内部上方，每个控制点的编号写在控制点的上方，括号中是该控制点对应的两个自由度编号。该设计域对应的等几何分析方程为

图4 由12个矩形单元、30个控制点和60个自由度描述的二维等几何设计域

(18)

基于SIMP法的等几何拓扑优化中，控制点的密度会随着迭代的进行向0或1靠近。假设在某次迭代之后，图4中影响e1、e2、e5和e6单元的所有控制点的控制点密度均更新为0，如图5所示。如果忽略赋给空单元的非常小的弹性模量Emin，那么单元e1、e2、e5和e6的单元刚度矩阵可视为零矩阵，控制点1影响的所有单元均为空单元。根据总体刚度矩阵的组装方法可推出，整体刚度矩阵中控制点1所对应的第1、31行和1、31列的元素均为0，同理，控制点2、7、8对应的2、32、7、37、8、38自由度所在行与列的元素也为0。仅考虑1、31自由度对应的等几何分析方程情况，可以将式(18)写为

图5 可缩减控制点示意图

(19)

鉴于总体刚度矩阵的第1、31行和1、31列均为零元素，可将这些行、列移除，同时移除U和F中的对应元素，从而得到结果与原方程相同的规模缩减方程。虽然缩减方程无法求解U1和U31的值，但是，由式(8)可知，无论控制点的位移如何，空单元e1、e2、e5和e6的单元灵敏度数值都为0，即U1和U31的数值不会影响优化求解结果，无需计算。同理，图5中的U2、U32、U7、U37、U8和U38也无需计算。因此，当一个控制点影响的所有单元均为空单元时，该控制点对应的自由度可以缩减。随着迭代的进行，将出现更多的空单元，即有更多的自由度可以缩减。

2.2 收敛加速方法

2.2.1基于设计变量缩减的收敛加速

受拓扑优化加速算法研究[16，19-20]的启发，本文引入了一种基于设计变量缩减的收敛加速方法，具体的收敛策略是：在等几何拓扑优化迭代中，当某设计变量连续多次迭代保持基本不变时，可以停止更新该设计变量，以加速收敛。

以单元规模为64×32的二维悬臂梁等几何拓扑优化为例，随机选取的几个控制点处的控制点密度随迭代次数的变化曲线见图6。从图6中可以看出，大多数控制点密度在接近0或1后即保持基本不变，因此在之后的迭代中可以不更新这些已保持稳定的设计变量。

图6 二维等几何拓扑优化在64×32规模下的几个控制点密度变化曲线

判断控制点密度是否需要更新的标准如下：

max(x(k)，x(k-1)，…，x(k-Q+2)，x(k-Q+1))<0.01

或

min(x(k)，x(k-1)，…，x(k-Q+2)，x(k-Q+1))>0.99

(20)

式中，x(k)为第k次迭代的设计变量值；Q为一个整数，本文将Q设置为5。

式(20)的判断条件表示，当某控制点的密度连续Q次迭代保持大于0.99或小于0.01时，即可在之后的迭代中不更新该变量。另外，对于不再更新的变量，将控制点密度大于0.99的赋值为1，而将控制点密度小于0.01的赋值为0，以加速收敛。

由图6可发现，索引为755和777的设计变量可能在迭代早期尚未稳定，为保证基于设计变量缩减的收敛加速方法的准确性，应在迭代的中后期使用该收敛加速方法。可以利用目标函数值的相对变化量来判断是否处在迭代早期，计算公式如下：

(21)

式中，cobj为目标函数的相对变化量；c(k)为第k次迭代的目标函数值；N为一个整数，本文将N设置为5。

当目标函数的相对变化量cobj大于设定的阈值δ时，可以认为优化过程仍处于迭代早期。只有当cobj<δ时，才可以进行基于设计变量缩减的迭代加速。

2.2.2基于灰度抑制的收敛加速

一个单元规模为128×64的二维悬臂梁等几何拓扑优化结果的控制点密度分布如图7所示。可见，即便使用了密度惩罚方法，最终的优化结果中仍有许多中间密度存在，这将影响等几何拓扑优化问题的收敛速度。

图7 控制点密度分布

针对此问题，参考GROENWOLD等[27]的研究，本文引入了基于灰度抑制的设计变量更新方法，在最优化准则方法的基础上继续进行灰度抑制。将sigmoid函数向x轴正向平移0.5单位长度可得

(22)

图8 不同参数a对应的平移后的sigmoid函数图像

本文中，迭代中的参数a由下式求得：

a=t/cobj

(23)

其中，t是一个常数。在等几何拓扑优化中，随着迭代的进行，cobj将逐渐减小，a逐渐增大，即加快设计变量二值化。如果常数t设置得太小，则加速不明显；如果t设置得过大，则可能影响到优化结果的准确性。经测试，本文将t设置为6。由于拓扑优化迭代早期的拓扑构型还不稳定，应根据式(21)判断何时引入基于灰度抑制的收敛加速方法。

3 多重加速流程的算法实现

本文阐述了两种等几何拓扑优化加速方法，分别是自由度缩减方法和迭代收敛加速方法，在等几何拓扑优化中同时使用多种加速方法可以获得显著的加速效果，具体的算法流程图见图9。图9中，加速流程从等几何拓扑优化进入迭代开始，到迭代收敛结束。图中黄色虚线框中是基于空单元的自由度缩减算法流程，蓝色虚线框中是基于位移变化的自由度缩减算法流程，绿色框中是基于设计变量缩减和灰度抑制的收敛加速算法。流程图中，L0是预设的迭代早期迭代次数，在经过至少L0次迭代之后开始进行基于位移变化的自由度缩减和收敛加速，一般设置L0>10。

图9 高效算法流程图

4 算例

4.1 二维悬臂梁等几何拓扑优化算例

二维悬臂梁是经典拓扑优化算例，该优化问题的设计域、边界条件和外加载荷如图10所示。本算例的优化目标是在保证体积约束的条件下，求得柔度最小的拓扑构型。本文在三个不同控制点网格规模下，分别对比加速前、自由度缩减后、收敛加速后和整体高效方法的优化时间与优化结果，进而研究各个算法的加速效果。

图10 二维悬臂梁优化问题的设计域、约束和外加载荷

本例中，悬臂梁的体积比约束VF=0.5，设计域划分为三种网格规模：128×64、256×128、512×256。二维悬臂梁等几何拓扑优化问题对应的控制点网格如图11所示，以10×5的单元网格(12×7的控制点网格)为例，控制点用绿色表示。

图11 二维悬臂梁等几何拓扑优化设计域的控制点网格

4.1.1基于自由度缩减的加速

本例中，自由度缩减算法的参数设置如下：误差阈值t= 0.04，位移变化阈值ε=0.0004。对比三种网格规模下等几何拓扑优化的等几何分析求解时间，前100次迭代的等几何分析求解时间-迭代次数曲线如图12所示。在迭代早期两条曲线基本重合，当引入基于自由度缩减的加速方法后，求解时间明显缩短，而且随着网格规模的增大，加速效果也愈加显著。图中蓝色曲线在迭代后期出现一些波动，这是由于触发了误差阈值限制，当误差值大于设定的阈值时，下一次迭代会取消基于单元位移的自由度缩减加速。可以通过调整误差阈值t来调整迭代后期的加速效率和时间曲线。

(a)网格规模128×64

由图12可知，基于自由度缩减的加速方法可以有效加速等几何分析求解速度，且随着迭代的进行，加速前后的求解时间差逐渐增大，说明基于自由度缩减的加速方法对大规模的等几何拓扑优化问题尤为有效。需要注意的是，基于自由度缩减的加速方法需要一些前处理，包括识别求解空单元对应的自由度、需要缩减的位移等，但是相较于自由度缩减所带来的加速效果，这些前处理时间的影响可以忽略。

4.1.2迭代收敛加速

本节通过对比使用迭代收敛加速算法前后，二维悬臂梁等几何拓扑优化问题的迭代次数、目标函数值等来讨论迭代收敛加速方法的效果。在所有算例中，优化问题的收敛条件均设置为：连续两次迭代的控制点密度的最大差值小于0.01。

迭代收敛加速的参数设置如下：目标函数改变的阈值δ=0.01，灰度抑制算法的参数a=6/cobj。在三种网格规模下，算例在收敛加速前后的迭代次数、目标函数值的对比如表1所示。由表1可知，迭代收敛加速算法可以显著加速等几何拓扑优化进程，减少大量的迭代次数，减小冗余计算量。得益于灰度抑制算法，加速后结果的柔度通常优于加速前的优化结果。

表1 二维悬臂梁等几何拓扑优化问题在收敛加速前后的迭代次数、柔度值对比

4.1.3高效等几何拓扑优化

综合利用基于自由度缩减的加速方法和迭代收敛加速方法，即可在保证结果准确、稳定的基础上获得良好的加速效果。本节通过二维悬臂梁算例对比高效等几何拓扑优化方法与传统等几何拓扑优化方法，来验证高效等几何拓扑优化方法的高效性。

高效等几何拓扑优化算法的参数设置与前述相同，以128×64、256×128、512×256三种网格规模划分设计域。图13所示为不同规模下使用高效等几何拓扑优化算法前后的优化结果模型，可以看出，使用加速算法前后，拓扑构型基本相同，说明提出的高效等几何拓扑优化具有可靠性。同时，得益于灰度抑制算法对中间单元的缩减，加速后的优化结果边界更加清晰、明显。

(a)128×64加速前 (b)128×64加速后

不同规模下高效算法应用前后的迭代次数、柔度和总优化时间对比如表2所示。可见，在基于自由度缩减的加速算法和收敛加速算法的作用下，二维悬臂梁的等几何拓扑优化求解速度得到了较大提高。总体优化的加速比分别为1.56、2.85和6.02，证明了该方法具有较好的加速效果。此外，前两个规模下加速后的柔度值下降，第三个规模下的柔度值仅增大了0.76%，处于可接受范围内，说明加速方法得到的目标函数依然准确。

表2 二维悬臂梁算例高效等几何拓扑优化加速效率对比

综上所述，高效等几何拓扑优化方法可以在保证优化结果准确的情况下大幅加速等几何拓扑优化问题的求解和迭代过程。

4.2 三维MBB梁等几何拓扑优化算例

MBB梁也是拓扑优化的经典基准算例，取它一半的设计域、边界条件和外加载荷如图14a所示。图中，左侧面的X、Z方向受位移约束，右下方的活动铰支座表示其Y、Z方向受位移约束。以体积比约束VF= 0.3的三维MBB梁等几何拓扑优化为例，验证自由度缩减和收敛加速的高效等几何拓扑优化方法在三维问题中的效果。

三维MBB梁等几何拓扑优化的目标是在保证体积约束的情况下获得柔度最小的拓扑结构，设计域对应的控制点网格如图14b所示，为了展示方便，图中的网格规模被缩小为16×9×4，绿色为控制点，红色为单元网格线。

(a)设计域与边界条件

本节在三种网格规模上(30×10×4，60×20×6和90×30×8)对比验证加速前后迭代次数、优化结果和优化时间。本例中，高效等几何拓扑优化方法的参数设置为：自由度缩减的误差阈值t= 0.04、位移变化阈值ε=0.005、目标函数改变的阈值δ=0.3、灰度抑制算法的参数a=6/cobj。如表3所示，加速后的优化迭代次数明显减少，目标函数柔度值也有所减小，整体迭代优化时间明显缩短，不同规模下的加速比分别为2.03、2.69和5.64。理论上，高效等几何拓扑优化算法在更大的网格规模下会获得更好的加速效果。

表3 三维MBB梁算例高效等几何拓扑优化加速效率对比

不同规模下加速前后的三维MBB梁等几何拓扑优化结果构型如图15所示。可见，加速前后的结果拓扑构型基本相同，因此高效等几何拓扑优化算法在加速计算的同时不影响最终的优化构型。由于引入了灰度抑制算法，使得图15右侧的加速后结果构型的中间灰度单元更少、边界更加清晰。

(a)30×10×4加速前 (b)30×10×4加速后

综上所述，本文所述的高效等几何拓扑优化方法在三维算例中依然有着较好表现，可大幅减少优化迭代次数、缩短优化时间，且随着问题规模的扩大，加速效果愈加显著。

5 结论

本文提出了一种高效等几何拓扑优化方法，包括两种自由度缩减方法和两种迭代收敛加速方法。其中，基于自由度缩减的加速方法通过缩减等几何分析方程的规模显著缩短每次迭代中等几何分析求解时间，迭代收敛加速方法通过加快设计变量的二值化，显著减少优化迭代次数。本文以二维悬臂梁和三维MBB梁两个经典优化问题为例，分别对比了迭代次数、目标函数值和迭代时间，验证了高效等几何拓扑优化的加速效果及其对优化结果的影响。结果表明，本文提出的高效等几何拓扑优化方法可以显著加速等几何拓扑优化，计算时间成本大幅缩减，对难以收敛的优化问题和大规模问题的加速效果尤其显著，验证了本文提出方法的有效性。另外，高效拓扑优化方法在处理柔度最小问题时采用了灰度抑制算法，减少了中间灰度单元的数量，可获得更加清晰的拓扑构型。将来，该方法可以与其他拓扑优化方法相结合，并推广到热、振动、流体等其他领域结构优化问题的高效求解。