CN109657397A - 基于频响函数的汽轮机叶片-转子系统稳定性的预测方法 - Google Patents

Abstract

本发明公开了基于频响函数的汽轮机叶片‑转子系统稳定性的预测方法，它解决了现有技术中转子系统无法准确稳定性预测的问题，具有能避免设计阶段无法准确预测机组稳定性造成的“先天不足”，将大幅提高设计效率降低设计成本的有益效果，其方案如下：基于频响函数的汽轮机叶片‑转子系统稳定性的预测方法，利用模态正交性求解大型多自由度汽轮机转子系统的频响函数，计算获得转子系统特征乘子，进而利用特征乘子与单位圆的关系判定系统稳定性。

Description

基于频响函数的汽轮机叶片-转子系统稳定性的预测方法

技术领域

本发明涉及非线性转子动力学领域，特别是涉及基于频响函数的汽轮机叶片-转子系统稳定性的预测方法。

背景技术

本部分的陈述仅仅是提供了与本公开相关的背景技术信息，不必然构成在先技术。

核电、热电联产等大型汽轮机组因绿色、高效等突出优势成为未来能源行业的主要发展趋势。转子系统一旦发生振动超标失稳，将对整个机组的安全“一票否决”，尤其大型汽轮机转子具有弹性长轴系、大柔性叶片等结构特点，且主要采用围带实现抑振增稳，在设计阶段准确预测转子系统稳定性对确保机组安全运行具有重要意义。

目前关于机组稳定性预测方法的报道多出现于论文中。文献中公开的大型机组稳定性预测方法主要分为以下两种：简化叶片的轴-轮盘稳定性预测方法和只考虑叶片的颤振预测方法。例如题为1000MW超超临界汽轮机汽流激振特性研究与应用、叶片-转子-轴承耦合系统的非线性动力学特性研究和高速透平机械稳定性评价与控制原理及方法研究等论文^[1-3]中忽略叶片柔性及几何非线性，将其简化为集中质量或单摆模型，无法考虑流场在叶片表面产生的压力脉动诱发的围带摩擦阻尼及其对机组稳定性的影响。题为二维叶栅气固耦合颤振分析、航空发动机叶片非线性动力学分析和基于流固耦合的叶片颤振分析等论文^[4-6]以及公开号为CN101908088A和CN101599104A等专利^[7,8]中仅给出了叶片的颤振预测方法，忽略了轴系结构，偏离机组实际结构和运行工况，无法真实预测整个机组转子系统的稳定性。

弹性长轴系、大柔性叶片及围带等结构特点及其与流场间的流固耦合作用使大型机组转子系统在不同运行工况下的动刚度、阻尼(尤其围带摩擦阻尼)具有强激励(含转速)依赖非线性，过度简化则无法求解。现有的转子系统稳定性预测及优化设计方法，不能综合考虑机组与流场的流固耦合作用以及轴系与叶片间变形耦合的影响，尤其对于采用围带结构增稳减振的转子系统，在设计阶段无法准确预测稳定性难以避免只凭经验设计的“先天不足”。

因此，需要对基于频响函数的汽轮机叶片-转子系统稳定性的预测方法进行新的研究设计。

发明内容

为了克服现有技术的不足，本发明提供了基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，有效避免设计阶段无法准确预测机组稳定性造成的“先天不足”，将大幅提高设计效率降低设计成本。

基于频响函数的汽轮机叶片-转子系统稳定性的预测方法的具体方案如下：

基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，建立考虑叶片柔性和围带摩擦阻尼的汽轮机叶片-转子系统运动微分方程，利用其模态正交性进行解耦，求解系统频响函数，计算系统特征乘子，进而利用特征乘子与单位圆的关系判定系统稳定性。该预测方法，丰富非线性转子动力学稳定性预测研究理论，更为在设计阶段准确预测稳定性以实现大型汽轮机组结构优化设计提供方法，可有效避免设计阶段无法准确预测机组稳定性造成的“先天不足”，降低相关设计成本。

