梁宽,付莉莉,张晓鹏,*,罗阳军

1. 大连理工大学 工业装备结构分析国家重点实验室, 大连 116024 2. 西安航天动力研究所,西安 710100

振动在航空航天薄壁结构中十分普遍,然而剧烈的振动将导致结构性能下降并带来安全隐患。通过优化技术对结构中材料与设备的合理布局以改善提高振动系统的动力学性能具有重要的意义。在各类优化方法中,结构动力学拓扑优化被认为是改善结构动力学设计的最有效途径,近年来受到学者的广泛关注。目前主流的结构拓扑优化方法,如变密度法、水平集法、渐进结构方法等,均已应用于结构动力学拓扑优化问题设计中。

相比于结构静力学拓扑优化问题,结构动力学拓扑优化问题面临着因局部最优解陷阱而导致的收敛困难、灵敏度分析复杂、局部模态等诸多难点,限制了结构动力学拓扑优化方法的发展和应用。针对上述困难,部分学者开展了长期的研究工作。例如,Du和Olhoff提出了广义增量频率(GIF)拓扑优化方法,从而克服了动力学优化中结构高频振动多峰值带来的收敛性和局部最优解难题。

结构动力学优化问题的研究目标主要分为3大类,包括：结构固有频率优化、结构稳态响应优化、结构瞬态动力学优化。经过多年的研究,固有频率的结构优化方法获得了长足的发展,并在结构重频率等难题中取得了重要的研究成果。近年来,学者们更为关注结构稳态与瞬态动力荷载下结构动力学拓扑优化的研究。在结构稳态动力学拓扑优化研究中,Rong等研究了能够有效抑制随机激励下动力学局部模态的结构拓扑优化方法；Zhu等研究了受到基础简谐振动的结构动力学拓扑优化方法,郑玲等则研究了基于渐进优化算法的车身结构减振降噪阻尼材料最优拓扑优化布局方法。在结构瞬态动力学优化中,Zhang和Kang研究了瞬态动力载荷下基于主动振动的压电智能结构拓扑优化方法；Yun和Youn研究了瞬态动力载荷下结构黏弹性阻尼材料最优拓扑布局方法。

已有研究虽然解决了部分结构稳态与瞬态动力学优化中的难题,但因结构动力学优化问题中的灵敏度分析与局部最优解陷阱等困难,使得结构动力学响应能够处理的目标函数类型有限,限制了结构动力学拓扑优化问题的应用。例如,在稳态工程结构动力学拓扑优化中,往往希望结构的动力学响应在整个工作频率区间内都具有较小的结构动力学响应。基于梯度的结构动力学拓扑优化方法,需要引入频率区间积分或者包络函数以实现目标函数的光滑化,从而可采用基于梯度的结构拓扑优化方法进行优化问题的求解。然而简谐激励下结构指定位置在不同频率下振幅变化往往十分剧烈,并具有“多峰值”特性,使得包络函数等呈现病态特性,导致优化效果不理想甚至是不收敛。采用基于非梯度优化算法的结构拓扑优化方法进行稳态与瞬态动力学载荷下结构拓扑优化是克服这些困难的可能途径,目前国内外尚未见相关研究。

最近提出的材料场级数展开(Material-Field Series-Expansion,MFSE)方法通过连续体拓扑表征和降维映射能够大幅度减少拓扑优化问题中设计变量个数,从而能够有效解决传统拓扑优化设计变量过多的难题,并结合序列Kriging代理模型算法,能够高效实现非梯度拓扑优化问题求解。这一优化策略具有不需要梯度信息、适合复杂非线性问题、具备一定的全局搜索能力、可充分利用并行计算等优势,十分适合求解复杂的结构稳态与瞬态动力学拓扑优化问题。

本文基于文献[21]中的结构拓扑表征和非梯度拓扑优化方法,给出了稳态载荷下结构指定频率区间内响应以及瞬态载荷下指定关心时间区间的结构动力学非梯度拓扑优化列式。采用有限元方法对结构的稳态与瞬态动力学响应进行求解,并基于材料场级数展开(MFSE)和序列Kriging代理模型优化算法,实现复杂动力学响应目标下结构拓扑优化问题的求解。

1 结构动力学响应优化问题描述

在外加动力载荷作用下,基于有限元离散的结构振动方程可以表示为

(1)

式中：、和为维矩阵,分别表示结构的总质量阵、总阻尼阵和总刚度阵,为结构总自由度个数;为结构瞬态位移响应;()为外加动载荷。

