CN111062177B - 一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法 - Google Patents

Abstract

本发明公开了一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法，通过围带摩擦阻尼计算模块、基于摄动修改动态优化围带摩擦阻尼模块、基于频响函数摄动修改围带参数优化稳定性模块、基于摄动修改动态优化供能/热效率模块和基于约束条件摄动修改工况参数优化模块依次进行数据处理，以最大稳定性的供能/热效率和对应工况下围带内阻尼最优为优化目标，以工况参数、围带参数为优化设计变量，综合优化达到对系统振动响应幅值的最大抑制效果。本发明操作流程清晰、计算效率高、结果准确，为汽轮机转子系统的优化设计提供理论依据和方法，对确保汽轮机安全运行具有重要理论意义和工程应用价值。

Description

一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法

技术领域

本公开涉及非线性转子动力学及稳定性动态优化，尤其涉及一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法。

背景技术

本部分的陈述仅仅是提到了与本公开相关的背景技术，并不必然构成现有技术。

具有长轴系、大柔性叶片及采用围带减振增稳等复杂结构特点的大型汽轮机组转子系统，在湿蒸汽两相流等流体激励下以及系统内阻尼等非线性因素引起的动静叶碰摩及叶片与轴的弯扭耦合等非线性振动和不稳定性是机组安全稳定运行的瓶颈限制，尤其随着大型机组逐渐向更高参数和轻量化方向发展，在设计阶段对整个运行周期的不同实际运行工况(包括升速启动、甩负荷停车、不同转速及不同流量等)复杂激励下转子系统稳定性进行准确的动态优化设计以避免凭经验设计的先天不足是十分必要的。

在实现本公开的过程中，发明人发现现有技术中存在以下技术问题：

目前现有技术中严重缺乏适合大型汽轮机组特点的转子系统稳定性动态优化设计方法，通常忽略叶片柔性和围带内阻尼等非线性特性的影响，优化设计变量大多只是转轴或轴承的特征尺寸，例如，现有文献1[黄晶晶.基于振动控制的柔性转子系统多目标优化技术研究.西安:西北工业大学,2012.]以转轴直径为设计变量、以质量最小为目标函数，对含轴承支承的阶梯轴转子系统采用遗传算法进行优化设计。现有文献2[崔立,郑建荣,周炜.考虑转子系统耦合影响的球轴承动态性能多目标优化设计.振动与冲击.2012,31(24):190-196.]采用自适应变异遗传优化算法(AMGA)对某转子系统过一阶临界转速时两转盘的最大振幅进行优化，但将优化对象简化为仅由转轴、刚性转盘和轴承组成的转子系统，忽略围带及柔性结构等非线性因素，且系统激励只考虑了转子所受的不平衡载荷。现有文献3[于瑾,高建,任朝辉.转子系统的动态优化设计研究.机械设计与制造.2012,4:241-243.]采用非劣排序遗传算法(NSGA-II)对转子系统的支承轴承额定动刚度、径向刚度和旋滚比进行优化，将优化对象简化为由刚性圆盘、弹性轴段和轴承构成的转子，考虑了转轴与轴承的耦合影响，但无法考虑柔性叶片、围带及其内阻尼等非线性影响，与实际结构相比偏差较大。

近年来，基于结构动力修改理论结合有限元软件在转子系统稳定性动态优化中得到了较广泛应用，如用模态分析结果进行修改的设计参数修改法、基于频响函数灵敏度分析的有限元模型修正法等。现有文献4[宋清华.高速旋转系统动态特性分析及优化方法研究.济南:山东大学,2008.]以考虑不平衡激励的柔性轴-弹性圆盘-刚性滑动支承旋转系统为优化对象，目标函数为能反映系统稳定性的频响函数实部和系统振动位移响应，以转轴和圆盘的结构尺寸为设计变量并分析其对频响函数的灵敏度，运用有限元模型的频响函数摄动修改法对系统结构进行多目标动力修改优化。

上述关于转子系统的动态优化均未考虑柔性叶片以及叶片振动诱发的围带摩擦产生的阻尼作用对系统稳定性的影响。现仅有文献5[刘锦.核电汽轮机转子叶片围带内阻尼定量计算与动态优化研究.济南:山东大学,2017.]提出定量计算围带摩擦阻尼的方法并采用摄动修改对叶片-围带-转子系统进行了动态优化，但优化目标函数仅为转速单一简谐激励下的围带摩擦阻尼最优，忽略了大型核/火电汽轮机组实际运行工况下湿蒸汽等复杂流体多频激励等非线性因素对系统稳定性的影响。而有围带汽轮机转子系统稳定性动态优化，是涉及多种转子结构非线性、围带阻尼激励依赖非线性、且需要在整个运行周期的多个运行工况(包括升速启动、甩负荷停车、不同转速及不同流量等)的复杂激励等多变量耦合影响下保证稳定运行的有约束的多目标动态优化问题，目前没有适合大型汽轮机组特点涉及多个复杂非线性变量耦合影响的优化算法，尚缺乏能在设计阶段进行汽轮机转子系统稳定性动态优化设计的方法。

发明内容

为了解决现有技术的不足，本公开提供了一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法；

第一方面，本公开提供了一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法；

一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法，包括以下步骤：

步骤(1)：获取汽轮机转子系统的运行工况参数；建立转子系统稳定性判定准则和稳定性极限阈值的求解方法；

步骤(2)：确定多目标动态优化的目标函数及其优化变量：

目标函数1，f₁(x)：稳定运行下汽轮机组供电/供热效率最大，涉及的优化变量为运行工况参数：转子运行转速Ω，相对进气流量G/G₀、进气速度V和进气角η；其中，G为实际进气流量，G₀为额定进气流量；

目标函数2，f₂(x)：转子系统运行稳定性极限阈值区间最优，即转子系统的围带阻尼最优，涉及的优化变量为围带参数：相邻围带间的摩擦系数μ、接触角β、接触刚度K_h和初始间隙e₀；

步骤(3)：设定约束条件：以供电/供热功率、效率曲线、各部件强度、疲劳寿命、无共振、强迫振动响应小于设定振幅、以及只供电或热电联产实际运行工况为约束条件；

步骤(4)：对转子系统稳定性进行有约束的多目标动态优化。

进一步地，在所述步骤(1)中，获取汽轮机转子系统的运行工况参数；具体包括：转子运行转速Ω、实际进气流量G、进气速度V、进气角η和额定进气流量G₀。

进一步地，在所述步骤(1)中，建立转子系统稳定性判定准则和稳定性极限阈值的求解方法，具体步骤如下：

步骤(1-1)：建立转子系统稳定性判定准则：基于Floquet理论，系统方程单值矩阵的特征乘子存在等于1的解(|μ_F|＝1)，系统处于临界稳定状态(即|μ_F|＝1为稳定性极限)，同时若其他特征乘子|μ_F|均小于1，则系统失稳；汽轮机组转子系统稳定性可依据特征乘子与单位圆的关系判定；

步骤(1-2)：建立稳定性极限阈值的求解方法。

进一步地，步骤(1-2)的具体步骤如下：

步骤(1-2-1)：基于Floquet理论，特征乘子μ_F可表示为：

μ_F＝exp(ΛT) (1)

其中，T＝Im(lnμ_F)/ω为主振动周期；Λ为系统特征方程的特征根。

将转子系统运动微分方程写为复坐标形式：

-m_rω_n ²+c_riω_n+k_r＝F₀exp(iωt) (2)

其中，ω为激励频率；ω_n为系统固有频率；i²＝-1；则转子系统第r阶频响函数H(ω)_r写为以下实部-虚部形式：

其中，定义频率比λ_r＝ω/ω_nr为系统所受激励频率与系统第r阶固有频率之比；阻尼比

