CN108710736B

CN108710736B - 一种考虑应力约束的宏微一体结构拓扑优化方法

Info

Publication number: CN108710736B
Application number: CN201810444607.7A
Authority: CN
Inventors: 邱志平; 王磊; 刘东亮; 夏海军
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2018-05-10
Filing date: 2018-05-10
Publication date: 2020-06-30
Anticipated expiration: 2038-05-10
Also published as: CN108710736A

Abstract

本发明公开了一种考虑应力约束的宏微一体结构拓扑优化方法，该方法针对同时具有宏观实体材料和微桁架的结构拓扑优化问题，采用代表体元法将微桁架等效表示为均值材料，并针对实体材料使用应力综合函数对其应力水平进行表征，然后采用伴随向量法求解了应力综合函数以及位移对设计变量的偏导值。最后，通过构建双层级的优化模型对考虑应力约束的宏微一体结构拓扑优化进行求解，其中，内层采用移动渐近线方法进行求解，实现结构的位移约束以及实体材料的应力约束，外层采用一元函数零点求解算法实现对微桁架的应力约束。所提供的实施例表明，拓扑优化方法能够有效控制宏微一体结构的应力水平，实现应力约束下的宏微一体结构拓扑优化。

Description

一种考虑应力约束的宏微一体结构拓扑优化方法

技术领域

本发明涉及连续体结构拓扑优化设计技术领域，特别涉及一种考虑应力约束的宏微一体结构拓扑优化方法。

背景技术

结构优化，特别是形状和拓扑优化，已被确定为结构设计中最具挑战性的任务之一。在过去二十年中，人们已经开发了各种技术和方法用于结构的拓扑优化。其中的代表性方法有水平集法、ESO法(evolutionary structural optimization)、SIMP法(solidisotropic material with penalization)等。随着拓扑优化技术的广泛使用以及3D打印成型技术的日益成熟，多材料结构也成为了可能，其设计也日益受到人们的重视，尤其是在航空航天结构中，由于多材料结构优异的单位质量性能，其使用能够很大程度上减小结构的质量，提高航空航天飞行器的承载能力。

值得注意的是，传统的多材料拓扑优化往往是基于各向同性材料的。但是，在当前的实际应用中，宏观实体材料和微桁架材料相结合的结构形式是最为广泛的，尤其是在卫星结构中的应用。然而，微桁架材料由于其单元微桁架的结构形式，难以等同于各向同性材料，常用的微桁架构型都是等同于正交各项异性材料的。因此，考虑材料的各向异性，并且能够反映微桁架特性的拓扑优化技术，有很大的应用前景和工程价值。

在实际应用中，结构优化的首要约束应当是应力，即使在多材料结构中也是如此，然而，由于多材料结构中应力约束的复杂性，在现有文献中还未有人能实现这一点。但是，考虑到宏微一体结构中实体材料和微桁架材料刚度性能差异较大，结构的主要载荷由实体材料承担，其应力较大的区域也将分布于实体材料中，因此，可以考虑将实体材料和微桁架的应力约束分离开，进行双层级的优化，内层实现实体材料结构的拓扑优化，外层优化通过调整实体材料的应力约束值，使得结构的微桁架材料应力满足要求。

发明内容

本发明要解决的技术问题是：克服现有技术的不足，提供一种考虑应力约束的宏微一体结构拓扑优化方法。本发明考虑实际工程中实体结构和微桁架应力分布的不同，以适用于实体材料和微桁架材料共存的双插值模型为基础，以结构的位移和应力为约束，最终优化得到的拓扑结果能够在满足结构位移和应力约束的基础上，充分利用两种不同材料的性能，获得传统单一材料所难以获得的结构性能。该方法考虑实体材料和微桁架材料同时在结构中存在的情况，采用代表体元法将微桁架等效表示为均值材料，并针对实体材料使用应力综合函数对其应力水平进行表征，最后，通过构建双层级的优化模型对考虑应力约束的宏微一体结构拓扑优化进行求解，其中，内层采用移动渐近线方法进行求解，实现结构的位移约束以及实体材料的应力约束，外层采用一元函数零点求解算法实现对微桁架的应力约束。最终实现应力约束下的宏微一体结构拓扑优化，为宏微一体结构设计提供一种新的设计方法。

本发明采用的技术方案为：一种考虑应力约束的宏微一体结构拓扑优化方法，可以对同时具有宏观实体材料和微桁架材料的结构进行应力约束拓扑优化，其实现步骤如下：

步骤一：基于应变能等效的原则，采用代表体元法将微桁架材料等效为均值材料：

其中，E₁、E₂和E₃为三向弹性模量，G₁₂、G₂₃和G₃₁为剪切模量，