进一步地，上述预测方法，包括如下内容：

1)建立考虑弹性转轴与柔性叶片的变形耦合非线性和流体激励诱发的围带摩擦阻尼激励依赖非线性的大型汽轮机叶片-转子系统运动微分方程；

2)获得系统运动微分方程中的流体激励F_f(t)，计算转子系统运动微分方程中的围带摩擦阻尼C_S(F,t)；

3)根据步骤2)获得的流体激励F_f(t)、围带摩擦阻尼C_S(F,t)，利用模态正交法解耦系统运动微分方程，求解综合考虑激励依赖围带阻尼、轴系叶片变形耦合影响的大型汽轮机叶片-转子系统频响函数；

4)通过叶片-转子系统频响函数计算系统特征乘子μ_F，并根据系统特征乘子μ_F与单位圆的关系判定系统在不同运行工况下的稳定性；

5)通过扫描所有运行工况，获得叶片-转子系统各运行工况的稳定性极限阈值，并将其组合绘制叶片-转子系统流致失稳三维极限曲面。

进一步地，所述步骤1)中所建立的系统运动微分方程是考虑围带摩擦阻尼的，运动微分方程表示为：

C(F,t)＝C_R+C₀+C_s(F,t)，K＝K_R+K_o (2)

式中，F(t)为系统的激励载荷矢量；Q(t)为位移矢量，M、C、K分别为系统质量、阻尼和刚度矩阵；其中，刚度矩阵中考虑了系统结构刚度K_R和支承刚度K_O；阻尼项中考虑了系统结构阻尼C_R、支承阻尼C_O及围带摩擦阻尼C_S(F,t)，且

F(t)＝F_f(t)+F_m(t) (4)

其中，F_f(t)为湿蒸汽非平衡凝结流场产生的流体激励，F_m(t)为转速激励和由不对中、不平衡等转子故障产生的机械激励。

进一步地，所述步骤1)所建立的系统运动微分方程中需辨识系统质量M、系统结构阻尼C_R、系统结构刚度K_R、系统支承阻尼C_O、系统支承刚度K_O，具体过程如下：

1-1)根据汽轮机组转轴的长度和直径、支承跨度、叶片参数(安装角、扭转角、展弦比)以及围带参数(接触角、接触间隙)，利用三维建模软件建立叶片围带接触角β、围带间初始间隙e₀的围带-叶片-转子系统几何模型，利用有限元分析软件，将系统离散为三维实体单元，得到式(1)-(2)中的系统质量、结构阻尼和结构刚度矩阵M、C_R、K_R以及单元形函数N；

1-2)将转子两端轴承支承近似为刚性支承，选择轴承参数(类型、载荷、润滑油型号和转速等)计算支承阻尼C_O和支承刚度K_O。

进一步地，所述步骤2)中利用数值方法获得系统运动微分方程中的流体激励F_f(t)，具体如下：

2-1)在计算流体软件中建立叶片流体计算域并划分网格，根据不同运行工况设置计算边界条件(进气速度V、进气角η)和转速Ω，选择适合的流体属性和湍流模型，进行非定常流场分析，获得叶片表面瞬态分布压力p_i(x,y,z,t)；

2-2)基于反距离加权插值法和虚功原理计算作用于叶片三维实体单元的流体激励F_f(t)：

F_f(t)＝N^TdF(t)^T (5)

式中，p_i(x,y,z,t)和d_i分别为流体计算域节点i(离散点)上的瞬态压力及其距叶片结构节点j(插值点)的距离；n为离散点的数量；N为叶片六面体单元型函数。

进一步地，所述步骤2)中围带摩擦阻尼C_S(F,t)的计算步骤如下：