为系统阻尼C_r与第r阶临界阻尼(C_rb＝2m_rω_nr)之比，阻尼C_r＝C_l+C_s，对于有围带的转子系统由于系统的结构阻尼C_l远小于围带摩擦阻尼C_s，为方便求解，忽略结构阻尼C_l，则有C_r＝C_s；kr为第r阶模态刚度。

转子系统运动微分方程(2)的特征方程为：

det{I-Λ[H(ω)]}＝0 (4)

则系统特征方程(4)的特征根Λ可写为

Λ＝H(ω)^-1 (5)

由频响函数(式(3))中所含各变量可知，对于具有几何非线性柔性叶片和激励依赖非线性围带阻尼的大型汽轮机转子系统，在不同运行工况下，其频率比

阻尼比

(各部件阻尼等)呈现不同的非线性特性。依据特征乘子与单位圆的关系判定转子系统稳定性,获得转子系统特征方程(式(4))的特征根Λ从而求解特征乘子|μ_F|(式(1))，由式(5)知特征根Λ为系统频响函数(式(3))的逆，则系统稳定性的判定可归结为求解系统在不同运行工况下的频响函数问题。

进一步地，据公式(3)，频响函数涉及系统各阶固有频率ω_nr、模态质量m_r、刚度k_r和阻尼C_r，其中阻尼C_r＝C_l+C_s，对于有围带的转子系统由于系统的结构阻尼C_l远小于围带阻尼C_s，为方便求解，忽略结构阻尼C_l，则有C_r＝C_s；分别求解各项以获得频响函数。

步骤(1-2-2)：计算各阶固有频率ω_nr、模态质量m_r和刚度k_r：

步骤(1-2-2-1)：建立几何模型并离散划分网格：根据设计的汽轮机组转轴的长度和直径、轴承支承跨度、叶片参数(安装角、扭转角、展弦比)以及围带参数(接触角β、接触间隙e₀)等建立转子系统三维几何模型，并将其导入有限元分析软件(以下简称FEMS)，在模态分析模块(Modal)中的网格离散模块(mesh)中将转子系统离散为三维实体单元；

步骤(1-2-2-2)：设置边界条件：在FEMS的模态分析(Modal)模块中的支承参数设置选项(Fix support)中利用滑动轴承参数计算软件(XROTOR)计算获得支承阻尼C_O和支承刚度K_O，在轴承支承处添加弹簧接触(connection-spring)，设置C_O和K_O；

步骤(1-2-2-3)：在FEMS的模态分析模块(Modal)中求解转子系统模态质量m_r和刚度k_r；再对整个转子系统施加转速(insert-Rotational Velocity)，求解(solution)得到忽略围带阻尼的固有频率ω_nr；

步骤(1-2-3)：计算获得围带阻尼C_s；

据多频激励下大型汽轮机转子围带摩擦阻尼等效模型和基于瞬态响应包络线的围带阻尼定量求解方法，多频激励F(t)下的围带阻尼C_s为：

其中，h＝1,2,3,......,p；p为正整数；

(F,t)为多频激励力F(t)中第h阶谐波激励分量作用下的围带阻尼：

其中，m_r是转子系统模态质量；ω_nr为忽略围带摩擦阻尼的系统固有频率，

和

分别为多频激励F(t)的第h阶谐波激励分量作用下转子系统最大振动位移响应节点的相邻位移响应峰值振幅，δ^h为对数减幅比；

步骤(1-2-3-1)：计算获得系统激励F(t)：运行在湿蒸汽等流体环境的大型汽轮机转子系统，所受的激励力F(t)由两部分组成：

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

其中，F_m(t)是相应运行转速下由动不平衡效应等引起的机械激励；F_f(t)是湿蒸汽等流体介质在叶片和围带表面的瞬态分布压力及脉动等流体多频激励；

步骤(1-2-3-1-1)：机械激励F_m(t)通过在FEMS的瞬态响应求解模块转速激励输入选项中选择运行转速作为F_m(t)激励输入；

步骤(1-2-3-1-2)：流体激励F_f(t)，利用基于流场仿真的流体激励求解方法得到如式(10)和(11)所示的流体激励F_f(t)：

F_f(t)＝N^Tdf(t)^T (10)

其中，p_i(x,y,z,t)和d_i分别为流场计算域中任意节点i(离散点)上在不同运行工况(进气流量G、进气速度V、进气角η和转速Ω)下的瞬态压力及其距叶片结构节点j(插值点)的距离；n为离散点的数量；N为叶片六面体单元形函数。

步骤(1-2-3-1-3)：利用基于FFT的多谐波平衡法确定该流体激励的前h阶频率成分(可取h≥5)以及各阶谐波频率下的谐波激励分量，将如式(10)和(11)表示的F_f(t)转化为一个常量和若干个不同频率成分谐波激励分量的叠加，如式(12)所示，可在FEMS的瞬态响应求解(Transient)模块中作为激励函数输入。

步骤(1-2-3-2)：获得系统多频激励F(t)下转子系统的最大振动位移响应，确定对数减幅比，具体步骤如下：

步骤(1-2-3-2-1)：设置边界条件：在FEMS的瞬态响应求解模块中的模型(Model)选项，在相邻围带接触面间添加摩擦接触(insert-connection-fraction)，设置摩擦系数μ和接触刚度K_h；

步骤(1-2-3-2-2)：施加激励并求解：在FEMS的瞬态响应求解模块施加F(t)，具体的，F_m(t)在转速激励输入选项(inertial)中选择运行转速作为激励输入；F_f(t)在载荷激励输入选项(insert-Force)中以

的形式施加由多阶谐波激励分量组成的多频激励力；在求解(solution)选项中求解第h阶谐波激励分量下的转子系统振动位移响应；重复上述步骤，依次分别求解出各阶谐波激励分量下的转子系统振动位移响应。

步骤(1-2-3-2-3)：确定第h阶谐波激励分量下的最大振动位移响应节点，输出该节点的最大振动位移响应曲线，读取相邻振动位移响应峰值的振幅

和

并依次代入式(10)求解出各阶振动位移响应对数减幅比δ^h；

步骤(1-2-3-2-4)：将步骤(1-2-2)所获转子系统质量m_r和忽略围带阻尼的固有频率ω_nr及步骤(1-2-3-2-3)所获的各阶对数减幅比δ^h代入式(7)计算围带阻尼

并将其代入式(6)求解得到多频激励F(t)下的围带阻尼C_s(F,t)。

步骤(1-2-4)：利用有限元分析软件求解转子系统频响函数并获取特征乘子：

据步骤(1-2-1)～(1-2-3)获得的转子系统各阶固有频率ω_nr、模态质量m_r、刚度k_r和围带阻尼C_s在FEMS的谐响应分析模块中求解获得转子系统的实部-虚部形式的频响函数(H(ω)_r)，并读取各阶特征根Λ_r(r为正整数)，将其代入式(1)即获得转子系统各阶特征乘子|μ_F|_r。

进一步地，在所述步骤(2)中，确定多目标动态优化的目标函数和优化设计变量：

目标函数1：f₁(x)：汽轮机组的供电/热总效率达到最大值，即

f₁(x)＝max[L_H ^s]＝max[BG_H ^s] (13)

其中，上标s表示当前运行工况是稳定的；L_H为该运行工况下输出的供电/热总效率，可表示为

L_H＝BG_H (14)

其中，G_H为实际用于供电/热的进气流量(kg/h)，按汽轮机组的两种不同用途划分为两大类运行工况：

1)常规供电：低压缸的相对进气流量G/G₀＝1(其中,G为实际进气流量、G₀为额定进气流量)为额定进气流量工况，G_H＝0；

2)热电联产：对于无备用低压缸的抽气式热电联产机组，供热季低压缸的相对进气流量G/G₀≤1(称小流量工况)，供热季实际用于供热的流量G_H＝G₀-G。B为单位进气流量G_H下考虑运行过程中能量转换消耗及效率等因素的产能效率系数。

