用可退化有限单元进行平面连续体拓扑优化

    1　引　言

    结构优化分为尺寸、形状和拓扑优化三个层次。前两个层次分别考虑杆件的横截面尺寸和结点位置,或板的厚度分布。而拓扑优化明显区别于前两者,难度更大,潜力也更大。结构拓扑优化研究对象分为杆系和连续体结构。杆系结构的拓扑优化考虑结点的多少、位置,以及各结点之间的连接关系,目前有有基结构法[1]和带随机搜索性质的优化方法[2]。基结构法是在有限的子空间内寻优,容易丢失最优解,另外还存在组合爆炸、解的奇异性[3]等问题。带随机搜索性质的方法克服了上述问题,但目前的计算量还过大,并不实用。连续体结构的拓扑优化要考虑连续体中空洞的多少,位置及大小,目前有均匀化方法[4]和比例函数方法[5]。均匀化方法理论严谨,但计算复杂,存在数字计算不稳定等问题[6];比例函数方法运算量小,更实用。本文采用的是比例函数方法。已经证明[7],在一定的条件下,连续体的拓扑优化结果就是桁架,从而使连续体结构的优化方法也可以应用于杆系结构。

    在以连续体为研究对象的结构拓扑优化方法中,一般将设计域划分成规则的矩形单元,通过删除单元实现结构拓扑优化。在形成单元时,无法预知最后形状,而单元一旦形成,就只能删除不能改变,这会给最后结果的边界描述带来困难,形成锯齿状边界。如果每次结构分析前都重新划分单元,其工作量巨大,尤其是三维复杂问题,重划分单元几乎是不可能的。这是以连续体为基结构的拓扑优化计算中所普遍存在的问题。本文根据计算结果将矩形单元退化为三角形单元,在不增加单元及计算自由度的情况下,使边界更为光顺。

    2　计算方法

    2.1　变厚度有限元刚度矩阵

    本文使用了矩形和三角形两种有限单元。矩形单元厚度根据四个节点的厚度,利用形函数构造:

其中:t_i为节点厚度,N_i为形函数。形函数N_i、应变矩阵Bi及弹性矩阵D在一般有限元书上都能查到,这里从略。本文推导出了变厚度单元刚度矩阵如下:

式中:,ξ= x/a,η= y/b,E为弹性模量,μ为横向变形系数,为单元在四个结点处厚度的平均值,亦即单元的平均厚度,其它参数为:

三角形单元的厚度取三个节点处单元厚度的平均值,其它与一般方法相同。

    2·2　变厚度板优化过程

    采用满应力法及HK技术。将初始设计域划分为规则的矩形单元,设各节点处初始厚度t0,由有限元分析得到各单元在节点处Misses应力,将围绕每一节点i的所有单元在节点i处应力的加权平均值作为节点i的应力σi,权重取单元面积。本文仅使用了矩形及三角形两种单元,且每种单元面积相同,三角形单元均为矩形单元的一半,故本文实际上矩形单元权重为2,三角形单元为1。采用下述迭代公式,确定下一轮计算时节点处的厚度: