改进双向渐进法的等效静载荷法动态响应拓扑优化方法

摘要

本发明涉及一种改进双向渐进法的等效静载荷法动态响应拓扑优化方法。该方法将改进的BESO方法应用于等效静载荷法中作为优化迭代方法，替换掉原始的变密度法，并建立新的收敛条件，从而改进动态响应拓扑优化流程。本发明提出的方法，能有效减少k大于1次设计循环阶段的优化时间，两个算例中优化所需CPU总时间比原始方法分别减少了21.8%、25.3%；该方法具有计算速度快、优化流程简单的特点；随着外载荷作用点的增多、优化规模的增大，优化所需的设计循环次数也会增加，该方法将会带来更大的优势。

说明书

技术领域

本发明属于结构动态响应拓扑优化领域，具体涉及一种改进双向渐进法的等效静载荷法动态响应拓扑优化方法。

背景技术

拓扑优化是指在满足给定约束下寻求材料的最优分布形式^[1]。常见拓扑优化方法有变密度法、渐进结构法以及均匀化法。Bend-soe>[2]于1988年提出了均匀化法，Bend-soe>[3]于1989年提出了变密度法。渐进结构优化法(ESO)是近年来兴起的一种解决各类结构优化问题的数值方法，由Steven等人^[4]于1993年提出，该方法通过删除结构中无效或低效的材料逐步去掉使剩下的结构渐渐趋于优化。2000年Querin等^[5]采用双向渐进法(BESO)利用Von>[6]于2009年指出，采用两个独立的参数RR和IR是不合理的，如果没有选择好这两个参数将得不到最优的结果。鉴于此，Huang和Xie等人^[7]针对刚度优化，提出一种基于灵敏度分析的改进BESO法，该方法解决了连续体拓扑优化的很多问题，如优化问题的恰当描述、棋盘格现象、网格依赖性以及解的收敛性，使其优化结果更加趋于稳定。PARK>[8-9]先后将改进的BESO法用于汽车发动机机罩的多目标优化，取得了理想的结果。

目前拓扑优化多用于静态载荷优化，而工程实际中结构受到的一般是动态载荷，直接进行动态响应优化分析^[10-11]虽然在理论上可行，但是计算很难收敛，同时计算过程繁冗，耗时长，未能在工程中得到应用。针对这些问题，Park于2002年提出了等效静载荷法^[12](ESLM),将动态载荷转化为等效静态载荷，并作为边界条件对结构进行动态优化，该方法最初只是用于尺寸优化。在随后的研究中，PARK证明了该方法获得的结果满足Karush-Kuhn-Tucker必要条件^[13]，并通过相应的算例验证了该方法对线性结构和非线性结构进行拓扑优化的可行性^[14-15]。该方法在国内亦得到广泛应用^[16-20]。相较于传统动力学优化法，等效静载荷方法在计算效率上已经得到很大提升，但是由于优化时需要若干个设计循环才能得到稳定结果，且每次设计循环都要进行动力学分析以及静态拓扑优化迭代，依然会耗费大量时间。

为了进一步提升等效静载荷法的计算效率，Kim E等人^[21]将系统约简法与等效静载荷法结合，能够更有效地求解大规模动力学优化问题。Lee>[15]针对非线性结构提出一种结点近似位移法来减少动态响应分析阶段所花费的时间，以此提高整体优化效率。

发明内容

本发明的目的在于解决当前等效静载荷法中单次设计循环耗时长，优化流程复杂等问题，提供一种改进双向渐进法的等效静载荷法动态响应拓扑优化方法，该方法将改进的BESO方法应用到等效静载荷法中作为优化迭代方法，替换掉原始的变密度法，从而提升计算效率，简化优化流程。

为实现上述目的，本发明的技术方案是：一种改进双向渐进法的等效静载荷法动态响应拓扑优化方法，将改进的BESO方法应用于等效静载荷法中作为优化迭代方法，替换掉原始的变密度法，并建立新的收敛条件，从而改进动态响应拓扑优化流程。

