近年来，同/反转双转子结构、中介轴承技术以及挤压油膜阻尼器(Squeeze film damper，SFD)在航空发动机上得到了广泛应用，如RB211，GE120，F119等。反向旋转双转子结构可以有效地改善燃油消耗、提高推重比以及降低陀螺力矩^[1]；中介轴承可以有效地减轻发动机的质量；SFD可以有效抑制、隔离转子振动。因此，针对含SFD、中介轴承的同/反转双转子系统开展振动响应特性研究对航空发动机的设计和制造具有重要意义。

国内外学者在多转子系统振动特性方面做了大量研究。胡绚等^[2]对中介轴承单独进行了建模和分析，对同向和反向旋转时中介轴承的寿命进行了对比分析，但是没有与双转子系统进行耦合分析。罗贵火等^[3-8]针对含滚动轴承的双转子系统进行了同/反转转子系统建模方法和加、减速响应特性的研究。冯国全等^[9-10]针对含中介轴承的反向旋转双转子系统，将中介轴承简化为弹簧和阻尼，进行了临界转速计算方法及支承不对中故障的研究。张大义等^[11]结合有限元法和振型筛选法，提出了适用于工程应用的临界转速求解方法。周海仑、陈果等^[12-13]以双转子航空发动机为研究对象，考虑中介轴承及SFD的非线性力，建立了双转子系统整机振动模型并进行了耦合振动分析。Chiang, Guskov^[14-15]将中介轴承简化为弹簧和阻尼，研究了中介轴承-双转子系统的动力特性。Philip等^[16-19]以含SFD的航空发动机多转子系统为研究对象进行了非线性不平衡响应计算方法的研究。虽然国内外众多学者在多转子系统振动特性方面进行了大量研究，但关于转速比影响的研究鲜有报道。

因此，本文在考虑中介轴承和SFD非线性力的情况下，建立双转子系统非线性动力学模型。利用该仿真模型对双转子系统的非线性响应特性进行分析，研究转速比对双转子系统非线性动力特性的影响并通过试验对仿真结果进行验证。

1 转子系统模型 1.1 动力学模型及求解方法

SFD-中介轴承-双转子系统结构简图见图 1。内、外转子分别模拟航空发动机低压转子和高压转子；盘1和盘3分别模拟低压和高压压气机，盘2和盘4分别模拟低压和高压涡轮；内外转子通过中介轴承Ⅳ连接；支承Ⅰ，Ⅱ，Ⅲ采用鼠笼+滚动轴承+挤压油膜阻尼器的形式。本文的研究中，认为轴承Ⅰ，Ⅱ和Ⅲ为刚性，考虑弹性支承的线性刚度、SFD和中介轴承的非线性力。转子系统转速比定义为$ \lambda \text{=}\frac{{{\mathit{\omega }}_{\text{2}}}}{{{\mathit{\omega }}_{1}}}$，其中ω₂为外转子转速，ω₁为内转子转速。λ＜0表示内外转子反向旋转；λ＞0表示内外转子同向旋转；忽略转子系统的外部阻尼和材料阻尼，考虑各圆盘处的陀螺力矩。

内外转子轴所用材料的密度为7 810 kg/m³，弹性模量为196 GPa，泊松比为0.298。Newmark-β法的计算参数α=0.25，β=0.5。转子系统几何尺寸见图 1；各弹性支承刚度见表 1；中介轴承参数见表 2；挤压油膜阻尼器参数见表 3；盘的质量属性及不平衡量配置见表 4；主模态矩阵Φ_k中保留模态的数量为40，前四阶主模态频率和振型见图 2。

转子系统运动方程可写为

$ \mathit{\boldsymbol{M\ddot u}} + \mathit{\boldsymbol{G\dot u}} + \mathit{\boldsymbol{Ku}} = {\mathit{\boldsymbol{F}}^n} + {\mathit{\boldsymbol{F}}^u} $

(1)

式中：F^u为不平衡力向量；Fⁿ为作用在转子系统上的非线性力向量，包括滚动轴承的非线性作用力和挤压油膜阻尼器的非线性作用力；M，K，G分别为惯性矩阵、刚度矩阵以及陀螺矩阵；u为节点位移列向量。

对于双转子系统，陀螺矩阵可写为

$ \mathit{\boldsymbol{G}} = {\rm{diag}}\left( {{\mathit{\boldsymbol{g}}_i}} \right) $

(2)

其中

