基于SIMP插值模型的棋盘格现象滤波控制算法研究

张 宇，罗陆锋

（天津职业技术师范大学 机械工程学院，天津 300222）

0 引言

棋盘格是连续体结构拓扑优化中较为常见的一种数值不稳定现象[1,2]。所谓棋盘格，是指在结构优化过程中某些区域出现的单元材质密度周期性分布的一种现象。棋盘格的出现与分析单元的选择有关，与单元的类型、单元的划分、单元弹性模量与单元密度之间的关系以及优化算法本身等都有关系。除了棋盘格以外，在拓扑优化的数值计算中一般还存在着网格依赖性、局部极值、多孔材料等问题。棋盘格和网格依赖两种现象一般同时出现在优化结果中，能够有效去除棋盘格的方法通常也能有效克服网格依赖现象[3]。

拓扑优化中数值计算不稳定性的消除非常重要，它关系到数值计算的收敛性和计算结果的可制造性问题。解决棋盘格除了采用高阶有限单元代替低阶单元方法外，一般还有周长约束、滤波法等。从算法本身来说，周长约束等于在原有的模型中添加了约束类型，可能造成收敛性和数值计算问题。滤波法是图像处理技术中的常用方法，其基本思想是将原函数与滤波函数进行卷积运算[3]，对原函数进行规整化处理，提高函数光滑性。该思想应用到拓扑优化中的具体表现，就是通过调整每次迭代后单元的灵敏度从而避免棋盘格的出现。在这种方法中，某一单元的灵敏度依赖于其自身以及滤波半径范围内相邻单元灵敏度的加权平均值。其优点是不需要在优化问题中加入额外约束，另

外容易实施，收敛性好，计算稳定。这种方法对消除棋盘格问题非常有效，也能一并解决网格依赖性问题。

1 基于SIMP插值模型的滤波算法

1.1 灵敏度分析

SIMP[1,4,5]方法的思想和前提是将离散单元内部的材料属性定义为常数，设计变量定义为离散单元的相对密度，尽量减少结构中间密度单元的数目，使结构单元密度尽可能为0或1[6]，从而用连续优化设计方法来近似离散优化设计[7]。

根据变密度方法的思想，单元材料为各向同性材料，假设单元内材料的宏观性质为各向同性，其泊松比与实体材料的泊松比相等。单元刚度矩阵与单元材料的等效性质相关[8]。为了控制整个结构位移，设目标函数为：

其中，C( ρ)为结构柔度，F和U分别为节点载荷阵和位移阵，K为刚度矩阵。如果我们以最大位移为约束条件，根据变密度方法的思想，当对连续体结构进行有限元离散以后，假设位移最大的节点的位移受到了限制，则约束条件如下：

式中u*为根据有限元分析的结果所给定的合理的位移上限，ρi代表单元密度，ε为人为给定的密度下限，i=j1, j2,…, jk指优化后密度保持不变的单元号。

有限元平衡方程为：

在插值前和插值后的材料弹性模量之间引入关系式：

其中，EP代表插值以后的弹性模量，ρiP代表迭代后的单元密度，E0和Emin分别代表固体和去除部分的弹性模量，弹性模量E( ρ)与实体弹性模量E0的比值同相对密度近似成指数关系，对上式进行简化：则目标函数可以修改为：对材料的平衡方程进行微分：该方程与目标函数相结合，可得：

根据以上各式，容易得到结构总体柔度的灵敏度等效形式：

1.2 滤波算法实现

滤波法本来是图像处理技术中的常用方法，其数学表达式为：

式中，F代表原函数，G代表滤波函数，该表达式的基本思想是将原函数与滤波函数进行卷积，对原函数进行规整化处理，提高函数光滑性。

将该思想引用到连续体拓扑优化中，针对本文变密度法模型，利用有限元法对设计区域进行离散化，每次迭代计算后对灵敏度场进行滤波，以滤波后的灵敏度场作为下次迭代的初始值，直至收敛。设j单元在某一次迭代后的灵敏度数值为，

现计算j单元滤波处理后的灵敏度值。在滤波范围内每个单元具有一个权重因子hi，hi的大小与该单元和当前计算单元的中心距离成反比，计算式如下：

式中，rf表示滤波特征半径，disk(i, j )表示单元i和当前计算单元j的距离。当disk(i, j )＞rf时，hi取值恒为0。

