CN104236705A - 汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法及系统 - Google Patents

Abstract

本发明属于汽轮发电机组轴系扭振监测技术领域，尤其涉及一种汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法及系统。该方法首先对机电网协调作用下电力系统大扰动导致的轴系扭振进行仿真计算，建立电力系统大扰动下发电机电磁力矩最大值与轴系疲劳寿命损耗之间的对应关系，制定损伤报警阈值和跳机保护阈值，然后根据发电机三相电流和电压计算发电机电磁力矩，在机组安全性受到威胁时及时地发出损伤报警或跳机保护信号，最后再对扭振故障进行轴系扭振动态响应计算，进而对轴系做出精确的损伤评估和剩余寿命评估。该方法可实现针对电力系统大扰动故障的扭振在线监测、分析与保护，同时可对机组轴系扭振疲劳寿命损耗进行精确计算。

Description

汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法及系统

技术领域

本发明属于汽轮发电机组轴系扭振监测技术领域，尤其涉及一种用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法及系统。

背景技术

汽轮发电机组轴系扭振根据它被激发机理的不同可以分为瞬态冲击类扭振和次同步振荡等共振类扭振。自动重合闸、非同期并网、瞬时性对称与不对称短路及线路开关切合操作等突发性机电扰动，将有可能产生短时间性冲击扭矩，形成短时间冲击性轴系扭振。这种扭振虽然时间不长，但在不利的配合条件下，仍有可能产生高到机组轴系无法承受的交变扭应力而造成一次性破坏，或者产生扭振疲劳损耗积累。

电力系统大扰动故障下汽轮发电机组轴系发生瞬态冲击类扭振时，其扭振振幅表现出突然增大后又逐渐衰减的特点，根据监测机组机头和机尾测速齿轮所得的轴系扭振相对角位移信号难以用带通滤波准确地滤出，因此不宜采用模态叠加法进行进一步分析计算。同时，电力系统大扰动故障下汽轮发电机组轴系瞬态冲击类扭振危害性很大，如果不能在发现该类扭振时作出及时的损伤报警或跳机保护动作，可能造成机组轴系极大的损伤。现有的监测与分析手段不能够快速、准确的对轴系扭振进行定位与分析，方法也需要进一步完善和改进。

根据瞬态冲击类扭振产生的机理，及其引发的机组轴系扭振所表现的特征，同时考虑实时在线监测、分析和保护的快速性与准确性，采用模型仿真法对轴系扭振响应进行求解，分析计算得到其对机组轴系的损伤程度。把该方法应用于汽轮发电机组轴系扭振在线监测与、分析与保护系统，得到很好的应用效果。

发明内容

针对上述问题，本发明提出了一种用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法及系统。

一种用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法，包括：

步骤1：通过电力系统仿真，模拟电力系统不同类型的大扰动故障下机组发生瞬态冲击类扭振时发电机的电磁力矩，仿真得到不同类型瞬态冲击类扭振故障下不同幅值的电磁力矩响应；

步骤2：根据仿真所得电磁力矩响应计算轴系扭振危险截面局部扭应力响应；

步骤3：对机组轴系扭振危险截面进行扭振疲劳寿命损耗计算；

步骤4：建立不同故障下机组电磁力矩最大值与轴系扭振危险截面扭振疲劳寿命损耗之间的对应关系；

步骤5：根据步骤4所建立的对应关系，确定机组轴系扭振损伤报警阈值和扭振跳机保护阈值；

步骤6：通过监测机组发电机三相电流和电压，实时在线计算相应电磁力矩；

步骤7：当发电机电磁力矩出现大幅跳变，电磁力矩值超过扭振损伤报警阈值时，对机组轴系发生的瞬态冲击类扭振进行报警；当电磁力矩值超过跳机保护阈值时，则判定机组轴系发生损伤程度较大的瞬态冲击类扭振，发出跳机保护信号。

所述大扰动故障包括：电力系统的两相短路故障、三相短路故障、120°非同期故障、180°非同期故障。

所述步骤2包括：首先建立机组轴系多段集中质量扭振模型；计算汽轮机蒸汽力矩，并根据汽轮机各级叶轮的出力分配情况，将蒸汽力矩按比例施加在扭振模型中汽轮机各级叶轮所在轮盘上；通过发电机三相电流电压计算机组在扭振故障下的电磁力矩响应，并将电磁力矩均匀分布在扭振模型中发电机绕组所在的轮盘上；

然后根据机组轴系所受蒸汽力矩和电磁力矩，采用Newmark-β法和传递矩阵法相结合，实时在线计算机组轴系扭振危险截面扭矩响应，然后利用轴系扭振危险截面的扭矩响应计算其名义扭应力响应，考虑结构对应力集中的影响，用理论应力集中系数对名义扭应力作修正，最终得到轴系扭振危险截面的局部扭应力响应。

所述计算汽轮机蒸汽力矩T_m的公式为：

$T_m=\frac{1}{{\mathrm{\ω}}_e}[(u_{\mathrm{ao}}{i}_{a0}+u_{b0}{i}_{b0}+u_{e0}{i}_{e0})+r(i_{a0}^2+i_{b0}^2+i_{c0}^2)]---\left(12\right)$

式中，r为电枢电阻，ωe为稳态发电机角速度，u_a0、u_b0、u_e0为扰动前发电机端三相电压测量值，i_a0、i_b0、i_c0为扰动前发电机端电流测量值。

所述计算机组轴系扭振危险截面扭矩响应的方法为：

首先建立轴系多段集中质量扭振模型，模型由N+1个单元组成，每个单元分别包括一个无弹性的刚性圆盘和一个无质量的弹性轴段；设I_i为第i个单元刚性圆盘的转动惯量，为第i个单元刚性圆盘的角加速度，C_i为第i个单元刚性圆盘的阻尼系数，为第i个单元刚性圆盘的角速度，T_Li为第i个单元刚性圆盘所受外力矩，T_i ^L和T_i ^R分别为第i个单元刚性圆盘左右两端扭矩，其中i＝1、2...N；

以第i个单元刚性圆盘为研究对象，由扭矩平衡条件得：

$I_i{\stackrel{..}{\mathrm{\φ}}}_i=-C_i{\stackrel{.}{\mathrm{\φ}}}_i+T_{\mathrm{Li}}+{T_i}^R-{T_i}^L---\left(13\right)$

将式(13)移项并改写为增量形式为：