基于FFT多谐波平衡法，将系统运动微分方程中的复杂多频激励项F(t)转换为不同频率成分谐波激励分量的叠加，假设每个谐波分量均激起叶片的简谐振动，则在第h阶谐波激励下，第i个叶片围带l侧上任一点在惯性坐标系(O-XYZ)下的振动位移响应q_h ⁱ(t)可写为式(8)：

式中，F₀为激励常量(直流量)，F_h、ω_h和φ_h分别为第h阶谐波激励分量及其频率和相位角，为相邻围带接触面间的相对运动幅值；

已知第h阶谐波激励下的围带摩擦阻尼为

δ＝ln(A1/A2) (10)

式中，是系统模态质量；ω_nr为忽略围带摩擦阻尼时系统的固有频率，δ为对数减幅率，A1和A2分别为系统振动响应在初始时刻和截止时刻的幅值；将各阶谐波分量下的等效摩擦阻尼叠加即得到复杂多频激励F(t)下的围带摩擦阻尼为

进一步地，所述步骤3)的具体过程如下：

根据步骤2)中求得的流体激励F_f(t)和求得的激励依赖围带摩擦阻尼C_S(F,t)代入系统运动微分方程中，即得综合考虑激励依赖围带阻尼、轴系叶片变形耦合影响的大型汽轮机叶片-转子系统的振动方程，其特征方程含有多(k)个特征根，即k个固有频率ω_n1，ω_n2，…ω_nr…ω_nk和相应的振型 根据振型正交性，有：

据此，解耦后的系统振动方程可写为：

式中，和分别为第r阶模态质量、模态刚度、模态阻尼和模态力，且有

则第r阶模态的频响函数(H(ω)_r)用实部和虚部表示为：

H(ω)_r＝ReH(ω)_r+ImH(ω)_r (14)

其中

式中，ω为系统激励角频率，ω_nr为系统第r阶模态固有频率，阻尼比ξ_r＝c_r/c_fr＝c_r/2m_rω_nr为系统阻尼c_r与第r阶临界阻尼系数c_fr之比；定义λr＝ω/ω_nr为对应于第r阶模态的频率比，代入式(13)-(15)中得大型多自由度汽轮机叶片-转子系统第r阶模态的频响函数(H(ω)_r)的实部和虚部可写为：

进一步地，所述步骤4)中系统特征乘子μ_F求解过程如下：

当系统以接近第r阶模态固有频率ω_nr的颤振频率ω_cr振动时，系统运动微分方程(式(13))的特征方程可写为：

det[I-diagΛ_rH(ω_cr)]＝0(r＝1，2，3，...，k) (19)

式中，I为k×k阶单位矩阵，diagΛr、H(ω_cr)为k×k阶对角阵。则第r阶模态的特征根为：

Λ^r＝[H(ω^cr)]^-1＝[H(ω)_r]^-1 (20)

据矩阵运算法则，复数矩阵(H(ω)_r)的逆矩阵可写为：

[H(ω)_r]^-1＝(A+iB)^-1

＝(A+BA^-1B)^-1-iA^-1B(A+BA^-1B)^-1 (21)

式中A、B代表频响函数的实部ReH(ω)_r与虚部ImH(ω)_r(式(18))；将式(17)-(18)代入式(20)-(21)中，即可计算获得系统第r阶模态特征根；

系统特征乘子(μ_F)为：

μ_F＝exp(Λ_rT) (22)

其中，T为主振周期，通常情况下T＝2πΩ/60。将第r阶模态特征根Λr代入(21)中即可获得大型多自由度汽轮机叶片-转子系统特征乘子。

进一步地，所述步骤4)的判定标准如下：

a)临界失稳：|μ_F|落在单位圆上；

b)Hopf分岔失稳：|μ_F|溢出单位圆的方式为共扼复数；

c)1周期分岔失稳：|μ_F|溢出单位圆的方式为正实数；

d)倍周期分岔失稳：|μ_F|溢出单位圆的方式为负实数；

若系统稳定，则返回步骤2)，进行下一组运行工况下系统稳定性判定；若系统失稳，则在三维坐标系下标记当前工况点。

