刘云龙，买买提明•艾尼

（新疆大学机械工程学院，新疆 乌鲁木齐 830047）

1 引言

根据不同的优化手段，可将结构优化设计分为尺寸优化、形状优化和拓扑优化[1]。而拓扑优化则是结构优化设计的最高层次，其主要思想是将寻求结构的最优拓扑问题转化为在给定的设计区域内寻求最优材料的分布问题，以达到结构的最大刚度或最大基频[2]。对于拓扑优化设计，众多的学者都做了许多工作。其中，文献[3]在1904年提出的桁架理论可以认为是拓扑优化设计的起源；随后，文献[4、5]引入了材料“微结构”的概念，指明了连续体结构拓扑优化的发展方向；1988年，文献[6]在微结构概念的基础上，提出了均匀化拓扑优化方法，这标志着连续体拓扑优化的发展进入黄金时期，而当前ANSYS等商业工程软件都有基于均匀化方法的拓扑优化模块。

ANSYS软件具有优良的建模功能和强大的分析能力，是进行拓扑优化设计的有力工具之一。但是基于均匀化方法的拓扑优化结果边缘往往具有棋盘格式的不规则和非线性特点，这种棋盘格式就给拓扑优化结果的CAD模型重构和CAE分析带来一定的困难[7、8]。棋盘格式的出现与所采用的材质设计变量无关，而是与有限元单元有密切关系，也就是说棋盘格式是连续体结构拓扑优化中依赖单元类型和数量的一种固有的现象。因此，在拓扑优化中，不论是采用均匀化设计方法还是采用密度法，只要用有限元法进行拓扑结构优化均会不可避免的出现棋盘格式。

棋盘格式问题的解决一直是很多学者关心的问题之一，并提出了一些切实有效的避免棋盘格式的求解策略。比如，采用后处理技术将棋盘格式过滤的方法，采用较为稳定的有限元模型或改变优化目标函数的泛函从而使优化过程趋于稳定等等。这些方法都能不同程度的缓解或减弱棋盘格式现象，但是不能完全消除，这依然为拓扑优化结构的实际工程应用带来了一定的困难。

这里则应用最小二乘法的线性拟合原理对边缘棋盘格式进行线性光滑处理，并根据边缘光滑处理的结果重构新的光滑拓扑优化结构，从而完全消除了棋盘格式。通过对新的光滑拓扑优化结构模型再进行网格划分和有限元分析验证了此方法的有效性。这种使用最小二乘法对棋盘格式进行线性光滑处理方法，可完全消除棋盘格式并建立基本保持拓扑结构的新的有效拓扑结构模型。本方法理论简单、容易编程实现、建模速度快、易于操作等优点，为拓扑优化设计提供了一种新的手段。

2 多边缘棋盘格式的最小二乘法描述

2.1 基于均匀化方法的拓扑优化方法

均匀化方法是一种经典的拓扑优化方法，有着严密的数学和力学理论基础，由于在基于均匀化方法的拓扑优化方面已有很多的文献和详细的报道，因此在此不再重复[6、9]。这里则是直接利用ANSYS中的基于均匀化方法的拓扑优化方法进行了建模和优化分析，得到了存在棋盘格式的小梁结构边缘。并针对这种具有棋盘格式的边缘，使用最小二乘法进行了光滑处理。

2.2 基于最小二乘法的多边缘棋盘格式光滑处理算法

最小二乘法是一个重要的数学与统计模型，在回归分析、数据拟合、图像处理等方面有着广泛的应用[10]。最小二乘法通过最小化误差的平方和找到一组数据的最佳函数匹配，是用最简单的方法求得一些绝对不可知的真值，而令误差平方之和最小[11]。利用最小二乘法，可以对平面内一系列非线性的散点进行线性或非线性处理。本研究考虑到拓扑优化结构的实际工程应用价值，尽量把小梁结构部分简化为线性的梁，因此本研究用线性的最小二乘法来近似小梁结构边缘棋盘格式的非线性散点（节点），并用一条直线方程y=ax+b来描述这一系列非线性的散点所组成的非线性曲线。对于复杂拓扑结构存在多个独立的小梁结构边缘棋盘格式域，并组成多个独立的散点域，其基于最小二乘法的多边缘棋盘格式光滑处理算法介绍如下。

设在平面直角坐标系内有l个散点域，每个散点域有nl个散点组成的散点集，i=1，2，3…n。在第l个散点域内用一个一次函数方程式来表达这n l个散点，则有：

式中：y l—xl的预估值。令S为这n个散点的实际值与预估值的差的平方和，即：

此时，令S(al，bl)=0，并在S(al，b l)中分别对al和b l求偏导，则有：

联立式（4）、式（5），且把a和b看作是未知数，则式（4）、式（5）的联立式是一个两元一次方程组，求解这个两元一次方程组，得：

令al为斜率，bl为截距，此时可以得到一条平面直角坐标系第l个散点域内的直线。通过求解所有散点域的光滑直线并连接相邻直线后即可得到完全消除棋盘格式的新的拓扑优化结构。

