基于广义模态截断扩增方法的结构频响拓扑变量灵敏度分析

周大为,陈飙松,李云鹏,张 盛

(大连理工大学工业装备结构分析国家重点实验室运载工程与力学学部工程力学系,大连 116024)

摘 要:基于结构拓扑优化设计中的变密度法,采用伴随法推导结构在简谐激励下的频响振幅对单元设计变量的解析灵敏度列式。针对灵敏度数值求解中传统模态位移法引发的低精度问题,通过引入计算成本较低的广义模态截断扩增方法提升灵敏度的计算精度。数值算例将该文方法与全局有限差分法和其他灵敏度计算方法进行了比较。结果证明了该文方法在不同的激励频率及有限元网格密度下高效求解高精度灵敏度的有效性。

关键词:结构拓扑优化;解析灵敏度列式;频响振幅;伴随法;模态位移法;广义模态截断扩增法

在结构拓扑优化的变密度法中,灵敏度是结构某状态变量对单元设计变量的导数,表征了某一单元设计变量对该状态变量的影响程度,是优化过程中各设计变量变化的重要依据。在灵敏度分析中,精度和效率是需要重点考虑的两个方面。首先,精度较低甚至错误的灵敏度会给出错误的优化方向,使优化结果不可信;其次,大规模的优化问题导致设计变量的数目呈指数级增长,低效率的灵敏度算法将十分耗时,使计算成本过高。因此,研究同时具有高精度和高效率的灵敏度算法是十分必要的。

另一方面,继结构自振频率响应[1-2]拓扑优化问题,作为工程中常见的动力响应,结构在简谐激励下的稳态频响拓扑优化问题已受到广泛的关注。Ma和Kikuchi等[3]给出了稳态频响的动柔顺性指标的表达式,并使用均匀化方法结合改进优化准则法对无阻尼频响问题进行了拓扑优化。Jog[4]推导了周期性荷载下结构稳态响应的灵敏度并进行了拓扑优化。Jensen[5]研究了结构受频带简谐激励的拓扑优化问题。Yoon[6]比较了不同缩减方法在灵敏度计算和拓扑优化中的表现。Du和Olhoff[7―8]讨论了在频响动柔顺性拓扑优化问题中可能出现的结构变“柔”的现象,并提出GIF方法。Kang等[9]研究了自由阻尼层板结构在简谐激励下的阻尼材料分布问题。叶红玲等[10]基于 ICM 方法研究了以频带内频响振幅响应为约束的拓扑优化问题。Liu和Zhang等[11]研究了在频响拓扑优化问题中,分别使用模态加速法、直接法和模态位移法对灵敏度计算和优化结果的影响。

本文针对稳态频响问题,推导了结构某点频响振幅的灵敏度列式,在计算中引入广义模态截断扩增法,指出了该方法对显著提高灵敏度精度并兼顾计算效率产生的关键作用。数值算例的对比验证说明了本文灵敏度算法的有效性。

1 结构某点频响振幅响应及其灵敏度

1.1 结构某点频响振幅

结构在外载荷下的线性动力方程为:

其中:MK分别是n阶质量矩阵和刚度矩阵;ft是时间相关外载向量;x是结构位移向量;fD是阻尼力向量,n是结构自由度数目。

阻尼作用机理十分复杂。在稳态频响问题中,出于简化计算的目的,常采用两类阻尼模型:瑞利阻尼模型和复刚度阻尼模型。瑞利阻尼模型假设结构的阻尼阵与其质量阵和刚度阵呈线性比例关系:

结构所受阻尼力为:

结构线性动力方程为:

复刚度阻尼模型则引入复刚度的概念:

其中:ς为复刚度阻尼系数;i是虚数单位。结构线性动力方程表示为:

稳态频响问题的外载荷向量及位移向量表示为:

将式(7)代入式(4),可以得到瑞利阻尼下频响方程:

将式(7)代入式(6),并定义复刚度阻尼阵S=ςK,可以得到复刚度阻尼下频响方程:

定义结构的动刚度矩阵:

频响位移U可以表示为:

