在航空发动机的研制过程中，人们为了提高其推重比及其结构效率，不断缩小发动机的转、静子间隙，与此同时，转子质量的不平衡、转子与静子的不对中、高温膨胀等因素也会导致转子系统振动大于转、静子间隙，进而引发大量的碰摩问题。根据资料统计，由碰摩故障引发的航空发动机转子系统运动失稳、轴系破坏、叶片折断直至机械失效等事故日益突出，在近年来我国航空发动机的重大事故中，有30%是由碰摩故障而引起的直接或间接振动故障所导致^[1-2]，该类型故障的特殊性和复杂性也日益受到国内外学者的关注。

国外，Childs^[3]研究了碰摩转子的拟周期现象，研究发现在碰摩转子中，倍周期的分叉进入到混动区域是必然的，同时还发现了稳定区域范围会随着系统的参数变化而发生改变。Muszynska^[4]基于动量守恒定律研究了系统的碰摩问题，并建立了一个较为完备的关于约束微分系统的碰摩力学模型。同时，在上述模型中引入了弹性恢复系统，用来表示碰摩时能量的损失，提出了对模型进行完善的方法^[5]。Beatty^[6]提出了一种通过分段线形刚度表示碰摩力的数学模型，并基于谐波响应实验数据，确定了良性接触与故障初期不稳定接触之间的界限。Adams等^[7]提出了一种数值仿真分析方法，并研究了转子系统的转速、不平衡量和径向间隙对碰摩系统动力学特性的影响和规律。Piccoli等^[8]提出了一种识别混沌运动的测试方法，并测试获得了庞加莱图、混沌相图、李雅普诺夫指数等参数。同时，还将转子和定子的偏心旋转影响引入到Jeffcott转子系统的碰摩模型中，研究了混沌运动对转子系统运动稳定性的影响。Groll等^[9]提出了一种基于谐波平衡法的数值算法，用以计算在周期激励下一个非线性转/定子接触系统的周期响应。Popprath等^[10]描述了一种与定子有间歇性接触的Jeffcott转子，其中的定子模型可作为一个附加的振动系统，并与转子模型通过非线性接触力相互作用。Banakh等^[11]研究了转子和浮动密封圈之间的振动和冲击影响，还考虑了转子和密封圈间隙之间的水动力以及密封圈和外壳之间的干摩擦。

国内，褚福磊等^[12]讨论了转子系统的转速及不平衡量的变化对碰摩转子系统的振动特性的影响，并建立了一个由油膜轴承支撑的碰摩转子动力学模型。王德友^[13]总结了旋转机械中，转、静子碰摩时的振动特性，并且指出碰摩故障是其他故障(如不平衡、偏心、干摩擦等)造成的二次效应，并且碰摩故障还会使结构的振动幅值大幅增加。戴兴建等^[14]考虑了限位器的碰摩，建立了碰摩转子系统的动力学方程，模型中对碰摩力采用了分段线弹性以及库伦摩擦模型的描述方法，对系统中的双频进动、局部碰摩、整圈碰摩以及拟整圈碰摩时的振动特性进行了研究。单颖春等^[15]针对现有的模型碰摩力模型存在的问题，首先应用ANSYS软件中的接触分析，计算了转静子间的非线性接触刚度，获得非线性碰摩力模型。同时，还利用时域波形、频谱分析和小波变换对响应进行了分析，得出相关的碰摩故障特征。何田等^[15]建立了双转子碰摩系统的动力学方程，利用分叉图、相图、Poincaré图等从时域研究了高、低速转子随转速变化时显示的周期、拟周期和混沌等复杂动力学现象及其变化规律。

虽然目前对航空发动机碰摩故障问题进行了大量研究，并取得了阶段性的研究成果，但主要是集中在单转子碰摩问题上。本文以某型号航空发动机的双转子模型实验台为研究对象，在合理简化的基础上，基于ADAMS多体动力学仿真软件，建立了双转子碰摩系统的动力学模型，同时对其动力学特性进行了仿真研究。

1 双转子碰摩系统ADAMS模型 1.1 多体系统的动力学原理

在ADAMS中，可根据多体动力学相关理论，建立系统的拉格朗日运动方程。对于每个刚体，可列出6个广义坐标，且根据带乘子的拉格朗日方程，获得相应的约束方程^[17]：

$ \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial \mathit{\boldsymbol{K}}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}}} \right) - \left( {\frac{{\partial \mathit{\boldsymbol{K}}}}{{\partial {\mathit{\boldsymbol{q}}_j}}}} \right) + \sum\limits_{i = 1}^n {\frac{{\partial {\mathit{\boldsymbol{\psi }}_i}}}{{\partial {\mathit{\boldsymbol{q}}_j}}}{\lambda _i}} = {\mathit{\boldsymbol{F}}_j} $

(1)

$ {\mathit{\boldsymbol{\psi }}_i} = 0\;\;\;\left( {i = 1,2, \cdots ,m} \right) $