若考虑结构承受瞬态荷载,求解其瞬态动力学响应的一般方法是,将方程(1)在时间域上进行离散,常见的求解方法为Newmark法、Wilson-法等。通过逐步积分,可计算出任意指定时刻的结构振动响应()。

当结构受到简谐激励()=ei(表示载荷向量幅值,为外激励圆频率)的作用,其稳态位移响应可表示为()=ei(表示结构振幅向量)。因此,结构稳态动力学方程可由式(1)可得到如下：

(-+iθ+)=

(2)

基于上述方程,提出适当的目标函数和约束函数,选取精确高效的结构拓扑描述模式,即可建立结构动力学拓扑优化模型。

对稳态结构动力学响应优化问题而言,结构动力学拓扑优化常选择结构动柔度为目标函数,结构动柔度(=||)可以表征结构整体的动响应特征。为了抑制结构稳态动响应多峰值而带来的数值困难,本文选择动柔度对数为目标函数,其表达如下：

(3)

对结构瞬态动力响应优化而言,其目标函数往往选取关心位置在指定时间段内的动力学响应的积分作为目标函数,考虑到振动瞬态响应会随时间产生正负号交替,常取振幅响应的平方(或绝对值)在时间上的积分。因此,本文选取结构指定位置瞬态的动力响应为目标函数：

(4)

结构动力学优化设计往往需要在满足指定体积约束下实现结构动响应最优设计。基于上述目标函数,结构稳态与瞬态动力学拓扑优化问题列式可以表达为

(5)

式中：和分别表征最优设计允许的材料用量和设计域满布材料时材料用量。需要指出的是,式(5)中表示材料的拓扑布局,利用不同的拓扑描述其表达方式有所不同。例如,变密度法往往选取单元或节点密度描述材料场拓扑布局,而水平集方法则利用水平集函数描述材料分布。本文将采用材料场级数展开策略对材料布局进行表征和降维,并在第2节中详述。

2 基于材料场级数展开策略的结构非梯度拓扑优化方法

2.1 基于材料场级数展开策略的结构拓扑表征与降维映射

基于非梯度优化算法的复杂结构拓扑优化问题求解的必要前提是用少量独立的拓扑设计变量实现复杂结构拓扑的描述,本文采用基于材料场级数展开(MFSE)策略对结构的拓扑布局进行描述,则实相材料和空相材料的分布可以通过设计域内观察位置的材料场函数()∈[-1, 1]来表示

(6)

与传统的基于单元密度的结构拓扑描述不同,在MFSE描述框架下相邻2点的材料分布一般不是完全独立的(这也符合结构中材料分布的一般规律),整个材料场的空间相关性仅依赖于观察点的空间相对位置。因此,需要引入空间中任意2点和的相关函数(,)来描述材料场的空间相关性,这一相关函数可以取为

(7)

在实际模拟分析中,需要将无限维材料场离散为观察点,通常在结构设计领域内均匀分布在个观察点上,因此相关矩阵可以构造为

=

(8)

因此材料场可以进一步表示为

(9)

式中：()=[(,)(,)…(,)];和是相关矩阵的第阶特征值和特征向量;(=1,2,…,)为其中的待定系数,可以理解为广义的设计变量,然而对于拓扑优化设计问题而言,观察点常分布于结构单元中心(或节点)处,其数量常达到数千甚至数万,这一量级的设计变量对于基于非梯度的优化拓扑优化设计问题是无法实现的。因此,在实际操作中,可将只截取前阶特征向量而忽略小特征对,因此材料场可以重新表示为

(10)

图1 材料场观察点分布、有界场拓扑描述与级数展开和降维表征示意图Fig.1 Distribution of observation points, topological description of bounded field, series expansion and dimensionality reduction of material field

为了将基于材料场展开系数描述的有界材料场映射到整个设计域的所有单元上,且防止结构动力学优化中的虚假局部模态问题,本文引入了(RAMP)Rational Approximation of Material Properties插值模型：

(11)

()=+()(-)

(12)

式中:和分别为单元的杨氏模量和密度;、和、分别表示实心和空心单元的杨氏模量和密度;为惩罚因子,单元相对密度变量()可以表示为

(13)

其中：表示第个单元的中心坐标。这里,为了获得更光滑的材料密度场分布可以使用Heaviside 函数进行过滤。

基于MFSE方法,式(5)中的稳态和瞬态动力学拓扑优化问题,可以进一步写成

(14)

2.2 基于非梯度序列Kriging优化算法的拓扑优化问题求解

