一种考虑结构稳定性的单材料结构拓扑优化方法和系统

技术领域

本发明属于结构拓扑优化领域，更具体地，涉及一种考虑结构稳定性的单材料结构拓扑优化方法和系统。

背景技术

在结构设计过程中，不仅要考虑结构的强度和刚度要求，同时也需要保证在正常作状态下结构的稳定性。在实际工程中，结构通常是由梁、杆、板壳等基本构件组合而成。如果这些构件受压，那么它们的稳定性将成为影响结构安全的一个重要因素。事实历史上也曾经发生过多起由于基本构件的失稳而引起的灾难性事故。因此，在结构概念设计阶段就考虑结构的稳定性是非常有必要的。而结构优化技术可以有效地提高设计质量，降低产品的制造和使用成本，因此在结构优化中考虑结构稳定性的结构拓扑优化设计更具有实际意义。

尽管连续体结构拓扑优化中对结构稳定性的研究已取得了一定的进展，但还存在一些问题。首先，在实际工程中结构的刚度和稳定性是设计人员必须考虑的两个重要因素，但已有的拓扑优化提法一般没有同时考虑结构刚度和稳定性；其次，目前所做研究大多基于SIMP(Solid Isotropic Microstructures with Penalization，带惩罚指数的固体各向同性微结构模型)方法，该方法存在以下缺点：优化出来的拓扑结构边界不够清晰，特别是当过滤半径比较大的时候。这些灰度区域没有物理意义，设计如果无后处理无法直接用于制造。Level-set(水平集)法的优点：用一个高纬度的水平集implicitly表达拓扑结构的边界，从而解决了SIMP法的灰度区域问题。拓扑结果边界清晰，无灰度区域，设计可以直接用于制造。缺点：由于设计变量间接与优化问题挂钩，中间涉及一些被水平集切割的有限单元的近似，从而影响优化精度。水平集方程需要用PDE方程来更新，中间还需要重置水平集方程来保证PDE的持续更新，从而大大降低优化收敛速度或者甚至无法收敛。PDE需要连续形状灵敏度来更新，相比SIMP的离散设计灵敏度更难。线弹性体的连续形状灵敏度已经发展很成熟，但是非线性结构连续形状灵敏度非常难求，需要很高的数学基础。现有技术中众多基于参数化水平集方法的结构拓扑优化方法均未考虑结构稳定性。

发明内容

针对现有技术的缺陷和改进需求，本发明提供了一种考虑结构稳定性的单材料结构拓扑优化方法和系统，其目的在于在考虑结构稳定性的同时，使得优化结果更加清晰，分支更加明显。

为实现上述目的，按照本发明的第一方面，提供了一种考虑结构稳定性的单材料结构拓扑优化方法，该方法包括以下步骤：

S1.构建用于描述设计域内单材料结构初始拓扑的初始水平集函数，确定插值初始水平集函数的扩展系数初值；

S2.构建以各扩展系数为设计变量，以单材料结构总柔度最小化为优化目标，以单材料结构的体积和屈曲载荷因子为优化约束条件的拓扑优化模型；

S3.在当前扩展系数和水平集函数下，基于优化拓扑模型计算单材料结构的结构总柔度、体积和屈曲载荷因子；

S4.根据单材料结构的结构总柔度、体积、屈曲载荷因子和拓扑优化模型，更新扩展系数和水平集函数；

S5.判断当前扩展系数与参考扩展系数的差值是否小于阈值，若是，则优化结束，根据当前水平集函数确定设计域内单材料结构拓扑并输出；否则，转入步骤S3。

参考扩展系数可以是上一轮的扩展系数，也可以是综合前几轮的扩展系数得到，也可以人为设定。

优选地，屈曲载荷因子用于单材料结构表征稳定性，其计算过程如下：

(1)根据当前水平集函数确定设计域内单材料当前结构拓扑；

(2)计算单材料当前结构的整体几何刚度矩阵和整体刚度矩阵；

(3)由整体几何刚度矩阵和整体刚度矩阵构建特征值方程；

(4)求解特征方程，特征值即为屈曲载荷因子。

优选地，使用插值矩阵方法计算整体几何刚度矩阵，过程如下：

(1)计算单材料应力σ＝DBU；

(2)基于单材料应力σ构建应力矩阵S，构建单元形函数相关的插值矩阵G；

(3)计算单元几何刚度矩阵k

k

(4)根据节点耦合方式，将单元几何刚度矩阵k

其中，D表示平面应力问题下的弹性矩阵，B表示应变矩阵，U表示单元节点位移场，σ

