基于混合元胞自动机法的结构拓扑优化软件平台开发

Software Platform Development for Structural Topology Optimization Based on Hybrid Cellular Automata

摘要: 元胞自动机结构拓扑优化的理论及其方法不断成熟和完善，逐渐应用于实际工程中。本文阐述了混合元胞自动机法的原理，分析了软件开发的技术架构和开发环境，并基于MATLAB与ANSYS等平台开发出面向产品优化的混合元胞自动机法结构拓扑优化设计软件平台，以降低设计人员采用混合元胞自动机法进行结构拓扑优化的难度。最后，通过经典的二维数值算例，验证软件的有效性。

Abstract: In this paper, using the excitation and de-excitation by impact of ion and free electron, we analyze the impact equilibrium condition of The theory and method of structural topology optimization of cellular automata are maturing and perfecting, and are gradually applied in practical engineering. This paper expounds the principle of hybrid cellular automata method, analyses the technical framework and development environment of software development, and develops a software platform for structural topology optimization design of hybrid cellular automata method based on MATLAB and ANSYS, so as to reduce the difficulty of structural topology optimization by designers using hybrid cellular automata method. Finally, the validity of the software is verified by a classical two-dimensional numerical example.

1. 引言

混合元胞自动机法(Hybrid Cellular Automaton Method，简称HCAM) [1] [2] [3] 是Tovar A.等受到骨骼结构重塑过程的启发，耦合了元胞自动机(Cellular Automaton Method，简称CA)和FEM (Finite Element Method，简称FEM)提出的一种拓扑优化仿生设计方法，是一种无梯度优化算法，成功解决了CA法没有解决场状态的全局分析问题，为连续体拓扑优化提供了新的研究方向。国内外学者对基于HCA的拓扑优化进行了研究。Zakhama等人 [4] 提出了二维连续介质拓扑优化问题的多网格加速计算方法。郭连水等 [5] 针对大变形非线性结构拓扑优化问题，提出了基于HCA法的多空间域连续体结构拓扑优化方法，之后还将基于应变、多域的拓扑优化算法引入混合元胞自动机法中，进行了大变形情况下的车辆防撞性研究 [6] 。杜义贤等 [7] 根据反馈控制原理，基于元胞自动机(CA)理论更新与拓扑设计相关的设计变量，从而抑制连续拓扑优化中的数值不稳定性。张一雨等 [8] 基于HCA算法的动态拓扑优化方法对叶片实现进一步减重，并验证了该动态优化方法的有效性。

与传统的方法相比，混合元胞自动机法不仅避免了在设计过程中目标函数敏度信息的计算，极大地提高了计算效率，还有效抑制了数值不稳定的问题，如棋盘格问题和网格依赖性问题；不仅适用于二维连续体结构拓扑优化，而且还可以直接应用于三维连续体结构拓扑优化中 [9] 。但现有的算法流程复杂，对仿真操作有较高的要求，增加了拓展实际工程运用领域的难度。因此，开发通用的混合元胞自动机法结构拓扑优化设计软件具有重要意义。本文基于MATLAB与ANSYS等平台，开发混合元胞自动机法结构拓扑优化设计软件，并完成案例分析以验证软件的有效性。

2. 混合元胞自动机法基本原理简介

组成HCA模型的元胞栅格是离散的规则网格单元，每个元胞的状态由当前时刻t状态及其邻域元胞的状态决定下一个时刻t +1该元胞的状态。元胞第J种的状态更新规则可以表示为：

(1)

式中： $\Delta N$ 为元胞的邻域元胞个数， $a_i^2\left(t\right),\cdots ,a_i^J\left(t\right)$ 为元胞i的邻域元胞， $R_i^K$ 为元胞i的局部更新规则，对于所有位置的元胞更新规则都相同。

使用HCAM进行拓扑优化迭代的过程中，每个元胞的设计变量 $x_i\left(t\right)$ 可定义为单元密度、几何尺寸、弹性模量等，每个元胞的状态场变量 $$ 可以定义为应力、应变、应变能、或者它们的函数。这时，元胞的状态可以用式(2)表示 [10] ：

$a_i\left(t\right)=\left[\begin{array}{l}{x}_i\left(t\right)\\ S_i\left(t\right)\end{array}\right]$ (2)

