基于变密度法的结构强度拓扑优化策略

丁 卯，耿 达，周明东, 来新民

(上海交通大学 机械与动力工程学院， 上海 200240)

复杂工况下的结构设计必须考虑结构在外力作用下应力的分布情况，以避免局部应力过高导致结构失效等问题[1].现代工业的不断发展对具有高比刚度及高比强度结构的需求愈发迫切.结构拓扑优化方法基于最优化理论，结合计算机辅助工程(CAE)技术，自动地设计符合性能要求的工程结构，近年来已作为新兴的结构设计方法在航空航天及汽车等领域的高比刚度、高比强度结构设计中得到了广泛的研究与应用[2-3].目前结构比刚度问题的拓扑优化已得到充分发展，而考虑结构强度的拓扑优化方法仍存在许多问题[4-6].

基于变密度法的结构拓扑优化方法优化过程中结构拓扑变化灵敏，被广泛应用于实际工程结构的优化设计.但传统基于变密度法的结构强度拓扑优化方法受优化问题非线性影响，结构近似最大单元应力难以准确计算，应力变化趋势难以精准评估[7-8],且优化过程对优化参数的选择较敏感，不易收敛.同时最终优化结果往往存在较多无实际物理意义的灰度单元，后处理前后结构强度性能偏差较大，常需反复试错才能得到符合强度设计需求的结构.Zhou等[9]采用基于单纯复形的结构拓扑优化方法，使用二阶单元划分结构网格，并选取较高应力指数系数以构建结构最大单元应力近似计算公式[10]，使优化结果的应力水平得以降低.Lian等[11]采用“结构拓扑开洞-结构形状优化-可变复型网格跟新”的迭代循环方法进行考虑结构应力的拓扑优化，优化结果不存在灰度单元，且优化过程中结构应力水平近似计算较准确.

在此基础上，本文针对基于变密度法考虑结构强度拓扑优化问题，采用过滤-投影结构参数化方法，通过研究主要优化参数对结构比强度优化问题优化过程的影响规律，提出一种新的优化策略，实现了在拓优化过程中对结构应力水平及其变化趋势的准确控制，同时保证优化过程的稳定性.最后通过典型结构的算例实验对该方法的合理性及实用性进行验证.

1 考虑结构强度的变密度法拓扑优化 模型

基于变密度法的结构强度拓扑优化研究一直以来是工程应用关注的对象和学界的研究热点.本文基于固体各向同性材料惩罚(SIMP)法，对以结构体积为约束，最小化结构最大应力的拓扑优化问题进行研究[12].其优化列式可表述为

s.t.KU=F

(1)

0≤ρe≤1,e=1,2,…,ne

式中：ρ={ρ1,ρ2,…,ρe}表征各单元的相对密度，其中ρe在0与1之间连续变化，ρe=0和ρe=1分别表示该位置为空单元及实体单元，0<ρe<1则表示此处单元具有中间密度材料；σmax为结构单元最大应力，表征结构整体应力水平，本文采用Von Mises屈服准则对单元应力进行表达(后文亦称σmax为最大单元Von Mises应力)；K、U及F分别为结构整体刚度矩阵、位移向量及外力向量；ve为对应单元的体积；V*为所允许的结构体积上限；ne为结构内单元总数量设计变量场.

(2)

式中：Ωe为单元e附近所有形心落在过滤区域内的单元集合；加权函数ωj为相应单元中心离距的线性函数；ρj为该对应单元的相对密度；单元e的过滤区域为以该单元形心为圆心，过滤半径为r0的圆形区域；rj为Ωe内单元j与单元e的中心离距.

本文采用Von Mises屈服准则进行单元应力计算[13-14]，并构建如下所示的类SIMP法单元应力插值函数，以解决单元应力奇异性问题：

(3)

(4)

e=1,2,…,ne

式中：σ11、σ22及σ12为单元e应力矢量的3个分量.

单元应力是局部物理量，优化过程中需对各单元应力进行控制以实现对结构整体应力水平的预测.而将每个单元应力均作为优化约束时，存在约束过多、难以求解的问题.因此，本文采用如下所示的p-norm法近似计算：

(5)

2 结合投影法的参数化方法

基于传统变密度拓扑优化方法所设计的结构中，中间密度单元使结构应力在对应位置处无法准确计算，而为获得单元密度非0即1的离散结构，一般对优化结果采用如下所示的投影后处理，以去除结构中间密度单元：

(6)

(7)

图1 投影函数曲线Fig.1 Curves of projection function

针对上述具有较高非线性的结构应力优化问题，本文采用移动渐近线(MMA)法[16]作为优化器对其进行寻优.优化器中设计变量变化阈值(ML)表征每次优化迭代中各设计变量的变化范围，亦即优化算法中的搜索步长[12]， 其同p-norm法的p和投影函数的β耦合，综合影响优化过程的稳定性和收敛性.本文基于式(1)的优化模型，通过控制变量法对上述3个主要优化参数的影响规律进行研究.

以图2所示的L型梁作为优化对象，图中R为半径.采用 6 400 个边长为1 mm的一阶四节点四边形单元对其进行网格划分，L型梁左上端固定，为避免产生应力集中，将竖直向下大小为5 N的力分解加载在结构右端扇形区域内的所有节点上，实体材料弹性模量E=1.0 MPa，泊松比ν=0.3，优化过程中结构体积约束值为0.3.

图2 二维L型梁有限元模型(mm)Fig.2 Finite element model of a L-shaped beam (mm)

