基于变密度法的分区加权敏度过滤方法*

高 涵,徐 雷,胡元昊,郭战岭

(四川大学 机械工程学院,四川 成都 610065)

0 引 言

作为一种有效的结构优化方法,结构拓扑优化方法可在设计初期创建最佳设计,找到合理的材料分布方式,获得更坚固、更轻盈的优化结构,可以减少所需的材料,大幅降低成本。

与尺寸优化和形状优化相比,拓扑优化拥有更多的设计自由度,可以克服结构参数化的局限性,获得更大的设计区域,是结构优化领域的一大研究热点[1]。目前,常见的拓扑优化方法主要有变密度法、均匀化法、渐进结构优化法和水平集法等[2]。

在拓扑优化的结果中,由于数值不稳定而导致结构[3]的不可制造,经常会出现边界扩散、棋盘格、网格依赖性和细长杆结构等。为解决拓扑优化过程中所出现的数值不稳定等问题,国内外学者进行了深入的研究。

SIGMUND O[4]采用引入最小过滤半径进行距离加权平均的敏度过滤方法,有效解决了棋盘格、网格依赖性等问题,但该方法仍存在优化效率低和边界扩散等问题。YIN F等人[5]建立了一种以质量和位移作为目标函数和约束,基于概率可靠性的高效拓扑优化模型,可以快速地获取清晰的拓扑优化边界;虽然其迭代速度较快,但个别区域仍存在边界扩散的问题。为了获得可制造性高的优化结构,ZHOU M等人[6]提出了一种在拓扑优化中达到最小长度尺寸的方法;但由于施加最小长度尺度不一定能保证稳定的优化性能,会导致出现一些灰度单元。DING M等人[7]采用基于过滤—投影的结构参数化方法,完成了对结构变化的精确控制,使中间密度单元比例不断降低;但该方法的迭代次数较多,效率较低。廉睿超和张国锋等人[8，9]提出了一种考虑分区的敏度过滤方法,可有效抑制边界扩散现象;但由于该方法的优化效率较低,同时可能出现细长杆或多孔结构,降低了结构的可制造性。

为此,针对Sigmund敏度过滤方法中边界扩散和优化效率低等问题,笔者提出一种基于变密度法的分区加权敏度过滤方法,将原敏度过滤区域一分为二,确定加权因子修正系数、惩罚因子和最小过滤半径的优化参数,引入加权因子对不同敏度过滤区域进行加权处理,使靠近中心的区域获得更高的敏度值,远离中心的区域获得较低的敏度值,以获取边界清晰的拓扑优化结构。

1 变密度法优化模型与求解算法

1.1 变密度法拓扑优化模型

拓扑优化是一种以给定的设计变量、约束条件和受力情况为依据,以获取更大的刚度、更轻的重量为优化目标,通过较少的先验决策提高结构优化效率,在给定区域内寻找材料的最优分布方式的轻量化设计的结构优化方法[10]。

结构优化可分为3种,即形状优化、尺寸优化和拓扑优化[11]。工程应用中运用最广泛的拓扑优化方法是变密度法,它以柔顺度最小作为优化目标进行拓扑优化[12],其数值表达式以相对密度作为设计变量,材料的弹性模量和密度满足的数学表达式为:

(1)

式中:Ei—第i个单元插值后的弹性模量,Pa;E0—固体部分材料的弹性模量,Pa;Emin—空洞部分材料的弹性模量,Pa;ρi—单元相对密度;P—惩罚因子。

其中:

(2)

拓扑优化过程中,笔者设定体积约束,以柔度值达到最小为优化目标,确定材料单元密度和材料属性二者的函数关系,使材料属性能够以材料单元密度函数的形式表达。

其数学模型为:

(3)

式中:ρ—单元相对密度的向量;C(ρ)—给定拓扑的最小柔度值;U—单元节点的结构位移向量位移;K—结构刚度矩阵;N—单元数目;ui—第i个元素的位移列向量;F—元素节点的施加载荷向量;V(ρ)—优化后的结构体积;f—预先设定的体积分数;V0—设计域的体积;vi—第i个元素的单元体积;ρi—第i个元素的相对密度值;ρmin—包含最低允许相对密度的向量。