$ {\mathit{\boldsymbol{g}}_i} = {\omega _l}\left[ {\begin{array}{*{20}{c}} 0&0&0&0\\ 0&0&0&0\\ 0&0&0&{ - {I_{pi}}}\\ 0&0&{{I_{pi}}}&0 \end{array}} \right]\;\;\;\;\;1 \le i \le 4 $

式中: I_pi为第i个圆盘的极转动惯量；ω_l表示圆盘所在转子的转速，i=1, 2时，ω_l=ω₁，i=3, 4时，ω_l=ω₂。

方程(1)为非线性微分方程组，其求解效率取决于所采用的求解方法以及模型的自由度数。为提高计算求解效率，本文采用有限元法和模态综合法建立如式(1)所示的转子系统模型。

根据固定界面模态综合法将系统自由度划分为内部自由度和界面自由度。转子支承处的自由度为界面自由度，其他为内部自由度。由此，式(1)可改写为

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{M}}_{{\rm{II}}}}}&{{\mathit{\boldsymbol{M}}_{{\rm{IJ}}}}}\\ {{\mathit{\boldsymbol{M}}_{{\rm{JI}}}}}&{{\mathit{\boldsymbol{M}}_{{\rm{JJ}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\ddot u}}}_{\rm{I}}}}\\ {{{\mathit{\boldsymbol{\ddot u}}}_{\rm{J}}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{G}}_{{\rm{II}}}}}&{{\mathit{\boldsymbol{G}}_{{\rm{IJ}}}}}\\ {{\mathit{\boldsymbol{G}}_{{\rm{JI}}}}}&{{\mathit{\boldsymbol{G}}_{{\rm{JJ}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\dot u}}}_{\rm{I}}}}\\ {{{\mathit{\boldsymbol{\dot u}}}_{\rm{J}}}} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{K}}_{{\rm{II}}}}}&{{\mathit{\boldsymbol{K}}_{{\rm{IJ}}}}}\\ {{\mathit{\boldsymbol{K}}_{{\rm{JI}}}}}&{{\mathit{\boldsymbol{K}}_{{\rm{JJ}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{u}}_{\rm{I}}}}\\ {{\mathit{\boldsymbol{u}}_{\rm{J}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}}_{\rm{I}}^u}\\ \mathit{\boldsymbol{0}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{0}}\\ {\mathit{\boldsymbol{F}}_{\rm{J}}^n} \end{array}} \right]} \end{array} $

(3)

式中：u_I为内部自由度合集，仅有不平衡力作用；u_J为界面自由度合集，轴承Ⅰ，Ⅱ和Ⅲ处的界面自由度存在弹性支承线性作用力和挤压油膜阻尼器的非线性作用力；轴承Ⅳ处的界面自由度仅受到中介轴承的非线性作用力。

根据固定界面模态综合法，物理坐标与模态坐标的变换关系为

$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{u}}_{\rm{I}}}}\\ {{\mathit{\boldsymbol{u}}_{\rm{J}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k}}&{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{I}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{q}}_k}}\\ {{\mathit{\boldsymbol{u}}_{\rm{J}}}} \end{array}} \right] = \mathit{\boldsymbol{T}}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{q}}_k}}\\ {{\mathit{\boldsymbol{u}}_{\rm{J}}}} \end{array}} \right] $

(4)

式中：Φ_k为质量归一化的主模态矩阵；Φ_c为质量归一化的约束模态矩阵；I为单位矩阵；q_k为主模态坐标；T为变换矩阵。

Φ_c可通过式(5)计算得到

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c} = - {\left( {{\mathit{\boldsymbol{K}}_{{\rm{II}}}}} \right)^{ - 1}}{\mathit{\boldsymbol{K}}_{{\rm{IJ}}}} $

(5)

为求解主模态，忽略轴段转动惯量以及轮盘陀螺力矩并约束所有界面自由度。然后对此系统进行模态分析，即求解对应式(6)的特征值问题，可获得质量归一化的主模态矩阵Φ_k和主模态频率向量Ω_k，则

$ {\mathit{\boldsymbol{M}}_{{\rm{II}}}}{{\mathit{\boldsymbol{\ddot u}}}_{\rm{I}}} + {\mathit{\boldsymbol{K}}_{{\rm{II}}}}{\mathit{\boldsymbol{u}}_{\rm{I}}} = 0 $

(6)

