汽轮机大尺寸动压轴承特性数值模拟研究*

(华中科技大学能源与动力工程学院 湖北武汉 430074)

2016年，我国核电累计发电量为2 105.19亿kW·h，约占全国累计发电量的3.56%，比2015年同期上升了25.07%[1]。按照我国电力行业发展规划，我国核电装机容量还将继续增加。随着核电在我国的加速布局，核电汽轮机也在不断向更高容量、更大尺寸和重载荷的方向发展。如东方电机股份有限公司某核电机组轴承直径已经达到800 mm，比压达到2.63 MPa[2]。

虽然用于轴承润滑性能研究的流体动力润滑理论及轴承计算的传统方法已经趋于成熟，但是随着工业实践对动压轴承结构的更高要求，按照传统动压理论假设设计的大尺寸动压轴承在运行中出现了一些问题，比如碰摩等不稳定现象时有发生，表明传统动压理论对于大尺寸动压轴承的设计可能存在某些不适用的地方。

传统上，国内外学者一般利用求解基于层流假设的雷诺方程的方法研究油膜轴承动静特性[3-5]。这种方法具有计算时间短的突出优点，但忽略了惯性项、油膜曲率等因素的影响。而且基于长轴承或者短轴承理论对滑动的求解也存在一定的局限性，比如忽略了油膜轴承的具体局部结构，例如油槽等。随着计算流体动力学技术的发展和计算机性能的提高，通过求解N-S方程的方法研究轴承动特性[6-8]已经越来越普遍了。采用CFD三维建模，可以考虑任意结构形式轴承及油膜内部空化效应的影响[9]。通过CFD软件提供的自定义函数实现轴颈的速度和位移扰动，便于分析动压轴承刚度阻尼特性[10]。存在的缺点主要是油膜圆周尺寸相对油膜间隙的大比例增加了高质量网格划分的难度，以及针对流动、空化、两相流不合理不匹配的模型所导致收敛性差的问题。

目前，针对核电汽轮机等行业的大尺寸动压轴承的研究分析还并不充分，文献[2]基于经典润滑理论的广义雷诺方程，结合流固耦合方法考虑轴瓦的热弹性变形，分析了直径为800 mm、1 750 MW级核电半速发电机试验用滑动轴承的润滑特性，轴瓦温度和油膜压力的仿真值与实验值偏差较小，最小油膜厚度和仿真值与测试值偏差15%～22%。文献[11]指出，在大尺寸油膜轴承的实际运行过程中，轴承发散区油楔会产生负压，更易导致润滑油气泡析出，发生空化效应并影响压力场的分布。文献[12]指出，油膜的空化效应将降低承载力并导致轴承结构损伤。

本文作者基于有限体积法，以某核电轴承结构为对象，考虑空化效应和温黏效应，对典型大尺寸动压轴承流动问题进行模拟研究。

1 基本方程和计算方法

1.1 空化理论

动压轴承中的油膜压力是不均匀的，存在高压区和低压区。当油膜压力低于一定值时，润滑油中溶解的空气会析出，该现象称为空化效应。因此在对流场进行求解之前，引入多相流模型。从计算效率和模型适用范围的角度考虑，文中采用Mixture模型来模拟润滑油和空气之间的两相耦合作用。

Singhal et al.模型，又称为全空化模型，是由SINGHAL等[13]于2001年提出。该模型考虑了相变、气泡动力学、湍流压力波动和不凝结气体的影响。两相连续方程为

液相：

(1)

气相：

(2)

混合：

(3)

式中：l表示液相；v表示气相;ρ表示混合密度，

ρ=αρv+(1-α)ρl

(4)

综合式(1)—(3)，混合密度(ρ)和气相体积分数(α)表示为

(5)

在流体和气泡之间具有零速度滑移的流动液体中，气泡动力学方程源于广义Rayleigh-Plesset方程：

(6)

式中：RB为气泡半径；σ为流体表面张力系数;ρl为流体密度；pB为气泡表面压力；p为当地远场压力。

忽略二阶项和表面张力，方程可简化为

(7)

该方程提供了将气泡动力学的效应引入空化模型的物理方法，也是混合物密度的方程。蒸汽体积分数可以用气泡数密度(n)和气泡半径(RB)表示为

