航空发动机双转子系统动力学特性分析

Analysis of Dynamic Characterization of Dual-Rotor System in Aero-Engine

摘要: 本文针对五点支承的复杂双转子系统，考虑了中介轴承的非线性，通过有限单元法建立了带中介轴承的双转子系统动力学模型。采用Timoshenko梁单元模拟轴段，刚性圆盘模拟转盘，分别得到刚度矩阵、陀螺矩阵和质量矩阵，并将各单元矩阵组合得到系统总体矩阵和系统动力学微分方程。采用Newmark-β法求解系统动力学响应，分析各转盘动力学特性。研究了转速比和不平衡量对系统动力学响应的影响。研究结果表明，双转子系统出现两个共振区；改变转速比影响高压转子基频和对应幅值，导致盘1出现拍振现象；改变不平衡量会引起各盘主频率成分的幅值增加，尤其在共振区内。

Abstract: In this paper, for the complex dual rotor system with five-point support, the nonlinearity of the inter-shaft bearing is taken into account, and the dynamics model of the dual rotor system with intermediate bearing is established by the finite unit method. The Timoshenko beam unit is used to simulate the shaft segment and the rigid disk is used to simulate the rotor, and the stiffness matrix, gyro matrix and mass matrix are obtained respectively, and the overall matrix of the system and the differential equations of the system dynamics are obtained by combining the matrices of each unit. The Newmark-β method is used to solve the system dynamic response and analyze the dynamic characteristics of each turntable. The effects of the speed ratio and the amount of unevenness on the system dynamic response are investigated. The results show that two resonance zones appear in the dual rotor system; changing the rotational speed ratio affects the fundamental frequency of the high-pressure rotor and the corresponding amplitude, which leads to the beat vibration phenomenon of disk 1; changing the amount of unevenness causes an increase in the amplitude of the main frequency components of each disk, especially in the resonance zone.

1. 引言

在航空发动机系统中，双转子系统低压转子和高压转子组成，是影响发动机的性能和安全运行的关键部件。在大多数航空发动机中，低压转子和高压转子通过中介轴承连接。这种设计能提高发动机效率和性能但也会引起系统的耦合振动。由于不平衡力、非线性轴间轴承力和复杂的运行条件，内转子和外转子之间的耦合振动问题日益复杂，给发动机的稳定性和安全性带来了巨大挑战。

大部分双转子系统采用双线轴形式，即低压转子和高压转子两根轴分别以不同转速旋转，同时，高压转子和低压转子通过中介轴承连接，导致系统的动力学响应更加复杂。因此，对于复杂双转子系统动力学模型的建立和动力学特性研究就具有重要的意义。转子系统的建模普遍采用传递矩阵法和集中质量法，得益于计算机硬件的发展，有着高精度和高效率等特点的有限单元法逐渐成为主流方法[1]-[6]。孙传宗等[1]建立了双转子系统三维实体有限元模型，通过Craig-Bampton模态综合法实现模型维度的缩减。凌文辉等[7]建立了考虑轴间碰摩的复杂双转子系统动力学模型，通过模态综合法和Newmark隐式积分法结合求解系统非线性响应。Chen等[8]提出一种改进的谐波平衡–交变频域/时域方法，用于分析双转子系统–轴承–机匣非线性动力学特性。Yang等[9]建立了基座松动和多级涡轮机叶片–机匣摩擦的双转子–轴承–双机匣系统有限元模型，详细研究了支座松动和碰摩对双转子系统动力学特性的潜在影响。徐伟文等[10]建立了考虑挤压油膜阻尼器静偏心的转子系统的有限元模型，发现挤压油膜阻尼器静偏心对转子系统临界转速和对应幅值都有影响。Wang等[11]建立了带中间轴承的双转子–叶片–机匣系统的有限元模型，研究并发现叶片–机匣碰摩会对机匣和转子产生冲击载荷，导致振动幅值急剧增大。

2. 双转子系统动力学建模

本节介绍了某航空发动机双转子系统的模型，主要由风扇转子、低压涡轮转子、高压转子和支承系统组成，图1为双转子系统模型简化图。风扇转子和低压转子通过套齿联轴器连接，忽略联轴器的非线性，将低压风扇段转子和低压涡轮段转子刚性连接，简化成一根轴，双转子系统的具体参数见表1。低压转子与高压转子通过中介轴承连接，低压转子由部分空心轴和部分实心轴组成，采用1-1-1支承形式，而高压转子由空心轴组成，其支承形式为1-0-1。除中介支承外其余轴承均考虑为线弹簧–阻尼形式，支点的组合刚度是轴承刚度与弹性鼠笼刚度串联组合的。