将式(4)代入式(3)并左乘T^T得到缩减后的系统运动方程为

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{I}}&{{{\mathit{\boldsymbol{\bar M}}}_{{\rm{IJ}}}}}\\ {{{\mathit{\boldsymbol{\bar M}}}_{{\rm{JI}}}}}&{{{\mathit{\boldsymbol{\bar M}}}_{{\rm{JJ}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\ddot q}}}_{\rm{k}}}}\\ {{{\mathit{\boldsymbol{\ddot u}}}_{\rm{J}}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\bar G}}}_{{\rm{II}}}}}&{{{\mathit{\boldsymbol{\bar G}}}_{{\rm{IJ}}}}}\\ {{{\mathit{\boldsymbol{\bar G}}}_{{\rm{JI}}}}}&{{{\mathit{\boldsymbol{\bar G}}}_{{\rm{JJ}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\dot q}}}_{\rm{k}}}}\\ {{{\mathit{\boldsymbol{\dot u}}}_{\rm{J}}}} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\bar K}}}_{{\rm{II}}}}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&{{{\mathit{\boldsymbol{\bar K}}}_{{\rm{JJ}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{q}}_{\rm{k}}}}\\ {{\mathit{\boldsymbol{u}}_{\rm{J}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k^{\rm{T}}\mathit{\boldsymbol{F}}_{\rm{I}}^u}\\ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c^{\rm{T}}\mathit{\boldsymbol{F}}_{\rm{I}}^u} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{0}}\\ {\mathit{\boldsymbol{F}}_{\rm{J}}^n} \end{array}} \right]} \end{array} $

(7)

其中

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\bar M}}}_{{\rm{JJ}}}} = {\mathit{\boldsymbol{M}}_{{\rm{JJ}}}} + {\mathit{\boldsymbol{M}}_{{\rm{JI}}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c} + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_{{\rm{II}}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c} + {\mathit{\boldsymbol{M}}_{{\rm{IJ}}}}} \right)\\ {{\mathit{\boldsymbol{\bar M}}}_{{\rm{IJ}}}} = {{\mathit{\boldsymbol{\bar M}}}_{{\rm{JI}}}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_{{\rm{II}}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c} + {\mathit{\boldsymbol{M}}_{{\rm{IJ}}}}} \right)\\ {{\mathit{\boldsymbol{\bar G}}}_{{\rm{II}}}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k^{\rm{T}}{\mathit{\boldsymbol{G}}_{{\rm{II}}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k}\\ {{\mathit{\boldsymbol{\bar G}}}_{{\rm{JJ}}}} = {\mathit{\boldsymbol{G}}_{{\rm{JJ}}}} + {\mathit{\boldsymbol{G}}_{{\rm{JI}}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c} + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c^{\rm{T}}\left( {{\mathit{\boldsymbol{G}}_{{\rm{II}}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c} + {\mathit{\boldsymbol{G}}_{{\rm{IJ}}}}} \right)\\ {{\mathit{\boldsymbol{\bar G}}}_{{\rm{IJ}}}} = {{\mathit{\boldsymbol{\bar G}}}_{{\rm{JI}}}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k^{\rm{T}}{\mathit{\boldsymbol{G}}_{{\rm{II}}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c} + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k^{\rm{T}}{\mathit{\boldsymbol{G}}_{{\rm{IJ}}}}\\ {{\mathit{\boldsymbol{\bar K}}}_{II}} = {\rm{diag}}\left( {\mathit{\Omega }_r^2} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;1 \le r \le n\\ {{\mathit{\boldsymbol{\bar K}}}_{{\rm{JJ}}}} = {\mathit{\boldsymbol{K}}_{{\rm{IJ}}}} + {\mathit{\boldsymbol{K}}_{{\rm{JI}}}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c} + {\mathit{\boldsymbol{k}}_{{\rm{JJ}}}} \end{array} \right. $

(8)

式中：Ω_r为r阶主模态的频率；n为保留主模态的阶数；F_I^u为作用在内部自由度上的不平衡力向量；F_Jⁿ为作用在界面自由度上的非线性作用力向量，且有

$ \mathit{\boldsymbol{F}}_{\rm{J}}^n = - {\mathit{\boldsymbol{k}}_{{\rm{JJ}}}}{\mathit{\boldsymbol{u}}_{\rm{J}}} + \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}}_{\rm{J}}^{{\rm{SFD}}}}\\ {\mathit{\boldsymbol{F}}_{\rm{J}}^{\rm{B}}} \end{array}} \right] $

(9)

式中，F_J^SFD=[f_1x^sfd, f_1y^sfd, …, f_Nx^sfd，f_Ny^sfd]^T为挤压油膜阻尼器非线性力向量。F_J^B=[f_x^b f_y^b -f_x^b -f_y^b]^T为中介轴承的非线性作用力向量。k_JJ为鼠笼式弹支的刚度矩阵，其表达式为

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\; {\mathit{\boldsymbol{k}}_{{\rm{JJ}}}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{k}}_{\rm{J}}^{\rm{B}}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}} \end{array}} \right]\\ \mathit{\boldsymbol{k}}_{\rm{J}}^{\rm{B}} = {\rm{diag}}\left( {\mathit{\boldsymbol{k}}_{{\rm{J}}i}^{\rm{B}}} \right)\;\;\;\;\;1 \le i \le N\\ \mathit{\boldsymbol{k}}_{{\rm{J}}i}^{\rm{B}} = \left[ {\begin{array}{*{20}{c}} {{k_{{\rm{J}}i}}}&0\\ 0&{{k_{{\rm{J}}i}}} \end{array}} \right]\;\;\;\;1 \le i \le N \end{array} $