1.2 拓扑优化的优化准则法

为了获得理想的拓扑优化结果,需要合适的数值求解算法。目前应用较广的拓扑优化数值求解算法有优化准则法、数学规划法和随机搜索法等[13]。优化准则法是通过将正的灵敏度值所对应的设计变量置零,保留负的灵敏度值,可以处理相应约束的最优准则,并据此建立优化迭代数学模型[14]。

优化准则算法的主要实现方式是通过联立目标函数和预设约束条件得到Lagrange方程,将约束条件与目标函数结合成无约束问题。

目标函数的优化准则法数学模型为:

(4)

式中:λ1,λ2,λ3,λ4—拉格朗日乘子;ρi—元素i的设计变量;ρmin—设计变量的取值下限;ρmax—设计变量的取值上限。

优化准则法迭代收敛快、计算量较小,不随结构复杂度及设计变量的增多而复杂化,在简单约束条件下的结构优化中得到了广泛应用[15],所以笔者采用优化准则法。

其设计变量迭代公式可表示为:

(5)

式中:ρe—单元e的相对密度;t—平移限度;η—阻尼系数;Be—中间变量;ξ—拉格朗日乘子。

其中:

(6)

2 敏度过滤

2.1 Sigmund敏度过滤方法

在基于变密度法的拓扑优化中,经常会出现网格依赖性和棋盘格等数值不稳定的现象,因此,Sigmund敏度过滤方法被提了出来。其主要步骤是先确定一个中心单元,然后设定一个最小过滤半径为rmin的敏度过滤区域,再引入加权因子,对各单元到中心单元的距离进行加权平均处理,使靠近中心单元的各个单元获得较高的敏度值;而远离中心单元处于敏度过滤边界的各个单元获得较低的敏度值,从而使过滤半径范围内各单元敏度的加权平均值代替中心单元的敏度值;从中心单元到边界的敏度值大体呈下降趋势,可防止出现局部设计区域内的单元密度值剧烈波动,有效抑制网格依赖性和棋盘格现象的出现。

过滤后的单元敏度的数学表达式为:

(7)

由式(7)可知:加权因子是根据单元e和单元i之间的距离进行线性取值,加权因子从中心单元到过滤边界所在单元呈线性递减趋势。这样虽然可以保证靠近中心的单元获得较高的权重值,但靠近过滤边界的低密度单元和与中心相隔较远距离的单元会对中心单元的敏度产生较大的影响,拓扑优化的边界会出现过度磨平的优化结果,容易产生边界扩散的问题。

采用Sigmund敏度过滤方法的拓扑优化结果如图1所示。

图1 Sigmund敏度过滤方法拓扑优化结果

2.2 分区加权敏度过滤方法

为了解决Sigmund敏度过滤中边界扩散的问题,笔者提出了一种分区加权敏度过滤方法。其本质就是将原敏度过滤区域划分为两个子区域,在不同的子区域采用不同的加权因子,以确保中心单元所处的内部子区域提高单元敏度影响,同时降低远离中心单元的外部子区域的敏度影响权重,有效抑制出现大量灰度单元的边界扩散现象。

敏度过滤划分区域示意图如图2所示。

图2 敏度过滤划分示意图

区域Ⅰ和区域Ⅱ为内部和外部子区域,通过改变区域Ⅰ的过滤半径rn和最小过滤半径rmin的大小,可以改变边界A和边界B的相对位置,进而改变各区域所包含的单元数目。

为保证靠近中心单元的敏度值不被区域Ⅱ内的单元影响,从而导致权重降低,可直接将区域Ⅰ内单元加权因子赋值为1,以提高区域Ⅰ单元对目标函数敏度的影响权重。

区域Ⅱ的加权因子由设定的指数函数来确定,在区域Ⅱ内该指数函数能实现靠近边界A时的权重缓慢降低,而靠近边界B时权重显著降低的效果,进一步弱化区域Ⅱ内单元敏度对中心单元的影响。

在保证Sigmund敏度过滤的优化稳定性的前提下,可以更加有效地抑制边界扩散问题,获取清晰的拓扑优化结构边界。

其加权因子的数学表达式为:

(8)

其中:

(9)

