CN103091107A - 汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法 - Google Patents

Abstract

本发明公开了旋转机械转子动力学技术领域中的一种汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法。该方法将逐步积分法的思想与传统传递矩阵法相结合，对非线性方程组中各量进行增量改造，将非线性微分方程转化为关于振动状态增量的线性微分方程；在对表达式进行增量改造基础上，建立起机组基于Riccati变换的弯扭耦合振动增量传递矩阵表达式，通过代入多次迭代计算的外力和外力矩，即可进行机组碰摩时弯扭耦合振动瞬态动态响应分析。本发明可适用于类似汽轮发电机组这样的复杂转子系统，分析结果更准确；也可用于强非线性系统的求解，拓展了传递矩阵法在非线性领域的应用。

Description

汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法

技术领域

本发明属于旋转机械转子动力学技术领域，尤其涉及一种汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法。

背景技术

随着汽轮发电机组向大容量高效率方向的发展，转子与静子的间隙越来越小，导致转静子的碰摩故障不断发生，碰摩故障已成为机组严重的主要故障模式之一。碰摩一旦发生，往往导致机组结构上的损坏，甚至导致灾难性的后果，因此研究大型旋转机械的动静摩擦的振动特性，对于充分利用振动信号诊断摩擦故障，防止转子摩擦向中、晚期过渡，确保设备的安全稳定运行，具有非常重要的意义。

关于转子系统碰摩问题的研究一直是转子动力学研究的重要课题，其振动动态特性响应分析的研究是转子动力学的主要内容之一。目前已有众多学者进行了该方面的研究，但这些研究大多针对形如Jeffcott转子的简单转子系统，而对复杂转子系统（如汽轮发电机组转子系统）的碰摩故障的研究则比较少见，复杂转子系统具有比简单转子系统更加复杂的动力学特性，所以有必要研究复杂转子系统的碰摩故障以揭示其动力学特性。分析复杂转子系统的方法主要有模态叠加法、直接时域积分法和传递矩阵法三类，对于汽轮发电机组这样链式结构的转子系统，采用传递矩阵法是一种简单有效的方法；由于机组碰摩故障发生时，轴系弯振和扭振参量相互耦合，机组轴系是一个高度复杂的非线性振动系统，对于这样的系统，采用线性振动理论中的一些求解方法，已不再适用。

发明内容

本发明针对现有技术的不足，提出一种汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法，将逐步积分法的思想与传统传递矩阵法相结合，对非线性方程组中各量进行增量改造，将非线性微分方程转化为关于振动状态增量的线性微分方程，建立起机组基于Riccati变换的弯扭耦合振动增量传递矩阵表达式，可用于汽轮发电机组弯扭耦合振动瞬态动态响应分析。

为了实现上述目的，本发明提出的技术方案是，一种汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法，其特征是所述方法包括：

步骤1：设定初始时刻t₀和步进长度Δt；

步骤2：采集时刻t₀机组各节点的位移和速度，计算时刻t₀机组各节点的加速度;

步骤3：令θ＝0；

步骤4：根据时刻t₀+θΔt机组各节点的位移、速度和加速度，计算时刻t₀+θΔt机组各轴段的外力/外力矩；

步骤5：将时刻t₀+θΔt机组各轴段的外力/外力矩分别作为时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩；

步骤6：根据时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩，利用基于Riccati变换的轴系弯扭耦合振动增量传递矩阵计算时刻t₀+(θ+1)Δt机组各节点的位移增量；

步骤7：根据时刻t₀+(θ+1)Δt机组各节点的位移增量，利用改进的wilson-θ增量表达式计算时刻t₀+(θ+1)Δt的位移、速度和加速度；

步骤8：根据时刻t₀+(θ+1)Δt的位移、速度和加速度，计算时刻t₀+(θ+1)Δt机组各轴段的外力/外力矩；

步骤9：判断是否满足设定条件，如果满足设定条件，则执行步骤11；否则，执行步骤10；

所述设定条件为：时刻t₀+(θ+1)Δt机组各轴段的外力/外力矩与时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩的差值绝对值小于设定值；

步骤10：令时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩等于时刻t₀+(θ+1)Δt机组各轴段的外力/外力矩，返回步骤6；

步骤11：令θ=θ+1，返回步骤5。

所述计算时刻t₀+θΔt机组各轴段的外力/外力矩具体为，当轴段为碰摩轴段时，计算时刻t₀+θΔt机组该轴段的摩擦力矩；当轴段为轴颈轴段时，计算时刻t₀+θΔt机组该轴段的油膜力；当轴段为动叶栅所在轴段时，计算时刻t₀+θΔt机组动叶栅受到的蒸汽弯曲力矩；当轴段为发电机绕组轴段时；计算时刻t₀+θΔt机组的发电机绕组轴段受到的电磁力矩。

所述利用基于Riccati变换的轴系弯扭耦合振动增量传递矩阵计算时刻t₀+(θ+1)Δt机组各节点的位移增量的计算公式为：

$\left\{\begin{array}{c}{f}_{i+1}^L=U_{11}{f_i}^L+U_{12,i}\·e_i^L+F_{f,i}\\ e_{i+1}^L=U_{21}{f_i}^L+U_{22,i}\·e_i^L+F_{e,i}\end{array}\right.,$

其中，f_i ^L为第i个轴段左侧力向量增量，

为第i个轴段左侧位移向量增量，U₁₁和U₂₁分别为基于Riccati变换的轴系弯扭耦合振动增量传递矩阵参数，U_12,i、U_22,i、F_f,i和F_e,i分别为第i个轴段的基于Riccati变换的轴系弯扭耦合振动增量传递矩阵参数；

$U_{11}=\left|\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ L& 0& 1& 0& 0\\ 0& L& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right|,$ L轴段的长度；

$U_{21}=\left|\begin{array}{ccccc}{\mathrm{\γ}}_3& 0& {\mathrm{\γ}}_2& 0& 0\\ 0& {\mathrm{\γ}}_3& 0& {\mathrm{\γ}}_2& 0\\ {\mathrm{\γ}}_2& 0& {\mathrm{\γ}}_1& 0& 0\\ 0& {\mathrm{\γ}}_2& 0& {\mathrm{\γ}}_1& 0\\ 0& 0& 0& 0& \frac{1}{K_t}\end{array}\right|,$ K_t为弹性轴段扭转刚度， ${\mathrm{\γ}}_1=\frac{L}{\mathrm{EJ}},$ ${\mathrm{\γ}}_2=\frac{L^2}{2\mathrm{EJ}},$

E为轴段的弹性模量，J为轴段等效抗弯刚度；

$U_{12,i}=\left|\begin{array}{ccccc}{A}_{11}& A_{12}& 0& 0& A_{13}\\ A_{22}& A_{21}& 0& 0& A_{23}\\ {\mathrm{LA}}_{11}& L{A}_{12}& A_{31}& A_{32}& L{A}_{13}\\ L{A}_{22}& L{A}_{21}& -A_{32}& A_{31}& L{A}_{23}\\ A_{51}& A_{52}& 0& 0& C_1\end{array}\right|,$

$U_{22,i}=\left|\begin{array}{ccccc}{\mathrm{\γ}}_3{A}_{11}+1& {\mathrm{\γ}}_3{A}_{12}& {\mathrm{\γ}}_2{A}_{31}+L& {\mathrm{\γ}}_2{A}_{32}& {\mathrm{\γ}}_3{A}_{13}\\ {\mathrm{\γ}}_3{A}_{22}& {\mathrm{\γ}}_3{A}_{21}+1& -{\mathrm{\γ}}_2{A}_{32}& {\mathrm{\γ}}_2{A}_{31}+L& {\mathrm{\γ}}_3{A}_{23}\\ {\mathrm{\γ}}_2{A}_{11}& {\mathrm{\γ}}_2{A}_{12}& {\mathrm{\γ}}_1{A}_{31}+1& {\mathrm{\γ}}_1{A}_{32}& {\mathrm{\γ}}_2{A}_{13}\\ {\mathrm{\γ}}_2{A}_{22}& {\mathrm{\γ}}_2{A}_{21}& -{\mathrm{\γ}}_1{A}_{32}& {\mathrm{\γ}}_1{A}_{31}+1& {\mathrm{\γ}}_2{A}_{23}\\ A_{51}/K_t& A_{52}/K_t& 0& 0& C_1/K_t+1\end{array}\right|,$

K_x和C_x分别为轴段等效圆盘在X方向的刚度和阻尼系数，m为轴段等效圆盘质量；

K_y和C_y分别为轴段等效圆盘在Y方向的刚度和阻尼系数；

A₁₂=A₂₁=0；A₁₃=meA₁，e为轴段等效圆盘的质量偏心距，A₂₃=meB₁， $\left\{\begin{array}{c}{A}_1\\ B_1\end{array}\right\}={\left\{\begin{array}{c}\frac{6\sin \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}+\frac{6\stackrel{\·}{\mathrm{\φ}}\cos \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δt}}+(\stackrel{\·\·}{\mathrm{\φ}}\cos \mathrm{\φ}-{\stackrel{\·}{\mathrm{\φ}}}^2\sin \mathrm{\φ})\\ \frac{6\stackrel{\·}{\mathrm{\φ}}\sin \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δt}}-\frac{6\cos \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}+(\stackrel{\·\·}{\mathrm{\φ}}\sin \mathrm{\φ}+{\stackrel{\·}{\mathrm{\φ}}}^2\cos \mathrm{\φ})\end{array}\right\}}_t,$ φ、

和

分别为轴段转子转过的角度、角速度和角加速度，t为当前时刻；

I_d为轴段等效圆盘转动惯量；

ω为轴段等效圆盘的转动角速度，I_p为轴段等效圆盘极惯矩；

$A_{51}=-\frac{6\mathrm{me}\sin \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2},$ $A_{52}=\frac{6\mathrm{me}\cos \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2},$ $C_1=\frac{6(J+{\mathrm{me}}^2)}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}+\frac{3C_t}{{\mathrm{\θ}}^3\mathrm{\Δ}t}-\mathrm{me}(\stackrel{\·\·}{x}\cos \mathrm{\φ}+\stackrel{\·\·}{y}\sin \mathrm{\φ}),$ J为轴段等效抗弯刚度，和

分别为轴段等效圆盘的形心在X方向和Y方向的振动加速度，C_t为扭振阻尼；

$F_{f,i}=\left[\begin{array}{c}{A}_2\\ B_2\\ L{A}_2+A_3\\ L{B}_2+B_3\\ C_2\end{array}\right],$

$\left\{\begin{array}{c}{A}_2\\ B_2\end{array}\right\}={\left\{\begin{array}{c}[\frac{6m}{{\mathrm{\θ}}^2\mathrm{\Δt}}-C_x(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2})]\stackrel{\·}{x}+[\frac{3m}{\mathrm{\θ}}-C_x\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}]\stackrel{\·\·}{x}\\ [\frac{6m}{{\mathrm{\θ}}^2\mathrm{\Δt}}-C_y(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2})]\stackrel{\·}{y}+[\frac{3m}{\mathrm{\θ}}-C_y\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}]\stackrel{\·\·}{y}\end{array}\right\}}_t$