(10)

式中：k_Ji为鼠笼式弹支的刚度系数，N为鼠笼式弹支的数量。

将式(7)展开得到

$ {{\mathit{\boldsymbol{\ddot q}}}_k} + {{\mathit{\boldsymbol{\bar G}}}_{{\rm{II}}}}{{\mathit{\boldsymbol{\dot q}}}_k} + {{\mathit{\boldsymbol{\bar K}}}_{{\rm{II}}}}{\mathit{\boldsymbol{q}}_k} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k^{\rm{T}}\mathit{\boldsymbol{F}}_{\rm{I}}^u - {{\mathit{\boldsymbol{\bar M}}}_{{\rm{IJ}}}}{{\mathit{\boldsymbol{\ddot u}}}_{\rm{J}}} - {{\mathit{\boldsymbol{\bar G}}}_{{\rm{IJ}}}}{{\mathit{\boldsymbol{\dot u}}}_{\rm{J}}} $

(11)

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\bar M}}}_{{\rm{JJ}}}}{{\mathit{\boldsymbol{\ddot u}}}_{\rm{J}}} + {{\mathit{\boldsymbol{\bar G}}}_{{\rm{JJ}}}}{{\mathit{\boldsymbol{\dot u}}}_{\rm{J}}} + {{\mathit{\boldsymbol{\bar K}}}_{{\rm{JJ}}}}{\mathit{\boldsymbol{u}}_{\rm{J}}} = }\\ {\mathit{\boldsymbol{F}}_{\rm{J}}^n + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c^{\rm{T}}\mathit{\boldsymbol{F}}_{\rm{I}}^u - {{\mathit{\boldsymbol{\bar M}}}_{{\rm{JI}}}}{{\mathit{\boldsymbol{\ddot q}}}_k} - {{\mathit{\boldsymbol{\bar G}}}_{{\rm{IJ}}}}{{\mathit{\boldsymbol{\dot q}}}_k}} \end{array} $

(12)

观察式(11, 12)可发现，式(11)中无非线性作用力，采用显式Newmark法求解；式(12)中存在挤压油膜阻尼器和中介轴承的非线性作用力，采用隐式Newmark法求解。

根据Newmark-β法的假设，在区间[t_n t_n+1]内有

$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot q}}_k^{n + 1}}\\ {\mathit{\boldsymbol{\dot u}}_{\rm{J}}^{n + 1}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot q}}_k^n}\\ {\mathit{\boldsymbol{\dot u}}_{\rm{J}}^n} \end{array}} \right] + \left\{ {\left( {1 - \beta } \right)\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot q}}_k^n}\\ {\mathit{\boldsymbol{\ddot u}}_{\rm{J}}^n} \end{array}} \right] + \beta \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot q}}_k^{n + 1}}\\ {\mathit{\boldsymbol{\ddot u}}_{\rm{J}}^{n + 1}} \end{array}} \right]} \right\}\Delta t $

(13)

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{q}}_k^{n + 1}}\\ {\mathit{\boldsymbol{u}}_{\rm{J}}^{n + 1}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{q}}_k^n}\\ {\mathit{\boldsymbol{u}}_{\rm{J}}^n} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot q}}_k^n}\\ {\mathit{\boldsymbol{\dot u}}_{\rm{J}}^n} \end{array}} \right]\Delta t + }\\ {\left\{ {\left( {0.5 - \alpha } \right)\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot q}}_k^n}\\ {\mathit{\boldsymbol{\ddot u}}_{\rm{J}}^n} \end{array}} \right] + \alpha \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot q}}_k^{n + 1}}\\ {\mathit{\boldsymbol{\ddot u}}_{\rm{J}}^{n + 1}} \end{array}} \right]} \right\}\Delta {t^2}} \end{array} $

(14)

式中t_n+1=t_n+Δt为时间步长，上标n和n+1分别表示t_n和t_n+1时刻。

