基于水力振源合理施加的水电机组轴系统振动分析

王 建， 宋志强， 刘昱杰

(西安理工大学 水利水电学院， 陕西 西安 710048)

1 研究背景

水电在国家经济、能源发展及环境安全等方面有着巨大的影响，随着水轮发电机容量和机组尺寸的增大，水电站机组及厂房的振动问题成为制约机组高效运行，甚至威胁其安全的关键。引发水电站机组及厂房振动的因素除了机组和厂房结构本身的非线性特性外，水力、机械、电磁振源的复杂性及其相互作用也是一个重要的方面。由于过流部件构造复杂，水轮机运行工况多变使得获知、表征和施加水力振源成为一大难题。转轮受到水流作用产生一种近似随机的变值荷载，在以往的轴系统分析中，只能将其忽略或简化为只考虑一种优势频率的简谐荷载，该荷载的无法准确模拟给轴系统振动分析造成了很大影响[1-3]。

关于水轮发电机组水力振源的研究目前已经取得了较大进展，文献[4]应用Newmark算法计算了不同工况条件下水轮机转轮部位应力，考虑了周期性的水力振源作用下结构的振动分析，但并未考虑脉动压力的非周期性部分；文献[5]运用转子动力学软件ARMD给予了机组轴系横向振动的评价，并准确预估了不同工况下水力激励力对机组瞬态响应的振动分析，提出了水力激励力对轴系稳定性研究的重要意义；文献[6]建立ANSYS主轴模型进行了轴系结构自振特性分析，分析了多种外激励共同作用下时振动响应。随着计算流体动力学(CFD)的应用及发展，合理获知机组过流部件表面脉动压力成为了可能[7]，在流体计算结果越来越准确的情况下，将CFD获知的脉动压力导入轴系结构进行水力振源的模拟，有了更为可靠的理论依据和试验参考，很多学者已经在CFD计算的基础上对流道内的复杂流体扰动进行了探究并结合模型试验验证了CFD流体计算的可行性[8-12]，在此基础上，Khare等[13]得出了模型试验所不能观测到的细节部位的流体运行情况，并分别计算了不同导叶开度和转轮转速条件下的流体运行特性；尤建锋等[14]建立了与试验模型一致的完整CFD流体计算模型并进行了多工况的数值模拟，对比分析了压力脉动特性；钱忠东等[15]通过改变导叶开度对比分析了不同工况条件下流道内各部位脉动压力分布规律。上述研究虽然准确获知了流道脉动压力，但更多是研究机组水力效率、转轮叶片疲劳破坏及空化等，合理施加转轮脉动压力进行机组整个轴系统结构振动分析的研究目前尚不多见。

机组在受到转轮处脉动压力作用的同时也受到定转子气隙不均匀产生的不平衡电磁力和转子-转轮偏心产生的离心力影响，本文通过CFD计算提取转轮表面脉动压力时程，转化为荷载文件，实现转轮脉动压力更为合理的模拟和施加，在此基础上综合电磁、机械振源进行了多振源耦合情况下水电机组三维转子-轴承系统结构的振动分析，并对转子-转轮不平衡力相位差等进行了敏感性分析。

2 脉动压力的CFD计算及提取

2.1 全流道流体域计算模型

本文以某电站混流式水轮机为研究对象，其转轮直径2.165 m，转轮叶片13个，固定导叶12个，活动导叶24个，导叶开度角23°，模拟工况条件为最高水头、最优负荷(H=77 m，P=102 MW)和额定水头、额定负荷(H=61 m，P=88 MW)下水轮机的内部流动特性，采用Unigraphics NX 8.0软件建立包括蜗壳、固定导叶、活动导叶、转轮和尾水管在内的全流道计算域模型。由于水轮机部件结构复杂，各部件之间相互关联，为确保模型的完整性，减小分析误差，采用ANSYS ICEM-CFD软件进行高精度的六面体网格划分，流道模型共有4 842 301个单元，4 880 919个节点，流道内各部件网格划分如图1所示。

图1 流体计算域网格划分示意图

2.2 控制方程

对三维不可压缩非稳态Navier-Stokes方程采用时间平均法，进行Reynolds时均法得到连续方程和动量方程分别为式(1)和(2)：

(1)

(2)

湍流模型采用剪切应力输运SST(Shear Stress Transport)模型，用以模拟旋转部件与静止部件的动静干涉，固体壁面采用标准壁面函数进行处理，控制方程为式(3)：

(3)

(4)