于是结构某一自由度的频响振幅为:

所关心的振幅方向可能并非该节点某个自由度的方向,因此可以定义结构某节点沿某方向的频响振幅A

其中,vn维行向量,定义了该节点的振幅关心方向。

1.2 基于伴随法的灵敏度推导

在式(14)两端对单元设计变量求导:

t,式(15)可以改写为:

下面用伴随法求解式(16)。对频响方程求导:

显然有:

其中,Λ是伴随向量。若令向量Λ满足伴随方程:

则有:

式(19)存在动刚度矩阵W的转置,通常情况下结构的质量阵和刚度阵是对称的,而根据本文中阻尼阵的定义可知,阻尼阵是对称的,因此动刚度矩阵也是对称的。伴随方程可以改写为:

式(20)可进一步化为单元级别的计算:

其中,weλeue分别是单元级别的动刚度阵、伴随向量和位移向量。设单元刚度阵、单元质量阵及单元阻尼阵的插值模型分别为:

式中:为单元初始刚度阵;为单元初始质量阵;为单元初始阻尼系数,则容易得到∂we/∂eρ,最后便可得到结构某点沿某方向的频响振幅的灵敏度,当向量v是单位向量时,即退化为某自由度频响振幅的灵敏度。

2 广义模态截断扩增法

由灵敏度列式可以看出求解灵敏度的关键在于求解伴随向量和频响位移。伴随式(21)本质上是一个频响方程,只有右端项不同。因此灵敏度求解的精度和效率最终取决于频响方程求解的精度和效率。

2.1 直接法和模态位移法

求解频响方程的方法主要有两类:直接法和模态位移法。直接法不对原方程进行任何处理,直接求解该复系数线性方程。该方法没有任何对解的近似,理论上可以得到问题的“真实”解,但计算成本过高。

模态位移法以结构的若干低阶模态的线性叠加近似结构的位移解,即前p阶模态截断:

其中,Φ是由前p阶模态组成的模态矩阵。

以瑞利阻尼下的频响方程为例,将式(24)代入频响方程式(8),方程两端左乘TΦ,方程的维数缩减为p阶,记ΦTΦTΦT,可以得到:

缩减后的方程规模大大降低,可以直接求解。得到q后再由式(24)得到近似的位移解。

2.2 广义模态截断扩增法

模态位移法虽然大幅降低了计算成本,但忽略高阶模态对位移的贡献导致了难以避免的精度损失。尤其对激励处于中高频段的问题,其结果与直接法相差较大。虽然可以取更多高阶模态参与缩减,但高阶模态计算成本较高,精度也难以保证。

一个重要的改进是由 Williams[12]于 1945年提出的模态加速法,该方法首先计算修正向量:

而位移的近似表达式则改写为:

Zhang[11]在频响拓扑优化中采用了该种方法,与模态位移法相比,灵敏度精度得到明显提升。另一种做法是将向量mφ连同模态向量一起参与缩减后再得到近似的位移解,即模态截断扩增法:

上述两种改进方法与模态位移法相比,虽然拥有显著的精度提升,但在更广的激励频率范围内,仍无法保证与直接法一致的精度。

基于模态截断扩增法,Rixen[13]提出了高阶修正的广义模态截断扩增法(GMTAM):

其中,j为修正阶数。将按式(29)计算的若干阶修正向量连同模态向量一同参与到缩减中,可以得到近似的位移解:

使用该方法求解可以在较广的频域内大幅提升求解精度,且修正向量的计算比较简单,不会显著增加计算量。当修正阶数j只取1时,该方法即退化为一般的模态截断扩增法。

值得注意的是,若结构所受载荷f满足线性关系,则在计算灵敏度时不必再计算对应的修正向量,只需重复使用f对应的修正向量即可。

3 数值算例

本节用本文方法计算频响振幅灵敏度,并与其他几种方法给出的结果进行比较。算例以有限差分直接法的灵敏度结果作为比较标准。设计变量初始值为0.5,增量为0.001。插值函数为

3.1 四边固定约束平板算例