3 多边缘棋盘格式结构的光滑处理

3.1 拓扑优化及其边界条件

在平面上有一宽为3000mm，高为1000mm的板，材料选用45号钢，即弹性模量为2.1×1011，泊松比为0.269，屈服强度355MPa[12]，如图1所示。平面板的左侧受固定约束；上端距左侧1000mm处受一方向垂直向下，大小为1000N的力；下端最右侧处受一方向垂直向上，大小为1000N的力。网格尺寸为50mm，采用四边形网格主导的四边形∕三边形网格。拓扑优化的约束条件设置为减少50%的面积。在上述条件下先在ANSYS Workbench中对模型进行求解。

图1 平面板的边界条件Fig.1 Boundary Conditions of Planar Plates

为了便于拓扑优化结果的提取，建立的项目图表，如图2所示。以便将ANSYS Workbench中的数据转换进入ANSYS APDL。之后右击C2项目图表，选择Edit in Mechanical APDL以进入ANSYS APDL环境。

图2 转换数据所用的项目图表Fig.2 Item Charts Used to Convert Data

进入ANSYS APDL后，利用APDL命令流对模型再进行一次拓扑优化计算，所需要的APDL命令流为：

即可得到约束条件为面积减少50%的拓扑优化结果，如图3所示。其中，红色部分为保留材料的部分，蓝色部分为去除材料的部分。可以看出，拓扑优化结果具有明显的多个带棋盘格式的梁结构。

图3 平面板面积减少50%的拓扑优化结果Fig.3 Topological Optimization Results for 50%Reduction of Planar Plate Area

3.2 拓扑优化结构的光滑处理

由第二节描述的最小二乘法可知，只要得到某条非线性曲线附近若干散点坐标，就可利用最小二乘法求得一条直线方程，并用这条直线方程来近似代替非线性曲线。这里就用以上方法对复杂拓扑优化多小梁边缘棋盘格式结构进行光滑处理。其具体的光滑处理过程如下。

3.2.1 拓扑优化结果的提取

通过ANSYS分析得到的拓扑优化结果，如图4所示。图中可以看出拓扑优化结果包含多个小梁结构及多边缘棋盘格式的结构。通过详细分析多边缘棋盘格式将各个边缘设置了总共22个编号，每一个编号代表一个边缘棋盘格式散点域，整个拓扑优化结构内总共有22个散点域。为了将这些散点域组成的棋盘格式设为一条线进行光滑处理，首先对这些散点域内的棋盘格式点（节点）进行坐标的提取。以图4中的10号棋盘格式散点域为例，利用ANSYS APDL中的Picked Entities+命令对10号棋盘格式散点域内的节点进行拾取，如图5所示。

图4 复杂拓扑优化多小梁结构的边缘棋盘格式编号及散点域Fig.4 Edge Chessboard Format Numbering and Spot Domain of Multitrabecular Structures with Complex Topology Optimization

图5 10号棋盘格式散点域内节点的拾取Fig.5 Picking up Nodes in Dispersion Domain of Chessboard Format No.10

在选择全部节点之后，单击拾取器中的OK即可得到这些节点的坐标。同理，利用相同的方法可以提取其余21条棋盘格式散点域内的节点，在此不一一赘述。

在节点拾取的过程中通过观察坐标可以发现，1、2、5、7号线已经是平面直角坐标系内的直线了，它们分别为：x=-1500、y=1000、x=1500、x=-1449.4。

3.2.2 利用最小二乘法进行光滑处理

由图5（a）拾取器中的信息可知，10号棋盘格式散点域内共有19个节点，它们的坐标分别为：（-160.57，499.72）（-110.25，499.68）（-65.786，449.70）（-116.00，449.88）（-256.12，549.52）（-206.07，549.57）（-155.88，549.57）（-19.418，400.02）（-70.070，399.95）（-302.41，599.44）（-252.38，599.44）（-22.903，350.30）（-350.02，649.36）（-300.03，649.36）（-400.01，699.36）（-350.02，699.36）（-450.01，749.36）（-400.01，749.36）（-450.01，799.36）。

将以上19个点带入式（6）和式（7）得：

a10≈-0.892024562

b10≈364.9228648

以上系数a10和b10可得10号棋盘格式散点域光滑处理后的一条直线方程y10=-0.892024562x+364.9228648。可用这条直线方程来近似代替图4中的10号棋盘格式散点域内散点所组成的曲线。同理，除1、2、5、7号之外，利用相同的方法均可得到代替剩余17个棋盘格式散点域内曲线的近似直线方程。

3.2.3 利用WPS表格简化计算过程

众所周知，WPS的表格软件拥有强大的函数计算功能。在求解图4中每条线的过程中，可以利用WPS表格来计算每条线的斜率al和截距bl。