式(1)～(4)中：ρ为流体的密度，kg/m3；μ为动力黏性系数；u为表征速度，m/s；τ为雷诺应力张量，MPa；Pk为湍流脉动动能k的生成项，k由k-ε模型求得；Pω为湍流脉动频率ω的生成项，由k-ω求得；Γk，Γω为k与ω的有效扩散系数；Yk，Yω为k与ω的耗散项；Dω为正交扩散项。

2.3 边界条件

计算域的进口为金属蜗壳的进口，流体垂直进口流入，进口边界条件设置为质量-流量边界；出口边界条件为压力出口，指定平均静压为0；固体壁面为无滑移、绝热壁面；交界面方式采用动静转子交界。

2.4 控制方程的离散和求解设置

模型采用有限体积法进行离散，建立多种运行条件下的原型全流道模型，利用大型有限元软件ANSYS中的CFX模块，采用Realizable 模型和SIMPLEC算法对水轮机过流部件的内部流动进行计算，从振源特性和水流特性分别进行流体计算。转轮转速为136.4 rad/min，转频为2.273 Hz，旋转周期为0.44 s，转轮瞬态分析中，将周期分为120个时间步进行脉动压力提取，时间步长取为Δt=0.00367 s，其中，转子-转轮轴系结构旋转角速度为ω=2πn=14.2838 rad/s，定常计算结果提取转轮旋转两个周期的脉动压力时程值。

2.5 脉动压力结果分析

为了直观地分析引水部件和转轮内部的压力脉动情况，分别在水轮机蜗壳进口、转轮进口、叶道间，转轮出口布置了4个压力观测点如图2所示。通过完成蜗壳进口脉动压力的提取并进行傅里叶变换后与模型试验实测数据对比如图3所示。可见CFD脉动压力计算结果可信，并在此基础上进行转轮部位的脉动压力提取。

水轮机转轮是整个水轮发电机组的核心部件，转轮中流态的好坏直接影响了机组的性能。对两种工况转轮处脉动压力提取，处理流体计算结果可知最优工况和额定工况下固定导叶和活动导叶通道内均无脱流，流线比较光顺，固定导叶出口角和活动导叶进口角配合很好，提取的脉动压力结果可以作为本文研究数据参考。现以最高水头H=77 m、最优工况条件下的转轮脉动压力提取结果为例，就转轮脉动压力的时域特性和频域特性进行简要分析。

脉动压力的近似随机性通过提取转轮进口、叶道间、转轮出口3个监测点脉动压力时域曲线进行表征，并就傅里叶变换后的频域曲线进行转轮部位脉动压力分析如图4所示。

由图4中压力脉动时域图(图4(a)、4(c)、4(e))和频域图(图4(b)、4(d)、4(f))可以看出，转轮进口脉动压力幅值较大，转轮进口压力幅值为0.423 MPa，叶道间脉动压力幅值为0.125 MPa，转轮出口脉动压力幅值为0.079 MPa；转轮进口和叶道间主频均为2.23 Hz的低频，该频率是转轮转频2.273 Hz的约1倍，在整个水轮机流道内传播，此外，转轮进口处还存在53.4 Hz的高频，是转轮叶片数(Zr=13)的4倍，该高频脉动是由于转轮和导叶的动静干涉作用产生的，转轮出口主频为接近1倍转频的低频，由各部位频域图对比可以看到叶道间振幅增加明显。

3 轴系统模型建立及荷载施加

3.1 轴系统三维有限元模型建立

利用有限元软件Ansys建立轴系统三维模型，模型由上导轴承、发电机转子、下导轴承、水导轴承、法兰、主轴及转轮组成。转轮采用solid187单元离散，主轴结构采用solid45单元离散。用combin14弹簧单元模拟上导、下导、水导轴承，轴系统单元总数为64 936个，节点数为113 904个，实际单元划分及结构模型如图5所示。

图2 流体计算压力观测点位置示意图 图3 蜗壳进口处脉动压力频域图结果对比

图4 转轮特征部位脉动压力时域图和频域图

导轴承的刚度和阻尼随着轴颈运动及导轴承间隙变化而非线性变化，而且刚度与阻尼之间存在一定的相互关系，即当导轴承的刚度发生变化时阻尼也随之相应变化[1]，如图6所示，其中λk=Kx/μnl(R/C)3，λc=ωCx/μnl(R/C)3,μ为轴承内油的黏度，N·s/cm2；n为机组额定转速，rad/s；l为轴瓦的轴向长度，m；R为导轴承轴瓦半径，m；C为导轴承间隙，m；C′为轴瓦在支承点处的间隙，m；ω为大轴旋转圆频率，rad/s；Kx为导轴承刚度，N/m。

以某实际水电站机组为例，取导轴承轴瓦半径R=0.77 m，转速n=136.4 r/min，ω=14.2838 rad/s，导轴承间隙C=0.02 cm，黏度μ=0.33×10-5N·s/cm2。代入上述数据后，Kx=1.709×106λk，Cx=0.12×106λc。

