基于动态高斯滤波的拓扑优化灵敏度过滤方法

王 伟,张岐良+,徐颖珊

(1.中山大学 系统科学与工程学院，广东 广州 510006；2.北京空天技术研究所，北京 100039)

0 引言

结构优化[1]是在给定的设计区域内，在满足约束条件和载荷作用下，实现结构的最佳性能。随着航空航天业以及建筑业的快速发展，越来越多的工程师和建筑师开始注重结构优化，都试图利用结构优化技术来得到最完美的结构。现如今，结构优化已渗入到工程中的各个方面，例如可用于多相材料结构设计[2-3]、柔性结构设计[4]、热结构优化[3,5-7]、结构轻量化[8]等。

结构优化分为尺寸优化、形状优化和拓扑优化，其中拓扑优化是目前公认最具有难度、最有挑战性的难题。根据拓扑结构表示模型的不同，拓扑优化方法可分为两个分支[9-10]：材料描述模型和边界描述模型。材料描述模型主要包括变密度法[11-13]和ESO(evolutionary structural optimization)法[14]，它是通过有限个密度单元离散设计域；边界描述模型主要包括水平集法(Level Set Method，LSM)[15]、移动可变形组件(Moving Morphable Components,MMC)方法[16]和移动可变形孔洞(Moving Morphable Void，MMV)方法[17]，它是用高维函数隐式或显式地表示结构边界的演化过程。

变密度法因原理简单、易于实现等优点被广泛使用，但数值不稳定现象如网格依赖性、棋盘格等常出现在拓扑结构中，为了解决此类现象，在文献[12]中，SIGMUND提出了传统的灵敏度过滤方法。传统的灵敏度过滤方法考虑过滤半径内所有单元对中心单元的影响，通过引入权重因子，对过滤半径内所有单元灵敏度进行距离加权平均，来修改中心单元的灵敏度。通过灵敏度过滤方法，达到解决网格依赖性、棋盘格等现象的目的，但由于这种“加权平均”效果使拓扑优化边界上出现中间密度单元，造成拓扑优化边界模糊，出现所谓的“边界扩散”现象。

为解决因灵敏度过滤造成的“边界扩散”现象，廉睿超等[18]对灰度单元进行分层双重惩罚，即对过滤后的单元灵敏度再进行惩罚；龙凯等[19]考虑密度梯度的灵敏度过滤方法，即在原有灵敏度过滤表达式中增加密度梯度权函数项；张志飞等[20]利用双重固体各向同性惩罚微结构模型(Solid Isotropic Microstructures with Penalization，SIMP)，首先利用带灵敏度过滤的SIMP方法进行拓扑优化，然后用不带灵敏度过滤的SIMP方法，从而有效抑制灰度单元的产生；陈垂福等[21]考虑变过滤半径的灵敏度方法，通过在优化过程中逐渐减小过滤半径来抑制中间密度单元的出现；匡兵等[22]对灵敏度过滤公式进行改进，考虑了过滤前中心单元灵敏度的影响；李家春等[23]基于空间扩张策略建立了双向插值函数的拓扑优化模型；SIGMUND[24]提出密度投影方法，利用Heaviside函数对中间密度投影来获得0/1解决方案；FU等[25]提出一种边缘光滑的材料分配策略，利用基于tanh函数表达式的Heaviside平滑函数对网格点密度进行投影，使密度值快速趋向0和1两端；HUANG[26]提出一种浮动投影拓扑优化(Floating Projection Topology Optimization，FPTO)方法，该方法基于tanh 函数利用参数调节对设计变量进行浮动投影。

当前方法存在优化迭代步数多、优化参数难以确定等问题，且优化结果不稳定，最终的拓扑结构仍可能出现中间密度单元。针对具有体积约束的结构柔度最小优化问题，采用变密度法来求解该优化问题，其主要步骤包括：建立模型、有限元计算、灵敏度过滤与设计变量更新。针对Sigmund的传统灵敏度过滤方法因“加权平均”造成的“边界扩散”现象，本文提出一种基于动态高斯滤波的灵敏度过滤方法，即在优化过程中逐渐减小高斯函数标准差的值，降低周围单元对中心单元的影响，对单元灵敏度进行动态滤波，以避免在优化中产生中间密度单元，消除边界模糊现象，并通过经典数值算例验证了该灵敏度过滤方法的可行性和稳定性。

