史峰源，李世强，刘志芳

（1. 太原理工大学 机械与运载工程学院，应用力学研究所, 山西, 太原 030024；2. 太原理工大学 山西省材料强度与结构冲击重点实验室, 山西, 太原 030024）

拓扑优化的目的是运用结构拓扑信息寻求结构材料在给定设计域内的最佳材料分布，从而满足一个或多个设计约束的同时最大化或最小化某些设计目标[1]. 拓扑优化是研究结构最优拓扑形式的优化问题，在优化的过程中去除多余或者影响小的材料从而构建满足约束条件的最优结构. 目前，比较成熟的连续体拓扑优化方法主要有均匀化方法[2−3]（homogenization method）、变密度法[4−5]（solid isotropic material penalization，SIMP）、水平集方法[6−7]（level set method）以及渐进结构优化方法[7−8]（evolutionary structural optimization，ESO）等. 渐进结构优化方法作为拓扑优化的一种重要方法，最早是由XIE 等[8]基于结构通过逐渐删除低应力材料而使结构进化到最优的概念的基础上提出的. 渐进结构优化方法不会产生中间密度，优化过程中会得到清晰的拓扑构型，但渐进结构优化方法只能不断删除低应力单元，QUERIN等[9−10]提出了双向渐进结构拓扑优化方法（BESO），不但可以逐渐删除低应力单元，同时也可以逐渐添加已经被删除的单元. 这样既能删减单元同时又可以添加单元，解决了最初只能删减单元的问题.HUANG 等[11]提出了解的收敛性与网格无关的一种BESO 方法，这种方法结合了灵敏度滤波和利用历史信息稳定优化方案，使优化过程更加高效平稳，在结构设计问题中表现出了良好的性能. 目前，双向渐进结构拓扑优化方法已经显示出高效和稳健的优化性能，成为一种无论是在学术研究还是工程应用中都广泛使用的设计方法[12].

从SIGMUND[13]开始将杨氏模量、导电率或导热率不同的两种材料同时引入设计域，在某些情况下用以显著提升结构特性，从此多相材料拓扑优化设计得到广大学者的关注[14]，复合材料由于多种材料特性之间的相互补充，从而在性能方面要比传统材料更加优良，具有比强度高、比模量大、可设计性等优点，并且能满足结构的一些特殊功能，如传热性[15]、导电性等.

目前拓扑优化方法多用于结构静态载荷优化，而工程实际中结构常受到动态载荷[16]，直接对结构进行非线性动力响应优化分析虽然在理论上可行，但是计算过程异常复杂、计算效率低，而且计算难以收敛，无法得到较为理想的拓扑构型[17]. 为了可以高效地进行结构非线性动态优化设计，PARK 等[18]提出了等效静力载荷优化方法，通过构造等效静力载荷的方法将动力优化设计问题转化成为静力优化设计问题[19]，并且验证了方法的可行性，提高了非线性动力优化的设计效率，这种方法也因此得到了广泛应用[20−21].

本文针对结构非线性动态响应拓扑优化优化流程复杂，计算效率低等问题，将双向渐进结构拓扑优化方法和等效静力载荷优化方法相结合，对非线性动态响应拓扑优化流程进行改进，采用内外双层循环的方式，外层循环实现对结构的非线性动态响应分析，依据结构的非线性动态响应构造等效静态载荷，内层循环依据构造的等效静态载荷，并将权值法引入其中，采用改进的双向渐进结构优化方法进行优化，简化优化流程，优化过程更加平稳，整体优化时间大幅减小，计算效率大幅提高. 采用双向渐进结构优化方法能够避免中间密度的产生，提高了工程可制备性，同时也将材料插值模型引入其中，实现冲击载荷下双相复合材料拓扑优化设计，在保证结构质量相同的情况下使结构具有更高的刚度. 最后通过两个算例来验证改进方法的可行性和高效性，在得到优化结果后对最终优化结构进行对比和动态响应分析，说明本文提出的方法更加适用于动载荷下结构拓扑优化设计.

1 等效静力载荷优化法

如果各个时刻下的结构响应通过积分计算的方法来获取，这会使得结构分析变得十分复杂并且计算效率低. 因此等效静力载荷优化方法提出了一种通过构造等效静力载荷可以高效的获得非线性结构响应的等效分析方法，这种方法要求线弹性静力分析在关心时刻下要与非线性动力分析具有相同的结构响应.