图3～6所示为ML、p及β取不同值时的结构拓扑、应力分布云图及结构最大应力.结构材料分布图下方的应力云图中，颜色从蓝到红表征单元应力由0逐渐增大至结构最大单元应力值，如图7所示.图3～6中，以相同p和ML作为一组算例，通过对比各组算例的优化结果可知，β越大，优化结果的中间密度单元比例越低，越能避免如图3(q)、3(r)、3(t) 等所示的结构断裂现象.而由各图中的(a)、(e)、(i)、(m)及(q)可知，随着p不断增大，优化结果的最大单元Von Mises应力有逐渐下降的趋势.当p取12及30等较高值时，优化过程受优化问题非线性影响较大，结构最大单元应力值有上升趋势.由此可知，在实际应用中，当p和β取较小值时，优化问题非线性较弱，结构内中间密度单元比例较大，此时设计变量变化阈值ML宜取较大值，以使结构拓扑变化更灵敏，从而更快地获得具有较低应力水平的优化结果.当为保证优化过程中对结构应力水平和最大应力变化趋势的控制精度而采用较高p值，同时为降低优化结果中间密度单元比例采用而较高β值时，应使用较低的设计变量变化阈值ML.

图3 β=2时优化结果Fig.3 Optimized results at β=2

图4 β=4时优化结果Fig.4 Optimized results at β=4

图5 β=8时优化结果Fig.5 Optimized results at β=8

图6 β=16时优化结果Fig.6 Optimized results at β=16

图7 优化结果应力分布图Fig.7 Stress distribution map of optimized results

3 优化策略

基于上一节所研究的主要优化参数对优化过程及优化结果的影响规律，本文提出了先结构拓扑优化，后近似形状优化的优化策略，以同时提升优化过程中结构应力水平的控制精度和优化过程的稳定性，并使最终设计结果的强度性能满足要求.首先在优化过程的拓扑优化阶段，通过采用较小β及p作为初始值，降低优化问题非线性影响，以实现结构拓扑的快速变化.待符合收敛条件后，增大p至30以提高优化过程中结构应力水平控制精度.满足收敛条件后增大β至64以大幅降低结构中间密度单元比例，使结构变化主要发生在边界，以构建近似形状优化的结构优化过程，最终采用后处理投影方法，获得非黑即白的设计结果.

为探明拓扑优化阶段p及整个优化过程中ML的选取方法，同时验证本文优化策略的有效性，以式(1)为优化模型，并以图2所示的L型梁为研究对象进行优化实验.为提高单元应力计算精度，降低非线性对优化问题的影响，采用二阶矩形单元对L型梁结构进行网格划分，网格数量、尺寸、荷载、边界条件及体积约束等设置与前文相同.在不同的初始p值(图中p′)及ML值下，最终设计结果的σmax变化曲线如图8～10所示.

图8 ML=0.05时优化结果的σmaxFig.8 σmax of optimized results at ML=0.05

图9 ML=0.1时优化结果的σmaxFig.9 σmaxof optimized results at ML=0.1

图10 ML=0.2时优化结果的σmaxFig.10 σmax of optimized results at ML=0.2

根据图8～10可知，当p′≤9时，优化结果的最大Von Mises应力较小.当p′>9时，优化问题的非线性程度较高，优化结果的最大Von Mises值较采用低p′值时更高，整体呈增大趋势，且存在较大振荡现象.当ML取0.1时，易于获得具有更低应力水平的最终设计结果，并保证结构拓扑的灵敏变化.优化策略如图11所示.

图11 优化策略流程Fig.11 Flowchart of optimization strategy

为了验证上述方法的合理性与一般性，以图12所示的Portal结构[11]作为验证实验的优化对象.采用边长为1 mm的二阶正方形单元对该结构进行网格划分，结构左下角5个节点的所有自由度固定，右下角5个节点仅固定竖直方向的自由度，大小为5 N的力平均作用在结构顶部中间5 mm范围的节点上.采取与L型梁算例相同优化参数设置构建对比实验，设计结果的最大单元Von Mises应力变化曲线如图13～15所示.

图12 Portal结构尺寸及边界条件(mm)Fig.12 Size and boundary condition of the portal structure (mm)

图13 ML=0.05时优化结果的σmaxFig.13 σmax of optimized results at ML=0.05

根据图中验证实验的优化结果可知，设计结果最大单元应力的变化趋势以及优化参数对设计结果强度的影响规律与L型梁实验一致，该系列优化结果以图16为例，表明上文所总结的优化参数影响规律及提出的参数选取策略具有一定的普适性及实用性，可用于指导工程实际应用.

图14 ML=0.1时优化结果的σmaxFig.14 σmax of optimized results at ML=0.1

图15 ML=0.2时优化结果的σmaxFig.15 σmaxof optimized results at ML=0.2

图16 Portal结构优化结果Fig.16 Optimized results of portal structure

5 结语

本文采用基于过滤-投影的参数化方法，实现在迭代优化过程中结构中间密度单元比例的不断下降，进而降低后处理前后优化结果应力水平的偏差，使其满足工程实际的设计需求.通过对经典L型梁结构的比强度优化实验结果的分析，研究了关键优化参数对优化过程及优化结果应力水平的影响规律，并在此基础上，提出了先结构拓扑优化，后近似形状优化的新优化策略，同时在优化过程中保证了对结构应力水平的控制精度以及优化结果的稳定性.在上述优化策略下对典型Portal结构进行优化设计，并对优化结果的应力水平进行了分析，阐明了该方法的合理性及实用性，可为实际工程中结构比强度拓扑优化设计提供指导.