经过基于材料场级数展开(MFSE)方法处理,使得拓扑优化设计问题中的设计变量个数大幅度减少,而且单一设计变量的改变不只是改变结构拓扑结构的某一局部,而是发生整体的改变。这使得传统结构拓扑优化问题中设计变量过多而无法使用非梯度优化算法的难题得到显著的改善。对于绝大多数的结构拓扑优化问题,经过MFSE策略对结构的材料场拓扑表征与降维,可以将结构拓扑优化问题的独立设计变量降到100个以内(部分问题可以降到50以下),然而这一设计变量维度对于现有的非梯度算法来说,依然具有计算量大和难以收敛的问题。

针对结构拓扑优化问题的特性,采用了基于自适应策略的序列Kriging算法求解拓扑优化问题。首先将优化问题转化为无约束优化问题

(15)

式中：为障碍乘子。这里障碍乘子并非一个固定值,与自适应障碍因子类似,其取值需要根据目标函数的大小而改变,实际算例中这一因子只需要保证约束能够满足即可,根据数值经验,障碍因子可以取为=10floor(1+lg|()|),floor(·)表示向下取整。

为了找到无约束优化问题的全局解,基于Kriging的算法采用了初始样本和序列采加点的过程。在代理模型的初始建模中,采用拉丁超立方抽样,以得样本在多维设计空间中均匀分布。但是采用标准的拉丁超立方采样方法的过程,会产生大量无效样本数据,这将使基于Kriging的算法很难找到最优解。因此,针对结构拓扑优化问题的基本性态相对稳定的特性,提出了设计变量采样范围的自适应策略,以提高拓扑优化问题的搜索效率和求解精度,这一自适应策略可以有效地防止不切实际的采样点的产生。基于自适应策略的序列Kriging算法拓扑优化问题求解过程如图2所示。其自适应流程如下：首先,以=为中心,确定初始子设计空间,在该设计空间中进行拉丁超立方采样和代理模型加点寻优。其次,使用当前最优点(即当前代理模型的性能最优点)以调整子设计空间的中心位置并逐渐减小其范围,构建新的子优化问题。随着子优化问题的更新和求解,最终得到优化解。

图2 基于自适应策略的序列Kriging算法求解过程示意Fig.2 Solution process diagram of sequential Kriging algorithm based on adaptive strategy

3 算例及结果讨论

3.1 指定频率区间结构稳态动力学优化算例

一四边固支的夹层方板结构如图3所示,夹层板的长度为1.2 m,分为3层。基底层的厚度为0.5 mm,其密度为2 700 kg/m,弹性模量=69×10Pa,泊松比为0.3。阻尼材料层的厚度为3.0 mm,密度为980 kg/m,弹性模量=22×10Pa,泊松比为0.49。约束层的厚度及材料属性与基底层的厚度、材料属性一致。方形夹层板的中心处受到一个集中简谐激励的作用,简谐激励()=ei中=200 N,=2π,载荷频率区间为[10, 20] Hz。本文假设夹层板结构层与层之间完美连接,另外由于基底层和约束层的刚度效应远大于阻尼材料层的刚度效应,所以在计算壳单元刚度阵时阻尼材料层的偏心效应被予以忽略。本算例中,优化问题的目标函数为结构加载点动柔度的对数,其表达式如式(3)所示。

首先将结构划分为100×100共10 000个单元的网格,考虑到该方形夹层板结构为对称结构,取1/4阻尼材料层结构作为优化设计域,观察点的个数为50×50,独立设计变量个数取=50,体积约束=05。应该指出,观察点的分布可以与有限元网格划分相互独立,通过对有限元网格加密,理论上将使得分析更为精确。在本文方法的具体实施过程中,为了方便,观察点位于每一个有限单元的中心,即观察点个数等于有限元单元个数。相关长度为0.36 m(设计域边长的0.3倍),当相邻2个子优化问题目标函数之间的相对误差小于0.000 1时,优化过程停止。从数值经验来看,这一停止准则的选取将一定程度影响优化问题最终的子迭代步个数和有限元分析总次数。当这一值取值过小时,在优化的中后段优化结果的变化十分微弱,因此只需选择合适的收敛准则即能有价值的优化解。

图3 简谐荷载作用下的固支夹层板Fig.3 A square plate with four edges clamped under a time-harmonic load