图1 为等效静力载荷的等效过程，在相对应的时间步中，等效静力载荷是能够与结构施加动载荷产生相同位移场的静载荷. 在等效时，将结构的非线性动力学分析按所需的时间步分开，每一个时间步的动态响应等效为静力载荷的一个工况，并且要求相对应时间步结构在动载荷下产生的位移和相对应等效静态载荷产生的位移相等[19].

图1 等效静态载荷的等效过程Fig. 1 Equivalent process of equivalent static load

非线性动力问题的控制方程为

2 基于等效静力载荷优化方法的双向渐进结构拓扑优化

2.1 多相材料的双向渐进结构优化方法

在进行单相材料双向渐进结构拓扑优化时，材料的弹性模量可以作为单元密度变量的函数进行插值：

在双向渐进结构优化方法中，可以使用离散设计变量进行结构优化，并且在这种方法中利用单个元素的灵敏度数值的相对排序来进行优化，由于采用“软杀”方法，为了防止刚度矩阵奇异的发生，用较小的相对密度来表示灵敏度数值低的待删除单元的相对密度，后续该单元仍参与有限元分析而不是直接删除该单元，同样在后续的结构优化中该单元如果具有较高应力时也可以恢复该单元.

由上式可以看出当有多种材料参与优化设计时，灵敏度数值完全不同于使用单元应变能为标准的原始的双向渐进结构优化方法.

2.2 多工况的双向渐进结构优化方法

原始的双向渐进结构优化方法适用于单一静态载荷工况，为了进行动态响应优化，每一个关心时刻的动态响应等效为等效静力载荷的一个工况，则需要在原始的双向渐进结构优化方法中引入权值法，即有多个载荷工况，每一个载荷工况赋予相应的权重，并且所有载荷工况下的权重系数相加为1，多工况的双向渐进拓扑优化“软杀”方法的优化问题可以表示为：

2.3 基于等效静力载荷优化方法的复合材料渐进结构拓扑优化

将等效静力载荷优化方法和双向渐进结构优化方法相结合，并将权值法和多相材料设计方法引入其中，从而高效解决非线性复合材料动态优化设计问题，由式(6)和式(10)可以得出基于等效静力载荷优化方法的复合材料双向渐进结构优化问题描述为

由式(9)和式(11)可以得出基于等效静力载荷优化方法的复合材料双向渐进结构优化单元灵敏度的计算公式为

为了保证所关心的各个时刻构造的等效静力载荷使结构在线弹性静力分析中的结构响应与初始动载荷结构在非线性动力分析中的结构响应相等，优化设计流程采用内外双层迭代方式进行结构的拓扑优化设计，其中，外层循环实现对结构的非线性动态响应分析，依据结构的非线性动态响应构造等效静态载荷，内层循环依据构造的等效静态载荷，并将权值法引入其中，采用改进的双向渐进结构优化方法进行优化.

3 基于等效静力载荷优化方法的双向渐进结构拓扑优化流程

优化流程如图2 所示，具体如下：

图2 优化流程Fig. 2 Optimized process

① 初始化设计域，使用有限元网格对整个设计域进行离散，设定设计变量初始值，定义目标体积分数、进化率、惩罚指数等参数.

② 基于当前设计进行动载荷下非线性动力响应分析并提取出各个关心时刻的结构响应. 本文在Abaqus 中进行结构的非线性动力响应分析并且在Matlab 中提取出各个关心时刻的结构响应.

③ 基于②提取的各个关心时刻的结构响应构造等效静力载荷. 由于关心优化设计多个时刻下的结构响应，所以需要构造多个等效静力载荷.

④ 基于③构造的等效静力载荷运用改进的多工况双向渐进结构优化方法进行拓扑优化设计(式(12)). 此时内层迭代由于考虑多个等效静力载荷，将每一个关心时刻构造的等效静力载荷作为一个工况，引入权值法来解决多工况静力优化设计问题.

⑤ 确定下一次迭代的目标体积，若当前体积Vi大于目标体积时，则由公式Vi+1=Vi(1−ER)计算下一次迭代的体积，若当前体积Vi小于或等于目标体积时，则令目标体积为下一次迭代的体积.

⑥ 更新设计变量，将所有单元按照计算的单元灵敏度进行排序，由目标材料体积和灵敏度排序可以得到灵敏度阈值 αth，当单元灵敏度 αi≤αth时，设计变量由1 转换为xmin， 当单元灵敏度 αi>αth时，设计变量由xmin转换为1.

