基于单元过滤的自支撑结构拓扑优化方法

王雷,闫素娜,赵强强,洪军

(1.西安交通大学机械工程学院,710049,西安;2.西安交通大学现代设计及转子轴承系统教育部 重点实验室,710049,西安;3.中国电子科技集团公司第二十七研究所,450047,郑州)

拓扑优化受到传统制造技术的限制,实际工程应用的拓扑优化结构仍然有一定的局限性,没有充分发挥拓扑优化的优势。近年来,增材制造技术得到了快速发展,与传统制造技术相比,它具有能够实现高效率、低成本、制造复杂结构的突出优势[1],为拓扑优化结构的制备提供了强大的工具,极大地拓宽了结构设计的空间。因此,将拓扑优化与增材制造结合,研究面向增材制造的创新设计方法具有广阔前景[2]。然而,增材制造技术存在一些工艺约束,制约着优化结构的制备,例如悬垂特征,打印时需要在其下方添加支撑防止坍塌,打印完成后再通过人工去除,不仅浪费打印材料和打印时间,而且去除支撑会对零件性能和表面质量造成损坏[3]。

为了减少使用支撑,一些学者从工艺参数方向入手,利用优化算法寻找能够使支撑体积最小的最佳成型方向。其中,Das等提出了一种利用多目标优化方法,在满足零件GD&T标准的同时,寻找能够使支撑体积最小的成型方向[4];Morgan等提出了一种利用优化算法寻找使支撑体积最小的成型方向的方法[5];Strano等提出了一种两步优化设计结构的方法,第一步寻找结构最优打印方向,第二步基于此方向优化设计多孔支撑结构,使支撑体积最小[6]。上述通过优化算法寻找最佳成型方向的方法只能使支撑体积达到最小,不能完全避免使用支撑。

另一些学者从结构设计方向入手,研究面向增材制造的自支撑结构设计方法。其中,Laery等提出了一种基于熔融沉积成型(FDM)工艺的自支撑条件,通过修改结构构型,使结构打印时无需支撑[7];Li等提出了一种面向激光选区熔化(SLM)工艺的零件轻量化/无支撑优化设计方法,第一步获得拓扑优化结构,第二步使用自支撑结构元素修改结构构型[8];桂馨提出了一种基于固体各向同性材料惩罚模型(SIMP)法和能量均匀化法的自支撑结构多尺度拓扑优化,在宏观结构优化阶段,利用SIMP法获取自由分布的宏观材料空间,从而划分微结构类型区域,在微结构设计阶段,进行自支撑微结构设计,最终获得自支撑结构[9]。上述通过修改结构构型来获得自支撑结构的方法会造成结构体积的增加,而且添加的材料会影响结构性能。

Garaigordobil等提出了一种使用轮廓检测算法提取结构边界,通过显示不等式约束条件控制边界法向量的倾斜角,从而获得自支撑结构的方法[10];Qian提出了一种基于密度梯度信息确定结构边界,通过显示约束条件控制边界倾斜角的方法[11];Zhang等提出了一种利用局部单元密度确定边界法向量的方法,通过控制法向量倾斜角来控制边界倾斜角和悬吊特征尺寸[12]。上述通过边界法向量控制边界倾斜角的方法需要通过密度梯度信息提取结构边界法向量信息,增加了优化时间。Zhao等提出了一种通过卷积算子判断单元的支撑性,通过显示约束条件控制非自支撑单元数量,最终获得自支撑结构的方法[13],此方法需要在传统拓扑优化的基础上确定成型方向,再进行自支撑结构优化设计;Kuo等提出了一种使用连续逻辑聚合函数评估单元支撑状态的方法,以显示约束形式加入到结构拓扑优化方法中,优化设计自支撑结构[14];邹君等提出了一种使用显示约束函数表征结构自支撑约束违反程度,并控制此函数值的方法优化设计自支撑结构[15]。上述通过显示约束条件控制非自支撑单元的数量的方法会增加优化问题求解难度。

