限界法克服连续体结构拓扑优化载荷病态问题

俞燎宏，荣见华

（1.宜春学院物理科学与工程技术学院，江西 宜春 336000；2.长沙理工大学汽车与机械工程学院，湖南 长沙410076；3.工程车辆轻量化与可靠性技术湖南省高校重点实验室，湖南 长沙 410076）

1 引言

自文献[1]提出均匀化方法以后，连续体结构拓扑优化方法迅速发展，出现了实体各向同性材料惩罚法（SIMP）[2]、渐进结构优化法（ESO/BESO）[3-4]、水平集法（Level Set Method）[5-7]、独立连续映射法（ICM）[8]和可移动变形组件法（MMC）[9]等。其中，SIMP方法目前应用最为广泛。在实际工程中，很多结构通常受到多种工况载荷的作用。比如主要承载机械的结构上可能需要安装传感器或操作平台等，故该结构除了强载荷作用外，还在偏离强载荷作用部位的某处存在小载荷作用。要想在优化过程中获得性能良好的拓扑结构，必须同时考虑所有的载荷工况。

由于不同工况载荷的大小、作用方向和作用位置不同，对结构的作用效果也不尽相同。那么，在基于传统优化方法的多工况载荷作用下结构拓扑优化过程中，可能出现没有材料支承小载荷的优化结构，这就是“载荷病态”问题。为了解决载荷病态问题，文献[10]提出了分层优化方法；文献[11]提出了ICM 应力全局化方法；文献[12]提出了分数模常法；文献[13]对上述方法进行了综合分析和评判，并提出了敏度分层过滤解决载荷病态的策略。该策略能够得到较为清晰的结构，可以克服荷载病态问题，但是优化结果依赖于两个关键系数。

多工况载荷作用下的结构柔顺度拓扑优化问题本质上是一个多目标优化问题，关于多目标优化问题的研究综述可见文献[14]。针对以多工况载荷下结构柔顺度最小为目标函数、结构体积为约束条件的优化问题，文献[15]采用限界法提出了一种新的求解思路，能以较少的工作量寻求更优的结构拓扑。限界法[16-17]非常适合于不可微的最小最大优化问题。文献[18]将限界法应用到振动连续体结构多特征频率和频率间隙最大化的拓扑优化设计。文献[19]利用限界法开展了制造约束下复合材料框架最大基频离散选材与结构拓扑优化。

经过对现有载荷病态处理方法分析和研究的基础上，提出采用限界法来克服载荷病态问题。首先，给出以最大的结构柔顺度最小化为目标，以体积为约束的拓扑优化模型。其次，引入限界变量将多个柔顺度目标函数转化为约束条件，结合变体积约束限技术，建立新的等效近似优化模型；针对同一工况内含有病态载荷的问题，在新建的优化模型中增加柔顺度小量变化约束。然后，对目标函数和约束函数进行灵敏度计算，采用移动渐近线方法（MMA）优化求解。最后，给出了二维和三维算例，验证了所提方法的可行性和有效性。

2 多工况载荷下的结构拓扑优化模型

2.1 材料插值与Heaviside过滤

通常，基于SIMP材料插值模型和Heaviside过滤技术，将结构的设计域离散成N个单元，每个单元赋予一个设计变量和一个相对密度变量。其中，第e号单元的相对密度变量x͂e(e=1，2，…，N)与它周围单元设计变量xi(i=1，2，…，N)之间存在如下关系：

当与Heaviside过滤技术一起使用时，常数加权可以更有效地减少了灰度，而线性加权则倾向于生成更平滑的拓扑构型。在文中采用常数加权形式。

为了获得0/1分布清晰的结构拓扑，采用文献[20]提出的光滑Heaviside函数，其表达式为：

式中：xˉe—单元物理变量；β—曲率参数。

通过式（1）和式（5）形成设计变量场x、过度变量场x͂和物理变量场xˉ。结构有限元分析时，用物理变量xˉe替代相对密度变量x͂e。

结构的单元体积和刚度矩阵通过公式（5）来获得：

式中：Ve和Κe—第e号单元的体积和刚度矩阵；—第e号单元的固有刚度矩阵单元体积的惩罚函数文中选取刚度矩阵的惩罚函数其中，E0—实体材料的弹性模量；Emin—空洞材料的弹性模量，为了避免刚度矩阵奇异，这里选取Emin=10-9E0；p—惩罚因子，文中选取p=3。

2.2 多工况载荷下结构柔顺度拓扑优化模型

在实际工程中，结构通常在复杂的工况载荷下工作，不同工况载荷对结构的要求也有差别。在多工况载荷作用下，以结构体积为约束条件，以最大的结构柔顺度最小化为目标的拓扑优化模型可以表示为：

式中：L—工况载荷总数；Fl—第l号工况载荷向量；Cl和Ul—第l号工况载荷作用下结构的柔顺度和位移向量；V0—整个结构的初始体积；f*—预先给定的体积比。

3 利用限界法克服载荷病态问题

3.1 基于限界法的拓扑优化模型