$\mathrm{\Δ}{T_i}^R=\mathrm{\Δ}{T_i}^L+I_i\mathrm{\Δ}\stackrel{..}{\mathrm{\φ}}+C_i\mathrm{\Δ}{\stackrel{.}{\mathrm{\φ}}}_i-\mathrm{\Δ}{T}_{\mathrm{Li}}---\left(14\right)$

其中，ΔI_i为第i个单元刚性圆盘的转动惯量在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘的角加速度在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘的角速度在单位时间变化量Δt时间段内的增量，ΔT_Li为第i个单元刚性圆盘所受外力矩在单位时间变化量Δt时间段内的增量，ΔT_i ^L和ΔT_i ^R分别为第i个单元刚性圆盘左右两端扭矩在单位时间变化量Δt时间段内的增量；

式(14)根据Newmark-β法推导得：

ΔT_i ^R＝ΔT_i ^L+A_iΔφ+B_i

$A_i=\frac{I_i}{{\mathrm{\β\Δt}}^2}+\frac{C_i\mathrm{\γ}}{\mathrm{\β\Δt}},B_i=-I_i(\frac{1}{\mathrm{\β\Δt}}{\stackrel{.}{\mathrm{\φ}}}_i+\frac{1}{2\mathrm{\β}}{\stackrel{..}{\mathrm{\φ}}}_i)-C_i[\frac{\mathrm{\γ}}{\mathrm{\β}}{\stackrel{.}{\mathrm{\φ}}}_i+(\frac{\mathrm{\γ}}{2\mathrm{\β}}-1){\stackrel{..}{\mathrm{\φ}}}_i\mathrm{\Δt}]-\mathrm{\Δ}{T}_{\mathrm{Li}}---\left(15\right)$

其中A_i与B_i为中间变量，Δt为单位时间变化量，Δφ为第i个单元刚性圆盘的角位移在单位时间变化量Δt时间段内的增量，参数γ取0.25，参数β取0.5，式(15)表示了第i个单元刚性圆盘左右两端的扭矩变化量关系；

然后设为第i个单元刚性圆盘右端角位移在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘左端角位移在单位时间变化量Δt时间段内的增量，第i个单元刚性圆盘左右两端的角位移变化量关系为：

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

以轴系扭振模型第i个单元弹性轴段为研究对象，推导其左右两端扭矩变化量和角位移变化量关系得：

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

其中为第i+1个单元刚性圆盘左端扭矩在单位时间变化量Δt时间段内的增量，为第i+1个单元刚性圆盘左端角位移在单位时间变化量Δt时间段内的增量，K_i为第i个单元弹性轴段的扭转刚度；

综合上述过程得轴系扭振模型第i+1个单元左右两端扭矩变化量和角位移变化量与第i个单元左右两端扭矩变化量和角位移变化量的传递关系为：