Figure 1. Schematic diagram of a dual rotor system

图1. 双转子系统示意图

Table 1. Parameters of each shaft section

表1. 各轴段的参数

参数	编号和大小
长度/mm	L1	L2	L3	L4	L5	L6	L7	L8
	89	100	289	718	80	120	300	205
	L9	L10	L11	L12
	813	1185	79	1902
宽度/mm	D1	D2	D3	D4	D5	D6	D7	D8
	56	70	45	70	65	50	56	70

2.1. 中介轴承非线性力模型

双转子系统中一般用中介轴承来连接高压转子和低压转子，通常为圆柱滚子轴承，其主要由滚动体、保持架和轴承内外圈组成，中介轴承示意图如图2所示。轴承的内外圈与滚珠之间的接触形式是点接触。图中虚线位置为内外圈未变形的位置，轴颈发生位移导致形变产生时，根据赫兹接触理论，滚子和滚道之间的接触变形会产生一个非线性恢复力。

中介轴承内圈与低压转子连接，轴承外圈与高压转子连接，滚子与滚道接触时在x方向和y方向上的非线性恢复力为$F_{bx}$ 和$F_{by}$ 的表达式为：

$\{\begin{array}{c}{F}_{bx}={\displaystyle \sum }_{j=1}^{N_b}{C}_b{\left(x\text{cos}{\beta }_j+y\text{sin}{\beta }_j-G_0\right)}^{10/9}H\left(x\text{cos}{\beta }_j+y\text{sin}{\beta }_j-G_0\right)\text{cos}{\beta }_j\\ F_{by}={\displaystyle \sum }_{j=1}^{N_b}{C}_b{\left(x\text{cos}{\beta }_j+y\text{sin}{\beta }_j-G_0\right)}^{10/9}H\left(x\text{cos}{\beta }_j+y\text{sin}{\beta }_j-G_0\right)\text{sin}{\beta }_j\end{array}$ (1)

$\{\begin{array}{l}{\beta }_j={\omega }_{cage}t+\frac{2\pi }{N_b\left(j-1\right)},j=1,2,\cdots ,N_b\hfill \\ {\omega }_{cage}=\frac{{\omega }_L{r}_i+{\omega }_H{r}_o}{\left(r_i+r_o\right)}\hfill \end{array}$ (2)

$H\left({\delta }_j\right)=\{\begin{array}{c}1,{\delta }_j>0\\ 0,{\delta }_j<0\end{array},{\delta }_j=x\cos {\beta }_j+y\sin {\beta }_j-G_0$ (3)

式中，$N_b$ 表示滚子的数量，共23个滚子；$C_b$ 是赫兹接触刚度，与接触材料和形状有关，取其刚度值为1.33 × 10⁹ N/m；${\beta }_j$ 是第j个滚子的转角位置；${\omega }_{cage}$ 是轴承保持架的角速度；$r_o$ 和$r_i$ 是轴承外圈半径和内圈半径，大小分别为50 mm和32.5 mm；G₀ = 0.1 μm表示轴承径向间隙。${\delta }_j$ 是第j个滚子与滚道的接触变形，$H\left({\delta }_j\right)$ 为Heaviside函数。${\omega }_L$ 是低压转子的转速，${\omega }_H$ 是高压转子的转速。值得注意的是，x和y表示滚子在径向上的相对振动位移，$x=x_i-x_o$ ，$y=y_i-y_o$ ，其中$x_i,x_o,y_i$ 和$y_o$ 分别表示中介轴承内外圈对应轴颈节点在x方向和y方向上的位移分量。

Figure 2. Schematic diagram of the intershaft bearing model

图2. 中介轴承模型示意图

2.2. 双转子系统有限元建模

图3为双转子系统离散化示意图，将低压转子离散成9个单元共10个节点，将高压转子分为5个单元共6个节点，故双转子系统共16个节点。转轴采用考虑转动惯量、陀螺效应及剪切变形的Timoshenko梁单元进行建模。图4(a)为梁单元示意图，每个梁单元包括两个节点，每个节点考虑4个自由度，即x_A、y_B、θ_xA、θ_yA和x_B、y_B、θ_xB、θ_yB表示节点A和节点B的两个平动自由度和两个旋转自由度，故任意梁单元的位移矢量为：

$q_e=\left[\begin{array}{cccccccc}{x}_A& y_A& {\theta }_{xA}& {\theta }_{yA}& x_B& y_B& {\theta }_{xB}& {\theta }_{yB}\end{array}\right]$ (4)

Figure 3. Schematic of discretization of dual rotor system

图3. 双转子系统离散化示意图

由拉格朗日方程推导得到轴单元在固定坐标系中的运动方程为：

$\left(M_e^t+M_e^r\right){\ddot{q}}_e-\omega G_e{\dot{q}}_e+K_e{q}_e=Q_e$ (5)

式中，$M_e^t$ 和$M_e^r$ 为轴单元的质量矩阵；$G_e$ 和$K_e$ 分别为轴单元的阻尼矩阵和刚度矩阵；$Q_e$ 表示外力。因为每个节点考虑四个自由度，一个梁单元有两个节点，故所有的矩阵都是8阶方阵，且都是关于对角线对称或反对称矩阵。

Figure 4. Shaft unit and disk unit finite element model

图4. 轴单元和盘单元有限元模型

各矩阵的具体表达式为：

质量矩阵：

$M_e^t=\frac{{\rho }_{l}l}{{\left(1+{\phi }_s\right)}^2}\left[\begin{array}{cccccccc}{M}_{t1}& & & & & & & \\ 0& M_{t1}& & & 对& 称& & \\ 0& -M_{t4}& M_{t2}& & & & & \\ M_{t4}& 0& 0& M_{t2}& & & & \\ M_{t3}& 0& 0& M_{t5}& M_{t1}& & & \\ 0& M_{t3}& -M_{t5}& 0& 0& M_{t1}& & \\ 0& M_{t5}& M_{t6}& 0& 0& M_{t4}& M_{t2}& \\ -M_{t5}& 0& 0& M_{t6}& -M_{t4}& 0& 0& M_{t2}\end{array}\right]$ (6)

$\begin{array}{l}{M}_{t1}=\frac{13}{35}+\frac{7}{10}{\phi }_s+\frac{1}{3}{\phi }_s^2，M_{t2}=\left(\frac{1}{105}+\frac{1}{60}{\phi }_s+\frac{1}{120}{\phi }_s^2\right)l^2，M_{t3}=\frac{9}{10}+\frac{3}{10}{\phi }_s+\frac{1}{6}{\phi }_s^2，\\ M_{t4}=\left(\frac{11}{210}+\frac{11}{120}{\phi }_s+\frac{1}{24}{\phi }_s^2\right)l，M_{t5}=\left(\frac{13}{420}+\frac{3}{40}{\phi }_s+\frac{1}{24}{\phi }_s^2\right)l，M_{t6}=-\left(\frac{1}{140}+\frac{1}{60}{\phi }_s+\frac{1}{120}{\phi }_s^2\right)l^2\end{array}$ (7)

$M_e^r=\frac{{\rho }_{l}I}{l{\left(1+{\phi }_s\right)}^{2}A}\left[\begin{array}{cccccccc}{M}_{r1}& & & & & & & \\ 0& M_{r1}& & & 对& 称& & \\ 0& -M_{r4}& M_{r2}& & & & & \\ M_{r4}& 0& 0& M_{r2}& & & & \\ -M_{r1}& 0& 0& -M_{r4}& M_{r1}& & & \\ 0& -M_{r1}& M_{r4}& 0& 0& M_{r1}& & \\ 0& -M_{r4}& M_{r3}& 0& 0& M_{r4}& M_{r2}& \\ M_{r4}& 0& 0& M_{r3}& -M_{r4}& 0& 0& M_{r2}\end{array}\right]$ (8)

$M_{r1}=\frac{6}{5}，M_{r2}=\left(\frac{2}{15}+\frac{1}{6}{\phi }_s+\frac{1}{3}{\phi }_s^2\right)l^2，M_{r3}=\left(-\frac{1}{30}-\frac{1}{6}{\phi }_s+\frac{1}{6}{\phi }_s^2\right)l^2，M_{r4}=\left(\frac{1}{10}-\frac{1}{2}{\phi }_s\right)l$ (9)

陀螺效应矩阵：

$G_e=\frac{{\rho }_{l}I}{15l{\left(1+{\phi }_s\right)}^{2}A}\left[\begin{array}{cccccccc}0& & & & & & & \\ G_1& 0& & & 反& 对& 称& \\ -G_2& 0& 0& & & & & \\ 0& -G_2& G_4& 0& & & & \\ 0& G_1& -G_2& 0& 0& & & \\ -G_1& 0& 0& -G_2& G_1& 0& & \\ -G_2& 0& 0& G_3& G_2& 0& 0& \\ 0& -G_2& -G_3& 0& 0& G_2& G_4& 0\end{array}\right]$ (10)

$M_{r1}=\frac{6}{5}，M_{r2}=\left(\frac{2}{15}+\frac{1}{6}{\phi }_s+\frac{1}{3}{\phi }_s^2\right)l^2，M_{r3}=\left(-\frac{1}{30}-\frac{1}{6}{\phi }_s+\frac{1}{6}{\phi }_s^2\right)l^2，M_{r4}=\left(\frac{1}{10}-\frac{1}{2}{\phi }_s\right)l$ (11)

刚度矩阵：

$K_e=\frac{EI}{l^3\left(1+{\phi }_s\right)}\left[\begin{array}{cccccccc}{K}_1& & & & & & & \\ 0& K_1& & & 对& 称& & \\ 0& -K_4& K_2& & & & & \\ K_4& 0& 0& K_2& & & & \\ -K_1& 0& 0& -K_4& K_1& & & \\ 0& -K_1& K_4& 0& 0& K_1& & \\ 0& -K_4& K_3& 0& 0& K_4& K_2& \\ K_4& 0& 0& K_3& -K_4& 0& 0& K_2\end{array}\right]$ (12)

$K_1=12，K_2=\left(4+{\phi }_s\right)l^2，K_3=\left(2-{\phi }_s\right)l^2，K_4=6l$ (13)

其中，ρ_l为材料密度；l为单元长度；I为截面惯性矩；A为单元横截面面积；E为材料弹性模量；${\phi }_s=12EI/G{A}_s{l}^2$ ，A_s为有效抗剪面积，与轴截面形状有关。

图4(b)为刚性圆盘单元示意图，x_d、y_d、θ_xd和θ_yd表示盘的两个平动振动位移和两个转动振动位移，并以集中质量叠加到对应节点上。则刚性圆盘的运动方程为：

$M_d{\ddot{q}}_d+\omega G_d{\dot{q}}_d=Q_d$ (14)

其中，M_d和G_d分别为圆盘质量矩阵和陀螺矩阵；Q_d为外力。各矩阵表达式为：

$M_d=\left[\begin{array}{cccc}{m}_d& & & \\ & m_d& & \\ & & J_d& \\ & & & J_d\end{array}\right]，G_d=\left[\begin{array}{cccc}0& & & \\ & 0& & \\ & & 0& -J_p\\ & & J_p& 0\end{array}\right]，q_d=\left[\begin{array}{c}{x}_d\\ y_d\\ {\theta }_{xd}\\ {\theta }_{yd}\end{array}\right]$ (15)

式中，m_d、J_d和J_p分别为圆盘的质量、直径转动惯量和极转动惯量。

将双转子系统各个梁单元矩阵对应组装，即将节点、刚性圆盘和轴承各个单位矩阵叠加得到总体矩阵形式。以刚度矩阵为例，系统总体刚度矩阵如图5所示。故双转子系统总体动力学方程为：

$M_{sys}\ddot{q}+\left(C_{sys}-\omega G_{sys}\right)\dot{q}+K_{sys}q=F_{sys}$ (16)

式中，M_sys，C_sys，G_sys和K_sys分别为系统总质量矩阵、阻尼矩阵、陀螺矩阵和刚度矩阵；q为广义位移矢量；F_sys为外激励矢量，由中介轴承非线性力矢量F_b和各个转盘的偏心力矢量F_u组成，具体表示为：

$F_b={\left[0,0,0,0,\cdots ,F_{bx},F_{by},0,0,\cdots ;-F_{bx},-F_{by},0;0;0;0;0;0\right]}^{{\rm T}}$ (17)

$\begin{array}{c}{F}_u=[0,0,0,0,m_{d1}{e}_1{\omega }_L^2\cos \left({\omega }_{L}t\right),m_{d1}{e}_1{\omega }_L^2\sin \left({\omega }_{L}t\right),0,0,\cdots ,\\ m_{d4}{e}_4{\omega }_H^2\cos \left({\omega }_{H}t\right),m_{d4}{e}_4{\omega }_H^2\sin \left({\omega }_{H}t\right),0,0,\cdots ,\\ {m_{d7}{e}_7{\omega }_H^2\cos \left({\omega }_{H}t\right),m_{d7}{e}_7{\omega }_H^2\sin \left({\omega }_{H}t\right),0,0,0,0,0,0]}^{\text{T}}\end{array}$ (18)

式中，$m_{di}$ 和$e_i\left(i=1,2,\cdots ,7\right)$ 代表第i个圆盘的集中质量和偏移量。$C_{sys}=\alpha M_{sys}+\beta K_{sys}$ 是瑞利阻尼，α和β是比例系数，具体表示为：

$\{\begin{array}{c}\alpha =\frac{2\left({\xi }_2/{\omega }_2-{\xi }_1/{\omega }_1\right)}{\left(1/{\omega }_1-1/{\omega }_2\right)}\\ \beta =\frac{2\left({\xi }_2{\omega }_2-{\xi }_1{\omega }_1\right)}{\left({\omega }_1^2-{\omega }_2^2\right)}\end{array}$ (19)

其中，ξ₁和ξ₂是阻尼系数，分别取0.02和0.04；ω₁和ω₂是系统的前两阶固有频率。

Figure 5. Schematic diagram of the total system stiffness matrix assembly

图5. 系统总体刚度矩阵组合示意图

3. 仿真结果与分析

在实际工作中，高压和低压转子分别旋转并保持一定的速比，将高低压转子设定为反比，速比为1.2。采用Newmark-β法来求解双转子系统非线性动力学响应，转速范围设定为1020~6660 r/min，转速间隔为120 r/min。研究转速对系统振动响应的影响，分别以双转子系统上每个转盘作为研究对象，取偏心量为256 g∙mm，得到各个转盘的幅频响应，如图6所示。

仿真结果可知双转子系统出现两个共振区，分别对应着系统前两阶固有频率。低压转子盘1和盘2的幅频响应的结果相似，都是位于低压转子前端，引起系统二阶共振，其中盘1在二阶共振区的振动更突出。低压转子盘3和高压转子盘6和盘7位于系统末端，引发系统一阶共振，盘7在一阶共振区的振动更明显。当旋转机械的转速接近系统的固有频率时，激励频率与固有频率重合，导致系统进入共振状态。随着外部激励能量的持续输入(比如不平衡激励力)，而系统阻尼无法及时耗散，能量积累使振动幅值显著增大，导致系统出现共振区。在双转子系统中，低压转子和高压转子通过中介轴承连接并相互耦合，这种强耦合效应会引起系统的共振现象，导致系统出现两个共振区。

Figure 6. The amplitude-frequency response curve of each disk of dual rotor system

图6. 双转子系统各盘幅频响应曲线

3.1. 转速比对系统动力学响应的影响

本节分析转速比η对系统动力学响应的影响，前文中设置高低压转子为反向旋转，这样能互相抵消陀螺力矩，提高系统不平衡振动响应的幅值。取速比η范围为−1~−1.4，其他仿真参数不变。图7和图8分别表示盘1在4380 r/min和盘7在1980 r/min时的时域轨迹、轴心轨迹和频谱图。

由图7和图8可知，速比直接影响着系统运行周期。当η为−1时，这种情况下高低压转子转速相等，在频谱图中的组合频也都是倍频，故只显示出单个低压转子频率f₁，转盘的轴心轨迹呈单个圆形或椭圆形。当η为−1.1~−1.4时，其中速比小数位为奇数时，轴心轨迹的花瓣相对于小数位为偶数时更加密集，系统的周期数变大。由于改变速比时固定内转子转速不变，改变的是外转子转速，导致外转子基频f₂增加，但随着速比的增加，f₂的幅值随之变小。改变转速比会直接改变高压转子的激励频率，影响两个转子间的能量传递与分配，当转速比接近−1.2时，能量在两个转子间高效传递，导致该速比下振动幅值最大。盘7的幅值明显表现出先增加后减小，除了上述原因，还因为改变转速比还会改变系统的临界转速。临界转速是系统固有频率与激励频率重合时的转速，据上文分析高压转子是引起系统一阶共振主要原因，当高压转子转速为2400 r/min时才会引起系统一阶共振。图8是盘7在低压转子为1980 r/min得到的振动响应，而改变速比会使得高压转子转速改变，使得以低压转子转速为研究对象时会发生一阶共振区偏移。当速比绝对值大于1.2时，共振频率减小，共振区偏向左移；当速比绝对值小于1.2时，共振频率增大，共振区发生右移。

Figure 7. Dynamic response of disk 1 at ω_L = 4380 r/min

图7. 盘1在ω_L = 4380 r/min时的动力学响应

Figure 8. Dynamic response of disk 7 at ω_L = 1980 r/min

图8. 盘7在ω_L = 1980 r/min时的动力学响应

3.2. 偏心量对系统动力学响应的影响

本节主要研究双转子系统非线性不平衡响应，以不平衡量U为参数，选取不平衡量U为156 g∙mm、256 g∙mm和356 g∙mm，其他参数保持不变。不同不平衡量下盘1和盘7的幅频响应如图9所示。观察发现随着不平衡量的增加，转盘的激励力增大，导致系统振动位移增大，尤其在共振区的振幅明显增大。所以偏心量的改变会引起系统不平衡激励力的变化。随着偏心量的增加，盘1的激励力增大导致二阶共振响应增强，进而通过中介轴承的强耦合效应导致盘7在二阶共振位移增大。同理可知，盘7在一阶共振区的振幅增大导致盘1在一阶共振区出现振幅增加现，故偏心量的增加会使得模态耦合效应增强。

Figure 9. The amplitude-frequency response of disk 1 and disk 7 under different unbalanced quantity

图9. 不同不平衡量下盘1和盘7的幅频响应

图10和图11分别表示盘1在4380 r/min和盘7在1980 r/min下的时域轨迹、轴心轨迹和频谱图。随着不平衡量的增加，时域图和轴心轨迹图出现的规律与前面分析的规律一致。观察频谱图可知，盘1在4380 r/min时的振动响应以低压转子基频为主要频率成分，此时低压转子通过中介轴承的强耦合作用使得高压转子在二阶共振区振动响应增强；而盘7在1980 r/min时的振动响应以高压转子基频为主要频率成分，其振动通过中介轴承迫使低压转子发生强振动。随着不平衡量U的增加，主要增加转盘主要频率成分的幅值。随着主要频率成分幅值的增加，导致两个频率之间的差距过大，减少系统出现拍振现象。但偏心量过大使得系统振幅增加会降低系统的稳定性，尤其是在高转速下。因此，在设计和运行中，需要通过动平衡技术减小偏心量，优化系统的动力学特性，以确保系统的稳定性和可靠性。

Figure 10. Dynamic response of disk 1 at ω = 4380 r/min

图10. 盘1在ω_L = 4380 r/min时的动力学响应

Figure 11. Dynamic response of disk 7 at ω = 1980 r/min

图11. 盘7在ω_L = 1980 r/min时的动力学响应

4. 总结

本文建立了考虑中介轴承的双转子系统动力学有限元模型，采用数值积分的方法求解系统动力学响应，研究了系统各个转盘、转速比和不平衡量对系统非线性动力学响应的影响。主要结论如下：

随着转速的增加，双转子系统出现两个共振区，分别对应着系统前两阶固有频率，低压转子盘1在二阶共振区的振幅最大，而盘7在一阶共振区的振幅更突出。

改变转速比会使得高压转子基频频率增大，转速比小位数为奇数时，系统轴心轨迹花瓣更密集，系统周期数变大。同时改变速比会影响转子间的模态耦合，甚至导致系统的临界转速发生偏移。

不平衡量的增加会使得转盘激励力变大，导致系统振动位移增大，转盘的主频率成分的幅值增大，特别是在共振区内更明显。不平衡量的增加也会使转子间通过中介轴承的耦合效应增强，导致系统非线性效应增强。

致 谢

本研究工作由“机械系统与振动全国重点实验室开放基金课题资助MSV202507”资助。