笔者通过设定η值的大小来改变内部和外部两个子区域分别所包含的单元数。

由式(8)可知:区域Ⅱ内加权因子的最大值可以通过调节修正系数k来改变。为避免边界A处过渡单元敏度值出现过大的差值,从而产生边界扩散现象,应尽量保证区域Ⅱ的权重因子最大值与区域Ⅰ保持一致,故笔者将k值设定为0.67可保证函数的连续性。

由于修正系数β影响指数函数的曲率变化程度,笔者以最小过滤半径rmin取2.5,区域Ⅰ的过滤半径rn取1.0为例,分析β取不同值时指数函数图像的变化趋势。

不同β值的加权因子函数图像如图3所示。

图3 不同β值的加权因子函数图像虚线—Sigmund敏度过滤方法的加权原理图像;实线—不同β值指数加权的分区敏度过滤方法的原理图像

图3中,当Δd(e,i)<1时,指数加权函数取值明显皆大于线性加权函数取值,可有效提高区域Ⅰ内单元敏度的影响;当Δd(e,i)>1时,在靠近边界A处函数缓慢变化,靠近边界B处加权因子取值已经显著减小接近于0,靠近边界B的单元敏度影响权重得到了进一步降低;

β值越小,加权因子靠近边界A处函数变化越缓慢,权重不至于快速下降,从而导致权重过小,靠近边界B处权重皆处于较低值。经综合考虑,为保证函数变化趋势满足优化合理性,笔者将修正系数β设定为1;

而修正系数设定并非固定不变,也需根据具体算例略加调整。

该处笔者采用的分区加权敏度过滤方法易于编程实现,且抑制边界扩散等的效果明显。

综合考虑加权因子对优化结果的影响,在式(7)的基础上,敏度过滤方法为:

(10)

2.3 敏度过滤方法的衡量指标

为了衡量拓扑优化的结果,笔者引入了离散率、灰度率、柔度值和迭代次数四大衡量指标。其中,离散率反映灰度单元密度偏离0和1的幅度,是判断拓扑优化结果是否收敛到离散解的主要依据,当离散率为0时,结构完全离散。

离散率数学表达式为:

(11)

上式中,初始单元密度设为0.5,此时Dh为100得到最大值,表示所有单元都未发生离散。而当所有单元密度为0或1时,Dh为0得到最小值,表示所有单元都已发生离散。故在一般情况下,在拓扑优化的过程中,要尽可能使其离散率达到最小值,以得到更优的离散单元结构。

通过灰度率可以量化优化结果中存在灰度单元的多少,反映拓扑优化结果中灰度单元的占比程度。

灰度率数学表达式为:

(12)

式中:ρj—密度处于0～1之间的灰度单元;v0—预设的体积分数。

当灰度率Gh越大时,灰度单元占比越大,边界扩散现象越明显。所以在优化过程中,要尽可能使灰度率减小,以有效弱化边界扩散的问题。

柔度值反映结构在受力时的稳定性。柔度值越小,刚度值越大,结构受力时抵抗弹性变形的能力越强,变形越小,结构的稳定性越好。

而迭代次数则反映拓扑优化过程的优化时间和优化效率。迭代次数越少,优化时间越短,优化效率就越高。所以针对柔度值和迭代次数,在优化过程中也是尽可能使其达到较小值,优化结果会达到更优。

2.4 敏度过滤方法不同优化参数的设定

2.4.1 惩罚因子P

基于变密度法的分区加权敏度过滤方法的模型中,为了获得更加清晰的优化结果,确保优化结果的可制造性和合理性,笔者引入合适的惩罚因子P来对中间密度进行惩罚。经过惩罚的材料单元密度值将会快速趋近0或1,使模型能更好地趋近基于离散变量的拓扑结构模型。

SIMP密度惩罚函数图如图4所示。

图4 SIMP密度惩罚函数图

由图4可知:P越大,中间密度单元越易于趋近0,趋近速度也越快,所以理论上P越大越好。

为了研究P值变化对各衡量指标的影响,并确定惩罚因子P的值,笔者对经典二维应力结构进行拓扑优化。假定修正参数k=0.67、β=1.0,优化体积为0.3,可得到P不同取值时在柔度值、迭代次数、离散率及灰度率方面的分析结果。