为泊松比，D^H为等效弹性矩阵，

与

分别为单元的平均应力张量与平均应变张量，

分别为三个方向的正应变，

和

分别为三个方向的剪应变，

分别为三个方向的正应力，

和

分别为三个方向的剪应力。

步骤二：基于传统的带惩罚因子的材料插值模型，构建适用于宏微一体结构的双材料插值模型：

其中

为插值后第i个单元的单元刚度矩阵，

为对应于实体材料的单元刚度矩阵，

为对应于微桁架等效材料的单元刚度矩阵，x_1,i和x_2,i分别为第i个单元的设计变量1和设计变量2，p(p>1)为惩罚因子；

步骤三：构建考虑应力约束的宏微一体结构拓扑优化数学模型，以最小化结构质量作为目标，以结构应力和位移为约束：

其中，M是结构质量，ρ₁和ρ₂分别为实体材料和微桁架材料的密度，v_1,i和v_2,i分别为第i个单元的实体材料和微桁架材料的体积，n为设计域划分的单元总数，

为实体材料的许用应力，σ_1,i和f_i分别为第i个单元的实体材料和微桁架材料的许用应力，u_j和u_j,targ分别为第j个位移及其约束，m为位移约束的个数，x_1,i为第i个单元的设计变量1，x_2,i为第i个单元的设计变量2；x ₁和x ₂分别为设计变量1和设计变量2的下界，

和

分别为设计变量1和设计变量2的上界；

步骤四：基于伴随向量法求解位移约束对设计变量的敏度：

其中K_1,j为实体材料所对应的第j个单元的刚度矩阵，K_2,j为微桁架材料所对应的第j个单元的刚度矩阵，λ_j为伴随向量；

步骤五：基于伴随向量法求解应力综合函数对设计变量的敏度：

其中，f_s为应力综合函数，D_1,i为第i个单元的弹性矩阵，B_i为第i个单元的应变矩阵，p_s为应力综合函数罚值，σ_s,i为第i个单元中心点的米塞斯应力。

步骤六：使用移动渐近线优化算法(Method of Moving Asymptotes)，以最小化结构质量为目标，以结构的位移和应力为约束，利用质量、位移和应力对设计变量的灵敏度进行迭代求解，在迭代过程中，如果当前设计不满足位移和应力约束，或前后两个迭代步之间设计变量的变化量的绝对值之和大于预设值ε时，则返回步骤二进行新一轮的迭代优化，否则，进行步骤七；

步骤七：如果当前设计满足位移和实体材料的应力约束，而且前后两个迭代步之间设计变量的变化量的绝对值之和小于预设值ε时，则内层优化迭代结束，得到宏微一体结构位移和实体材料应力约束下质量最小的结构构型。

步骤八：采用微单元应力强度校核公式对微单元的应力进行校核：

其中，n_N和m_N分别为无量纲的轴力和弯矩，f为校核函数。

步骤九：如果微桁架应力满足要求(f<0)，则优化结束，否则采用二分法调整实体结构的应力许用值，并返回步骤二进行新一轮的迭代。

本发明与现有技术相比的优点在于：

本发明提供了一种考虑应力约束的宏微一体结构拓扑优化的新思路，通过构建适用于宏微一体结构的双材料插值模型，将实体材料的应力和微桁架材料的应力划分为两个不同的优化问题，最终实现了位移和应力约束下的宏微一体结构拓扑优化，可以获得更优的多材料结构拓扑形式。解决了传统多材料拓扑优化只适用于各向同性材料，或者无法进行应力约束的难题。该方法在确保结构位移和应力满足约束条件的同时，通过实体材料和微桁架的组合，有效提升了结构的力学性能，为航空航天飞行器的结构设计提供了理论方法。

附图说明

图1是本发明考虑应力约束的宏微一体结构拓扑优化流程图；

图2是本发明计算微桁架截面应力的示意图，其中，图2(a)是完全弹性阶段，图2(b)是完全塑性阶段，图2(c)是截面轴力，图2(d)是截面弯矩，图2(e)是求解弯矩的等效方法，图2(f)是求解弯矩的等效方法；

图3是本发明的应力约束宏微一体拓扑优化实施例的初始模型示意图；

图4是本发明的应力约束宏微一体拓扑优化实施例的所采用的微桁架单胞示意图；

图5是本发明针对接头结构的应力约束宏微一体拓扑优化结果示意图，其中，图5(a)是优化结果正视图，图5(b)是优化结果后视图；

图6是本发明实施例优化得到的接头拓扑结构设计域应力分布云图；

图7是本发明实施例优化得到的接头拓扑结构的微桁架等效单元应力分布云图。

