一种位移和应力混合约束下的连续体结构非概率可靠性拓扑优化方法与流程

本发明涉及含连续体结构的拓扑优化设计领域，特别涉及考虑材料属性和载荷环境的不确定性对结构的刚度和强度的影响以及基于位移和应力的非概率可靠度指标约束下连续体结构的不确定拓扑优化方案的制定。

背景技术：

随着科学技术和生产力的日益发展，人类伸向太空、海洋的触角不断延伸，结构优化研究的范围越来越宽泛。因资源的有限、工程技术的激烈竞争和环境需要保护等问题，使得结构优化设计变得越来越重要。结构优化设计分为三个层次：尺寸优化、形状优化和拓扑优化。与尺寸优化和几何优化相比，结构拓扑优化不仅待确定的参数更多，而且拓扑变量对优化目标的影响更大，因而取得经济效益更大。因此，针对连续体结构的拓扑优化研究具有重要的理论意义和工程实用价值。

然而，随着科技水平的不断进步，工程结构日趋精密和复杂，通常涉及多个领域中大量信息的交换。此外，结构服役的环境也越来越恶化，材料的制造加工工艺造成的材料属性的分散性也不可避免，这些不确定因素会对结构的工作性能产生重要影响。拓扑优化作为结构优化设计的概念设计阶段，其优化设计结果对最终的结构形式有着决定性的影响，因此在拓扑优化设计阶段考虑不确定性的影响是十分必要的。实际结构中应力和位移约束是十分重要的，不考虑它们的设计是不能付诸工程使用的。因此，研究位移和应力混合约束下的连续体结构非概率可靠性优化设计方法具有重大意义。

当前，国内外学者与工程技术人员对考虑连续体结构的拓扑优化分析与设计研究主要集中在两个方面：(1)以结构柔顺度为目标函数，受体积约束的确定性拓扑优化问题；(2)以结构重量为目标函数，以结构位移为约束的非概率可靠性拓扑优化问题。上述工作一定程度上丰富了连续体结构的拓扑优化设计研究。但是上述工作并没有考虑位移和应力两种约束下的连续体结构非概率可靠性拓扑优化设计方法，而且已经提出的非概率可靠性拓扑优化设计方法使得结构的安全余量过大，结构的经济效益受损，大大限制了其理论的工程实用化进程。

由于实际工程中贫信息、少数据的情况时有发生，建立以非概率理论框架为基础的位移和应力混合约束下的连续体结构拓扑优化设计方法具有显著的现实意义。目前，相关研究工作尚不成熟，现有连续体结构的拓扑优化设计方案经常无法严格满足所需的应用要求，亦或是安全冗余度过大，造成严重的资源浪费与时间成本损耗。

技术实现要素：

本发明要解决的技术问题是：克服现有技术的不足，提供一种位移和应力混合约束下的连续体结构非概率可靠性拓扑优化方法。本发明充分考虑实际工程问题中普遍存在的不确定性因素，以提出的非概率可靠性度量指标作为优化模型的约束条件，所得到的设计结果更加符合真实情况，工程适用性更强。

本发明采用的技术方案实现步骤如下：一种位移和应力混合约束下的连续体结构非概率可靠性拓扑优化方法，该方法包括如下步骤：

第一步：采用带罚因子的固体各向同性微结构/材料插值模型(SIMP模型)，考虑结构材料属性和载荷环境的不确定性并以区间量来度量这种不确定性，运用区间模型来描述位移不确定量以及应力不确定量，利用非概率可靠性度量指标，以结构的最小体积作为优化目标，以结构关键部位的位移或应力作为约束，建立相应的非概率可靠性拓扑优化模型如下：

其中，V是优化区域的体积，ρ_i和V_i分别为第i个单元的相对密度和体积，N为优化区域划分的单元总数，ρ为单元相对密度的下限。K为单元的总体刚度矩阵，u为单元的总体位移列向量，F为总体载荷列向量。是第j个位移约束点的实际位移区间值，是第j个位移约束的容许位移区间值，m为位移约束的个数。是第k个应力约束点的实际应力区间值，是第k个应力约束的许用应力区间值，n为应力约束的个数。R_s是非概率可靠性指标，是第j个位移约束对应的目标非概率可靠度，是第k个应力约束对应的目标非概率可靠度。对于SIMP模型，单元的弹性模量是材料相对密度的函数

其中P>1是惩罚因子，用于实现对中间密度单元的惩罚。按照经验，一般取P＝3，E₀是完全实心材料的弹性模量。

第二步：以区间量来表征材料的弹性模量和载荷的不确定性，以材料的弹性模量和载荷为变量，采用顶点组合法，对每个材料弹性模量和载荷的不确定量的顶点组合进行有限元计算，得到相应的位移约束点的位移和应力约束单元中心点的应力，然后进行比较，其中位移(应力)约束点的位移(应力)的最大值作为位移约束点位移(应力)的上界，位移(应力)约束点的位移(应力)的最小值作为位移约束点位移(应力)的下界。