进一步地，所述工况为进气速度V_i+ΔV、进气角η_i+Δη、转速Ω_i+ΔΩ，其中，ΔV、Δη和ΔΩ分别为各工况参数V_i、η_i和Ω_i相应的微增量。

与现有技术相比，本发明的有益效果是：

1)本发明通过预测方法的给出，解决了传统方法无法综合考虑大型机组转子系统因长轴系(弹性)、大柔性叶片及围带等结构特点及其与流场间的流固耦合作用导致的强激励(含转速)依赖非线性对机组稳定性的影响的缺陷，有效提高了多自由度非线性叶片-转子系统稳定性预测的精度。

2)本发明不仅丰富非线性转子动力学稳定性预测研究理论，更为在设计阶段准确预测稳定性以实现大型汽轮机组结构优化设计提供方法，可有效避免设计阶段无法准确预测机组稳定性造成的“先天不足”，将大幅提高设计效率降低设计成本。

附图说明

构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。

图1为本发明的基于频响函数的汽轮机叶片-转子系统稳定性预测方法流程图；

图2为本发明的一个具体实施例的某1000MW汽轮机叶片-转子系统示意图；

图3为本发明的一个具体实施例的末级叶片流体计算域示意图；

图4为本发明的一个具体实施例的末级叶片流场表面分布压力示意图；

图5为本发明的一个具体实施例的三维稳定性极限曲面图。

具体实施方式

应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

正如背景技术所介绍的，现有技术中存在的不足，为了解决如上的技术问题，本申请提出了基于频响函数的汽轮机叶片-转子系统稳定性的预测方法。

本申请的一种典型的实施方式中，如图2所示，基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，建立考虑叶片柔性和围带摩擦阻尼的汽轮机叶片-转子系统运动微分方程，利用其模态正交性进行解耦，求解系统频响函数，计算系统特征乘子，进而利用特征乘子与单位圆的关系判定系统稳定性。

上述方法，如图1所示，包括如下内容：

1)以某1000MW汽轮机组为具体实施例，如图2所示，建立考虑弹性转轴与柔性叶片的变形耦合非线性和流体激励诱发的围带摩擦阻尼激励依赖非线性的大型汽轮机叶片-转子系统运动微分方程：

C(F，t)＝C_R+C₀+C_s(F,t),K＝K_R+K₀ (2)

式中，F(t)为系统的激励载荷矢量；Q(t)为位移矢量，M、C、K分别为系统质量、阻尼和刚度矩阵；其中，刚度项中考虑了系统结构刚度K_R和支承刚度K_O；阻尼项中考虑了系统结构阻尼C_R、支承阻尼C_O及围带摩擦阻尼C_S(F,t)，且

F(t)＝F_f(t)+F_m(t) (4)

其中，F_f(t)为湿蒸汽非平衡凝结流场产生的流体激励，F_m(t)为转速激励和不对中、不平衡等转子故障产生的机械激励；

2)辨识系统运动微分方程(1)中的系统质量M、系统结构阻尼C_R、系统结构刚度K_R、系统支承阻尼C_O、系统支承刚度K_O，具体过程如下：

i)根据1000MW汽轮机转轴的长度和直径、支承跨度、叶片参数(安装角、扭转角、展弦比)以及围带参数(接触角、接触间隙)，利用三维建模软件建立末级叶片围带接触角β、围带间初始间隙e₀的围带-叶片-转子系统几何模型，利用有限元分析软件ANSYS Workbench中的网格划分功能，将系统离散为三维实体单元(Solid 185)，得到式(1)-(2)中的系统质量、结构阻尼和结构刚度矩阵M、C_R、K_R以及单元形函数N；

ii)将转子两端的可倾瓦滑动轴承支承近似为刚性支承，选择轴承参数(类型、载荷、润滑油型号和转速等)计算得支承阻尼C_O和支承刚度K_O为

