非比例阻尼结构谐响应分析的自适应模型降阶方法①

2023-07-08李玉韦柏汉松

固体火箭技术 2023年3期

李玉韦,魏晓,曹航,柏汉松,王博,郝鹏

(1.中国航发沈阳发动机研究所,沈阳 110015;2. 大连理工大学工程力学系,大连 116024)

0 引言

为了在服役过程中控制结构的振动水平,通常需要获得结构在简谐激励下的稳态响应。工程中常用的谐响应分析方法包括直接频响法、模态叠加法和模态加速度法等。直接频响法因其计算精度高,常被应用于计算结构的稳态响应[1],但该方法计算效率较低,不适合处理激励频率较为密集或自由度数目较大的工程问题。模态叠加法通过引入模态正交化条件解耦动力学方程,提高了复杂结构的谐响应分析效率,被广泛应用于复杂工程结构[2-3]。由于模态叠加法忽略高阶模态的影响,导致该方法计算精度降低。为提高模态叠加法的计算精度,需要选择更多模态参与谐响应分析,但对于如何选择模态以及选择多少模态参与计算还缺乏科学的准则。为了弥补模态叠加法忽略高阶模态导致计算精度降低的问题,模态加速度法引入拟静力响应来补偿截断模态的贡献,从而达到用较少模态获得高精度响应的效果[4-5]。模态叠加法或模态加速度法计算结构稳态响应的高效性体现在结构阻尼矩阵满足模态正交化条件,从而可以对动力学方程进行解耦,提高响应分析的效率。但是,对于由不同材料和构件组成的非比例阻尼结构,阻尼矩阵不再满足模态正交化条件,此时,模态叠加法或模态加速度法不再适用,否则会导致错误的计算结果[6-7]。

Krylov子空间法[8]通过构建一组标准正交基来实现模型降阶,提高结构响应计算效率。这种模型降阶方法数学理论完善、算法稳定,亦不要求结构阻尼矩阵满足模态正交化条件,因此该方法被广泛应用于工程结构的模型降阶[9-10]。一阶Krylov子空间法应用于结构模型降阶时,需要将二阶动力学方程转换为一阶状态方程,再对状态方程降阶[11]。这种线性化处理方式的优点是收敛速度快,但线性化过程破坏了原结构的矩阵特征,计算规模是原动力学方程的两倍,计算量显著增加。为此,基于二阶Krylov子空间作为前射子空间,柏兆俊和苏仰峰提出了二阶Arnoldi(Second-order Arnoldi Reduction,SOAR)算法[12-14],该算法不会破坏结构动力学方程的二阶特性和结构矩阵特点,在保证分析精度的前提下,可以有效减少数值分析的复杂度,提高计算效率。

针对SOAR方法需要使用昂贵的矩阵反演生成正交投影矩阵且无法保持原系统稳定性的问题,MALIK 等[15]用高斯核加权的插值点在期望的频率范围内生成减缩基矩阵。TAMRI等[16]提出基于数值秩性能系数的概念来自动选择最优的降阶模型。本文提出一种简单的自适应模型降阶方法,利用交叉验证和二分法策略确定展开频点数和正交基阶数,自适应地建立在目标频段内具有更高计算精度的降阶模型,提高非比例阻尼结构谐响应分析效率。

1 谐响应分析方法

简谐激励载荷下结构的频响方程为

(K+jωC-ω2M)U(ω)=F

(1)

直接频响法将一系列离散频点直接代入式(1)得到结构在频点ω处的响应:

U(ω)=(K+jωC-ω2M)-1F

(2)

可以看出,直接频响法需要得到每个频点ω处动刚度阵的逆矩阵,计算规模较大,不适用于激励频率密集或自由度数目较多的问题。

1.1 模态叠加法

模态叠加法是把对应无阻尼结构的模态为空间基底,经过坐标变换使得原动力学方程解耦,通过求解n个独立的方程获得模态位移,进而通过叠加各阶模态的贡献求得结构的响应。