(2)

式中：K为系统动能，q_j描述系统的广义坐标，ψ_i为统的约束方程，F_j为在广义坐标方向的广义力，λ_j为m×1的拉格朗日乘子列阵。

物体B_i(i=1, 2, …, N)的动能可表示为

$ {\mathit{\boldsymbol{T}}_i} = \sum\limits_{k = 1}^l {\frac{1}{2}m_i^k{{\mathit{\boldsymbol{\dot r}}}^{{k^{\rm{T}}}}}\mathit{\boldsymbol{\dot r}}_i^k} $

(3)

式中：l为物体B_i的节点数；m_i^k为节点k的质量；r_i^k为节点k在惯性基中的坐标阵，r_i^k=B_i^kv_i。

ADAMS中柔性体模块是基于弹性小变形假设建立的。其基本思想是赋予柔性体一个模态集，采用模态展开法，利用模态坐标来表示弹性位移。通过计算每一时刻物体的弹性位移来描述其变形运动。设柔性体的模态坐标是q_f={q₁, q₂, …, q_n}^T，则柔性体的广义坐标可表示为

$ {\mathit{\boldsymbol{r}}_i} = \mathit{\boldsymbol{x}} + \mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{S}}_i} + {\mathit{\boldsymbol{\varphi }}_i}{\mathit{\boldsymbol{q}}_f}} \right) $

(4)

式中：x为从整体坐标系原点到局部坐标系原点的位置矢量，A为局部坐标系相对整体坐标系的方向余弦矩阵，S_i为i点在未变形时在局部坐标系的位置，φ_i为点i的移动自由度的模态矩阵子块。

物体B_i(刚-柔系统)上任一点i在整体坐标系的位置向量可以表示为

$ \mathit{\boldsymbol{\xi }} = {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{x}}&\mathit{\boldsymbol{\psi }}&{{\mathit{\boldsymbol{q}}_f}} \end{array}} \right]^{\rm{T}}} = {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{x}}&\mathit{\boldsymbol{y}}&\mathit{\boldsymbol{z}}&\mathit{\boldsymbol{\psi }}&\mathit{\boldsymbol{\varphi }}&{{\mathit{\boldsymbol{q}}_{i,j = 1,n}}} \end{array}} \right]^{\rm{T}}} $

(5)

刚-柔系统的动能可以表示为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{E}} = \frac{1}{2}\int {\rho \mathit{\boldsymbol{\dot r}}_i^{\rm{T}}\mathit{\boldsymbol{\dot r}}{\rm{d}}V} = \frac{1}{2}\sum\limits_{i = 1}^N {\left( {{m_i}{{\mathit{\boldsymbol{\dot r}}}^{\rm{T}}}\mathit{\boldsymbol{\dot r}} + \mathit{\boldsymbol{\omega }}_i^{\rm{T}}{\mathit{\boldsymbol{I}}_i}{\mathit{\boldsymbol{\omega }}_i}} \right)} = }\\ {\frac{1}{2}{{\mathit{\boldsymbol{\dot \xi }}}^{\rm{T}}}\mathit{\boldsymbol{M}}\left( \mathit{\boldsymbol{\xi }} \right)\mathit{\boldsymbol{\dot \xi }}} \end{array} $

(6)

式中：m_i为结点i的模态质量，I_i为结点i的模态惯量，ω为物体坐标系的角速度向量。

该系统的势能包括两部分，其一为物体B_i所固有的重力势能和物体弹性变形势能，另外与系统的有势力有关，如线弹簧力元的势能等。关于物体B_i(i=1, 2, …, N)固有的势能，记为V_i，它包括重力势能和物体弹性变形势能，分别记为V_i^g和V_i^u，即有

$ {\mathit{\boldsymbol{V}}_i} = \mathit{\boldsymbol{V}}_i^g + \mathit{\boldsymbol{V}}_i^u $

(7)

令重力势能的参考点为惯性基的坐标原点，物体B_i的重力势能为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{V}}_i^g = - \sum\limits_{k = 1}^l {m_i^k\mathit{\boldsymbol{g}}_i^{\rm{T}}\mathit{\boldsymbol{r}}_i^k} = }\\ { - \mathit{\boldsymbol{g}}_i^{\rm{T}}\sum\limits_{k = 1}^l {m_i^k\left( {{\mathit{\boldsymbol{r}}_i} + {\mathit{\boldsymbol{A}}^i}\rho {'}_{i0}^k + {\mathit{\boldsymbol{A}}^i}\mathit{\boldsymbol{ \boldsymbol{\varPhi} '}}_i^k{\mathit{\boldsymbol{a}}_i}} \right)} } \end{array} $

(8)

式中g为重力在惯性基中的坐标阵。

令K_aⁱ为物体B_i的模态刚度阵，则其弹性变形势能V_i^u为