1 变密度法拓扑建模与求解方法

变密度法[11-13]将设计域的材料分布问题视为各单元的密度值问题，引入连续设计变量x作为单元的相对密度值，通过SIMP模型材料属性的合理近似模型(Rational Approximation of Material Properties，RAMP)引入惩罚因子，并建立材料弹性模量与单元相对密度之间的函数关系，对中间密度值的单元进行有限度的惩罚，尽量减少中间密度单元数，使单元密度值尽可能趋向0或1。

1.1 密度插值模型

本文采用SIMP模型，建立的材料弹性模量与单元相对密度之间的函数关系如下[13]：

(1)

式中:E0为固体材料的弹性模量；Emin为空洞部分的弹性模量，其值很小，通常取E0/109，以此来防止整体刚度矩阵奇异；每个单元e被赋予一个密度xe；Ee为单元e的弹性模量；p为惩罚因子。

基于SIMP模型的弹性模量插值公式，结构的整体刚度矩阵表达式以及整体刚度矩阵灵敏度信息表达式分别如下[13]：

(2)

(3)

1.2 基于SIMP的拓扑优化模型

基于密度插值模型，建立以结构柔度最小为目标，具有体积约束的优化问题[13]：

minc(x)=UTKU=

s.t.

KU=F,

0≤x≤1。

(4)

式中:F为施加的载荷矢量，c为结构柔度，U为整体位移矢量，K为整体刚度矩阵，ue是单元位移矢量，ke为单元刚度矩阵，x为设计变量，N为用于离散设计域的单元数目，p为惩罚因子，V(x)为材料体积，V0为设计域体积，f为规定的体积分数。

1.3 有限元分析与计算

在结构设计领域，有限元分析是一种常用的近似数值计算方法，它是用若干个有限大小的单元体离散化连续体。有限元计算主要包括以下步骤：

(1)将连续体变换为离散体。将连续体离散成有限个单元体的集合。

(2)对单元进行分析。求出单元的结点信息、结点力信息以及刚度矩阵等。

(3)进行整体分析。利用单元刚度矩阵组装整体刚度矩阵，并建立结点平衡方程组，求解出单元各结点的位移大小。

在变密度法中，有限元求解是其中的一个步骤，通过调用有限元求解器更新结构的柔度值与单元灵敏度，进而更新设计变量的值，然后根据新的设计变量值再次调用有限元更新变量，直至达到拓扑优化的收敛标准。

1.4 优化模型求解方法

以结构柔度最小为目标，具有体积约束的拓扑优化问题本质上是一个二次规划问题，而求解这样的二次规划问题主要有两种方法：一种是优化准则法(Optimality Criteria，OC)[1]，一种是移动渐进线法(Method of Moving Asymptotes,MMA)[27]。其中OC法主要用来解决单约束问题，MMA法主要用来解决多约束问题，本文建立的模型为单约束柔度最小化问题，因此采用OC法求解更为方便，效率也更高。

OC法[1]通过引入拉格朗日乘子，将带约束的目标函数转换成拉格朗日函数，如下所示：

L=c+λ1(V(x)-fV0)+λ2(KU-F)+

(5)

(6)

拓扑优化迭代终止条件由设计变量的最大变化值决定, 当变化量小于等于设定的容许误差值时，拓扑优化达到收敛标准，从而跳出循环，数学表达式如下[13]：

max(max|xk+1-xk|)≤ε。

(7)

式中:k为迭代步数，xk+1为第(k+1)步设计变量的值，xk为第k步设计变量的值，ε为设定的容许误差。

2 灵敏度过滤方法

基于本文的拓扑优化模型，目标函数c对单元密度xe的灵敏度为[13]：

(8)

2.1 传统灵敏度过滤方法

Sigmund的传统灵敏度过滤方法(简称传统方法)[13]如下：

(9)

Hei=max(0,rmin-dist(i,e))，

(10)

Ne={i|dist(i,e)}≤rmin。

(11)

式中:dist(i,e)为单元i的中心与单元e的中心之间的距离，rmin为过滤器的过滤半径。