根据式(13, 14)可得

$ \begin{array}{l} \mathit{\boldsymbol{\ddot q}}_k^{n + 1} = a\mathit{\boldsymbol{q}}_k^{n + 1} - \mathit{\boldsymbol{A}}_q^n\\ \mathit{\boldsymbol{\dot q}}_k^{n + 1} = b\mathit{\boldsymbol{q}}_k^{n + 1} - \mathit{\boldsymbol{B}}_q^n \end{array} $

(15)

$ \begin{array}{l} \mathit{\boldsymbol{\ddot u}}_{\rm{J}}^{n + 1} = a\mathit{\boldsymbol{u}}_{\rm{J}}^{n + 1} - \mathit{\boldsymbol{A}}_{\rm{J}}^n\\ \mathit{\boldsymbol{\dot u}}_{\rm{J}}^{n + 1} = b\mathit{\boldsymbol{u}}_{\rm{J}}^{n + 1} - \mathit{\boldsymbol{B}}_{\rm{J}}^n \end{array} $

(16)

其中

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}}_q^n = \frac{1}{{\alpha {{\left( {\Delta t} \right)}^2}}}\mathit{\boldsymbol{q}}_k^n + \frac{1}{{\alpha \Delta t}}\mathit{\boldsymbol{\dot q}}_k^n + \left( {\frac{1}{{2\alpha }} - 1} \right)\mathit{\boldsymbol{\ddot q}}_k^n}\\ {\mathit{\boldsymbol{A}}_{\rm{J}}^n = \frac{1}{{\alpha {{\left( {\Delta t} \right)}^2}}}\mathit{\boldsymbol{u}}_{\rm{J}}^n + \frac{1}{{\alpha \Delta t}}\mathit{\boldsymbol{\dot u}}_{\rm{J}}^n + \left( {\frac{1}{{2\alpha }} - 1} \right)\mathit{\boldsymbol{\ddot u}}_{\rm{J}}^n}\\ {\mathit{\boldsymbol{B}}_q^n = \frac{\beta }{{\alpha \Delta t}}\mathit{\boldsymbol{q}}_k^n + \left( {\frac{\beta }{\alpha } - 1} \right)\mathit{\boldsymbol{\dot q}}_k^n + \left( {\frac{\beta }{\alpha } - 2} \right)\mathit{\boldsymbol{\ddot q}}_k^n}\\ {\mathit{\boldsymbol{B}}_{\rm{J}}^n = \frac{\beta }{{\alpha \Delta t}}\mathit{\boldsymbol{u}}_{\rm{J}}^n + \left( {\frac{\beta }{\alpha } - 1} \right)\mathit{\boldsymbol{\dot u}}_{\rm{J}}^n + \left( {\frac{\beta }{\alpha } - 2} \right)\mathit{\boldsymbol{\ddot u}}_{\rm{J}}^n}\\ {a = \frac{1}{{\alpha {{\left( {\Delta t} \right)}^2}}}\;\;\;\;b = \frac{\beta }{{\alpha \Delta t}}} \end{array} $

(17)

将式(13—17)代入式(11, 12)得

$ \mathit{\boldsymbol{q}}_k^{n + 1} = \mathit{\boldsymbol{S}}_q^{ - 1}\left( {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k^{\rm{T}}\mathit{\boldsymbol{F}}_{\rm{I}}^{u,n + 1} - {\mathit{\boldsymbol{V}}_q}\mathit{\boldsymbol{u}}_{\rm{J}}^{n + 1} + {\mathit{\boldsymbol{W}}_q}} \right) $

(18)

$ \begin{array}{*{20}{c}} {\left( {{\mathit{\boldsymbol{S}}_{\rm{J}}} - {\mathit{\boldsymbol{V}}_{\rm{J}}}\mathit{\boldsymbol{S}}_q^{ - 1}{\mathit{\boldsymbol{V}}_q}} \right)\mathit{\boldsymbol{u}}_{\rm{J}}^{n + 1} = \mathit{\boldsymbol{F}}_{\rm{I}}^{n,n + 1} + }\\ {\left( {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_c^{\rm{T}} - {\mathit{\boldsymbol{V}}_{\rm{J}}}\mathit{\boldsymbol{S}}_q^{ - 1}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k^{\rm{T}}} \right)\mathit{\boldsymbol{F}}_I^{u,n + 1} - {\mathit{\boldsymbol{V}}_{\rm{J}}}\mathit{\boldsymbol{S}}_q^{ - 1}{\mathit{\boldsymbol{W}}_q} + {\mathit{\boldsymbol{W}}_{\rm{J}}}} \end{array} $