假定结构特征值和模态分别是Λ和Φ:

Φ=[φ1,φ2,…,φn]

(3)

其中,Λ和Φ满足特征值方程和正交性条件:

KΦ=ΛMΦ,ΦTKΦ=Λ,ΦTMΦ=I

(4)

若阻尼为瑞利阻尼,则阻尼阵C满足模态正交化条件:

ΦTCΦ=2ξdiag(ω1,ω2,…,ωn)

(5)

式中ξ为瑞利阻尼比。

利用式(4)和式(5)将式(1)解耦为

(6)

可以看出,模态叠加法通过模态的线性组合来近似结构的位移,避免了动刚度阵的分解,提高了计算效率。实用中,模态叠加法一般是通过模态截断选取有限数目的模态近似表示结构的位移:

(7)

式中l为选取的模态数,l<

模态叠加法忽略了高阶模态,使得该方法在高频范围内的求解精度降低。

1.2 模态加速度法

对于无阻尼结构,式(1)可以改写为

U(ω)=K-1F+ω2K-1MU(ω)

=K-1F+ω2K-1MΦ(Λ-ω2I)-1ΦTF

(8)

将式(4)代入式(8)可得

U(ω)=K-1F+ω2K-1MU(ω)

=K-1F+ω2K-1MΦ(Λ-ω2I)-1ΦTF

(9)

忽略高阶模态,式(9)可变换为如下形式:

(10)

可以看出,模态加速度法通过引入拟静力响应来补偿截断模态的贡献,从而达到用较少模态获得高精度响应的效果。

2 基于SOAR的自适应模型降阶方法

2.1 SOAR算法

在任一频点ωa处,频响方程(1)的等价形式如下:

(11)

其中,

(12)

通过二阶Krylov子空间对式(12)进行降阶,对应减缩基矩阵如下:

Tk(ωa)=span{r0(ωa),r1(ωa),…,rk-1(ωa)}

(13)

其中,r0(ωa),r1(ωa),…,rk-1(ωa)代表在频点ωa处展开的k个标准正交基向量:

(14)

建立减缩基矩阵Tk(ωa)后,式(11)可写为

(15)

其中,

(16)

可以看出,基于二阶Krylov子空间建立降阶模型与展开频点数及标准正交基阶数有关,展开频点数和标准正交基阶数越多,降阶模型的计算精度越高。为平衡降阶模型的计算精度和效率,本文提出了自适应确定展开频点数和标准正交基阶数的降阶策略。

2.2 自适应降阶策略

结构动柔顺度[17]误差常被用于衡量降阶模型的计算精度:

(17)

式(17)需要获得所有频点的位移来计算结构的动柔顺度误差,这使得建立降阶模型过程中需要浪费大量的时间去评估降阶模型的计算精度,失去了利用降阶模型提高计算效率的意义。本文提出利用交叉验证的策略来评估降阶模型的预测精度:

(18)

交叉验证策略仅需计算展开频点的位移精确解,大大提高了建立降阶模型的效率。为了在保证降阶模型计算精度的前提下,尽可能减少展开频点的数目,本文建立了基于二分法的自适应降阶策略,具体流程如图1所示。

图1 基于SOAR算法的自适应模型降阶流程Fig.1 Flowchart of the SOAR-based adaptive model order reduction method

第一步:设置每个频点的初始标准正交基阶数r和最大正交基数目Maxorder,正交基阶数增加的步长p;设置结构动柔顺度误差阈值Tol;初始化集合S1={ω1}和集合S2={ω2},其中ω1和ω2分别为目标频段的上下限。

第二步:利用SOAR算法得到展开频点ωa处的r个标准正交基,并组装成减缩基矩阵Tr(ωa):

ωa=max{S1}