传统方法中，考虑了过滤半径内所有单元对中心单元的影响，通过引入权重因子对过滤半径内的所有单元灵敏度进行距离加权平均，来修改中心单元的灵敏度。过滤半径越大，对中心单元造成影响的周围单元越多，消除棋盘格和网格依赖性效果越好，但是过大的过滤半径会导致拓扑优化边界的不清晰，出现“边界扩散”现象；过小的过滤半径虽然能减少灰度单元的产生，但容易出现棋盘格和网格依赖性现象。因此，传统方法不能在避免棋盘格和网格依赖性的基础上，消除“边界扩散”现象，这样会影响拓扑结构的可制造性。

2.2 动态高斯滤波的灵敏度过滤方法

在拓扑优化初始时，传统方法利用过滤半径内的单元对中心单元的距离加权平均影响，来消除棋盘格和网格依赖性现象；随着拓扑优化的进行，拓扑结构慢慢形成，这时棋盘格和网格依赖性现象已经基本消除，但由于过滤半径的存在，中心单元仍会受到周围单元的影响，出现“边界扩散”现象。为此，本文提出一种基于动态高斯滤波的灵敏度过滤方法，在拓扑优化迭代的初始阶段，为避免出现棋盘格和网格依赖性现象，将高斯函数的标准差设定为最大值，考虑过滤半径内其他单元对中心单元的影响；随着优化的进行，逐渐减小标准差的值，防止过滤方法的“加权平均”效果产生灰度单元，从而避免出现“边界扩散”现象。其过滤公式如下：

(12)

对于高斯核函数来说，标准差σ越大，高斯核函数W的平滑程度就越好，过滤半径内的其他单元对中心单元的影响就越大；标准差σ越小，高斯核函数W越尖锐，滤波效果就越弱。其中，均值为0，不同σ值对应的高斯核函数图像如图1所示。

当标准差σ小到一定程度时，高斯滤波就相当于对单元灵敏度没有进行过滤。因此，在拓扑优化初始时，将标准差σ设定为最大值，考虑过滤半径内其他单元对中心单元的影响；随着优化的进行，逐渐减小标准差σ的值，避免出现“边界扩散”的现象。

为了衡量最终拓扑结构中单元密度的离散程度，构建拓扑优化过程的离散性能指标[19]：

(13)

式中:s为结构的离散率，N为设计域中所有离散单元的数目，xi是第i个单元的相对密度。s为0～1之间的数，s越小，则更多单元的密度趋向0与1两端，中间密度单元越少，拓扑结构越清晰；当s为0时，拓扑结构构型完全为离散0-1矩阵的形式。

高斯函数的标准差σ在拓扑优化过程中是动态变化的，本文利用离散率s对高斯参数σ进行更新，公式如下：

σ=σmin+(σmax-σmin)×s。

(14)

式中:σmin为标准差σ的最小值，σmax为标准差σ的最大值，通过该公式更新标准差σ，可使σ从最大值σmax逐渐变化到最小值σmin。

基于动态高斯灵敏度过滤的结构拓扑优化具体步骤如图2所示。

3 数值算例

将本文动态高斯滤波方法与变密度法结合，并用数值算例来说明基于动态高斯滤波的灵敏度过滤方法(简称本文方法)的可行性与稳定性。在所有算例中，单位都是统一的，并假定以下参数：E0=1，Emin=10-9，材料的泊松比ν=0.3，rmin=1.5，ε=0.01，本文优化算法程序在MATLAB R2012a上实现[13]，实验环境为：CPU:11th Gen Intel Core i7-11700，GPU:NVIDIA GeForce GTX 1650 SUPER和32.0 GB运行内存。

为了评价最终优化结构的清晰程度，量化结构中灰度单元的数目，引入结构灰度率这一性能指标，其公式如下：

(15)

式中:g为结构的灰度率，xmin取0.001，h(xi)为拓扑结构中单元密度在xmin～1之间的单元数目，N为设计域中单元的总数目。灰度率g越大，优化结果存在的灰度单元越多；灰度率g越小，灰度单元越少，优化结构边界越清晰。

3.1 过滤参数分析

如图3所示，在简支梁上部中间位置施加竖直向下，大小为1 kN的集中力，将设计域划分为60×20的网格，以结构柔度最小为目标，目标体积为50%建立拓扑优化模型，因为简支梁的结构和所受载荷均对称，所以取简支梁的一半，并设置对称约束，进行结构拓扑优化。

σmin和σmax取值的不同组合会影响拓扑优化结果，为了评估其影响，本开展了参数研究。以简支梁模型为例，选取惩罚因子p=4，对参数σmin与σmax进行参数选优，建立实验表如表1所示。拓扑优化最终数据如表2和图4～图7所示，最终结构如图8所示，依据迭代步数、柔度值、结构离散率和灰度率来确定合适的参数组合值。