Δf_x为轴段在X方向上的外力增量，Δf_y为轴段在Y方向上的外力增量，

和

分别为轴段等效圆盘的形心在X方向和Y方向的振动速度，

和分别为轴段等效圆盘的形心在X方向和Y方向的振动加速度；

$\left\{\begin{array}{c}{A}_3\\ B_3\end{array}\right\}=I_d{\left\{\begin{array}{c}\frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}{\stackrel{\·}{\mathrm{\α}}}_x+\frac{3}{\mathrm{\θ}}{\stackrel{\·\·}{\mathrm{\α}}}_x\\ \frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}{\stackrel{\·}{\mathrm{\α}}}_x+\frac{3}{\mathrm{\θ}}{\stackrel{\·\·}{\mathrm{\α}}}_x\end{array}\right\}}_t-\mathrm{\ω}\left[\begin{array}{cc}0& I_p\\ {-I}_p& 0\end{array}\right]{\left\{\begin{array}{c}(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2}){\stackrel{\·}{\mathrm{\α}}}_x+\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}{\stackrel{\·\·}{\mathrm{\α}}}_x\\ (\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2}){\stackrel{\·}{\mathrm{\α}}}_x+\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}{\stackrel{\·\·}{\mathrm{\α}}}_x\end{array}\right\}}_t-\left\{\begin{array}{c}{\mathrm{\ΔL}}_x\\ {\mathrm{\ΔL}}_y\end{array}\right\},$ ${\stackrel{\·}{\mathrm{\α}}}_x$ 和

分别为弯矩和剪力共同作用下与轴段在X方向上的角速度和角加速度，

和

分别为弯矩和剪力共同作用下与轴段在Y方向上的角速度和角加速度，ω为轴段等效圆盘的转动角速度，ΔL_x为轴段等效圆盘在X方向所受的外部弯矩增量，ΔL_y为轴段等效圆盘在Y方向所受的外部弯矩增量；

$F_{e,i}=\left[\begin{array}{c}{\mathrm{\γ}}_3{A}_2+{\mathrm{\γ}}_2{A}_3\\ {\mathrm{\γ}}_3{B}_2+{\mathrm{\γ}}_2{B}_3\\ {\mathrm{\γ}}_2{A}_2+{\mathrm{\γ}}_1{A}_3\\ {\mathrm{\γ}}_2{B}_2+{\mathrm{\γ}}_1{B}_3\\ C_2/K_t\end{array}\right],$

$C_2=-(J+{\mathrm{me}}^2){(\frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}\stackrel{\·}{\mathrm{\φ}}+\frac{3}{\mathrm{\θ}}\stackrel{\·\·}{\mathrm{\φ}})}_t+C_t{[(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2})\stackrel{\·}{\mathrm{\φ}}+\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}\stackrel{\·\·}{\mathrm{\φ}}]}_{\text{t}}$

$+\mathrm{me}{[\sin \mathrm{\φ}(\frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}\stackrel{\·}{x}+\frac{3}{\mathrm{\θ}}\stackrel{\·\·}{x})-\cos \mathrm{\φ}(\frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}\stackrel{\·}{y}+\frac{3}{\mathrm{\θ}}\stackrel{\·\·}{y})]}_t-{\mathrm{\ΔT}}_L;$

和

分别为轴段等效圆盘的形心在X方向和Y方向的振动速度，

和

分别为轴段等效圆盘的形心在X方向和Y方向的振动加速度，ΔT_L为扭矩增量。

所述利用改进的wilson-θ增量表达式计算时刻t₀+(θ+1)Δt的位移、速度和加速度的计算公式为：

$\left\{\begin{array}{c}\mathrm{\Δ}{\stackrel{\·\·}{q}}_{t+\mathrm{\θt}}={\stackrel{\·\·}{q}}_{t+(\mathrm{\θ}+1)\mathrm{\Δt}}-{\stackrel{\·\·}{q}}_{t+\mathrm{\θ\Δt}}=\frac{6}{{\mathrm{\θ}}^3{\mathrm{\Δt}}^2}(q_{t+(\mathrm{\θ}+1)\mathrm{\Δt}}-q_{t+\mathrm{\θ\Δt}})-\frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}{\stackrel{\·}{q}}_{t+\mathrm{\θ\Δt}}-\frac{3}{\mathrm{\θ}}{\stackrel{\·\·}{q}}_{t+\mathrm{\θ\Δt}}\\ \mathrm{\Δ}{\stackrel{\·}{q}}_{t+\mathrm{\θ\Δt}}={\stackrel{\·}{q}}_{t+(\mathrm{\θ}+1)\mathrm{\Δt}}-{\stackrel{\·}{q}}_{t+\mathrm{\θ\Δt}}=\frac{3}{{\mathrm{\θ}}^3\mathrm{\Δt}}(q_{t+(\mathrm{\θ}+1)\mathrm{\Δt}}-q_{t+\mathrm{\θ\Δt}})+(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2}){\stackrel{\·}{q}}_{t+\mathrm{\θ\Δt}}+\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}{\stackrel{\·\·}{q}}_{t+\mathrm{\θ\Δt}}\\ \mathrm{\Δ}{q}_{t+\mathrm{\θ\Δt}}=q_{t+(\mathrm{\θ}+1)\mathrm{\Δt}}-q_{t+\mathrm{\θ\Δt}}=\frac{1}{{\mathrm{\θ}}^3}(q_{t+(\mathrm{\θ}+1)\mathrm{\Δt}}-q_{t+\mathrm{\θ\Δt}})+(1-\frac{\mathrm{\Δt}}{{\mathrm{\θ}}^2}){\stackrel{\·}{q}}_{t+\mathrm{\θ\Δt}}+\frac{{\mathrm{\Δt}}^2(\mathrm{\θ}-1)}{2\mathrm{\θ}}{\stackrel{\·\·}{q}}_{t+\mathrm{\θ\Δt}}\end{array}\right.$

其中，q_t+θΔt、

和

分别为时刻t₀+θΔt机组各节点的位移、速度和加速度，q_t+(θ+1)Δt、

和

分别为时刻t₀+(θ+1)Δt机组各节点的位移、速度和加速度，Δq_t+θΔt、

和

分别为时刻t₀+(θ+1)Δt机组各节点的位移增量、速度增量和加速度增量。

所述设定值为0.2%。

本发明利用改进的wilson-θ增量表达式计算位移加速度，通过基于Riccati变换的轴系弯扭耦合振动增量传递矩阵计算位移增量，其实质是将逐步积分法的思想与传统传递矩阵法相结合，对非线性方程组中各量进行增量改造，可用于强非线性系统的求解，拓展了传递矩阵法在非线性领域的应用；另外，在计算外力/外力矩时，不但能分析其弯振动态响应特性，同时也能进行扭振动态响应的分析，更全面地揭示复杂转子系统的碰摩故障动力学特性；最后，利用前一时刻的位移、速度和加速度，通过迭代运算计算当前时刻的外力/外力矩，其计算结果更接近真实值。

附图说明

图1是汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法流程图；

图2是汽轮发电机组轴段碰摩分析示意图。

具体实施方式

下面结合附图，对优选实施例作详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。

本发明的目的在于，将逐步积分法的思想与传统传递矩阵法相结合，对非线性方程组中各量进行增量改造，使轴系弯振和扭振参量相互耦合，以解决现有技术存在的不足。为此，在实施本发明前，需要进行一些准备，即将逐步积分法的思想与传统传递矩阵法相结合，对非线性方程组中各量进行增量改造。

首先，对wilson-θ表达式进行增量改造。

wilson-θ法是逐步积分法的一种，其表达式为：

${\stackrel{\·\·}{q}}_{t+\mathrm{\Δt}}=\frac{6}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}(q_{t+\mathrm{\θ\Δt}}-q_t)-\frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}{\stackrel{\·}{q}}_t+(1+\frac{3}{\mathrm{\θ}}){\stackrel{\·\·}{q}}_t---\left(1\right)$

${\stackrel{\·}{q}}_{t+\mathrm{\Δt}}=\frac{3}{{\mathrm{\θ}}^3\mathrm{\Δ}t}(q_{t+\mathrm{\θ\Δt}}-q_t)+(1+\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2}){\stackrel{\·}{q}}_t+\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}{\stackrel{\·\·}{q}}_t---\left(2\right)$

$q_{t+\mathrm{\Δt}}=q_t+\frac{1}{{\mathrm{\θ}}^3}(q_{t+\mathrm{\θ\Δt}}-q_t)+(1-\frac{\mathrm{\Δt}}{{\mathrm{\θ}}^2}){\stackrel{\·}{q}}_t+\frac{\mathrm{\Δ}{t}^2(\mathrm{\θ}-1)}{2\mathrm{\θ}}{\stackrel{\·\·}{q}}_t---\left(3\right)$

式中，q为转子的广义坐标。上式说明，如果已知瞬时t某节点的位移q_t、速度和加速度且还知道t+θΔt时刻该结点的位移q_t+θΔt，则可由此式求出t+Δt时刻该结点的位移q_t+Δt、速度

和加速度对于t+θΔt时刻轴系各节点的位移q_t+θΔt，可用下述基于Riccati变换的增量传递矩阵法求得。

将上式改造为增量表达形式：

$\mathrm{\Δ}{\stackrel{\·\·}{q}}_t={\stackrel{\·\·}{q}}_{t+\mathrm{\Δt}}-{\stackrel{\·\·}{q}}_t=\frac{6}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}(q_{t+\mathrm{\θ\Δt}}-q_t)-\frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}{\stackrel{\·}{q}}_t-\frac{3}{\mathrm{\θ}}{\stackrel{\·\·}{q}}_t---\left(4\right)$

${\mathrm{\Δ}{\stackrel{\·}{q}}_t=\stackrel{\·}{q}}_{t+\mathrm{\Δt}}-{\stackrel{\·}{q}}_t=\frac{3}{{\mathrm{\θ}}^3\mathrm{\Δ}t}(q_{t+\mathrm{\θ\Δt}}-q_t)+(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2}){\stackrel{\·}{q}}_t+\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}{\stackrel{\·\·}{q}}_t---\left(5\right)$

$\mathrm{\Δ}{q}_t=q_{t+\mathrm{\Δt}}-q_t+\frac{1}{{\mathrm{\θ}}^3}(q_{t+\mathrm{\θ\Δt}}-q_t)+(1-\frac{\mathrm{\Δt}}{{\mathrm{\θ}}^2}){\stackrel{\·}{q}}_t+\frac{\mathrm{\Δ}{t}^2(\mathrm{\θ}-1)}{2\mathrm{\θ}}{\stackrel{\·\·}{q}}_t---\left(6\right)$

公式（4）-（6）即为改进的wilson-θ增量表达式，其用于根据时刻t的位移、速度和加速度以及时刻t+Δt的位移增量，计算时刻t+Δt的位移、速度和加速度。

其次，建立基于Riccati变换的轴系弯扭耦合振动增量传递矩阵。