确定f₁(x)的优化设计变量：优化设计变量是可控可调的初始变量，在设计阶段已知功率、效率曲线约束条件下，选取整个运行周期不同运行工况的工况参数：相对进气流量G/G₀及其相应的进气速度V和进气角η，转子运行转速Ω作为优化设计变量；

目标函数2：f₂(x)：转子系统运行稳定性极限阈值区间最优，由步骤(1-2)中稳定性极限阈值求解方法可知，对于有围带起主要减振增稳作用的大型汽轮机组，可通过优化围带阻尼进而优化转子系统稳定性极限阈值-特征乘子μ_F，故取优化目标函数2为各稳定运行工况下的围带阻尼最优，即：

其中，上标s表示稳定的运行工况；C_s为该运行工况下转子系统围带阻尼；μ、K_h、β和e₀依次分别为相邻围带接触面的摩擦系数、接触刚度系数、接触角和初始间隙。

确定f₂(x)的优化设计变量：在设计阶段，转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数初步设计确定的前提下，目标函数2的优化设计变量为围带参数：μ，K_h，β和e₀。

进一步地，在所述步骤(3)中：以机组整个运行周期各运行工况的供电/供热功率和效率曲线、各部件强度/刚度、疲劳寿命、无共振、强迫振动响应小于设定振幅为约束条件如下：

a.依据所设计的机组在各运行工况的供电/供热功率和效率曲线进行转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数的初步设计确定。并在优化设计过程中，各优化设计变量及其摄动修改增量的选取均应满足机组的功率和效率曲线；

b.为避免强度失效，各部件在不同运行工况载荷下的各类应力/应变不超过材料设定值,即满足：σ_imax(x,y,z)≥[σ(x,y,z)]；

c.无共振发生：运行转速Ω避开0.75-1.25倍的各阶临界转速Ω_r(r＝1,2,3,.....p,p为正整数)区间，即满足：Ω≤0.75Ω_r或Ω≥1.25Ω_r；

d.为避免动静件碰摩故障，转子系统各节点的各向最大振动位移响应q_imax(x,y,z)不超过设计设定振幅值(参照相应的设计手册),即满足：q_imax(x,y,z)≤[q(x,y,z)]；

e.为适应实际工况，转子系统的转速Ω小于等于常规运行转速的1.2倍；系统的发电/热功率P大于等于额定发电/热功率。

进一步地，在所述步骤(4)中，对转子系统稳定性进行有约束的多目标动态优化，具体步骤如下：

步骤(4-1)：确定满足目标函数2的围带参数：给定第i组运行工况参数(转速Ω_i，实际进气流量G_i(V_i,η_i))，并根据步骤(1-2)中的围带阻尼求解方法求解确定该运行工况下围带阻尼最优的第j组围带参数(摩擦系数μ_j,接触刚度K_hj,接触角β_j,初始间隙e_0j)，满足目标函数2；

步骤(4-2)：利用步骤(1)中的稳定性判定准则判定第i组运行工况参数(Ω_i，G_i(V_i,η_i))下具有第j组围带参数(μ_j，K_hj，β_j，e_0j)的转子系统的稳定性，若系统失稳，则进入步骤(4-2-1)；若稳定，则进入步骤(4-2-2)；

步骤(4-2-1)：在不同运行工况的三维坐标系(Ω,V,η)下标记当前工况点(Ω_i,V_i,η_i)为具有第j组围带参数的转子系统在该工况下的失稳区；

步骤(4-2-2)：判断第i组运行工况(Ω_i,G_i(V_i,η_i))下的供电/热效率L^s _H是否满足目标函数1(公式(13)～公式(14))，若不满足，则辨识确定敏感优化设计变量，在满足步骤(3)约束条件下选取一组微摄动变量(ΔΩ，ΔV，Δη)，摄动修改第i组运行工况参数(其中，1)纯供电机组，微摄动变量ΔV＝0和Δη＝0，ΔΩ≠0；2)热电联产机组供热季:微摄动变量ΔV、Δη、ΔΩ≠0)为第i+1组运行工况(V_i+1＝V_i+ΔV，η_i+1＝η_i+Δη，Ω_i+1＝Ω_i+ΔΩ)，并返回步骤(4-1)计算在第i+1组运行工况下围带阻尼最优的第j+1组围带参数(μ_j+1,K_hj+1,β_j+1,e_0j+1)和步骤(4-2)判定第i+1组运行工况下是否稳定；若不稳定，进入步骤(4-2-1)；若稳定，则进行步骤(4-2-3)；

步骤(4-2-3)：第j+1组围带参数(μ_j+1，K_hj+1，β_j+1，e_0j+1)即为第i+1组运行工况下满足两个目标函数的最优围带参数；记录第j+1组围带参数以及对应的最优供电/热效率和稳定性极限阈值；

步骤(4-3)：重复步骤(4-1)到步骤(4-2)扫描所有运行工况，获得动态优化后的转子系统在各运行工况的稳定性极限阈值，并将其组合绘制得到转子系统的三维稳定性极限阈值曲线。

与现有技术相比，本公开的有益效果是：

1.本公开提出的基于围带阻尼的汽轮机转子系统稳定性动态优化方法，适合大型汽轮机转子系统包含多非线性变量及其耦合影响等特点，克服了现有技术中对柔性叶片、围带等结构及不同运行工况下的复杂激励等过度简化导致严重偏离实际的弊端，可大幅提高含多非线性变量的大型汽轮机组转子系统的优化设计精度，为在设计阶段对转子系统稳定性进行有约束的多目标动态优化设计提供了方法，可保证机组在整个运行周期的不同运行工况(包括升速启动、甩负荷停车、不同转速及不同流量等)复杂激励下稳定运行，避免凭经验设计的先天不足。

2.本公开提出的基于围带阻尼的汽轮机转子系统稳定性动态优化方法，以最大稳定性的供能/热效率和稳定性极限阈值最优为优化目标，采用微量摄动修改多优化变量(围带参数和运行工况参数)，通过优化围带阻尼特性提升系统的稳定运行极限阈值、有效扩大稳定运行区域，从本公开实施例可看出，优化后的稳定运行区间扩大40％，其中，进气速度的稳定性极限阈值增大55％-77％，对有流量变化机组的稳定运行优化效果十分显著。

3.本公开提出的基于围带阻尼的转子系统稳定性动态优化方法，从对转子系统减振增稳起关键作用的围带阻尼入手，利用基于围带内阻尼的频响函数来预测判定系统稳定性，通过敏感变量辨识和采用微量摄动修改法对多变量进行动态优化，不仅提高了优化设计精度，还可大幅度缩减含多非线性优化变量的动态优化问题的计算规模、简化计算过程，方便工程应用。

附图说明

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

图1为本公开的一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法流程图；

图2为本公开的一个具体实施例的汽轮机转子系统有限元模型；

图3为本公开的一个具体实施例的汽轮机转子系统优化后的三维(Ω-η-V)稳定性极限曲面；

图4为本公开的一个具体实施例的汽轮机转子系统优化前的三维(Ω-η-V)稳定性极限曲面；

图5(a)-图5(c)为本公开的一个具体实施例的汽轮机转子系统优化后不同进气角η下的转速-进气速度稳定性极限。

图6为本公开的系统最大位移响应曲线。

具体实施方式

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

实施例一，如图1所示，本实施例提供了一种基于围带阻尼的汽轮机转子系统动态优化方法，包括：

步骤(1)：获取汽轮机转子系统的运行工况参数；建立转子系统稳定性判定准则和稳定性极限阈值的求解方法；

步骤(2)：确定多目标动态优化的目标函数及其优化设计变量：