Gaynor等提出了一种通过映射单元设计变量的方法实现结构最小特征尺寸和最大悬垂角度的控制[16],但其扇形支撑区域包含太多单元,计算复杂。Langelaar提出了一种过滤方法,模拟增材制造成型过程逐层过滤单元进行自支撑结构的优化设计,并将结构优化、支撑结构布局和成型方向3个目标结合,寻求最优结构设计和成型方向[17-18];Pellens等提出了4种不同形式的过滤器,以不同顺序组合进行过滤,将特征尺寸和悬垂角度约束合并到结构拓扑优化问题中进行结构优化设计[19];Wang等提出了一种基于分层过滤策略的SIMP改进方法和一种基于悬垂灵敏度分析的双向渐进结构拓扑优化(BESO)改进方法以减少优化设计中悬垂结构的产生[20];杜义贤等通过建立每个单元的可打印密度与相应支撑区域内单元的可打印密度之间的关系来过滤掉不符合约束条件的结构,优化设计自支撑结构[21];Zou等提出了一种基于max/min函数识别非自支撑单元,通过不等式约束控制非自支撑单元数量,得到自支撑结构的方法[22]。上述过滤方法中采用近似max/min函数通过判断上下相邻单元的密度关系,确定单元的自支撑性,此判断条件不符合实际逐层堆积成型原理,会出现自下向上密度逐渐增加的现象。

Luo等提出了一种自支撑封闭空腔结构拓扑优化方法,其采用非线性虚拟温度法(N-VTM)识别封闭空腔,通过多重过滤识别悬垂面,基于对数函数的约束控制封闭空腔的最小悬垂角[23];Guo等提出了一种基于移动可变形组件的自支撑结构设计方法,通过控制组件的倾斜角度保证组件自支撑,最终获得自支撑结构[24];Emiel等将设计域分层,以设计域内每个点材料到达时间为判断依据判断支撑性,优化设计自支撑结构[25];王亚光等提出了一种基于参数化水平集法的悬垂约束处理方法[26];韩永生等提出了一种基于BESO法的无需支撑(自支撑)拓扑优化方法[27]。

本文模拟增材制造逐层堆积材料的成型过程,提出一种基于Heaviside函数的自支撑单元过滤法则,并基于此构建面向增材制造的自支撑结构拓扑优化方法,旨在解决传统拓扑优化结构在增材制造过程中需要使用支撑的问题。最终,通过数值案例验证本文所提方法的有效性,结果表明本文所提方法能够在4个方向上对二维结构进行优化设计,实现打印过程自支撑,并且在保证结构刚度的前提下,节省打印材料和打印时间。

图1 增材制造悬垂特征模型示意图Fig.1 Schematic diagram of overhang features in additive manufacturing

1 自支撑单元过滤法则

1.1 增材制造悬垂特征

图1给出了增材制造悬垂特征模型示意图,悬垂特征模型主要包括倾斜角α小于自支撑临界角的悬垂面、悬吊面和悬吊线。

为了避免使用支撑,对于悬吊面和悬吊线特征,设计阶段可以在其下方添加特征连接至打印平台;对于悬垂面特征,Thomas研究了SLM工艺自支撑临界角为45°[28],则可以在设计阶段使倾斜角α大于45°,从而无需添加支撑。

1.2 自支撑单元过滤法则

图2给出了悬垂面特征模型及其有限元模型,倾斜角为45°,采用长宽比为1∶1的四节点矩形单元离散此模型,构建如图2b所示悬垂面特征有限元模型。

(a)悬垂面特征模型 (b)悬垂面特征有限元模型图2 悬垂面特征模型及其有限元模型 Fig.2 Overhang model and its finite element model in additive manufacturing

取单个单元作为研究对象,图3给出了自支撑单元过滤法则示意图,成型方向为竖直向上,目标单元xij与其支撑单元集3个单元形成最小倾斜角为45°的倾斜面,其中i和j分别表示单元水平和垂直方向的位置。

图3 自支撑单元过滤法则示意图Fig.3 Schematic diagram of filtering rule of self-supporting element

目标单元是否自支撑与支撑单元集中3个单元密度有关,若支撑单元集中至少有一个单元密度大于0,则目标单元为自支撑单元。基于上述条件可得自支撑单元判断法则为:若支撑单元集中单元平均密度大于0,则目标单元为自支撑单元,反之目标单元为非自支撑单元。对于设计域全局单元,将逐层判断是否自支撑。

1.3 自支撑过滤特殊单元

设计域中存在一些特殊单元在优化算法中要特殊处理,包括边界单元和底层单元,自支撑过滤特殊单元如图4所示。

(a)边界单元 (b)底层单元图4 自支撑过滤特殊单元Fig.4 Special element of self-supporting filtering