对汽轮发电机组转子系统采用多段集中质量模型，系统被模化为多达数百个集中单元轴段，每个集中单元轴段由刚性薄圆盘和弹性轴段两部分组成。

A、建立轴段弯振增量传递矩阵表达式。

以一集中单元轴段为研究对象，分别建立刚性薄圆盘和弹性轴段两部分增量传递矩阵表达式。

1）刚性薄圆盘。

刚性薄圆盘左右截面的剪力及弯矩分别为Q^L、Q^R和M^L、M^R，由D’Alembert原理可得：

$m\left\{\begin{array}{c}{\stackrel{\·\·}{x}}_c\\ {\stackrel{\·\·}{y}}_c\end{array}\right\}={\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^L-{\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^R-\left[K\right]{\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i-\left[C\right]{\left\{\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·}{y}\end{array}\right\}}_i+\left\{\begin{array}{c}{f}_x\\ f_y\end{array}\right\}---\left(7\right)$

$I_d\left\{\begin{array}{c}{\stackrel{\·\·}{\mathrm{\α}}}_x\\ {\stackrel{\·\·}{\mathrm{\α}}}_y\end{array}\right\}={\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^R-{\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^L-\mathrm{\ω}\left[\begin{array}{cc}0& I_p\\ -I_p& 0\end{array}\right]\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\α}}}_x\\ {\stackrel{\·}{\mathrm{\α}}}_y\end{array}\right\}+\left\{\begin{array}{c}{L}_x\\ L_y\end{array}\right\}---\left(8\right)$

同时，

${\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i^R={\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i^L,$ ${\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_i^R={\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_i^L---\left(9\right)$

式中，m为圆盘质量（轴段等效圆盘的质量），x_c、

和

分别为圆盘质量中心在X方向的振动位移、速度和加速度，y_c、

和

分别为圆盘质量中心在Y方向的振动位移、速度和加速度，Q_x和Q_y分别为圆盘在X方向和Y方向所受的剪切力，K为圆盘刚度，[K]为圆盘刚度矩阵，C为圆盘阻尼系数，[C]为圆盘阻尼系数矩阵，x和

分别为圆盘形心在X方向的振动位移和速度，y和

分别为圆盘形心在Y方向的振动位移和速度，f_x和f_y分别为圆盘在X方向和Y方向所受的外力，L和R分别代表圆盘的左截面和右截面，i代表圆盘所处的轴段，I_d为圆盘转动惯量，I_p为圆盘极惯矩，α_x、和

分别为弯矩和剪力共同作用下与X方向的转角、角速度和角加速度，α_y、和

分别为弯矩和剪力共同作用下与Y方向的转角、角速度和角加速度，M_x和M_y分别为圆盘左右截面在X方向和Y方向的弯矩，ω为圆盘转动角速度，L_x和L_y分别为圆盘在X方向和Y方向所受的外部弯矩。

由于 $\left\{\begin{array}{c}{\stackrel{\·\·}{x}}_c\\ {\stackrel{\·\·}{y}}_c\end{array}\right\}$ 为圆盘质心坐标，它可表示为

$\left\{\begin{array}{c}{\stackrel{\·\·}{x}}_c\\ {\stackrel{\·\·}{y}}_c\end{array}\right\}=\left\{\begin{array}{c}\stackrel{\·\·}{x}\\ \stackrel{\·\·}{y}\end{array}\right\}-e\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\φ}}}^2\cos \mathrm{\φ}+\stackrel{\·\·}{\mathrm{\φ}}\sin \mathrm{\φ}\\ {\stackrel{\·}{\mathrm{\φ}}}^2\sin \mathrm{\φ}-\stackrel{\·\·}{\mathrm{\φ}}\cos \mathrm{\φ}\end{array}\right\}---\left(10\right)$

式（10）中φ为转子转过的角度，满足φ=φ₁+β，其中φ₁为不计扭角时转过的角度，β为扭角,e为轴段等效圆盘的质量偏心距。稳定旋转时，φ₁=ωt+φ₀，加速时

A为升速率，B为初转速。将式（10）代入（7），并整理式（7）、（8）和（9）得：

${\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^R={\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^L-m{\left\{\begin{array}{c}\stackrel{\·\·}{x}\\ \stackrel{\·\·}{y}\end{array}\right\}}_i-\left[K\right]{\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i-\left[C\right]{\left\{\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·}{y}\end{array}\right\}}_i+\mathrm{me}\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\φ}}}^2\cos \mathrm{\φ}+\stackrel{\·\·}{\mathrm{\φ}}\sin \mathrm{\φ}\\ {\stackrel{\·}{\mathrm{\φ}}}^2\sin \mathrm{\φ}-\stackrel{\·\·}{\mathrm{\φ}}\cos \mathrm{\φ}\end{array}\right\}+\left\{\begin{array}{c}{f}_x\\ f_y\end{array}\right\}---\left(11\right)$

${\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^R={\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^L+I_d\left\{\begin{array}{c}{\stackrel{\·\·}{\mathrm{\α}}}_x\\ {\stackrel{\·\·}{\mathrm{\α}}}_y\end{array}\right\}+\mathrm{\ω}\left[\begin{array}{cc}0& I_p\\ -I_p& 0\end{array}\right]\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\α}}}_x\\ {\stackrel{\·}{\mathrm{\α}}}_y\end{array}\right\}-\left\{\begin{array}{c}{L}_x\\ L_y\end{array}\right\}---\left(12\right)$

${\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i^R={\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i^L,$ ${\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_i^R={\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_i^L---\left(13\right)$

上式为非线性方程组，各量用增量Δq=q(t+Δt)-q(t)代替，即可将非线性微分方程转化为关于振动状态增量的线性微分方程，从而可以采用适合于线性振动系统的任何一种方法进行分析。增量线性化具体方法为：

令方程组中某参量为f(x₁,x₂,...,x_n)，则其增量为：

Δf=f[x₁(t+Δt),x₂(t+Δt),...,x_n(t+Δt)]-f[x₁(t),x₂(t),...,x_n(t)]

将等式右侧第一项在t时刻按泰勒级数展开，且略去二阶及以上无穷小项，即可得增量表达式：

$\mathrm{\Δf}=\frac{\∂f(x_1,x_2,...,x_n)}{\∂x_1}{|}_t\·\mathrm{\Δ}{x}_1+\frac{\∂f(x_1,x_2,...{,x}_n)}{\∂x_2}{|}_t\·\mathrm{\Δ}{x}_2+...+\frac{\∂f(x_1,x_2,...{,x}_n)}{\∂x_n}{|}_t\·\mathrm{\Δ}{x}_n$

下面以式(11)右端项

为例说明增量线性化具体方法。令 $f(\mathrm{\φ},\stackrel{\·}{\mathrm{\φ}},\stackrel{\·\·}{\mathrm{\φ}})=\mathrm{me}({\stackrel{\·}{\mathrm{\φ}}}^2\cos \mathrm{\φ}+\stackrel{\·\·}{\mathrm{\φ}}\sin \mathrm{\φ}),$ 则

$f(\mathrm{\φ}(t+\mathrm{\Δt}),\stackrel{\·}{\mathrm{\φ}}(t+\mathrm{\Δt}),\stackrel{\·\·}{\mathrm{\φ}}(t+\mathrm{\Δt}))-f(\mathrm{\φ}\left(t\right),\stackrel{\·}{\mathrm{\φ}}\left(t\right),\stackrel{\·\·}{\mathrm{\φ}}\left(t\right))$ 为其增量，将式中第一项在t时刻按泰勒级数展开，且略去二阶及以上无穷小，可得：

将式（11）-（13）用增量方式表达，同时方程组中的

均用Δq及t时刻各已知量替换，即将式（4）、（5）代入，经整理得：

${\left\{\begin{array}{c}\mathrm{\Δ}{Q}_x\\ \mathrm{\Δ}{Q}_y\end{array}\right\}}_i^R={\left\{\begin{array}{c}\mathrm{\Δ}{Q}_x\\ \mathrm{\Δ}{Q}_y\end{array}\right\}}_i^L+\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]{\left\{\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right\}}_i+\left\{\begin{array}{c}{A}_{13}\\ A_{23}\end{array}\right\}\mathrm{\Δ\φ}+\left\{\begin{array}{c}{A}_2\\ B_2\end{array}\right\}---\left(14\right)$

${\left\{\begin{array}{c}\mathrm{\Δ}{M}_x\\ \mathrm{\Δ}{M}_y\end{array}\right\}}_i^R={\left\{\begin{array}{c}\mathrm{\Δ}{M}_x\\ \mathrm{\Δ}{M}_y\end{array}\right\}}_i^L+\left[\begin{array}{cc}{A}_{31}& A_{32}\\ A_{41}& A_{42}\end{array}\right]\left\{\begin{array}{c}\mathrm{\Δ}{\mathrm{\α}}_x\\ \mathrm{\Δ}{\mathrm{\α}}_y\end{array}\right\}+\left\{\begin{array}{c}{A}_3\\ B_3\end{array}\right\}---\left(15\right)$

${\left\{\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right\}}_i^R={\left\{\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right\}}_i^L,$ ${\left\{\begin{array}{c}{\mathrm{\Δ\α}}_x\\ {\mathrm{\Δ\α}}_y\end{array}\right\}}_i^R={\left\{\begin{array}{c}{\mathrm{\Δ\α}}_x\\ {\mathrm{\Δ\α}}_y\end{array}\right\}}_i^L---\left(16\right)$

式（14）-（16）即为刚性薄圆盘的弯振增量传递方程，式中， $A_{11}=-(K_x+C_x\frac{3}{{\mathrm{\θ}}^3\mathrm{\Δt}}+\frac{6m}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}),$ $A_{22}=-(K_y+C_y\frac{3}{{\mathrm{\θ}}^3\mathrm{\Δt}}+\frac{6m}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}),$ A₁₂=A₂₁=0，A₁₃=meA₁，A₂₃=meB₁，

A₄₁=-A₃₂，A₄₂=A₃₁，K_x和C_x分别为轴段等效圆盘在X方向的刚度和阻尼系数，m为轴段等效圆盘质量，β为轴段等效圆盘扭角，K_y和C_y分别为轴段等效圆盘在Y方向的刚度和阻尼系数；A₁₂=A₂₁=0；A₁₃=meA₁，e为轴段等效圆盘的质量偏心距。

$\left\{\begin{array}{c}{A}_1\\ B_1\end{array}\right\}={\left\{\begin{array}{c}\frac{6\sin \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}+\frac{6\stackrel{\·}{\mathrm{\φ}}\cos \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δt}}+(\stackrel{\·\·}{\mathrm{\φ}}\cos \mathrm{\φ}-{\stackrel{\·}{\mathrm{\φ}}}^2\sin \mathrm{\φ})\\ \frac{6\stackrel{\·}{\mathrm{\φ}}\sin \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δt}}-\frac{6\cos \mathrm{\φ}}{{\mathrm{\θ}}^3\mathrm{\Δ}{t}^2}+(\stackrel{\·\·}{\mathrm{\φ}}\sin \mathrm{\φ}+{\stackrel{\·}{\mathrm{\φ}}}^2\cos \mathrm{\φ})\end{array}\right\}}_t,$

