抽水蓄能机组推力轴承静力学特性分析

王雪梅，孙建伟，刘 斌，徐立建，孙 洁

（1．河海大学能源与电气学院，江苏省南京市 211100；2．南京市江宁区赵村水库管理所，江苏省南京市211199；3．江苏省骆运水利工程管理处，江苏省宿迁市 223899）

0 引言

随着我国电力产业的不断发展，为保证供电系统的安全、稳定、经济运行[1，2]，发展建设抽水蓄能电站的作用日益凸显，但是随着抽水蓄能电站投入运行，出现了一系列亟待解决的问题[3，4]。推力轴承作为水轮机机组的重要部件之一，承担了水轮发电机组的轴向载荷。由于抽水蓄能机组经常频繁启停的工作特性，推力轴承的工作环境更加复杂，而且推力轴承性能不但关系到机组的出力和效率，还会影响到整个机组的安全稳定运行[5]。

保证良好的润滑是推力轴承正常运行的关键[6]，在瓦面形成稳定的油膜才能有效防止出现烧瓦事故。由于油膜厚度较薄，实际测量其特性难以实现，因此大部分学者采用数值模拟研究推力轴承的油膜特性。聂赛等[7]研究了推力轴承油膜在不同转速和不同压力进口的特性，分析了油膜在不同工况下的压力和速度分布规律。张艳琴等[8]为研究油膜温度问题，数值模拟了不同转速和不同载荷工况下的油膜温度分布规律，并通过油膜温度场实验验证了数值模拟结果的准确性。李永海等[9]数值模拟了喷油润滑式可倾瓦推力轴承温度分布，研究了推力瓦在喷油润滑工况下的热变形。刘广东等[10]研究了双距形腔多油垫静压推力轴承油膜在高速重载工况下的温度场和压力场分布，发现偏载距离进一步增加，推力轴承会出现失压现象，发生摩擦学失效。

江西洪屏抽水蓄能电站在启动过程中出现了高压注油泵出口压力过低的情况，当注油泵出口压力过低不足以在推力瓦间形成可靠稳定的油膜，就会导致烧瓦事故。大部分学者关注油膜在高速、重载下的特性，在此基础上对推力轴承进行优化，但是关于推力瓦进口压力造成的影响研究较少，本文基于江西洪屏抽水蓄能电站出现的情况，建立推力轴承模型，基于流固耦合理论，采用数值模拟计算高压注油泵压力报低所对应转速下的油膜特性，进而探究在不同高压注油泵出口压力下油膜对推力轴承的作用机理。

1 模型的建立

1.1 推力轴承油膜的建立

一般情况下推力轴承由推力头、推力瓦、镜板、支撑系统、供油系统组成。如图1所示为水轮机组推力轴承主要结构示意图。其中推力瓦为静止部件，推力头和镜板连接随着机组一同旋转，在机组启动时高压注油泵通过推力瓦进油孔注油，形成静压油膜。当机组到达稳定转速高压注油泵不再注油，润滑油通过瓦面的楔形斜面达到润滑效果，此时为动压油膜[11]。

1.2 固体模型的建立

本文主要研究油膜对推力轴承的作用机理，为便于分析，将整体分为两大部分：流体区域和固体区域[12]。以江西洪屏抽水蓄能电站的推力轴承为研究对象，建立固体三维模型。其中江西洪屏抽水蓄能电站推力轴承由12块推力瓦组成，采用弹簧簇支撑[13]为了简化计算，认为每一块推力瓦瓦间油膜分布情况和运动规律基本一致，且润滑油流过轴瓦与镜板之间对推力轴承产生的作用具有周期性，因此采用一块推力瓦为研究对象，相应地推力头也采用其1/12进行建模。如图2所示为固体部分三维模型。

图2 固体部分三维模型Figure 2 Solid part 3D model

2 数值模拟计算

2.1 流固耦合基本理论

对推力轴承进行流固耦合计算时需先计算得到油膜的流体计算结果，在流体计算时，认为所计算的流体为连续介质，满足如下雷诺方程[14]：

式中：r——柱坐标系径向方向；

θ——柱坐标系圆周方向；

Ω——机组转速；

p——油膜压力；

h——油膜厚度；

μ——油的动力黏度。

润滑油的黏度受温度影响，润滑油和温度之间的关系可用如下黏温方程表示[15]：

式中：μ0——润滑油在温度T0下的动力黏度；

T——绝对温度。

在进行流固耦合分析时，认为固体是弹性的，在此基础上结合结构动力学方程和流体动量方程分析，得到如下方程[16]：

式中：[M] ——质量矩阵；

[K] ——刚度矩阵；

