唐腾飞, 高隆隆, 廖利辉, 习毅, 李宝仁
(华中科技大学 机械科学与工程学院 FESTO气动中心,湖北,武汉 430074)
结构拓扑优化是一种针对给定结构设计域的初始优化设计. 在给定的约束条件下,它通过在设计域内寻找合适的材料分布或者结构形式的分布,使得目标函数达到最优. 20世纪80年代前的结构优化研究主要集中在设计部件的尺寸、形状优化,自从1988年Bends∅e针对拓扑优化创造性的研究[1],将微结构复合材料引入结构优化,解决了各种优化方法用于拓扑结构优化的难题,随后的几十年内,拓扑优化方法得到了迅猛发展. Xie等[2]提出了渐进结构优化方法ESO,属于结构拓扑优化方法的一种. ESO方法通过逐步删除设计域结构中具有低效单元灵敏度的材料,最终得到优化结构[3]. 因其方法原理简单且易于理解,并在工程实践中具有良好的应用性,得到了科研人员与工程师的广泛关注. 目前,ESO方法已经应用于车辆结构设计、蝶阀结构设计等方面并获得很好的效果[4-6].
在ESO方法中,删除率rRR和进化率rER控制着优化结果的优化趋势,属于关键的影响参数. 删除率rRR和进化率rER两参数设置过大,会导致构不稳定甚至优化失败;参数设置过小,会使得优化过程中删除单元量增加,从而引起可能的结优化过程十分漫长,效率低且计算量增加.
本文针对传统渐进结构优化法中存在的上述问题,提出了优化过程中参数自适应的方法,并引入约束罚函数,达到优化过程快与优化结果稳定目标.
为了不失普遍性,以体积作为约束条件,本文对基于刚度的拓扑优化问题进行研究. 问题可以表示为
(1)
(xi=0或1)
式中:f为外载荷向量;u为位移向量;C为平均柔度值,其值越小,反映出刚度越大;V为设计域的体积,Vi为单元i的体积,i为单元序号,N为单元总数;设计变量xi表示单元的存在于删除状态.
基于均匀应变能设计准则的渐进结构优化方法,其单元应力能作为敏度值,为了有效删除具有低效敏度值的单元,选取余弦函数作为基础的递减函数[7],进化率rER自适应的单元删除方法如下表示:
(2)
式中:A为幅值参数;αmin为最小单元应力能值;αmax为最大单元应力能值;rER为进化率;m与p值为罚值. 通常来说,A依据经验选择为0.01,本文选择0.015作为A值,使得进化率rER取值范围为[0,0.015],为了表述方便,将0.01作为一个标准值,因此进化率rER取值范围为[0,1.5]的标准值中变化. 罚值m与p的适当选取,有助于优化过程快速收敛且减小收敛过程中的震荡问题.
为了避免不合理结构引起的个别单元的单元应力能极端变化,从而无法有效反映出结构体的真实刚度水平,本文选取了如下策略:当前设计域中,按单元应力能值的大小排序,分成相等的两部分,单元应力能值较大的集合与单元应力能值较小的集合. 最小应力能αmin是单元应力能值较小的集合的平均值,最大应力能αmax是单元应力能值较大的集合的平均值. 因其定义,αmin/αmax的取值范围为(0,1). 本文选取了4种结果进行分析,分别是: ① 传统固定定进化率rER,rER=0.01; ② 取p=0,m=1的自适应进化率rER1; ③ 取p=2,m=1的自适应进化率rER2; ④ 取p=3,m=1.5的自适应进化率rER3. 如图1所示,随着单元应力能比值的增大,采用自适应进化率rER1,rER2,rER3值随之减小,各进化率的变化形式不一致,rER1单调减小,rER2,rER3属于‘z’型变化趋势.
在拓扑优化的求解过程中,因为结构设计域的有限元单元采用低阶单元,经常会引起棋盘格的问题,这种问题产生的主要原因是对单元之间采用对角连接的结构的刚度的扩大导致. 这种问题导致优化后的结构无法很好地应用于工程实际中,因此需要抑制棋盘格结构的产生,采用灵敏度平均过滤算法属于其中一种有效的方法[8]. 在不考虑计算资源的情况下,拓扑优化方法形成的优化结果的局部可以形成非常细致的微小结构[9],实际中则需平衡计算资源与优化结果的工程可行性. 为减少计算所消耗的资源,本文采用了一种棋盘格过滤算法[10],这种算法可以表示为如下形式.
首先,找到单个单元节点所连接的所有单元,通过所计算的单元敏度值,将节点所连接的所有单元敏度值进行平均,这里,每个单元的体积被认为是相同. 以上过程采用公式表示为
(3)
式中:αnode为节点的敏度值,通过周边相连接的单元的敏度平均值获得;αj为每个单元的单元应力能;M为节点连接单元数量.
其次,获得每个节点的敏度值后,找到每个单元所对应的节点集合,计算节点敏度值的平均值,作为单元的应力能值,也就是单元的敏度值. 表示为
(4)
式中:αele为单元的应力能值;αnode为上一步获得的节点敏度值;K为单元连接的节点数.
参数自适应渐进结构优化方法流程的约束是体积约束,其流程为:
① 离散化设计域结构建立有限元结构网格;
② 有限元软件获得单元应力能,通过式(4)计算出过滤后的单元应力能;
③ 通过删除率dRR移除具有最低敏度值的单元集合;
④ 如果达到稳定状态,依据式(2)改变进化率rER值;
⑤ 计算删除率dRR的值;
⑥ 重复②~⑤直到达到约束要求.
短悬臂梁结构的设计域及载荷边界约束的模型如图2所示. 短悬臂梁的厚度为0.1 mm. 网格尺寸选择为5 mm,因此设计域的有限元单元有32×20个. 设计域左侧端面被固定,在右侧中部设有垂直向下的载荷,载荷P=1 000 N被用于优化过程中. 为了方便,材料属性在不同仿真例子中均相同,材料的杨氏模量E=100 GPa,泊松比为ν=0.3,初始删除率rRR=0.01. 设计域的体积约束设置V=0.4,意味着40%的设计域最终被保留. 本文的ESO优化方法采用的是“硬删除”方法,所有被移除的单元被直接移除.
采用如图1所示的4种方法得到的结果如图3所示,4种不同的优化结果具有相似的优化结构,但又各具不同之处,主要在于主体结构中的“X”型支持处于不同位置.
优化过程中,采用自适应进化率rER1并在最后仿真时出现结构过度删除,属于不合理现象,取值的体积为42%.
采用固定进化率与不同自适应进化率的4种方法所获得的优化结果的最大应力值与柔度值如表1所示. 需要注意的是,因为优化目标针对的是结构刚度值,具有越小柔度值的结果越好;最大应力值只是作为仿真结果仅作为参考. 从4种结果的柔度值来看,具有参数自适应进化率γER1和γER3两种结果优于固定进化率rER的结果,体现出参数自适应可以获得更优方案的潜力.
表1 固定与自适应进化率的最大应力值与平均柔度值
由图4反映出的设计域体积变化规律可知,采用固定进化率rER的方法经过了41步得到最后结果,可以明显发现,在35步迭代后,采用固定进化率每次删除的单元数都较多,导致在最后迭代时的体积迅速下降,从优化的角度来说容易产生不稳定现象;采用参数自适应进化率rER3的方法经过35步有限元仿真即得到最后结果,更快地完成了优化进程,而且在最后优化迭代时,体积降低间隔小,产生的结构更具有稳定性. 采用参数自适应进化率rER1和rER2的方法,所用的迭代步数分别为43次和58次,说明随着优化进程的进行,进化率降低至0.01以下,导致了更多的迭代步数.
综合以上分析结果,可以发现,具有罚值p=3,m=1.5的自适应进化率rER3的方法,从目标的刚度值与优化迭代次数,明显优于传统的固定进化率rER方法;从结构应力分布和以及可参考的最大应力来看,也具有更好的结构形式与效果.
图5为不同方法在优化过程中的进化率ER值的变化规律,自适应进化率在前期大于固定的进化率0.01,可以移除更多的单元数量,减少仿真迭代步数,随着仿真的进行,自适应进化率逐步减小;当优化进入结束过程时,自适应进化率小于固定的进化率0.01,从而减少了优化结束阶段单元的删除量,保证了优化结构的稳定性,因此,采用自适应进化率实际上对设计流程中的单元删除数量有较好的调节作用.
刚度是指材料或结构在受力时抵抗弹性变形的能力,调节阀刚度越高,承受压力的变形更高,同时会改变阀体结构固有频率等参数,具有更好的安全性. 原始调节阀的模型如图6所示,调节阀口径为DN100. 通过流道提取得到如图7所示的初始设计域. 设计域的边界条件如图所示,为两端固定,流道区域施加3 MPa压力.
设计采用Abaqus有限元软件,单元类型选择为C3D10,共18 346个单元. 材料属性为杨氏模量E=207 GPa,泊松比为ν=0.3. 优化的约束为体积约束,约束体积为0.3的刚度设计,选取罚值p=3,m=1.5的参数自适应渐进结构优化方法.
图8为优化结果,与现有结构相比,目前的结构已经近似最优解外,有5处存在具有特点之处:入口端A处壁厚较厚,出口B处壁厚较薄,这与传统的等壁厚设计不一致;D处拐角弧度较大,C处弧度较小,E处圆柱体厚形成梯形设计. 限于制造工况,最终得如图9所示的优化模型.
优化后调节阀模型有限元仿真结果表明调节阀的柔度从1.305×1010N·m减小1.192×1010N·m,相当于刚度增加10%,在相同压力下形变量更小. 结构应力采用Von Mises应力表示,从63 MPa略微减小至58 MPa,主要原因是最大应力值所处的拐角处为非设计区域. 结果表明,自适应的渐进结构优化方法可以很好地运用于三维的结构模型优化,对阀类设计提供参考.
针对传统的渐进结构优化方法过程中参数值固定引起优化收敛速度缓慢,本文提出了一种基于动态进化率的参数自适应渐进结构优化方法,引入了自适应进化率函数以及罚值,设计了新的优化流程. 通过典型算例,得到罚值较优的自适应参数进结构优化方法. 通过应用本文方法于调节阀设计,获得了提升10%刚度的调节阀结构,为阀类结构设计提供参考.