如图 1所示四边固定平板模型,平板尺寸为3 m×2 m,板厚0.02 m。划分为四边形板单元网格。材料属性为:弹性模量70 GPa,泊松比0.3,密度2700 kg/m3。材料具有复刚度阻尼,复刚度阻尼系数为ς=0.5。板中心受竖直向下简谐激励,激励频率分别为20 Hz及300 Hz,幅值为100 kN。计算板中心点处法线方向振幅响应的灵敏度。将平板划分为较稀疏的30×20个及较致密的60×40个四边形板单元,分别对比各方法的计算结果。

图1 四边固定平板模型
Fig.1 Four-edge-clamped flat plate model

3.1.1 稀疏网格下的灵敏度

根据对称性,截取1/4稀疏网格模型,如图2。测试单元号为311、316、411、416、511及516。

图2 1/4平板模型
Fig.2 A quarter of the plate model

首先考察激励频率为20 Hz时的灵敏度。各方法给出的灵敏度如表1(a)和表1(b)所示。图3展示了模态位移法与有限差分直接法的灵敏度的相对误差。

表1(a)稀疏网格下各测试单元灵敏度结果对比表—广义模态截断扩增法—20 Hz
Table 1(a)Comparison of sensitivity results of each test element under thin mesh-GMTAM-20 Hz

表1(b)稀疏网格下各测试单元灵敏度结果对比表—模态位移法—20 Hz
Table 1(b)Comparison of sensitivity results of each test element under thin mesh-MDM-20 Hz

图3 稀疏网格下灵敏度结果相对误差—模态位移法—20 Hz
Fig.3 Relative differences of sensitivity results under thin mesh-MDM-20 Hz

20 Hz低于该结构一阶自振频率 33 Hz,是低频。由表1(a)可以看出,对于各测试单元,本文方法、解析直接法与有限差分直接法三种方法给出的灵敏度吻合很好。对于GMTAM,模态截断阶数p仅需取4,修正阶数r仅需取2就可以得到精度很高的结果。

在表1(b)中,当模态截断阶数p增加时,模态位移法给出的灵敏度并未收敛:多数单元的灵敏度出现“震荡”,图 3直观地表明了这一现象。对于311号单元,即使p取80,与有限差分直接法给出的灵敏度相比仍有27%的差别。

现在考察激励频率为300 Hz时的灵敏度。本文方法给出的灵敏度如表2所示。

表2 稀疏网格下各测试单元灵敏度结果对比表—广义模态截断扩增法—300 Hz
Table 2 Comparison of sensitivity results of each test element under thin mesh-GMTAM-300 Hz

300 Hz是中高频。由表2可以看出,本文方法、解析直接法与有限差分直接法给出的灵敏度基本一致。对于GMTAM,模态截断阶数取30,修正阶数取4可以得到精度很高的结果,这是因为激振频率增高,结构的更多高阶模态被激起参与到振动中。图4表明模态位移法的灵敏度收敛性仍然很差。

图4 稀疏网格下灵敏度结果相对误差—模态位移法—300 Hz
Fig.4 Relative differences of sensitivity results under thin mesh-MDM-300 Hz

表3 致密网格下各测试单元灵敏度结果对比表—广义模态截断扩增法—20 Hz
Table 3 Comparison of sensitivity results of each test element under fine mesh-GMTAM-20 Hz

表4 致密网格下各测试单元灵敏度结果对比表—广义模态截断扩增法—300 Hz
Table 4 Comparison of sensitivity results of each test element under fine mesh-GMTAM-300 Hz

图5和图6则进一步表明,提高模态位移法的模态截断阶数无法保证灵敏度精度的提升。在致密网格下,模态位移法给出的灵敏度的相对误差明显增加,这是因为网格加密后,结构出现更多运动细节,更多的高阶模态具有对结构位移的贡献,从而收敛速度减慢。

GMTAM在保证精度的同时,计算成本对比直接法有显著下降,见表5。可以看到,本文方法比直接法节省了约一半的时间,仅比模态位移法多消耗了2 s,由此可证明使用GMTAM计算灵敏度是高效的。

