CN111062177B - 一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法 - Google Patents
一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法 Download PDFInfo
- Publication number
- CN111062177B CN111062177B CN201911326462.1A CN201911326462A CN111062177B CN 111062177 B CN111062177 B CN 111062177B CN 201911326462 A CN201911326462 A CN 201911326462A CN 111062177 B CN111062177 B CN 111062177B
- Authority
- CN
- China
- Prior art keywords
- shroud
- damping
- rotor system
- excitation
- under
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000013016 damping Methods 0.000 title claims abstract description 140
- 238000000034 method Methods 0.000 title claims abstract description 89
- 238000013461 design Methods 0.000 claims abstract description 77
- 238000005457 optimization Methods 0.000 claims abstract description 72
- 230000004044 response Effects 0.000 claims abstract description 63
- 238000005316 response function Methods 0.000 claims abstract description 31
- 230000004048 modification Effects 0.000 claims abstract description 13
- 238000012986 modification Methods 0.000 claims abstract description 13
- 230000008569 process Effects 0.000 claims abstract description 10
- 238000004364 calculation method Methods 0.000 claims abstract description 8
- 230000005284 excitation Effects 0.000 claims description 111
- 238000006073 displacement reaction Methods 0.000 claims description 34
- 230000001052 transient effect Effects 0.000 claims description 27
- 239000012530 fluid Substances 0.000 claims description 24
- 238000004458 analytical method Methods 0.000 claims description 14
- 230000009467 reduction Effects 0.000 claims description 12
- 238000010248 power generation Methods 0.000 claims description 8
- 230000000694 effects Effects 0.000 claims description 7
- 230000009471 action Effects 0.000 claims description 6
- 238000006243 chemical reaction Methods 0.000 claims description 6
- XXXSILNSXNPGKG-ZHACJKMWSA-N Crotoxyphos Chemical compound COP(=O)(OC)O\C(C)=C\C(=O)OC(C)C1=CC=CC=C1 XXXSILNSXNPGKG-ZHACJKMWSA-N 0.000 claims description 4
- 230000001687 destabilization Effects 0.000 claims description 4
- 239000000463 material Substances 0.000 claims description 4
- 239000005364 simax Substances 0.000 claims description 4
- 238000004088 simulation Methods 0.000 claims description 4
- 230000003068 static effect Effects 0.000 claims description 4
- 230000001419 dependent effect Effects 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 238000000605 extraction Methods 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 230000003094 perturbing effect Effects 0.000 claims description 3
- 230000001629 suppression Effects 0.000 abstract 1
- 238000002715 modification method Methods 0.000 description 6
- 230000008878 coupling Effects 0.000 description 5
- 238000010168 coupling process Methods 0.000 description 5
- 238000005859 coupling reaction Methods 0.000 description 5
- 230000001133 acceleration Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 230000007547 defect Effects 0.000 description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000002068 genetic effect Effects 0.000 description 2
- 238000009434 installation Methods 0.000 description 2
- 230000010349 pulsation Effects 0.000 description 2
- 206010010356 Congenital anomaly Diseases 0.000 description 1
- 230000003416 augmentation Effects 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000013078 crystal Substances 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000010206 sensitivity analysis Methods 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 230000005514 two-phase flow Effects 0.000 description 1
- 239000009891 weiqi Substances 0.000 description 1
Images
Landscapes
- Turbine Rotor Nozzle Sealing (AREA)
- Control Of Turbines (AREA)
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,f1(x):稳定运行下汽轮机组供电/供热效率最大,涉及的优化变量为运行工况参数:转子运行转速Ω,相对进气流量G/G0、进气速度V和进气角η;其中,G为实际进气流量,G0为额定进气流量;
目标函数2,f2(x):转子系统运行稳定性极限阈值区间最优,即转子系统的围带阻尼最优,涉及的优化变量为围带参数:相邻围带间的摩擦系数μ、接触角β、接触刚度Kh和初始间隙e0;
步骤(3):设定约束条件:以供电/供热功率、效率曲线、各部件强度、疲劳寿命、无共振、强迫振动响应小于设定振幅、以及只供电或热电联产实际运行工况为约束条件;
步骤(4):对转子系统稳定性进行有约束的多目标动态优化。
进一步地,在所述步骤(1)中,获取汽轮机转子系统的运行工况参数;具体包括:转子运行转速Ω、实际进气流量G、进气速度V、进气角η和额定进气流量G0。
进一步地,在所述步骤(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)/ω为主振动周期;Λ为系统特征方程的特征根。
将转子系统运动微分方程写为复坐标形式:
-mrωn 2+criωn+kr=F0exp(iωt) (2)
其中,ω为激励频率;ωn为系统固有频率;i2=-1;则转子系统第r阶频响函数H(ω)r写为以下实部-虚部形式:
其中,定义频率比λr=ω/ωnr为系统所受激励频率与系统第r阶固有频率之比;阻尼比为系统阻尼Cr与第r阶临界阻尼(Crb=2mrωnr)之比,阻尼Cr=Cl+Cs,对于有围带的转子系统由于系统的结构阻尼Cl远小于围带摩擦阻尼Cs,为方便求解,忽略结构阻尼Cl,则有Cr=Cs;kr为第r阶模态刚度。
转子系统运动微分方程(2)的特征方程为:
det{I-Λ[H(ω)]}=0 (4)
则系统特征方程(4)的特征根Λ可写为
Λ=H(ω)-1 (5)
由频响函数(式(3))中所含各变量可知,对于具有几何非线性柔性叶片和激励依赖非线性围带阻尼的大型汽轮机转子系统,在不同运行工况下,其频率比阻尼比(各部件阻尼等)呈现不同的非线性特性。依据特征乘子与单位圆的关系判定转子系统稳定性,获得转子系统特征方程(式(4))的特征根Λ从而求解特征乘子|μF|(式(1)),由式(5)知特征根Λ为系统频响函数(式(3))的逆,则系统稳定性的判定可归结为求解系统在不同运行工况下的频响函数问题。
进一步地,据公式(3),频响函数涉及系统各阶固有频率ωnr、模态质量mr、刚度kr和阻尼Cr,其中阻尼Cr=Cl+Cs,对于有围带的转子系统由于系统的结构阻尼Cl远小于围带阻尼Cs,为方便求解,忽略结构阻尼Cl,则有Cr=Cs;分别求解各项以获得频响函数。
步骤(1-2-2):计算各阶固有频率ωnr、模态质量mr和刚度kr:
步骤(1-2-2-1):建立几何模型并离散划分网格:根据设计的汽轮机组转轴的长度和直径、轴承支承跨度、叶片参数(安装角、扭转角、展弦比)以及围带参数(接触角β、接触间隙e0)等建立转子系统三维几何模型,并将其导入有限元分析软件(以下简称FEMS),在模态分析模块(Modal)中的网格离散模块(mesh)中将转子系统离散为三维实体单元;
步骤(1-2-2-2):设置边界条件:在FEMS的模态分析(Modal)模块中的支承参数设置选项(Fix support)中利用滑动轴承参数计算软件(XROTOR)计算获得支承阻尼CO和支承刚度KO,在轴承支承处添加弹簧接触(connection-spring),设置CO和KO;
步骤(1-2-2-3):在FEMS的模态分析模块(Modal)中求解转子系统模态质量mr和刚度kr;再对整个转子系统施加转速(insert-Rotational Velocity),求解(solution)得到忽略围带阻尼的固有频率ωnr;
步骤(1-2-3):计算获得围带阻尼Cs;
据多频激励下大型汽轮机转子围带摩擦阻尼等效模型和基于瞬态响应包络线的围带阻尼定量求解方法,多频激励F(t)下的围带阻尼Cs为:
步骤(1-2-3-1):计算获得系统激励F(t):运行在湿蒸汽等流体环境的大型汽轮机转子系统,所受的激励力F(t)由两部分组成:
F(t)=Ff(t)+Fm(t) (9)
其中,Fm(t)是相应运行转速下由动不平衡效应等引起的机械激励;Ff(t)是湿蒸汽等流体介质在叶片和围带表面的瞬态分布压力及脉动等流体多频激励;
步骤(1-2-3-1-1):机械激励Fm(t)通过在FEMS的瞬态响应求解模块转速激励输入选项中选择运行转速作为Fm(t)激励输入;
步骤(1-2-3-1-2):流体激励Ff(t),利用基于流场仿真的流体激励求解方法得到如式(10)和(11)所示的流体激励Ff(t):
Ff(t)=NTdf(t)T (10)
其中,pi(x,y,z,t)和di分别为流场计算域中任意节点i(离散点)上在不同运行工况(进气流量G、进气速度V、进气角η和转速Ω)下的瞬态压力及其距叶片结构节点j(插值点)的距离;n为离散点的数量;N为叶片六面体单元形函数。
步骤(1-2-3-1-3):利用基于FFT的多谐波平衡法确定该流体激励的前h阶频率成分(可取h≥5)以及各阶谐波频率下的谐波激励分量,将如式(10)和(11)表示的Ff(t)转化为一个常量和若干个不同频率成分谐波激励分量的叠加,如式(12)所示,可在FEMS的瞬态响应求解(Transient)模块中作为激励函数输入。
步骤(1-2-3-2):获得系统多频激励F(t)下转子系统的最大振动位移响应,确定对数减幅比,具体步骤如下:
步骤(1-2-3-2-1):设置边界条件:在FEMS的瞬态响应求解模块中的模型(Model)选项,在相邻围带接触面间添加摩擦接触(insert-connection-fraction),设置摩擦系数μ和接触刚度Kh;
步骤(1-2-3-2-2):施加激励并求解:在FEMS的瞬态响应求解模块施加F(t),具体的,Fm(t)在转速激励输入选项(inertial)中选择运行转速作为激励输入;Ff(t)在载荷激励输入选项(insert-Force)中以的形式施加由多阶谐波激励分量组成的多频激励力;在求解(solution)选项中求解第h阶谐波激励分量下的转子系统振动位移响应;重复上述步骤,依次分别求解出各阶谐波激励分量下的转子系统振动位移响应。
步骤(1-2-3-2-3):确定第h阶谐波激励分量下的最大振动位移响应节点,输出该节点的最大振动位移响应曲线,读取相邻振动位移响应峰值的振幅和并依次代入式(10)求解出各阶振动位移响应对数减幅比δh;
步骤(1-2-3-2-4):将步骤(1-2-2)所获转子系统质量mr和忽略围带阻尼的固有频率ωnr及步骤(1-2-3-2-3)所获的各阶对数减幅比δh代入式(7)计算围带阻尼并将其代入式(6)求解得到多频激励F(t)下的围带阻尼Cs(F,t)。
步骤(1-2-4):利用有限元分析软件求解转子系统频响函数并获取特征乘子:
据步骤(1-2-1)~(1-2-3)获得的转子系统各阶固有频率ωnr、模态质量mr、刚度kr和围带阻尼Cs在FEMS的谐响应分析模块中求解获得转子系统的实部-虚部形式的频响函数(H(ω)r),并读取各阶特征根Λr(r为正整数),将其代入式(1)即获得转子系统各阶特征乘子|μF|r。
进一步地,在所述步骤(2)中,确定多目标动态优化的目标函数和优化设计变量:
目标函数1:f1(x):汽轮机组的供电/热总效率达到最大值,即
f1(x)=max[LH s]=max[BGH s] (13)
其中,上标s表示当前运行工况是稳定的;LH为该运行工况下输出的供电/热总效率,可表示为
LH=BGH (14)
其中,GH为实际用于供电/热的进气流量(kg/h),按汽轮机组的两种不同用途划分为两大类运行工况:
1)常规供电:低压缸的相对进气流量G/G0=1(其中,G为实际进气流量、G0为额定进气流量)为额定进气流量工况,GH=0;
2)热电联产:对于无备用低压缸的抽气式热电联产机组,供热季低压缸的相对进气流量G/G0≤1(称小流量工况),供热季实际用于供热的流量GH=G0-G。B为单位进气流量GH下考虑运行过程中能量转换消耗及效率等因素的产能效率系数。
确定f1(x)的优化设计变量:优化设计变量是可控可调的初始变量,在设计阶段已知功率、效率曲线约束条件下,选取整个运行周期不同运行工况的工况参数:相对进气流量G/G0及其相应的进气速度V和进气角η,转子运行转速Ω作为优化设计变量;
目标函数2:f2(x):转子系统运行稳定性极限阈值区间最优,由步骤(1-2)中稳定性极限阈值求解方法可知,对于有围带起主要减振增稳作用的大型汽轮机组,可通过优化围带阻尼进而优化转子系统稳定性极限阈值-特征乘子μF,故取优化目标函数2为各稳定运行工况下的围带阻尼最优,即:
其中,上标s表示稳定的运行工况;Cs为该运行工况下转子系统围带阻尼;μ、Kh、β和e0依次分别为相邻围带接触面的摩擦系数、接触刚度系数、接触角和初始间隙。
确定f2(x)的优化设计变量:在设计阶段,转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数初步设计确定的前提下,目标函数2的优化设计变量为围带参数:μ,Kh,β和e0。
进一步地,在所述步骤(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.为避免动静件碰摩故障,转子系统各节点的各向最大振动位移响应qimax(x,y,z)不超过设计设定振幅值(参照相应的设计手册),即满足:qimax(x,y,z)≤[q(x,y,z)];
e.为适应实际工况,转子系统的转速Ω小于等于常规运行转速的1.2倍;系统的发电/热功率P大于等于额定发电/热功率。
进一步地,在所述步骤(4)中,对转子系统稳定性进行有约束的多目标动态优化,具体步骤如下:
步骤(4-1):确定满足目标函数2的围带参数:给定第i组运行工况参数(转速Ωi,实际进气流量Gi(Vi,ηi)),并根据步骤(1-2)中的围带阻尼求解方法求解确定该运行工况下围带阻尼最优的第j组围带参数(摩擦系数μj,接触刚度Khj,接触角βj,初始间隙e0j),满足目标函数2;
步骤(4-2):利用步骤(1)中的稳定性判定准则判定第i组运行工况参数(Ωi,Gi(Vi,ηi))下具有第j组围带参数(μj,Khj,βj,e0j)的转子系统的稳定性,若系统失稳,则进入步骤(4-2-1);若稳定,则进入步骤(4-2-2);
步骤(4-2-1):在不同运行工况的三维坐标系(Ω,V,η)下标记当前工况点(Ωi,Vi,ηi)为具有第j组围带参数的转子系统在该工况下的失稳区;
步骤(4-2-2):判断第i组运行工况(Ωi,Gi(Vi,ηi))下的供电/热效率Ls H是否满足目标函数1(公式(13)~公式(14)),若不满足,则辨识确定敏感优化设计变量,在满足步骤(3)约束条件下选取一组微摄动变量(ΔΩ,ΔV,Δη),摄动修改第i组运行工况参数(其中,1)纯供电机组,微摄动变量ΔV=0和Δη=0,ΔΩ≠0;2)热电联产机组供热季:微摄动变量ΔV、Δη、ΔΩ≠0)为第i+1组运行工况(Vi+1=Vi+ΔV,ηi+1=ηi+Δη,Ωi+1=Ωi+ΔΩ),并返回步骤(4-1)计算在第i+1组运行工况下围带阻尼最优的第j+1组围带参数(μj+1,Khj+1,βj+1,e0j+1)和步骤(4-2)判定第i+1组运行工况下是否稳定;若不稳定,进入步骤(4-2-1);若稳定,则进行步骤(4-2-3);
步骤(4-2-3):第j+1组围带参数(μj+1,Khj+1,βj+1,e0j+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,f1(x):稳定运行下汽轮机组供电/供热效率最大,涉及的优化变量主要为运行工况参数:转子运行转速Ω,相对进气流量G/G0(G为实际进气流量,G0为额定进气流量)及其相应的进气速度V和进气角η,等;
目标函数2,f2(x):转子系统运行稳定性极限阈值区间最优,即转子系统的阻尼(主要包括围带阻尼)最优,涉及的优化变量主要为围带参数:相邻围带间的摩擦系数μ、接触角β、接触刚度Kh和初始间隙e0,等;
步骤(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)/ω为主振动周期;Λ为系统特征方程的特征根。
将转子系统运动微分方程写为复坐标形式:
-mrωn 2+criωn+kr=F0 exp(iωt) (2)
其中,ω为激励频率;ωn为系统固有频率;i2=-1;则转子系统第r阶频响函数H(ω)r可写为以下实部-虚部形式:
其中,定义频率比λr=ω/ωnr为系统所受激励频率与系统第r阶固有频率之比;阻尼比为系统阻尼Cr与第r阶临界阻尼(Crb=2mrωnr)之比,其中,阻尼Cr=Cl+Cs,对于有围带的转子系统由于系统的结构阻尼Cl远小于围带摩擦阻尼Cs,为方便求解,忽略结构阻尼Cl,则有Cr=Cs;kr为第r阶模态刚度。
转子系统运动微分方程(2)的特征方程为:
det{I-Λ[H(ω)]}=0 (4)
则系统特征方程(4)的特征根Λ可写为:
Λ=H(ω)-1 (5)
由频响函数(式(3))中所含各变量可知,对于具有几何非线性柔性叶片和激励依赖非线性围带阻尼的大型汽轮机转子系统,在不同运行工况下,其频率比阻尼比(各部件阻尼等)呈现不同的非线性特性,使系统的频响特性十分复杂。依据特征乘子与单位圆的关系判定转子系统稳定性,关键是获得转子系统特征方程(式(4))的特征根Λ从而求解特征乘子|μF|(式(1)),由式(5)知特征根Λ为系统频响函数(式(3))的逆,则系统稳定性的判定可归结为求解系统在不同运行工况下的频响函数问题。
进一步地,据式(3),频响函数涉及系统各阶固有频率ωnr、模态质量mr、刚度kr和阻尼Cr,其中阻尼Cr=Cl+Cs,对于有围带的转子系统由于系统的结构阻尼Cl远小于围带阻尼Cs,为方便求解,忽略结构阻尼Cl,则有Cr=Cs。需分别求解各项以获得频响函数。
步骤(1-2-2):计算各阶固有频率ωnr、模态质量mr和刚度kr:
步骤(1-2-2-1):建立几何模型并离散划分网格:根据设计的汽轮机组转轴的长度和直径、轴承支承跨度、叶片参数(安装角、扭转角、展弦比)以及围带参数(接触角β、接触间隙e0)等建立转子系统三维几何模型,并将其导入有限元分析软件(以下简称FEMS),在模态分析模块(Modal)中的网格离散模块(mesh)中将转子系统离散为三维实体单元;
步骤(1-2-2-2):设置边界条件:在FEMS的模态分析(Modal)模块中的支承参数设置选项(Fix support)中利用滑动轴承参数计算软件(XROTOR)计算获得支承阻尼CO和支承刚度KO,在轴承支承处添加弹簧接触(connection-spring),设置CO和KO;
步骤(1-2-2-3):在FEMS的模态分析模块(Modal)中求解转子系统模态质量mr和刚度kr;再对整个转子系统施加转速(insert-Rotational Velocity),求解(solution)得到忽略围带阻尼的固有频率ωnr;
步骤(1-2-3):计算获得围带阻尼Cs;
据多频激励下大型汽轮机转子围带摩擦阻尼等效模型和基于瞬态响应包络线的围带阻尼定量求解方法,多频激励F(t)下的围带阻尼Cs为:
步骤(1-2-3-1):计算获得系统激励F(t):运行在湿蒸汽等流体环境的大型汽轮机转子系统,所受的激励力F(t)由两部分组成:
F(t)=Ff(t)+Fm(t) (9)
其中,Fm(t)是相应运行转速下由动不平衡效应等引起的机械激励;Ff(t)是湿蒸汽等流体介质在叶片和围带表面的瞬态分布压力及脉动等流体多频激励;
步骤(1-2-3-1-1):机械激励Fm(t)可通过在FEMS的瞬态响应求解模块转速激励输入选项中选择运行转速作为Fm(t)激励输入;
步骤(1-2-3-1-2):流体激励Ff(t),作为实施方式之一,可利用基于流场仿真的流体激励求解方法得到如式(10)和(11)所示的流体激励Ff(t):
Ff(t)=NTdf(t)T (10)
其中,pi(x,y,z,t)和di分别为流场计算域中任意节点i(离散点)上在不同运行工况(进气流量G、进气速度V、进气角η和转速Ω)下的瞬态压力及其距叶片结构节点j(插值点)的距离;n为离散点的数量,具体实施例中可取n=4;N为叶片六面体单元形函数。
步骤(1-2-3-1-3):利用基于FFT的多谐波平衡法确定该流体激励的前h阶频率成分(可取h≥5)以及各阶谐波频率下的谐波激励分量,将如式(10)和(11)表示的Ff(t)转化为一个常量和若干个不同频率成分谐波激励分量的叠加,如式(12)所示,可在FEMS的瞬态响应求解(Transient)模块中作为激励函数输入。
步骤(1-2-3-2):获得系统多频激励F(t)下转子系统的最大振动位移响应,确定对数减幅比,具体步骤如下:
步骤(1-2-3-2-1):设置边界条件:在FEMS的瞬态响应求解模块中的模型(Model)选项,在相邻围带接触面间添加摩擦接触(insert-connection-fraction),设置摩擦系数μ和接触刚度Kh;
步骤(1-2-3-2-2):施加激励并求解:在FEMS的瞬态响应求解模块施加F(t),具体的,Fm(t)在转速激励输入选项(inertial)中选择运行转速作为激励输入;Ff(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)所获转子系统质量mr和忽略围带阻尼的固有频率ωnr及步骤(1-2-3-2-3)所获的各阶对数减幅比δh代入式(7)计算围带阻尼并将其代入式(6)求解得到多频激励F(t)下的围带阻尼Cs。
步骤(1-2-3):利用有限元分析软件求解转子系统频响函数并获取特征乘子:
据步骤(1-2-1)-(1-2-2)获得的转子系统各阶固有频率ωnr、模态质量mr、刚度kr和围带阻尼Cs在FEMS的谐响应分析模块中求解获得转子系统的实部-虚部形式的频响函数(H(ω)r),并读取各阶特征根Λr,将其代入式(1)即获得转子系统各阶特征乘子|μF|r。
进一步地,在所述步骤(2)中,确定多目标动态优化的目标函数和优化设计变量:
目标函数1:f1(x):汽轮机组的供电/热总效率达到最大值,即
f1(x)=max[LH s]=max[BGH s] (13)
其中,上标s表示当前运行工况是稳定的;LH为该运行工况下输出的供电/热总效率,可表示为:
LH=BGH (14)
其中,GH为实际用于供电/热的进气流量(kg/h),按汽轮机组的两种不同用途划分为两大类运行工况:1)常规供电:低压缸的相对进气流量G/G0=1(其中,G为实际进气流量、G0为额定进气流量)为额定进气流量工况,GH=0;2)热电联产:对于无备用低压缸的抽气式热电联产机组,供热季低压缸的相对进气流量G/G0≤1(称小流量工况),供热季实际用于供热的流量GH=G0-G。B为单位进气流量GH下考虑运行过程中能量转换消耗及效率等因素的产能效率系数。
确定f1(x)的优化设计变量:优化设计变量一定是可控可调的初始变量,在设计阶段已知功率、效率曲线等约束条件下,选取整个运行周期不同运行工况的工况参数:相对进气流量G/G0及其相应的进气速度V和进气角η,转子运行转速Ω作为优化设计变量;
目标函数2:f2(x):转子系统运行稳定性极限阈值区间最优,由步骤(1-2)中稳定性极限阈值求解方法可知,对于有围带起主要减振增稳作用的大型汽轮机组,可通过优化围带阻尼进而优化转子系统稳定性极限阈值-特征乘子μF,故取优化目标函数2为各稳定运行工况下的围带阻尼最优,即:
f2(x)=max[Cs s]=max[f(μs,Ks h,βs,e0 s)] (15)
其中,上标s表示稳定的运行工况;Cs为该运行工况下转子系统围带阻尼;μ、Kh、β和e0依次分别为相邻围带接触面的摩擦系数、接触刚度系数、接触角和初始间隙。
确定f2(x)的优化设计变量:在设计阶段,转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数初步设计确定的前提下,目标函数2的优化设计变量主要为围带参数:μ,Kh,β和e0。
进一步地,在所述步骤(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.为避免动静件碰摩故障,转子系统各节点的各向最大振动位移响应qimax(x,y,z)不超过设计设定振幅值(参照相应的设计手册),即满足:qimax(x,y,z)≤[q(x,y,z)];
e.为适应实际工况,转子系统的转速Ω小于等于常规运行转速的1.2倍;系统的发电/热功率P大于等于额定发电/热功率。
进一步地,在所述步骤(4)中,采用微量摄动修改法对转子系统稳定性进行有约束的多目标动态优化,具体步骤如下:
步骤(4-1):确定满足目标函数2的围带参数:给定第i组运行工况参数(转速Ωi,实际进气流量Gi(Vi,ηi)),并根据步骤(1-2)中的围带阻尼求解方法求解确定该运行工况下围带阻尼最优的第j组围带参数(摩擦系数μj,接触刚度Khj,接触角βj,初始间隙e0j),满足目标函数2;
步骤(4-2):利用步骤(1)中的稳定性判定准则判定第i组运行工况参数(Ωi,Gi(Vi,ηi))下具有第j组围带参数(μj,Khj,βj,e0j)的转子系统的稳定性,若系统失稳,则进入步骤(4-2-1);若稳定,则进入步骤(4-2-2);
步骤(4-2-1):在不同运行工况的三维坐标系(Ω,V,η)下标记当前工况点(Ωi,Vi,ηi)为具有第j组围带参数的转子系统在该工况下的失稳区;
步骤(4-2-2):判断第i组运行工况(Ωi,Gi(Vi,ηi))下的供电/热效率Ls H是否满足目标函数1(式(13)-(14)),若不满足,则辨识确定敏感优化设计变量,在满足步骤(3)约束条件下选取一组微摄动变量(ΔΩ,ΔV,Δη),摄动修改第i组运行工况参数(其中,1)纯供电机组,微摄动变量ΔV=0和Δη=0,ΔΩ≠0;2)热电联产机组供热季:微摄动变量ΔV、Δη、ΔΩ≠0)为第i+1组运行工况(Vi+1=Vi+ΔV,ηi+1=ηi+Δη,Ωi+1=Ωi+ΔΩ),并返回步骤(4-1)计算在第i+1组运行工况下围带阻尼最优的第j+1组围带参数(μj+1,Khj+1,βj+1,e0j+1)和步骤(4-2)判定第i+1组运行工况下是否稳定;若不稳定,进入步骤(4-2-1);若稳定,则进行步骤(4-2-3);
步骤(4-2-3):第j+1组围带参数(μj+1,Khj+1,βj+1,e0j+1)即为第i+1组运行工况下满足两个目标函数的最优围带参数;记录第j+1组围带参数以及对应的最优供电/热效率和稳定性极限阈值;
步骤(4-3):通过上述步骤扫描所有运行工况,获得动态优化后的转子系统在各运行工况的稳定性极限阈值,并将其组合绘制可得到转子系统的三维稳定性极限阈值曲线。
结合具体实例对本公开进行进一步说明:
以某1000MW级大型汽轮机组为例,该汽轮机组是无备用低压转子的热电联产机组,换热季需抽取低压缸部分进气GH用以供热以满足热电联产需求,低压缸的实际进气流量G=G0-GH,即低压缸的相对进气流量G/G0≤1,为小流量工况,由于供热季实际用于供热的流量GH决定机组的供热效率,故对于目标函数1,低压缸的进气流量G越小越优化。依据该机组供电/热功率和效率曲线初步设计的低压末级转子结构如图2所示,该转子采用平行四边形围带,初始围带参数设置为:摩擦系数μ=0.05,接触刚度系数Kh=0.2,接触角β=15°,初始间隙e0=1mm;常规转速Ω=3000rpm。采用本公开提出的基于围带阻尼的汽轮机转子系统稳定性动态优化方法对该汽轮机低压末级转子系统进行动态优化设计。
步骤(1)获取汽轮机转子系统的运行工况参数;建立转子系统稳定性判定准则和稳定性极限阈值的求解方法:
步骤(1-1)获取汽轮机转子系统的运行工况参数;选择据特征乘子与单位圆关系判定转子系统稳定性的准则。
步骤(1-2):建立稳定性极限阈值的求解方法:1)依据该机组供电/热功率和效率曲线初步设计的转子系统结构建立转子系统三维几何模型如图2,并进一步离散为三维实体单元,设置边界条件,求解转子系统模态质量mr和刚度kr;再对整个转子系统施加转速(insert-Rotational Velocity),求解(solution)得到忽略围带阻尼的固有频率ωnr;2)进一步采用基于流场仿真的流体激励求解方法求解转子系统多频激励F(t);3)根据1)、2)利用围带阻尼定量求解方法,求解多频激励F(t)下的围带阻尼Cs;4)根据1)、2)、3)结果求解频响函数获取步骤(1-1)中判定准则所需要的特征乘子。
步骤(2)确定多目标动态优化的目标函数和优化设计变量:
目标函数1:f1(x):汽轮机组的供电/热总效率达到最大值,即
f1(x)=max[LH s]=max[BGH s]
LH=BGH
其中,GH=G0-G,换热季低压缸的相对进气流量G/G0≤1,为小流量工况,实际用于供热的流量GH决定机组的供热效率,满足优化目标1的条件是在保证稳定运行前提下低压缸的流量G越小即机组供电/热功率越大;初设产能效率系数B=90%。
确定f1(x)的优化设计变量:初步设计阶段,在已知功率、效率曲线等约束条件下,选取整个运行周期各运行工况(包括升速启动、甩负荷停车、不同转速及不同流量等)参数:相对进气流量G/G0及相应的进气速度V、进气角η和转子运行转速Ω作为优化设计变量;作为一种实施方式之一,各变量及取值范围采用多因素多水平正交法结合响应曲面法确定。
目标函数2:f2(x):转子系统运行稳定性极限阈值区间最优:
f2(x)=max[Cs]=max[f(μs,Ks h,βs,e0 s)]
其中,上标s表示稳定的运行工况;Cs为该运行工况下转子系统围带阻尼。确定f2(x)的优化设计变量:在设计阶段,转子系统(包括各轴段、叶片、围带及轴承支承等)的结构及相关参数初步设计确定的前提下,目标函数2的优化设计变量主要为围带参数:相邻围带接触面的摩擦系数μ,接触刚度系数Kh,接触角β、e0、μ、Kh、β和初始间隙e0。各变量及取值范围采用多因素多水平正交法结合响应曲面法确定。
步骤(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.为避免动静件碰摩故障,转子系统各节点的各向最大振动位移响应qimax(x,y,z)不超过设计设定振幅值(参照相应的设计手册),即满足:qimax(x,y,z)≤[q(x,y,z)];
e.为适应实际工况,转子系统的转速Ω小于等于常规运行转速的1.2倍;系统的发电/热功率P大于等于额定发电/热功率。
据上述约束条件,初步确定各优化设计变量的优化取值范围为:
运行工况参数:Ω优化取值范围:0<Ω≤4000rpm;额定工况转速为Ω=3000rpm;G/G0优化取值范围:0.15<G/G0≤1,其中相应的V优化取值范围:45m/s<V≤300m/s,η优化取值范围:90°<η≤160°。围带参数:μ优化取值范围:0<μ≤0.2,Kh优化取值范围:0<Kh≤1,β优化取值范围:0°≤β≤20°,e0优化取值范围:0<e0≤1.2。
步骤(4),采用微量摄动修改法对转子系统稳定性进行有约束的多目标动态优化,具体步骤如下:
步骤(4-1)设定初始计算的(i=1)运行工况为该转子系统的额定运行工况(G/G0=1(Vi=300m/s,ηi=90°),Ω=3000rpm),据步骤(1-2)的围带阻尼计算方法从初步设计的(j=1)围带参数((μ1,Kh1,β1,e01))开始计算在(i=1)运行工况下的各组围带参数的围带阻尼,确定该(i=1)运行工况下围带阻尼最优的第j组围带参数(μj,Khj,βj,e0j),并记录对应的围带阻尼最优值
步骤(4-2)据步骤(1-2)稳定性极限阈值计算方法,在有限元分析软件(ANSYS)的谐响应分析模块中求解第(i=1)运行工况下的具有第j组围带参数(μj,Khj,βj,e0j)的转子系统的频响函数,读取各阶特征根,并计算该工况下的系统特征乘子据Floquet理论判定该运行工况下第j组围带参数转子系统的稳定性;
若系统失稳,进入步骤(4-2-1),在不同运行工况的三维坐标系(Ω,V,η)下标记当前(i=1)工况点(Ωi,Gi(Vi,ηi))为具有第j组围带参数的转子系统在该工况下的失稳区;
若稳定,则进入步骤(4-2-2),判断该(i=1)运行工况(Ωi,Gi(Vi,ηi))下的供电/热效率Ls H是否满足目标函数1(式(13)-(14)),若不满足,则辨识确定敏感优化设计变量,热电联产机组供热季:各微摄动变量ΔV、Δη、ΔΩ≠0,且按相对进气量(G/G0)逐渐减小的原则,在满足步骤(3)约束条件下选取一组微摄动变量(ΔΩ,ΔG(ΔV,Δη)),将当前(i=1)运行工况参数摄动修改为第i+1组运行工况(Ωi+1=Ωi+ΔΩ,Gi+1=(Gi+ΔG)(Vi+1=Vi+ΔV,ηi+1=ηi+Δη)),并返回步骤(4-1)计算在第i+1组运行工况下围带阻尼最优的第j+1组围带参数(μj+1,Khj+1,βj+1,e0j+1)和步骤(4-2)判定第i+1组运行工况下是否稳定;若不稳定,进入步骤(4-2-1);
若稳定,则进行步骤(4-2-3):第j+1组围带参数(μj+1,Khj+1,βj+1,e0j+1)即为第i+1组运行工况下满足两个目标函数的最优围带参数;记录第j+1组围带参数以及对应的最优供电/热效率和稳定性极限阈值;
(4-3)通过上述步骤扫描所有运行工况,获得动态优化后的转子系统在各运行工况的稳定性极限阈值,得到本算例下的最优围带参数组合为(μ=0.23,Kh=0.62,β=8.54°,e0=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)/ω为主振动周期;Λ为系统特征方程的特征根;
将转子系统运动微分方程写为复坐标形式:
-mrωn 2+criωn+kr=F0exp(iωt) (2)
其中,ω为激励频率;ωn为系统固有频率;i2=-1;则转子系统第r阶频响函数H(ω)r写为以下实部-虚部形式:
其中,定义频率比λr=ω/ωnr为系统所受激励频率与系统第r阶固有频率之比;阻尼比ξr=Cr/2mrωnr为系统阻尼Cr与第r阶临界阻尼之比,阻尼Cr=Cl+Cs,对于有围带的转子系统由于系统的结构阻尼Cl远小于围带摩擦阻尼Cs,为方便求解,忽略结构阻尼Cl,则有Cr=Cs;kr为第r阶模态刚度;
转子系统运动微分方程(2)的特征方程为:
det{I-Λ[H(ω)]}=0 (4)
则系统特征方程(4)的特征根Λ可写为
Λ=H(ω)-1 (5)
由频响函数中所含各变量可知,对于具有几何非线性柔性叶片和激励依赖非线性围带阻尼的大型汽轮机转子系统,在不同运行工况下,其频率比阻尼比呈现不同的非线性特性;依据特征乘子与单位圆的关系判定转子系统稳定性,获得转子系统特征方程的特征根Λ从而求解特征乘子|μF|,由式(5)知特征根Λ为系统频响函数的逆,则系统稳定性的判定可归结为求解系统在不同运行工况下的频响函数问题;
据公式(3),频响函数涉及系统各阶固有频率ωnr、模态质量mr、刚度kr和阻尼Cr,其中阻尼Cr=Cl+Cs,对于有围带的转子系统由于系统的结构阻尼Cl远小于围带阻尼Cs,为方便求解,忽略结构阻尼Cl,则有Cr=Cs;分别求解各项以获得频响函数;
步骤(1-2-2):计算各阶固有频率ωnr、模态质量mr和刚度kr;
步骤(1-2-3):计算获得围带阻尼Cs;
步骤(1-2-4):利用有限元分析软件求解转子系统频响函数并获取特征乘子:
据步骤(1-2-1)~(1-2-3)获得的转子系统各阶固有频率ωnr、模态质量mr、刚度kr和围带阻尼Cs在FEMS的谐响应分析模块中求解获得转子系统的实部-虚部形式的频响函数H(ω)r,并读取各阶特征根Λr,r为正整数,将其代入式(1)即获得转子系统各阶特征乘子|μF|r;其中,步骤(1-2-3)的具体步骤包括:
据多频激励下大型汽轮机转子围带摩擦阻尼等效模型和基于瞬态响应包络线的围带阻尼定量求解方法,多频激励F(t)下的围带阻尼Cs为:
步骤(1-2-3-1):计算获得系统激励F(t):运行在湿蒸汽流体环境的大型汽轮机转子系统,所受的激励力F(t)由两部分组成:
F(t)=Ff(t)+Fm(t) (9)
其中,Fm(t)是相应运行转速下由动不平衡效应引起的机械激励;Ff(t)是湿蒸汽流体介质在叶片和围带表面的瞬态分布压力及脉动流体多频激励;
步骤(1-2-3-2):获得系统多频激励F(t)下转子系统的最大振动位移响应,确定对数减幅比;
步骤(2):确定多目标动态优化的目标函数及其优化变量:
目标函数1,f1(x):稳定运行下汽轮机组供电/供热效率最大,涉及的优化变量为运行工况参数:转子运行转速Ω,相对进气流量G/G0、进气速度V和进气角η;其中,G为实际进气流量,G0为额定进气流量;
目标函数2,f2(x):转子系统运行稳定性极限阈值区间最优,即转子系统的围带阻尼最优,涉及的优化变量为围带参数:相邻围带间的摩擦系数μ、接触角β、接触刚度Kh和初始间隙e0;
步骤(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组运行工况下的供电/热效率Ls 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计算获得支承阻尼CO和支承刚度KO,在轴承支承处添加弹簧接触connection-spring,设置CO和KO;
步骤(1-2-2-3):在FEMS的模态分析模块Modal中求解转子系统模态质量mr和刚度kr;再对整个转子系统施加转速insert-Rotational Velocity,求解solution得到忽略围带阻尼的固有频率ωnr。
3.如权利要求1所述的方法,其特征是,步骤(1-2-3-1)的具体步骤包括:
步骤(1-2-3-1-1):机械激励Fm(t)可通过在FEMS的瞬态响应求解模块转速激励输入选项中选择运行转速作为Fm(t)激励输入;
步骤(1-2-3-1-2):流体激励Ff(t),利用基于流场仿真的流体激励求解方法得到如式(10)和(11)所示的流体激励Ff(t):
Ff(t)=NTdf(t)T (10)
其中,pi(x,y,z,t)和di分别为流场计算域中任意节点i上在不同运行工况下的瞬态压力及其距叶片结构节点j的距离;n为离散点的数量;N为叶片六面体单元形函数;
步骤(1-2-3-1-3):利用基于FFT的多谐波平衡法确定该流体激励的前h阶频率成分以及各阶谐波频率下的谐波激励分量,将如式(10)和(11)表示的Ff(t)转化为一个常量和若干个不同频率成分谐波激励分量的叠加,如式(12)所示,可在FEMS的瞬态响应求解Transient模块中作为激励函数输入;
4.如权利要求1所述的方法,其特征是,步骤(1-2-3-2)的具体步骤包括:
步骤(1-2-3-2-1):设置边界条件:在FEMS的瞬态响应求解模块中的模型Model选项,在相邻围带接触面间添加摩擦接触,设置摩擦系数μ和接触刚度Kh;
步骤(1-2-3-2-2):施加激励并求解:在FEMS的瞬态响应求解模块施加F(t),具体的,Fm(t)在转速激励输入选项inertial中选择运行转速作为激励输入;Ff(t)在载荷激励输入选项insert-Force中以的形式施加由多阶谐波激励分量组成的多频激励力;在求解solution选项中求解第h阶谐波激励分量下的转子系统振动位移响应;重复上述步骤,依次分别求解出各阶谐波激励分量下的转子系统振动位移响应;
步骤(1-2-3-2-3):确定第h阶谐波激励分量下的最大振动位移响应节点,输出该节点的最大振动位移响应曲线,读取相邻振动位移响应峰值的振幅和并依次代入式(10)求解出各阶振动位移响应对数减幅比δh;
5.如权利要求1所述的方法,其特征是,在所述步骤(2)中,确定多目标动态优化的目标函数和优化设计变量:
目标函数1:f1(x):汽轮机组的供电/热总效率达到最大值,即
f1(x)=max[LH s]=max[BGH s](13)
其中,上标s表示当前运行工况是稳定的;LH为该运行工况下输出的供电/热总效率,表示为
LH=BGH (14)
其中,GH为实际用于供电/热的进气流量,单位:kg/h,按汽轮机组的两种不同用途划分为两大类运行工况:
1)常规供电:低压缸的相对进气流量G/G0=1,低压缸的相对进气流量为额定进气流量工况,GH=0;其中,G为实际进气流量、G0为额定进气流量;
2)热电联产:对于无备用低压缸的抽气式热电联产机组,供热季低压缸的相对进气流量G/G0≤1,供热季实际用于供热的流量GH=G0–G;B为单位进气流量GH下考虑运行过程中能量转换消耗及效率因素的产能效率系数;
确定f1(x)的优化设计变量:优化设计变量是可控可调的初始变量,在设计阶段已知功率、效率曲线约束条件下,选取整个运行周期不同运行工况的工况参数:相对进气流量G/G0及其相应的进气速度V和进气角η,转子运行转速Ω作为优化设计变量;
目标函数2:f2(x):转子系统运行稳定性极限阈值区间最优,由步骤(1-2)中稳定性极限阈值求解方法可知,对于有围带起主要减振增稳作用的大型汽轮机组,可通过优化围带阻尼进而优化转子系统稳定性极限阈值-特征乘子μF,故取优化目标函数2为各稳定运行工况下的围带阻尼最优,即:
f2(x)=max[Cs s]=max[f(μs,Ks h,βs,e0 s)] (15)
其中,上标s表示稳定的运行工况;Cs为该运行工况下转子系统围带阻尼;μ、Kh、β和e0依次分别为相邻围带接触面的摩擦系数、接触刚度系数、接触角和初始间隙;
确定f2(x)的优化设计变量:在设计阶段,转子系统的结构及相关参数初步设计确定的前提下,目标函数2的优化设计变量为围带参数:μ,Kh,β和e0。
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.为避免动静件碰摩故障,转子系统各节点的各向最大振动位移响应qimax(x,y,z)不超过设计设定振幅值,即满足:qimax(x,y,z)≤[q(x,y,z)];
e.为适应实际工况,转子系统的转速Ω小于等于常规运行转速的1.2倍;系统的发电/热功率P大于等于额定发电/热功率。
Applications Claiming Priority (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN2018116315781 | 2018-12-29 | ||
CN201811631578 | 2018-12-29 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111062177A CN111062177A (zh) | 2020-04-24 |
CN111062177B true CN111062177B (zh) | 2021-09-03 |
Family
ID=70301305
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201911326462.1A Expired - Fee Related CN111062177B (zh) | 2018-12-29 | 2019-12-20 | 一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111062177B (zh) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111859730B (zh) * | 2020-06-05 | 2023-10-13 | 合肥通用机械研究院有限公司 | 一种燃料电池离心压缩机转子构形优化设计方法 |
CN112861291B (zh) * | 2021-03-17 | 2022-03-25 | 浙江理工大学 | 一种基于阻尼比分析的多级离心泵环形密封设计方法 |
CN113312730B (zh) * | 2021-06-25 | 2022-07-22 | 内蒙古京泰发电有限责任公司 | 一种双驱汽轮机转子应力监测方法 |
CN113742931B (zh) * | 2021-09-13 | 2024-01-26 | 中国电子信息产业集团有限公司第六研究所 | 一种区块链边缘安全检测方法、系统、电子设备 |
CN118094784B (zh) * | 2024-04-28 | 2024-07-09 | 东北大学 | 一种联合激励下双自由度附件系统中联合响应的确定方法 |
Family Cites Families (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101122541B (zh) * | 2007-08-03 | 2010-08-25 | 东方电气集团东方汽轮机有限公司 | 汽轮机叶片振动试验方法及装置 |
JP5591152B2 (ja) * | 2011-02-28 | 2014-09-17 | 三菱重工業株式会社 | タービン動翼 |
CN105422187A (zh) * | 2015-12-17 | 2016-03-23 | 东方电气集团东方汽轮机有限公司 | 660mw等级汽轮机的给水系统及其给水泵汽轮机末级动叶片 |
CN106528982B (zh) * | 2016-10-26 | 2019-08-23 | 西安交通大学 | 一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法 |
CN106503375B (zh) * | 2016-10-28 | 2020-01-10 | 山东大学 | 一种基于cn群理论确定汽轮机转子临界转速的方法及系统 |
CN108087048B (zh) * | 2017-12-05 | 2020-01-31 | 东方电气集团东方汽轮机有限公司 | 热电联产汽轮发电机组的运行方式 |
-
2019
- 2019-12-20 CN CN201911326462.1A patent/CN111062177B/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN111062177A (zh) | 2020-04-24 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111062177B (zh) | 一种基于围带阻尼的汽轮机转子系统稳定性动态优化方法 | |
CN109657397B (zh) | 基于频响函数的汽轮机叶片-转子系统稳定性的预测方法 | |
Piollet et al. | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37 | |
Qin et al. | Multi-objective optimization design on high pressure side of a pump-turbine runner with high efficiency | |
CN110929419B (zh) | 基于围带零阻尼的汽轮机转子系统失稳极限快速预测方法 | |
Sun et al. | Investigation of tip clearance flow effects on an open 3d steam turbine flutter test case | |
Diener et al. | Multi-Disciplinary optimization of a mixed-flow compressor impeller | |
Chelabi et al. | Analysis of the threedimensional accelerating flow in a mixed turbine rotor | |
Schuff et al. | Influence of the steady deformation on numerical flutter prediction for highly loaded and flexible fan blades | |
Naung et al. | An Experimental and Numerical Study on the Aerodynamic Performance of Vibrating Wind Turbine Blade with Frequency- Domain Method | |
Колодяжна et al. | Aeroelastic Characteristics of Rotor Blades of Last Stage of a Powerful Steam Turbine | |
Li et al. | Research on aerodynamic damping of bladed disk with random mistuning | |
CN114065423B (zh) | 快速评估航空发动机风扇叶片颤振的方法 | |
Schobeiri et al. | Endwall contouring using continuous diffusion: a breakthrough method and its application to a three-stage high pressure turbine | |
Zheng et al. | Aeroelastic vibration analysis of a 1.5 stage compressor | |
Kim et al. | Non-axisymmetric endwall profile optimization of a high-pressure transonic turbine using approximation model | |
Fruth et al. | Influence of the Blade Count Ratio on Aerodynamic Forcing: Part II—High Pressure Transonic Turbine | |
Xie et al. | Strength design and numerical analysis of radial inflow turbine impeller for a 100kw microturbine | |
Zhang et al. | Mistuning effects on aero-elastic stability of contra-rotating turbine blades | |
EP3536949A1 (en) | Wind farm, and operation method, control device, and operation control program for wind farm | |
Choi et al. | Resonant response of mistuned bladed disks including aerodynamic damping effects | |
Fleeter et al. | Fatigue life prediction of turbomachine blading | |
Naung | Computational method for aerodynamic and aeromechanical analysis of offshore wind turbines | |
Popov et al. | Axial compressor optimization method | |
Yang et al. | Numerical research of the ram-rotor with different geometric parameters |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210903 |
|
CF01 | Termination of patent right due to non-payment of annual fee |