具体实施方式

下面结合附图以及具体实施例进一步说明本发明。

如图1所示，本发明提出了一种考虑应力约束的宏微一体结构拓扑优化方法，包括以下步骤：

(1)基于应变能等效的原则，采用代表体元法将微桁架材料等效为均值材料：

微桁架材料可由典型单胞在空间周期性排布形成，通过取出一定数量的单胞，施加特定位移及载荷边界条件，在应变能等效的意义下使之等效为均质各向异性材料，进而可确定材料的等效宏观性能，这是代表体元法的思想。微桁架的宏观应力和宏观应变可以通过在体积上均匀化代表体元的应力和应变张量来得到，继而基于能量等效可计算材料等效弹性性能。在体积上进行均匀化后的平均应力和平均应变可表示为：

其中，σ_ij和ε_ij分别为代表体元内真实应力张量与应变张量，V_RVE为代表体元体积。

等效材料中的应变能为：

代表体元中的应变能为：

将(2)与(3)相减，得：

根据高斯积分，(4)中的体积积分可化为表面积分：

其中，S为外表面，n_j为单位外法向量，在S上有：

因此，有：

U′-U＝0 (7)

至此，可证明通过体积平均化应力与应变所获取的等效材料是能够保证与代表体元在应变能上是等效的。等效材料的本构方程可表示为：

为泊松比，D^H为等效弹性矩阵，

与

分别为平均应力张量与平均应变张量。

对于代表体元进行有限元分析可得到其内部的应力场和应变场，对应的平均值若直接通过体积平均求解，其过程在单元数较多的情况下将较为繁琐，利用高斯积分可转化为利用代表体元在表面边界上的位移求解，有：

至此，通过对微桁架单胞施加不同典型加载方式与周期性边界条件，可通过有限元技术求解得到相应的平均应力与平均应变，代入(8)即可得到特定微桁架材料的等效弹性性能。经过以上过程，代表体元的等效弹性模量(E₁、E₂、E₃、G₁₂、G₂₃和G₃₁)可通过有限元分(FEM)获取。

(2)基于传统的带惩罚因子的材料插值模型，构建适用于宏微一体结构的双材料插值模型：

其中

为插值后第i个单元的单元刚度矩阵，

为对应于实体材料的单元刚度矩阵，

(3)构建考虑应力约束的宏微一体结构拓扑优化数学模型，以最小化结构质量作为目标，以结构应力和位移为约束：

其中，M是结构质量，ρ_1,i和ρ_2,i分别为实体材料和微桁架材料的密度，M为位移约束的个数，n为设计域划分的单元总数，

和

分别为实体材料和微桁架材料的许用应力，x_1,i为第i个单元的设计变量1，x_2,i为第i个单元的设计变量2；x ₁和x ₂分别为设计变量1和设计变量2的下界，

和

分别为设计变量1和设计变量2的上界；

(4)基于伴随向量法求解位移约束对设计变量的敏度：

由于通过差分的形式求解约束函数对设计变量的偏导数

和

将会导致很大的计算量，因此需要通过间接的方式求解，构造如下约束函数的拉格朗日增广函数：

其中，λ_j(j＝1,2,…,m)为拉格朗日乘子向量，即伴随向量，向量F为结构载荷向量，K为刚度矩阵，u为位移向量。由平衡方程F-Ku＝0得

因此式(10)对设计变量x_k,i(k＝1,2)求全导数得：

其中，

上式对任意λ均成立，因此通过调整λ的值，可以使得du/dx_i,k所在项的系数为零，即令：

由于刚度矩阵具有对称性，上式可以改写为：

式(14)表明，通过给结构施加虚拟载荷

所求得的位移即为伴随向量λ。得到λ的值后，结构位移对设计变量x_k,i(k＝1,2)的偏导数则由下式给出：

其中λ_j、K_j、u_j分别为对应于u_j的伴随向量、总体单元刚度矩阵和位移向量。由于载荷向量F与设计变量无关，即dF/dx_k,i＝0，所以式(15)可以写为：

由宏微一体结构的双材料插值模型，可得dK_j/dx_k,i为：

所以式(16)最终写为：

对于第i(i＝1,2,…,n)个单元，考虑对应力进行一定程度上的松弛，设其单元应力分量为：

设应力综合函数为：

其中，σ_0,1为材料1的应力许可值，p_s为应力综合函数罚值，一般取p_s＝30或更高。

应力综合函数f_s对设计变量x_i的偏导数为：

将应力计算公式代入式(21)，得：

式(22)的上下两部分可分别可整理为：

其中N为结构自由度数，a_1,k和a_2,k(k＝1,2,…,N)为系数，由式(22)决定。则有：