⑦ 循环迭代②～⑥，直到所有材料当前体积达到目标体积，并且得到的优化结果满足收敛条件，在目标函数的变化中定义以下收敛准则：

当所有材料体积达到目标体积后，运用式(14)判断得到的优化结果是否满足收敛条件，如果满足收敛条件，优化结束，输出最终优化结果. 如果不满足收敛条件则进入第二步，进行当前结构的动态响应分析.

4 冲击载荷下结构优化与动态响应分析

4.1 简单冲击载荷下结构优化与动态响应分析

如图3 所示，设计域为长240 mm，高40 mm 的二维固支梁结构，在其上边中点处施加一个正弦冲击载荷F=3 000sin(5πt)的作用力，作用时间为0.2 ms，将其离散为240×40 个4 节点四边形单元的有限元模型.

图3 固支梁设计区域Fig. 3 Clamped beam design area

此算例分别进行4 种工况下的优化：

工况1：作为对比进行静态载荷下单相材料优化，静态载荷取脉冲载荷峰值F=3 000 N，材料特性为弹性模量E=200 GPa ，泊松比 µ=0.3；

工况2：冲击载荷下单相材料优化，材料特性为弹性模量E=200 GPa ，泊松比 µ=0.3；

工况3：冲击载荷下复合材料优化，两种材料弹性模量分别为E1=200 GPa，E2=20 GPa，泊松比µ=0.3；

工况4：冲击载荷下复合材料优化，两种材料弹性模量分别为E1=200 GPa，E2=40 GPa，泊松比µ=0.3.

其余优化参数均为：体积约束V=0.5，过滤半径rmin=3 mm， 惩罚指数p=3.0， 体积进化率ER=0.02.

最终的优化结果如图4 所示，从最终的优化结果可以看出，单相材料进行优化时，由于只有一种材料承受所施加的载荷，结构各个部分没有双相复合材料的优化结果光滑. 复合材料优化时由于有软材料的加入，使得优化结果中硬材料结构具有一定的弧形，且随着优化时结构中软材料弹性模量的增加，更加能独立承担所受载荷时，宏观上相应的硬材料结构中的支撑会变少且硬材料结构会更加平滑，如图4（c）、4（d）对比所示.

图4 优化结果Fig. 4 Results of optimization

原始的“软杀”双向渐进结构拓扑优化方法在进行结构优化时，为了防止刚度矩阵奇异的发生，用较小的相对密度来表示灵敏度数值低的待删除单元的相对密度，后续该单元仍参与有限元分析而不是直接删除该单元，优化时会赋予这些待删除单元一个很低的弹性模量，但这并不意味着这种方法也能够进行复合材料结构优化，采用原始的“软杀”方法进行复合材料优化时，设计域内的各种材料分布无法连接成结构从而所得到的柔顺度或者应变能会十分异常，当有多种材料参与优化设计时，灵敏度数值完全不同于使用单元应变能为标准的原始双向渐进结构优化方法. 所以本文在内层循环中引入复合材料优化方法(式（12）（13）)，在内层循环迭代中可以合理、平稳且高效地进行复合材料优化，如图6、图7所示.

图5、图6 为单相材料(E=200 GPa)与双相复合材料(E1=200 GPa、E2=20 GPa)优化历程图，可以看出，单相材料在进行结构优化时，从结构优化历程来看，目标函数没有复合材料优化时平稳，主要原因是因为在优化过程中，一些重要的连接杆件被打断，从而造成结构传力路径有所改变，结构的应力分布也会随之改变，结构会重新分配单元，高应力的地方添加单元，低应力的地方删减单元，所以优化历程会存在一些起伏；而复合材料在进行结构优化时，由于有软硬两种材料在结构中起承载作用，当硬材料结构由于优化过程中体积分数逐渐减小而造成重要的硬材料连接杆件被打断时，软材料单元会立即填充硬材料单元删除部分，结构传力路径和结构应力分布不会发生间断性突变，所以优化历程会很平稳，并且更容易收敛，优化迭代次数和优化时间会减少，降低时间成本. 算例中单相材料优化(如图5)进行了43 次迭代，复合材料优化(如图6)进行了41 次迭代，由于算例设计域选取的是二维固支梁，所以单相材料优化和复合材料优化迭代次数相差不多，当设计域和边界条件更加复杂时，复合材料优化在优化历程平稳程度和优化效率等方面相比于单相材料优化优势会更加明显. 所以与单相材料优化相比，双相复合材料优化既能有效降低结构柔顺度，又能平稳且高效地进行结构优化.