第三步：根据得到的位移的上下界以及应力的上下界，结合非概率可靠性模型，分别得到位移和应力约束相应的非概率可靠度为

其中，R^I为位移(应力)的容许区间值，S^I为位移(应力)的实际区间值。R为位移(应力)的容许区间值的下界，为位移(应力)的容许区间值的上界。S为位移(应力)的实际区间值的下界，为位移(应力)的实际区间值的上界。

第四步：采用优化特征位移替代非概率可靠性指标来改善问题的收敛性。优化特征位移定义为实际失效平面到目标失效平面的移动位移，而目标失效平面是与原实际失效平面平行，其可靠度等于目标非概率可靠度的失效平面。利用优化特征位移可以将原优化模型改写为

其中，d为优化特征位移。

第五步：运用伴随向量法得到位移上下界以及应力上下界对设计变量的灵敏度，然后利用复合函数的求导法则，先求解优化特征位移关于位移(应力)上下界的灵敏度，然后再求解位移(应力)上下界关于设计变量的灵敏度，将两者相乘即得到位移和应力的优化特征位移对设计变量的灵敏度。

第六步：利用得到的优化特征位移值以及位移和应力的优化特征位移对设计变量的灵敏度构造原问题的近似模型，设置相关的经验参数，运用MMA算法对优化问题进行求解，得到新的设计变量。

第七步：迭代过程中，如果当前设计不满足可靠度约束的许用值，或者尽管满足可靠度约束，但相较于上一个可行解，目标函数的相对变化百分比大于预设值ξ时，则继续运用MMA算法进行迭代计算，将已经完成迭代次数的值增加一，并返回第二步，否则，迭代过程结束。

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

本发明提供了位移和应力混合约束下的连续体结构非概率可靠性拓扑优化新思路，在对连续体结构进行拓扑优化设计时，可以充分考虑不确定性对结构性能的影响，在保证结构刚度和强度满足可靠度约束的前提下可大大降低结构重量，提高性能的同时，降低设计周期和经济成本。与传统的考虑单一约束(柔顺度、位移)的拓扑优化方法相比，本方法考虑了结构刚度和强度两种约束下的拓扑优化设计方法，本方法和实际工程结合更加紧密，具有更加重大的意义。

附图说明

图1是本发明针对位移和应力混合约束下的连续体结构非概率可靠性拓扑优化设计流程图；

图2是本发明提出的非概率可靠性模型中的应力-强度非概率集合干涉模型示意图；

图3是应力-强度干涉模型的标准化空间示意图；

图4是本发明提出的优化特征位移的临界斜率示意图；

图5是本发明实施例中拓扑优化设计区域以及边界和载荷条件示意图；

图6是本发明实施例中确定性拓扑优化和不同可靠度约束下拓扑优化的最终构型图，其中，图6(a)为确定性优化结果，图6(b)中R_s＝0.90，图6(c)中R_s＝0.95；

图7是发明实施例中确定性拓扑优化和不同可靠度约束下拓扑优化的目标函数的迭代历程曲线。

具体实施方式

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

如图1所示，本发明提出了一种位移和应力混合约束下的连续体结构非概率可靠性拓扑优化方法，包括以下步骤：

(1)采用带罚因子的固体各向同性微结构/材料插值模型(SIMP模型)，考虑结构材料属性和载荷环境的不确定性，运用区间模型来描述位移不确定量以及应力不确定量，利用非概率可靠性度量指标，以结构的最小体积作为优化目标，以结构关键部位的位移或应力作为约束，建立相应的非概率可靠性拓扑优化模型如下：

其中，V是优化区域的体积，ρ_i和V_i分别为第i个单元的相对密度和体积，N为优化区域划分的单元总数，ρ为单元相对密度的下限。K为单元的总体刚度矩阵，u为单元的总体位移列向量，F为总体载荷列向量。是第j个位移约束点的实际位移区间值，是第j个位移约束的容许位移区间值，m为位移约束的个数。是第k个应力约束点的实际应力区间值，是第k个应力约束的许用应力区间值，n为应力约束的个数。R_s是非概率可靠性指标，是第j个位移约束对应的目标非概率可靠度，是第k个应力约束对应的目标非概率可靠度。对于SIMP模型，单元的弹性模量是材料相对密度的函数：

其中P>1是惩罚因子，用于实现对中间密度单元的惩罚。按照经验，一般取P＝3，E₀是完全实心材料的弹性模量。