目标函数1，f₁(x)：稳定运行下汽轮机组供电/供热效率最大，涉及的优化变量主要为运行工况参数：转子运行转速Ω，相对进气流量G/G₀(G为实际进气流量，G₀为额定进气流量)及其相应的进气速度V和进气角η，等；

目标函数2，f₂(x)：转子系统运行稳定性极限阈值区间最优，即转子系统的阻尼(主要包括围带阻尼)最优，涉及的优化变量主要为围带参数：相邻围带间的摩擦系数μ、接触角β、接触刚度K_h和初始间隙e₀，等；

步骤(3)：设定约束条件：以供电/供热功率、效率曲线、各部件强度、疲劳寿命、无共振、强迫振动响应小于设定振幅、以及只供电或热电联产等实际运行工况为约束条件；

步骤(4)：对转子系统稳定性进行有约束的多目标动态优化。

进一步地，在所述步骤(1)中，建立转子系统稳定性判定准则和稳定性极限阈值的求解方法，具体步骤如下：

步骤(1-1)：建立转子系统稳定性判定准则：基于Floquet理论，系统方程单值矩阵的特征乘子存在等于1的解(|μ_F|＝1)，系统处于临界稳定状态(即|μ_F|＝1为稳定性极限)，同时若其他|μ_F|均小于1，则系统失稳。汽轮机组转子系统稳定性可依据特征乘子与单位圆的关系判定。

步骤(1-2)：建立稳定性极限阈值的求解方法；

进一步地，步骤(1-2)的具体步骤如下：

步骤(1-2-1)：基于Floquet理论，特征乘子μ_F可表示为：

μ_F＝exp(ΛT) (1)

其中，T＝Im(lnμ_F)/ω为主振动周期；Λ为系统特征方程的特征根。

将转子系统运动微分方程写为复坐标形式：

-m_rω_n ²+c_riω_n+k_r＝F₀ exp(iωt) (2)

其中，ω为激励频率；ω_n为系统固有频率；i²＝-1；则转子系统第r阶频响函数H(ω)_r可写为以下实部-虚部形式：

其中，定义频率比λ_r＝ω/ω_nr为系统所受激励频率与系统第r阶固有频率之比；阻尼比

为系统阻尼C_r与第r阶临界阻尼(C_rb＝2m_rω_nr)之比，其中,阻尼C_r＝C_l+C_s，对于有围带的转子系统由于系统的结构阻尼C_l远小于围带摩擦阻尼C_s，为方便求解，忽略结构阻尼C_l，则有C_r＝C_s；k_r为第r阶模态刚度。

转子系统运动微分方程(2)的特征方程为：

det{I-Λ[H(ω)]}＝0 (4)

则系统特征方程(4)的特征根Λ可写为：

Λ＝H(ω)^-1 (5)

由频响函数(式(3))中所含各变量可知，对于具有几何非线性柔性叶片和激励依赖非线性围带阻尼的大型汽轮机转子系统，在不同运行工况下，其频率比

阻尼比

(各部件阻尼等)呈现不同的非线性特性，使系统的频响特性十分复杂。依据特征乘子与单位圆的关系判定转子系统稳定性,关键是获得转子系统特征方程(式(4))的特征根Λ从而求解特征乘子|μ_F|(式(1))，由式(5)知特征根Λ为系统频响函数(式(3))的逆，则系统稳定性的判定可归结为求解系统在不同运行工况下的频响函数问题。

进一步地，据式(3)，频响函数涉及系统各阶固有频率ω_nr、模态质量m_r、刚度k_r和阻尼C_r，其中阻尼C_r＝C_l+C_s,对于有围带的转子系统由于系统的结构阻尼C_l远小于围带阻尼C_s，为方便求解，忽略结构阻尼C_l，则有C_r＝C_s。需分别求解各项以获得频响函数。

步骤(1-2-2)：计算各阶固有频率ω_nr、模态质量m_r和刚度k_r：

步骤(1-2-2-1)：建立几何模型并离散划分网格：根据设计的汽轮机组转轴的长度和直径、轴承支承跨度、叶片参数(安装角、扭转角、展弦比)以及围带参数(接触角β、接触间隙e₀)等建立转子系统三维几何模型，并将其导入有限元分析软件(以下简称FEMS)，在模态分析模块(Modal)中的网格离散模块(mesh)中将转子系统离散为三维实体单元；

步骤(1-2-2-2)：设置边界条件：在FEMS的模态分析(Modal)模块中的支承参数设置选项(Fix support)中利用滑动轴承参数计算软件(XROTOR)计算获得支承阻尼C_O和支承刚度K_O，在轴承支承处添加弹簧接触(connection-spring)，设置C_O和K_O；

步骤(1-2-2-3)：在FEMS的模态分析模块(Modal)中求解转子系统模态质量m_r和刚度k_r；再对整个转子系统施加转速(insert-Rotational Velocity)，求解(solution)得到忽略围带阻尼的固有频率ω_nr；

步骤(1-2-3)：计算获得围带阻尼C_s；

据多频激励下大型汽轮机转子围带摩擦阻尼等效模型和基于瞬态响应包络线的围带阻尼定量求解方法，多频激励F(t)下的围带阻尼C_s为：

其中，

为多频激励力F(t)中第h阶谐波激励分量作用下的围带阻尼：

其中，m_r是转子系统模态质量；ω_nr为忽略围带摩擦阻尼的系统固有频率，

和

分别为多频激励F(t)的第h阶谐波激励分量作用下转子系统最大振动位移响应节点的相邻位移响应峰值振幅，δ^h为对数减幅比；

步骤(1-2-3-1)：计算获得系统激励F(t)：运行在湿蒸汽等流体环境的大型汽轮机转子系统，所受的激励力F(t)由两部分组成：

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

其中，F_m(t)是相应运行转速下由动不平衡效应等引起的机械激励；F_f(t)是湿蒸汽等流体介质在叶片和围带表面的瞬态分布压力及脉动等流体多频激励；

步骤(1-2-3-1-1)：机械激励F_m(t)可通过在FEMS的瞬态响应求解模块转速激励输入选项中选择运行转速作为F_m(t)激励输入；

步骤(1-2-3-1-2)：流体激励F_f(t)，作为实施方式之一，可利用基于流场仿真的流体激励求解方法得到如式(10)和(11)所示的流体激励F_f(t)：

F_f(t)＝N^Tdf(t)^T (10)

其中，p_i(x,y,z,t)和d_i分别为流场计算域中任意节点i(离散点)上在不同运行工况(进气流量G、进气速度V、进气角η和转速Ω)下的瞬态压力及其距叶片结构节点j(插值点)的距离；n为离散点的数量，具体实施例中可取n＝4；N为叶片六面体单元形函数。

步骤(1-2-3-1-3)：利用基于FFT的多谐波平衡法确定该流体激励的前h阶频率成分(可取h≥5)以及各阶谐波频率下的谐波激励分量，将如式(10)和(11)表示的F_f(t)转化为一个常量和若干个不同频率成分谐波激励分量的叠加，如式(12)所示，可在FEMS的瞬态响应求解(Transient)模块中作为激励函数输入。

步骤(1-2-3-2)：获得系统多频激励F(t)下转子系统的最大振动位移响应，确定对数减幅比，具体步骤如下：

步骤(1-2-3-2-1)：设置边界条件：在FEMS的瞬态响应求解模块中的模型(Model)选项，在相邻围带接触面间添加摩擦接触(insert-connection-fraction)，设置摩擦系数μ和接触刚度K_h；

步骤(1-2-3-2-2)：施加激励并求解：在FEMS的瞬态响应求解模块施加F(t)，具体的，F_m(t)在转速激励输入选项(inertial)中选择运行转速作为激励输入；F_f(t)在载荷激励输入选项(insert-Force)中以

的形式施加由多阶谐波激励分量组成的多频激励力；在求解(solution)选项中求解第h阶谐波激励分量下的转子系统振动位移响应；重复上述步骤，依次分别求解出各阶谐波激励分量下的转子系统振动位移响应。