(19)

其中

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{W}}_q} = {\mathit{\boldsymbol{M}}_{{\rm{II}}}}\mathit{\boldsymbol{A}}_q^n + {\mathit{\boldsymbol{G}}_{{\rm{II}}}}\mathit{\boldsymbol{B}}_q^n + {\mathit{\boldsymbol{M}}_{{\rm{IJ}}}}\mathit{\boldsymbol{A}}_{\rm{J}}^n + {\mathit{\boldsymbol{G}}_{{\rm{IJ}}}}\mathit{\boldsymbol{B}}_{\rm{J}}^n\\ {\mathit{\boldsymbol{W}}_{\rm{J}}} = {\mathit{\boldsymbol{M}}_{{\rm{JJ}}}}\mathit{\boldsymbol{A}}_{\rm{J}}^n + {\mathit{\boldsymbol{G}}_{{\rm{JJ}}}}\mathit{\boldsymbol{B}}_{\rm{J}}^n + {\mathit{\boldsymbol{M}}_{{\rm{JI}}}}\mathit{\boldsymbol{A}}_q^n + {\mathit{\boldsymbol{G}}_{{\rm{JI}}}}\mathit{\boldsymbol{B}}_q^n\\ {\mathit{\boldsymbol{S}}_q} = a{\mathit{\boldsymbol{M}}_{{\rm{II}}}} + b{\mathit{\boldsymbol{G}}_{{\rm{II}}}} + {\mathit{\boldsymbol{K}}_{{\rm{II}}}}\\ {\mathit{\boldsymbol{S}}_J} = a{\mathit{\boldsymbol{M}}_{{\rm{JJ}}}} + b{\mathit{\boldsymbol{G}}_{{\rm{JJ}}}} + {\mathit{\boldsymbol{K}}_{{\rm{JJ}}}}\\ {\mathit{\boldsymbol{V}}_q} = a{\mathit{\boldsymbol{M}}_{{\rm{IJ}}}} + b{\mathit{\boldsymbol{G}}_{{\rm{IJ}}}}\\ {\mathit{\boldsymbol{V}}_J} = a{\mathit{\boldsymbol{M}}_{{\rm{JI}}}} + b{\mathit{\boldsymbol{G}}_{{\rm{JI}}}}\\ \mathit{\boldsymbol{F}}_{\rm{J}}^{n,n + 1} = \mathit{\boldsymbol{F}}_{\rm{J}}^{n,n}\left( {\mathit{\boldsymbol{u}}_{\rm{J}}^{n + 1},\mathit{\boldsymbol{\dot u}}_{\rm{J}}^{n + 1}} \right) \end{array} \right. $

(20)

将式(16)代入式(20)中的F_J^{n, n+1}=F_J^{n, n}(u_Jⁿ⁺¹, $ \mathit{\boldsymbol{\dot{u}}}_{\text{J}}^{\mathit{n}+1}$)，随后代入式(19)可以得到关于u_Jⁿ⁺¹的非线性代数方程组。求解(19)得到u_Jⁿ⁺¹后，将其代入式(16)得到q_kⁿ⁺¹。随后通过式(15, 16)得到$ \mathit{\boldsymbol{\dot{q}}}_{\mathit{k}}^{\mathit{n}+1}, \mathit{\boldsymbol{\ddot{q}}}_{\mathit{k}}^{\mathit{n}+1}, \mathit{\boldsymbol{\dot{u}}}_{\text{J}}^{\mathit{n}+1}, \mathit{\boldsymbol{\ddot{u}}}_{\text{J}}^{\mathit{n}+1}$。然后将该时刻的响应值作为下一时间区间的初始条件，开始新的计算。当所有计算完成之后，物理坐标下的响应可通过式(3)变换得到。

非线性方程组(14)的维数等于非线性力处的自由度个数，而数值计算的效率主要取决于数值算法求解非线性方程组的效率。所以，在每个时刻非线性方程组的计算规模仅与非线性自由度个数有关，进而提高计算效率。

1.2 非线性力模型

根据短轴承假设和Reynolds边界条件^[12]，式(9)中挤压油膜阻尼器非线性力为