图4为优化迭代历史曲线,从迭代历史曲线可以看出,在优化的每个子迭代步中包含2个阶段：① 随机采样阶段,在此阶段通过拉丁超立方采样构建Kriging代理模型,在此阶段迭代历史出现因随机采样而产生的“振荡”现象,在本例中每个子迭代步进行3×次随机采样；② 利用EI(Expectation Improvement)加点准则对于构建的代理模型进行寻优,在此阶段优化问题的目标函数将平稳下降,本算例中对每个子迭代步最多进行60次加点。在使用Kriging代理模型优化算法求解了大约15个子优化问题之后,目标函数收敛。目标函数值从初始设计中的365.61 N·m·s稳定的下降至优化设计中的105.68 N·m·s,降幅比较明显。结构的最终拓扑优化结果如图5所示,图6为结构参考设计与优化设计在0～40 Hz频域内的扫频曲线,与参考设计相比,结构在经过优化设计后,共振峰从给定的频率区间内(18.6 Hz)移至区间之外(22.1 Hz),实现了将结构的共振峰避开特殊频率区间的优化目的。

图4 目标函数f1的迭代历史Fig.4 Iterative history of objective function f1

当选取不同的给定频域区间时,阻尼材料层的拓扑优化结果会发生变化,图7给出了指定频率区间为[50, 80] Hz的拓扑优化结果,与指定频率区间为[10, 20] Hz的拓扑优化结果相比,阻尼材料层的拓扑布局变得复杂,这是因为高频率的激励荷载激发了结构更为复杂的局部模态,从而改变了阻尼材料层的拓扑布局。图8为结构参考设计与优化设计在40～100 Hz频域内的扫频曲线,可以看出夹层板结构在经过优化布局后,其共振峰也从给定的频率区间内移至给定区间之外。表1给出了考虑载荷频率区间为[50, 80] Hz时的优化结果和参考解的前10阶固有频率对比列表,可以看出相比于仅考虑[10, 20] Hz的动力学优化相比,考虑[50, 80] Hz频率区间的优化问题包含了结构更高阶的固有频率,其频率区间内的动力学优化减振体现了多个振动模态下的综合减振效应。与已有阻尼减振结构拓扑优化相比,针对同样的问题本文的优化结果与已有方法十分相似,说明本文方法计算结果具有正确性。

图5 阻尼材料层优化结果(频率区间[10, 20] Hz)Fig.5 Optimization result of damping layer (frequency range[10, 20] Hz)

图6 结构指定位置的扫频曲线优化解(频率区间[10, 20] Hz)与参考解(ρe=0.5)的比较Fig.6 Comparison between optimal solution (frequency range[10, 20] Hz) and reference solution (ρe=0.5) of sweep curve at specified position of structure

图7 阻尼材料层优化结果(频率区间[50, 80] Hz)Fig.7 Optimization result of damping layer (frequency range[50, 80] Hz)

图8 结构指定位置的扫频曲线优化解(给定频率区间[50, 80] Hz)与参考解(ρe=0.5)的比较Fig.8 Comparison between optimal solution (given frequency range[50, 80] Hz) and reference solution (ρe=0.5) of sweep curve at specified position of structure

表1 优化结构与参考结构各阶固有频率比较

在运载航天器系统中,载荷适配器可以为载荷与航天器之间提供一个缓冲截面。通过在载荷适配器表面铺设阻尼层可以减小航天器与卫星之间的振动。本文载荷适配器结构由圆柱段、圆锥台段与法兰3部分构成,其上部的直径为500 mm,底部的直径为700 mm,如图9所示。圆柱段的高度和厚度分别为200 mm和1.5 mm,圆锥台段的高度与厚度分别为200 mm和0.8 mm,法兰的外直径为740 mm,厚度为2 mm。结构的圆锥台段铺设阻尼层,厚度为1 mm,结构基底与阻尼的材料属性与算例1一致。

图9 简谐激励下载荷适配器示意图Fig.9 Payload adapter subjected to four harmonic concentrated forces

结构整体划分为240×50个有限元单元,对于此周向循环结构,取阻尼层的周向1/4作为优化设计域,观察点的个数为60×30,设计变量个数为50,相关长度取为阻尼层母线长度的0.3。圆柱段的顶部受到固定约束,4个简谐激励等距地加载在法兰底部内侧的纵向,其幅值为=10 kN,载荷频率区间为[360, 420] Hz, 设计要求在荷载频域内无共振峰值。设计域体积约束上限为=05,优化问题的目标函数取为4个加载点动柔度之和。优化结果如图10所示,图11 为参考解(=05)与优化解的扫频曲线比较。

图10 阻尼材料层拓扑优化结果Fig.10 Optimization result of damping layer

图11 结构指定位置的扫频曲线优化解与参考解(ρe=0.5)比较Fig.11 Comparison between optimal solution and reference solution (ρe=0.5) of sweep curve at specified position of structure

3.2 瞬态动力载荷结构拓扑优化算例