在本发明一实施例中，所述改进的BESO方法具体如下，

在对单元处理之前，应事先给出每次迭代的单步目标体积V_k+1；初始体积V_k大于或小于整体的体积约束V*；因此，体积进化迭代式可表示为：

V_k+1＝V_k(1±ER),(k＝1,2,3...)>

式(1)中k为当前迭代次数，ER为体积进化率；

若目标体积V_k+1小于当前体积V_k，则必须将单元按灵敏度大小排列并删除灵敏度低的单元，在删除单元之前，空单元的灵敏度有可能大于待删单元灵敏度，因此要先考虑增添单元，这里引入增删阀值α^th，α^th定义为待删单元灵敏度中的最大值，有：

α_i>α^th>

式(2)中α_i为增添单元的灵敏度；当增添单元的体积大于体积增添率AR时，需要重新计算α^th直到满足增添单元的体积小于或等于AR为止；然后删除灵敏度小于该值的单元；显然，删除单元的体积V_del可以表示为：

V_del＝V_k-V_k+1+V_add>

式(3)中V_add为增添单元的体积；

若目标体积V_k+1大于当前体积V_k，则必须将单元按灵敏度大小排列并增添灵敏度高的单元，在增添单元之前，实单元的灵敏度有可能小于待增单元灵敏度，因此要先考虑删除单元，这里引入增删阀值β^th，β^th定义为待增单元灵敏度中的最小值，有：

β_i＜β^th>

其中β_i为删除单元的灵敏度；当删除单元的体积大于体积增添率AR时，需要重新计算β^th直到满足删除单元的体积小于或等于AR为止；然后增添灵敏度大于该值的单元；显然，增添单元的体积V_add可以表示为：

V_add＝V_k-V_k+1+V_del>

重复进行有限元分析与单元增删的过程，直至迭代循环满足两个终止条件：一是体积约束，即当前进化体积已满足整体约束体积V*；二是收敛准则，考虑目标函数的变化量，收敛判别式规定为：

式(6)中k为当前迭代步数，τ为内部许用收敛因子；M是一整数，用以限制并稳定柔顺度在多次连续迭代中的平均变化量；计入局部多个M值，能够使目标函数波动较小，曲线更为平滑，有效防止局部收敛。

在本发明一实施例中，所述M取值为5。

在本发明一实施例中，以最小化柔顺度为目标函数，以体积为约束条件，以单元本身为设计变量，基于灵敏度分析改进的BESO刚度优化数学表达式可描述为：

式(7)中F和U分别为外载荷向量和结点位移向量，C为结构柔顺度，V_i是个体单元的体积，V*为体积约束，N为设计域中单元总数量；x_i为双向赋值的非连续设计变量，即材料单元的密度；K为结构总体刚度矩阵，x_min为材料单元的密度的一个小值。

在本发明一实施例中，所述改进的BESO方法中，删除单元的方法为软杀法，即将灵敏度低的待删单元的密度人为减至一个小值x_min，以防止刚度矩阵奇异，可令x_min＝0.001，但后续仍参与有限元分析，而不是直接删除；因此，采用单元软杀法的BESO方法的单元灵敏度计算式可归结为：

式(8)中K_i⁰表示实单元的单元刚度矩阵，u_i为单元节点位移向量,P为惩罚因子，可令P＝3。

在本发明一实施例中，将权值法引入到BESO法中，多载荷工况的刚度拓扑优化问题可以描述为：

式(9)中，m是载荷工况的总数，ω_s、C_s分别是工况s的预设权重系数和平均柔顺度，并有

由于各载荷工况的位移场相互独立，带权重的目标函数设计变量x_i求偏导数可得：

于是灵敏度计算式可以归结为：

在本发明一实施例中，在等效静态载荷法中，等效静态载荷被定义为能够产生与动态载荷相同位移场的一组载荷集；可以看出，第i个等效静载荷作用下产生的位移响应与第i个时间步上的动态响应完全相等；

已知动态分析控制方程：