[C] ——阻尼矩阵；

{F}——外荷载向量；

{u}——位移向量。

在进行流固耦合计算中，需要通过流固耦合面进行流体和固体的数据传导，数据传导需要满足压力、位移、热流以及温度守恒[15]，如式（4）所示：

式中：f——流体；

s——固体；

τ——流体与固体应力；

u——位移；

q——热流量；

T——温度。

2.2 网格划分及约束设置

流体部分计算油膜的网格划分最为重要，采用结构网格划分，其余部分采用非结构网格，经网格无关性验证油膜厚度为5层，流体域计算总网格数为65万个。固体区域网格通过ANSYS19.2中Static Structural，设置推力瓦网格尺寸为0.008m，推力头网格尺寸为0.01m，生成非结构网格。

对推力轴承进行流固耦合计算时，需要给结构体添加相应约束，图3（a）为推力瓦所受约束，其中A为对推力瓦顶部的固定约束（fixed support），固定住推力瓦瓦面避免产生位移。B处为重力作用约束（standard earth gravity），推力瓦在工作时受到地心引力的作用，设置重力加速度为9.8066 m/s2。图3（b）为推力头所设置得约束条件，其中A为角加速度（rotational velocity），推力头在旋转过程中受到离心力作用，设定旋转速度为275r/min。B为施加集中应力（force），此处设定为推力头的一个面上，所设置的值为537N，该值为计算工况荷载的1/12加载在模型上[17]。C和D均为圆柱约束（cylindrical support），设定推力头的切向转动。

图3 约束条件设置Figure 3 Constraint setting

2.3 计算工况设置

推力轴承的油膜是通过高压注油泵建立的，洪屏抽水蓄能电站在抽水调相过程中，当机组转速到达额定转速（500r/min）的55%时，高压注油泵压力都会降低至6MPa，直至转速上升至额定转速，高压注油泵持续低压运行。探究高压注油泵出口压力降低过程中对推力轴承结构造成的影响，由于注油泵和推力轴承进油压力是一一对应的，为简化计算取额定转速的55%，设置不同推力轴承进口压力（4～11MPa）进行数值模拟计算。

2.4 数值计算方法

推力轴承在运行过程中与油膜接触处因为油膜压力分布差异会引起一定的位移量，为探究油膜对推力轴承产生的影响，可通过ANSYS19.2中Static Structural进行计算[18]。首先在Fluent中获得油膜的数值计算结果，通过建立流固耦合面进行数据传导，将油膜表面压力通过镜板面加载到推力轴承上，之后在Static Structural中进行稳态结构计算，获得各个网格节点的位移变化。然后根据耦合面位移变化进行流体区域相应面网格的变形及体网格重构，接下来再进行下一轮的流体计算[19]，下一轮计算重复前面所述过程，经若干轮计算和传递后停止计算，此时的计算结果应满足设定的收敛条件。如图4所示为单向流固耦合计算流程图。

图4 单向流固耦合计算流程Figure4 Fluid-structure coupling calculation process

3 计算结果及分析

3.1 推力轴承油膜压力分析

随着注油泵出口压力的降低，油膜表面的压力也会不可避免地受到一些影响。由于机组在55%额定转速6MPa时出现注油泵出口压力报低情况，因此主要关注6MPa左右油膜压力分布云图，如图5所示为55%额定转速，进口压力为5.5MPa油膜表面压力分布云图。油膜最高压力位于瓦面中心靠近主轴位置，其余位置的压力以该位置为中心递减，在油膜表面最外侧出现负压；随着进口压力的增加，各个油膜表面压力分布情况和进口压力为5.5MPa油膜表面压力分布大致一致，如图6所示为推力轴承不同进口压力和油膜表面最大压力关系图，随着进口压力的增加，油膜表面压力近似呈线性增加。可见进口压力对油膜的分布情况影响较小，主要影响的是油膜表面压力值大小。

图5 5.5MPa压力进口油膜表面压力分布Figure 5 Surface pressure distribution of oil film at 5.5MPa pressure inlet

图6 不同进口压力油膜表面最大压力变化情况Figure6 The maximum pressure change of the oil film surface at different inlet pressures

3.2 推力瓦静力学分析

如图7所示为转速275r/min，不同进口压力的推力瓦的应力分布情况。应力集中于推力瓦中心，随着进口压力的减小最大应力值减少，且应力分布逐渐不均匀，出现明显的梯度分布，在4～6MPa这一区间的进口压力变化中，推力瓦瓦面伴随着明显的压力梯度变化，也逐渐出现应力集中从瓦面中心转移至出油口左上角位置处的现象，说明进口压力降低到一定程度时，会出现多个应力集中区域；虽然在较低压力进口推力瓦瓦面出现多个区域应力集中，但是整体应力值随着进口压力的降低而降低，出现应力集中的区域应力值并不大。

