郑伟伟，周克民

（华侨大学 土木工程学院，福建 厦门361021）

结构拓扑优化的研究包括均匀化方法［1］、演化优化算法［2］、水平集方法［3］、ICM 法［4］等.为了获得清晰的结构，抑制中间密度单元，这会引起了数值不稳定现象，如棋盘格现象［5］，需要格外的技术解决这个问题［6-8］.Michell［9］揭示了拓扑优化结构是非均匀质各向异性类桁架连续体，即由无限密的杆件构成，这一特性启发了从类桁架连续体中获取最优拓扑的想法［10-11］.工程中更多地使用带孔等厚板，通过剔除密度小于一定值的节点形成孔洞.演化优化算法［12-13］有采用过这种做法，但是由于演化优化算法采用的是各向同性材料，且设计变量为单元的密度，因此，该算法得到的结果会出现棋盘格现象且与网格密度有关.本文通过优化孔洞的分布形成带孔等厚板的优化算法，在优化迭代过程中逐步由类桁架连续体演化为带孔板，使优化结果更接近工程需要.

1 类桁架连续体材料模型

1.1 弹性矩阵

在类桁架连续体材料模型中，材料模型是由无限多的非均匀连续分布的杆件构成，假设任意点都是由两组正交的杆件构成，它们在任意点的密度用t1和t2表示.假设应力σi和应变εi的线性关系为

式（1）中：E是弹性模量.

式（2）中：diag［·］表示对角矩阵.

在这种假设下，当t1＝t2时，式（2）的各项为同性材料.若假设结构坐标轴到材料坐标轴的转角为α，那么结构坐标系中的弹性矩阵可由坐标转换矩阵获得，即

式（3）中：T（α）为坐标转换矩阵，即

弹性矩阵（3）可以表示为

式（5）中：sb，r和gr（α）分别为常数矩阵和函数矩阵分量，Ar为常数矩阵，有

节点j（j＝1，2，…，J）处的杆件密度t1，j，t2，j和方向αj作为设计变量，J是节点的总数，设计变量节点j处的弹性矩阵可以表示为

有限元内部任意一点的弹性矩阵可以由该单元节点位置的弹性矩阵插值得到

式（7）中：Nj（ξ，η）是型函数；ξ，η是局部坐标；Se是属于单元e的节点的集合.

将式（6）带入式（7），得到

1.2 刚度矩阵

将单元刚度矩阵（8）带入刚度矩阵，有

2 优化方法

将节点处的杆件密度和方向作为设计变量，材料的体积作为目标函数，每一节点施加应力约束，优化问题的列式为

式（12）中：σb，j为节点处杆件b的应力；σp为允许应力.

迭代过程有如下5个步骤.

步骤1 设计域被划分为有限个单元，单元采用4节点矩形单元，杆件初始密度和方向设置为

上标0代表迭代次数，初始密度剔除标准设为0；

步骤2 通过有限元方法，计算节点处的主应力方向和主应力方向上的应变的大小；

步骤3 利用满应力准则，优化杆件的密度，杆件的方向与主应力的方向一致，即

步骤4 为了防止刚度矩阵奇异，设置低于密度剔除标准的杆件密度为t0（这里取10-6），而不是直接设置为0，高于平均密度的杆件的密度将设置为

步骤5 密度剔除标准.随迭代次数而增加，tic＋1＝tic＋Δt，Δt是密度剔除标准的增量.

通过大量的数值计算，取10-4比较合适.太小则孔洞形成的较慢，需要更多步的迭代；太大则形成的拓扑与最优拓扑相差较大.重复步骤2～5，直到没有密度的节点数占总节点数的60%，剔除材料比例可以由用户根据需要自行决定.

图1 悬臂梁尺寸Fig.1 Dimension of cantilever beam

3 数值算例

运用文中的算法，计算两个经典的算例，两个算例中：弹性模量E＝210 GPa；允许应力σp＝160 MPa.

第1个算例是一个右端中点受集中荷载的悬臂梁，尺寸如图1所示.图1中：厚度（h）为0.01 m；设计域被划分为64×40矩形单元.所有节点密度（两个方向的杆件密度之和）的等值线填充图，如图2，3所示.图2，3中：等值线的间距为最大密度的5%.图2是迭代31次后的结果，设计域开始出现一些孔.图3是迭代62次后的结果，孔不断扩展、溶合，结构拓扑趋于简单.

图2 迭代31次后的密度等值线图Fig.2 Isogram of density after 31 iterations

图3 迭代62次后的密度等值线图Fig.3 Isogram of density after 62 iterations

大部分区域的密度分布较均匀，只有支座和荷载作用点附近区域密度较大，但没有出现严重的应力集中现象.体积比的定义为优化后结构总体积占初始体积的百分比.体积比（Rv）和保留的节点数百分比（Rn）随着迭代次数（N）的变化曲线，如图4所示.