$ \left\{ \begin{array}{l} f_x^{{\rm{std}}} = - \frac{{\mu R{L^3}}}{{{c^2}{{\left( {{x^2} + {y^2}} \right)}^{1/2}}}}\left[ {x\left( {\dot \varepsilon {I_2} + \varepsilon \dot \psi {I_1}} \right) - } \right.\\ \;\;\;\;\;\;\;\left. {y\left( {\dot \varepsilon {I_1} + \varepsilon \dot \psi {I_3}} \right)} \right]\\ f_y^{{\rm{sfd}}} = - \frac{{\mu R{L^3}}}{{{c^2}{{\left( {{x^2} + {y^2}} \right)}^{1/2}}}}\left[ {y\left( {\dot \varepsilon {I_2} + \varepsilon \dot \psi {I_1}} \right) + } \right.\\ \;\;\;\;\;\;\;\left. {x\left( {\dot \varepsilon {I_1} + \varepsilon \dot \psi {I_3}} \right)} \right] \end{array} \right. $

(21)

其中

$ \left\{ \begin{array}{l} \varepsilon = \sqrt {{x^2} + {y^2}} /c\\ \dot \varepsilon = \left( {x\dot x + y\dot y} \right)/\left( {c\sqrt {{x^2} + {y^2}} } \right)\\ \dot \psi = \left( {y\dot x - x\dot y} \right)/\left( {{x^2} + {y^2}} \right)\\ \tan \psi = y/x \end{array} \right. $

(22)

式中：x和y分别为挤压油膜阻尼器所在轴径处的水平和垂直方向位移；I_j(j=1, 2, 3)为Sommerfeld积分；R和L分别为挤压油膜阻尼器的半径和轴向承载宽度；μ和c分别为滑油动力黏度系数和阻尼器的径向间隙。

基于Hertz接触理论和纯滚动假设，式(9)中中介轴承的非线性力可表述为^[7]

$ \left\{ \begin{array}{l} f_x^{\rm{b}} = {k_{\rm{n}}}\sum\limits_{j = 1}^{{N_{\rm{b}}}} {u_{{\theta _j}}^\xi H\left( {{u_{{\theta _j}}}} \right)\sin {\theta _j}} \\ f_y^{\rm{b}} = {k_{\rm{n}}}\sum\limits_{j = 1}^{{N_{\rm{b}}}} {u_{{\theta _j}}^\xi H\left( {{u_{{\theta _j}}}} \right)\cos {\theta _j}} \end{array} \right. $

(23)

其中

$ \left\{ \begin{array}{l} {\theta _j} = \frac{{2{\rm{ \mathsf{ π} }}\left( {j - 1} \right)}}{{{N_{\rm{b}}}}} + {\omega _c}t\\ {u_{{\theta _j}}} = \left( {{x^{{\rm{ir}}}} - {x^{{\rm{or}}}}} \right)\cos {\theta _j} + \left( {{y^{{\rm{ir}}}} - {y^{{\rm{or}}}}} \right)\sin {\theta _j} - \frac{\gamma }{2}\\ H\left( {{\theta _j}} \right) = \left\{ \begin{array}{l} 0\;\;\;\;\;\;{u_{{\theta _j}}} \le 0\\ {u_{{\theta _j}}}\;\;\;\;{u_{{\theta _j}}} > 0 \end{array} \right.\\ {\omega _c} = \left( {{\omega _{{\rm{in}}}}r + {\omega _{{\rm{out}}}}R} \right)/\left( {R + r} \right) \end{array} \right. $

(24)

式中：上标ir和or分别表示轴承的内环和外环；k_n表示赫兹接触刚度，与互相接触的材料和形状有关；N_b表示滚珠的个数；对于本文中的中介轴承而言$ \xi =\frac{3}{2}$；θ_j表示t时刻第j个滚珠的角度；γ和ω_c分别表示轴承的径向间隙和保持架的旋转速度；u_θj表示第j个滚珠的弹性变形；r和R分别表示轴承内圈和外圈的半径；ω_in和ω_out分别表示轴承内圈和外圈的转速。

2 转子系统的非线性响应特性分析

利用上述转子数据和式(1—24)，对不同转速比情况下转子系统的非线性响应特性进行了计算和分析。计算转速范围为4~400 rad/s，计算步长为2 rad/s。计算过程中，利用前一渐进稳态响应的结果作为下一个渐进稳态计算的初始值。不失一般性，本文以盘2和盘4的响应为例分析内、外转子的不平衡响应。

图 3为λ=±1.65和λ=±2情况下盘2的三维频谱图。图 4为λ=±1.65和λ=±2情况下盘4的三维频谱图。图 5为不同转速比情况下盘2响应幅值随转速变化的曲线。图 6为不同转速比情况下盘2响应的分岔图。由图 3—6可以看出：(1) 由于中介轴承的存在，内、外转子的响应互相耦合，不同转速比情况下有着相似的三维频谱图：盘2以及盘4处的响应在整个转速范围内均存在内、外转子的自转频率ω₁和ω₂，分别对应内、外转子的不平衡激励频率，这两者的贡献在系统响应中占绝对优势，此现象称之为交叉激振现象，见图 3，4。

(2) 对比图 3和图 4可以很明显看出盘2和盘4的响应规律类似：交叉激振频率为主，组合频率成分几乎完全相同，因此后续仅对盘2的响应进行分析。

(3) 由于SFD和中介轴承的非线性作用力，盘2和盘4的响应中除了ω₁和ω₂的频率成分之外还出现了两者的组合频率。2ω₁+ω₂, ω₁+2ω₂两者最为明显，在所有转速比下均能较明显地观察到。较低的组合频率ω₁+ω₂仅在λ=±2情况下出现；较高的组合频率3ω₁+2ω₂, 2ω₁+3ω₂则仅在λ=±1.65情况下能观察到，见图 3, 4。

(4) 转速比绝对值相同的情况下，同向旋转的各阶临界转速均不小于反向旋转；增大转速比绝对值会使得内转子为主激励的临界转速上升，外转子为主激励的临界转速下降，见图 3和表 5。从表 5中可以看到当λ=±2时在0~400 rad/s范围内存在两阶以外转子为主激励的临界转速，而λ=±1.65时仅存在一阶。这说明增大转速比可能导致工作转速范围内临界转速阶数增加，不利于转子系统的启动和停止。此外，从图 5可以看出，总的来说，反向旋转情况下振幅略小于同向旋转。

(5) 从图 5可以看到，同转情况下临界转速附近盘2响应中ω₁和ω₂频率成分的幅值均小于反转情况。

(6) 转速比对于转子系统轴心轨迹和运动的周期性存在较大影响：λ=±1.65情况下转子系统处于4周期的次谐响应状态，轴心轨迹呈对称的花瓣状；λ=±2情况下转子系统则处于单周期运动状态，轴心轨迹呈一定宽度的不规则闭合曲线带，见图 6—8。此外，从图 7可以看到相同转速下λ=-1.65时庞加莱图上所形成的4个点放大后均为扭曲、封闭的图形，而λ=1.65时则形成了较规则的圆形，说明反向旋转导致转子系统响应更加复杂，进而导致系统在相空间中的运动轨迹具有更复杂的几何结构。

3 试验验证

为了验证本文建模及计算结果的正确性和准确性，通过试验验证了λ=±1.65情况下的数值计算结果。双转子试验器几何尺寸、轴段材料参数、支承刚度、阻尼器参数、轮盘参数及不平衡量配置见图 1及表 4。试验过程中，待转子加速至目标转速并稳定后通过东方振动所的数据采集卡采集传感器信号。4个电涡流位移传感器分别测量盘2和盘4互相垂直的两方向振动响应。实验结束后通过配套的DASP软件分析振动信号并导出分析结果。由于盘2和盘4响应规律类似，故仅列出盘2的试验结果。

图 9, 10为λ=-1.65情况下盘2数值计算结果与试验结果的对比。图 11, 12为λ=1.65情况下盘2数值结果与试验结果的对比。通过图 9—12可以看出本文的计算结果与试验结果吻合较好，说明了本文建模和计算的准确性。

4 结论

以航空发动机双转子系统为研究对象，建立了含SFD和中介轴承的双转子系统耦合动力学模型，模型考虑了SFD和中介轴承非线性力的影响。通过数值模拟分析了不同转速比下双转子系统的非线性响应特性，并通过试验进行验证。本文主要结论如下：

(1) 本文方法可以方便快捷地建立双转子系统的非线性动力学模型，且改进的Newmark算法求解效率主要取决于转子系统中非线性力处的自由度数，因而具有较高的计算效率。

(2) 受陀螺力矩的影响，转速比绝对值相同的情况下，同向旋转双转子系统的临界转速不小于反向旋转，且反转情况下盘2的响应幅值略小于同转情况。增大转速比绝对值会使得内转子为主激励的临界转速上升，外转子为主激励的临界转速下降。

(3) 中介轴承使得内、外转子的振动互相耦合：转子系统的响应中有较明显的交叉激振现象，且响应中除了内外转子的不平衡激励频率之外还出现了两者的组合频率，但是不同转速比情况下的频率组合情况并不一样。

(4) λ=±1.65情况下，转子系统处于4周期的次谐响应状态，轴心轨迹呈对称的花瓣状；λ=±2情况下，转子系统处于单周期运动状态，轴心轨迹呈一定宽度的不规则闭合曲线带；反向旋转导致转子系统响应更加复杂，进而导致相空间中的运动轨迹更加复杂。