Tr(ωa)=span{r0(ωa),r1(ωa),…,rr-1(ωa)}

(19)

ωb=min{S2}

(20)

第四步:通过SOAR算法扩充展开频点ωa处的标准正交基数目,并组装成减缩基矩阵Tr+p(ωa):

Tr+p(ωa)=span{r0(ωa),r1(ωa),…,rr+p-1(ωa)}

(21)

第五步:利用减缩基矩阵Tr+p(ωa)计算展开频点处的动柔顺度误差:

(22)

S1=S1∪{ωc}

(23)

S2=S2∪{ωa}

S1=S1{ωa}

(24)

若移除频点ωa后,集合S1仍为非空集合,则返回至第二步,否则迭代结束。

迭代结束后,将集合S2中所有频点的标准正交基进行正交化得到减缩基矩阵T:

T=orth[Tr+p(ω1),Tr+p(ω2),…,Tr+p(ωi)]

(25)

其中,orth表示对矩阵进行正交化。本文采用奇异值分解技术(Singular Value Decomposition,SVD)进行正交化。

3 自适应模型降阶方法验证

3.1 阻尼涂层板的谐响应分析

本节采用阻尼涂层板验证自适应降阶策略有效性,阻尼涂层板结构长300 mm,宽25 mm ,厚3.5 mm,其中阻尼涂层厚2 mm,如图2所示。基体材料的弹性模量E=71 GPa,密度ρ=2700 kg/m3,泊松比μ=0.3,阻尼比ξ=0.02。阻尼涂层材料的弹性模量E=50 GPa,密度ρ=1400 kg/m3,泊松比μ=0.3,阻尼比ξ=0.05。基体结构右端固支,左端施加简谐激励,激励频率的范围为[0, 800 Hz], 分析步长为2 Hz。

图2 阻尼涂层平板Fig.2 A plate with free-layer damping

选取初始标准正交基阶数r=3、步长p=1,最大正交基数目Maxorder=20,根据图1给出的自适应模型降阶流程,在目标频段内确定了5个展开频点以及每个频点处所需标准正交基的数目,如表 1所示。将5个展开频点生成的标准正交基进行SVD正交化得到降阶模型的减缩基矩阵T。

表1 自适应确定的展开频点数和标准正交基阶数Table 1 Expansion frequency points and its related orders of orthonormal basis frequency points

图3 (a)分别给出了基于全阶模型(FOM)和降阶模型(ROM_T)得到的动柔顺度曲线。可以看出,降阶模型与全阶模型的计算结果十分吻合,相对误差最大值不超过1×10-7(见图3(b))中黑色实线)。

(a) Dynamic compliance

(b) Relative error图3 阻尼涂层平板的动柔顺度及相对误差Fig.3 Dynamic compliance and relative error of the free-layer damping plate

为了验证本文提出的交叉验证策略有效性,图3 (b)分别给出了由单个展开频点的正交基建立的降阶模型预测误差。图中括号内的数字表示展开频点,下标表示标准正交基的阶数,例如,T4(200)表示展开频点为200 Hz,该展开频点处的标准正交基阶数为4。图中竖直灰色实线代表展开频点,水平灰色虚线代表阈值Tol。可以看出,降阶模型(ROM_T4(200))在频点200 Hz附近频段内的预测精度较高,在远离200 Hz的频段内的预测精度较差。同样,降阶模型(ROM_T6(400))在频点400 Hz附近频段内的预测精度较高,在远离400 Hz频段内的预测精度较差。其余降阶模型具有类似规律。

提取表 1中5个降阶模型(ROM_T20(0)、ROM_T3(100)、ROM_T4(200)、ROM_T6(400) 、ROM_T5(800))计算动柔顺度误差的下包络,如图3(b)中的红色实线所示。可以看出,红色实线完全位于水平灰色虚线的下方,表明展开频点间的预测误差亦满足精度要求,这也证明了本文利用交叉验证策略来评估降阶模型的预测精度是合理有效的。

