本发明属于多相材料拓扑优化设计领域，更具体地，涉及一种考虑结构稳定性的多相材料拓扑优化设计方法和系统。
背景技术：
在受压载荷的作用下，机械结构部件的主要失效形式通常有两种:材料的失效或失稳和结构的失稳。对于结构的失稳也就是通常称为的结构屈曲现象，即在没有超过结构的强度限制，结构也会由于屈曲变形，导致失稳破坏。对于结构的失稳问题可主要分为两种情况：线性屈曲分析和非线性屈曲分析，其中非线性屈曲分析由于计算结果较为准确，但由于其计算复杂及求解效率低，一般在结构拓扑优化中不被考虑；然而线性屈曲分析虽然计算结果的准确性不如非线性屈曲分析的好，但由于其计算简单及求解效率高而被广泛应用于结构拓扑优化设计中。在结构拓扑优化设计中，线性屈曲分析主要起到预测作用，它可以保证结构可以在给定的载荷下不发生结构屈曲破坏。对于连续体结构的拓扑优化，当结构的限制体积很小时，优化结果易产生细长的杆状结构，如果该结构受压，那么结构的失稳就会成为其主要的失效形式。因此，对于这种情况，如果不考虑结构的屈曲失稳就很可能导致不理想的设计结果。
与传统的单材料结构相比，多相材料结构可以充分发挥每种材料的性能优势，如耐腐蚀、耐磨、重量轻等，尤其在工程应用中，如飞机等，已经成为拓扑优化领域的研究热点。随着增材制造技术的发展，对于包含多种不同材料的零件结构的3d打印也已经成为了一种新的趋势，因此研究多相材料结构的稳定性优化问题具有重要意义。
在拓扑优化领域中，多材料拓扑优化是在给定的设计域内寻找多种不同材料的最优分布的设计。现有对于多相材料结构稳定性设计问题研究，由于其设计变量多，优化约束多，导致多相材料结构的屈曲分析变得更加复杂，且求解效率低。而且目前已有的多材料稳定性研究是基于密度方法进行研究，使用该方法拓扑优化后的结构的边界会存在大量灰度单元，导致优化后的多材料结构边界模糊。
因此找到一种边界光滑、优化效率高的多相材料结构设计方法与稳定性问题相结合是当下亟待解决的问题。
技术实现要素：
针对现有技术的以上缺陷或改进需求，本发明提供了一种考虑结构稳定性的多相材料拓扑优化设计方法和系统，其目的在于实现边界光滑、优化效率高且稳定性提高的多相材料拓扑。
为实现上述目的，本发明提供了一种考虑结构稳定性的多相材料拓扑优化设计方法，包括：
s1.以扩展系数作为结构设计变量，通过插值隐式水平集函数构建不同材料相的结构，求取各相材料的体积；
s2.构建以各扩展系数为设计变量，以多相材料结构总柔度最小化为优化目标，以各相材料的体积和失稳载荷因子为优化约束条件的拓扑优化模型；
s3.在当前扩展系数和水平集函数下，计算多相材料整体结构失稳载荷因子；
s4.计算目标函数与优化约束条件对设计变量的敏感度，根据计算结果更新扩展系数和拓扑优化模型；
s5.判断拓扑优化模型是否满足收敛条件，若是，则输出多相材料结构的最优拓扑结构；若否，则返回到步骤s3。
进一步地，所述拓扑优化模型的表达式为：
find:αk＝[α1,α2,…,αm]t
subjectto:gi(φ)＝∫ωχi(φ)dω≤vi,max
αk,min≤αk≤αk,max
其中，m、αk和n分别为水平集的数量、扩展系数和材料相数量，j(u,φ)和gi(φ)分别表示结构的总柔度，及材料相i的体积，λ*和vi,max分别表示屈曲约束及材料相i的体积约束，λp表示在屈曲模态集j中第p阶失稳载荷因子，αk,max和αk,min分别表示设计变量的上下限，fi(u,u)代表插值前材料相i的应变能，χi(φ)表示第i种材料的特征函数。
进一步地，失稳载荷因子的计算公式为：
(k λkkg)φk＝0
其中，k和kg分别为单材料结构在外力荷载作用下由位移场计算得到的整体刚度矩阵和应力场计算得到的整体几何刚度矩阵，λk为单材料结构的第k阶失稳载荷因子，φk为对应的第k阶屈曲模态，即在第k阶极限载荷fk作用下单材料结构的位移。
进一步地，整体几何刚度矩阵的计算过程为：
01.计算多材料单元应力σ＝χi(φ)dbu；
02.基于多材料单元应力σ构建应力矩阵s，构建单元形函数相关的插值矩阵g；
03.计算单元几何刚度矩阵kg：
kg＝∫χi(φ)gtsgdv
04.根据节点耦合方式，将单元几何刚度矩阵kg组装成整体几何刚度矩阵kg；
其中，d表示平面应力问题下的弹性矩阵，b表示应变矩阵，u表示单元节点位移场，σx、σy、τxy分别表示单材料x方向正应力、y方向正应力、xy平面切应力，j表示雅克比矩阵，mi表示中间变量，pi表示单元形函数，ξ、η表示单元坐标系。
进一步地，基于自伴随方法和伴随变量法计算目标函数与优化约束条件对设计变量的敏感度。
进一步地，失稳载荷因子对设计变量的灵敏度为：
λ为失稳载荷因子，ε为应变场；t表示矩阵的转置，u表示结构场的位移，d是弹性矩阵，σ为结构的单元应力，w为伴随变量，δ表示dirac函数，为heaviside函数的一阶微分，
进一步地，采用mma移动渐近线算法更新扩展系数。
进一步地，引入全局径向基函数gsrbf对水平集函数进行插值。
总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果。
(1)本发明针对多相材料，将参数化水平集方法与稳定性问题相结合，可以得到边界光滑、优化效率高且具有较高稳定性的多相材料分布结构，由于研究对象是复合多相材料，具有更加优越的材料性能和组合功能，可以满足航空航天、火箭、导弹等军工制造中的结构性能设计需求。
(2)本发明引入全局径向基函数gsrbf对水平集函数进行插值，克服其在优化过程中需要重新初始化、速度扩展、迭代步长需满足逆风差分格式cfl条件及无法与许多成熟优化算法结合等数值缺陷。
附图说明
图1是本发明提供的考虑结构稳定性的多相材料拓扑优化设计方法的流程示意图；
图2是本发明实例悬臂梁初始设计域及初始化孔洞位置；
图3是本发明实例悬臂梁优化结果图；
图4是目标函数(柔度)迭代图；
图5是材料体积分数迭代图；
图6是第一阶失稳载荷因子迭代图；
图7是失稳载荷因子-柔度图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
请参阅图1，本发明提供的考虑结构稳定性的多相材料拓扑优化设计方法包括以下步骤：
(1)采用扩展系数插值隐式水平集函数，构建不同材料相的多材料结构；接着，采用max函数将所有材料相的水平集函数集成到一个水平集函数中，以得到集成有所有材料的多材料水平集函数。
具体地，通过隐式水平集函数构建不同材料相的结构，以插值隐式水平集函数的扩展系数为结构设计变量，在n个固定水平集结点处的结构隐式水平集函数为：
其中，x＝x1,x2,...,xn表示所有插值节点坐标，即水平集节点；n表示节点总数；αl,q表示在第l个水平集函数节点q处的扩展系数；φ表示组合材料1和材料2的水平集函数集，该水平集函数由高斯径向基函数φl,q(x)插值；φl,q(x)表示第l个高斯径向基函数，公式为：
其中，c是形状参数，等于水平集网格面积或者体积的倒数；xq表示水平集函数的第q个节点的坐标；||x-xq||是用于计算当前采样点x到xq节点距离的欧几里得范数。
进一步地，使用m个子水平集函数可以最多表示n＝2m个相。因此，可以定义设计空间：
其中，φ1(x)和φ2(x)分别表示两个不同的水平集函数，用来组合成结构中不同的材料相，u表示求和符号，ωk表示第k个材料相。
接着，使用max函数将所有函数集成到一个水平集函数φ(x,y)中：
φ(x,y)＝max(φ1,φ2,...,φm)
其中，m为所有水平集函数的数量。
(2)通过各个材料相的水平集函数的不同组合得到不同材料相的拓扑模型；接着，基于heaviside函数获得不同材料相连续变化的单元密度，并进一步求取各材料的体积。
具体地，通过水平集多材料理论构造结构中实际材料相和空相等材料相。本实施方式中以三相材料结构为例说明，但本发明所提供的方法不限于三相材料。通过不同的水平集函数的不同组合，分别得到不同材料相的拓扑模型如下：
其中，φ1(x)和φ2(x)分别表示两个不同的水平集函数。
接着；利用heaviside函数获得不同材料相连续变化的单元密度，不同材料相的特征函数χr(r＝1,2,3,4)为：
其中，φ1表示材料1的水平集函数，φ2表示材料2的水平集函数；h表示heaviside函数：
接着，各材料的体积计算公式如下：
(3)定义材料的物理参数、位移载荷和力载荷，并基于不同材料相连续变化的单元密度，将多材料结构中每种材料的弹性模量以特定方式进行插值处理，以得到多材料结构整体结构的等效弹性模量；
具体地，定义材料1的弹性模量为e1，定义材料2的弹性模量为e2，定义空相的弹性模量为e3，三种材料的泊松比均为ν，并将多材料结构的节点i的位移置为0，即
uix＝0,uiy＝0
其中，uix和uiy分别为节点i的x方向位移和y方向位移。
接着，将结构中各材料相的弹性模量互相插值得到结构的等效弹性模量e:
其中，er(r＝1,2,3,4)表示第r种材料相的弹性模量；h表示heaviside函数。
(4)基于参数化水平集的多材料理论建立考虑结构稳定性的多材料结构最小柔度拓扑优化模型，并基于上一步得到的等效弹性模量在结构设计域中通过有限元分析求解整体结构的位移场和应力场，根据得到的位移场计算多材料结构最小柔度拓扑优化模型的目标函数，根据应力场计算失稳载荷因子；接着，基于自伴随方法与伴随变量法对结构的设计变量进行灵敏度分析，并采用mma移动渐近线算法更新全局设计变量，继而确定满足稳定性约束的不同材料相的最优分布。
具体地，所述的基于参数化水平集理论建立的考虑结构稳定性的多材料最小柔度拓扑优化模型的表达式为：
find:αk＝[α1,α2,…,αm]t
subjectto:gi(φ)＝∫ωχi(φ)dω≤vi,max
αk,min≤αk≤αk,max
其中，m、αk、n和n分别为子水平集的数量、采用gsrbf插值时的扩展系数、有限元网格节点数量和材料相数量，j(u,φ)和gi(φ)分别表示结构的总柔度，及材料相i的体积，λ*和vi,max分别表示屈曲约束及材料相i的体积约束，λp表示在屈曲模态集j中第p阶失稳载荷因子，αk,max和αk,min分别表示设计变量的上下限。
在上式中，fi(u,u)代表插值前材料相i的应变能：
其中e(i)为在每种材料i的弹性模量。
χi(φ)表示第i种材料的特征函数，公式如下：
其中，hk是φk的heaviside函数h(φk)，
χ3(φ)＝h11h21＝h(φ1)h(φ2)
基于虚功原理，针对有限元平衡方程弱形式进行计算，对应的弱形式如下：
其中，a表示双线性能量式；l表示单线性荷载形式；dω为结构设计域的积分算子；h表示heaviside函数，用于表征结构形式的特征函数；ε为应变场；t表示矩阵的转置；u表示结构场的位移；v表示在动力学上允许的位移空间u中的一个虚拟位移；τ表示应用在边界
构建考虑结构稳定性的多材料最小柔度拓扑优化模型具体包括以下步骤：
(4.1)初始化设计参数，定义各相材料体积上限vi,max及整体结构失稳载荷因子下限λ*；
(4.2)基于各相材料物理参数插值得到的弹性模量e计算结构单元刚度矩阵，并组装成整体刚度矩阵k，接着计算节点位移u和结构单元应力σ，并计算整体几何刚度矩阵kg；
具体地，计算整体几何刚度矩阵公式如下：
(1)计算多材料单元应力σ＝χi(φ)dbu；
(2)基于多材料应力σ构建应力矩阵s，构建单元形函数相关的插值矩阵g；
(3)计算单元几何刚度矩阵kg：
kg＝∫χi(φ)gtsgdv
(4)根据节点耦合方式，将单元几何刚度矩阵kg组装成整体几何刚度矩阵kg；
其中，d表示平面应力问题下的弹性矩阵，b表示应变矩阵，u表示单元节点位移场，σx、σy、τxy分别表示单材料x方向正应力、y方向正应力、xy平面切应力，j表示雅克比矩阵，mi表示中间变量，pi表示单元形函数，ξ、η表示单元坐标系。
本发明首次基于参数化水平集方法求解失稳载荷因子，在求解过程中，首次使用多材料特征函数χi插值求解单元应力，进一步求解单元几何刚度矩阵，再根据节点耦合的方式将单元几何刚度矩阵组装成整体几何刚度矩阵，最后通过构建特征方程求解多材料结构失稳载荷因子，通过这样的方式求解，可以使多材料结构的失稳载荷因子计算的更加准确。
(4.3)第三步根据上一步计算得到的位移u计算结构整体的应变能，即目标函数j(u,φ)；接着使用结构单元应力σ计算失稳载荷因子λ。
具体地，计算多材料结构目标函数公式如下：
其中，ε为应变场；t表示矩阵的转置；u表示结构场的位移；v表示在动力学上允许的位移空间u中的一个虚拟位移；e是多材料结构的等效弹性模量；dω为结构设计域的积分算子；
接着，计算失稳载荷因子公式如下：
(k λkkg)φk＝0
其中，k和kg分别为单材料结构在外力荷载作用下由位移场计算得到的整体刚度矩阵和应力场计算得到的整体几何刚度矩阵，λk为单材料结构的第k阶失稳载荷因子，φk为对应的第k阶屈曲模态，即在第k阶极限载荷fk作用下单材料结构的位移。
(4.4)基于自伴随方法和伴随变量法求解目标函数与约束函数针对设计变量进行灵敏度分析，并采用mma移动渐近线算法更新全局设计变量，判断模型是否满足收敛条件，若否，则返回到步骤(4.2)，若是，则输出多材料结构的最优拓扑结构。
具体地，根据链式求导法则计算目标函数与约束函数针对优化设计变量的一阶微分，计算如下：
其中，其中，m、αk和n分别为子水平集的数量、采用gsrbf插值时的扩展系数和材料相数量，j(u,φ)和gi(φ)分别表示结构的总柔度及材料相i的体积，
其中，多材料失稳载荷因子针对设计变量的灵敏度求解过程具体为：
由于失稳载荷因子λ是针对多材料结构整体的，因此在优化求解过程中不用写成每个材料相累加的形式，在求解过程中，本发明引入关于失稳载荷因子的平衡方程弱形式公式：
a(u,v)＝λ·b(u,v)(3-43)
其中:
a(u,v,φ)＝∫ωεt(u)dε(v)χi(φ)dω(3-44)
b(u,v)＝∫ωσuvdω(3-45)
其中，ε为结构在受外载荷作用下的应变场，d为由组成结构的材料参数计算出的结构弹性矩阵，σ为结构的单元应力。
a(u,v,φ)和l(v,φ)对时间t求导得：
将平衡方程弱形式两边对时间t求偏导得到：
提取所有含
将式(3-46)、(3-47)和(3-49)代入式(3-48)可得：
使用伴随变量方法避免直接对u的求导，使用伴随变量w替换公式(3-50)中的v，可得
伴随变量w可通过求解下列双线性形式的伴随方程获得，
将式(3-53)代入公式(3-51)中，得：
将水平集函数速度场vn代入公式(3-54)可得：
通过链式法则对屈曲约束直接求其关于时间变量t的偏导数可得：
对比公式(3-55)和(3-56)可得屈曲约束关于扩展系数αk的导数：
本发明首次提出基于参数化水平集方法的多材料失稳载荷因子灵敏度的分析，在本发明中，通过平衡方程弱形式和伴随变量法求解多材料失稳载荷因子灵敏度，这样求解出的灵敏度的梯度信息更加准确，可以使多材料结构在拓扑优化过程中更快的收敛，提高优化效率。
实例
在对结构进行有限元分析时，其包含的单元通常有两种状态：受拉伸和受压缩，这也就是屈曲分析中出现正/负屈曲载荷因子的原因，当结构受力相反时会得到负的屈曲载荷因子，但由于在整个优化过程中，荷载的加载方向始终保持不变，因此，在本算例中，应该将负的失稳载荷因子从所提取的屈曲载荷集中排除，只对第一个正屈曲载荷因子施加约束限制。
本发明给出一个二维具有三相材料的短悬臂梁的算例来说明采用本研究方法的有效性。如图2中(a)所示，设计域被建模为单位厚度的矩形空间，并且通过宽度和高度均为单位大小的四节点等参元离散成n＝20*40的有限元网格(本研究考虑平面应力问题)。梁的左边缘是固定的，在右边缘的中点处施加一单位集中力。本小节所有优化结果均是基于图2中(b)所示初始孔洞位置。三种材料的弹性模量分别设置为e1＝10-6，e2＝1和e3＝2，其中e1表示空相，此外，对于不同材料的泊松比均设置为0.3。无量纲材料特性的设置是为了方便不同设计之间的比较，以下算例中均采用这种策略。每种材料的体积分数、颜色和从以最大化刚度为目标的多材料优化结果中获得的第一阶失稳载荷因子如表1所示。
表1多相材料参数表
受屈曲约束的悬臂梁拓扑优化结果如图3所示，从图中可以看出，随着第一阶失稳载荷因子的增大，多材料悬臂梁的强材料会分布在其受压部位，并且在悬臂梁的支撑边缘的材料相比非屈曲拓扑形状会变得越来越接近，此外，为避免发生屈曲变形，压杆尺寸往往大于拉杆的尺寸。从结果图中还可以看到，基于参数化水平集方法的多材料稳定性优化的结果边界更加光滑，更重要的是通过引入径向基函数、小波变换等数值处理方法，使得优化效率更高。
通过对比图4和5可知，施加屈曲约束前后，结构的柔度和各相材料的体积分数在迭代过程中变化平稳，并无明显不同。图6为施加屈曲约束前后第一阶失稳载荷因子的迭代曲线，从图中可以看出，当施加屈曲约束后，第一阶失稳载荷因子在前期的数值波动较大，这是由于基于参数化水平集的多材料稳定性方法在优化过程中与初始孔洞位置相关，由于初始孔洞位置分布均匀，因此在水平集函数演变从对称结构演变成强材料下移的过程中，容易出现结构突变，因此在前期会出现数值较大波动，当大部分强材料下移结束后，结构变化趋于稳定，因此在迭代后期可以看到第一阶失稳载荷因子呈现上升趋势，最终趋于平稳。
图7为悬臂梁第一阶失稳载荷因子和柔度变化关系图，如图所示，结构的稳定性与柔度变化规律相似，随着结构的第一阶失稳载荷因子增大，结构的柔度也随之增大，即结构的刚度越来越小，从中也可以发现，结构的稳定性与刚度是呈现负相关特性。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
1.一种考虑结构稳定性的多相材料拓扑优化设计方法，其特征在于，包括：
s1.以扩展系数作为结构设计变量，通过插值隐式水平集函数构建不同材料相的结构，求取各相材料的体积；
s2.构建以各扩展系数为设计变量，以多相材料结构总柔度最小化为优化目标，以各相材料的体积和失稳载荷因子为优化约束条件的拓扑优化模型；
s3.在当前扩展系数和水平集函数下，计算多相材料整体结构失稳载荷因子；
s4.计算目标函数与优化约束条件对设计变量的敏感度，根据计算结果更新扩展系数和拓扑优化模型；
s5.判断拓扑优化模型是否满足收敛条件，若是，则输出多相材料结构的最优拓扑结构；若否，则返回到步骤s3。
2.根据权利要求1所述的一种考虑结构稳定性的多相材料拓扑优化设计方法，其特征在于，所述拓扑优化模型的表达式为：
find:αk＝[α1,α2,…,αm]t
subjectto:gi(φ)＝∫ωχi(φ)dω≤vi,max
αk,min≤αk≤αk,max
其中，m、αk和n分别为水平集的数量、扩展系数和材料相数量，j(u,φ)和gi(φ)分别表示结构的总柔度，及材料相i的体积，λ*和vi,max分别表示屈曲约束及材料相i的体积约束，λp表示在屈曲模态集j中第p阶失稳载荷因子，αk,max和αk,min分别表示设计变量的上下限，fi(u,u)代表插值前材料相i的应变能，χi(φ)表示第i种材料的特征函数。
3.根据权利要求2所述的一种考虑结构稳定性的多相材料拓扑优化设计方法，其特征在于，失稳载荷因子的计算公式为：
(k λkkg)φk＝0
其中，k和kg分别为单材料结构在外力荷载作用下由位移场计算得到的整体刚度矩阵和应力场计算得到的整体几何刚度矩阵，λk为单材料结构的第k阶失稳载荷因子，φk为对应的第k阶屈曲模态，即在第k阶极限载荷fk作用下单材料结构的位移。
4.根据权利要求3所述的一种考虑结构稳定性的多相材料拓扑优化设计方法，其特征在于，整体几何刚度矩阵的计算过程为：
01.计算多材料单元应力σ＝χi(φ)dbu；
02.基于多材料单元应力σ构建应力矩阵s，构建单元形函数相关的插值矩阵g；
03.计算单元几何刚度矩阵kg：
kg＝∫χi(φ)gtsgdv
04.根据节点耦合方式，将单元几何刚度矩阵kg组装成整体几何刚度矩阵kg；
其中，d表示平面应力问题下的弹性矩阵，b表示应变矩阵，u表示单元节点位移场，σx、σy、τxy分别表示单材料x方向正应力、y方向正应力、xy平面切应力，j表示雅克比矩阵，mi表示中间变量，pi表示单元形函数，ξ、η表示单元坐标系。
5.根据权利要求1-4任一项所述的一种考虑结构稳定性的多相材料拓扑优化设计方法，其特征在于，基于自伴随方法和伴随变量法计算目标函数与优化约束条件对设计变量的敏感度。
6.根据权利要求5所述的一种考虑结构稳定性的多相材料拓扑优化设计方法，其特征在于，失稳载荷因子对设计变量的灵敏度为：
λ为失稳载荷因子，ε为应变场；t表示矩阵的转置，u表示结构场的位移，d是弹性矩阵，σ为结构的单元应力，w为伴随变量，δ表示dirac函数，为heaviside函数的一阶微分，
7.根据权利要求1-6任一项所述的一种考虑结构稳定性的多相材料拓扑优化设计方法，其特征在于，采用mma移动渐近线算法更新扩展系数。
8.根据权利要求1-7任一项所述的一种考虑结构稳定性的多相材料拓扑优化设计方法，其特征在于，引入全局径向基函数gsrbf对水平集函数进行插值。
9.一种考虑结构稳定性的多相材料拓扑优化设计系统，其特征在于，包括：计算机可读存储介质和处理器；
所述计算机可读存储介质用于存储可执行指令；
所述处理器用于读取所述计算机可读存储介质中存储的可执行指令，执行权利要求1至8任一项所述的考虑结构稳定性的多相材料拓扑优化设计方法。
技术总结