表1 不同参数组合值对应的实验表

表2 不同参数组合值简支梁优化数据

由表2和图8可知，本文方法得到的拓扑结构离散率s的值很小，结构中单元密度趋于两极化，且灰度率g很小，拓扑结构边界清晰，在消除棋盘格现象的同时，可以解决“边界扩散”问题，最终的拓扑结构几乎无灰度单元。由图4可知，σmin位于0.01～0.1之间时，优化迭代步数较少，最少迭代步数仅为39次，σmax增大时迭代步数会有所增加，但继续减小σmin值无助于减少迭代步数；当σmin取0.5或1时迭代步数大幅增加，最大迭代步数为200次，是最少迭代步数的5倍左右。由图5知，当σmin一定时，σmax的值越大，相应的结构柔度值越小；当σmax一定时，随着σmin减小，结构柔度值整体会有一个下降的趋势，且σmin=0.1是一个分界线，σmin≤0.1时的结构柔度值与σmin>0.1时的结构柔度值相差较大。由图6可知，当σmax一定时，σmin越小，结构离散率s的值越小，单元密度越趋于两极化；当σmin≤0.1时，s的值处于一个很小的范围，最小s值仅为4.6807e-13。由图7可知，当σmin在0.01～0.1之间取值时，结构的灰度率g极小，最小值为0，拓扑结构中几乎无灰度单元，结构边界清晰；当σmin在0.5～1之间取值时，灰度率会相应增大，结构边界处会出现灰度单元，最大灰度率达到2.75%。由图6和图7可看出，σmax越大，结构相应的离散率和灰度率会越大。同时考虑优化的迭代步数、结构柔度、结构的离散率与灰度率，数值实验表明σmin取值范围应在0.01～0.1区间内。

在简支梁算例中，当σmin分别取0.01、0.05、0.1，图8的(1)～(5)、(6)～(10)、(11)～(15)分别对应不同σmax值的拓扑结构。当σmin分别取0.01与0.05，σmax取5和10时，拓扑结构右下角有不连通区域；当σmin取0.1，σmax取5时，结构右下角有不连通区域，σmax取25时，结构内部有微小孔洞区域，对制造工业技术要求较高。数值实验表明，当σmax取值过大，结构柔度值会随之降低，但拓扑结构内部会出现微小孔洞区域；当σmax取值过小，拓扑结构会出现不连通区域；在充分考虑拓扑结构可制造性的前提下，σmax取值应在15～20之间。

由图8可知，当σmin取值范围在0.01～1，σmax取值范围在5～25时，得到的拓扑结构相似，说明过滤参数的取值对拓扑结构构型的影响小，从而证明了本文方法的适用性与有效性。

3.2 简支梁优化算例

依据上述参数分析，选择σmin=0.01，σmax=20，选取简支梁的一半结构(如图3)进行拓扑优化。为了检验本文方法的有效性与高效性，将本文方法与传统方法以及其他文献提出的方法进行对比实验。实验相关数据和最终拓扑结构如表3所示。

由表3可知，传统方法的结构灰度率达到了34.2500%，拓扑结构边界模糊，存在“边界”扩散现象。本文方法以及文献[19-21]、文献[24]提出的方法均降低了结构的柔度值，对灰度单元起到了抑制作用，所得到的最终拓扑结构与传统方法相似。文献[19]提出的方法对灰度单元抑制的效果不太理想，灰度率g的值为12.2500%，结构内部存在着灰度单元。文献[20]和文献[24]得到的拓扑结构内部存在微小孔洞区域，限制了结构的可制造性，且文献[24]的优化迭代次数达到了422次，是传统方法的4.5倍左右。文献20的方法和本文方法的结构灰度率g均为0，结构中几乎无灰度单元，消除了由于“加权平均”效果导致的在拓扑优化边界上出现中间密度单元造成的结构边界模糊，但本文方法的迭代次数仅为45次，是所有方法中最少的，比文献[21]的迭代次数减少了32.8358%。本文方法得到的结构离散率最低，离散率的值仅为4.6807e-11%，结构内部单元密度趋向0～1两端；本文方法得到的结构柔度值为195.9681，略高于其他方法的结构柔度值，但相对于传统方法降低了3.5554%。综上可知，本文方法在减少迭代次数和抑制“边界扩散”现象上都体现出了一定的优势。