(8)

由式(7)和(8)可得

(9)

由式(7)和(9)，得相变速率R

(10)

根据气泡半径重建相变速率表达式

(11)

式(11)精确表达了空化中从液相到气相(气泡的生长或气化)的质量传递，单位体积的质量交换率与气相密度(ρv)、液相密度(ρl)和混合密度(ρ)相关。该式通常作为气泡破裂建模的第一近似值。方程右边采用绝对压力差来修正并被视为源汇相。在不考虑非凝结气体、质量传输和黏性阻尼时，气泡压力(pB)等于饱和蒸汽压力(pv)。SINGHAL等[13]提出的模型计入了紊流和不凝结气体的影响。控制方程为

(12)

式中：fv为气相质量分数；fg为不凝结气体质量分数；Γ为扩散系数。

当p≤pv，

(13)

当p>pv，

(14)

通过对湍流压力波动的局部值的估计来校正饱和压力

(15)

式中：Fvap=0.02，Fcond=0.01。

1.2 刚度阻尼计算方法

油膜的刚度阻尼是动压轴承的重要特性参数。转子在平衡位置受外载荷作用，对油膜产生位移扰动和速度扰动，轴承油膜力将发生变化。由摄动理论，将油膜力线性化处理后，油膜力的变化与转子扰动位移和扰动速度之间的关系如下：

(16)

(17)

采用动网格技术，令转子分别在水平方向和垂直方向移动微小位移ΔX、ΔY，并计算出有无扰动位移前后油膜在水平和垂直方向的油膜力的变化ΔFX、ΔFY，通过式(16)、(17)可以求得KXX、KYX、KXY和KYY。

(18)

(19)

2 计算域网格与边界条件

文中分析的油膜轴承是核电汽轮机动压轴承，相比抽油机、小型油膜轴承，具有更大的尺寸和相对更小的半径间隙。轴承直径D=800 mm，宽度B=680 mm，相对顶隙比ψ=0.15%，载荷为F=1 374.65 kN。入口直径d=160 mm，方形槽宽度为轴承宽度的0.7倍。润滑油为20#汽轮机油，分析中计入温黏效应。入口油压为0.2 MPa，出口油压为0.1 MPa。空化模型中空化压力采用文献[14]的推荐值29 158 Pa。

油膜流体域网格在ICEM CFD中采用六面体结构化网格划分，进口处采用O-BLOCK结构划分，采用Determinant 2×2×2标准检查全局网格质量大于0.6。采取50万、100万、200万、400万4种网格数量方案，如图1和图2所示为最大压力和承载力随网格数量变化曲线，可知当网格数量大于200万后，油膜最大压力以及油膜承载力变化较小，已经满足计算精度的要求。兼顾计算精度、计算时间，后续分析时网格数量统一取为200万。

图1 最大压-网格数量曲线

图2 承载力-网格数量曲线

3 结果与分析

3.1 流动模型和空化模型的比较讨论

表1 各转速工况下雷诺数

表2示出了轴承在转速为1 500 r/min、偏心率为0.7时，分别在层流和紊流realizable k-epsilon模型下的静态参数。可见，层流模型计算的数值与经验值偏离较远，紊流模型的结果与经验数值更为接近。同时，采用紊流模型计算的最大温度值、压力值与文献[2]同直径的大尺寸轴承试验得到的最大温度值82.8 ℃、压力值6.65 MPa也较为接近。因此，采用湍流模型比较适合所研究的轴承的油膜流动状态。

表2 不同流动模型分析结果

空化现象常存在于轴承油膜中，常用的数值空化模型有3种：Singhal et al.模型，Zwart-Gerber-Belamri模型和Schnerr and Sauer模型。文献[17]论述了各个空化模型的特点，对于旋转机械的空化现象模拟，该文献建议使用后2种空化模型。而在文中的计算分

析中，通过调整初始化方式和松弛因子，Zwart-Gerber-Belamri模型和Schnerr and Sauer模型收敛异常困难，epsilon和气相的残差很难降低至10-2量级，而Singhal et al.模型收敛精度高，各残差均在10-5量级以下。因此文中采用Singhal et al.模型模拟油膜中空化效应。

3.2 基于力平衡的油膜参数求解