$\left\{\begin{array}{c}{A}_2\\ B_2\end{array}\right\}={\left\{\begin{array}{c}[\frac{6m}{{\mathrm{\θ}}^2\mathrm{\Δt}}-C_x(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2})]\stackrel{\·}{x}+[\frac{3m}{\mathrm{\θ}}-C_x\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}]\stackrel{\·\·}{x}\\ [\frac{6m}{{\mathrm{\θ}}^2\mathrm{\Δt}}-C_y(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2})]\stackrel{\·}{y}+[\frac{3m}{\mathrm{\θ}}-C_y\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}]\stackrel{\·\·}{y}\end{array}\right\}}_t$

$\left\{\begin{array}{c}{A}_3\\ B_3\end{array}\right\}=I_d{\left\{\begin{array}{c}\frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}{\stackrel{\·}{\mathrm{\α}}}_x+\frac{3}{\mathrm{\θ}}{\stackrel{\·\·}{\mathrm{\α}}}_x\\ \frac{6}{{\mathrm{\θ}}^2\mathrm{\Δt}}{\stackrel{\·}{\mathrm{\α}}}_x+\frac{3}{\mathrm{\θ}}{\stackrel{\·\·}{\mathrm{\α}}}_x\end{array}\right\}}_t-\mathrm{\ω}\left[\begin{array}{cc}0& I_p\\ {-I}_p& 0\end{array}\right]{\left\{\begin{array}{c}(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2}){\stackrel{\·}{\mathrm{\α}}}_x+\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}{\stackrel{\·\·}{\mathrm{\α}}}_x\\ (\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\θ}}^2}){\stackrel{\·}{\mathrm{\α}}}_x+\frac{\mathrm{\Δt}(\mathrm{\θ}-3)}{2\mathrm{\θ}}{\stackrel{\·\·}{\mathrm{\α}}}_x\end{array}\right\}}_t-\left\{\begin{array}{c}{\mathrm{\ΔL}}_x\\ {\mathrm{\ΔL}}_y\end{array}\right\},$

Δf_x为轴段在X方向上的外力增量，Δf_y为轴段在Y方向上的外力增量，

和

分别为轴段等效圆盘的形心在X方向和Y方向的振动速度，

和

分别为轴段等效圆盘的形心在X方向和Y方向的振动加速度，t为当前时刻，

和分别为弯矩和剪力共同作用下与轴段在X方向上的角速度和角加速度，

和分别为弯矩和剪力共同作用下与轴段在Y方向上的角速度和角加速度，ω为轴段等效圆盘的转动角速度，ΔL_x为轴段等效圆盘在X方向所受的外部弯矩增量，ΔL_y为轴段等效圆盘在Y方向所受的外部弯矩增量。

2）弹性轴段。

设该轴段编号为i，两端截面的编号分别为i及i+1，则针对第i轴段（只有弹性、无质量）的推导与一般的弯振场矩阵完全相同，有：