(2)以区间量来表征材料的弹性模量和载荷的不确定性，以材料的弹性模量和载荷为变量，采用顶点组合法，进行有限元计算，得到位移约束点的位移和应力约束单元中心点的应力，进行比较得到位移约束点位移的上下界和应力约束点应力的上下界。

首先介绍计算位移上下界的方法。在有限元位移法中，静力平衡方程为：

Ku＝F (3)

考虑结构材料和载荷的不确定性，在平衡方程中，将刚度矩阵K和载荷向量f考虑为区间变量K^I和f^I。则形成的区间有限元控制方程为：

K^Iu^I＝f^I (4)

其中u^I为位移区间向量。由于控制方程是线性的，因此可以采用顶点组合法来求解u^I中任意分量的上下界。

顶点组合法：

如果f(x₁,x₂,…,x_n)对自变量x_i(i＝1,2,…,n)是单调的，将自变量考虑为区间变量时，即：

由函数的单调性可知，f的取值范围为：

其中r为顶点(区间两端点)组合序数，k_i＝1,2，i＝1,2,…,n；r＝1,2,…,2ⁿ。

所以，根据式(6)，得到第j个约束对应的位移的取值区间为：

其中其中下标corj表示位移区间向量u^I中对应于第j个位移约束的分量；上标k_i＝1,2，当k_i＝1时表示对应值取下界，当k_i＝2时表示对应值取上界，即i＝1,2,…,N。

然后介绍应力上下界的求解方法。为了表述方便，记：

σ_r＝h(u^e) (8)

式中，u^e表示应力约束单元所对应的结点的位移分量组成的位移列向量。它与总体位移列向量有如下关系：

u^e＝s_eu (9)

式中s_e为单元位移提取矩阵。对于矩形双线性单元，有：

(u^e)_8×1＝(s_e)_8×nu_n×1 (10)

其中矩阵(s_e)_8×n的每一行为某个元素为1，其余元素为0的行向量。

根据顶点组合法，得到单元应力的取值区间为：

其中其中上标k_i＝1,2，当k_i＝1时表示对应值取下界，当k_i＝2时表示对应值取上界，即(K^-1)²＝K^-1，，i＝1,2,…,N。

(3)根据得到的位移的上下界以及应力的上下界，结合非概率可靠性模型，得到位移和应力约束的非概率可靠度为：

其中，R^I为位移(应力)的容许区间值，S^I为位移(应力)的实际区间值。

(4)在运用MMA算法求解拓扑优化问题时，第2节所述非概率可靠度的梯度信息存在缺陷，即在设计域中存在梯度信息为零的区域(可靠度为0或1)，会造成一定的数值困难。下面引入优化特征位移d来解决这一困难。

优化特征位移d的定义为：原失效平面到目标失效平面的移动位移。其中目标失效平面是与原失效平面平行的平面，并且其可靠度为目标值。由于可靠度一般接近1，故目标失效平面一般位于不确定域的右下角，图4为两种临界情况。

首先计算临界情况下失效平面的斜率，设η为目标可靠度。对于k₁，有(2×2/k₁×1/2)/4＝1-η，解得k₁＝1/2(1-η)，同理可得k₂＝2(1-η)，针对原失效平面斜率k取值的不同情况，使用直线间的距离公式，并定义原失效平面在目标失效平面上方的距离为正，反之为负，给出距离d的表达式为：

当d>0时，失效平面在与目标非概率可靠度η对应的目标失效平面上方，此时由于安全区域的面积小于目标值，对应的非概率可靠度R_s＜η，不满足要求。当d≤0时，失效平面在与目标非概率可靠度η对应的目标失效平面下方，此时由于安全区域的面积大于等于目标值，对应的非概率可靠度R_s≥η，满足设计要求。

(5)运用伴随向量法得到位移上下界以及应力上下界对设计变量的灵敏度，然后利用复合函数的求导法则，先求解优化特征位移关于位移(应力)上下界的灵敏度，然后再求解位移(应力)上下界关于设计变量的灵敏度，将两者相乘即得到位移和应力的优化特征位移对设计变量的灵敏度。

伴随向量法：

考虑优化问题中任意一种非概率可靠性指标约束的情况。设约束函数为f，其对单个设计变量x_i的全导数为：

其中上式中用表示仅f显含x_i的部分对x_i进行求导。和通常可以直接计算，为了避免直接计算使计算效率更高，构造如下的约束函数的增广拉格朗日函数：

其中，λ为与平衡方程关联的任意的乘子向量，也被称为伴随向量。由于F-Ku＝0，故式(15)对设计变量x_i求全导数得：

式(16)对任意λ均成立，因此可以选取适当的λ使得du/dx_i所在项的系数为零，即令：

利用刚度矩阵的对称性可以将上式改为：

对比式(3)可知，可以仿造有限元位移求解得过程，将视为虚拟载荷，则可将虚拟位移λ求出，这可以利用有限元程序来求解。求解出λ后，约束函数的灵敏度则由下式给出：