则单元滤波后的灵敏度值可以表示为：

在本文所采用的拓扑优化算法中，每一次迭代后，有一部分单元会被删除，也就是说这一部分已被删除的单元不应该再在滤波过程中起作用，为了忽略这部分单元的影响，将模型中加入惩罚系数p，则滤波函数可以修改为：

当单元存在时，pi取值为1，否则为零。利用这种方法在每一次迭代后对设计区域中所有设计单元进行滤波，形成滤波后的灵敏度场，作为下一次迭代的初始值，如此循环直到得到最优拓扑形式。

从该算法的滤波半径设置以及滤波公式可以看出，在滤波半径范围之内的单元对单元灵敏度更新有影响，并且当滤波半径趋近于0时，修改后的灵敏度值趋于不变，当滤波半径趋近于无限大时，各单元灵敏度相等。通过对每次迭代后全局灵敏度进行修改，可以有效克服数值计算中的棋盘格现象和网格依赖。

2 数值模拟分析

添加了滤波算法后的基于SIMP的变秘密度法拓扑优化方法的算法流程图如图1所示。

采用数值模拟方法加以验证。分析模型采用单边（左侧）固支悬臂梁，右边中点加载集中载荷。梁的弹性模量为210×109，泊松比为0.3，厚度为2mm，右侧中点加载向下集中载荷为1000N，位移约束为-0.18×10-3m。有限元划分网格24×16个。

图1 加入滤波算法后的变密度拓扑优化流程

图2和图3分别给出了未加入滤波算法的结构拓扑优化结果和加入了滤波算法之后的新的拓扑优化结果。表1是加入滤波算法前后的各项数据指标。

图2 未加入滤波算法前的拓扑结构

图3 加入滤波算法后的拓扑结构

从图2到图3的结构拓扑优化结果可以看出，在加入滤波算法后，棋盘格现象消除；由表1可以看出，滤波后所得结构满足约束要求，结构最终质量略有变动但变动不大。棋盘格结构中存在的应力集中现象也得到了一定的控制，整体应力分布更加均匀。

表1 加入滤波算法前后的各项指标

3 结束语

基于工程中常用的SIMP变密度函数模型，进行了灵敏度分析，提出了以结构最小柔度为目标函数、位移为约束的变密度插值模型为基础的滤波思想，并推导了具体的计算方法。通过算例分析，明显消除了棋盘格现象，同时对网格依赖性有所控制。证明了这种滤波方法是切实有效的。

在滤波算法算例中，设置的滤波半径为1.8，也就是说采取的是二阶滤波，对于滤波算法来讲，选取合适的滤波半径也是非常重要的。滤波半径过小，则灵敏度变化不大，不能起到很好的消除棋盘格式的效果；反之，则会导致灵敏度变化过于剧烈，所得结构完全脱离原来拓扑模型。同时，滤波法难以对全局数值变量进行控制，而且没有严格的理论依据，有时会导致边界扩散化，产生模糊的边界。这也是滤波方法的不足之一，是在具体实施过程中需要密切注意的问题。

[1] 张丽, 饶华球, 昌俊康. 对连续体结构拓扑优化的一点认识[J]. 制造业自动化. 2010, 32(2): 93-96.

[2] 左孔天. 连续体结构拓扑优化理论与应用研究[D]. 武汉:华中科技大学, 2004.

[3] 郭中泽, 张卫红, 陈裕泽. 结构拓扑优化设计综述[J], 机械设计, 2007, 24(8): 1-4.

[4] Mlejnek H. P, Schirrmacher R..An engineer's approach to optimal material distribution and shape finding[J]. Comp.Meth. Appl. Mech. Engrg, 1993, 106: 1-26.

[5] M le jned H p, S chinm ascher R. An engineers approach to optimal material distribution and shape finding[J]. Comput Method Appl Mech Engrg 1993, 106(1-2): 1-26.

[6] 汤颖颖. 基于变密度法的连续体拓扑优化设计[D], 西安:长安大学, 2008.

[7] 邹春江, 左孔天, 向宇. 基于SIMP方法微电容加速度计结构固有频率拓扑优化[J], 科学技术与工程,2011,11(29): 7088-7091.

[8] 牟淑志, 杜春江. 基于ANSYS二次开发的结构拓扑优化[J], 计算机应用与软件, 2010, 27(2): 228-230.