3)利用数值方法获得系统运动微分方程(1)中的流体激励F_f(t)，具体过程如下：

i)在计算流体软件CFX中建立末级叶片流体计算域(如图3)并划分网格，根据不同运行工况设置计算边界条件(进气速度V、进气角η和转速Ω)，选择适合的流体属性和湍流模型，进行非定常流场分析，获得叶片表面瞬态分布压力p_i(x,y,z,t)，如图4所示；

ii)基于反距离加权插值法和虚功原理计算作用于叶片三维实体单元的流体激励F_f(t)：

F_f(t)＝N^TdF(t)^T (5)

式中，p_i(x,y,z,t)和d_i分别为流体计算域节点i(离散点)上的瞬态压力及其距叶片结构节点j(插值点)的距离；n为离散点的数量，具体实施例中取n＝4；N为叶片六面体单元形函数。

4)计算系统运动微分方程(1)中的围带摩擦阻尼C_S(F,t)，具体过程如下：

基于FFT多谐波平衡法，将式(1)中的复杂多频激励项F(t)转换为不同频率成分谐波激励分量的叠加(如式(7))，假设每个谐波分量均激起叶片的简谐振动，则在第h阶谐波激励下，第i个叶片围带l侧上任一点在惯性坐标系(O-XYZ)下的振动位移响应q_h ⁱ(t)可写为式(8)：

式中，F₀为激励常量(直流量)，因其不随时间而变，且对系统动态特性影响甚微，可忽略不计。F_h、ω_h和φ_h分别为第h阶谐波激励分量及其频率和相位角。为相邻围带接触面间的相对运动幅值。

已知第h阶谐波激励下的围带摩擦阻尼为

δ＝ln(A1/A2) (10)

式中，是系统模态质量；ω_nr为忽略围带摩擦阻尼时系统的固有频率，δ为对数减幅率，A1和A2分别为系统振动响应在初始时刻和截止时刻的幅值。将各阶谐波分量下的等效摩擦阻尼叠加即得到复杂多频激励F(t)下的围带摩擦阻尼为

工程中一般取前三阶谐波分量(n＝3)即可满足要求。

步骤4)中，利用有限元分析软件ANSYS Workbench中的模态分析模块，定义支承刚度K_O和支承阻尼C_o，施加转速等效离心载荷激励，得到不同工况下系统前r阶模态质量m_r和固有频率ω_nr；利用ANSYS中的瞬态分析模块，通过在围带接触面添加接触的方式定义摩擦系数μ和接触刚度系数K_h；定义支承刚度K_O和支承阻尼C_o；施加转速和各阶谐波激励载荷，得到系统振动响应曲线；在系统的最大振动位移响应曲线读取出系统振动响应在初始时刻和截止时刻的幅值A1和A2，代入式(10)计算对数减幅率δ，再将δ、m_r和ω_nr代入式(9)，即可求解出不同工况复杂多频激励诱发的围带摩擦阻尼C_s(F,t)。

5)求解综合考虑激励依赖围带阻尼、轴系叶片变形耦合影响的大型多自由度汽轮机叶片-转子系统频响函数，具体过程如下：

将步骤2)求得的M、C_R、K_R、C_O、K_O，步骤3)求得的流体激励F_f(t)和步骤4)求得的激励依赖围带摩擦阻尼C_S(F,t)代入系统运动微分方程(式(1))中，即得综合考虑激励依赖围带阻尼、轴系叶片变形耦合影响的大型汽轮机叶片-转子系统的振动方程。其特征方程含有多(k)个特征根，即k个固有频率ω_n1，ω_n2，…ω_nr…ω_nk和相应的振型根据振型正交性，有：

据此，解耦后的系统振动方程可写为：

式中，和分别为第r阶模态质量、模态刚度、模态阻尼和模态力。且有

则第r阶模态的频响函数(H(ω)_r)用实部和虚部表示为：

