靳立涵JIN Li-han

（重庆交通大学机电与车辆工程学院，重庆 400074）

0 引言

结构优化共有三种优化方式分别为：尺寸优化、形状优化和拓扑优化。其中拓扑优化是一种在有限的空间中根据设计要求确定材料最优布局的数学方法。近三十年来，随着计算机技术的发展，拓扑优化技术也取得了巨大的进步。相关的从业人员和研究者提出了许多的拓扑优化算法。

拓扑优化算法中有三个主流方法分别是：变密度方法、双向渐进结构法[1-2]和水平集方法[3]。对于变密度法，Sigmund[4]在MATLAB 平台上介绍了基于固体各向同性材料惩罚（SIMP）的99 行二维拓扑优化程序，简洁易读，很好地阐述了密度法。Andreassen[5]等在99 行代码的基础上提出了88 行拓扑优化代码，88 行算法在优化相同算例时比99 行算法快100 倍。Liu 和Tovar[6]还在MATLAB 平台上公开了一个169 行拓扑优化程序，在保证代码简单和优化效率高的同时实现了三维目标的优化。Ferrari 等人[7]提出了一种新的99 行和125 行拓扑优化代码，分别处理二维和三维拓扑优化问题，是目前MATLAB 平台上最高效的密度法拓扑优化代码。

本文将MATLAB 与ABAQUS 脚本接口进行集成用以拓扑优化，该集成方法的目的是充分利用MATLAB 平台的矩阵计算能力和商业有限元软件友好的可视化界面。该方法通过一段Python 代码和一段MATLAB 代码实现。这两个代码基于变密度法开发，通过MATLAB 程序调用有商业限元软件实现了对目标设计域在体积约束下的最小合规性的拓扑优化。

1 SIMP 拓扑优化方法理论

本文关于ABAQUS-MATLAB 拓扑优化方法研究的目的，为了得到一个结合两个软件优点并且能适用于基于FEM（有限元法）的多种拓扑优化方法的集成框架。基于变密度方法（SIMP）的特点：发展时间长较为成熟、易拓展和灰度显示。因此本文选择变密度法作为集成框架的理论基础。

SIMP 法拓扑优化模型：

为了将离散型问题转化为连续性优化问题，变密度法提出了一种引入中间密度单元的方法进行计算。但在实际制造过程中，一个单元内只有填充材料和不填充材料两种状态，中间密度的存在使得模型无法制造。为此变密度法还提出了一个通过对中间密度单元进行惩罚的方式来避免产生中间密度单元，使中间密度单元在迭代计算后趋向于实体单元或非实体单元。

在本文中讨论的SIMP 表达式为：

其中xi是编号为i 的单元的相对密度。

SIMP 法优化的数学模型为：

本方案是求解在体积或者质量约束下的优化模型的最小柔度，最终得到在约束条件下的最大刚度的优化模型。其中X 是单元相对密度矢量为单元设计变量，C 是结构柔度，F 和U 分别是载荷矢量和位移矢量，ki为单元刚度矩阵，k0为初始单元刚度矩阵，vi为单元体积，f 是保留的体积分数，V0是初始体积，xmin是设计变量的下限值，xmax是设计单元的上限值，m 为最大单元数。

2 ABAQUS-MATLAB 集成方法及实现

图1 是ABAQUS-MATLAB 集成框架示意图，其中MATLAB 程序完成拓扑优化模型求解，包括灵敏度过滤和设计变量更新；Python 程序由MATLAB 调用后会执行ABAQUS 有限元计算内核，实现每个迭代步骤中有限元求解。在运行之前，CAE 文件、MATLAB 代码和Python 代码需要放在同一个文件夹中，其余文件由程序生成。

图1 ABAQUS-MATLAB 集成方案

根据结构拓扑优化设计要求，在ABAQUS/CAE GUI中建立有限元模型，模型包括网格单元、材料属性、边界条件和载荷应用，其中材料属性可以被更新。MATLAB 代码完成求解初始化，包括有限元模型的单元总数和迭代的相关参数。随后，在循环迭代开始时进行拓扑优化，Python 程序读取MATLAB 代码生成的各单元设计变量，根据（1）式的材料插值模型创建各单元对应的杨氏模量，并更新有限元模型的材料属性。

Python 程序调用ABAQUS 内核，求解有限元模型，并从结果文件中提取相关参数。MATLAB 代码根据优化相关参数完成灵敏度过滤和设计变量更新，在获取新的灵敏度后根据收敛准则判断是否收敛，如果不收敛则进行下一次循环迭代，直到收敛。*.XLS 文件是MATLAB 和ABAQUS之间传输数据的介质，该文件是通过代码调用Microsoft Excel 生成的。