边界单元是指处于设计域两侧边界的单元,如图4a所示,它的特殊性在于支撑单元集中只包含2个单元。底层单元是指处于设计域底部的单元,如图4b所示,与成型方向有关,它的特殊性在于直接与打印平台连接,无论如何都是自支撑单元。

2 自支撑结构拓扑优化模型

2.1 自支撑单元过滤

采用自支撑单元过滤的方法判断单元是否自支撑,保留自支撑单元,删除非自支撑单元,自支撑单元过滤表达式如下

(1)

(2)

其中β为近似参数,T为阈值,Me为单元e支撑单元集中的单元平均自支撑密度,表达式如下

(3)

图5 Heaviside阈值函数曲线Fig.5 Curve of Heaviside threshold function

2.2 单元密度过滤

棋盘格式解是指结构中出现实体单元与空白单元交替的现象,形成一种棋盘状的构型,棋盘格式解如图6所示,此结构在实际工程中无法应用。

图6 棋盘格式解Fig.6 Checkerboard solution

为了消除棋盘格式解现象,本文采用单元密度过滤法,其核心思想是利用过滤后的单元密度替代单元实际密度[29],单元物理密度表达式如下

(4)

Ne={i|‖Xi-Xe‖≤rmin}

(5)

其中Xi为单元i的坐标,Xe为单元e的坐标,rmin为过滤半径。

Hei表达式如下

Hei=max(0,rmin-Δ(e,i))

(6)

式中:Δ(e,i)为单元i和e中心距。密度过滤区域如图7所示,是以目标单元e为圆心、过滤半径rmin为半径形成的圆形区域。

图7 密度过滤区域Fig.7 Density filtering area

2.3 自支撑结构拓扑优化模型

基于SIMP变密度法,以最小化柔度(最大化刚度)为目标函数,以结构体积分数为约束条件,构建自支撑结构拓扑优化模型,数学模型如下

(7)

式中:c(x)为结构柔度;v为约束条件;x为单元设计变量;V(x)为结构体积;V0为设计域总体积;f为结构体积分数。

3 基于MMA的拓扑优化模型求解

本文所提的自支撑结构拓扑优化模型采用MMA法求解[30],更新单元设计变量x。

3.1 有限元分析

采用长宽比为1∶1的四节点矩形单元离散设计域,初始单元刚度矩阵为k0,基于SIMP模型计算单元e的单元刚度矩阵ke,表达式如下

(8)

SIMP模型单元密度与弹性模量之间的关系可表示为

(9)

惩罚因子P是为了惩罚灰色单元密度,使其趋向0或1[31]。图8给出了SIMP模型指数函数曲线,不同曲线代表不同P值,P值越大函数值越趋近于0或1。

图8 SIMP模型指数函数曲线Fig.8 Exponential function curve of SIMP model

通过组装单元刚度矩阵求得整体刚度矩阵K,建立刚度方程求解整体位移矩阵U,刚度方程表达式如下

KU=F

(10)

式中:K为整体刚度矩阵;U为整体位移矩阵;F为已知载荷矩阵。

最终,通过U和K计算目标函数,结构柔度表达式如下

(11)

式中:n为离散化设计域单元的数量;ue为单元位移矩阵。

3.2 灵敏度分析

MMA求解需要分析目标函数和约束条件对设计变量的灵敏度,灵敏度信息决定设计变量的更新。其中,目标函数对设计变量的灵敏度表达式如下

(12)

式中:c为目标函数。

约束条件对设计变量的灵敏度表达式如下

(13)

式中:xf为单元自支撑密度。

3.3 自支撑结构拓扑优化流程

本文首先采用自支撑单元过滤,从更新的单元中过滤得到自支撑单元,保证设计域内只留下自支撑单元,再采用密度过滤消除棋盘格式现象,并且可以通过密度过滤避免结构中出现细小特征。然而,对于相反的过滤顺序,通过试验可知同样可以得到最优的自支撑结构,但是优化结构中会出现细小特征,优化结构中的细小特征如图9所示。

图9 优化结构中的细小特征Fig.9 Small features in optimal structures

图10给出了自支撑结构拓扑优化流程图,具体步骤如下:

(1)前处理,定义结构几何模型,添加边界条件,施加载荷,设置优化参数;

(2)采用四节点矩形单元离散设计域;

(3)初始化单元设计变量x;

(4)自支撑单元过滤,单元密度由单元设计变量x转换到单元自支撑密度xf;

(5)密度过滤,单元密度由单元自支撑密度xf转换到单元物理密度xp;