有益效果：本发明使用插值矩阵方法计算几何刚度矩阵，使得求解出的几何刚度矩阵更加准确，对结构的稳定性计算更加准确，使最终的结构拓扑产生更多分支。

优选地，所述拓扑优化模型的表达式如下：

Find：α＝[α

其中，α

有益效果：第一个约束的意义：结构拓扑优化的本质是为了在保证结构力学性能的前提下减少结构的质量，通过设置该约束可以控制拓扑优化后结构的质量，起到减重的目的，如果没有磁约束，那么将不能使结构的拓扑优化起到减重目的。

第二个约束的意义：屈曲载荷因子是衡量结构稳定性的指标，本发明通过限制最小的屈曲载荷因子大于指定值，以此使结构的稳定性满足工程要求，当该约束未得到满足时，则优化出的结构的稳定性达不到要求。

第三个约束的意义：限制扩展系数，使基于参数化水平集方法的结构拓扑优化更加稳定，否则容易造成优化不收敛等问题。

第四个约束的意义：使用径向基函数对水平集函数进行插值，通过插值可以有效解决传统水平集函数在优化过程中速度场需要扩展到整个设计域和无法与许多成熟算法结合等数值求解问题。

优选地，步骤S4包括以下步骤：

(1)分别对拓扑优化模型中的目标函数、体积约束条件、屈曲载荷因子约束条件对扩展系数求偏导，作为结构总柔度对扩展系数的敏感度、体积对对扩展系数的敏感度、屈曲载荷因子对扩展系数的敏感度；

(2)将上述敏感度、目标函数值、体积值和屈曲载荷因子数值同时代入移动渐近线方法，求解得到新扩展系数，进一步确定新水平集函数。

有益效果：本发明采用MMA对整个结构进行优化分析，该方法求解多约束问题，尤其是屈曲约束问题，能够保证优化的准确性和收敛性。

优选地，计算屈曲载荷因子约束对设计变量的敏感度，过程如下：

(1)构建关于屈曲载荷因子的平衡方程弱形式公式：

a(u,v)＝λ·b(u,v)

其中：

a(u,v)＝∫

b(u，v)＝∫

其中，a(u，v)表示能量双线性函数，l(v)表示载荷线性函数，v表示在动力学上允许的位移场中的一个虚拟位移域，ε为单元应变场，D为结构整体的弹性矩阵，σ为单元应力场；

(2)根据形状导数的定义及其引理，推导能量双线性形式和载荷线性形式关于时间变量t的导数

其中，v

(3)将平衡方程弱形式两边对时间t求偏导；

(4)提取所有含

(5)将步骤(1)-(3)得到的式子代入步骤(4)的式子，得到

(6)使用伴随变量w替换步骤(5)的式子中的v，得到

伴随变量w通过求解下列双线性形式的伴随方程获得：

a(u，w)＝λ·b(u，w)

∫

(7)求得屈曲载荷因子关于t的导数

(8)通过链式法则对屈曲约束直接求其关于时间变量t的偏导数得到：

进一步得到

有益效果：本发明对屈曲载荷因子对设计变量的敏感度求解方式，使求解出的扩展系数更加准确，使最终的结构拓扑产生更多分支。

优先地，水平集函数采用高斯径向基函数进行插值。

有益效果：本发明通过使用高斯径向基函数插值，可以保证屈曲载荷因子的计算更加准确。

优选地，使用DWT压缩优化过程中整体几何刚度矩阵和整体刚度矩阵。

有益效果：本发明采用DWT减轻全插值矩阵造成的计算成本。

为实现上述目的，按照本发明的第二方面，提供了一种考虑结构稳定性的单材料结构拓扑优化系统，包括：计算机可读存储介质和处理器；

所述计算机可读存储介质用于存储可执行指令；

所述处理器用于读取所述计算机可读存储介质中存储的可执行指令，执行第一方面所述的考虑结构稳定性的单材料结构拓扑优化方法。

总体而言，通过本发明所构思的以上技术方案，能够取得以下有益效果：

相对于现有技术中基于水平集函数的结构拓扑优化，本发明构建优化模型时引入基于扩展系数的屈曲载荷因子约束，通过扩展系数插值水平集函数，水平集函数决定了伪密度，伪密度决定了几何刚度矩阵和刚度矩阵，而屈曲载荷因子是用它们两个构成的特征方程求解出来的，因此，扩展系数将影响最后的屈曲载荷因子数值，从而在优化过程中将单材料结构的稳定性纳入考虑，使得优化结果更加清晰，分支更加明显。

附图说明

图1是本发明提供的一种考虑结构稳定性的单材料结构拓扑优化方法的流程示意图；

图2是本发明实例提供的悬臂梁的示意图，其中，(a)为初始设计域，(b)为初始化孔洞位置；

图3是本发明实例提供的悬臂梁优化结果图；

图4是本发明提供的目标函数(柔度)和体积分数迭代图；

图5是本发明提供的屈曲约束与屈曲载荷因子关系图；

图6是本发明提供的是屈曲载荷因子-柔度图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

首先，对本发明涉及的术语进行解释：

单材料，是指只具有一种材料属性的材料，如泊松比、弹性模量等。

单材料结构，是指只由一种材料构成的结构，适用于工程中所有单材料结构，例如航空航天所使用的结构。

如图1所示，本发明提供一种考虑结构稳定性的单材料结构拓扑优化方法，包括以下步骤：

(1)基于水平集方法的单材料结构稳定性拓扑优化设计，将其与矩阵压缩技术相结合，以此提高整体稳定性分析效率；

具体地，通过隐式水平集函数构建不同材料相的结构，以插值隐式水平集函数的扩展系数为结构设计变量，在N个固定水平集节点处的结构隐式水平集函数为：

其中，x＝x

其中，c是形状参数，等于水平集网格面积或者体积的倒数；x

进一步地，为提高优化效率，DWT是减轻全插值矩阵造成的计算成本的关键。将原始插值矩阵A转换为其相同大小的小波形式

假设预定义的t＝0时的水平集函数为Φ

如上所述，将矩阵压缩技术引入到基于高斯径向基函数的参数化水平集优化设计中，形成了一种新PLSM，不仅计算成本低，而且其优化设计的性能显著提高。在优化过程中可以发现仅仅计算了一次

(2)基于应力插值矩阵及单元形函数插值矩阵求解单元几何刚度矩阵，进一步通过构造线性特征值方程，计算结构屈曲载荷因子。

具体地，结构的屈曲载荷因子通过求解线性特征值公式可得到，其表达形式为：

(K+λ

其中，K和K

F

其中，F为施加在结构上的外载荷，F

由线性特征值公式可知，若想求得结构的屈曲载荷因子，关键是计算结构的整体刚度矩阵K和整体几何刚度矩阵K

单元刚度矩阵k

k

其中，B为应变矩阵，D为平面应力问题下的弹性矩阵。

单元几何刚度矩阵k

k

其中，S和G分别为应力矩阵和与单元形函数相关的插值矩阵，计算公式如下：

G＝T.[M

其中，U和ε分别表示单元节点位移和应变，N表示为四边形单元形函数。

(3)对于水平集函数的Heaviside函数H(φ

(4)基于参数化水平集的单材料结构理论建立考虑结构稳定性的单材料结构最小柔度拓扑优化模型，并在结构设计域中通过有限元分析求解整体结构的位移场，根据得到的位移场计算单材料结构最小柔度拓扑优化模型的目标函数；接着，基于自伴随方法与伴随变量法对结构的设计变量进行灵敏度分析，并采用MMA移动渐近线算法更新全局设计变量，继而确定满足稳定性约束的结构中材料的最优分布。

具体地，所述的基于参数化水平集理论建立的考虑结构稳定性的单材料结构最小柔度拓扑优化模型的表达式为：

Find：α＝[α

Subject to：G(Φ)＝∫

a(u，v)＝l(v)

α

其中，N为有限元网格节点数量，J(u，Φ)和G(Φ)分别表示结构的总柔度及材料的体积，λ

在目标函数公式中，f(u，u)代表结构的应变能：

其中，ε表示为结构应变场，D为材料的弹性矩阵。

基于虚功原理，针对有限元平衡方程弱形式进行计算，对应的弱形式如下：

其中，a表示双线性能量式；l表示单线性荷载形式；dΩ为结构设计域的积分算子；H表示Heaviside函数，用于表征结构形式的特征函数；ε为应变场；T表示矩阵的转置；u表示结构场的位移；ν表示在动力学上允许的位移空间U中的一个虚拟位移；τ表示应用在边界

构建考虑结构稳定性的单材料结构最小柔度拓扑优化模型具体包括以下步骤：

(3.1)第一步为整个优化中的初始化过程，首先需要定义优化过程中约束，即结构中材料体积上限V

(3.2)第二步是限元分析，在有限元分析中，首先需要定义材料物理参数，其次分别计算节点位移U、整体刚度矩阵K和整体几何刚度矩阵Kg，并代入求解线性特征值方程(K+λKg)Φ＝0，计算屈曲载荷因子λ；

(3.3)第三步根据上一步计算得到的位移U计算结构整体的应变能，即目标函数J(u，Φ)；

(3.4)第四步中完成对目标函数敏度

(3.5)第五步对上一步得到的敏度进行敏度过滤及对敏度矩阵小波变换，采用适用于多约束优化问题的优化算法——移动渐近线方法(MMA)对整个结构进行优化分析，并根据MMA优化算法得到的扩展系数更新水平集函数，实现结构的演化；最后通过计算本次与前一步的目标函数差值是否大于1x10-6来决定是否进行下一次优化求解。

具体地，根据链式求导法则计算目标函数与约束函数针对优化设计变量的一阶微分，计算如下：

目标函数与约束函数针对时间t的一阶微分：

能量双线性形式a(u，v，Φ)及载荷线性形式l(v，Φ)关于时间变量t的导数分别为：

通过链式法则直接对目标函数J(u，Φ)求其关于时间变量t的偏导数：

得到目标函数J(u，Φ)关于扩展系数α的导数：

同样原理，可以得到各相材料体积约束关于设计变量α的敏度：

结构的屈曲优化问题是一个典型的非自伴随问题，因此，通过使用伴随变量方法避免直接对u的求导，使用伴随变量w替换公式中的v，可得

将水平集函数速度场v

通过以上求解即可获得

更新设计变量并判断优化模型是否满足收敛条件，若否，则返回到步骤(3.2)；若是，则输出单材料结构的最优拓扑结构。由于以上分析是多约束问题，因此本发明采用MMA移动渐近线方法进行优化求解。

实施例

本实施例给出一个短悬臂梁的稳定性优化算例来说明所提出的研究方法的有效性，图2是本发明实例提供的悬臂梁的示意图，其中，(a)为初始设计域，(b)为初始化孔洞位置。如图2中(a)所示，将悬臂梁的左端固定，在梁的右端中点处施加如下图所示的垂直向下的集中载荷F＝1，并通过宽度和高度均为单位大小的四节点等参单元离散成N＝40*80的有限元网格(本研究仅考虑平面应力问题)。用均匀分布的圆孔对短悬臂梁结构进行初始化，结果如图2中(b)所示。在优化过程中材料的弹性模量设置为E＝1，空白材料的弹性模量设置为E0＝10-6，泊松比均为0.3。无量纲材料特性的设置是为了方便不同设计之间的比较，以下算例中均采用这种策略。

从图3中可清晰看出，随着屈曲载荷因子的下限逐渐增大，上半部分受拉构件逐渐向下偏移，下部受压构件开始逐渐变粗并开始出现分支结构，以此来提高结构的稳定性，从优化结果中还可以看出，基于PLSM的考虑结构稳定性的拓扑优化设计结果边界更加清晰，且分支结构可以明显显示，并不会像基于SIMP方法是以大部分中间密度来支撑，以此增加结构稳定性，这也是基于PLSM的考虑结构稳定性的拓扑优化设计方法的优点所在。从图3中还可以发现，随着屈曲约束的增加，结构边界的边界会出现微小的波动，这在基于SIMP的结构稳定性分析中也可以看到，这是由于屈曲约束的施加后导致数值波动造成的，这种波动会随着屈曲约束的增加而变得越来越明显，在今后的研究中可通过改变优化算法来进行改进。

图4给出了体积约束为0.2时的目标函数和结构体积的迭代曲线图，从迭代曲线中可以看出，在迭代初期，迭代曲线存在一定的波动，这是由于材料用量少时，基于PLSM方法的结构拓扑在演化初期存在结构突变导致的，但这并不影响后续的优化进程，从图中可以看出，在后期的迭代中，目标函数和体积分数的迭代曲线是非常平稳的。

如图5所示，表示当体积约束为0.2时第一阶屈曲载荷因子和屈曲约束的关系曲线，从图5中可以看出，当到达某一极限值时，随着屈曲约束的增加，屈曲载荷因子并非一直都能满足约束条件，而是有一个极限值，但从曲线图中可以观察到，尽管屈曲约束没有得到满足，但是相比极限值之前的结构屈曲载荷因子，在极限值之后的屈曲载荷因子也会有可能变大。

如图6所示，表示当体积约束为0.2时第一阶屈曲载荷因子和柔度的关系曲线，从图6中得出以下结论：在材料体积给定的情况下，在极限值之前，随着屈曲载荷因子的提高，拓扑结构的柔度会增大，这说明结构的稳定性的提高是以牺牲结构的刚度为代价的。从图6中的第一阶屈曲载荷因子和柔度的数值可以发现，当材料的体积多时，其刚度和稳定性均要高于低材料体积。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。