2.1 MATLAB 程序的实现

本文所编写的MATLAB 程序结构如图2 所示。当开始优化时，初始化设计模块对设计域内的设计变量和优化参数进行初始化。初始化后的设计变量经设计变量输出模块以电子表格的形式输出到工作目录。使用系统命令调用Python 程序执行有限元分析。待有限元分析结束后使用有限元结果处理模块将有限元结果组装成灵敏度信息。使用网格过滤子程序对灵敏度信息进行过滤。过滤后的灵敏度信息经过设计变量更新子程序的处理得到新的设计变量。将新旧设计变量通过收敛准则计算是否收敛，如果收敛则输出优化结果否则程序就回到设计变量输出模块继续执行。

图2 MATLAB 程序结构

2.2 Python 程序的实现

Python 程序的结构如图3 所示。该程序可以分为有限元前处理、有限元分析和有限元后处理三部分。当接到MATLAB 程序发出的系统命令时，该程序开始执行。在前处理阶段程序的操作对象为.CAE 文件，使用MATLAB 输出的设计变量对模型中所有单元的材料属性进行更新。使用有限元分析模块将更新后的模型提交有限元分析并等待分析结束。分析结束后使用有限元后处理模块，将有限元分析结果中的弹性应变能以电子表格的形式输出到根目录。如果处于第一次迭代中Python 程序还会输出模型中所有单元的中心坐标。

图3 Python 程序结构

3 优化算例验证

根据图4（a）中模型示意图，创建一个长度为60 宽为20 的二维悬臂梁结构。梁的左侧施加对称约束而梁的右下角施加完全固定约束。在梁的右上角施加一个方向向下大小为1 的力。在材料方面，创建的材料属性中将弹性属性中的杨氏模量值设定为1 泊松比为0.33。在划分网格时，设定每个网格为网格尺寸的边长为1 形状为四边形的四节点网格，因此在梁中共有1200 个网格单元如图4（b）中所示。完成上述步骤后，将该模型的名称、工作目录、模型数据库名称、网格单元数等信息输入拓扑优化集成程序中开始进行拓扑优化，进过68 次迭代后最终的结果如图4（c）所示。

图4 二维算例的设计域与优化结果

A-M 平台与单ABAQUS 平台运算对比：

在经过修改后，上述ABAQUS-MATLAB 拓扑优化平台（简称A-M 平台）可以实现双向渐进结构法进行拓扑优化。有学者曾使用Python 程序在ABAQUS 平台中实现了双向渐进结构优化方法[8]。在同一硬件条件下，使用同一优化对象与相同的优化参数对两个平台的运算时间进行了对比，所得到优化结果如图5 所示，其中（a）为A-M 平台的优化结果（b）则是ABAQUS 平台的优化结果。

图5 A-M 平台与ABAQUS 平台悬臂梁优化结果

从表1 中可以直观地看到两个平台在同一条件下，对同一对象优化所耗费的时间。在总用时方面A-M 平台远小于ABAQUS 平台，反映在单次迭代的平均用时上也是A-M 平台耗时最少约为ABAQUS 平台的三分之一左右。造成这种结果的原因是因为ABAQUS 平台的计算是基于Python 语言，由于该语言在处理矩阵计算时需要完全遍历的特点导致在网格过滤阶段ABAQUS 平台的耗时要远远大于A-M 平台。

表1 两个平台各项用时对比

4 轮毂的拓扑优化

使用A-M 集成优化方法对轮毂进行拓扑优化。如图6（a）所示优化前的轮毂结构由白色的内外环非设计区域和中间青色的设计区域组成，其中设计区域使用铝合金材料非设计区域使用刚体建模，各部件的尺寸在图中标注。

图6 轮毂结构优化示意图

围绕外环施加四个等大的切向力，内环内侧则完全固定，三个部件之间使用Tie 约束。设计区域的网格单元共7360 个如图6（b）所示。优化目标体积分数为20%，过滤半径为2.0。经过159 次迭代之后得到了图6（c）所示的优化结果。该优化结果满足预设的体积约束，材料分布在集中力施加点与内环周围且以中轴线对称。

5 结论

提出和实现了一种基于变密度方法的ABAQUSMATLAB 集成平台。通过经典算例验证了该集成平台的可行性。经过修改后该平台可以用双向渐进结构法进行拓扑优化。与另一种拓扑优化平台对比得到了本集成平台耗时更少的结果，证明了A-M 集成方法的高效性。使用该集成方法以轮毂为对象进行拓扑优化，得到了结构合理造型美观的优化结果。