(6)有限元分析,基于单元物理密度xp和初始单元刚度矩阵k0计算整体刚度矩阵K,求解整体位移矩阵U,计算目标函数和约束条件;

(7)灵敏度分析;

(8)采用MMA法求解优化模型,更新单元设计变量x;

(9)根据收敛条件判断结果,如果不收敛则从步骤4继续循环,直至收敛;如果收敛则输出优化结果,并结束计算。

图10 自支撑结构拓扑优化流程图 Fig.10 Flow chart of self-supporting structures topology optimization

4 自支撑结构拓扑优化

应用本文所提的自支撑结构拓扑优化方法解决MBB梁拓扑优化问题。图11给出了MBB梁对称设计域,长宽比为3∶1,设计域左侧添加滚动支撑,右下角添加滚动支撑,左上角施加沿-Y方向的单位集中载荷F。将设计域离散为180×60个单元,并初始化单元设计变量x为0.5。

图11 MBB梁对称设计域Fig.11 Symmetric design domain of MBB beam

4.1 优化参数确定

在本文所提出的自支撑结构拓扑优化方法中,惩罚因子P决定着对灰色单元密度的惩罚程度,P值越大,单元密度越趋近0或1,但是会导致函数非线性程度增加;阈值T决定单元自支撑条件,T值越大,要求支撑集单元密度越大,单元自支撑条件越严苛。因此为了确定一组合适的P值和T值,本文以Y方向的MBB梁优化设计为例,选择不同的参数P和T进行优化结果对比,最终根据优化结构的迭代步数和性能确定优化参数。

不同参数P和T优化结果对比如表1所示,可知惩罚因子P和阈值T对自支撑结构优化有重要的影响。T对结构优化的影响从T为0.1～0.4可知,迭代步数增大,结构柔度增大,特别在T为0.4时出现优化失败的结果,所以T越大,会导致结构性能下降,迭代收敛难度增加,甚至出现优化失败;P对结构优化的影响从P为2～5可知,结构中灰色单元明显减少,结构边界逐渐清晰,但是迭代步数增大,所以P越大,结果中灰色单元越少,但是会导致迭代收敛难度增加。

综合分析,本文旨在保证结构中不存在大量灰色单元的前提下,使迭代步数减小,结构性能增加,得到最优的自支撑结构。从表1可知,在参数P为3、T为0.1的组合下,得到的结构柔度最小,迭代步数最小,并且没有大量灰色单元存在。所以本文选择P为3、T为0.1作为算例的优化参数组合,最终MBB梁拓扑优化参数取值如表2所示。

本数值案例由Matlab R2020b(64位)、Windows 10、Intel(R)Core(TM)i5-8400cpu@2.8GHz CPU、8GB RAM实现。

4.2 MBB梁拓扑优化设计

分别应用传统拓扑优化方法和本文所提的面向增材制造的自支撑结构拓扑优化方法优化设计MBB梁。图12给出了传统拓扑优化MBB梁,迭代步数为565,柔度为194.8。

利用Autodesk Inventor Professional 2020提取优化结果数据,建立MBB梁三维模型,厚度设置为5 mm。利用Autodesk公司开发的3D打印辅助分析软件Autodesk Netfabb Ultimate 2020沿Y、-Y、X、-X4个方向进行成型仿真分析,成型仿真分析参数如表3所示。

表1 不同参数P和T优化结果对比

表2 MBB梁拓扑优化参数

图12 传统拓扑优化MBB梁Fig.12 MBB beams with traditional topology optimization

表3 成型仿真分析参数

图13给出了传统拓扑优化MBB梁分析结果,蓝色部分为支撑。表4给出了传统拓扑优化MBB梁成型数据统计。

表4 传统拓扑优化MBB梁成型数据统计

(a)Y方向分析结果

(b)-Y方向分析结果

(c)X方向分析结果 (d)-X方向分析结果图13 传统拓扑优化MBB梁分析结果Fig.13 Analysis results of MBB beam with traditional topology optimization

图14～17给出了Y、-Y、X、-X4个方向优化自支撑MBB梁及分析结果。图14a给出了最终的优化结果,迭代步数为616,柔度为233.6;图14b～17b给出了分析结果,无需支撑结构。

(a)Y方向优化结果

(b)Y方向分析结果图14 Y方向优化自支撑MBB梁及分析结果 Fig.14 Optimized self-supporting MBB beam in Y direction and analysis results