仍以图4中10号棋盘格式散点域为例，如图6所示。将在ANSYS APDL中提取的19个点的X坐标粘贴进A列，Y坐标粘贴进B列；C2单元格输入=A2*A2，并下拉单元格至C20；D2单元格输入=A2*B2，并下拉单元格至D20；G1单元格输入=SUM（C2：C20）；G2单元格输入=SUM（A2：A20）；G3单元格输入=SUM（D2：D20）；G4单元格输入=SUM（B2：B20）；G5单元格输入=（19*G3-G2*G4）∕（19*G1-G2*G2）；G6单元格输入=（G1*G4-G2*G3）∕（19*G1-G2*G2）。此时，WPS表格软件即可自动计算出10号棋盘格式散点域内所对应的直线斜率al和截距bl，它们分别为G5和G6中的数值，如图6所示。

图6 利用WPS表格软件计算直线斜率和截距Fig.6 Using WPS to Calculate Slope and Interception of Lines

使用相同方法，加之之前已得到的1、2、5、7号直线，可得22条近似直线。其中，4号直线的斜率a4为1.39167×10-33，截距b4为-2.23408×10-30，由于数值非常小，因此可以认为4号直线为：y=0。

综上，可得图4中1、2、4、5、7号线的近似直线方程，如表1所示。而对于图4中剩余17条线的近似直线方程，可以求得其截距al和斜率bl，如表2所示。

表1 1，2，4，5，7号线的近似直线的方程Tab.1 Approximate Linear Equations for Line 1，2，4，5 and 7

表2 剩余17条近似直线的斜率a l和截距b lTab.2 Slope A and Intercept B of the Remaining 17 Approximate lines

4 光滑结果的重构与新拓扑结构分析

根据拓扑优化得到的非线性非光滑结果，将3.2.3中所得的（1～5）号直线方程、（6～10）号直线方程、（11～13）号直线方程、（14～16）号直线方程、（17～19）号直线方程、（20～22）号直线方程分别联立，求得它们的交点，即可得重构线性光滑模型所需的22个关键点，如表3所示。

表3 重构线性光滑模型所需的22个关键点坐标Tab.3 22 Key Points Coordinates for Reconstructing Linear Smooth Model

根据表2中的22个关键点坐标，即可在ANSYS Workbench中对拓扑优化且进行了光滑处理的平面板模型进行重构。对重构的拓扑优化模型进行网格划分，如图7所示。网格尺寸为10mm，采用四边形网格主导的四边形∕三边形网格。

图7 重构的模型Fig.7 Reconstructed Model

将3.1的边界条件再次施加到重构的光滑模型上，模型的左侧受固定约束；上端距左侧1000mm处受一方向垂直向下，大小为1000N的力；下端最右侧处受一方向垂直向上，大小为1000N的力。并设置模型的弹性模量为2.1×1011，泊松比为0.269。通过有限元分析，分别对其进行位移、应力和应变分析，可得的计算结果，如图8～图10所示。

图8 重构模型的位移分析Fig.8 Displacement Analysis of Reconstructed Model

图10 重构模型的应变分析Fig.10 Strain Analysis of Reconstructed Model

由图8可知，重构模型的最大变形量为0.69344mm，最大变形分布在重构拓扑优化模型的最右端。但是，从图9和图10中可以看出，应变和应力都分布的比较均匀，基本消除了应力集中。最大应变为0.00020234，非常小，最大应变发生在重构模型中17、18号直线的交点处，如图9所示。同样，最大应力为38.078MPa，最大应力也发生在重构模型中17、18号直线的交点处，这是由于交接点出现锐角而未进行圆滑处理的缘故，虽然这样最大应力仍比45号钢的屈服强度355MPa小得多，如图10所示。综上，对利用最小二乘法光滑处理后所得到的重构拓扑优化模型进行有限元分析结果表明，应力应变分布均匀，表明最小二乘法光滑处理后重构所得到的拓扑优化模型基本保持了原拓扑结构，同时验证了本方法的有效性。

图9 重构模型的应力分析Fig.9 Stress Analysis of Reconstructed Model

5 结论

（1）建立了多小梁边缘棋盘格式复杂拓扑结构的最小二乘法光滑处理算法，并利用WPS对拓扑优化的结果进行计算，实现了光滑处理，为实际工程应用提供了一种简便且易于操作的多小梁边缘棋盘格式的复杂拓扑结构光滑处理方法；（2）通过AN⁃SYS Workbench和ANSYS APDL软件实现了拓扑优化结果的提取，并将光滑处理后的拓扑优化结果在ANSYS Workbench中进行了重构，实现了完全消除棋盘格式结构的光滑拓扑优化结构，并对其进行了相应的有限元分析，验证了本方法的有效性；（3）通过新光滑拓扑优化结构有限元分析可知，在给定的边界条件下和面积减少50%的拓扑优化约束条件下，利用最小二乘法进行边界光滑处理得到的重构拓扑优化模型的应力应变分布均匀，消除了应力集中，光滑处理后的重构拓扑优化模型基本保持了原拓扑结构形状。