步骤(1-2-3-2-3)：确定第h(h＝1,2,3,…p,其中，p为正整数)阶谐波激励分量下的最大振动位移响应节点，输出该节点的最大振动位移响应曲线，读取相邻振动位移响应峰值的振幅

和

并依次代入式(10)求解出各阶振动位移响应对数减幅比δ^h；

步骤(1-2-3-2-4)：将步骤(1-2-2)所获转子系统质量m_r和忽略围带阻尼的固有频率ω_nr及步骤(1-2-3-2-3)所获的各阶对数减幅比δ^h代入式(7)计算围带阻尼

并将其代入式(6)求解得到多频激励F(t)下的围带阻尼C_s。

步骤(1-2-3)：利用有限元分析软件求解转子系统频响函数并获取特征乘子：

据步骤(1-2-1)-(1-2-2)获得的转子系统各阶固有频率ω_nr、模态质量m_r、刚度k_r和围带阻尼C_s在FEMS的谐响应分析模块中求解获得转子系统的实部-虚部形式的频响函数(H(ω)_r)，并读取各阶特征根Λ_r，将其代入式(1)即获得转子系统各阶特征乘子|μ_F|_r。

进一步地，在所述步骤(2)中，确定多目标动态优化的目标函数和优化设计变量：

目标函数1：f₁(x)：汽轮机组的供电/热总效率达到最大值，即

f₁(x)＝max[L_H ^s]＝max[BG_H ^s] (13)

其中，上标s表示当前运行工况是稳定的；L_H为该运行工况下输出的供电/热总效率，可表示为：

L_H＝BG_H (14)

其中，G_H为实际用于供电/热的进气流量(kg/h)，按汽轮机组的两种不同用途划分为两大类运行工况：1)常规供电：低压缸的相对进气流量G/G₀＝1(其中,G为实际进气流量、G₀为额定进气流量)为额定进气流量工况，G_H＝0；2)热电联产：对于无备用低压缸的抽气式热电联产机组，供热季低压缸的相对进气流量G/G₀≤1(称小流量工况)，供热季实际用于供热的流量G_H＝G₀-G。B为单位进气流量G_H下考虑运行过程中能量转换消耗及效率等因素的产能效率系数。

确定f₁(x)的优化设计变量：优化设计变量一定是可控可调的初始变量，在设计阶段已知功率、效率曲线等约束条件下，选取整个运行周期不同运行工况的工况参数：相对进气流量G/G₀及其相应的进气速度V和进气角η，转子运行转速Ω作为优化设计变量；

目标函数2：f₂(x)：转子系统运行稳定性极限阈值区间最优，由步骤(1-2)中稳定性极限阈值求解方法可知，对于有围带起主要减振增稳作用的大型汽轮机组，可通过优化围带阻尼进而优化转子系统稳定性极限阈值-特征乘子μ_F，故取优化目标函数2为各稳定运行工况下的围带阻尼最优，即：

f₂(x)＝max[C_s ^s]＝max[f(μ^s,K^s _h,β^s,e₀ ^s)] (15)

其中，上标s表示稳定的运行工况；C_s为该运行工况下转子系统围带阻尼；μ、K_h、β和e₀依次分别为相邻围带接触面的摩擦系数、接触刚度系数、接触角和初始间隙。

确定f₂(x)的优化设计变量：在设计阶段，转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数初步设计确定的前提下，目标函数2的优化设计变量主要为围带参数：μ，K_h，β和e₀。

进一步地，在所述步骤(3)中：以机组整个运行周期各运行工况的供电/供热功率和效率曲线、各部件强度/刚度、疲劳寿命、无共振、强迫振动响应小于设定振幅等为约束条件如下：

a.依据所设计的机组在各运行工况的供电/供热功率和效率曲线进行转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数的初步设计确定。并在优化设计过程中，各优化设计变量及其摄动修改增量的选取均应满足机组的功率和效率曲线；

b.为避免强度失效，各部件在不同运行工况载荷下的各类应力/应变不超过材料设定值,即满足：σ_imax(x,y,z)≥[σ(x,y,z)]；

c.无共振发生：运行转速Ω尽量避开0.75-1.25倍的各阶临界转速Ω_r(r＝1,2,3,.....p,p为正整数)区间，即满足：Ω≤0.75Ω_r或Ω≥1.25Ω_r；

d.为避免动静件碰摩故障，转子系统各节点的各向最大振动位移响应q_imax(x,y,z)不超过设计设定振幅值(参照相应的设计手册),即满足：q_imax(x,y,z)≤[q(x,y,z)]；

e.为适应实际工况，转子系统的转速Ω小于等于常规运行转速的1.2倍；系统的发电/热功率P大于等于额定发电/热功率。

进一步地，在所述步骤(4)中，采用微量摄动修改法对转子系统稳定性进行有约束的多目标动态优化，具体步骤如下：

步骤(4-1)：确定满足目标函数2的围带参数：给定第i组运行工况参数(转速Ω_i，实际进气流量G_i(V_i,η_i))，并根据步骤(1-2)中的围带阻尼求解方法求解确定该运行工况下围带阻尼最优的第j组围带参数(摩擦系数μ_j,接触刚度K_hj,接触角β_j,初始间隙e_0j)，满足目标函数2；

步骤(4-2)：利用步骤(1)中的稳定性判定准则判定第i组运行工况参数(Ω_i，G_i(V_i,η_i))下具有第j组围带参数(μ_j，K_hj，β_j，e_0j)的转子系统的稳定性，若系统失稳，则进入步骤(4-2-1)；若稳定，则进入步骤(4-2-2)；

步骤(4-2-1)：在不同运行工况的三维坐标系(Ω,V,η)下标记当前工况点(Ω_i,V_i,η_i)为具有第j组围带参数的转子系统在该工况下的失稳区；

步骤(4-2-2)：判断第i组运行工况(Ω_i，G_i(V_i,η_i))下的供电/热效率L^s _H是否满足目标函数1(式(13)-(14))，若不满足，则辨识确定敏感优化设计变量，在满足步骤(3)约束条件下选取一组微摄动变量(ΔΩ,ΔV,Δη)，摄动修改第i组运行工况参数(其中，1)纯供电机组，微摄动变量ΔV＝0和Δη＝0，ΔΩ≠0；2)热电联产机组供热季:微摄动变量ΔV、Δη、ΔΩ≠0)为第i+1组运行工况(V_i+1＝V_i+ΔV，η_i+1＝η_i+Δη，Ω_i+1＝Ω_i+ΔΩ)，并返回步骤(4-1)计算在第i+1组运行工况下围带阻尼最优的第j+1组围带参数(μ_j+1,K_hj+1,β_j+1,e_0j+1)和步骤(4-2)判定第i+1组运行工况下是否稳定；若不稳定，进入步骤(4-2-1)；若稳定，则进行步骤(4-2-3)；

步骤(4-2-3)：第j+1组围带参数(μ_j+1，K_hj+1，β_j+1，e_0j+1)即为第i+1组运行工况下满足两个目标函数的最优围带参数；记录第j+1组围带参数以及对应的最优供电/热效率和稳定性极限阈值；

步骤(4-3)：通过上述步骤扫描所有运行工况，获得动态优化后的转子系统在各运行工况的稳定性极限阈值，并将其组合绘制可得到转子系统的三维稳定性极限阈值曲线。

结合具体实例对本公开进行进一步说明：