式(12)中，M、C、K分别为结构质量、阻尼和刚度矩阵，均为设计变量x的函数；和u(t)分别为加速度、速度和位移向量，f(t)为外部施加的动载荷；

对式(12)式进行变换有：

式(13)中，t＝1,2,3…m；

再用静态变量s替换f_eq(t)中的时间变量t，即可算出等效静载荷：

f_eq(s)＝K(x)u(s),s＝1,2,3...m>

在本发明一实施例中，该改进双向渐进法的等效静载荷法动态响应拓扑优化方法，具体实现如下，

S1：进行有限元建模，设定初始设计变量x(k)＝x(0)，初始设计循环次数k＝0；

S2：对以x(k)作为密度的有限元模型进行动力学分析，求其时间域上的节点位移响应u(t)；

S3：将时间域上的连续响应u(t)离散为m个时间点上的响应u(s)，并计算m个时刻的等效静载荷f_eq(s)，计算公式如下：

f_eq(s)＝K(x)u(s)；s＝1,2,3...m>

S4：将等效静载荷输入静态多目标拓扑优化数学模型中，每一个时刻的等效载荷为一种工况，每个工况赋予ω_k＝1/l的权重因子，进行优化迭代直至收敛；

S5：当k＝0时，令k＝k+1并进入步骤S2；当k≠0，如果满足收敛条件，优化结束；如果不满足收敛条件，依然令k＝k+1并进入步骤S2；至此完成一次设计循环；收敛条件可以表示为：

该方法将渐进结构法替换变密度法作为内循环方法，第一次外循环时，初始设计域为满设计域，此时按照V_k+1＝V_k(1±ER)，其中k＝1,2,3…，计算每次迭代的目标体积，从第二次外循环开始，均满足V_k+1＝V*，即设计域同时包含0、1单元，此时有V_del＝V_add＝AR，这种情况下，渐进结构法的优化速度要比满设计域时快。

相较于现有技术，本发明具有以下有益效果：本发明方法将改进的BESO法引入到等效静载荷优化中，改进了动态响应拓扑优化流程，并建立了新的收敛条件；通过matlab软件对2个算例进行了优化分析，经比较可以看出，本发明提出的方法，能有效减少k大于1次设计循环阶段的优化时间，两个算例中优化所需CPU总时间比原始方法分别减少了21.8％、25.3％。该方法具有计算速度快、优化流程简单的特点；随着外载荷作用点的增多、优化规模的增大，优化所需的设计循环次数也会增加，该方法将会带来更大的优势。

附图说明

图1为本发明静态载荷等效过程图。

图2为本发明改进双向渐进法的等效静载荷法动态响应拓扑优化方法流程图。

图3为算例一。

图4(a)为针对算例一的原始方法优化结果图。

图4(b)为针对算例一的本申请方法优化结果图。

图5为算例二。

图6(a)为针对算例二的原始方法优化结果图。

图6(b)为针对算例二的本申请方法优化结果图。

具体实施方式

下面结合附图，对本发明的技术方案进行具体说明。

本发明的一种改进双向渐进法的等效静载荷法动态响应拓扑优化方法，将改进的BESO方法应用于等效静载荷法中作为优化迭代方法，替换掉原始的变密度法，并建立新的收敛条件，从而改进动态响应拓扑优化流程。

所述改进的BESO方法具体如下，

在对单元处理之前，应事先给出每次迭代的单步目标体积V_k+1；初始体积V_k大于或小于整体的体积约束V*；因此，体积进化迭代式可表示为：

V_k+1＝V_k(1±ER),(k＝1,2,3...)>

式(1)中k为当前迭代次数，ER为体积进化率；

若目标体积V_k+1小于当前体积V_k，则必须将单元按灵敏度大小排列并删除灵敏度低的单元，在删除单元之前，空单元的灵敏度有可能大于待删单元灵敏度，因此要先考虑增添单元，这里引入增删阀值α^th，α^th定义为待删单元灵敏度中的最大值，有：

α_i>α^th>