表3 不同方法的简支梁优化结果对比

图9和图10分别是采用不同方法得到的结构体柔度与体积比随迭代步数变化的曲线。

如图9所示，本文方法和文献[20]提供的结构优化方法在迭代步数为21步左右，结构柔度开始收敛，早于其他方法，可见本文方法可以加快结构柔度的收敛。如图10所示，采用本文方法在确保得到目标体积的条件下，在优化后期结构体积收敛于目标体积，而其他方法得到的结构体积一直在目标体积附近小幅波动，本文方法相对于其他方法具有更强的稳定性。

3.3 L型梁优化算例

L型梁的尺寸、边界条件与受力情况如图11所示。L型梁的顶部被完全约束，右侧上角位置处受到竖直向下，大小为1 kN的集中力，以结构柔度最小为目标，目标体积为30%建立拓扑优化模型，根据参数选优，选取σmin=0.01，σmax=20，设计区域用90×90个离散单元进行划分，分别用不同方法对L型梁进行拓扑优化，最终拓扑相关数据和拓扑结构如表4所示。

由表4可知，传统方法的结构灰度率为12.6420%，结构边界模糊，存在中间密度单元。本文方法结构的灰度率为1.2346e-2%，结构中无明显灰度单元，拓扑结构边界清晰，无“边界扩散”现象。文献[19]、文献[21]与本文方法的迭代次数接近，但本文方法的结构柔度值比文献[19]和文献[21]的结构柔度值低，且拓扑结构无明显的微小孔洞区域，文献[21]的拓扑结构右上角有明显的孔洞区域，而文献[19]的结构柔度值较大，为138.759 9。文献[20]和文献[24]得到的结构柔度值都较低，相对传统方法分别降低了4.0160%和4.5222%，但文献[20]所得结构存在孔洞区域，文献[24]的优化迭代次数最多，为374次，优化收敛速度慢。综上可知，本文方法不仅可以避免出现“边界扩散”现象，还可以加快拓扑优化的收敛速度。

表4 不同方法的L型梁优化结果对比

3.4 悬臂梁优化算例

为了研究本文方法的网格依赖性问题，以悬臂梁为例，如图12所示。悬臂梁左侧被完全固定，右侧中间位置受到竖直向下，大小为1 kN的集中力，以结构柔度最小为目标，目标体积为50%建立拓扑优化模型，研究不同网格数量下的拓扑结构，并与传统方法进行对照。选取参数σmin=0.01，σmax=20，两种方法的最终拓扑结构如图13和图14所示。

图13是不同网格下采用传统方法得到的拓扑结构，可以看出拓扑结构出现了一定的网格依赖性现象。当网格划分越细，拓扑结构会出现越多的细支结构，结构变得越来越复杂，这会加大结构实际制造的难度、成本与时间，不利于实际工程制造。由图14可知，不同网格下用本文方法得到的拓扑结构基本未出现网格依赖性、棋盘格等数值不稳定现象，网格划分的越细，拓扑结构边界越光滑，更有利于实际工程的制造；图13的各个拓扑结构都有中间密度单元，造成拓扑边界不清晰，而图14的各个拓扑结构基本没有中间密度单元，边界清晰。由此可知本文方法不仅无棋盘格、网格依赖性等数值不稳定现象，还解决了拓扑优化结构的“边界扩散”现象。

4 结束语

本文提供一种基于动态高斯滤波的灵敏度过滤方法。随着拓扑优化的进行，减小高斯函数标准差的值，防止“加权平均”效果产生中间密度单元，避免了“边界扩散”现象，提高了拓扑优化结果的可行性和稳定性，同时更快、更容易达到收敛。通过对不同参数组合值进行实验，确定了参数值的最佳取值范围，并通过经典算例证实了本文方法不仅保持了传统方法的优点，还解决了拓扑结构的“边界扩散”现象。

在实际工程中，结构的中间密度单元并不存在，而本文方法可以有效避免中间密度单元的产生，得到具有清晰边界的拓扑优化结果，也使得最终的拓扑结构更容易导入CAD软件中进行后处理，故本方法对在实际工程中应用拓扑优化具有现实意义。但本方法对过滤参数的选取具有一定的局限性，如何自动化地进行参数选优将是后续研究的重点。