如图12所示,悬臂矩形夹层板的尺寸为0.6 m×0.4 m,由厚度为1.0 mm的阻尼材料层与厚度均为0.5 mm的基底层和约束层组成。该夹层板结构各层的材料属性与算例1夹层板结构各层的材料属性一致。夹层板的自由边中点受到冲击荷载()的作用,其时程曲线如图13 所示。在该算例中目标函数所考虑的时间段终止时刻为=0.05 s,时域积分方法的时间步长Δ=2×10s。

图12 冲击荷载作用下的悬臂夹层板Fig.12 A cantilever plate under impact load

整个结构被划分为60×40个单元,相关长度为0.12 m(设计域短边的0.3倍),设计变量个数取=50,阻尼材料的体积约束分数取=05。优化问题的目标函数取为自由边中点振幅的平方在时间段[0,]内的积分。本算例中,目标函数为结构指定位置瞬态的动响应,其表达式如式(4)所示。

图13 悬臂夹层板冲击荷载时程曲线Fig.13 Impact force applied to cantilever plate

优化迭代曲线如图14所示,优化过程在经过15个子迭代步,共2 749步迭代后收敛,目标函数从初始设计中的10.97 m·s平稳的下降至4.78 m·s,下降了56.43%。图15为最终的阻尼材料层拓扑优化布局。初始设计和优化设计的目标位置时程位移响应如图16所示,可以发现结构在经过优化设计后,其振动幅值有明显的降低,振动频率也低于结构初始设计的振动频率,这说明拓扑优化解良好地实现了结构振动控制。

图14 目标函数f2的迭代历史Fig.14 Iterative history of objective function f2

图15 阻尼材料层拓扑优化结果(T=0.05 s)Fig.15 Optimization result of damping layer (T=0.05 s)

修改时间段终止时刻为0.2 s,其优化结果如图17所示。可以看出,当改变目标函数所考虑的时间段终止时刻,拓扑优化结果发生了明显的变化。图18为初始设计与优化设计的目标位置位移时程响应,可以发现结构在经过优化设计后,在时间历程前期其振动衰减得更加迅速,在时间历程后期其振动趋于平稳。

图16 结构指定位置位移时程响应的优化解(T=0.05 s)与参考解(ρe=0.5)比较Fig.16 Comparison between optimal solution (T=0.05 s) and reference solution (ρe=0.5) for displacement time history response of structures at specified positions

图17 阻尼材料层拓扑优化结果(T=0.2 s)Fig.17 Optimization result of damping layer (T=0.2 s)

图18 结构指定位置位移时程响应的优化解(T=0.2 s)与参考解(ρe=0.5)比较Fig.18 Comparison between optimal solution (T=0.2 s) and reference solution (ρe=0.5) for displacement time history response of structures at specified positions

4 结 论

本文研究了基于材料场级数展开策略的结构动力学非梯度拓扑优化方法,该方法能够利用少量拓扑场描述变量实现复杂结构的拓扑表征,从而大幅降低结构拓扑优化变量个数,进一步采用自适应序列Kriging代理模型算法求解拓扑优化问题,实现了基于非梯度优化算法的大规模结构动力学拓扑优化问题的求解。所提方法因不需要灵敏度信息,使得优化问题目标函数和约束函数的选取与以往动力学优化问题相比具有更大的自由度,更适用于复杂动力学优化问题的建模与求解。数值算例表面,所提优化问题能很好地处理结构在稳态与瞬态动力学荷载下的拓扑优化问题。

本文所处理的2类优化问题设计目标,属于常规结构动力学优化目标,可以通过灵敏度分析和梯度算法进行求解和设计,但考虑频率区间或时间段内的振动优化问题中,由于原目标函数不连续,常常需要引入包络函数等技术,使得灵敏度更为复杂,部分优化问题采用基于梯度的优化算法因多局部解效应导致收敛性较差。同时,该方法因不需要灵敏度信息,易与各种仿真分析软件结合,未来在复杂物理场和非线性结构拓扑优化问题(如多物理场、非线性优化、接触与碰撞)中,具有较大的应用前景,目前已经在声学与光学超材料设计、大变形结构拓扑优化、流固耦合问题优化等方面进行了初步尝试。

值得注意的是,此前已有文献给出了10万以上自由度的工程结构算例,本文所提出的基于材料场级数展开的结构拓扑优化方法需要的总计算次数(有限元分析)与传统梯度类优化算法相比较多。然而,随着计算机硬件技术和并行计算能力的飞速发展,本文在未来大规模动力学优化问题仍然有较大的应用潜力。