式(2)中α_i为增添单元的灵敏度；当增添单元的体积大于体积增添率AR时，需要重新计算α^th直到满足增添单元的体积小于或等于AR为止；然后删除灵敏度小于该值的单元；显然，删除单元的体积V_del可以表示为：

V_del＝V_k-V_k+1+V_add>

式(3)中V_add为增添单元的体积；

若目标体积V_k+1大于当前体积V_k，则必须将单元按灵敏度大小排列并增添灵敏度高的单元，在增添单元之前，实单元的灵敏度有可能小于待增单元灵敏度，因此要先考虑删除单元，这里引入增删阀值β^th，β^th定义为待增单元灵敏度中的最小值，有：

β_i＜β^th>

其中β_i为删除单元的灵敏度；当删除单元的体积大于体积增添率AR时，需要重新计算β^th直到满足删除单元的体积小于或等于AR为止；然后增添灵敏度大于该值的单元；显然，增添单元的体积V_add可以表示为：

V_add＝V_k-V_k+1+V_del>

重复进行有限元分析与单元增删的过程，直至迭代循环满足两个终止条件：一是体积约束，即当前进化体积已满足整体约束体积V*；二是收敛准则，考虑目标函数的变化量，收敛判别式规定为：

式(6)中k为当前迭代步数，τ为内部许用收敛因子；M是一整数，用以限制并稳定柔顺度在多次连续迭代中的平均变化量，通常取5；计入局部多个M值，能够使目标函数波动较小，曲线更为平滑，有效防止局部收敛。

以最小化柔顺度为目标函数，以体积为约束条件，以单元本身为设计变量，基于灵敏度分析改进的BESO刚度优化数学表达式可描述为：

式(7)中F和U分别为外载荷向量和结点位移向量，C为结构柔顺度，V_i是个体单元的体积，V*为体积约束，N为设计域中单元总数量；x_i为双向赋值的非连续设计变量，即材料单元的密度；K为结构总体刚度矩阵，x_min为材料单元的密度的一个小值。

所述改进的BESO方法中，删除单元的方法为软杀法，即将灵敏度低的待删单元的密度人为减至一个小值x_min，以防止刚度矩阵奇异，可令x_min＝0.001，但后续仍参与有限元分析，而不是直接删除；因此，采用单元软杀法的BESO方法的单元灵敏度计算式可归结为：

式(8)中K_i⁰表示实单元的单元刚度矩阵，u_i为单元节点位移向量,P为惩罚因子，可令P＝3。

将权值法引入到BESO法中，多载荷工况的刚度拓扑优化问题可以描述为：

式(9)中，m是载荷工况的总数，ω_s、C_s分别是工况s的预设权重系数和平均柔顺度，并有

由于各载荷工况的位移场相互独立，带权重的目标函数设计变量x_i求偏导数可得：

于是灵敏度计算式可以归结为：

如图1所示，在等效静态载荷法中，等效静态载荷被定义为能够产生与动态载荷相同位移场的一组载荷集；可以看出，第i个等效静载荷作用下产生的位移响应与第i个时间步上的动态响应完全相等；

已知动态分析控制方程：

式(12)中，M、C、K分别为结构质量、阻尼和刚度矩阵，均为设计变量x的函数；和u(t)分别为加速度、速度和位移向量，f(t)为外部施加的动载荷；

对式(12)式进行变换有：

式(13)中，t＝1,2,3…m；

再用静态变量s替换f_eq(t)中的时间变量t，即可算出等效静载荷：

f_eq(s)＝K(x)u(s),s＝1,2,3...m>

如图2所示，本发明的改进双向渐进法的等效静载荷法动态响应拓扑优化方法，具体实现如下，

S1：进行有限元建模，设定初始设计变量x(k)＝x(0)，初始设计循环次数k＝0；

S2：对以x(k)作为密度的有限元模型进行动力学分析，求其时间域上的节点位移响应u(t)；

S3：将时间域上的连续响应u(t)离散为m个时间点上的响应u(s)，并计算m个时刻的等效静载荷f_eq(s)，计算公式如下：

f_eq(s)＝K(x)u(s)；s＝1,2,3...m>