其中，

和

由式(17)确定。根据线弹性结构的位移叠加原理可知，

和

可以通过一次伴随向量求出。

(6)使用移动渐近线优化算法(Method of Moving Asymptotes)，以最小化结构柔度为目标，以结构相对质量分数为约束，利用柔度和相对质量分数对涉及变量的灵敏度进行迭代求解，在迭代过程中，如果当前设计不满足相对质量约束，或前后两次迭代的设计变量变化绝对值之和大于预设值ε时，则返回步骤二进行新一轮的迭代优化，否则，进行步骤六；

(7)如果当前设计不满足相对质量约束，或前后两次迭代的设计变量变化绝对值之和大于预设值ε时，则内层优化迭代结束，得到宏微一体结构实体应力约束下质量最小的结构构型。

(8)采用微单元应力强度校核公式对微单元的应力进行校核：

假设材料为理想塑性材料，屈服强度为σ_s。对于直径d＝2r的圆截面均质梁，假设该截面形成塑性铰时，认为其发生了屈服，此时，中性轴的偏移假设为h，该截面上的应力分布如图2所示，其中F为半圆截面CGD的质心，E为截面ABCD的质心。图2(b)为截面塑性铰形成时截面的应力分布，其可以分解为截面轴力F_N和截面弯矩M_N，分别如图2(c)和图2(d)所示。

截面轴力F_N计算较为简单，即：

F_N＝2·S_ABCD·σ_s (27)

式中截面ABCD的面积S_ABCD为：

对于截面弯矩M_N，如图2(d)、(e)、(f)解释的等效求解方法，即：

M_N＝(2·S_CDG·y_F-2·S_ABCD·y_E)σ_s (29)

为了简便，可对截面轴力和弯矩进行下面无量纲化处理：

通过整理式(27)～(31)，得到截面发生屈服时，截面轴力F_N和截面弯矩M_N的如下关系式：

进一步，可得到截面的屈服函数：

代入截面的轴力和弯矩，当f<0时，截面不发生屈服，当f>0时，截面屈服，微单元失效，其中，n_N和m_N分别为无量纲的轴力和弯矩，f为校核函数。

(9)如果微桁架应力满足要求，则优化结束，否则采用二分法调整实体结构的应力许用值，并返回步骤二进行新一轮的迭代。

实施例：

了更充分地了解该发明的特点及其对工程实际的适用性，本发明针对如图3所示的某型飞机舱门的接头结构实施宏微一体结构应力约束拓扑优化。接头结构的初始设计域为图3中灰色部分，深灰色网格部位为非设计域。接头的上端施加了如图3所示的y方向载荷F＝45000kN，下端四个圆孔固支。设计域所采用的两种材料分别为：铝合金(弹性模量E＝70GPa，泊松比μ＝0.3，密度ρ＝2700kg/m³)，铝合金粉末3D打印成型的微桁架(弹性系数为E_x＝E_y＝563Mpa、E_z＝4595Mpa、μ_xy＝0.8625、μ_yz＝μ_xz＝0.0574、G_xy＝3016Mpa、G_yz＝G_xz＝2952Mpa，密度ρ＝315.09kg/m³)。非设计域采用了铝合金材料。微桁架的单胞构型如图4所示，单胞的尺寸为4×4×4mm，单胞杆径为0.3mm。铝合金材料的许用应力为280MPa，由该许用应力计算，当微桁架等效单元的许用应力小于42.15MPa时，微桁架应力校核函数f<0。接头的位移约束为1.5mm。

经过所提的应力约束宏微一体结构拓扑优化，得到结果如图5～图7所示，图中设计域内的深色部位是微桁架，灰色部位为铝合金实体结构。从结果的拓扑构型上看，实体材料位于力的主要传递路径上，微桁架则起到辅助支撑的作用，结构的整体布局较为合理。优化后得到的接头结构加载点的位移为1.497mm，满足设计的约束要求。优化后结构设计域的最大米塞斯应力为233.913MPa，小于280MPa的需用位移。微桁架的最大米塞斯应力为36.37MPa，小于应力校核函数f<0所对应的最小米塞斯应力(42.15MPa)，所以微桁架的应力也满足设计约束。因此，本算例中考虑应力约束的宏微一体结构拓扑优化方法实现了实体材料和微桁架材料应力约束以及结构位移约束条件下的拓扑优化，在降低结构质量的同时，提升了结构的性能。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制；其可扩展应用于宏微一体结构应力约束拓扑优化设计领域，凡采用等同变换或者等效替换而形成的技术方案，均落在本发明权利保护范围之内。

本发明未详细阐述部分属于本领域技术人员的公知技术。