${\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_{i+1}^L={\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^R,$ ${\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_{i+1}^L={\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^R+L_i{\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^R---\left(17\right)$

${\left\{\begin{array}{c}x\\ y\end{array}\right\}}_{i+1}^L={\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i^R+L_i{\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_i^R+\frac{{L_i}^2}{2\mathrm{EJ}}{\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^R+\frac{{L_i}^3}{6\mathrm{EJ}}{\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^R---\left(18\right)$

${\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_{i+1}^L={\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_i^R+\frac{L_i}{\mathrm{EJ}}{\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^R+\frac{{L_i}^2}{6\mathrm{EJ}}{\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^R---\left(19\right)$

同样写成增量表达式，即为：

${\left\{\begin{array}{c}\mathrm{\Δ}{Q}_x\\ {\mathrm{\ΔQ}}_y\end{array}\right\}}_{i+1}^L={\left\{\begin{array}{c}{\mathrm{\ΔQ}}_x\\ \mathrm{\Δ}{Q}_y\end{array}\right\}}_i^R,$ ${\left\{\begin{array}{c}{\mathrm{\ΔM}}_x\\ {\mathrm{\ΔM}}_y\end{array}\right\}}_{i+1}^L={\left\{\begin{array}{c}\mathrm{\Δ}{M}_x\\ {\mathrm{\ΔM}}_y\end{array}\right\}}_i^R+L_i{\left\{\begin{array}{c}\mathrm{\Δ}{Q}_x\\ {\mathrm{\ΔQ}}_y\end{array}\right\}}_i^R---\left(20\right)$

${\left\{\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right\}}_{i+1}^L={\left\{\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right\}}_i^R+L_i{\left\{\begin{array}{c}{\mathrm{\Δ\α}}_x\\ {\mathrm{\Δ\α}}_y\end{array}\right\}}_i^R+{{\mathrm{\γ}}_2\left\{\begin{array}{c}{\mathrm{\ΔM}}_x\\ {\mathrm{\ΔM}}_y\end{array}\right\}}_i^R+{{\mathrm{\γ}}_3\left\{\begin{array}{c}\mathrm{\Δ}{Q}_x\\ {\mathrm{\ΔQ}}_y\end{array}\right\}}_i^R---\left(21\right)$

${\left\{\begin{array}{c}{\mathrm{\Δ\α}}_x\\ \mathrm{\Δ}{\mathrm{\α}}_y\end{array}\right\}}_{i+1}^L={\left\{\begin{array}{c}\mathrm{\Δ}{\mathrm{\α}}_x\\ {\mathrm{\Δ\α}}_y\end{array}\right\}}_i^R+{{\mathrm{\γ}}_1\left\{\begin{array}{c}{\mathrm{\ΔM}}_x\\ \mathrm{\Δ}{M}_y\end{array}\right\}}_i^R+{{\mathrm{\γ}}_2\left\{\begin{array}{c}{\mathrm{\ΔQ}}_x\\ {\mathrm{\ΔQ}}_y\end{array}\right\}}_i^R---\left(22\right)$

式中， ${\mathrm{\γ}}_1=\frac{L}{\mathrm{EJ}},$ ${\mathrm{\γ}}_2=\frac{L^2}{2\mathrm{EJ}},$ ${\mathrm{\γ}}_3=\frac{L^3}{6\mathrm{EJ}},$ E为轴段弹性模量，J为轴段等效抗弯刚度，L_i=L为第i个轴段的长度。

B、建立轴段扭振增量传递矩阵表达式。

1）刚性薄圆盘。

刚性薄圆盘承受的扭矩有惯性矩、弯振惯性力所产生扭矩、阻力矩、左右截面扭矩差以及外矩，由扭矩所满足的条件，可得：

$J\stackrel{\·\·}{\mathrm{\φ}}=m{\stackrel{\·\·}{x}}_{c}e\sin \mathrm{\φ}-m{\stackrel{\·\·}{y}}_{c}e\cos \mathrm{\φ}-C_t\stackrel{\·}{\mathrm{\φ}}+\mathrm{TL}\left(t\right)+{T_i}^R{-T_i}^L---\left(23\right)$

C_t为扭振阻尼，TL(t)为表示外部扭矩的参量，T_i ^L和T_i ^R分别为第i个圆盘左截面和右截面的扭矩，即

${T_i}^R={T_i}^L+(J+{\mathrm{me}}^2)\stackrel{\·\·}{\mathrm{\φ}}-\mathrm{me}(\stackrel{\·\·}{x}\sin \mathrm{\φ}-\stackrel{\·\·}{y}\cos \mathrm{\φ})+C_t\·\stackrel{\·}{\mathrm{\φ}}-\mathrm{TL}\left(t\right)---\left(24\right)$

改写成增量形式为：

$\mathrm{\Δ}{T_i}^R=\mathrm{\Δ}{T_i}^L+(J+{\mathrm{me}}^2)\mathrm{\Δ}\stackrel{\·\·}{\mathrm{\φ}}+C_t\mathrm{\Δ}\stackrel{\·}{\mathrm{\φ}}$

$-\mathrm{me}[\sin \mathrm{\φ}\·\mathrm{\Δ}\stackrel{\·\·}{x}-\cos \mathrm{\φ}\·\mathrm{\Δ}\stackrel{\·\·}{y}+(\stackrel{\·\·}{x}\cos \mathrm{\φ}+\stackrel{\·\·}{y}\sin \mathrm{\φ})\mathrm{\Δ\φ})]-\mathrm{\Δ}{T}_L---\left(25\right)$

将式（4）-（6）代入，经整理可得

$\mathrm{\Δ}{T_i}^R=\mathrm{\Δ}{T_i}^L+C_1\·\mathrm{\Δ\φ}+A_{51}\mathrm{\Δx}+A_{51}\mathrm{\Δx}+A_{52}\mathrm{\Δy}+C_2---\left(26\right)$

式中， $C_1=\frac{6(J+{\mathrm{me}}^2)}{{\mathrm{\β}}^3\mathrm{\Δ}{t}^2}+\frac{{3C}_t}{{\mathrm{\β}}^3\mathrm{\Δt}}-\mathrm{me}(\stackrel{\·\·}{x}\cos \mathrm{\φ}+\stackrel{\·\·}{y}\sin \mathrm{\φ}),$ $A_{51}=-\frac{6\mathrm{me}\sin \mathrm{\φ}}{{\mathrm{\β}}^3\mathrm{\Δ}{t}^2},$ $A_{52}=\frac{6\mathrm{me}\cos \mathrm{\φ}}{{\mathrm{\β}}^3\mathrm{\Δ}{t}^2}.$

$C_2=-(J+{\mathrm{me}}^2){(\frac{6}{{\mathrm{\β}}^2\mathrm{\Δt}}\stackrel{\·}{\mathrm{\φ}}+\frac{3}{\mathrm{\β}}\stackrel{\·\·}{\mathrm{\φ}})}_t+C_t{[(\frac{\mathrm{\Δt}}{2}-\frac{3}{{\mathrm{\β}}^2})\stackrel{\·}{\mathrm{\φ}}+\frac{\mathrm{\Δt}(\mathrm{\β}-3)}{2\mathrm{\β}}\stackrel{\·\·}{\mathrm{\φ}}]}_t$

$+\mathrm{me}{[\sin \mathrm{\φ}(\frac{6}{{\mathrm{\β}}^2\mathrm{\Δt}}\stackrel{\·}{x}+\frac{3}{\mathrm{\β}}\stackrel{\·\·}{x})-\cos \mathrm{\φ}(\frac{6}{{\mathrm{\β}}^2\mathrm{\Δt}}\stackrel{\·}{y}+\frac{3}{\mathrm{\β}}\stackrel{\·\·}{y})]}_t-{\mathrm{\ΔT}}_L$

同时，有

和

分别为第i个圆盘右截面和左截面转过的角度，ΔT_L为扭矩增量。

其增量表达式为：

${\mathrm{\Δ\φ}}_i^R={\mathrm{\Δ\φ}}_i^L---\left(27\right)$

对于弹性轴段，由于无质量等截面弹性轴段的转动惯量忽略不计，它满足下述关系

$\left\{\begin{array}{c}{T}_{i+1}^L={T_i}^R\\ {\mathrm{\φ}}_{i+1}^L={\mathrm{\φ}}_i^R+\frac{{T_i}^R}{K_t}\end{array}\right.---\left(28\right)$

其增量表达式为：

$\left\{\begin{array}{c}{\mathrm{\ΔT}}_{i+1}^L=\mathrm{\Δ}{T_i}^R\\ {\mathrm{\Δ\φ}}_{i+1}^L={\mathrm{\Δ\φ}}_i^R+\frac{\mathrm{\Δ}{T_i}^R}{K_t}\end{array}\right.---\left(29\right)$

K_t为弹性轴段扭转刚度。

C、建立轴段耦合振动增量传递矩阵表达式。

令转子各截面的力向量增量为

位移向量增量为

则刚性薄圆盘的弯扭耦合振动增量传递方程可写成如下矩阵形式：

${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i^R=\left[\begin{array}{cc}{D}_{11}& D_{12}\\ D_{21}& D_{22}\end{array}\right]\·{\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i^L+{\left\{\begin{array}{c}{D}_{13}\\ D_{23}\end{array}\right\}}_i---\left(30\right)$

式中各矩阵元素分别为：D₁₁=D₂₂=I，I表示单元对角矩阵D₂₁=0，D₂₃=0， $D_{12}=\left|\begin{array}{ccccc}{A}_{11}& A_{12}& 0& 0& A_{13}\\ A_{22}& A_{21}& 0& 0& A_{23}\\ 0& 0& A_{31}& A_{32}& 0\\ 0& 0& -A_{32}& A_{31}& 0\\ A_{51}& A_{52}& 0& 0& C_1\end{array}\right|,$ $D_{13}=\left|\begin{array}{c}{A}_2\\ B_2\\ A_3\\ B_3\\ C_2\end{array}\right|,$ 弹性轴段的弯扭耦合振动增量传递方程可写成如下矩阵形式：

${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}^L=\left[\begin{array}{cc}{B}_{11}& B_{12}\\ B_{21}& B_{22}\end{array}\right]{\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i^R+\left\{\begin{array}{c}{B}_{13}\\ B_{23}\end{array}\right\}---\left(31\right)$

式中各矩阵元素分别为：B₁₂=0,B₁₃=B₂₃=0， $B_{11}=\left|\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ L& 0& 1& 0& 0\\ 0& L& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right|,$ $B_{22}=\left|\begin{array}{ccccc}1& 0& L& 0& 0\\ 0& 1& 0& L& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right|,$ $B_{21}=\left|\begin{array}{ccccc}{\mathrm{\γ}}_3& 0& {\mathrm{\γ}}_2& 0& 0\\ 0& {\mathrm{\γ}}_3& 0& {\mathrm{\γ}}_2& 0\\ {\mathrm{\γ}}_2& 0& {\mathrm{\γ}}_1& 0& 0\\ 0& {\mathrm{\γ}}_2& 0& {\mathrm{\γ}}_1& 0\\ 0& 0& 0& 0& \frac{1}{K_t}\end{array}\right|.$

将刚性圆盘和弹性轴段组合为一部件，有：

${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}^L=\left[B\right]\left(\left[D\right]{\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i^L\begin{array}{c}\end{array}+\left\{\begin{array}{c}{D}_{13}\\ D_{23}\end{array}\right\}\right)+\left\{\begin{array}{c}{B}_{13}\\ B_{23}\end{array}\right\}$

$=\left[B\right]\left[D\right]{\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i^L+\left[B\right]\begin{array}{c}\left\{\begin{array}{c}{D}_{13}\\ D_{23}\end{array}\right\}+\left\{\begin{array}{c}{B}_{13}\\ B_{23}\end{array}\right\}\end{array}$

[B]和[D]为两个系数矩阵，上式写成

${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}^L=\left[\begin{array}{cc}{U}_{11}& U_{12}\\ U_{21}& U_{22}\end{array}\right]{\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i^L+\left\{\begin{array}{c}{F}_f\\ F_e\end{array}\right\}---\left(32\right)$

式中，U₁₁=B₁₁，U₂₁=B₂₁。

$U_{12}=\left|\begin{array}{ccccc}{A}_{11}& A_{12}& 0& 0& A_{13}\\ A_{22}& A_{21}& 0& 0& A_{23}\\ {\mathrm{LA}}_{11}& L{A}_{12}& A_{31}& A_{32}& L{A}_{13}\\ L{A}_{22}& L{A}_{21}& -A_{32}& A_{31}& L{A}_{23}\\ A_{51}& A_{52}& 0& 0& C_1\end{array}\right|,$

$U_{22}=\left|\begin{array}{ccccc}{\mathrm{\γ}}_3{A}_{11}+1& {\mathrm{\γ}}_3{A}_{12}& {\mathrm{\γ}}_2{A}_{31}+L& {\mathrm{\γ}}_2{A}_{32}& {\mathrm{\γ}}_3{A}_{13}\\ {\mathrm{\γ}}_3{A}_{22}& {\mathrm{\γ}}_3{A}_{21}+1& -{\mathrm{\γ}}_2{A}_{32}& {\mathrm{\γ}}_2{A}_{31}+L& {\mathrm{\γ}}_3{A}_{23}\\ {\mathrm{\γ}}_2{A}_{11}& {\mathrm{\γ}}_2{A}_{12}& {\mathrm{\γ}}_1{A}_{31}+1& {\mathrm{\γ}}_1{A}_{32}& {\mathrm{\γ}}_2{A}_{13}\\ {\mathrm{\γ}}_2{A}_{22}& {\mathrm{\γ}}_2{A}_{21}& -{\mathrm{\γ}}_1{A}_{32}& {\mathrm{\γ}}_1{A}_{31}+1& {\mathrm{\γ}}_2{A}_{23}\\ A_{51}/K_t& A_{52}/K_t& 0& 0& C_1/K_t+1\end{array}\right|,$

$F_f=\left[\begin{array}{c}{A}_2\\ B_2\\ L{A}_2+A_3\\ L{B}_2+B_3\\ C_2\end{array}\right],$ $F_e=\left[\begin{array}{c}{\mathrm{\γ}}_3{A}_2+{\mathrm{\γ}}_2{A}_3\\ {\mathrm{\γ}}_3{B}_2+{\mathrm{\γ}}_2{B}_3\\ {\mathrm{\γ}}_2{A}_2+{\mathrm{\γ}}_1{A}_3\\ {\mathrm{\γ}}_2{B}_2+{\mathrm{\γ}}_1{B}_3\\ C_2/K_t\end{array}\right].$

D)对耦合振动增量传递矩阵表达式进行Riccati变换。

将式（32）写为：

$\left\{\begin{array}{c}{f}_{i+1}^L=U_{11}{f_i}^L+U_{12,i}\·e_i^L+F_{f,i}\\ e_{i+1}^L=U_{21}{f_i}^L+U_{22,i}\·e_i^L+F_{e,i}\end{array}\right.$

其中，f_i ^L为第i轴段左侧力向量增量；U_12,i和U_22,i分别为第i轴段矩阵系数；为第i轴段左侧位移向量增量，F_f,i和F_e,i分别为第i轴段参量，且U_12,i=U₁₂，U_22,i=U₂₂，F_f,i=F_f，F_e,i=F_e。

引入Riccati变换，令f_i=S_ie_i+P_i，S_i和P_i分别为第i轴段待定的Riccati传递矩阵，则得到S_i和P_i的递推公式及e_i的递推公式

$S_{i+1}=[U_{11}\·S+U_{12}]\·{[U_{21}S+U_{22}]}_i^{-1}---\left(33\right)$

$e_i={[U_{21}\·S+U_{22}]}_i^{-1}\·e_{i+1}-{[U_{21}\·S+U_{22}]}^{-1}{[U_{21}\·P+\mathrm{Fe}]}_i---\left(35\right)$

由轴段左端边界条件可得S₁=0，P₁=0，利用式（33）、（34）则可顺次得到S₂，P₂，S₃，P₃…….，S_N+1，P_N+1，对于右端截面N+1，有

f_N+1=S_N+1·e_N+1+P_N+1     （36）

且f_N+1=0（边界条件），故

$e_{N+1}=-S_{N+1}^{-1}{P}_{N+1}---\left(37\right)$

然后利用式（35），可从右到左算出各e_i(i＝N,N-1,N-2,......,1)，相应也可算出各截面f_i，所求结果即为各截面在t+θΔt时刻各状态矢量增量。

最后，还要确定加速度的计算方法。

在t₀时刻，各节点的位移和速度是运动的初始条件，一般是已知的，但初始加速度则需另外求出，根据式（11）、（12）和（24）可整理出有关加速度的方程组：

$m\stackrel{\·\·}{x}-\mathrm{me}\sin \mathrm{\φ}\stackrel{\·\·}{\mathrm{\φ}}=(Q_x^L-Q_x^R)-k_{x}x-c_x\stackrel{\·}{x}+\mathrm{me}{\stackrel{\·}{\mathrm{\φ}}}^2\cos \mathrm{\φ}+f_x---\left(38\right)$

$m\stackrel{\·\·}{y}+\mathrm{me}\cos \mathrm{\φ}\stackrel{\·\·}{\mathrm{\φ}}=(Q_y^L-Q_y^R)-k_{y}y-c_y\stackrel{\·}{y}+\mathrm{me}{\stackrel{\·}{\mathrm{\φ}}}^2\sin \mathrm{\φ}+f_y---\left(39\right)$

$-\mathrm{me}\sin \mathrm{\φ}\stackrel{\·\·}{x}+\mathrm{me}\cos \mathrm{\φ}\stackrel{\·\·}{y}+(J+{\mathrm{me}}^2)\stackrel{\·\·}{\mathrm{\φ}}=(T_i^R-T_i^L)-C_t\·\stackrel{\·}{\mathrm{\φ}}+\mathrm{TL}---\left(40\right)$

$I_d{\stackrel{\·\·}{\mathrm{\α}}}_x=M_x^R-M_x^L-{\mathrm{\ωI}}_p{\stackrel{\·}{\mathrm{\α}}}_y+L_x---\left(41\right)$

$I_d{\stackrel{\·\·}{\mathrm{\α}}}_y=M_y^R-M_y^L+{\mathrm{\ωI}}_p{\stackrel{\·}{\mathrm{\α}}}_x+L_y---\left(42\right)$

其中，k_x、c_x、k_y和c_y为滑动轴承动力系数，当轴段为非轴颈轴段时，k_x、c_x、k_y和c_y均为0，f_x和f_y分别为圆盘（轴段）在X方向和Y方向所受的外力。对于轴段i有：

${\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^L-{\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}}_i^R={\left(k_1\right)}_{i-1}{\left\{\begin{array}{c}x\\ y\end{array}\right\}}_{i-1}+{\left(k_2\right)}_{i-1}{\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_{i-1}-[{\left(k_1\right)}_{i-1}+{\left(k_1\right)}_i]{\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i+$

$[{\left(k_2\right)}_{i-1}-{\left(k_2\right)}_i]{\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_i+{\left(k_1\right)}_i{\left\{\begin{array}{c}x\\ y\end{array}\right\}}_{i+1}-{\left(k_2\right)}_i{\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_{i+1}---\left(43\right)$

${\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^R-{\left\{\begin{array}{c}{M}_x\\ M_y\end{array}\right\}}_i^L={-\left(k_2\right)}_{i-1}{\left\{\begin{array}{c}x\\ y\end{array}\right\}}_{i-1}{-\left(k_3\right)}_{i-1}{\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_{i-1}+[{\left(k_2\right)}_{i-1}-{\left(k_2\right)}_i]{\left\{\begin{array}{c}x\\ y\end{array}\right\}}_i-$

$[\underset{r=i-1}{\overset{i}{\mathrm{\Σ}}}{(L{k}_2-k_3)}_r{]\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_i+{\left(k_2\right)}_i{\left\{\begin{array}{c}x\\ y\end{array}\right\}}_{i+1}-{\left(k_3\right)}_i{\left\{\begin{array}{c}{\mathrm{\α}}_x\\ {\mathrm{\α}}_y\end{array}\right\}}_{i+1}---\left(44\right)$

其中， ${\left(k_1\right)}_i={\left(\frac{12\mathrm{EJ}}{L^3}\right)}_i,$ ${\left(k_2\right)}_i={\left(\frac{6\mathrm{EJ}}{L^2}\right)}_i,$ ${\left(k_3\right)}_i={\left(\frac{2\mathrm{EJ}}{L}\right)}_i$

对于扭振，则由方程（29），可容易得到

$T_i^R-T_i^L=k_i({\mathrm{\φ}}_{i+1}^L-{\mathrm{\φ}}_i^R)-k_{i-1}({\mathrm{\φ}}_i^L-{\mathrm{\φ}}_{i-1}^R)---\left(45\right)$

k_i为第i个轴段的刚度。将式（43）、（44）和（45）代入式（38）-（42）中，即可得到初始加速度。

在完成上述准备工作后，可以按照图1所示的步骤，进行汽轮发电机组碰摩故障的弯扭耦合振动特性分析。图1是汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法流程图。图1中，本发明提供的方法包括：

步骤1：设定初始时刻t₀和步进长度Δt。

步骤2：采集时刻t₀机组各节点的位移和速度，计算时刻t₀机组各节点的加速度。

由于已知时刻时刻t₀机组各节点的位移和速度，因此可以利用上述公式（38）-（45）计算时刻t₀机组各节点的加速度。

步骤3：令θ＝0。

步骤4：根据时刻t₀+θΔt机组各节点的位移、速度和加速度，计算时刻t₀+θΔt机组各轴段的外力/外力矩。

当获知某时刻机组各节点的位移、速度和加速度后，即可计算该时刻机组各轴段的外力/外力矩。

本发明采用多段集中质量模型对汽轮发电机组轴系进行模化，根据模型中每个轴段受外力和外力矩的不同，可将轴段分为四类：碰摩轴段、轴颈轴段、动叶栅所在轴段和发电机绕组轴段。计算得到四类轴段外力和外力矩后，即可得到基于Riccati变换的轴系弯扭耦合振动增量传递矩阵。

（1）当轴段为碰摩轴段时。

碰摩发生时轴段坐标和碰摩模型如图2所示。图中，o为静子中心，o₁为单圆盘形心初始位置，o₂为单圆盘形心位置，c为单圆盘质心，φ为转过的角度，δ为转静子间隙圆半径，x-o₁-y是建立在系统静止时转子形心的固定坐标系，F_N为碰摩正压力，F_T为切向摩擦力。系统静止时，转子形心与静子形心存在一定的不对中量（对应图2中的oo₁距离），其值为δ₀，旋转过程中，转静子形心距为当r＞δ时，碰摩发生，碰摩力可表示为：

$\left\{\begin{array}{c}{F}_x\\ F_y\end{array}\right\}=-\frac{(r-\mathrm{\δ})k_s}{r}\left[\begin{array}{cc}1& -{\mathrm{\Θ}}_{\mathrm{\ψ}}\mathrm{\μ}\\ {\mathrm{\Θ}}_{\mathrm{\ψ}}\mathrm{\μ}& 1\end{array}\right]\left\{\begin{array}{c}x\\ y-{\mathrm{\δ}}_0\end{array}\right\}---\left(46\right)$

式中，F_x为发生碰摩时沿转子形心X方向的摩擦力，F_x为发生碰摩时沿转子形心Y方向的摩擦力，x为发生碰摩时沿转子形心X方向的位移，y为发生碰摩时沿转子形心Y方向的位移，k_s为定子的径向刚度，μ为转子与静子间的摩擦系数，摩擦力方向由符号Θ_ψ决定：

${\mathrm{\&Theta;}}_{\mathrm{\&psi;}}=\left\{\begin{array}{cc}-1& \mathrm{\&psi;}<0\\ 0& \mathrm{\&psi;}=0\\ 1& \mathrm{\&psi;}>0\end{array}\right.---\left(47\right)$

ψ是动静碰摩点的线速度，它与转子自转角速度和涡动角速度有关，表达式为：

$\mathrm{\&psi;}=\frac{x\stackrel{\&CenterDot;}{y}-(y-{\mathrm{\&delta;}}_0)\stackrel{\&CenterDot;}{x}}{r}+R_0\mathrm{\&omega;}---\left(48\right)$

式中，r₀为碰摩点到转子形心的距离，为发生碰摩时沿转子形心X方向的速度，

为发生碰摩时沿转子形心Y方向的速度。摩擦力的存在，将使定子在碰摩点对转子形心产生摩擦力矩M_F，有：

M_F=-F_TR₀=-Θ_ψμF_NR₀    （49）

（2）当轴段为轴颈轴段时。

机组轴颈轴段在滑动轴承内高速旋转，将受到油膜力的作用。通常油膜力按非线性力处理，本文通过建立非线性油膜力数据库的方法计算非线性油膜力。以椭圆轴承为例作介绍。

椭圆轴承的油膜力是由上瓦和下瓦所受力的合成，即：

$\left\{\begin{array}{c}{F}_X(x,y,\stackrel{\&CenterDot;}{x},\stackrel{\&CenterDot;}{y})=\sqrt{{(1-{2\mathrm{\&theta;}}_1^{\&prime;})}^2+{4\mathrm{\&epsiv;}}_1^{\&prime;2}}{f}_{x10}({\mathrm{\&epsiv;}}_1,{\mathrm{\&theta;}}_1,{\mathrm{\&alpha;}}_1)+\sqrt{{(1-{2\mathrm{\&theta;}}_2^{\&prime;})}^2+{4\mathrm{\&epsiv;}}_2^{\&prime;2}}{f}_{x20}({\mathrm{\&epsiv;}}_2,{\mathrm{\&theta;}}_2,{\mathrm{\&alpha;}}_2)\\ F_Y(x,y,\stackrel{\&CenterDot;}{x},\stackrel{\&CenterDot;}{y})=\sqrt{{(1-{2\mathrm{\&theta;}}_1^{\&prime;})}^2+{4\mathrm{\&epsiv;}}_1^{\&prime;2}}{f}_{y10}({\mathrm{\&epsiv;}}_1,{\mathrm{\&theta;}}_1,{\mathrm{\&alpha;}}_1)+\sqrt{{(1-{2\mathrm{\&theta;}}_2^{\&prime;})}^2+{4\mathrm{\&epsiv;}}_2^{\&prime;2}}{f}_{y20}({\mathrm{\&epsiv;}}_2,{\mathrm{\&theta;}}_2,{\mathrm{\&alpha;}}_2)\end{array}\right.---\left(50\right)$

式中，F_x10、f_y10、f_x20和f_y20即为数据库中存储的当量油膜力分量，它由参数(ε,θ，α)决定；

为轴颈轴段的油膜力在上瓦的分力，

为轴颈轴段的油膜力在下瓦的分力，x和y分别为轴段轴心在X和Y方向的位移，

和

分别为轴段轴心在X和Y方向的速度，ε₁、θ₁、ε’₁、θ’₁和α₁分别是轴段轴心相对于下瓦中心o₁在极坐标系(ε,θ)下的运动参量，ε₂、θ₂、ε’₂、θ’₂和α₂分别是轴段轴心相对于上瓦中心o₂在极坐标系(ε,θ)下的运动参量。

对于运行在一定工况下的，具有一定宽径比及瓦包角的椭圆轴承，根据轴心的运动状态、轴承的椭圆度和公式(50)，即可直接计算椭圆轴承的非线性油膜力分量。

（3）当轴段为动叶栅所在轴段时。

动叶栅所在轴段是指高、中、低压缸内各级处的转子轴段，当机组带一定负荷运行时，动叶栅将会受到蒸汽的弯曲力矩，该力矩即为动叶栅所在轴段所受的外力矩，其值为：

$M=\frac{N_{\mathrm{el}}}{{\mathrm{\&omega;}}_e}---\left(51\right)$

N_el为汽轮机的机械功率输出，ω_e为转子转动的角频率，有：

$N_{\mathrm{el}}=N_e\frac{{\mathrm{\&Delta;H}}_{\mathrm{tl}}}{{\mathrm{\&Delta;H}}_t}\frac{D_l}{D}\frac{{\mathrm{\&eta;}}_{\mathrm{oel}}}{{\mathrm{\&eta;}}_{\mathrm{oe}}}---\left(52\right)$

式中，N_e为机组额定功率，η_oe为机组相对内效率，机组工况变化不大时，通常保持不变，取为常数。η_oel为偏离额定工况时的机组相对内效率，ΔH为机组理想比焓降，有：

$\frac{{\mathrm{\&Delta;H}}_{t1}}{{\mathrm{\&Delta;H}}_t}=\frac{T_{01}[1-{(P_{21}/P_{01})}^{\frac{R-1}{R}}]}{T_0\left[{1-(P_2/P_0)}^{\frac{R-1}{R}}\right]}$

T₀、P₀和P₂分别为级或级组在额定工况下蒸汽入口温度、压力和蒸汽排汽压力，T₀₁、P₀₁和P₂₁分别为级或级组在偏离额定工况时的蒸汽入口温度、压力和蒸汽排汽压力，R为气体常数，D为蒸汽流量，有：

$\frac{D_l}{D}=\sqrt{\frac{P_{H01}^2-P_{H21}^2}{P_{H0}^2-P_{H2}^2}}\sqrt{\frac{T_{H0}}{T_{H01}}}$

T_H0、P_H0和P_H2分别为额定工况下主蒸汽入口温度、压力和乏汽排汽压力，T_H01、P_H01和P_H21分别为偏离额定工况时的主蒸汽入口温度、压力和乏汽排汽压力。

（4）当轴段为发电机绕组轴段时。

发电机绕组轴段受到电磁力矩M_e的作用，其表达式为：

M_e=ψ_di_q-ψ_qi_d=(-x_di_d+x_afI_f+x_aDI_D)i_q-(-x_qi_q+x_aQI_Q)i_d   （53）

i_d和i_q分别为定子绕组纵轴和橫轴的电流，x_d为定子绕组d轴（纵轴）自感系数，x_q为定子绕组q轴（横轴）自感系数，x_aD为定子绕组d轴互感系数，x_aQ为定子绕组q轴互感系数，x_af为定子绕组励磁互感系数，I_f为定子绕组励磁电流，ψ_d和ψ_q分别为定子绕组d轴和q轴的磁链。根据Park变换，有dqo基本方程：

$\left[\begin{array}{c}{i}_d\\ i_q\end{array}\right]=\frac{2}{3}\left[\begin{array}{ccc}\cos \mathrm{\&gamma;}& \cos (\mathrm{\&gamma;}-\frac{2}{3}\mathrm{\&pi;})& \cos (\mathrm{\&gamma;}+\frac{2}{3}\mathrm{\&pi;})\\ -\sin \mathrm{\&gamma;}& -\sin (\mathrm{\&gamma;}-\frac{2}{3}\mathrm{\&pi;})& -\sin (\mathrm{\&gamma;}+\frac{2}{3}\mathrm{\&pi;})\end{array}\right]\left[\begin{array}{c}{I}_a\\ I_b\\ I_c\end{array}\right]---\left(54\right)$

式中，γ为电机转子d轴与定子a相绕组轴线的夹角，I_a、I_b和I_c分别为a、b、c三相电流，有：

${\mathrm{\&gamma;}}^t={\&Integral;}_0^t\mathrm{\&omega;dt}+{\mathrm{\&gamma;}}_0---\left(55\right)$

I_D和I_Q分别为定子绕组和阻尼绕组纵轴和橫轴的电流，ω为发电机每一时步的角速度。

步骤5：将时刻t₀+θΔt机组各轴段的外力/外力矩分别作为时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩。

由于计算某一时刻各个轴段的位移、速度和加速度需要用到外力/外力矩，而相关外力/外力矩获得比较困难，因此本发明使用该时刻的前一时刻的位移、速度和加速度计算出的外力和外力矩（即前一时刻的外力和外力矩）近似作为该时刻的外力/外力矩，再使用迭代过程逼近该时刻的外力/外力矩的真实值。

步骤6：根据时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩，利用基于Riccati变换的轴系弯扭耦合振动增量传递矩阵计算时刻t₀+(θ+1)Δt机组各节点的位移增量。

当获得某时刻各个轴段的外力/外力矩后，可以利用公式 $\left\{\begin{array}{c}{f}_{i+1}^L=U_{11}{f}_i^L+U_{12,i}\&CenterDot;e_i^L+F_{f,i}\\ e_{i+1}^L=U_{21}{f}_i^L+U_{22,i}\&CenterDot;e_i^L+F_{e,i}\end{array}\right.$ 结合公式（33）-（37）计算该时刻机组各节点的位移增量。

步骤7：根据时刻t₀+(θ+1)Δt机组各节点的位移增量，利用改进的wilson-θ增量表达式计算时刻t₀+(θ+1)Δt的位移、速度和加速度。

计算得到位移增量后，利用前一时刻的位移、速度、加速度和该时刻的位移增量，通过公式（4）-（6）即可计算出该时刻的位移、速度和加速度。此时需要注意的是，前一时刻为t₀+θΔt，当前时刻为t₀+(θ+1)Δt，因此公式（4）-（6）变为：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;t}}={\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}=\frac{6}{{\mathrm{\&theta;}}^3{\mathrm{\&Delta;t}}^2}(q_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-q_{t+\mathrm{\&theta;\&Delta;t}})-\frac{6}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}{\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}-\frac{3}{\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}={\stackrel{\&CenterDot;}{q}}_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-{\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}=\frac{3}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;t}}(q_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-q_{t+\mathrm{\&theta;\&Delta;t}})+(\frac{\mathrm{\&Delta;t}}{2}-\frac{3}{{\mathrm{\&theta;}}^2}){\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}+\frac{\mathrm{\&Delta;t}(\mathrm{\&theta;}-3)}{2\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}\\ \mathrm{\&Delta;}{q}_{t+\mathrm{\&theta;\&Delta;t}}=q_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-q_{t+\mathrm{\&theta;\&Delta;t}}=\frac{1}{{\mathrm{\&theta;}}^3}(q_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-q_{t+\mathrm{\&theta;\&Delta;t}})+(1-\frac{\mathrm{\&Delta;t}}{{\mathrm{\&theta;}}^2}){\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}+\frac{{\mathrm{\&Delta;t}}^2(\mathrm{\&theta;}-1)}{2\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}\end{array}\right..$

步骤8：根据时刻t₀+(θ+1)Δt的位移、速度和加速度，计算时刻t₀+(θ+1)Δt机组各轴段的外力/外力矩。

利用步骤4中的方法，使用步骤7中计算出的位移、速度和加速度，计算新的机组各轴段的外力/外力矩。

由于之前得到的临时外力/临时外力矩只是外力/外力矩真实值的近似值，因此本发明通过迭代的方法不断更新外力/外力矩，以逼近其真实值。迭代的过程是使用本步骤计算的外力/外力矩与临时外力/临时外力矩做比较，从而确定本步骤计算的外力/外力矩是否接近真实值。

步骤9：判断是否满足设定条件，如果满足设定条件，则执行步骤11；否则，执行步骤10。

设定条件为：时刻t₀+(θ+1)Δt机组各轴段的外力/外力矩与时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩的差值绝对值小于设定值。

将每一次迭代获得的临时外力/临时外力矩（步骤5）和计算的外力/外力矩（步骤8）作差，判断其差值的绝对值是否小于设定值，如果其差值的绝对值小于设定值，说明计算的外力/外力矩已经非常接近真实值，此时可以结束迭代过程。设定值一般取0.2%。

步骤10：临时外力/临时外力矩和计算的外力/外力矩的差值绝对值不小于设定值，说明此时计算的外力/外力矩还不够接近真实值，此时令机组各轴段的临时外力/临时外力矩等于步骤8计算的机组各轴段的外力/外力矩，返回步骤6，继续迭代过程。

步骤11：当临时外力/临时外力矩和计算的外力/外力矩的差值绝对值小于设定值时，说明步骤8计算的外力/外力矩已经非常接近真实值。此时，可以将步骤8计算的外力/外力矩作为该时刻的外力/外力矩真实值，且将步骤8中计算外力/外力矩时所用到的各轴段的位移、速度和加速度作为该时刻各轴段的位移、速度和加速度的真实值。

令θ=θ+1，返回步骤5，继续循环执行上述过程，则可以得到t₀+Δt，t₀+2Δt，t₀+3Δt，...等时刻的外力/外力矩、位移、速度和加速度的值。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (5)

1.一种汽轮发电机组碰摩故障的弯扭耦合振动特性分析方法，其特征是所述方法包括：

步骤1：设定初始时刻t₀和步进长度Δt；

步骤2：采集时刻t₀机组各节点的位移和速度，计算时刻t₀机组各节点的加速度;

步骤3：令θ＝0；

步骤4：根据时刻t₀+θΔt机组各节点的位移、速度和加速度，计算时刻t₀+θΔt机组各轴段的外力/外力矩；

步骤5：将时刻t₀+θΔt机组各轴段的外力/外力矩分别作为时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩；

步骤6：根据时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩，利用基于Riccati变换的轴系弯扭耦合振动增量传递矩阵计算时刻t₀+(θ+1)Δt机组各节点的位移增量；

步骤7：根据时刻t₀+(θ+1)Δt机组各节点的位移增量，利用改进的wilson-θ增量表达式计算时刻t₀+(θ+1)Δt的位移、速度和加速度；

步骤8：根据时刻t₀+(θ+1)Δt的位移、速度和加速度，计算时刻t₀+(θ+1)Δt机组各轴段的外力/外力矩；

步骤9：判断是否满足设定条件，如果满足设定条件，则执行步骤11；否则，执行步骤10；

所述设定条件为：时刻t₀+(θ+1)Δt机组各轴段的外力/外力矩与时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩的差值绝对值小于设定值；

步骤10：令时刻t₀+(θ+1)Δt机组各轴段的临时外力/临时外力矩等于时刻t₀+(θ+1)Δt机组各轴段的外力/外力矩，返回步骤6；

步骤11：令θ=θ+1，返回步骤5。

2.根据权利要求1所述的方法，其特征是所述计算时刻t₀+θΔt机组各轴段的外力/外力矩具体为，当轴段为碰摩轴段时，计算时刻t₀+θΔt机组该轴段的摩擦力矩；当轴段为轴颈轴段时，计算时刻t₀+θΔt机组该轴段的油膜力；当轴段为动叶栅所在轴段时，计算时刻t₀+θΔt机组动叶栅受到的蒸汽弯曲力矩；当轴段为发电机绕组轴段时；计算时刻t₀+θΔt机组的发电机绕组轴段受到的电磁力矩。

3.根据权利要求1所述的方法，其特征是所述利用基于Riccati变换的轴系弯扭耦合振动增量传递矩阵计算时刻t₀+(θ+1)Δt机组各节点的位移增量的计算公式为：

$\left\{\begin{array}{c}{f}_{i+1}^L=U_{11}{f_i}^L+U_{12,i}\&CenterDot;e_i^L+F_{f,i}\\ e_{i+1}^L=U_{21}{f_i}^L+U_{22,i}\&CenterDot;e_i^L+F_{e,i}\end{array}\right.,$

其中，f_i ^L为第i个轴段左侧力向量增量，为第i个轴段左侧位移向量增量，U₁₁和U₂₁分别为基于Riccati变换的轴系弯扭耦合振动增量传递矩阵参数，U_12,i、U_22,i、F_f,i和F_e,i分别为第i个轴段的基于Riccati变换的轴系弯扭耦合振动增量传递矩阵参数；

$U_{11}=\left|\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ L& 0& 1& 0& 0\\ 0& L& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right|,$ L轴段的长度；

$U_{21}=\left|\begin{array}{ccccc}{\mathrm{\&gamma;}}_3& 0& {\mathrm{\&gamma;}}_2& 0& 0\\ 0& {\mathrm{\&gamma;}}_3& 0& {\mathrm{\&gamma;}}_2& 0\\ {\mathrm{\&gamma;}}_2& 0& {\mathrm{\&gamma;}}_1& 0& 0\\ 0& {\mathrm{\&gamma;}}_2& 0& {\mathrm{\&gamma;}}_1& 0\\ 0& 0& 0& 0& \frac{1}{K_t}\end{array}\right|,$ K_t为弹性轴段扭转刚度， ${\mathrm{\&gamma;}}_1=\frac{L}{\mathrm{EJ}},$ ${\mathrm{\&gamma;}}_2=\frac{L^2}{2\mathrm{EJ}},$

E为轴段的弹性模量，J为轴段等效抗弯刚度；

$U_{12,i}=\left|\begin{array}{ccccc}{A}_{11}& A_{12}& 0& 0& A_{13}\\ A_{22}& A_{21}& 0& 0& A_{23}\\ {\mathrm{LA}}_{11}& L{A}_{12}& A_{31}& A_{32}& L{A}_{13}\\ L{A}_{22}& L{A}_{21}& -A_{32}& A_{31}& L{A}_{23}\\ A_{51}& A_{52}& 0& 0& C_1\end{array}\right|,$

$U_{22,i}=\left|\begin{array}{ccccc}{\mathrm{\&gamma;}}_3{A}_{11}+1& {\mathrm{\&gamma;}}_3{A}_{12}& {\mathrm{\&gamma;}}_2{A}_{31}+L& {\mathrm{\&gamma;}}_2{A}_{32}& {\mathrm{\&gamma;}}_3{A}_{13}\\ {\mathrm{\&gamma;}}_3{A}_{22}& {\mathrm{\&gamma;}}_3{A}_{21}+1& -{\mathrm{\&gamma;}}_2{A}_{32}& {\mathrm{\&gamma;}}_2{A}_{31}+L& {\mathrm{\&gamma;}}_3{A}_{23}\\ {\mathrm{\&gamma;}}_2{A}_{11}& {\mathrm{\&gamma;}}_2{A}_{12}& {\mathrm{\&gamma;}}_1{A}_{31}+1& {\mathrm{\&gamma;}}_1{A}_{32}& {\mathrm{\&gamma;}}_2{A}_{13}\\ {\mathrm{\&gamma;}}_2{A}_{22}& {\mathrm{\&gamma;}}_2{A}_{21}& -{\mathrm{\&gamma;}}_1{A}_{32}& {\mathrm{\&gamma;}}_1{A}_{31}+1& {\mathrm{\&gamma;}}_2{A}_{23}\\ A_{51}/K_t& A_{52}/K_t& 0& 0& C_1/K_t+1\end{array}\right|,$

K_x和C_x分别为轴段等效圆盘在X方向的刚度和阻尼系数，m为轴段等效圆盘质量；

K_y和C_y分别为轴段等效圆盘在Y方向的刚度和阻尼系数；

A₁₂=A₂₁=0；A₁₃=meA₁，e为轴段等效圆盘的质量偏心距，A₂₃=meB₁， $\left\{\begin{array}{c}{A}_1\\ B_1\end{array}\right\}={\left\{\begin{array}{c}\frac{6\sin \mathrm{\&phi;}}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;}{t}^2}+\frac{6\stackrel{\&CenterDot;}{\mathrm{\&phi;}}\cos \mathrm{\&phi;}}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;t}}+(\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}\cos \mathrm{\&phi;}-{\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}^2\sin \mathrm{\&phi;})\\ \frac{6\stackrel{\&CenterDot;}{\mathrm{\&phi;}}\sin \mathrm{\&phi;}}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;t}}-\frac{6\cos \mathrm{\&phi;}}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;}{t}^2}+(\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}\sin \mathrm{\&phi;}+{\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}^2\cos \mathrm{\&phi;})\end{array}\right\}}_t,$ φ、

和

分别为轴段转子转过的角度、角速度和角加速度，t为当前时刻；

I_d为轴段等效圆盘转动惯量；

ω为轴段等效圆盘的转动角速度，I_p为轴段等效圆盘极惯矩；

$A_{51}=-\frac{6\mathrm{me}\sin \mathrm{\&phi;}}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;}{t}^2},$ $A_{52}=\frac{6\mathrm{me}\cos \mathrm{\&phi;}}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;}{t}^2},$ $C_1=\frac{6(J+{\mathrm{me}}^2)}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;}{t}^2}+\frac{3C_t}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;}t}-\mathrm{me}(\stackrel{\&CenterDot;\&CenterDot;}{x}\cos \mathrm{\&phi;}+\stackrel{\&CenterDot;\&CenterDot;}{y}\sin \mathrm{\&phi;}),$ J为轴段等效抗弯刚度，

和

分别为轴段等效圆盘的形心在X方向和Y方向的振动加速度，C_t为扭振阻尼；

$F_{f,i}=\left[\begin{array}{c}{A}_2\\ B_2\\ L{A}_2+A_3\\ L{B}_2+B_3\\ C_2\end{array}\right],$

$\left\{\begin{array}{c}{A}_2\\ B_2\end{array}\right\}={\left\{\begin{array}{c}[\frac{6m}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}-C_x(\frac{\mathrm{\&Delta;t}}{2}-\frac{3}{{\mathrm{\&theta;}}^2})]\stackrel{\&CenterDot;}{x}+[\frac{3m}{\mathrm{\&theta;}}-C_x\frac{\mathrm{\&Delta;t}(\mathrm{\&theta;}-3)}{2\mathrm{\&theta;}}]\stackrel{\&CenterDot;\&CenterDot;}{x}\\ [\frac{6m}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}-C_y(\frac{\mathrm{\&Delta;t}}{2}-\frac{3}{{\mathrm{\&theta;}}^2})]\stackrel{\&CenterDot;}{y}+[\frac{3m}{\mathrm{\&theta;}}-C_y\frac{\mathrm{\&Delta;t}(\mathrm{\&theta;}-3)}{2\mathrm{\&theta;}}]\stackrel{\&CenterDot;\&CenterDot;}{y}\end{array}\right\}}_t$

Δf_x为轴段在X方向上的外力增量，Δf_y为轴段在Y方向上的外力增量，

和分别为轴段等效圆盘的形心在X方向和Y方向的振动速度，

和

分别为轴段等效圆盘的形心在X方向和Y方向的振动加速度；

$\left\{\begin{array}{c}{A}_3\\ B_3\end{array}\right\}=I_d{\left\{\begin{array}{c}\frac{6}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_x+\frac{3}{\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&alpha;}}}_x\\ \frac{6}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_x+\frac{3}{\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&alpha;}}}_x\end{array}\right\}}_t-\mathrm{\&omega;}\left[\begin{array}{cc}0& I_p\\ {-I}_p& 0\end{array}\right]{\left\{\begin{array}{c}(\frac{\mathrm{\&Delta;t}}{2}-\frac{3}{{\mathrm{\&theta;}}^2}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_x+\frac{\mathrm{\&Delta;t}(\mathrm{\&theta;}-3)}{2\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&alpha;}}}_x\\ (\frac{\mathrm{\&Delta;t}}{2}-\frac{3}{{\mathrm{\&theta;}}^2}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_x+\frac{\mathrm{\&Delta;t}(\mathrm{\&theta;}-3)}{2\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&alpha;}}}_x\end{array}\right\}}_t-\left\{\begin{array}{c}{\mathrm{\&Delta;L}}_x\\ {\mathrm{\&Delta;L}}_y\end{array}\right\},$ ${\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_x$ 和

分别为弯矩和剪力共同作用下与轴段在X方向上的角速度和角加速度，

和

分别为弯矩和剪力共同作用下与轴段在Y方向上的角速度和角加速度，ω为轴段等效圆盘的转动角速度，ΔL_x为轴段等效圆盘在X方向所受的外部弯矩增量，ΔL_y为轴段等效圆盘在Y方向所受的外部弯矩增量；

$F_{e,i}=\left[\begin{array}{c}{\mathrm{\&gamma;}}_3{A}_2+{\mathrm{\&gamma;}}_2{A}_3\\ {\mathrm{\&gamma;}}_3{B}_2+{\mathrm{\&gamma;}}_2{B}_3\\ {\mathrm{\&gamma;}}_2{A}_2+{\mathrm{\&gamma;}}_1{A}_3\\ {\mathrm{\&gamma;}}_2{B}_2+{\mathrm{\&gamma;}}_1{B}_3\\ C_2/K_t\end{array}\right],$

$C_2=-(J+{\mathrm{me}}^2){(\frac{6}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}\stackrel{\&CenterDot;}{\mathrm{\&phi;}}+\frac{3}{\mathrm{\&theta;}}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}})}_t+C_t{[(\frac{\mathrm{\&Delta;t}}{2}-\frac{3}{{\mathrm{\&theta;}}^2})\stackrel{\&CenterDot;}{\mathrm{\&phi;}}+\frac{\mathrm{\&Delta;t}(\mathrm{\&theta;}-3)}{2\mathrm{\&theta;}}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}]}_{\text{t}}$