拟取导轴承刚度为Kx=2×108N/m，λk=1.171，由图查得C′/C=0.79，得λc=2.25，导轴承阻尼系数为Cx=2.69×107N/m，为使计算简便，本文计算模型拟取3个导轴承刚度系数和阻尼系数均为Kx和Cx。

3.2 荷载计算及施加

3.2.1 水力振源 作用在转轮上的水力分为周期性和非周期性两种，分别对应水力不平衡力和脉动压力，周期性力即水力不平衡力，水力不平衡力主要是由转轮部分不圆或机组大轴摆度使迷宫间隙发生周期性变化产生的迷宫压力脉动，取值参照式(5)计算[16]。

F1=K水ε

(5)

式中：K水为水的不平衡系数；ε为转轮处摆度，m；算例取F1=39.12 kN。

图6 导轴承刚度与阻尼系数的关系

近似随机的非周期性水力振源一般难以用传统方法考虑。非周期性脉动压力即上述CFD计算所得转轮表面实时脉动压力。模拟过程是在Ansys软件中建立以坐标原点为基准点的转轮模型，进行网格划分，提取转轮区域各节点编号及其坐标，写入节点文件，同时建立以坐标原点为基准点的全流道计算模型，保证结构体与流道模型的位置对应，在此基础上应用CFD-ICEM进行流体域网格的划分，完成CFD流体非定常计算后，进入CFX-Post后处理中编写节点导入命令，提取各节点部位的脉动压力时程后，利用APDL语言将脉动压力施加到结构中，至此完成脉动压力的合理模拟及施加。

3.2.2 机械不平衡力 假设只存在动偏心，各部位轴心在初始位置都是原始位置零点，系统偏心是由于系统转动所引起的旋转部件质量分布不均匀引起的，所产生的离心力为：

F2=me0ω2

(6)

式中：m为转动部件的质量，kg；ω为转动角频率，rad/s。

离心力主要产生于转子和转轮处，根据已知参数转子质量为mr=6.5×106kg，转轮质量为mt=3.5×106kg，求得F2r=856.913 kN、F2t=583.51 kN ，不平衡力作为外力施加在相应的结构上。偏心距参考水轮机制造标准：转子偏心er=6.3/ω，转轮偏心et=40/ω，取er=0.44 mm，et=1.2 mm。

3.2.3 不平衡电磁拉力 由于制造、安装偏差及运行过程中的磨损、变形等会造成定子和转子间产生不均匀间隙，进而产生不平衡电磁拉力，当转子偏心为e时，不平衡磁拉力为[17]：

(7)

式中:L、D分别为转子高度和直径，m；取L=1 134 cm；D=275 cm;B为磁通量密度，取2393×10-4T；β为系数，一般取值0.3左右；δ为正常间隙，取值为1.6 cm，e=0.44 cm。本文不平衡电磁拉力取F3=152.7 kN。

为了保证数值模拟更接近机组实际运行条件，不平衡电磁拉力和离心力均由互相垂直的正弦荷载和余弦荷载来模拟。

4 基于水力振源合理施加的轴系振动分析

4.1 水力振源对轴系振动的影响

4.1.1 脉动压力对轴系结构的影响分析 将CFD提取的脉动压力结果进行处理后单独施加到轴系结构中进行动力响应分析，并提取转轮、转子周期内的轴心偏移情况如表1，并根据XY向(水平向)位移进行偏移合成绘制转轮、转子轴心偏移图如图7。

结合表1和图7可以看出，在两种工况条件下，导入脉动压力，转轮、转子轴心运行均受到影响。当水头H=61 m时，转轮轴心偏移为0.005～0.238 mm，转子轴心偏移为0.001～0.01 mm；当水头H=77 m时，转轮轴心偏移为0.009～0.233 mm，转子轴心偏移为0.001～0.007 mm。水头H=61 m时转轮、转子位移偏移较大，原因是由于此工况条件下转轮部位脉动压力振幅最大，转轮处于部分负荷状态，振动更加明显，转轮部位脉动压力较大，对整个轴系结构有明显的影响。

表1 两种水头工况下转轮、转子偏移情况

4.1.2 振源耦合条件下轴系统振动分析 通过在转子处施加不平衡电磁拉力模拟电磁振源分布、加载离心力的方式模拟转子和转轮处的机械振源以及导入实时脉动压力时程模拟水力振源，研究整个转轮-转子轴系结构在机械振源、电磁振源、水力振源3种振源耦合情况下的振动响应。图8为轴系结构在只导入机械振源和电磁振源作用下的转轮、转子轴心轨迹图，图9为两种运行工况下，轴系结构在3种振源耦合情况下的振动分析，并根据响应结果绘制的轴系统关键部位轴心轨迹图。