设计变量根据材料模型确定，本文将单元相对密度作为设计变量，场变量 $S_i\left(t\right)$ 的定义取决于要解决的优化问题。在一离散区域，结构整体的应变能可以表示为：

$U={\displaystyle \underset{i=1}{\overset{N}{\sum }}{U}_i{v}_i}$ (3)

式中： $U_i$ 为CA中元胞单元i内由于变形而储存的应变能密度， $v_i$ 为元胞单元i的体积。

$U_i=\frac{1}{2}{\sigma }^{\text{T}}\epsilon$ (4)

式中： $\epsilon ={\left[{\epsilon }_x,{\epsilon }_y,{\epsilon }_z,{\gamma }_{xy},{\gamma }_{yz},{\gamma }_{zx}\right]}^{\text{T}}$ 、 $\sigma ={\left[{\sigma }_x,{\sigma }_y,{\sigma }_z,{\tau }_{xy},{\tau }_{xz},{\tau }_{yz}\right]}^{\text{T}}$ ， ${\epsilon }_x$ 、 ${\epsilon }_y$ 、 ${\epsilon }_z$ 、 ${\gamma }_{xy}$ 、 ${\gamma }_{yz}$ 、 ${\gamma }_{xz}$ 分别为节点的x向应变、y向应变、z向应变、切应变， ${\sigma }_x$ 、 ${\sigma }_y$ 、 ${\sigma }_z$ 分别为x向、y向、z向正应力， ${\tau }_{xy}$ Math_35##Math_36#为

xy，yz，xz三个面上的剪应力。

应用HCAM进行结构拓扑优化设计，应变能密度 $U_i\left(t\right)$ 的定义可由具体优化目标函数确定，优化的目的是使元胞的应变能密度均值 ${\bar{U}}_i$ 与状态场变量设定值的差值最小化。因此，优化问题可以表示为：