以某1000MW级大型汽轮机组为例，该汽轮机组是无备用低压转子的热电联产机组，换热季需抽取低压缸部分进气G_H用以供热以满足热电联产需求，低压缸的实际进气流量G＝G₀-G_H，即低压缸的相对进气流量G/G₀≤1，为小流量工况，由于供热季实际用于供热的流量G_H决定机组的供热效率，故对于目标函数1，低压缸的进气流量G越小越优化。依据该机组供电/热功率和效率曲线初步设计的低压末级转子结构如图2所示，该转子采用平行四边形围带，初始围带参数设置为：摩擦系数μ＝0.05,接触刚度系数K_h＝0.2,接触角β＝15°,初始间隙e₀＝1mm；常规转速Ω＝3000rpm。采用本公开提出的基于围带阻尼的汽轮机转子系统稳定性动态优化方法对该汽轮机低压末级转子系统进行动态优化设计。

步骤(1)获取汽轮机转子系统的运行工况参数；建立转子系统稳定性判定准则和稳定性极限阈值的求解方法：

步骤(1-1)获取汽轮机转子系统的运行工况参数；选择据特征乘子与单位圆关系判定转子系统稳定性的准则。

步骤(1-2)：建立稳定性极限阈值的求解方法：1)依据该机组供电/热功率和效率曲线初步设计的转子系统结构建立转子系统三维几何模型如图2，并进一步离散为三维实体单元，设置边界条件，求解转子系统模态质量m_r和刚度k_r；再对整个转子系统施加转速(insert-Rotational Velocity)，求解(solution)得到忽略围带阻尼的固有频率ω_nr；2)进一步采用基于流场仿真的流体激励求解方法求解转子系统多频激励F(t)；3)根据1)、2)利用围带阻尼定量求解方法，求解多频激励F(t)下的围带阻尼C_s；4)根据1)、2)、3)结果求解频响函数获取步骤(1-1)中判定准则所需要的特征乘子。

步骤(2)确定多目标动态优化的目标函数和优化设计变量：

目标函数1：f₁(x)：汽轮机组的供电/热总效率达到最大值，即

f₁(x)＝max[L_H ^s]＝max[BG_H ^s]

L_H＝BG_H

其中，G_H＝G₀-G，换热季低压缸的相对进气流量G/G₀≤1，为小流量工况，实际用于供热的流量G_H决定机组的供热效率，满足优化目标1的条件是在保证稳定运行前提下低压缸的流量G越小即机组供电/热功率越大；初设产能效率系数B＝90％。

确定f₁(x)的优化设计变量：初步设计阶段，在已知功率、效率曲线等约束条件下，选取整个运行周期各运行工况(包括升速启动、甩负荷停车、不同转速及不同流量等)参数：相对进气流量G/G₀及相应的进气速度V、进气角η和转子运行转速Ω作为优化设计变量；作为一种实施方式之一，各变量及取值范围采用多因素多水平正交法结合响应曲面法确定。

目标函数2：f₂(x)：转子系统运行稳定性极限阈值区间最优：

f₂(x)＝max[C_s]＝max[f(μ^s,K^s _h,β^s,e₀ ^s)]

其中，上标s表示稳定的运行工况；C_s为该运行工况下转子系统围带阻尼。确定f₂(x)的优化设计变量：在设计阶段，转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数初步设计确定的前提下，目标函数2的优化设计变量主要为围带参数：相邻围带接触面的摩擦系数μ，接触刚度系数K_h，接触角β、e₀、μ、K_h、β和初始间隙e₀。各变量及取值范围采用多因素多水平正交法结合响应曲面法确定。

步骤(3)：设定约束条件如下：

a.依据所设计的机组在整个运行周期各运行工况的供电/供热功率和效率曲线进行转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数的初步设计确定。并在优化设计过程中，各优化设计变量及其摄动修改增量的选取均应满足机组的功率和效率曲线；

b.为避免强度失效，各部件在不同运行工况载荷下的各类应力/应变不超过材料设定值,即满足：σ_imax(x,y,z)≥[σ(x,y,z)]；

c.无共振发生：运行转速Ω尽量避开0.75-1.25倍的各阶临界转速Ω_r(r＝1,2,3,.....)区间，即满足：Ω≤0.75Ω_r或Ω≥1.25Ω_r；

d.为避免动静件碰摩故障，转子系统各节点的各向最大振动位移响应q_imax(x,y,z)不超过设计设定振幅值(参照相应的设计手册),即满足：q_imax(x,y,z)≤[q(x,y,z)]；

e.为适应实际工况，转子系统的转速Ω小于等于常规运行转速的1.2倍；系统的发电/热功率P大于等于额定发电/热功率。

据上述约束条件，初步确定各优化设计变量的优化取值范围为：

运行工况参数：Ω优化取值范围：0＜Ω≤4000rpm；额定工况转速为Ω＝3000rpm；G/G₀优化取值范围:0.15＜G/G₀≤1,其中相应的V优化取值范围：45m/s＜V≤300m/s，η优化取值范围：90°＜η≤160°。围带参数：μ优化取值范围:0＜μ≤0.2，K_h优化取值范围:0＜K_h≤1，β优化取值范围:0°≤β≤20°，e₀优化取值范围:0＜e₀≤1.2。

步骤(4)，采用微量摄动修改法对转子系统稳定性进行有约束的多目标动态优化，具体步骤如下：

步骤(4-1)设定初始计算的(i＝1)运行工况为该转子系统的额定运行工况(G/G₀＝1(V_i＝300m/s,η_i＝90°)，Ω＝3000rpm)，据步骤(1-2)的围带阻尼计算方法从初步设计的(j＝1)围带参数((μ₁,K_h1,β₁,e₀₁))开始计算在(i＝1)运行工况下的各组围带参数的围带阻尼，确定该(i＝1)运行工况下围带阻尼最优的第j组围带参数(μ_j,K_hj,β_j,e_0j)，并记录对应的围带阻尼最优值

步骤(4-2)据步骤(1-2)稳定性极限阈值计算方法，在有限元分析软件(ANSYS)的谐响应分析模块中求解第(i＝1)运行工况下的具有第j组围带参数(μ_j,K_hj,β_j,e_0j)的转子系统的频响函数，读取各阶特征根，并计算该工况下的系统特征乘子

据Floquet理论判定该运行工况下第j组围带参数转子系统的稳定性；

若系统失稳，进入步骤(4-2-1)，在不同运行工况的三维坐标系(Ω,V,η)下标记当前(i＝1)工况点(Ω_i,G_i(V_i,η_i))为具有第j组围带参数的转子系统在该工况下的失稳区；

若稳定，则进入步骤(4-2-2)，判断该(i＝1)运行工况(Ω_i,G_i(V_i,η_i))下的供电/热效率L^s _H是否满足目标函数1(式(13)-(14))，若不满足，则辨识确定敏感优化设计变量，热电联产机组供热季:各微摄动变量ΔV、Δη、ΔΩ≠0,且按相对进气量(G/G₀)逐渐减小的原则，在满足步骤(3)约束条件下选取一组微摄动变量(ΔΩ，ΔG(ΔV，Δη))，将当前(i＝1)运行工况参数摄动修改为第i+1组运行工况(Ω_i+1＝Ω_i+ΔΩ，G_i+1＝(G_i+ΔG)(V_i+1＝V_i+ΔV，η_i+1＝η_i+Δη))，并返回步骤(4-1)计算在第i+1组运行工况下围带阻尼最优的第j+1组围带参数(μ_j+1,K_hj+1,β_j+1,e_0j+1)和步骤(4-2)判定第i+1组运行工况下是否稳定；若不稳定，进入步骤(4-2-1)；

若稳定，则进行步骤(4-2-3)：第j+1组围带参数(μ_j+1，K_hj+1，β_j+1，e_0j+1)即为第i+1组运行工况下满足两个目标函数的最优围带参数；记录第j+1组围带参数以及对应的最优供电/热效率和稳定性极限阈值；

(4-3)通过上述步骤扫描所有运行工况，获得动态优化后的转子系统在各运行工况的稳定性极限阈值，得到本算例下的最优围带参数组合为(μ＝0.23,K_h＝0.62,β＝8.54°,e₀＝0.35)，通过利用本公开的稳定性动态优化方法，得到优化后的转子系统三维(Ω-η-V)稳定性极限曲面。