由于约束函数为位移或应力，不显含设计变量x_i，而且载荷F不随设计变量变化，即dF/dx_i＝0，则式(19)可以改写为：

位移上下界灵敏度分析：

对于位移上界灵敏度分析，则构造位移上界对应的拉格朗日增广函数为：

式中，为位移上界对应的顶点组合中的载荷列向量，为位移上界对应的顶点组合中的刚度矩阵，为位移上界对应的顶点组合下的位移列向量。

运用伴随向量法，根据式(20)可得到：

式中，满足：

对于位移下界灵敏度分析，可采取同样的方法，得到：

式中，为位移下界对应的顶点组合中的刚度矩阵，为位移下界对应的顶点组合下的位移列向量。满足：

应力上下界灵敏度分析：

对于应力上界灵敏度分析，构造应力上界对应的拉格朗日增广函数为：

式中，为应力上界对应的顶点组合中的载荷列向量，为应力上界对应的顶点组合中的刚度矩阵，为应力上界对应的顶点组合下的位移列向量。

运用伴随向量法，根据式(20)可得到：

式中，满足：

式(28)中，的求解可利用复合函数的求导法则进行。

式(29)中，为应力上界对应的顶点组合下的应力约束单元的结点位移列向量。

由式(8)和式(9)，式(29)可转化为：

式(30)中，计算如下：

其中，

同理，对于应力下界灵敏度的分析，可得：

式中，为应力下界对应的顶点组合中的刚度矩阵，为应力下界对应的顶点组合下的位移列向量。满足：

式(33)中，为应力下界对应的顶点组合下的应力约束单元的结点位移列向量。

优化特征位移对设计变量的灵敏度求解方法：

根据复合函数的求导法则，其约束函数d对设计变量的灵敏度为：

其中：

式中而A,B中的X为

对于位移约束，将式(34)与式(22)及式(24)式结合，有：

对于应力约束，将式(34)与式(27)及式(32)结合，有：

(6)利用得到的优化特征位移值以及位移和应力的优化特征位移对设计变量的灵敏度构造原问题的近似模型，并运用MMA算法对优化问题进行求解，得到新的设计变量。

(7)迭代过程中，如果当前设计不满足可靠度约束的许用值，或者尽管满足可靠度约束，但相较于上一个可行解，目标函数的相对变化百分比大于预设值ξ时，则继续运用MMA算法进行迭代计算，将已经完成迭代次数的值增加一，并返回第二步，否则，迭代过程结束。

实施例：

为了更充分地了解该发明的特点及其对工程实际的适用性，本发明针对如图5所示的矩形区域进行拓扑优化设计。其中，矩形域的尺寸为25mm×20mm，厚度为0.25mm，划分为100×80个单元。材料为钢，取弹性模量E＝205Gpa，泊松比μ＝0.3。矩形区域的左端固定，在右端最下面的A点加上F＝20N的力，不考虑重力的影响，约束A点的位移和应力，使得u_A<0.55mm，σ_A<450Mpa，弹性模量E、载荷F及A点的位移约束u_A和应力约束σ_A相对名义值的波动为：即E＝[189,231]Gpa，F＝[18,22]N，u_A＝[0.5,0.6]mm，σ_A＝[400,500]Mpa。初始材料的相对密度为1，对该结构进行拓扑优化。

该实例给出了不考虑不确定性影响的拓扑优化结果以及采用两种不同水平的可靠度设计许用值进行约束，即R_s分别为0.9和0.95，图6给出了确定性拓扑优化和不同可靠度约束下拓扑优化的最终构型图。图7(a)给出了不考虑不确定性影响下目标函数的迭代历程曲线。图7(b)(c)给出了考虑不确定性影响下目标函数的迭代历程曲线。可以看出，可靠性拓扑优化的构型相对于确定性的拓扑优化的更加粗实一些，体积分数更高，结构更加安全。确定性的拓扑优化构型体积分数较低，经数据分析，确定性拓扑优化结果的可靠性只有45％左右，无法满足工程的实际要求。

综上所述，本发明提出了一种位移和应力混合约束下的连续体结构非概率可靠性拓扑优化方法。首先，根据Qiu的非概率可靠性模型建立位移和应力混合约束下的连续体结构非概率可靠性拓扑优化模型；其次，运用顶点组合法得到位移约束点位移的上下界和应约束点应力的上下界；运用伴随向量法求解位移上下界的灵敏度和应力上下界的灵敏度；采用优化特征位移代替原来的可靠性指标来改善问题的收敛性，并利用复合函数求导法则得到优化特征位移关于设计变量的灵敏度；最后，利用已有信息构造原问题的近似模型，采用MMA算法进行迭代计算，直至满足约束条件和收敛条件。

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

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