为探究上述进口压力为4～6MPa时，推力瓦表面出现的多处应力集中对推力瓦造成的影响，以5.5MPa进口压力推力瓦瓦面应变分布云图为例，如图9所示，可以看出图7（c）中应力集中的区域对应图7中发生较大形变的位置，说明应力很大程度上影响着推力瓦瓦面的形变情况，其余进口压力的应变分布云图和图7中推力瓦应力分布云图相对应。进而可以说明，进口压力从6MPa开始出现多个最大形变中心。如图9所示为推力瓦在不同进口压力情况下的应变情况。结合图8中可以看出，发生最大形变的位于瓦面中心的位置，随着压力的减小该位置处的最大形变不断减小，瓦面最大和最小形变差异最大的发生在较高进口压力11MPa处，推力瓦瓦面发生的形变不到推力瓦外径的0.01%，对机组造成的影响较小。

图7 不同进口压力推力瓦应力分布Figure7 Stress distribution of thrust pads with different inlet pressures

图8 5.5MPa进口压力推力瓦应变分布（m）Figure8 Strain distribution of 5.5MPa inlet pressure thrust pad

图9 不同进口压力推力瓦应变变化情况Figure 9 Strain changes of thrust pads with different inlet pressures

3.3 推力头静力学分析

如图10所示为不同进口压力下推力头的应力变化情况图，可以看出推力头表面最大应力随着压力减低而降低，4MPa对应的最大应力和11MPa对应的最大应力相差近60%；但是不同进口压力之间的应力平均值相差不大，这可能是因为推力头表面最大应力集中区域较小，从而对推力头表面的平均值影响较小。如图11所示为5.5MPa进口压力推力头表面应力分布图，可见推力头表面出现应力集中的区域较小并且位于推力瓦形状发生突变位置，应考虑对推力头进行局部应力集中处理或者对应力集中处进行结构优化。

图10 不同进口压力推力头应力变化情况Figure10 Stress changes of thrust heads with different inlet pressures

图11 5.5MPa进口压力推力头应力分布（Pa）Figure 11 Stress distribution of thrust head with 5.5MPa inlet pressure

如图12所示为推力头在275r/min不同进口压力下的形变变化图，总体上看随着进口压力的减小推力头的最大形变逐渐增大，进口压力降低到2MPa时，最大形变增加了原来的1.5倍，相比较油膜对推力瓦形变产生的影响，推力头的形变比推力瓦的形变高出一个数量级，说明油膜对推力头的影响更大一些。从图13中可以看出推力头在不同进口压力情况下，最大形变均在推力头外沿，且从靠近主轴处形变逐渐增加直至推力头外沿，呈现外延上翘的形态。由于油膜较薄，而推力瓦产生的形变远远大于油膜的厚度，因此推力头受油膜影响产生的形变在一定程度上也影响了油膜的均匀分布。因此进口压力应控制在一定范围内，如果出现突然的压力降低会在短时间内出现较大的形变变化，进一步影响油膜的均匀分布，可能对机组运行产生不利影响。

图12 不同进口压力推力头应变变化情况Figure12 Strain changes of thrust heads with different inlet pressures

图13 5.5MPa进口压力推力头应变分布（m）Figure13 Strain distribution of thrust head with 5.5MPa inlet pressure

4 结论

本文通过分别建立推力瓦和推力头三维模型，基于流固耦合理论，对不同进口压力下的推力瓦和推力头进行了静力学分析，得出以下结论：

（1）油膜表面最大压力随进口压力的增加而增加，且呈线性增长趋势；不同进口压力对油膜的压力分布影响较小，主要影响油膜表面压力值大小。

（2）推力瓦表面应力随着进口压力的减小而减小，但是随着进口压力的减小会出现多处应力集中现象；相应地不同进口压力下推力瓦瓦面发生的形变不到推力瓦外径的0.01%，对机组造成的影响较小。

（3）推力头的最大形变随着进口压力的降低而减小，4MPa对应的最大应力和11MPa对应的最大应力相差近60%；其表面出现应力集中的区域位于推力瓦形状发生突变位置，应考虑对推力头进行局部应力集中处理或者对应力集中处进行结构优化。

（4）整体上推力头形变随着进口压力的降低而增加，2MPa对应的最大形变增加了原来的1.5倍，进口压力应控制在一定范围内，如果出现突然的压力降低会在短时间内出现较大的形变变化，进一步影响油膜的均匀分布，可能对机组运行产生不利影响。