S4：将等效静载荷输入静态多目标拓扑优化数学模型中，每一个时刻的等效载荷为一种工况，每个工况赋予ω_k＝1/l的权重因子，进行优化迭代直至收敛；

S5：当k＝0时，令k＝k+1并进入步骤S2；当k≠0，如果满足收敛条件，优化结束；如果不满足收敛条件，依然令k＝k+1并进入步骤S2；至此完成一次设计循环；收敛条件可以表示为：

该方法将渐进结构法替换变密度法作为内循环方法，第一次外循环时，初始设计域为满设计域，此时按照V_k+1＝V_k(1±ER)，其中k＝1,2,3…，计算每次迭代的目标体积，从第二次外循环开始，均满足V_k+1＝V*，即设计域同时包含0、1单元，此时有V_del＝V_add＝AR，这种情况下，渐进结构法的优化速度要比满设计域时快。

传统变密度法如果直接对含有中间密度单元的模型进行动力学分析会引起的网格畸变^[22]，因此需要在下一次外循环开始前运用EEM(eliminate>

下面根据原始等效静载荷法和本发明提出的方法计算两个算例。

单个动载荷条件下的优化：

如图3，设计域为160mm*100mm的二维悬臂梁结构，将其划分为40*25个单元的有限元模型。模型右边中点受到一个周期的正弦动载荷f1＝1000sin(5πt)[N]作用。每个时间点的载荷工况权值为ω_s＝1/l。算例其他优化参数如表1所示。

表1拓扑优化参数

两种方法的优化结果如图4所示。

结果数据对比如表2所示。

表2结果对比

从优化结果图4(a)和图4(b)可看出，对于受到单个动载荷的悬臂梁结构，本发明研究所提出的方法得到的优化结果与原始方法相似，并且目标函数值也非常接近，可以认为该方法是可行的。同时本发明研究提出的方法在第二次设计循环阶段的优化时间为28.08s，明显快于原始方法在第二次设计循环阶段的56.20s，两种方法的最终优化时间分别为90.48s、115.74s。因此从整体的优化效率上看，本发明提出的方法要高于原始方法。

多个动载荷条件下的优化：

如图5，设计域为200mm*100mm的二维悬臂梁结构，将其划分为50*25个单元的有限元模型。模型右上方受到一个周期的正弦动载荷f1＝10000sin(5πt)[N]作用，下方中点受到一个周期的余弦动载荷f2＝20000cos(5πt)[N]作用。

每个时间点的载荷工况权值为ω_s＝1/l。算例其他优化参数如表3所示。

表3拓扑优化参数

两种方法的优化结果如图6所示。

结果数据对比如表4所示。

表4结果对比

从优化结果图6(a)和图6(b)可看出，对于受到两个动载荷的悬臂梁结构，得到的结论与算例1类似，两种方法下优化结果与目标函数非常接近，可以认为BESO法作为优化迭代方法可行。同时随着载荷作用点数的增加，设计循环次数增加到3次，由表4的优化结果中可以看出，本发明提出的方法在第二、三次设计循环阶段的优化时间分别为43.14s和40.60s,均明显快于原始方法在第二、三次设计循环阶段的71.14s和75.87s，两种方法的最终优化时间分别为171.27s、229.16s。因此从整体的优化效率上看，本发明提出的方法要高于原始方法。

参考文献：

[1]BENDSOE M P,SIGMUND O.Topology optimization:Theory,Methods,and Applications[M].New York:Methods and Applications,2003.

[2]BENDSOE M P,KIKUCHI N.Generating optimal t-opologies in structural design using a homogenization method[J].Computer Methods in Applied Mechanics and Engineering,1988,71(2):197-224.

[3]BENDSOE M P.Optimal Shape Design as a Materi-al Distribution Problem[J].Structural Optimization,1989,1(4):193-202.

[4]Xie Y M,and StevenG P.A Simple Evolutionary Procedure for Structural Optimization[J].Computers and Structures,1993,49(5):885-896.

[5]Querin O M,Young V,Steven G P,Xie Y M.Computational efficiency and validation of bi-directional evolutionary structural optimization[J].Comput.Meth.Appl.Engng,2000(189):559-573.

[6]Rozvany G I.A critical review of established metho-ds of structural topology optimization[J].Structural and Multidisciplinary Optimization,2009,37(3):217-237.

[7]Huang X,Xie M.Evolutionary topology optimization of continuum structures:methods and applications[M].John Wiley&Sons,2010.

[8]Choi S H,Kim S R,Park J Y,Han S Y.Multi-Objective Optimization Of The Inner Reinforcement For A Vehicle’S Hood Considering Static Stiffness And Natural Frequency[J].International Journal Of Automotive Technology,2007,8(3):337-342.

[9]Choi K H,Park J Y,Ryu S P,Han S Y.Reliabili ty-Based Topology Optimization Based On Bidirectional Evolutionary Structural Optimization Using Multi-Objective Sensitivity Numbers[J].International Journal Of Automotive Technology,2011,12(6):849-856.

[10]Afimiwala KA,Mayne RW.Optimal Design of an Impact Absorber[J].Journal of Engineering for Industry,Transactions of the ASME 96:124-130.

[11]Fox RL,Kapoor MP.Structural Optimization in the Dynamic Regime[J].A Computational Approach.AIAA Journal 8:1798-1804.

[12]Choi W S,Park G J.Structural Optimization Using Equivalent Static Loads at All the Time Intervals[J].Computer Methods in Applied Mechanics and Engin-eering,2002,191(19-20):2105-2122.

[13]PARK G J,KANG B S.Validation of a Structural Optimization Algorithm Transforming Dynamic Loadsinto Equivalent Static Loads[J].Journal Of Optimizat-ion Theory And Applications,2003,118(1):191-200.

[14]Jang H H,Lee H A,Park G J.Preliminary Study on Linear Dynamic Response Topology OptimizationUsing Equivalent Static Loads[J].Transactions of theKorean Society of Mechanical Engineers,2009,33(12):1357-1493.

[15]Lee H A,Park G J.Nonlinear dynamic response top-ology optimization using the equivalent static loads method[J].Computer Methods in Applied Mechanics and Engineering,2015,283(1):956-970.

[16]苏义脑,唐雪平.初弯曲纵横弯曲梁的等效载荷法及其应用[J].力学与实践,20040.

Su Yinao,Tang Xueping.The equivalent load method and its application to the bending and transverse beam of the initial bending beam[J].Mechanics and Practice,20040.

[17]黄武龙.基于等效静态载荷方法的大型复杂结构的轻量化设计[D].广州:广东工业大学,2013(5)17.

Huang Wulong.Based on lightweight design of equivalent static load method of large complex structures[D].Guangzhou:Guangdong University of technology,2013(5)17.

[18]赵礼辉.ESL法在汽车结构优化设计中的应用[D].上海:上海交通大学,2009(2).

Zhao Lihui.Application of ESL method in automobile structure optimization design[D].Shanghai:Shanghai Jiao Tong University,2009(2).

[19]杨志军.基于等效静态载荷原理的高速机构结构拓扑优化方法[J].机械工程学报,2011,9(47):119-126.

Yang Zhijun.Topology optimization method of high speed mechanism based on equivalent static load principle[J].Journal of mechanical engineering,2011,9(47):119-126.

[20]芮强,王红岩.基于等效静态载荷法的结构动态优化[J].汽车工程,2014(1).

Rui Qiang,Wang Hongyan.Dynamic optimization of structure based on equivalent static load method[J].Automotive Engineering,2014(1).

[21]Kim E,Kim H,Baek S,Cho M.Effective structural o-ptimization based on equivalent static loads combined with system reduction method[J].Structural and Mul-tidisciplinary Optimization,2014,50:775-786.

[22]Lee J Y,Park G J.Dynamic Response Topology Optimization in the Time Domain Using Equivalent Static Loads[J].AIAA Journal,2012,50(1):226-234.。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。