综上，基于本公开方法对所述实施例转子系统进行动态优化，得到上述实例下的最优围带参数及其优化后的转子系统三维(Ω-η-V)稳定性极限曲面(图3)，图3中三维Ω-η-V极限曲面以下是稳定区，曲面以上是不稳定区，对比优化前的转子系统三维(Ω-η-V)稳定性极限曲面(图4)可以看出:优化后的转子系统在不同转速、不同进气工况下的增稳能力较优化前有大幅度提升，转子系统稳定运行区间显著增大。据优化结果进一步分析不同进气角η下的转速-进气速度稳定性极限(如图5(a)-图5(c)所示)，对于小进气角工况(如η＝90°)，系统在额定转速(3000rpm)下的进气速度V的稳定性极限阈值由优化前的330m/s(图5(a)中虚线)升高到450m/s(图5(a)中实线)，稳定运行区间增大了40.0％，对于进气角η较大的工况，如η＝125°和160°，优化后系统在额定转速(3000rpm)下的进气速度V的稳定性极限阈值分别增大了56.6％和76.4％，本公开提出的基于围带阻尼的稳定性动态优化方法可得到较好的优化效果。图6为本公开的系统最大位移响应曲线。

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

Claims (6)

1.一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法，其特征是，包括以下步骤：

步骤(1)：获取汽轮机转子系统的运行工况参数；建立转子系统稳定性判定准则和稳定性极限阈值的求解方法；具体步骤为：

步骤(1-1)：建立转子系统稳定性判定准则：基于Floquet理论，系统方程单值矩阵的特征乘子存在等于1的解，系统处于临界稳定状态，同时若其他特征乘子|μ_F|均小于1，则系统失稳；汽轮机组转子系统稳定性可依据特征乘子与单位圆的关系判定；

步骤(1-2)：建立稳定性极限阈值的求解方法，具体步骤为：

步骤(1-2-1)：基于Floquet理论，特征乘子μ_F可表示为：

μ_F＝exp(ΛT) (1)

其中，T＝Im(lnμ_F)/ω为主振动周期；Λ为系统特征方程的特征根；

将转子系统运动微分方程写为复坐标形式：

-m_rω_n ²+c_riω_n+k_r＝F₀exp(iωt) (2)

其中，ω为激励频率；ω_n为系统固有频率；i²＝-1；则转子系统第r阶频响函数H(ω)_r写为以下实部-虚部形式：

其中，定义频率比λ_r＝ω/ω_nr为系统所受激励频率与系统第r阶固有频率之比；阻尼比ξ_r＝C_r/2m_rω_nr为系统阻尼C_r与第r阶临界阻尼之比，阻尼C_r＝C_l+C_s，对于有围带的转子系统由于系统的结构阻尼C_l远小于围带摩擦阻尼C_s，为方便求解，忽略结构阻尼C_l，则有C_r＝C_s；k_r为第r阶模态刚度；

转子系统运动微分方程(2)的特征方程为：

det{I-Λ[H(ω)]}＝0 (4)

则系统特征方程(4)的特征根Λ可写为

Λ＝H(ω)^-1 (5)

由频响函数中所含各变量可知，对于具有几何非线性柔性叶片和激励依赖非线性围带阻尼的大型汽轮机转子系统，在不同运行工况下，其频率比

阻尼比

呈现不同的非线性特性；依据特征乘子与单位圆的关系判定转子系统稳定性,获得转子系统特征方程的特征根Λ从而求解特征乘子|μ_F|，由式(5)知特征根Λ为系统频响函数的逆，则系统稳定性的判定可归结为求解系统在不同运行工况下的频响函数问题；

据公式(3)，频响函数涉及系统各阶固有频率ω_nr、模态质量m_r、刚度k_r和阻尼C_r，其中阻尼C_r＝C_l+C_s，对于有围带的转子系统由于系统的结构阻尼C_l远小于围带阻尼C_s，为方便求解，忽略结构阻尼C_l，则有C_r＝C_s；分别求解各项以获得频响函数；

步骤(1-2-2)：计算各阶固有频率ω_nr、模态质量m_r和刚度k_r；

步骤(1-2-3)：计算获得围带阻尼C_s；

步骤(1-2-4)：利用有限元分析软件求解转子系统频响函数并获取特征乘子：

据步骤(1-2-1)～(1-2-3)获得的转子系统各阶固有频率ω_nr、模态质量m_r、刚度k_r和围带阻尼C_s在FEMS的谐响应分析模块中求解获得转子系统的实部-虚部形式的频响函数H(ω)_r，并读取各阶特征根Λ_r，r为正整数，将其代入式(1)即获得转子系统各阶特征乘子|μ_F|_r；其中，步骤(1-2-3)的具体步骤包括：

据多频激励下大型汽轮机转子围带摩擦阻尼等效模型和基于瞬态响应包络线的围带阻尼定量求解方法，多频激励F(t)下的围带阻尼C_s为：

其中，h＝1,2,3,......,p；p为正整数；

为多频激励力F(t)中第h阶谐波激励分量作用下的围带阻尼：

其中，m_r是转子系统模态质量；ω_nr为忽略围带摩擦阻尼的系统固有频率，

和

分别为多频激励F(t)的第h阶谐波激励分量作用下转子系统最大振动位移响应节点的相邻位移响应峰值振幅，δ^h为对数减幅比；

步骤(1-2-3-1)：计算获得系统激励F(t)：运行在湿蒸汽流体环境的大型汽轮机转子系统，所受的激励力F(t)由两部分组成：

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

其中，F_m(t)是相应运行转速下由动不平衡效应引起的机械激励；F_f(t)是湿蒸汽流体介质在叶片和围带表面的瞬态分布压力及脉动流体多频激励；

步骤(1-2-3-2)：获得系统多频激励F(t)下转子系统的最大振动位移响应，确定对数减幅比；

步骤(2)：确定多目标动态优化的目标函数及其优化变量：

目标函数1，f₁(x)：稳定运行下汽轮机组供电/供热效率最大，涉及的优化变量为运行工况参数：转子运行转速Ω，相对进气流量G/G₀、进气速度V和进气角η；其中，G为实际进气流量，G₀为额定进气流量；

目标函数2，f₂(x)：转子系统运行稳定性极限阈值区间最优，即转子系统的围带阻尼最优，涉及的优化变量为围带参数：相邻围带间的摩擦系数μ、接触角β、接触刚度K_h和初始间隙e₀；

步骤(3)：设定约束条件：以供电/供热功率、效率曲线、各部件强度、疲劳寿命、无共振、强迫振动响应小于设定振幅、以及只供电或热电联产实际运行工况为约束条件；

步骤(4)：对转子系统稳定性进行有约束的多目标动态优化，具体步骤如下：

步骤(4-1)：确定满足目标函数2的围带参数：给定第i组运行工况参数，并根据步骤(1-2)中的围带阻尼求解方法求解确定该运行工况下围带阻尼最优的第j组围带参数，满足目标函数2；

步骤(4-2)：利用步骤(1)中的稳定性判定准则判定第i组运行工况参数下具有第j组围带参数的转子系统的稳定性，若系统失稳，则进入步骤(4-2-1)；若稳定，则进入步骤(4-2-2)；

步骤(4-2-1)：在不同运行工况的三维坐标系(Ω,V,η)下标记当前工况点为具有第j组围带参数的转子系统在该工况下的失稳区；