由图8、9转轮和转子轴心轨迹图对比可知，在脉动压力的非线性影响下，结构各部位轴心运行轨迹呈现不规则变化，从图9轴心轨迹图可以看出转轮、下导和水导受脉动压力影响较为剧烈，当水头H=61 m时，导入脉动压力后转子摆动幅值为0.102～0.107 mm，转轮摆动幅值为0.586～0.725 mm，转子摆动幅值增加了4.91%，转轮摆动幅值增加了23.7%，上导轴承、下导轴承和水导轴承的摆动幅值分别为0.072、0.059和0.063 mm；当水头H=77 m时，转子摆动幅值为0.102～0.103 mm，转轮摆动幅值为0.586～0.664 mm，转子摆动幅值增加了0.98%，转轮摆动幅值增加了13.2%，上导轴承、下导轴承和水导轴承的摆动幅值分别为0.071、0.058和0.062 mm。轴系结构各部位偏移情况对比可知H=61 m水头作用下轴系结构振动较H=77m水头更为剧烈。

图7 两种水头工况下转轮、转子受脉动压力影响的轴心偏移图

图8 电磁振源和机械振源耦合作用下转轮、转子轴心轨迹图

图9 3种振源耦合作用下轴系结构关键部位轴心轨迹图

综上可得，脉动压力对整个结构的影响主要体现在转轮上，转轮轴心运行受脉动压力非线性影响较为剧烈，转子由于距离脉动压力施加部位较远，受到的影响较小；在额定水头，额定工况条件下，脉动压力对转轮振动影响高于最大水头，对于最优工况条件下的转轮振动，其脉动压力对转轮结构的振动有着重要影响，为使得分析更为符合实际，脉动压力的合理模拟和施加尤为重要。

4.2 转子与转轮荷载相位差对系统振动的影响

图10 不同相位差条件下各部位轴心偏移对比图

由图10(a)、(b)可以看出不考虑转轮脉动压力影响下，转轮、转子在电磁振源和机械振源耦合作用下振动平稳，转轮、转子位移幅值均出现在θ=180°时，分别为0.707和0.116 mm，可见当θ=180°时为荷载的最不利组合；由图10(c)可以看出，在考虑脉动压力的振源耦合作用下，在最不利荷载θ=180°时，转轮偏移幅值为0.757 mm，且当结构在振源作用0.07 s之后，转轮开始剧烈振动，当θ=90°时，转轮轴心运行轨迹变动幅度最大，振动幅值为0.751 mm，略小于θ=180°时的位移幅值，θ=0°和θ=45°时转轮运行轨迹变动范围较小，且θ=0°时，转轮运行轨迹范围最小，但仍在0.41～0.66 mm之间摆动；由图10(d)可以看出，转子受相位角θ影响较小，但θ=90°时，转子轴心运行轨迹变幅大于其他3种相位差，运行范围在0.031～0.117 mm之间，可见相位角θ大小对转轮运行稳定有着重要影响，对转子影响相对较小。 此外由各图对比可知，在θ=180°的时候，转轮位移幅值由0.707 mm变为0.757 mm，增长幅值为7.08%，转子位移幅值由0.116 mm变为0.117 mm，变化相对较小，原因是因为转轮受到脉动压力直接作用变化较明显，而转子距离转轮位置较远，受脉动压力影响较小。

5 结 论

(1) 建立水轮机全流道湍流模型，基于CFD数值模拟获取转轮叶片表面脉动压力时程，通过Ansys命令流将脉动压力导入轴系统结构的振动分析中，可实现水力脉动振源的合理模拟和施加，在此基础上，综合电磁和机械振源进行了多振源耦合作用下轴系统的振动分析。

(2) 轴系统考虑非周期水力即脉动压力作用后，转子振动变化较小，转轮轴心偏移幅值较大。由于各点脉动压力的方向是垂直该点所在表面的，转轮叶片沿圆周循环对称且本身具有的空间扭曲性，所以非周期性脉动压力对转轮某一个方向上的振幅影响不会很大，但使转轮运行轨迹趋于复杂，振动趋于剧烈。轴系统在非周期性脉动压力作用下的非线性变化明显，合理模拟和施加转轮叶片表面脉动压力是必要的。

(3)转子和转轮荷载相位差为90°时，转子受非周期性脉动压力影响较大，相位差为180°时，转轮受非周期性脉动压力影响较大。导轴承支撑条件和荷载相位差对轴系统振动有着重要影响，在轴系统动力特性分析及评价中值得重视。