图5 冲击载荷单相材料优化历程，E=200 GPaFig. 5 Evolutionary histories of single phase materials under impact load,E=200 GPa

图6 冲击载荷复合材料优化历程，E 1=200 GPa ，E2=20 GPaFig. 6 Evolutionary histories of composite materials under impact load,E1=200 GPa ,E2=20 GPa

图4（a）为直接运用双向渐进结构优化方法对设计域施加恒定的静态载荷进行单相材料静态优化得到的结构，图4（b）和4（c）分别为运用本文改进的冲击载荷下结构优化方法对设计域施加正弦脉冲载荷进行单相材料和双相复合材料动态优化得到的结构，为了对比这些结构在冲击载荷下的动力响应特性，将3 种结果与同质量同体积下的实心结构在Abaqus中施加相同的正弦脉冲载荷进行动态响应分析，得到位移-时间曲线如图8 所示.

图8 动态响应位移-时间曲线Fig. 8 Dynamic response displacement-time curve

由图8 可以看出，运用本文改进的优化方法得到的结构要比经典的静载荷优化的结构振幅小，4 种结构向上的最大位移分别为0.33，0.37，0.79 和0.23 mm，向下的最大位移分别为0.34，0.41，0.81 和0.24 mm. 以单相材料优化为例，如图9 所示，在相同正弦脉冲载荷作用下，利用改进的优化方法优化的单相材料结构向上（向下）的最大位移是经典静力学优化结果的89.18%（82.92%），是同体积同质量的实心结构的41.77%（41.97%）. 说明运用本文动态载荷优化方法优化的结构在动态载荷下承载能力更好，优化方法合理高效.

由图8 和图9 也可以看出，优化后的双相复合材料/结构振幅最小，说明双相复合材料结构比单相材料结构在动态载荷下具有更好的承载能力.

图9 动态响应结果对比Fig. 9 Comparison of dynamic response results

4.2 多个冲击载荷下结构优化与动态响应分析

如图10 所示，设计域为长80 mm，宽40 mm 的二维悬臂梁结构，在其上边中点处施加一个三角脉冲载荷F1=−20 000t+2 000(N)的作用，同时在其右边最下端施加一个三角脉冲载荷F2=−5 000t+1 000(N)的作用，作用时间均为0.2 ms，将设计域离散为80×40 个4 节点四边形单元的有限元模型.

图10 悬臂梁设计区域Fig. 10 Cantilever beam design area

此算例分别进行3 种工况下的优化：

工况1：作为对比进行静态载荷下单相材料优化，静态载荷分别取三角脉冲载荷峰值F1=2 000 N，F2=1 000 N， 材料特性为弹性模量E=200 GPa，泊松比 µ=0.3， 密度 ρ=7 800 kg/m3；

工况2：动载荷下单相材料优化，材料特性为弹性模量E=200 GPa ，泊松比 µ=0.3， 密度 ρ=7 800 kg/m3；

工况1 和工况2 体积约束均为V=0.5，工况3 为保证与前面两种工况结构质量相等，进而与前两种优化进行对比，此时硬材料体积V1=0.49，软材料体积V2=0.51. 其余优化参数均为：过滤半径rmin=3 mm，惩罚指数p=3.0， 体积进化率ER=0.02.

最终的优化结果如图11 所示.

图11 优化结果Fig. 11 Results of optimization

图12、图13 分别为冲击载荷下单相材料和复合材料结构优化历程图，可以看出优化历程都很平稳.说明本文提出的冲击载荷下结构的优化方法无论是单相材料优化还是复合材料优化均平稳高效，同时在保证质量相等的前提下，双相复合材料优化所得结构的柔顺度要比单相材料小. 复合材料结构中由于存在两种材料，在优化平稳程度和结构具有更小的柔顺度等方面会更加具有一些优势.

图12 冲击载荷下单相材料优化历程,E=200 GPaFig. 12 Evolutionary histories of single phase materials under impact load,E=200 GPa