图15a给出了最终的优化结果,此结果出现振荡。

(a)-Y方向优化结果

(b)-Y方向分析结果图15 -Y方向优化自支撑MBB梁及分析结果 Fig.15 Optimized self-supporting MBB beam in -Y direction and analysis results

图16a给出了最终的优化结果,迭代步数为1 297,柔度为194.4。

(a)X方向优化结果 (b)X方向分析结果图16 X方向优化自支撑MBB梁及分析结果 Fig.16 Optimized self-supporting MBB beam in X direction and analysis results

图17a给出了最终的优化结果,迭代步数为1 487,柔度为194.5。表5给出了自支撑MBB梁成型数据统计。

(a)-X方向优化结果 (b)-X方向分析结果图17 -X方向优化自支撑MBB梁及分析结果 Fig.17 Optimized self-supporting MBB beam in -X direction and analysis results

表5 自支撑MBB梁成型数据统计

4.3 对比分析

Langelaar模拟增材制造逐层堆积的成型过程,采用过滤方法优化设计自支撑结构的方法[18]受到了广泛的认可,将本文所提方法与Langelaar方法和传统拓扑优化方法进行对比分析,总结本文所提方法的优点和不足。

Langelaar方法在文章中从北、东、南、西4个方向优化设计同样的MBB梁,将设计域离散为180×60个单元。表6汇总了应用3种方法的优化结果,cMBB为MBB梁柔度,cref为传统MBB梁柔度。对比可得基于本文所提方法得到的优化结构:

(1)在结构柔度方面,与传统拓扑优化方法相比,Y方向增大了19.9%,X方向减小了0.2%,-X方向减小了0.2%,与Langelaar优化方法相比,Y方向增大了13.1%,X方向减小了0.2%,-X方向减小了1.1%;

(2)在所需材料方面,与传统拓扑优化方法相比,Y方向体积减小了20.6%,-Y方向体积减小了18.0%,X方向体积减小了18.1%,-X方向体积减小了20.0%;

(3)在打印时间方面,与传统拓扑优化方法相比,Y方向减小了16.6%,-Y方向减小了14.5%,X方向减小了10.5%,-X方向减小了11.5%;

(4)在迭代步数方面,与传统拓扑优化方法相比,Y方向增大了9%,X方向增大了129.6%,-X方向增大了163.2%。

通过对比分析,可知本文所提方法能够沿4个方向对MBB梁进行优化设计,其中在X和-X两个方向上的结构柔度均优于传统拓扑优化方法和Langelaar提出的拓扑优化方法,虽然在某些方向上优化效率有所损失,但是能够在保证结构刚度的前提下,无需使用支撑,减少打印材料和打印时间的浪费。

本文针对45°临界角,基于SIMP变密度法,构建了面向增材制造的平面自支撑结构拓扑优化方法。基于本文所提方法未来可向:①采用可变长宽比四节点矩形单元构建临界角可变的自支撑单元过滤方法,研究适应于临界角可变的自支撑结构拓扑优化方法;②采用三维自支撑单元过滤方法判断单元自支撑性,此时目标单元为1个,支撑单元为5个,将本文所提方法拓展至三维自支撑结构优化设计;③采用基于Heaviside函数的自支撑单元过滤方法,基于变密度法多材料插值模型或非均质材料插值模型,构建面向增材制造的多材料或非均质材料自支撑结构拓扑优化方法。

表6 MBB梁优化结果对比

5 结 论

本文模拟增材制造逐层堆积材料的成型过程,提出了一种基于Heaviside函数的自支撑单元过滤法则。基于SIMP变密度法构建了面向增材制造的自支撑结构拓扑优化方法,并将本文所提方法应用到MBB梁拓扑优化中。通过数值案例可得出如下结论:

(1)本文所提出的一种基于Heaviside函数的自支撑过滤法则能够有效判断单元自支撑性,并逐层保留自支撑单元,删除非自支撑单元;

(2)本文所提的面向增材制造的自支撑结构拓扑优化方法能够有效解决传统拓扑优化结构在增材制造过程中使用支撑的问题,并且在不损失结构刚度的前提下,节省打印材料和打印时间;

(3)本文所提的面向增材制造的自支撑结构拓扑优化方法能够在4个方向上对平面结构进行优化设计,并且可得成型方向对结构构型、刚度和优化效率有很大的影响,可依据结构使用要求选择最优设计方案。