H(ω)_r＝ReH(ω)_r+ImH(ω)_r (14)

其中

式中，ω为系统激励角频率，ω_nr为系统第r阶模态固有频率，阻尼比ξ_r＝c_r/c_fr＝c_r/2m_rω_nr为系统阻尼c_r与第r阶临界阻尼系数c_fr之比；定义λr＝ω/ω_nr为对应于第r阶模态的频率比，代入式(13)-(15)中得大型多自由度汽轮机叶片-转子系统第r阶模态的频响函数(H(ω)_r)的实部和虚部可写为：

6)求解大型多自由度汽轮机叶片-转子系统特征乘子μ_F，具体过程如下：

当系统以接近第r阶模态固有频率ω_nr的颤振频率ω_cr振动时，系统运动微分方程(式(13))的特征方程可写为：

det[Ι-diagΛ_rH(ω_cr)]＝0(r＝1,2,3,...,k) (19)

式中，I为k×k阶单位矩阵，diagΛr、H(ω_cr)为k×k阶对角阵。则第r阶模态的特征根为：

Λ_r＝[H(ω_cr)]^-1＝[H(ω)_r]^-1 (20)

据矩阵运算法则，复数矩阵(H(ω)_r)的逆矩阵可写为：

[H(ω)_r]^-1＝(A+iB)^-1

＝(A+BA^-1B)^-1-iA^-1B(A+BA^-1B)^-1 (21)

式中，A、B代表频响函数的实部ReH(ω)_r(式(17))与虚部ImH(ω)_r(式(18))。将式(17)-(18)代入式(20)-(21)中，即可计算获得系统第r阶模态的特征根。

系统特征乘子(μ_F)为：

μ_F＝exp(Λ_rT) (22)

其中，T为主振周期，通常情况下T＝2πΩ/60。将第r阶模态特征根Λr代入(22)中即可获得大型多自由度汽轮机叶片-转子系统特征乘子。

7)据系统特征乘子μ_F与单位圆的关系判定系统在任一运行工况(进气速度V_i、进气角η_i和转速Ω_i)(i＝1,2,3,…..n)下的稳定性，具体判定准则为：

a)临界失稳：|μ_F|落在单位圆上；

b)Hopf分岔失稳：|μ_F|溢出单位圆的方式为共扼复数；

c)1周期分岔失稳：|μ_F|溢出单位圆的方式为正实数；

d)倍周期分岔失稳：|μ_F|溢出单位圆的方式为负实数。

若系统稳定，则返回步骤(3)，进行下一组运行工况(进气速度(V_i+ΔV)、进气角(η_i+Δη)、转速(Ω_i+ΔΩ))(其中，ΔV、Δη和ΔΩ分别为各工况参数V_i、η_i和Ω_i相应的微增量)下系统稳定性判定；若系统失稳，则在三维(进气速度、进气角和转速)坐标系下标记当前工况点(V_i、η_i和Ω_i)。

8)重复步骤7)，通过扫描所有运行工况(V、η、Ω)，获得叶片-转子系统各运行工况的稳定性极限阈值，并将其组合绘制叶片-转子系统流致失稳三维(V-η-Ω)极限曲面，如图5所示。

以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (10)

1.基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，建立考虑叶片柔性和围带摩擦阻尼的汽轮机叶片-转子系统运动微分方程，利用其模态正交性进行解耦，求解系统频响函数，计算系统特征乘子，进而利用特征乘子与单位圆的关系判定系统稳定性。

2.根据权利要求1所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，包括如下内容：

1)建立考虑弹性转轴与柔性叶片的变形耦合非线性和流体激励诱发的围带摩擦阻尼激励依赖非线性的大型汽轮机叶片-转子系统运动微分方程；

2)获得系统运动微分方程中的流体激励F_f(t)，计算系统运动微分方程中的围带摩擦阻尼C_S(F,t)；