$+\mathrm{me}{[\sin \mathrm{\&phi;}(\frac{6}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}\stackrel{\&CenterDot;}{x}+\frac{3}{\mathrm{\&theta;}}\stackrel{\&CenterDot;\&CenterDot;}{x})-\cos \mathrm{\&phi;}(\frac{6}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}\stackrel{\&CenterDot;}{y}+\frac{3}{\mathrm{\&theta;}}\stackrel{\&CenterDot;\&CenterDot;}{y})]}_t-{\mathrm{\&Delta;T}}_L;$

和

分别为轴段等效圆盘的形心在X方向和Y方向的振动速度，和

分别为轴段等效圆盘的形心在X方向和Y方向的振动加速度，ΔT_L为扭矩增量。

4.根据权利要求1所述的方法，其特征是所述利用改进的wilson-θ增量表达式计算时刻t₀+(θ+1)Δt的位移、速度和加速度的计算公式为：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;t}}={\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}=\frac{6}{{\mathrm{\&theta;}}^3{\mathrm{\&Delta;t}}^2}(q_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-q_{t+\mathrm{\&theta;\&Delta;t}})-\frac{6}{{\mathrm{\&theta;}}^2\mathrm{\&Delta;t}}{\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}-\frac{3}{\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}={\stackrel{\&CenterDot;}{q}}_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-{\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}=\frac{3}{{\mathrm{\&theta;}}^3\mathrm{\&Delta;t}}(q_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-q_{t+\mathrm{\&theta;\&Delta;t}})+(\frac{\mathrm{\&Delta;t}}{2}-\frac{3}{{\mathrm{\&theta;}}^2}){\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}+\frac{\mathrm{\&Delta;t}(\mathrm{\&theta;}-3)}{2\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}\\ \mathrm{\&Delta;}{q}_{t+\mathrm{\&theta;\&Delta;t}}=q_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-q_{t+\mathrm{\&theta;\&Delta;t}}=\frac{1}{{\mathrm{\&theta;}}^3}(q_{t+(\mathrm{\&theta;}+1)\mathrm{\&Delta;t}}-q_{t+\mathrm{\&theta;\&Delta;t}})+(1-\frac{\mathrm{\&Delta;t}}{{\mathrm{\&theta;}}^2}){\stackrel{\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}+\frac{{\mathrm{\&Delta;t}}^2(\mathrm{\&theta;}-1)}{2\mathrm{\&theta;}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{t+\mathrm{\&theta;\&Delta;t}}\end{array}\right.$

其中，q_t+θΔt、

和

分别为时刻t₀+θΔt机组各节点的位移、速度和加速度，q_t+(θ+1)Δt、

和

分别为时刻t₀+(θ+1)Δt机组各节点的位移、速度和加速度，Δq_t+θΔt、

和分别为时刻t₀+(θ+1)Δt机组各节点的位移增量、速度增量和加速度增量。

5.根据权利要求1所述的方法，其特征是所述设定值为0.2%。

Images