从图4可知：由于荷载比较小，优化后纤维密度较初始密度小很多，优化后的体积只占初始体积很小一部分.保留的节点百分比指密度大于临界值（这里取2t0）的节点占总节点数的百分比，反映了类桁架连续体演化为带孔板的进程，当其等于40%时，迭代停止.满应力优化准则可以在少量的迭代下使结构体积达到最小，通过密度限值约束，去除小密度杆件，结构离散得越来越简单，但体积只增加了一点.

第2个算例是跨中受集中荷载的简支梁，尺寸如图5所示.图5中：厚度（h）为0.01 m；设计域被矩形单元划分为60×40.

图4 体积比及节点数比的变化Fig.4 Variation of volume ratio and nodes ratio

未经密度限值约束迭代100次得到的优化结果，如图6所示.图6中：为了使荷载集中点的密度不至于大其他区域的密度太多倍，仍施加了密度上限为平均密度的约束.从图6可知：未经密度限值约束的优化结果形成连续材料场，而不能形成杆件.

增加密度限值约束后迭代56次的结果，如图7所示.对比图6，7可以看出：增加密度限值约束后，在设计域的中间部分形成了几根杆件，而且与圆拱的密度基本一致，比较均匀.

图5 简支梁尺寸Fig.5 Dimension of simply supported beam

图6 没有密度限值约束的优化结果Fig.6 Optimal result without density restriction

图7 经密度限值约束的优化结果Fig.7 Optimal result with density restriction

体积比和保留节点数百分比这两项指标随迭代次数的变化，及施加密度限值约束前后的对比，如图8所示.从图8可知：施加密度限值约束后的体积比会比施加前稍微大一些，且体积比会随着迭代次数逐步提高；对比施加密度限值约束前后节点数比，非零密度（大于2t0）节点比在施加密度限值约束后比施加前减少快很多.随着迭代次数的变化线形，这一指标决定于密度剔除标准的确定方法，因此，可以通过改变密度剔除标准的函数形式（线性或非线性）来控制类桁架连续体向带孔板演化的过程.

图8 施加密度限值约束前后指标变化及对比Fig.8 Indexes variation and contrast before and after density restriction

4 结束语

研究了基于类桁架连续体材料模型的带孔板拓扑优化方法，通过在每一迭代步中剔除小于一定值的密度形成孔洞，限定大于一定值的密度达到密度均匀的目的，最终得到带孔板.下一步工作是对连续体施加最小宽度的约束，从而可以得到孔数量不同的连续体拓扑优化结果.

［1］ BENDSØE M P，KIKUCHI N.Generating optimal topologies in structural design using a homogenization method［J］.Comput Methods Appl Mech Eng，1988，71（2）：197-224.

［2］ XIE Yin-min，STEVEN G P.A simple evolutionary procedure for structural optimization［J］.Comput Struct，1993，49（5）：885-896.

［3］ SETHIAN J A，WIEGMANN A.Structural boundary design via level set and immersed interface methods［J］.J Comput Phys，1999，163（2）：489-528.

［4］ 隋允康，彭细荣.结构拓扑优化ICM 方法的改善［J］.力学学报，2005，37（2）：190-198.

［5］ SIGMUND O.Numerical instabilities in topology optimization：A survey on procedures dealing with checkerboards，mesh-dependencies and local minima［J］.Struct Optim，1998，16（1）：68-75.

［6］ HABER R B，JOG C S.BENDSφE M P.A new approach to variable-topology design using a constraint on the per-imeter［J］.Structural and Multidisciplinary Optimization，1996，11（1）：1-12.

［7］ PETERSSON J，SIGMUND O.Slope constrained topology optimization［J］.International Journal for Numerical Methods in Engineering，1998，41（8）：1417-1434.

［8］ Sigmund O.Materials with prescribed constitutive parameters：An inverse homogenization problem［J］.International Journal of Solids and Structures，1994，31（17）：2313-2329.

［9］ MICHELL A G M.The limits of economy of material in frame structure［J］.Phil Mag，1904，8（6）：589-597.

［10］ ZHOU Ke-ming，LI Xia.Topology optimization of structures under multiple load cases using a fiber-reinforced composite material model［J］.Comput Mech，2005，38（2）：163-170.

［11］ ZHOU Ke-ming，LI Xia.Topology optimization for minimum compliance under multiple loads based on continuous distribution of members［J］.Struct Multidisc Optim，2008，37（1）：49-56.

［12］ XIE Yin-min，STEVEN G P.Optimal design of multiple load case structures using an evolutionary procedure［J］.Engineering Computations，1994，11（4）：295-302.

［13］ CHU D N，XIE Yin-min，HIRA A，et al.Evolutionary structural optimization for problems with stiffness constraints［J］.Finite Elements in Analysis and Design，1996，21（4）：239-251.