3)根据步骤2)获得的流体激励F_f(t)、围带摩擦阻尼C_S(F,t)，利用模态正交法解耦系统运动微分方程，求解综合考虑激励依赖围带阻尼、轴系叶片变形耦合影响的大型汽轮机叶片-转子系统频响函数；

4)通过叶片-转子系统频响函数计算系统特征乘子μ_F，并根据系统特征乘子μ_F与单位圆的关系判定系统在不同运行工况下的稳定性；

5)通过扫描所有运行工况，获得叶片-转子系统各运行工况的稳定性极限阈值，并将其组合绘制叶片-转子系统流致失稳三维极限曲面。

3.根据权利要求2所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，所述步骤1)中所建立的系统运动微分方程是考虑围带摩擦阻尼的，运动微分方程表示为：

C(F，t)＝C_R+C_O+C_S(F，t)，K＝K_R+K_O (2)

式中，F(t)为系统的激励载荷矢量；Q(t)为位移矢量，M、C、K分别为系统质量、阻尼和刚度矩阵；其中，刚度矩阵中考虑了系统结构刚度K_R和支承刚度K_O；阻尼项中考虑了系统结构阻尼C_R、支承阻尼C_O及围带摩擦阻尼C_S(F,t)，且

F(t)＝F_f(t)+F_m(t) (4)

其中，F_f(t)为湿蒸汽非平衡凝结流场产生的流体激励，F_m(t)为转速激励和由不对中、不平衡等转子故障产生的机械激励。

4.根据权利要求3所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，所述步骤1)所建立的系统运动微分方程中需辨识系统质量M、系统结构阻尼C_R、系统结构刚度K_R、系统支承阻尼C_O、系统支承刚度K_O，具体过程如下：

1-1)根据汽轮机组转轴的长度和直径、支承跨度、叶片参数(安装角、扭转角、展弦比)以及围带参数(接触角、接触间隙)，利用三维建模软件建立叶片围带接触角β、围带间初始间隙e₀的围带-叶片-转子系统几何模型，利用有限元分析软件，将系统离散为三维实体单元，得到式(1)-(2)中的系统质量、结构阻尼和结构刚度矩阵M、C_R、K_R以及单元形函数N；

1-2)将转子两端轴承支承近似为刚性支承，选择轴承参数(类型、载荷、润滑油型号和转速等)计算支承阻尼C_O和支承刚度K_O。

5.根据权利要求2所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，所述步骤2)中利用数值方法获得系统运动微分方程中的流体激励F_f(t)，具体如下：

2-1)在计算流体软件中建立叶片流体计算域并划分网格，根据不同运行工况设置计算边界条件(进气速度V、进气角η、转速Ω)，选择适合的流体属性和湍流模型，进行非定常流场分析，获得叶片表面瞬态分布压力p_i(x,y,z,t)；

2-2)基于反距离加权插值法和虚功原理计算作用于叶片三维实体单元的流体激励F_f(t)：

F_f(t)＝N^TdF(t)^T (5)

式中，p_i(x,y,z,t)和d_i分别为流体计算域节点i(离散点)上的瞬态压力及其距叶片结构节点j(插值点)的距离；n为离散点的数量；N为叶片六面体单元型函数。

6.根据权利要求3所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，所述步骤2)中围带摩擦阻尼C_S(F,t)的计算步骤如下：

基于FFT多谐波平衡法，将系统运动微分方程中的复杂多频激励项F(t)转换为不同频率成分谐波激励分量的叠加，假设每个谐波分量均激起叶片的简谐振动，则在第h阶谐波激励下，第i个叶片围带l侧上任一点在惯性坐标系(O-XYZ)下的振动位移响应q_h ⁱ(t)可写为式(8)：

式中，F₀为激励常量(直流量)，F_h、ω_h和φ_h分别为第h阶谐波激励分量及其频率和相位角，为相邻围带接触面间的相对运动幅值；

已知第h阶谐波激励下的围带摩擦阻尼为