图4进一步讨论了增加展开频点数量及正交基阶数的必要性。图4中纵轴表示动柔顺度相对误差的最大值,横轴表示正交基的阶数,红色矩形框表示增加展开频点的位置,水平灰色虚线表示误差阈值。

图4 动柔顺度最大相对误差与标准正交基数目的关系Fig.4 Maximum relative error of dynamic compliance with respect to the order of orthonormal basis

由图4可以看出,当仅有1个展开频点时,随着正交基阶数的增加,降阶模型的计算精度逐渐提高,但当正交基数超过6时,降阶模型的计算精度趋于稳定,这意味着只增加正交阶数而不增加展开频点数量不能进一步提高降阶模型的计算精度;而在增加展开频点后,降阶模型的计算精度迅速提升,并且随着展开频点数目的增加,降阶模型的计算精度越来越高。当展开频点数目增加到5个时,降阶模型的最大动柔顺度误差最大值为1×10-7。因此,仅在单个展开频点生成多个标准正交基不能满足对降阶模型计算精度的要求,本文提出的自适应增加展开频点和标准正交基的策略可建立满足精度要求的降阶模型。

3.2 网格加筋筒壳的谐响应分析

本节利用网格加筋筒结构进一步验证自适应模型降阶方法的有效性。结构模型如图5所示,几何尺寸如下:圆柱筒直径D=5000 mm,长度L=16 000 mm,蒙皮厚度ts=5 mm,环筋数目21,纵筋数目30,筋条高度100 mm,筋条厚度tc=5 mm。蒙皮材料弹性模量E=71 GPa,筋条材料弹性模量为E=50 GPa,密度ρ=2700 kg/m3,泊松比μ=0.3,阻尼比为ξ=5%。边界条件为一端固支一端自由,在自由侧施加简谐激励,频率范围[0, 100 Hz],分析步长为1 Hz。

图5 正置正交网格加筋筒结构Fig.5 Orthogonal grid stiffened cylinder

选取初始标准正交基的阶数r=3、步长p=1,最大正交基数目Maxorder=20。根据图1给出的自适应模型降阶流程,在目标频段内确定了3个展开频点以及每个频点处所需标准正交基的数目,如表 2所示。将3个展开频点生成的标准正交基进行SVD正交化得到降阶模型的减缩基矩阵T。利用降阶模型计算结构的动柔顺度及相对误差如图6所示。

表2 自适应确定的展开频点和标准正交基阶数Table 2 Expansion frequency points and its related orders of orthonormal basis Frequency points

(a) Dynamic compliance

(b) Relative error图6 网格加筋筒结构的动柔顺度及误差Fig.6 Dynamic compliance and relative error of the orthogonal grid stiffened cylinder

图6中给出了选取结构前100阶模态时模态叠加法(Modal Superposition Method,MSM)和模态加速度法(Modal Acceleration Method,MAM)的计算结果。由图6可以看出,本文建立的降阶模型计算结果与直接频响法十分吻合,最大相对误差为1×10-7;模态叠加法和模态加速度法在低频段内计算精度能够满足要求,但在高频段内计算误差较大,模态叠加法的在高频段内预测误差最大值为421.97%,模态加速度法的预测误差最大值为46.63%。这是由于本文中的网格加筋筒壳结构是由两种材料组成的非比例阻尼结构,其阻尼矩阵不再满足模态正交化条件,导致模态叠加法和模态加速度法的求解精度降低。

表 3统计了直接频响法、模态叠加法法、模态加速度法和本文提出的自适应模型降阶法的计算时间,可以看出,自适应模型降阶法的计算时间为25.80 s,分别是直接频响法计算时间的0.78%,模态叠加法计算时间的37.52%,模态加速度法计算时间的35.47%。综上可以证明本文提出的自适应降阶模型表现出更优异的计算精度和效率。