$ \mathit{\boldsymbol{V}}_i^u = \frac{1}{2}{\mathit{\boldsymbol{a}}^{\rm{T}}}\mathit{\boldsymbol{K}}_a^i\mathit{\boldsymbol{a}} $

(9)

如果系统中有势力力元存在，则系统中存在这些力元的势能，记力元势能为V_k^e。如对于弹簧刚度为k的线弹簧力元，该力元的势能为

$ \mathit{\boldsymbol{V}}_k^e = \frac{1}{2}k{\left( {\left| \mathit{\boldsymbol{d}} \right| - \left| {{\mathit{\boldsymbol{d}}_0}} \right|} \right)^2} = \frac{1}{2}k{\left( {1 - \frac{{\left| {{\mathit{\boldsymbol{d}}_0}} \right|}}{{\sqrt {{\mathit{\boldsymbol{d}}^{\rm{T}}}\mathit{\boldsymbol{d}}} }}} \right)^2}{\mathit{\boldsymbol{d}}^{\rm{T}}}\mathit{\boldsymbol{d}} $

(10)

式中：|d|－|d₀|为力元的变形，|d₀|为弹簧原长。

由式(1)、(2)可以获得刚/柔混合体的动力学微分-代数方程。对于这种代数-微分方程的求解方法有多种，如可将二阶微分方程降阶为一阶微分方程来求解，或直接对二阶微分方程进行积分求解，ADAMS就采用了前一种方法。

1.2 双转子碰摩系统ADAMS模型的建立

航空发动机转子结构中，多数为双转子结构，因此在转、静件发生碰摩时，由于碰摩力的非线性作用，其振动响应的特征与单转子系统对应的振动响应特征有着非常明显的差别。针对如图 1所示的双转子碰摩试验系统，建立了双转子碰摩的多体动力学模型，其中，系统的弹性模量E为2.09×10¹¹Pa，泊松比υ为0.295，密度ρ为7.87×10³ kg/m³，转子系统的基本参数如表 1所示。

进行转子系统建模时，采取如下简化处理：1)两端弹性支承结构相同，将其均视为线性，支承刚度由轴承刚度和鼠笼结构刚度串联得到，分别在横向水平和垂直方向上，采用“spring”模拟支承，刚度、阻尼设置同表 1计算值；2)转动约束为理想约束；3)转轴视为刚体，不考虑转轴变形对动力学特性的影响；4)通过在高、低压转盘上各附加一质量块来模拟双转子不平衡；5)在距离低压转盘一定间隙处设置碰摩板，通过调整间隙、碰摩参数来模拟不同碰摩工况。

在建立转子系统几何模型的基础上，定义材料属性、合理设置约束与载荷，具体建模过程为：1)将转轴、转盘的材料参数根据材料属性进行设置；2)低压转子系统两端用“spring”单元来模拟轴承，施加在两端支承位置与大地之间；在输入轴支承处施加绕x轴旋转的点驱动；3)高压转子系统两端用“spring”单元来模拟轴承，一端施加在支承位置与大地之间；另一端施加在支承位置与转轴之间；4)动力学模型校对，设置好约束与载荷后，利用“Model Verify”对模型约束检查，以确保模型的正确性。最终建立的ADAMS双转子碰磨多体动力学模型如图 2所示。

2 双转子碰摩系统ADAMS模型分析

根据1.2节建立的模型，模拟实际工况进行仿真，初始不平衡量1 500 g·mm，转速3 000 r/min (50 Hz)，初始碰摩参数如表 2所示。通过调整碰摩参数，模拟不同的碰摩状态，进行碰摩响应仿真分析。

首先，通过选取不同碰摩刚度10⁴、10⁶和10⁸ N/m模拟双转子轻微、较重和严重碰摩工况，其对应碰摩阻尼分别为10、100和1 000 N·s/m。获得三种状态下低压转盘中心的时域、频域和轨迹分别如图 3~5所示。

由图 3~5可以看出，随着碰摩刚度以及阻尼系数的增加，振动响应的频谱变得越来越复杂，出现了高次谐波、组合谐波和次谐波成分，倍频频率成分越来越丰富。当碰摩刚度以及阻尼较大时，碰摩对于系统的振动响应影响越大，进而加剧碰摩。

为了研究不同刚度和阻尼对转子系统碰摩力的影响，图 6给出了不同碰摩刚度条件下的碰摩力。对碰摩力的变化趋势进行分析可知，随着碰摩刚度的增大，碰摩力将随着增大，进而加剧了碰摩程度。

3 结论

1) 通过对碰摩转子系统多体动力学模型进行仿真分析可知，在转静件碰摩从轻到重的整个过程中，其振动频谱成分与碰摩参数有着密切的关系。

2) 随着碰摩刚度以及阻尼系数的增加，振动响应的频谱变得越来越复杂，出现了高次谐波、组合谐波和次谐波成分。

3) 当碰摩刚度和阻尼系数较大时，随着碰摩刚度的增大，碰摩力将随着增大，进而加剧碰摩。