图13 冲击载荷复合材料优化历程，E1=200 GPa，E2=1.97 GPaFig. 13 Evolutionary histories of composite materials under impact load,E1=200 GPa ,E2=1.97 GPa

图11（a）为直接运用经典的双向渐进结构优化方法对设计域施加恒定的多个静态载荷进行单相材料优化得到的结构，图11（b）为运用本文改进的优化方法对设计域施加两个冲击载荷进行单相材料动态优化得到的结构，对比结果能够看出最终得到的优化结构区别明显. 为了对比这些结构在冲击载荷下的动力响应特性，将两种结果和同质量同体积下的实心结构在Abaqus 中施加两个脉冲载荷进行动态响应分析，得到位移-时间曲线如图14 所示.

图14 动态响应位移-时间曲线Fig. 14 Dynamic response displacement-time curve

由图14 可以看出，运用本文改进的优化方法获得的结构要比另外两种结构振幅小，3 种结构向下的最大位移分别为0.22，0.85 和0.73 mm，向上的最大位移分别为0.29，0.91 和0.73 mm，所以在2 个冲击载荷作用下，运用改进后的方法优化的结构向下（向上）的最大位移是经典方法优化结构的25.88%（31.87%），是同体积同质量的实心结构的30.14%（39.73%），如图15 所示.

图15 动态响应结果对比Fig. 15 Comparison of dynamic response results

由图15 也可以看出，双相复合材料结构比单相材料结构在受到相同两个冲击载荷作用时振幅小，说明双相复合材料结构比单相材料结构在受到多个冲击载荷作用时具有更好的承载能力.

对比两个算例的优化结果，在单个动态载荷作用下，运用本文提出的优化方法所得结构和经典静态载荷优化方法所得的结构较为相似(如图4(a)和4(b))，而在更加复杂冲击载荷作用下（如：同时作用两个冲击载荷），运用本文提出优化方法的优化结构与经典静态载荷优化方法的优化结构有明显区别(如图11(a)和11(b))，并且两种结构在同时施加多个脉冲载荷作用时，运用改进方法优化所得结构的承载能力远好于其他结构的承载能力. 值得注意的是，在受到多个冲击载荷作用时，运用经典静态载荷优化的结构比同质量同体积实心结构的承载能力差(如图14、图15)，这说明利用经典的静态载荷优化方法所得到的结构不能适应复杂的冲击载荷作用，必须构建相应的冲击载荷下结构的优化流程和方法.

本文为了说明所提出方法的有效性和高效性，进行了两个算例的优化加以对比说明，为了更加直观清晰，文中选取的动态载荷均能以一定的函数形式表示. 若结构所受到的动态载荷呈现无规律的情况，通过该方法依然可以对结构进行拓扑优化. 只需要保证在等效时，将结构的非线性动力学分析按所关心的时间步分开，每一个时间步的动态响应等效为静力载荷的一个工况，并且要求在相对应的时间步，结构在动载荷下产生的位移和相对应等效静态载荷产生的位移相等，即可将无规律的动态载荷等效为静力载荷.

5 结 论

本文对冲击载荷下结构拓扑优化流程和方法进行改进，将双向渐进结构拓扑优化方法和等效静力载荷优化方法相结合，采用内外双层循环的方式，外层循环实现对结构的非线性动态响应分析，依据结构的非线性动态响应构造等效静态载荷，内层循环依据构造的等效静态载荷，将权值法引入其中，采用多工况下的双向渐进结构优化方法进行优化，简化优化流程. 同时又将材料插值模型引入其中，实现双相复合材料/结构在冲击载荷作用下的结构优化设计. 通过两个算例可以看出，本文建立的方法在进行冲击载荷加载下的单相材料和双相复合材料/结构拓扑优化设计时，整体优化历程十分平稳，结构清晰合理，计算简洁高效.

将优化后的结构在Abaqus 中进行动态响应对比分析，结果表明运用本文建立的优化方法获得的结构在受到冲击载荷作用时具有更好的承载能力，随着动载荷工况复杂程度的提高，这种方法具有更明显的优势，优化后的结构在冲击载荷作用下承载能力更好. 通过算例与经典的静态载荷优化方法进行对比，不但证实所提出方法的有效性和高效性，更加说明利用经典的静态载荷优化方法所得到的结构不能适应复杂的冲击载荷作用，所以构建冲击载荷下结构的优化流程和方法，对于结构设计具有重要的理论和工程意义.