$\{\begin{array}{l}\min \left|e_i\right|=\left|{\bar{U}}_i-U^{*}\right|\\ \text{s}\text{.t}\text{.}0<x_{\min }\le x_i\le 1,i=1,2,\cdots ,N\end{array}$ (5)

式中：e_i为误差信号，U*为应变能密度的设定值。

式中有下限x_min的要求，从而可以避免有限元分析过程中产生奇异矩阵，本文中取x_min = 10⁻³。元胞单元i的应变能密度的均值 ${\bar{U}}_i\left(t\right)$ 可由元胞及其邻域元胞的状态场变量确定。设计域中的材料分布由局部控制规则确定，它是HCAM算法的基础。两位置控制、线性控制、积分控制和派生控制等局部控制规则已经在混合元胞自动机算法中实现，用来寻求场变量设定值与场变量平均值之间的差值减小到最小。本文利用分段常数函数表示设计变量的变化量，采用了两位置控制策略。为抑制优化中的数值不稳定性，应设置最大变化量限制每次迭代中设计变量的变化。

3. 混合元胞自动机法结构拓扑优化设计软件开发

为规范设计流程，也给使用HCAM开展结构拓扑优化设计的研究人员在改进以及优化设计上提供相关的参考信息，在HCAM结构拓扑优化设计生产中实现应用，进一步优化结构设计及产品的更新换代，开发了HCAM结构拓扑优化设计软件。

3.1. 软件分析流程

在MATLAB图形用户界面(GUIDE)调用APDL文件建立通用的有限元模型并进行求解计算。在建立有限元模型之前，在MATLAB图形用户界面定义模型的设计区域、材料性能、载荷、边界条件和网格密度等建模分析中使用的参数，并在建模分析中用参数名称表示相应的数据。MATLAB图形用户界面接受用户输入的参数，启动ANSYS批处理程序，调用APDL文件模块执行建模与分析过程，MATLAB程序将ANSYS对基结构进行静力分析后生成的各种结果输出到图形用户界面供用户启动后处理程序查看。ANSYS分析得到的结果数据存放入平台的数据库管理系统中，可进行动态更新、查询等功能。平台工作流程图如图1所示。

Figure 1. Platform flowchart

图1. 平台工作流程图

3.2. 技术构架

本文采用MATLAB的程序设计语言进行平台交互界面开发，将复杂的ANSYS命令和操作过程封装在MATLAB编程语言的集成开发环境中，用户只需在主界面中选择子界面并在子界面中输入各种参数，通过参数化建模与分析接口程序，平台程序可以调用后端的ANSYS命令进行分析计算。分析完成后，用户可以在主界面直接选择查看各种结果，方便非工程人员使用。运行ANSYS程序得到的分析结果将作为单独的文件存在。同时为利于管理，将数据放入数据库中，用MATLAB编程结合SQL语言开发数据库管理系统，方便用户对所有结果数据(包括数值，动画和图片)的查询、显示、更新等管理操作。

3.3. 开发环境

本软件平台以MATLAB为工具，运用APDL语言在Windows环境下编写和调试，把ANSYS软件分析过程封装在后台，并将ANSYS有限元分析作为子模块。在CA的优化程序中调用该模块对结构进行有限元计算。编译生成的目标文件可以直接下载到当地主机的存储空间，启动软件后用户只需要输入相关参数即可获得相应的优化分析结果，并可以对结果数据进行查询管理。

3.4. 软件主界面

HCAM结构拓扑优化设计软件界面主要包括参数子界面选择栏、结果显示选择栏、结果显示区和工作区等，结果显示区显示用户选择查看的结果，平台主界面、材料参数和边界条件设置界面如图2所示。

3.5. 软件功能实现

HCAM结构拓扑优化设计软件平台的模块可划分为四个部分：一是参数输入模块，包括样机结构参数输入、环境参数输入、边界条件输入、MATLAB文件调用等；二是分析模块，相应的设置完成后调用后台ANSYS程序进行优化分析计算；三是结果查询模块，进行结果的查询分析；四是数据管理模块，进行所有分析数据包括结果等在内的管理和更新。

Figure 2. Platform interface

图2. 平台界面

4. 案例分析

本文以一个简单的悬臂结构为例，验证混合元胞自动机法结构拓扑优化设计软件的功能。在参数输入模块设置平面板的基本结构为60 mm × 30 mm，厚度为6.0 mm，材料弹性模量E₀为210 GPa，泊松比μ为0.25，密度ρ为 ${\text{7850kg/m}}^{\text{3}}$ 。在边界条件子界面中设置结构左端边界固定约束，在右端自由边界中心位置处，施加向下的集中荷载 $F=\text{2}000\text{N}$ ，如图3所示，求其重量最轻的拓扑。在分析设置子界面设置FEM计算时选择高阶2维8节点的PLANE183单元，邻胞选择von Neumann型(邻胞个数为4个)，初始设计变量 $x_i\left(0\right)=1$ ，收敛误差取 $e={10}^{-3}$ ，采用反向控制规则。参数定义完成后用户点击“查看动画”按钮即可看到每一代拓扑优化构形图，如图4(a)所示为最终拓扑优化构形图。用户点击“应力云图”可查看图5所示的最终拓扑优化构形的应力云图。从图5中可见最大应力出现在集中载荷F加载处，应力值最低的蓝色区域与拓扑优化结果中得到的需去除区域保持一致，同时图4(a)与图4(b)中文献 [9] 的HCA 算法求得的拓扑构形相似，从而验证了本文软件进行结构拓扑优化的有效性和可靠性。

Figure 3. TOCS standard test component base structure

图3. TOCS标准测试件基结构

(a) HCAM的拓扑优化构形 (b) 文献 [9] 的拓扑优化构形

Figure 4. Results comparison of different optimization algorithms

图4. 不同优化算法的结果对比

Figure 5. Stress nephogram of final topological optimization configuration

图5. 最终拓扑优化构形的应力云图

5. 结论

本文在HCAM理论研究基础上，开发了HCAM结构拓扑优化设计软件平台，并通过经典案例验证了软件功能的有效性和可靠性，有效地减少了工程技术人员的使用难度及工作量，为连续体结构拓扑优化设计生产提供了有益的参考。

基金项目

本文受到国家自然科学基金项目(51605253)、国家级大学生创新创业训练计划项目(201811488004)、浙江省博士后择优资助项目(zj20180077)和浙江省基础公益研究计划项目(LGG18E050014)资助。

NOTES

^*通讯作者。