图5 致密网格下灵敏度结果相对误差—模态位移法—20 Hz
Fig.5 Relative differences of sensitivity results under fine mesh-MDM-20 Hz

图6 致密网格下灵敏度结果相对误差—模态位移法—300 Hz
Fig.6 Relative differences of sensitivity results under fine mesh-MDM-300 Hz

表5 灵敏度计算时间对比表
Table 5 Comparison of computation time using different sensitivity analysis methods

3.2 两端固定约束圆柱壳算例

如图 7(a)和图 7(b)所示圆柱壳模型,半径为1 m,高2 m,厚0.02 m。划分为40×40个四边形壳单元网格。材料属性与上一算例相同。圆柱壳在z方向对称面内的四个节点上受径向简谐激励,激励频率分别为80 Hz及200 Hz,幅值为100 kN。考察其中一个节点处径向(y轴方向,见图8)振幅响应的灵敏度。

图7 两端固定圆柱壳结构模型
Fig.7 Cylindrical shell Structure with two end fixed

根据对称性,取1/16结构如图8所示。测试单元如图8所示。为方便结果显示,将单元311、313、315、471、473、475、631、633、635、791、793、795依次重新编号为1号~12号单元。

图8 1/16圆柱壳模型
Fig.8 One sixteenth of the cylindrical shell model

图9和图10分别给出了激励频率为80 Hz及200 Hz的灵敏度对比图。

对于该算例,无论激励频率是高频或是低频,本文方法与有限差分直接法所得的灵敏度均保持一致,具有很好的收敛性。

图9 各测试单元灵敏度结果—广义模态截断扩增法—80 Hz
Fig.9 Sensitivity results of each test element-GMTAM-80 Hz

图10 各测试单元灵敏度结果—广义模态截断扩增法—200 Hz
Fig.10 Sensitivity results of each test element-GMTAM-200 Hz

4 结论

本文在变密度法结构拓扑优化框架下,研究了结构受简谐激励的频响振幅响应的灵敏度算法。指出在计算中使用广义模态截断扩增方法可以在保证计算效率的同时大幅提升计算精度。数值算例验证了该算法的实用性。研究工作表明,本文算法在激励频率为低频时只需很低的计算成本便可得到高精度的灵敏度结果;对中高激励频率在致密网格下的灵敏度计算同样有理想的效果;直接法计算精度高,但计算成本高昂;模态位移法则无法保证灵敏度的计算精度。

参考文献:

[1]薛开,雷寰兴,王威远.一种新的周长约束方法在阻尼频率拓扑优化中的应用[J].工程力学,2013,30(6):275-280.Xue Kai,Lei Huanxing,Wang Weiyuan.An application of a new perimeter constraint method in topological optimization for damped frequency [J].Engineering Mechanics,2013,30(6): 275-280.(in Chinese)

[2]李东泽,于登云,马兴瑞.基频约束下的桁架结构半定规划法拓扑优化[J].工程力学,2011,28(2): 181-185.Li Dongze,Yu Dengyun,Ma Xingrui.Truss topology optimization with fundamental frequency constraints via semi-definite programming [J].Engineering Mechanics,2011,28(2): 181-185.(in Chinese)

[3]Ma Z D,Kikuchi N,Hagiwara I.Structural topology and shape optimization for a frequency response problem [J].Computational Mechanics,1993,13(3): 157-174.

[4]Jog C S.Topology design of structures subjected to periodic loading [J].Journal of Sound and vibration,2002,253(3): 687-709.

[5]Jensen J S.Topology optimization of dynamic problems with Padé approximants [J].International Journal Numerical Methods in Engineering,2007,72(13):1605-1630.

[6]Yoon G H.Structural topology optimization for frequency response problem using model reduction schemes [J].Computer Methods in Applied Mechanics and Engineering,2010,199(25-28): 1744-1763.

[7]Olhoff Niels,Du Jianbin.Topological design for minimum dynamic compliance of structures under forced vibration [M]// Rozvany G I N,Lewiński T.Topology Optimization in Structural and Continuum Mechanics,CISM International Centre for Mechanical Sciences,vol 549.Vienna: Springer,2014: 325-329.