不同P值的计算结果如表1所示。由表1和优化结果可知:当P≤2时,虽然迭代次数较少,但会产生大量多孔和细长杆的不可制造结构,不符合拓扑优化的优化要求;当2

表1 不同P值的计算结果

不同P值对衡量指标的影响折线图如图5所示。

图5 不同P值对衡量指标的影响折线图

由图5可看出:当2

当P>4时,拓扑优化结果与前面差异显著,出现较多的中间密度单元,迭代次数和柔度值也大幅增加。

综上所述,惩罚因子P设定为3.5比较合理,可以满足大部分算例的优化要求,特殊情况下也可适当调整P值。

2.4.2 最小过滤半径rmin

由图2的划分示意图可知:当区域Ⅰ的过滤半径rn确定后,最小过滤半径rmin的大小直接决定区域Ⅱ包含的单元个数,从而影响敏度过滤的优化结果;

rmin值过小时,区域Ⅰ相对于区域Ⅱ来说包含过多单元,无法有效弱化单元距离加权平均的效果;

rmin值过大时,区域Ⅰ包含相对较少的单元,迭代过程中极易出现数值不稳定问题。

因此,为研究rmin值变化对各衡量指标的影响,并确定rmin的值,笔者同样对算例1对应的二维应力结构进行拓扑优化。

假定惩罚因子P为3.5,rmin的取值对于各衡量指标结果的影响,即不同rmin的衡量指标折线图,如图6所示。

图6 不同rmin的衡量指标折线图

由图6可知:当rmin>4时,迭代次数虽变化并不明显,但离散率、柔度值和灰度率都呈增大趋势,会导致优化结果缺乏可制造性;而当rmin为3时,离散率、柔度值和灰度率都处于最小值,且迭代次数也处于较小值。

综上所述,最小过滤半径rmin设定为3比较合理,在保证优化效率的前提下,可以获得更加理想的拓扑优化结构。

3 数值算例分析与验证

在MATLAB-2019a软件环境下,笔者采用拓扑优化算例,对分区加权敏度过滤方法的可行性和有效性进行验证,针对算例模型均采用平面四节点四边形单元进行离散化,对模型的材料属性和结构尺寸进行无量纲化处理[16]。

假定实体材料的弹性模量为1,泊松比为0.3,惩罚因子P为3,最小过滤半径为3.5。

3.1 矩形板模型及优化结果分析

算例1的二维平面应力结构设计区域设为120 mm×40 mm,两侧的中间节点处采用固定约束,一集中载荷F=1 N作用于下边界中心位置。

优化区域被划分为120×40个矩形四节点单元,优化结构许用材料体积分数设为0.3。

矩形板结构设计区域如图7所示。

图7 矩形板结构设计区域

根据算例1为单工况、两个位移约束的情况,笔者通过将Sigmund原敏度过滤方法以及各文献中对敏度过滤改进后的方法,和分区加权敏度过滤方法作对比,验证该方法的可行性和有效性。

不同方法对矩形板结构的优化结果如图8所示。

图8 不同方法对矩形板结构的优化结果

由图8可以看出:4种方法的最终优化拓扑结构大致相似,但在边界清晰度上,除了Sigmund敏度过滤方法的边界有明显的边界扩散现象外,其他3种方法的边界扩散现象均呈现大幅度的减弱;

方法一的优化结果会出现细长杆结构和孔洞结构,缺乏可制造性,方法一的优化结果和该方法的优化结构相近。

不同方法对矩形板结构的优化数据对比如表2所示。

表2 不同方法对矩形板结构的优化数据对比

由表2可知:相对于Sigmund敏度过滤方法,其他3种方法在离散率、灰度率和柔度值方面都大幅降低;在迭代次数、离散率和灰度值方面,相对于方法一和方法二,分区敏度过滤方法都更低一些,尤其是迭代次数和灰度率方面尤为显著,在提高优化效率和抑制灰度单元方面有一定的优势。

通过对比可知,该方法在有效获得清晰拓扑结构的同时,具有较快的收敛速度和较低的离散率,对边界扩散现象的控制较好,能获得更佳的优化结构。