为了得到额定转速1 500 r/min及给定外载荷1 374.65 kN下的轴承最小油膜厚度及相关参数，建立了不同偏心量的油膜轴承模型并进行试算，如表3和图3所示。当某偏心率在该转速下刚好满足外荷载等于承载力时，此时油膜偏心率即为所求偏心率。其他静态特性及相关参数在此偏心基础上进行计算求取。

表3 不同偏心下承载力计算

图3 承载力-偏心率曲线

试算结果表明，当偏心率为0.69时轴承的承载力为1 365.00 kN，与给定外载荷仅偏差约0.73%。此时计算的最小油膜厚度约为186 μm。图4所示为此时轴颈压力分布，最大压力约为6.35 MPa。图5所示为流体域中速度场的分布，最大线速度为62.74 m/s，入口方形槽中速度较低但存在旋流引起的涡。

图4 轴颈侧油膜压力分布

图5 速度场分布

图6所示为轴瓦温度分布云图，最大温度为81.52 ℃，在接近入口油槽处。图7所示为轴颈和轴瓦侧油膜的气穴体积分数云图，其中轴颈表面的气穴体积分数更高，空化效应较轴瓦侧更加明显。气穴体积分数最大处在接近入口油槽处，最大约为0.81。

图6 轴瓦侧温度分布云图

图7 轴瓦侧(左)和轴颈侧(右)气穴体积分数

3.3 转速波动下的轴承静特性

工程应用中，除了考察额定转速下的静特性外，也关注其他各重要转速如1 000、1 200、1 750 r/min时轴承的最小油膜厚度、温度、压力偏位角等静态特性参数。同样采用基于承载力和外载荷平衡的原则进行多转速工况下，多偏心尺寸的油膜轴承试算。表4给出了各转速下满足外载荷要求的最小油膜厚度值。由图8可知，在该转速范围内，转子转速越高轴承最小油膜厚度越大。转速过低时，如1 000 r/min时，最小油膜厚度仅为132 μm。

表4 定载荷各转速下的最小油膜厚度

图8 不同转速下的最小油膜厚度曲线

Fig 8 Minimum oil film thickness curve at different rotation speeds

图9示出了1 000～1 750 r/min转速范围内偏位角的变化趋势，转子转速越高轴承的偏位角越小，偏位角范围为40°～55°。图10示出了各转速下油膜中心周向的压力分布，转子转速越高，压力峰值越低，压力峰值的位置出现得越早。在各转速下压力峰值的范围为6～8 MPa。图11示出了各转速下油膜中心周向的温度分布，温度趋势均随角度增大而增大，在接近入口方形槽温度急速降低，最大温度在接近入口油槽处。

图9 不同转速下的偏位角曲线

图10 轴颈中心压力曲线

图11 轴瓦中心温度分布曲线

3.4 油膜动态特性求解

图12示出了各转速下油膜的刚度系数，图13示出了各转速下的阻尼系数。因计算刚度和阻尼采用相同的扰动速度和扰动位移，因此刚度和阻尼在随转速的变化趋势上相同。在1 000～1 750 r/min转速范围内，刚度系数的量级为109～1010，阻尼系数的量级为107～108。

图12 各转速下油膜刚度系数

图13 各转速下油膜阻尼系数

4 结论

(1)针对大尺寸动压轴承油膜流场，可根据临界雷诺数准则判断流动模态，而不应直接采用动压轴承设计手册的层流假设。针对文中研究对象在1 000、1 250、1 500、1 750 r/min时采用realizable k-epsilon模型计算更为准确。

(2)针对文中研究对象，在低于额定转速1 500 r/min以下工况时，润滑油温度由常温到较高温度转变时，流动模态有可能发生转换，此时轴承油膜的承载力和其他参数均有可能发生改变，易于导致轴承运行不稳，值得引起关注。

(3)在对动压轴承进行数值空化模拟收敛困难时，空化模拟中可采用Singhal et al.模型，该模型拥有更高的收敛精度和更快的收敛速度。

(4)计算得到了典型尺寸轴承在不同转速工况条件下油膜厚度规律、压力分布、温度分布、滑油流量等静态特性和轴承刚度、阻尼动态特性等数据，可为类似大型动压油膜轴承的设计和轴承的载荷分配提供参考。