[8]Olhoff Niels,Du Jianbin.Generalized incremental frequency method for topological design of continuum structures for minimum dynamic compliance subject to forced vibration at a prescribed low or high value of the excitation frequency [J].Structural and Multidisciplinary Optimization,2016,54(5): 1113-1141.

[9]Kang Z,Zhang X,Jiang S,et al.On topology optimization of damping layer in shell structures under harmonic excitations [J].Structural and Multidisciplinary Optimization,2012,46(1): 51-67.

[10]叶红玲,李耀明,张颜明,等.基于对数型 Heaviside近似函数作为过滤函数的动力响应结构拓扑优化ICM方法[J].工程力学,2014,31(6): 13-20.Ye Hongling,Li Yaoming,Zhang Yanming,et al.Structural topology optimization of frequency response problem by approximately logarithmic Heaviside function based on ICM method [J].Engineering Mechanics,2014,31(6): 13-20.(in Chinese)

[11]Hu Liu,Weihong Zhang,Tong Gao.A comparative study of dynamic analysis methods for structural topology optimization under harmonic force excitations [J].Structural and Multidisciplinary Optimization,2015,51(6): 1321-1333.

[12]Williams D E.Dynamic loads in aeroplanes under given impulsive loads with particular reference to landing and gust loads on a large flying boat [R].UK: Great Britain Royal Aeronautical Establishment Reports.SME,1945:3309-3316.

[13]Daniel Rixen.Generalized mode acceleration methods and modal truncation augmentation [C].Seattle,WA,U.S.A.: 42nd AIAA/ASME/ASCE/AHS/ASC Structures,Structrual Dynamics and Material Conference and Exhibition,AIAA Paper 2001-1300,2001.

SENSITIVITY ANALYSIS OF STRUCTURAL TOPOLOGY DESIGN VARIABLES UNDER HARMONIC EXCITATIONS BASED ON GENERALIZED MODAL TRUNCATION AUGMENTATION METHOD

ZHOU Da-wei,CHEN Biao-song,LI Yun-peng,ZHANG Sheng
(State Key Laboratory of Structural Analysis for Industrial Equipment,Department of Engineering Mechanics,Facxlty of Vehicle Engineering and Mechanics,Dalian University of Technology,Dalian 116024,China)

Abstract:Based on the variable density method for structural topology optimization,the analytical sensitivity formulation of the frequency response displacement amplitude of structures under harmonic excitations is proposed using the adjoint method.The generalized modal truncation augmentation method is introduced to obtain high accuracy without high computational cost in contrast to the poor accuracy of the traditional modal displacement method in sensitivity computation.Numerical examples are presented to compare the proposed method with the global finite difference method and other computation methods.The computational results demonstrate the effectiveness of the proposed method in computing accurate sensitivities with high efficiency under different excitation frequencies and different densities of finite element meshes.

Key words:structural topology optimization; analytical sensitivity formulation; frequency response displacement amplitude; adjoint method; modal displacement method; generalized modal truncation augmentation method

张 盛(1976―),男,吉林人,博士,从事计算力学及其软件开发研究(E-mail: zhangs@dlut.edu.cn).

李云鹏(1970―),男,辽宁人,硕士,从事计算力学及其软件开发研究(E-mail: lyp@dlut.edu.cn);

周大为(1990―),男(满族),吉林人,博士生,从事结构优化及计算力学研究(E-mail: Zhoudw@mail.dlut.edu.cn);

作者简介:

通讯作者:陈飙松(1973―),男,广东人,教授,博士,从事结构优化及计算力学软件开发研究(E-mail: Chenbs@dlut.edu.cn).

基金项目:国家重点研发计划项目(2016YFB0200702);国家自然科学基金项目(1372064,11761131005);高等学校学科创新引智计划项目(B14013)

收稿日期:2017-08-11;修改日期:2018-03-29

文章编号:1000-4750(2018)11-0001-07

doi:10.6052/j.issn.1000-4750.2017.08.0618

文献标志码:A

中图分类号:O327;TB535.1