${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L=\left[\begin{array}{cc}1& A_i\\ \frac{1}{K_i}& 1+\frac{A_i}{K_i}\end{array}\right]{\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L+\left\{\begin{array}{c}{B}_i\\ \frac{B_i}{K_i}\end{array}\right\}---\left(18\right)$

式(18)为矩阵形式，其中 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L$ 为第i+1个单元刚性圆盘左端扭矩变化量和角位移变化量写成矩阵表达形式， ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L$ 为第i个单元刚性圆盘左端扭矩变化量和角位移变化量写成的矩阵表达形式；

为方便推导计算，设中间变量f＝ΔT，e＝Δφ，U₁₁＝1,U₁₂＝A_i， F_f＝B_i，则式(18)写为如下矩阵表达形式：

${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}=\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+\left\{\begin{array}{c}{F}_f\\ F_e\end{array}\right\}---\left(19\right)$

其中 ${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}$ 为 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L,{\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i$ 为 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L;$

根据Riccati传递矩阵法，设Riccati变换为：

f_i＝S_ie_i+P_i        (20)

其中f_i、S_i、e_i、P_i为第i个单元的中间变量，f_i+1、S_i+1、e_i+1、P_i+1为第i+1个单元的中间变量(i＝1,2,3......N)，N代表轴系扭振模型单元数量，第N+1个单元为假设单元；

将式(20)带入式(19)，得到S_i+1、e_i和P_i+1的递推关系表达式如下：

$S_{i+1}=\frac{U_{11}{S}_i+U_{12}}{U_{21}{S}_i+U_{22}}---\left(21\right)$

P_i+1＝(U₁₁P_i+F_f)-S_i+1(U₂₁P_i+F_e)       (22)

$e_i=\frac{e_{i+1}}{U_{21}{S}_i+U_{22}}-\frac{U_{21}{P}_i+F_e}{U_{21}{S}_i+U_{22}}---\left(23\right)$

当轴系扭振模型由N个单元组成时，由于模型两端为自由端，为方便计算设模型由N+1个单元组成，且第N+1个单元为不存在单元；由模型两端为自由端得f₁＝0，e₁≠0，S₁＝0，P₁＝0，根据式(21)和式(22)计算得S₂，P₂，S₃，P₃…，S_N+1，P_N+1；对于模型第N+1个单元，由模型两端自由得f_N+1＝0，根据式(20)得根据式(23)依次计算得到e_i(i＝N，N-1，N-2，......，1)，相应算出各截面f_i，所求结果即为各单元状态参数在单位时间变化量Δt时间段内的增量。

所述轴系扭振危险截面的局部扭应力响应的计算方法为：

设机组轴系危险截面扭矩为T，轴系危险截面抗扭截面系数为W，危险截面应力集中系数为K，则机组轴系危险截面局部扭应力为求解得出轴系危险截面扭矩响应后，经计算得轴系危险截面扭应力响应，同时用应力集中系数对扭应力进行修正，得到危险截面局部扭应力响应。

所述步骤3根据轴系扭振危险截面材料扭转S-N曲线和局部扭应力响应，采用雨流法和线性累积法对轴系扭振危险截面进行扭振疲劳寿命损耗计算。

所述步骤4中的对应关系是指不同的电磁力矩最大标幺值条件下分别对应不同大扰动故障所引起的危险截面A和危险截面B分别所对应的轴系的扭振疲劳寿命损耗的百分比值。

所述步骤5中确定机组轴系扭振损伤报警阈值的方法为：随着电磁力矩最大标幺值增大，查找首次出现不同大扰动故障所引起的危险截面A或危险截面B分别所对应的轴系的扭振疲劳寿命损耗的百分比值为非零的情况，该情况所对应的电磁力矩最大标幺值为机组轴系扭振损伤报警阈值。

确定机组轴系扭振跳机保护阈值的方法为：随着电磁力矩最大标幺值增大，查找首次出现不同大扰动故障所引起的的危险截面A或危险截面B分别所对应的轴系的扭振疲劳寿命损耗达到1％的情况，该情况所对应的电磁力矩最大标幺值为机组轴系扭振损伤报警阈值。

一种用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的系统，包括：仿真模块、扭应力计算模块、疲劳寿命损耗计算模块、数据存储模块、阈值设定模块、电磁力矩监控模块、报警模块；

仿真模块依次通过扭应力计算模块、疲劳寿命损耗计算模块、数据存储模块、阈值设定模块、电磁力矩监控模块和报警模块相连；

其中，仿真模块的功能是：通过电力系统仿真，模拟电力系统大扰动故障下机组发生瞬态冲击类扭振时发电机的电磁力矩，仿真得到不同类型瞬态冲击类扭振故障下不同幅值的电磁力矩响应；

扭应力计算模块的功能是：根据仿真所得电磁力矩响应计算轴系扭振危险截面局部扭应力响应；

疲劳寿命损耗计算模块的功能是：对机组轴系扭振危险截面进行扭振疲劳寿命损耗计算；

数据存储模块的功能是：建立不同故障下机组电磁力矩最大值与轴系扭振危险截面扭振疲劳寿命损耗之间的对应关系；

阈值设定模块的功能是：确定机组轴系扭振损伤报警阈值和扭振跳机保护阈值；

电磁力矩监控模块的功能是：通过监测机组发电机三相电流和电压，实时在线计算相应电磁力矩；

报警模块的功能是：当发电机电磁力矩出现大幅跳变，电磁力矩值超过扭振损伤报警阈值时，对机组轴系发生的瞬态冲击类扭振进行报警；当电磁力矩值超过跳机保护阈值时，则判定机组轴系发生损伤程度较大的瞬态冲击类扭振，发出跳机保护信号。

本发明的有益效果在于：首先对机电网协调作用下电力系统大扰动导致的轴系扭振进行仿真计算，建立电力系统大扰动下发电机电磁力矩最大值与轴系疲劳寿命损耗之间的对应关系，制定损伤报警阈值和跳机保护阈值，然后根据发电机三相电流和电压计算发电机电磁力矩，在机组安全性受到威胁时及时地发出损伤报警或跳机保护信号，最后再对扭振故障进行轴系扭振动态响应计算，进而对轴系做出精确的损伤评估和剩余寿命评估；该方法和系统可实现针对电力系统大扰动故障的扭振在线监测、分析与保护，同时可对机组轴系扭振疲劳寿命损耗进行精确计算。

附图说明

图1为用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法流程图；

图2为汽轮发电机组轴系多段集中质量扭振模型示意图；

图3为用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的系统结构图。

具体实施方式

下面结合附图，对优选实施例作详细说明。

一种用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法，如图1所示，包括：

步骤1：通过电力系统仿真，模拟电力系统不同类型的大扰动故障下机组发生瞬态冲击类扭振时发电机的电磁力矩，仿真得到不同类型瞬态冲击类扭振故障下不同幅值的电磁力矩响应；

步骤2：根据仿真所得电磁力矩响应计算轴系扭振危险截面局部扭应力响应；

步骤3：对机组轴系扭振危险截面进行扭振疲劳寿命损耗计算；

步骤4：建立不同故障下机组电磁力矩最大值与轴系扭振危险截面扭振疲劳寿命损耗之间的对应关系；

步骤5：根据步骤4所建立的对应关系，确定机组轴系扭振损伤报警阈值和扭振跳机保护阈值；

步骤6：通过监测机组发电机三相电流和电压，实时在线计算相应电磁力矩；

步骤7：当发电机电磁力矩出现大幅跳变，电磁力矩值超过扭振损伤报警阈值时，对机组轴系发生的瞬态冲击类扭振进行报警；当电磁力矩值超过跳机保护阈值时，则判定机组轴系发生损伤程度较大的瞬态冲击类扭振，发出跳机保护信号。

步骤1：通过电力系统仿真，模拟电力系统大扰动故障下机组发生瞬态冲击类扭振时发电机的电磁力矩，仿真得到不同类型瞬态冲击类扭振故障下不同幅值的电磁力矩响应。

瞬态冲击类扭振电磁力矩仿真方法

由于瞬态冲击类扭振下的电磁力矩需要考虑系统的强非线性，同步发电机模型用派克方程描述的5阶微分方程表示，采用X_ad基值系统，则同步电机电磁力矩瞬时值T_e可表示为：

T_e＝ψ_di_q-ψ_qi_d＝(-x_di_d+x_qfI_f+x_aDI_D)i_q-(-x_qi_q+x_aQI_Q)i_d     (1)

式中，x_d为定子绕组d轴自感系数；x_q为定子绕组q轴自感系数；x_ad为d轴互感系数；x_aq为q轴互感系数；I_f励磁电流i_d、i_q、I_D、I_Q分别为定子绕组和阻尼绕组纵轴和横轴的电流，且x_af＝x_aD＝x_ad，x_aQ＝x_aq对于汽轮发电机有x_d＝x_q。

大型同步发电机绕组均是星形接法，零轴电流i₀＝0，则根据派克变换有dq0基本方程：

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

离散化后有：

$i_d^t=\frac{2}{3}\cos {\mathrm{\γ}}^t\·I_a^t+\frac{2}{3}\cos ({\mathrm{\γ}}^t-\frac{2}{3}\mathrm{\π})\·I_b^t+\frac{2}{3}\cos ({\mathrm{\γ}}^t+\frac{2}{3}\mathrm{\π})\·I_c^t---\left(3\right)$

$i_q^t=-[\frac{2}{3}\sin {\mathrm{\γ}}^t\·I_a^t+\frac{2}{3}\sin ({\mathrm{\γ}}^t-\frac{2}{3}\mathrm{\π})\·I_b^t+\frac{2}{3}\sin ({\mathrm{\γ}}^t+\frac{2}{3}\mathrm{\π})\·I_c^t]---\left(4\right)$

式中，t＝nΔt，Δt为计算步长，为三相电流采样值；γ为电机转子d轴与定子a相绕组轴线的夹角，有：

${\mathrm{\γ}}^t={\∫}_0^t\mathrm{\ωdt}+{\mathrm{\γ}}_0---\left(5\right)$

ω为每一时步的角速度，ω＝ω_s+ω_n，即同步角速度与扭转角速度之和。

对于P极电机，如果已知由响应计算求得的发电机轴段扭角δ，对应的电角度会变化Pδ。式(5)变成下面形式：

${\mathrm{\γ}}^t={\∫}_0^t{\mathrm{\ω}}_s\mathrm{dt}+{\mathrm{\γ}}_0+\mathrm{P\δ}---\left(6\right)$

因为γ₀为初始角，根据下式可求出初始角γ₀为：

式中，I_m为相电流幅值，U_m为相电压幅值，R_a为定子绕组内电阻，为功率因数角，可由实测的a相电压和电流采样值求出。

到此可计算出i_d、i_q，下面讨论一下如何求I_1d、I_1q。

稳态时，阻尼绕组电流很小，可忽略，即认为I_1d、I_1q初始值为零。由于阻尼绕组D、和阻尼绕组Q的端电压恒为零，进入暂态时，有下式：

$U_{1d}=-x_{\mathrm{ad}}\frac{{\mathrm{di}}_d}{\mathrm{dt}}+x_{1\mathrm{fd}}\frac{{\mathrm{dI}}_f}{\mathrm{dt}}+x_{11d}\frac{{\mathrm{dI}}_D}{\mathrm{dt}}+r_D{I}_D=0---\left(8\right)$

$U_{1q}=-x_{\mathrm{aq}}\frac{{\mathrm{di}}_d}{\mathrm{dt}}+x_Q\frac{{\mathrm{dI}}_Q}{\mathrm{dt}}+r_Q{I}_Q=0---\left(9\right)$

式中，x_ad＝x_d-x₁，x_1fd＝x_ad， $x_{11d}=x_{D1}+x_{\mathrm{ad}}=\frac{{(x_d^{\′}-x_1)}^2}{x_d^{\′}-x_d^{\′\′}}-x_d^{\′}+x_1+x_{\mathrm{ad}},$

$r_D=\frac{(x_{D1}+1/(x_d^{\′}-x_1))}{{\mathrm{\ω}}_B{T}_{d0}^{\′\′}},$ x_aq＝x_q-x₁， $x_Q=\frac{x_{\mathrm{ad}}^2}{x_q-x_q^{\′\′}},r_Q=\frac{x_Q}{{\mathrm{\ω}}_B{T}_{q0}^{\′\′}}$ 为同步电机基本参数。

若设 $\frac{{\mathrm{dI}}_D}{\mathrm{dt}}=\frac{I_D^{t+\mathrm{\Δt}}-I_D^t}{\mathrm{\Δt}},I_D=\frac{I_D^{t+\mathrm{\Δt}}+I_D^t}{2}$ 可得：

$({2x}_{11d}+\mathrm{\Δ}{\mathrm{tr}}_D)I_D^{t+\mathrm{\Δt}}=2x_{\mathrm{ad}}(i_d^{t+\mathrm{\Δt}}-i_d^t)-{2x}_{1\mathrm{fd}}(I_f^{t+\mathrm{\Δt}}-I_f^t)+({2x}_{11d}-\mathrm{\Δ}{\mathrm{tr}}_D)I_D^t---\left(10\right)$

$({2x}_Q+\mathrm{\Δ}{\mathrm{tr}}_Q)I_Q^{t+\mathrm{\Δt}}=2x_{\mathrm{aq}}(i_d^{t+\mathrm{\Δt}}-i_d^t)+({2x}_Q-\mathrm{\Δ}{\mathrm{tr}}_Q)I_Q^t---\left(11\right)$

通过上面两个式子可知，只需把求出的i_d、i_q代入，并令I_D、I_Q初始值为零，即可求解I_D、I_Q。把以上的计算结果代入式(1)可求出电磁力矩瞬时值T_e。

步骤2、计算机组轴系扭振危险截面局部扭应力响应。

首先，建立机组轴系多段集中质量扭振模型；计算汽轮机蒸汽力矩，并根据汽轮机各级叶轮的出力分配情况，将蒸汽力矩按比例施加在扭振模型中汽轮机各级叶轮所在轮盘上；通过发电机三相电流电压计算机组在扭振故障下的电磁力矩响应，并将电磁力矩均匀分布在扭振模型中发电机绕组所在的轮盘上。

然后，根据机组轴系所受蒸汽力矩和电磁力矩，采用Newmark-β法和传递矩阵法相结合，实时在线计算机组轴系扭振危险截面扭矩响应，然后利用轴系扭振危险截面的扭矩响应计算其名义扭应力响应，考虑结构对应力集中的影响，用理论应力集中系数对名义扭应力做修正，最终得到轴系扭振危险截面的局部扭应力响应。

具体过程包括：

1)机组蒸汽力矩计算方法

汽轮发电机组处于稳态运行时，汽轮机机械的蒸汽扭矩与电磁扭矩平衡，所以可由电磁扭矩求出蒸汽力矩T_m为：

$T_m=\frac{1}{{\mathrm{\ω}}_e}[(u_{\mathrm{ao}}{i}_{a0}+u_{b0}{i}_{b0}+u_{e0}{i}_{e0})+r(i_{a0}^2+i_{b0}^2+i_{c0}^2)]---\left(12\right)$

式中，r为电枢电阻，ωe为稳态发电机角速度，u_a0，u_b0，u_e0，i_a0，i_b0，i_c0为扰动前发电机端三相电压、电流测量值。

2)机组轴系扭振危险截面扭矩响应计算方法

首先建立轴系多段集中质量扭振模型，如图2所示，模型由N+1个单元组成，每个单元分别包括一个无弹性的刚性圆盘和一个无质量的弹性轴段；设I_i为第i个单元刚性圆盘的转动惯量，为第i个单元刚性圆盘的角加速度，C_i为第i个单元刚性圆盘的阻尼系数，为第i个单元刚性圆盘的角速度，T_Li为第i个单元刚性圆盘所受外力矩，T_i ^L和T_i ^R分别为第i个单元刚性圆盘左右两端扭矩，其中i＝1、2...N；以第i个单元刚性圆盘为研究对象，由扭矩平衡条件得：

$I_i{\stackrel{..}{\mathrm{\φ}}}_i=-C_i{\stackrel{.}{\mathrm{\φ}}}_i+T_{\mathrm{Li}}+{T_i}^R-{T_i}^L---\left(13\right)$

将式(13)移项并改写为增量形式为：

$\mathrm{\Δ}{T_i}^R=\mathrm{\Δ}{T_i}^L+I_i\mathrm{\Δ}\stackrel{..}{\mathrm{\φ}}+C_i\mathrm{\Δ}{\stackrel{.}{\mathrm{\φ}}}_i-\mathrm{\Δ}{T}_{\mathrm{Li}}---\left(14\right)$

其中，ΔI_i为第i个单元刚性圆盘的转动惯量在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘的角加速度在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘的角速度在单位时间变化量Δt时间段内的增量，ΔT_Li为第i个单元刚性圆盘所受外力矩在单位时间变化量Δt时间段内的增量，ΔT_i ^L和ΔT_i ^R分别为第i个单元刚性圆盘左右两端扭矩在单位时间变化量Δt时间段内的增量；

式(14)根据Newmark-β法推导得：

ΔT_i ^R＝ΔT_i ^L+A_iΔφ+B_i

$A_i=\frac{I_i}{{\mathrm{\β\Δt}}^2}+\frac{C_i\mathrm{\γ}}{\mathrm{\β\Δt}},B_i=-I_i(\frac{1}{\mathrm{\β\Δt}}{\stackrel{.}{\mathrm{\φ}}}_i+\frac{1}{2\mathrm{\β}}{\stackrel{..}{\mathrm{\φ}}}_i)-C_i[\frac{\mathrm{\γ}}{\mathrm{\β}}{\stackrel{.}{\mathrm{\φ}}}_i+(\frac{\mathrm{\γ}}{2\mathrm{\β}}-1){\stackrel{..}{\mathrm{\φ}}}_i\mathrm{\Δt}]-\mathrm{\Δ}{T}_{\mathrm{Li}}---\left(15\right)$

其中A_i与B_i为中间变量，Δt为单位时间变化量，Δφ为第i个单元刚性圆盘的角位移在单位时间变化量Δt时间段内的增量，参数γ取0.25，参数β取0.5，式(15)表示了第i个单元刚性圆盘左右两端的扭矩变化量关系；

然后设为第i个单元刚性圆盘右端角位移在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘左端角位移在单位时间变化量Δt时间段内的增量，第i个单元刚性圆盘左右两端的角位移变化量关系为：

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

以轴系扭振模型第i个单元弹性轴段为研究对象，推导其左右两端扭矩变化量和角位移变化量关系得：

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

其中为第i+1个单元刚性圆盘左端扭矩在单位时间变化量Δt时间段内的增量，为第i+1个单元刚性圆盘左端角位移在单位时间变化量Δt时间段内的增量，K_i为第i个单元弹性轴段的扭转刚度；

综合上述过程得轴系扭振模型第i+1个单元左右两端扭矩变化量和角位移变化量与第i个单元左右两端扭矩变化量和角位移变化量的传递关系为：

${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L=\left[\begin{array}{cc}1& A_i\\ \frac{1}{K_i}& 1+\frac{A_i}{K_i}\end{array}\right]{\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L+\left\{\begin{array}{c}{B}_i\\ \frac{B_i}{K_i}\end{array}\right\}---\left(18\right)$

式(18)为矩阵形式，其中 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L$ 为第i+1个单元刚性圆盘左端扭矩变化量和角位移变化量写成矩阵表达形式， ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L$ 为第i个单元刚性圆盘左端扭矩变化量和角位移变化量写成的矩阵表达形式；

为方便推导计算，设中间变量f＝ΔT，e＝Δφ，U₁₁＝1,U₁₂＝A_i， F_f＝B_i，则式(18)写为如下矩阵表达形式：

${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}=\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+\left\{\begin{array}{c}{F}_f\\ F_e\end{array}\right\}---\left(19\right)$

其中 ${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}$ 为 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L,{\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i$ 为 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L;$

根据Riccati传递矩阵法，设Riccati变换为：

f_i＝S_ie_i+P_i      (20)

其中f_i、S_i、e_i、P_i为第i个单元的中间变量，f_i+1、S_i+1、e_i+1、P_i+1为第i+1个单元的中间变量(i＝1,2,3......N)，N代表轴系扭振模型单元数量，第N+1个单元为假设单元；

将式(20)带入式(19)，得到S_i+1、e_i和P_i+1的递推关系表达式如下：

$S_{i+1}=\frac{U_{11}{S}_i+U_{12}}{U_{21}{S}_i+U_{22}}---\left(21\right)$

P_i+1＝(U₁₁P_i+F_f)-S_i+1(U₂₁P_i+F_e)      (22)

$e_i=\frac{e_{i+1}}{U_{21}{S}_i+U_{22}}-\frac{U_{21}{P}_i+F_e}{U_{21}{S}_i+U_{22}}---\left(23\right)$

当轴系扭振模型由N个单元组成时，由于模型两端为自由端，为方便计算设模型由N+1个单元组成，且第N+1个单元为不存在单元；由模型两端为自由端得f₁＝0，e₁≠0，S₁＝0，P₁＝0，根据式(21)和式(22)计算得S₂，P₂，S₃，P₃…，S_N+1，P_N+1；对于模型第N+1个单元，由模型两端自由得f_N+1＝0，根据式(20)得根据式(23)依次计算得到e_i(i＝N，N-1，N-2，......，1)，相应算出各截面f_i，所求结果即为各单元状态参数在单位时间变化量Δt时间段内的增量。

综上所述，轴系扭振危险截面扭矩响应计算流程如下：

a计算模型各单元参数I_i，C_i，K_i，ΔT_Li和Δt(i＝1,2,3......N)；给定轴系扭振模型各单元刚性圆盘初始时刻t₀时的φ_i，和(i＝1,2,3......N)。

b计算模型中间变量A_i，B_i，U₁₁，U₁₂，U₂₁，U₂₂，F_f，F_e。

c由模型两端为自由端可得f₁＝0，e₁≠0(为任意值)，可得S₁＝0，P₁＝0，根据式(21)和式(22)计算可得S₂，P₂，S₃，P₃……，S_N+1，P_N+1。

d由模型两端自由可得f_N+1＝0，根据式(20)可得根据式(23)可依次计算得到e_i(i＝N，N-1，N-2，......，1)，相应也可算出各截面f_i。

e由时刻t₀时计算得到的e_i和f_i计算时刻t₀+Δt时轴系扭振模型各单元刚性圆盘的φ_i和扭矩T_i(i＝1,2,3......N)。

f根据Newmark-β法计算时刻t₀+Δt时轴系扭振模型各单元刚性圆盘的和((i＝1,2,3......N)。

g以初始时刻为t₀+Δt，重复步骤b计算时刻t₀+2Δt时轴系扭振模型各单元刚性圆盘的T_i，φ_i，和((i＝1,2,3......N)。如此反复，直到计算得到要求时间段轴系扭振危险截面扭矩响应。

3)机组轴系危险截面扭应力计算方法

设机组轴系危险截面扭矩为T，轴系危险截面抗扭截面系数为W，危险截面应力集中系数为K，则机组轴系危险截面局部扭应力为

如上所述求解得出轴系危险截面扭矩响应后，经计算可得轴系危险截面扭应力响应，同时用应力集中系数对扭应力进行修正，得到危险截面局部扭应力响应。

步骤3：对机组轴系扭振危险截面进行扭振疲劳寿命损耗计算；

根据轴系扭振危险截面材料扭转S-N曲线和局部扭应力响应，采用雨流法和线性累积法对轴系扭振危险截面进行扭振疲劳寿命损耗计算。

1)雨流计数法

把应变-时间历程数据记录转过90°，时间坐标轴竖直向下，雨流计数法对载荷的时间历程进行计数。雨流计数法的基本计数规则为：

(1)雨流依次从载荷时间历程的峰值位置的内侧沿着斜坡往下流；

(2)雨流从某一个峰值点开始流动，当遇到比其起始峰值更大的峰值时要停止流动；

(3)雨流遇到上面流下的雨流时，必须停止流动；

(4)取出所有的全循环，记下每个循环的幅度；

(5)将第一阶段计数后剩下的发散收敛载荷时间历程等效为一个收敛发散型的载荷时间历程，进行第二阶段的雨流计数。计数循环的总数等于两个计数阶段的计数循环之和。

2)线性累积法

线性疲劳累积操作理论是指在循环载荷作用下，疲劳损伤是可以线性地累积的，各个应力之间相互独立和互不相关，当累积的损伤达到某一数值时，试件或构件就发生疲劳破坏。

表1  线性疲劳累积理论

Miner理论指出：

(1)一个循环造成的损伤：

$D=\frac{1}{N}$

式中，N为对应于当前载荷水平S的疲劳水平；

(2)等幅载荷下，n个循环造成的损伤：

$D=\frac{n}{N}$

变幅载荷下，n个循环造成的损伤：

$D=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{N_i}$

式中，Ni为对应于当前载荷水平Si的疲劳寿命。

步骤4：根据仿真所得电磁力矩响应计算轴系扭振危险截面局部扭应力响应，并根据轴系扭振危险截面材料扭转S-N曲线进行扭振疲劳寿命损耗评估，建立不同故障下机组电磁力矩最大值与轴系扭振危险截面扭振疲劳寿命损耗之间的对应关系；建立不同的电磁力矩最大标幺值条件下分别对应不同大扰动故障所引起的危险截面A和危险截面B分别所对应的轴系的扭振疲劳寿命损耗的百分比值的对应关系。例如某600MW机组在电磁力矩瞬态冲击下轴系的扭振疲劳寿命损耗如表2所示。

表2  某600MW机组在电磁力矩瞬态冲击下轴系的疲劳寿命损耗(％)

步骤5：根据步骤4所建立的对应关系，确定机组轴系扭振损伤报警阈值和扭振跳机保护阈值。具体方法为：随着电磁力矩最大标幺值增大，查找首次出现不同大扰动故障所引起的危险截面A或危险截面B分别所对应的轴系的扭振疲劳寿命损耗的百分比值为非零的情况，该情况所对应的电磁力矩最大标幺值为机组轴系扭振损伤报警阈值；

确定机组轴系扭振跳机保护阈值的方法为：随着电磁力矩最大标幺值增大，查找首次出现不同大扰动故障所引起的的危险截面A或危险截面B分别所对应的轴系的扭振疲劳寿命损耗达到1％的情况，该情况所对应的电磁力矩最大标幺值为机组轴系扭振损伤报警阈值。

由表2中数据可知，瞬态冲击类扭振故障下，机组轴系扭振疲劳寿命损耗主要取决于机组发电机电磁力矩的最大值，电磁力矩最大值相同的不同类型瞬态冲击类扭振故障所造成的轴系扭振疲劳寿命损耗较为接近。由上表可知，当该机组电磁力矩最大值为3倍标幺值时，轴系可能承受超出其疲劳极限的扭转应力。因此，当该机组出现超过3倍标幺值的电磁力矩时，该机组轴系有出现扭振疲劳寿命损耗的风险，故而把电磁力矩的3倍标幺值作为该机组的扭振损伤报警阈值。由上表还可得，当该机组电磁力矩最大值为8～9倍标幺值时，轴系扭振疲劳寿命损耗将达到1％左右，经计算机组电磁力矩最大值为8.5倍标幺值时，轴系扭振寿命损耗达到1％，故而把电磁力矩的8.5倍标幺值作为该机组的扭振跳机保护阈值。

步骤6：通过监测机组发电机三相电流和电压，实时在线计算相应电磁力矩。

步骤7：当发电机电磁力矩出现大幅跳变，电磁力矩值超过损伤报警阈值时，则对机组轴系发生的瞬态冲击类扭振进行报警；当电磁力矩值超过跳机保护阈值时，则判定机组轴系发生损伤程度较大的瞬态冲击类扭振，发出跳机保护信号。

一种用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的系统，如图3所示，包括：仿真模块、扭应力计算模块、疲劳寿命损耗计算模块、数据存储模块、阈值设定模块、电磁力矩监控模块、报警模块；

仿真模块依次通过扭应力计算模块、疲劳寿命损耗计算模块、数据存储模块、阈值设定模块、电磁力矩监控模块和报警模块相连；

其中，仿真模块的功能是：通过电力系统仿真，模拟电力系统大扰动故障下机组发生瞬态冲击类扭振时发电机的电磁力矩，仿真得到不同类型瞬态冲击类扭振故障下不同幅值的电磁力矩响应；

扭应力计算模块的功能是：根据仿真所得电磁力矩响应计算轴系扭振危险截面局部扭应力响应；

疲劳寿命损耗计算模块的功能是：对机组轴系扭振危险截面进行扭振疲劳寿命损耗计算；

数据存储模块的功能是：建立不同故障下机组电磁力矩最大值与轴系扭振危险截面扭振疲劳寿命损耗之间的对应关系；

阈值设定模块的功能是：确定机组轴系扭振损伤报警阈值和扭振跳机保护阈值；

电磁力矩监控模块的功能是：通过监测机组发电机三相电流和电压，实时在线计算相应电磁力矩；

报警模块的功能是：当发电机电磁力矩出现大幅跳变，电磁力矩值超过扭振损伤报警阈值时，对机组轴系发生的瞬态冲击类扭振进行报警；当电磁力矩值超过跳机保护阈值时，则判定机组轴系发生损伤程度较大的瞬态冲击类扭振，发出跳机保护信号。

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

Claims (10)

1.一种用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的方法，其特征在于，包括：

步骤1：通过电力系统仿真，模拟电力系统不同类型的大扰动故障下机组发生瞬态冲击类扭振时发电机的电磁力矩，仿真得到不同类型瞬态冲击类扭振故障下不同幅值的电磁力矩响应；

步骤2：根据仿真所得电磁力矩响应计算轴系扭振危险截面局部扭应力响应；

步骤3：对机组轴系扭振危险截面进行扭振疲劳寿命损耗计算；

步骤4：建立不同故障下机组电磁力矩最大值与轴系扭振危险截面扭振疲劳寿命损耗之间的对应关系；

步骤5：根据步骤4所建立的对应关系，确定机组轴系扭振损伤报警阈值和扭振跳机保护阈值；

步骤6：通过监测机组发电机三相电流和电压，实时在线计算相应电磁力矩；

步骤7：当发电机电磁力矩出现大幅跳变，电磁力矩值超过扭振损伤报警阈值时，对机组轴系发生的瞬态冲击类扭振进行报警；当电磁力矩值超过跳机保护阈值时，则判定机组轴系发生损伤程度较大的瞬态冲击类扭振，发出跳机保护信号。

2.根据权利要求1所述的方法，其特征在于，所述大扰动故障包括：电力系统的两相短路故障、三相短路故障、120°非同期故障、180°非同期故障。

3.根据权利要求1所述的方法，其特征在于，所述步骤2包括：首先建立机组轴系多段集中质量扭振模型；计算汽轮机蒸汽力矩，并根据汽轮机各级叶轮的出力分配情况，将蒸汽力矩按比例施加在扭振模型中汽轮机各级叶轮所在轮盘上；通过发电机三相电流电压计算机组在扭振故障下的电磁力矩响应，并将电磁力矩均匀分布在扭振模型中发电机绕组所在的轮盘上；

然后根据机组轴系所受蒸汽力矩和电磁力矩，采用Newmark-β法和传递矩阵法相结合，实时在线计算机组轴系扭振危险截面扭矩响应，然后利用轴系扭振危险截面的扭矩响应计算其名义扭应力响应，考虑结构对应力集中的影响，用理论应力集中系数对名义扭应力作修正，最终得到轴系扭振危险截面的局部扭应力响应。

4.根据权利要求1或3所述的方法，其特征在于，所述计算汽轮机蒸汽力矩T_m的公式为：

$T_m=\frac{1}{{\mathrm{\ω}}_e}[(u_{\mathrm{ao}}{i}_{a0}+u_{b0}{i}_{b0}+u_{e0}{i}_{e0})+r(i_{a0}^2+i_{b0}^2+i_{c0}^2)]---\left(12\right)$

式中，r为电枢电阻，ωe为稳态发电机角速度，u_a0、u_b0、u_e0为扰动前发电机端三相电压测量值，i_a0、i_b0、i_c0为扰动前发电机端电流测量值。

5.根据权利要求1或3所述的方法，其特征在于，所述计算机组轴系扭振危险截面扭矩响应的方法为：

首先建立轴系多段集中质量扭振模型，模型由N+1个单元组成，每个单元分别包括一个无弹性的刚性圆盘和一个无质量的弹性轴段；设I_i为第i个单元刚性圆盘的转动惯量，为第i个单元刚性圆盘的角加速度，C_i为第i个单元刚性圆盘的阻尼系数，为第i个单元刚性圆盘的角速度，T_Li为第i个单元刚性圆盘所受外力矩，T_i ^L和T_i ^R分别为第i个单元刚性圆盘左右两端扭矩，其中i＝1、2...N；

以第i个单元刚性圆盘为研究对象，由扭矩平衡条件得：

$I_i{\stackrel{..}{\mathrm{\φ}}}_i=-C_i{\stackrel{.}{\mathrm{\φ}}}_i+T_{\mathrm{Li}}+{T_i}^R-{T_i}^L---\left(13\right)$

将式(13)移项并改写为增量形式为：

$\mathrm{\Δ}{T_i}^R=\mathrm{\Δ}{T_i}^L+I_i\mathrm{\Δ}\stackrel{..}{\mathrm{\φ}}+C_i\mathrm{\Δ}{\stackrel{.}{\mathrm{\φ}}}_i-\mathrm{\Δ}{T}_{\mathrm{Li}}---\left(14\right)$

其中，ΔI_i为第i个单元刚性圆盘的转动惯量在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘的角加速度在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘的角速度在单位时间变化量Δt时间段内的增量，ΔT_Li为第i个单元刚性圆盘所受外力矩在单位时间变化量Δt时间段内的增量，ΔT_i ^L和ΔT_i ^R分别为第i个单元刚性圆盘左右两端扭矩在单位时间变化量Δt时间段内的增量；

式(14)根据Newmark-β法推导得：

ΔT_i ^R＝ΔT_i ^L+A_iΔφ+B_i

$A_i=\frac{I_i}{{\mathrm{\β\Δt}}^2}+\frac{C_i\mathrm{\γ}}{\mathrm{\β\Δt}},B_i=-I_i(\frac{1}{\mathrm{\β\Δt}}{\stackrel{.}{\mathrm{\φ}}}_i+\frac{1}{2\mathrm{\β}}{\stackrel{..}{\mathrm{\φ}}}_i)-C_i[\frac{\mathrm{\γ}}{\mathrm{\β}}{\stackrel{.}{\mathrm{\φ}}}_i+(\frac{\mathrm{\γ}}{2\mathrm{\β}}-1){\stackrel{..}{\mathrm{\φ}}}_i\mathrm{\Δt}]-\mathrm{\Δ}{T}_{\mathrm{Li}}---\left(15\right)$

其中A_i与B_i为中间变量，Δt为单位时间变化量，Δφ为第i个单元刚性圆盘的角位移在单位时间变化量Δt时间段内的增量，参数γ取0.25，参数β取0.5，式(15)表示了第i个单元刚性圆盘左右两端的扭矩变化量关系；

然后设为第i个单元刚性圆盘右端角位移在单位时间变化量Δt时间段内的增量，为第i个单元刚性圆盘左端角位移在单位时间变化量Δt时间段内的增量，第i个单元刚性圆盘左右两端的角位移变化量关系为：

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

以轴系扭振模型第i个单元弹性轴段为研究对象，推导其左右两端扭矩变化量和角位移变化量关系得：

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

其中为第i+1个单元刚性圆盘左端扭矩在单位时间变化量Δt时间段内的增量，为第i+1个单元刚性圆盘左端角位移在单位时间变化量Δt时间段内的增量，K_i为第i个单元弹性轴段的扭转刚度；

综合上述过程得轴系扭振模型第i+1个单元左右两端扭矩变化量和角位移变化量与第i个单元左右两端扭矩变化量和角位移变化量的传递关系为：

${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L=\left[\begin{array}{cc}1& A_i\\ \frac{1}{K_i}& 1+\frac{A_i}{K_i}\end{array}\right]{\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L+\left\{\begin{array}{c}{B}_i\\ \frac{B_i}{K_i}\end{array}\right\}---\left(18\right)$

式(18)为矩阵形式，其中 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L$ 为第i+1个单元刚性圆盘左端扭矩变化量和角位移变化量写成矩阵表达形式， ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L$ 为第i个单元刚性圆盘左端扭矩变化量和角位移变化量写成的矩阵表达形式；

为方便推导计算，设中间变量f＝ΔT，e＝Δφ，U₁₁＝1,U₁₂＝A_i， F_f＝B_i，则式(18)写为如下矩阵表达形式：

${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}=\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+\left\{\begin{array}{c}{F}_f\\ F_e\end{array}\right\}---\left(19\right)$

其中 ${\left\{\begin{array}{c}f\\ e\end{array}\right\}}_{i+1}$ 为 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_{i+1}^L,{\left\{\begin{array}{c}f\\ e\end{array}\right\}}_i$ 为 ${\left\{\begin{array}{c}\mathrm{\ΔT}\\ \mathrm{\Δ\φ}\end{array}\right\}}_i^L;$

根据Riccati传递矩阵法，设Riccati变换为：

f_i＝S_ie_i+P_i       (20)

其中f_i、S_i、e_i、P_i为第i个单元的中间变量，f_i+1、S_i+1、e_i+1、P_i+1为第i+1个单元的中间变量(i＝1,2,3......N)，N代表轴系扭振模型单元数量，第N+1个单元为假设单元；

将式(20)带入式(19)，得到S_i+1、e_i和P_i+1的递推关系表达式如下：

$S_{i+1}=\frac{U_{11}{S}_i+U_{12}}{U_{21}{S}_i+U_{22}}---\left(21\right)$

P_i+1＝(U₁₁P_i+F_f)-S_i+1(U₂₁P_i+F_e)      (22)

$e_i=\frac{e_{i+1}}{U_{21}{S}_i+U_{22}}-\frac{U_{21}{P}_i+F_e}{U_{21}{S}_i+U_{22}}---\left(23\right)$

当轴系扭振模型由N个单元组成时，由于模型两端为自由端，为方便计算设模型由N+1个单元组成，且第N+1个单元为不存在单元；由模型两端为自由端得f₁＝0，e₁≠0，S₁＝0，P₁＝0，根据式(21)和式(22)计算得S₂，P₂，S₃，P₃…，S_N+1，P_N+1；对于模型第N+1个单元，由模型两端自由得f_N+1＝0，根据式(20)得根据式(23)依次计算得到e_i(i＝N，N-1，N-2，......，1)，相应算出各截面f_i，所求结果即为各单元状态参数在单位时间变化量Δt时间段内的增量。

6.根据权利要求1或3所述的方法，其特征在于，所述轴系扭振危险截面的局部扭应力响应的计算方法为：

设机组轴系危险截面扭矩为T，轴系危险截面抗扭截面系数为W，危险截面应力集中系数为K，则机组轴系危险截面局部扭应力为求解得出轴系危险截面扭矩响应后，经计算得轴系危险截面扭应力响应，同时用应力集中系数对扭应力进行修正，得到危险截面局部扭应力响应。

7.根据权利要求1所述的方法，其特征在于，所述步骤3根据轴系扭振危险截面材料扭转S-N曲线和局部扭应力响应，采用雨流法和线性累积法对轴系扭振危险截面进行扭振疲劳寿命损耗计算。

8.根据权利要求1所述的方法，其特征在于，所述步骤4中的对应关系是指不同的电磁力矩最大标幺值条件下分别对应不同大扰动故障所引起的危险截面A和危险截面B分别所对应的轴系的扭振疲劳寿命损耗的百分比值。

9.根据权利要求1所述的方法，其特征在于，所述步骤5中确定机组轴系扭振损伤报警阈值的方法为：随着电磁力矩最大标幺值增大，查找首次出现不同大扰动故障所引起的危险截面A或危险截面B分别所对应的轴系的扭振疲劳寿命损耗的百分比值为非零的情况，该情况所对应的电磁力矩最大标幺值为机组轴系扭振损伤报警阈值；

确定机组轴系扭振跳机保护阈值的方法为：随着电磁力矩最大标幺值增大，查找首次出现不同大扰动故障所引起的的危险截面A或危险截面B分别所对应的轴系的扭振疲劳寿命损耗达到1％的情况，该情况所对应的电磁力矩最大标幺值为机组轴系扭振损伤报警阈值。

10.一种用于汽轮发电机组轴系大扰动瞬态冲击扭振监测的系统，其特征在于，包括：仿真模块、扭应力计算模块、疲劳寿命损耗计算模块、数据存储模块、阈值设定模块、电磁力矩监控模块、报警模块；

仿真模块依次通过扭应力计算模块、疲劳寿命损耗计算模块、数据存储模块、阈值设定模块、电磁力矩监控模块和报警模块相连；

其中，仿真模块的功能是：通过电力系统仿真，模拟电力系统大扰动故障下机组发生瞬态冲击类扭振时发电机的电磁力矩，仿真得到不同类型瞬态冲击类扭振故障下不同幅值的电磁力矩响应；

扭应力计算模块的功能是：根据仿真所得电磁力矩响应计算轴系扭振危险截面局部扭应力响应；

疲劳寿命损耗计算模块的功能是：对机组轴系扭振危险截面进行扭振疲劳寿命损耗计算；

数据存储模块的功能是：建立不同故障下机组电磁力矩最大值与轴系扭振危险截面扭振疲劳寿命损耗之间的对应关系；

阈值设定模块的功能是：确定机组轴系扭振损伤报警阈值和扭振跳机保护阈值；

电磁力矩监控模块的功能是：通过监测机组发电机三相电流和电压，实时在线计算相应电磁力矩；

报警模块的功能是：当发电机电磁力矩出现大幅跳变，电磁力矩值超过扭振损伤报警阈值时，对机组轴系发生的瞬态冲击类扭振进行报警；当电磁力矩值超过跳机保护阈值时，则判定机组轴系发生损伤程度较大的瞬态冲击类扭振，发出跳机保护信号。