3.2 MBB梁模型及优化结果分析

算例2的几何尺寸为240 mm×40 mm,顶端中部承受1 N的外载荷,目的是测试该方法在不同网格密度情况下是否会出现棋盘格、多孔材料和网格依赖性等数值不稳定现象。

为了简化优化的过程,笔者对1/2MBB梁结构进行网格划分,得到了单元总数量为120×40的有限元模型优化结构,许用材料体积分数设为0.3。

MBB梁设计区域如图9所示。

图9 MBB梁设计区域

笔者通过选用不同参数,分析及验证该方法在不同网格密度时是否具有可行性。

不同网格划分对MBB梁的优化结果如图10所示。

图10 不同网格划分对MBB梁的优化结果

由图10可以看出:在对MBB梁的优化结果分区方面,敏度过滤方法同样是优于Sigmund敏度过滤方法的;在不同网格划分的情况下,对MBB梁进行优化,均可得到理想的拓扑优化结构,拥有较为清晰的结构边界,且均没有求解不稳定现象的出现。

不同网格划分对MBB梁的优化数据对比,如表3所示。

表3 不同网格划分对MBB梁的优化数据对比

由表3可以看出,采用不同网格划分进行优化时,虽然在数值上各大衡量指标有所波动,但相对于Sigmund敏度过滤方法都有所减小,均体现出对灰度单元的抑制作用。

3.3 Michelle结构模型及优化结果分析

算例3的结构为两端固支的传统Michelle型结构,设计区域为120 mm×40 mm,左下角和右下角节点处采用固定约束,在图示结构上部1/4和3/4节点处皆受到F1=2 N的竖直向下载荷的作用;同时,结构下端中间节点处受到F2=1 N的竖直向下的载荷作用。

网格划分成120×40个四节点等尺寸的平面应力单元,优化结构许用材料体积分数设为0.3。

Michelle结构设计区域如图11所示。

图11 Michelle结构设计区域

为了进一步分析及验证该方法在多载荷情况下的可行性,笔者将采用分区敏度过滤方法得到的拓扑优化结果,与采用Sigmund敏度过滤方法得到的结果进行对比。对比采用体积比为0.3的约束条件,过滤半径取rmin=3.5。

不同方法对Michelle结构的优化结果如图12所示。

图12 不同方法对Michelle结构的优化结果

由图12可知:在多载荷情况下,该方法同样具有良好的可行性,可以抑制边界扩散的现象,避免产生细长杆等不可制造结构。

采用不同方法对Michelle结构的优化数据对比,如表4所示。

表4 不同方法对Michelle结构的优化数据对比

由表4可知:与Sigmund敏度过滤方法相比,虽然分区敏度过滤方法的迭代次数增加了,但其拓扑优化离散率、灰度率和柔度值都大幅降低,优化结构的柔度收敛值较小,能有效控制其离散率及灰度率指标。

4 结束语

在基于变密度法的结构拓扑优化中,优化结果常存在边界扩散和细长杆结构等数值不稳定的问题,为此,笔者提出了一种基于变密度法的分区加权敏度过滤方法,即通过引入加权因子对划分的内、外两个敏度过滤区域分别进行了加权处理,解决了拓扑优化中存在的灰度单元和细长杆结构等数值不稳定问题,使拓扑优化结构具备较好的可制造性。

研究结果如下:

(1)相对于Sigmund敏度过滤方法,该方法在离散率和灰度率方面提升显著,减小率达到80%以上,可以有效地解决边界扩散和细长杆结构的问题;

(2)该方法在迭代次数和柔度值方面有一定的提升,而且相对于其他敏度过滤方法,解决了它们只能改善边界扩散,但优化效率低的问题,不仅提高了拓扑优化的迭代效率,还提升了其结构刚度,有利于拓扑优化的后处理操作;

(3)各个典型算例分别从多个角度(单一载荷、多载荷、不同网格划分的约束和受载条件)进行了结构拓扑优化,验证了该方法具有一定的可行性。

但是,该方法目前只适用于二维平面结构的拓扑优化。因此,在后期的研究中,笔者将进一步探究该方法在三维结构拓扑优化中的可行性和有效性。