步骤(4-2-2)：判断第i组运行工况下的供电/热效率L^s _H是否满足目标函数1，若不满足，则辨识确定敏感优化设计变量，在满足步骤(3)约束条件下选取一组微摄动变量，摄动修改第i组运行工况参数，并返回步骤(4-1)计算在第i+1组运行工况下围带阻尼最优的第j+1组围带参数和步骤(4-2)判定第i+1组运行工况下是否稳定；若不稳定，进入步骤(4-2-1)；若稳定，则进行步骤(4-2-3)；

步骤(4-2-3)：第j+1组围带参数即为第i+1组运行工况下满足两个目标函数的最优围带参数；记录第j+1组围带参数以及对应的最优供电/热效率和稳定性极限阈值；

步骤(4-3)：通过步骤(4-1)～步骤(4-2)扫描所有运行工况，获得动态优化后的转子系统在各运行工况的稳定性极限阈值，并将其组合绘制得到转子系统的三维稳定性极限阈值曲线。

2.如权利要求1所述的方法，其特征是，步骤(1-2-2)具体步骤包括：

步骤(1-2-2-1)：建立几何模型并离散划分网格：根据设计的汽轮机组转轴的长度和直径、轴承支承跨度、叶片参数以及围带参数建立转子系统三维几何模型，并将其导入有限元分析软件，在模态分析模块Modal中的网格离散模块mesh中将转子系统离散为三维实体单元；

步骤(1-2-2-2)：设置边界条件：在FEMS的模态分析Modal模块中的支承参数设置选项Fix support中利用滑动轴承参数计算软件XROTOR计算获得支承阻尼C_O和支承刚度K_O，在轴承支承处添加弹簧接触connection-spring，设置C_O和K_O；

步骤(1-2-2-3)：在FEMS的模态分析模块Modal中求解转子系统模态质量m_r和刚度k_r；再对整个转子系统施加转速insert-Rotational Velocity，求解solution得到忽略围带阻尼的固有频率ω_nr。

3.如权利要求1所述的方法，其特征是，步骤(1-2-3-1)的具体步骤包括：

步骤(1-2-3-1-1)：机械激励F_m(t)可通过在FEMS的瞬态响应求解模块转速激励输入选项中选择运行转速作为F_m(t)激励输入；

步骤(1-2-3-1-2)：流体激励F_f(t)，利用基于流场仿真的流体激励求解方法得到如式(10)和(11)所示的流体激励F_f(t)：

F_f(t)＝N^Tdf(t)^T (10)

其中，p_i(x,y,z,t)和d_i分别为流场计算域中任意节点i上在不同运行工况下的瞬态压力及其距叶片结构节点j的距离；n为离散点的数量；N为叶片六面体单元形函数；

步骤(1-2-3-1-3)：利用基于FFT的多谐波平衡法确定该流体激励的前h阶频率成分以及各阶谐波频率下的谐波激励分量，将如式(10)和(11)表示的F_f(t)转化为一个常量和若干个不同频率成分谐波激励分量的叠加，如式(12)所示，可在FEMS的瞬态响应求解Transient模块中作为激励函数输入；

4.如权利要求1所述的方法，其特征是，步骤(1-2-3-2)的具体步骤包括：

步骤(1-2-3-2-1)：设置边界条件：在FEMS的瞬态响应求解模块中的模型Model选项，在相邻围带接触面间添加摩擦接触，设置摩擦系数μ和接触刚度K_h；

步骤(1-2-3-2-2)：施加激励并求解：在FEMS的瞬态响应求解模块施加F(t)，具体的，F_m(t)在转速激励输入选项inertial中选择运行转速作为激励输入；F_f(t)在载荷激励输入选项insert-Force中以

的形式施加由多阶谐波激励分量组成的多频激励力；在求解solution选项中求解第h阶谐波激励分量下的转子系统振动位移响应；重复上述步骤，依次分别求解出各阶谐波激励分量下的转子系统振动位移响应；

步骤(1-2-3-2-3)：确定第h阶谐波激励分量下的最大振动位移响应节点，输出该节点的最大振动位移响应曲线，读取相邻振动位移响应峰值的振幅

和

并依次代入式(10)求解出各阶振动位移响应对数减幅比δ^h；

步骤(1-2-3-2-4)：将步骤(1-2-2)所获转子系统质量m_r和忽略围带阻尼的固有频率ω_nr及步骤(1-2-3-2-3)所获的各阶对数减幅比δ^h代入式(7)计算围带阻尼

并将其代入式(6)求解得到多频激励F(t)下的围带阻尼C_s(F,t)。

5.如权利要求1所述的方法，其特征是，在所述步骤(2)中，确定多目标动态优化的目标函数和优化设计变量：

目标函数1：f₁(x)：汽轮机组的供电/热总效率达到最大值，即

f₁(x)＝max[L_H ^s]＝max[BG_H ^s](13)

其中，上标s表示当前运行工况是稳定的；L_H为该运行工况下输出的供电/热总效率，表示为

L_H＝BG_H (14)

其中，G_H为实际用于供电/热的进气流量，单位：kg/h，按汽轮机组的两种不同用途划分为两大类运行工况：

1)常规供电：低压缸的相对进气流量G/G₀＝1，低压缸的相对进气流量为额定进气流量工况，G_H＝0；其中，G为实际进气流量、G₀为额定进气流量；

2)热电联产：对于无备用低压缸的抽气式热电联产机组，供热季低压缸的相对进气流量G/G₀≤1，供热季实际用于供热的流量G_H＝G₀–G；B为单位进气流量G_H下考虑运行过程中能量转换消耗及效率因素的产能效率系数；

确定f₁(x)的优化设计变量：优化设计变量是可控可调的初始变量，在设计阶段已知功率、效率曲线约束条件下，选取整个运行周期不同运行工况的工况参数：相对进气流量G/G₀及其相应的进气速度V和进气角η，转子运行转速Ω作为优化设计变量；

目标函数2：f₂(x)：转子系统运行稳定性极限阈值区间最优，由步骤(1-2)中稳定性极限阈值求解方法可知，对于有围带起主要减振增稳作用的大型汽轮机组，可通过优化围带阻尼进而优化转子系统稳定性极限阈值-特征乘子μ_F，故取优化目标函数2为各稳定运行工况下的围带阻尼最优，即：

f₂(x)＝max[C_s ^s]＝max[f(μ^s,K^s _h,β^s,e₀ ^s)] (15)

其中，上标s表示稳定的运行工况；C_s为该运行工况下转子系统围带阻尼；μ、K_h、β和e₀依次分别为相邻围带接触面的摩擦系数、接触刚度系数、接触角和初始间隙；

确定f₂(x)的优化设计变量：在设计阶段，转子系统的结构及相关参数初步设计确定的前提下，目标函数2的优化设计变量为围带参数：μ，K_h，β和e₀。

6.如权利要求1所述的方法，其特征是，在所述步骤(3)中：以机组整个运行周期各运行工况的供电/供热功率和效率曲线、各部件强度/刚度、疲劳寿命、无共振、强迫振动响应小于设定振幅为约束条件如下：

a.依据所设计的机组在各运行工况的供电/供热功率和效率曲线进行转子系统的结构及相关参数的初步设计确定；并在优化设计过程中，各优化设计变量及其摄动修改增量的选取均应满足机组的功率和效率曲线；

b.为避免强度失效，各部件在不同运行工况载荷下的各类应力/应变不超过材料设定值,即满足：σ_imax(x,y,z)≥[σ(x,y,z)]；

c.无共振发生：运行转速Ω尽量避开0.75-1.25倍的各阶临界转速Ω_r区间，即满足：Ω≤0.75Ω_r或Ω≥1.25Ω_r；

d.为避免动静件碰摩故障，转子系统各节点的各向最大振动位移响应q_imax(x,y,z)不超过设计设定振幅值,即满足：q_imax(x,y,z)≤[q(x,y,z)]；

e.为适应实际工况，转子系统的转速Ω小于等于常规运行转速的1.2倍；系统的发电/热功率P大于等于额定发电/热功率。

Images