δ＝ln(A1/A2) (10)

式中，是系统模态质量；ω_nr为忽略围带摩擦阻尼时系统的固有频率，δ为对数减幅率，A1和A2分别为系统振动响应在初始时刻和截止时刻的幅值；将各阶谐波分量下的等效摩擦阻尼叠加即得到复杂多频激励F(t)下的围带摩擦阻尼为

7.根据权利要求2所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，所述步骤3)的具体过程如下：

根据步骤2)中求得的流体激励F_f(t)和求得的激励依赖围带摩擦阻尼Cs(F,t)代入系统运动微分方程中，即得综合考虑激励依赖围带阻尼、轴系叶片变形耦合影响的大型汽轮机叶片-转子系统运动微分方程，其特征方程含有多(k)个特征根，即k个固有频率ω_n1，ω_n2，…ω_nr…ω_nk和相应的振型 根据振型正交性，有：

据此，解耦后的系统运动微分方程可写为：

式中，和分别为第r阶模态质量、模态刚度、模态阻尼和模态力，且有

则第r阶模态的频响函数(H(ω)_r)用实部和虚部表示为：

H(ω)_r＝ReH(ω)_r+ImH(ω)_r (14)

其中

式中，ω为系统激励角频率，ω_nr为系统第r阶模态固有频率，阻尼比ξ_r＝c_r/c_fr＝c_r/2m_rω_nr为系统阻尼c_r与第r阶临界阻尼系数c_fr之比；定义λr＝ω/ω_nr为对应于第r阶模态的频率比，代入式(13)-(15)中得大型多自由度汽轮机叶片-转子系统第r阶模态的频响函数(H(ω)_r)的实部和虚部可写为：

8.根据权利要求2所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，所述步骤4)中系统特征乘子μ_F求解过程如下：

当系统以接近第r阶模态固有频率ω_nr的颤振频率ω_cr振动时，系统运动微分方程(式(13))的特征方程可写为：

det[Ι-diagΛ_rH(ω_cr)]＝0(r＝1,2,3,...,k) (19)

式中，I为k×k阶单位矩阵，diagΛr、H(ω_cr)为k×k阶对角阵。则第r阶模态的特征根为：

Λ_r＝[H(ω_cr)]^-1＝[H(ω)_r]^-1 (20)

据矩阵运算法则，复数矩阵(H(ω)_r)的逆矩阵可写为：

[H(ω)_r]^-1＝(A+iB)^-1

＝(A+BA^-1B)^-1-iA^-1B(A+BA^-1B)^-1 (21)

式中A、B代表频响函数的实部ReH(ω)_r与虚部ImH(ω)_r(式(18))；将式(17)-(18)代入式(20)-(21)中，即可计算获得系统第r阶模态特征根；

系统特征乘子(μ_F)为：

μ_F＝exp(Λ_rT) (22)

其中，T为主振周期，通常情况下T＝2πΩ/60。将第r阶模态特征根Λr代入(21)中即可获得大型多自由度汽轮机叶片-转子系统特征乘子。

9.根据权利要求2所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，所述步骤4)的判定标准如下：

a)临界失稳：|μ_F|落在单位圆上；

b)Hopf分岔失稳：|μ_F|溢出单位圆的方式为共扼复数；

c)1周期分岔失稳：|μ_F|溢出单位圆的方式为正实数；

d)倍周期分岔失稳：|μ_F|溢出单位圆的方式为负实数；

若系统稳定，则返回步骤2)，进行下一组运行工况下系统稳定性判定；若系统失稳，则在三维坐标系下标记当前工况点。

10.根据权利要求9所述的基于频响函数的汽轮机叶片-转子系统稳定性的预测方法，其特征在于，所述工况为进气速度V_i+ΔV、进气角η_i+Δη、转速Ω_i+ΔΩ，其中，ΔV、Δη和ΔΩ分别为各工况参数V_i、η_i和Ω_i相应的微增量。