参考文献[14-16]的限界法，引入一个限界变量z(z＞0)将优化模型式（6）中的多个目标函数转化为约束条件，即将多目标优化问题转换成单目标优化问题。为了协调Heaviside变换引起的结构体积冲突，同时确保每一迭代步的结构特性函数稳健变化，引入变体积限约束[21]。其转换后的近似等效优化模型为：

式中：Cmax—所有工况载荷作用下的最大结构柔顺度；VU(k)—设定的当前步结构体积限，其更新表达式为：

式中：V*—目标体积；V(k-1)—上一步得到的结构总体积；γ—体积限的经验参数，其取值范围为[0.005，0.025]；Q—Heaviside映射曲率参数β变化后的前三个迭代步的迭代步编号的集合。

3.2 柔顺度小量变化约束策略

文献[11]将载荷病态分为三类：第一类是多工况间有载荷病态，但工况内无载荷病态；第二类是仅在工况内部有载荷病态；第三类是多工况间有载荷病态，同时某工况内也有载荷病态。

为了便于描述，将引起载荷病态问题的工况载荷称为“病态载荷”。通过数值算例表明，优化模型式（7）可以解决第一类载荷病态问题，但是不能解决同一工况内有载荷病态的问题，即不能解决第二类和第三类载荷病态问题。

为了解决同一工况内有载荷病态的问题，借鉴文献[21]提出的方法：将柔顺度小量变化要求作为约束条件，建立新的近似优化模型。文中将其称为柔顺度小量变化约束策略，其实现过程如下：

（1）将同一工况内的病态载荷单独取出来作为新的虚载荷，并将其独立加载在结构上。

（2）将虚载荷作用下的结构柔顺度按照表达式（10）形成新的约束条件。

（3）将柔顺度小量变化约束表达式（10）加入优化模型式（7），形成新的优化模型式（11）。

据董事长王茂莆介绍，公司创业之初，国内喷灌机产业面临外资机构大规模并购的尴尬局面，但为了实现行业民族品牌的发展，华雨公司迎难而上，凭借“踏实、拼搏、责任”的企业精神和“诚信、共赢、开创”的经营理念，通过不懈努力，终于生产出具有我国自主知识产权的大型喷灌机设备，为赢得国产喷灌机市场定价权起到了积极的促进作用。王茂莆说，我国灌溉产业是个存量很大的市场，随着十八届三中全会的召开，“调整产业结构，转变发展方式”已经被纳入国家体制改革的宏观计划，未来新兴灌溉产业将会有更大的发展空间，因此把握产业命脉的关键是打好基本功，夯实自身优势。

如果工况内有多个病态载荷，则添加对应的多个柔顺度小量变化约束表达式。

4 灵敏度计算及优化求解

4.1 目标函数和约束函数的灵敏度计算

在优化模型式（11）中，目标函数fo(x)对物理变量的导数为：

结构柔顺度Cl(x)对物理变量xˉe的导数为：

式中：Ul，e—仅与结构第e号单元自由度相关的Ul的分量构成的矢量。

柔顺度小量变化约束函数ϕs(x)对物理变量的导数为：式中：—第s个病态载荷作用下的位移向量仅与结构第i号单元自由度相关的的分量构成的矢量。

结构体积V(x)对物理变量的导数为：

由前文可知，结构的任意性能函数f(x)对设计变量xi的导数为：

4.2 MMA优化求解

优化模型式（7）和式（11）采用移动渐近线方法（MMA）求解，具体求解过程参见文献[22]。对于载荷病态优化问题，由于不同工况载荷的作用效果相差较大，优化过程中结构柔顺度会发生剧烈振荡。为了抑制振荡行为，同时协调Heaviside变换引起的急剧变化，参考文献[23]，对MMA算法进行适当的修改。（1）修改前两次迭代上、下移动渐近线表达式（19）中的移动渐近系数s0，修改后的s0如式（20）；（2）将式（21）引入到MMA算法中，通过和xml来控制设计变量xi的变化。

式中：曲率参数β的初始值β(0)=1，每迭代50 步按β(m)=1.5*β(m-1)进行更新。

5 优化算例

5.1 算例1

文献[11]采用的基结构，尺寸为2000mm×400mm，厚度为9mm，弹性模量E=68.89GPa，泊松比v=0.3，左右两端由铰链支承，如图1所示。受到两组工况载荷作用，工况1：两个集中载荷=200N分别作用于上边界距离梁端100mm处的A点和C点；工况2：一个集中载荷F2=20000N 作用于梁上边界中点B。结构划分为100×20个矩形网格单元，体积比f*=0.5。

图1 基结构Fig.1 Basic Structure

本算例工况间有载荷病态但工况内无载荷病态，属于第一类载荷病态问题。本算例过滤半径rmin=2.0Δmi（nΔmin是网格单元边长的最小尺寸），收敛经验参数ε=0.002，采用所提方法获得的优化历程和优化结果，如图2所示。文献[11]方法获得的优化结果，如图3所示。

图2 优化历程和优化结果Fig.2 Optimization Process and Results

图3 文献[11]方法获得的优化结果Fig.3 Optimization Results Obtained by Reference[11]

其中图3（a）为反演前的拓扑值分布图（红色单元的拓扑值为1，其它颜色单元拓扑值介于0与1之间），图3（b）为反演后的最优拓扑图形。从图2和图3可见，优化结果整体构型基本相似，但支撑病态载荷传力路径的构件有所不同。采用所提方法不需要进行反演就能够获得0/1分布清晰的拓扑结构，如图2（d）所示。离散度指标md=2.63%，病态载荷的传力路径清晰。所提方法和文献[11]方法优化得到最大工况（工况2）的最小柔顺度（最小应变能）分别为20.04N·m和22.37N·m。因此，所提方法可获得柔顺度更小的优化结果，更符合柔顺度最小化的优化目标。该算例表明，对于工况间有载荷病态同时工况内无载荷病态的拓扑优化问题，采用所提方法是可行且有效的。

5.2 算例2

文献[13]采用的初始结构，为300mm×200mm的长方形平面区域，左右两边受固定支撑，如图4所示。单工况内受两个载荷作用，集中载荷F1=1000N作用于上边界中点A，方向向下；集中载荷F2=1N作用于下边界中点B，方向向上。结构划分为60×40个矩形网格单元，体积比f*=0.5。本算例仅在工况内有载荷病态，属于第二类载荷病态问题。对于这一类问题，将病态载荷F2作为虚载荷，并独立加载到结构中的B点。本算例中柔顺度小量变化经验参数ϕ*=0.002，采用所提方法获得的优化结果，如图5所示。如果初始结构仅受到集中载荷F1的作用，该问题不属于载荷病态问题，优化获得的结果，如图6所示。

图4 初始结构Fig.4 Initial Structure

图5 优化结果Fig.5 Optimization Result

图6 仅载荷F1作用下的优化结果Fig.6 Optimization Result Only Under F1

由图5和图6可知，优化结果中集中载荷F1的传力路径基本相似，图5仅比图6多了病态载荷F2的传力路径。可见，对于仅在工况内有载荷病态的问题，按提出的柔顺度小量变化约束策略不会改变主要载荷的传力路径，而且得到了清晰的病态载荷传力路径，表明该方法是可行的。

采用文献[13]方法，不同参数a和b对应的优化结果，如图7所示。对应的最小柔顺度分别为4975.556N·m、4993.970N·m、4940.627N·m和4966.611N·m。由图5和图7可知，载荷F1的传力路径基本相似，载荷F2的传力路径略有差异；图5结构0/1分布更加清晰，离散度指标md=0.15%。采用所提方法计算得到载荷F1作用下的最小柔顺度为4739.244N·m，比文献[13]方法计算得到的柔顺度值更小。该算例表明，对于同一工况内有载荷病态的拓扑优化问题，采用所提方法是可行且有效的。

图7 文献[13]方法优化结果Fig.7 Optimization Results Obtained by Reference[13]

5.3 算例3

初始结构为600mm×400mm×40mm 的长方体三维区域，左右两端面受固定支撑，如图8所示。

图8 初始结构Fig.8 Initial Structure

受到2组工况载荷作用，工况1：集中载荷F1=10000N竖直向下作用于结构顶面中心A点，F′1=10N竖直向上作用于结构底面中心B点；工况2：集中载荷F2′=F2′′=10N竖直向下分别作用于结构顶面C点和D点。划分60×40×4个矩形网格单元，体积比f*=0.4。

本算例工况间有载荷病态同时某工况内也有载荷病态，属于第三类载荷病态问题。对于这一类问题，将工况1中的病态载荷F1′作为虚载荷，并独立加载到结构中的B点，柔顺度小量变化经验参数ϕ*=0.005。采用所提方法获得的优化结果，如图9 所示。其中，X-Z平面视图，如图9（a）所示；30°三维视图，如图9（b）所示。由图9 可知，优化结果单元0/1 分布清晰，离散度指标md=3.7%，病态载荷F′1、F2′和F2′′的传力路径被保留并清晰展现。可见，按照提出的方法能很好地克服工况间有载荷病态同时某工况内也有载荷病态的拓扑优化问题，即能解决第三类载荷病态问题。该算例也表明，所提方法在三维空间结构优化中也是可行且有效的。

图9 三维结构优化结果Fig.9 Optimization Results of 3D Structure

6 结论

针对连续体结构拓扑优化中的载荷病态问题，通过限界法结合柔顺度小量变化约束策略和变体积约束限技术，基于移动渐近线方法（MMA），提出了一种克服载荷病态问题的新的优化求解方法。数值算例结果表明：

（1）所提方法可以清晰的展现病态载荷的传力路径，能够有效克服三种不同类型的载荷病态问题；（2）在相同优化条件下，与文献[11-13]方法相比，所提方法可高效地获得0/1分布更清晰、性能更优的结构拓扑；（3）该方法不仅适用于二维平面结构，同时也适用于三维空间结构，为载荷病态问题的优化求解提供了新的思路。