以超临界二氧化碳为介质的径向箔片轴承静特性分析*

车国铚 杨启超 滕 斌 魏志国 高志成

(1.青岛科技大学机电工程学院 山东青岛 266061；2.中国船舶重工业集团公司第719研究所湖北武汉 430064；3.广东智空动力科技有限公司 广东佛山 528216)

超临界二氧化碳布雷顿发电循环是以超临界状态的CO2作为循环工质的热力循环，二氧化碳的临界压力为7.38 MPa,临界温度为31.1 ℃[1]，在靠近临界点的区域，S-CO2表现出密度大、黏度小的特点，这一物性特点使得占涡轮机输出功比重较大的压缩功减少，提高了S-CO2布雷顿循环的效率[2]。在同等条件(压力8～20 MPa和温度450～650 ℃)下，S-CO2布雷顿循环的效率在应用较多的蒸汽朗肯循环和氦气布雷顿循环效率之上[3]。除此之外，S-CO2具有更高的介质密度，减小了系统设备的尺寸，使得系统结构紧凑[4]，在经济性方面具有更好的潜力。因此超临界二氧化碳布雷顿循环在太阳能、火电、核电等能源领域具有良好的发展前景[5]。美国桑迪亚(Sandia)国家实验室[6]、韩国能源技术研究院(KIER)[7]、日本东京工业大学(Tokyo Institute of Technology，TIT)[8]等研究机构均率先对S-CO2布雷顿循环进行了理论及实验研究，并在此基础上开发了不同功率的循环系统。

压缩机是S-CO2布雷顿循环系统中一个关键设备，S-CO2循环工质通过压缩机提高至涡轮机的入口压力，最终高温高压的S-CO2流体在涡轮机中做功。为了最大化系统效率，压缩机转子的转速通常较高，可达到100 000 r/min[9],因此用于支撑高转速转子的轴承技术是超临界二氧化碳布雷顿循环压缩机的关键技术之一。美国Sandia国家实验室曾选用滚珠轴承进行轴向力的承载，但轴承的运行寿命较短，其滚珠轴承寿命有20～2 000 h不等，后期采用气体箔片轴承来进行长时间的测试[10]。FLEMING等[11]曾根据S-CO2循环的不同功率等级提出了S-CO2透平机械以及其中关键部件的推荐路线。众研究机构设计实验的超临界二氧化碳循环功率多在100～300 kW范围[12]，根据推荐路线，轴承推荐选择气体轴承。并且美国Sandia国家实验室、TIT等在实验测试样机的轴承选择上均选用以S-CO2直接润滑气体轴承，成功地实现稳定、长时间的转子力的承载，验证了气体轴承的技术可靠性[13]。

以上众研究机构在设计S-CO2压缩机时，压缩机的进口压力设计在7.6～9 MPa范围，进口温度设计在305～330 K范围。该进口温度及压力处在临界点附近，S-CO2物性在临界点附近呈现出剧烈的非线性变化。在利用S-CO2的特殊物性减少压缩机耗功、提高循环效率的同时，也应该考虑物性的剧烈变化对轴承性能以及转子稳定性的影响。传统动压径向气体轴承的研究主要以空气为工质，发展至今，在数值模拟、理论分析及试验研究方面都取得了较多的成果。2004年虞烈[14]通过引入弹性箔片静、动变形，联立求解气体润滑Reynolds方程，给出了以空气为润滑介质的箔片动压轴承的完全气弹润滑解。2006年杨利花等[15]曾建立气体轴承性能测试试验台，并对径向箔片气体轴承进行承载力和起飞转速的实验研究。2016年贾晨辉等[16]采用流体动力学软件，以动压气体轴承为研究对象，研究气体轴承瞬态非线性动力学行为，探讨轴承-转子系统的稳定性。2017年冯凯等人[17]建立了箔片径向气体轴承静态特性和温度特性的测试试验台，通过变载荷变转速方式对起飞转速、摩擦力矩和温度进行测量。但目前对于S-CO2润滑介质气体轴承的研究还较少，近临界区S-CO2密度、黏度的大幅变化，以及高压下S-CO2内部流动的湍流状态，这都将增加对预测、计算径向轴承性能的难度。另外，基于空气工质的气体轴承研究及设计方法是否适用的问题也是研究学者关注的方向。

2017年国内温建全[18]在连续性方程和动量方程的基础上推导层流状态下适用于S-CO2径向轴承的雷诺方程，并考虑湍流速度脉动的影响，以修正系数的方式加入雷诺方程中。计算时S-CO2物性采用差值方法获得，并用有限差分数值计算方法初步计算了径向气体轴承的静态性能，但未分析工况参数及结构参数对径向气体轴承静特性的影响。

2013年国外CONBOY[19]建立了一个弹性流体动力学模型，考虑润滑层内的流体速度场、流体动压和摩擦损失,计算了S-CO2润滑的推力轴承中的压力场，结果表明，计算域内所有离散点均为湍流区。除此之外，进行了不同转速以及不同轴承外半径下，空气止推轴承与S-CO2止推轴承的承载能力对比。

此外，2020年中国科学院大学BI等[20]为研究S-CO2动压径向气体轴承的动态特性，考虑了S-CO2湍流状态以及非线性热力学性质，在传统扰动法的基础上，引入扰动密度、扰动黏度、扰动雷诺数和扰动湍流系数，扩展了频率扰动法，用数值方法计算了包含所有动力学变量的频率相关刚度和阻尼系数。

本文作者考虑S-CO2的实际物性以及湍流状态的影响，运用有限差分法数值求解雷诺方程，重点分析近临界区S-CO2的物性突变对气膜、承载力、摩擦力矩等的影响；另外分析工况参数以及轴承结构参数对径向箔片气体轴承静特性的影响，为S-CO2布雷顿循环压缩机中径向气体轴承的研究提供理论依据。

1 箔片动压径向气体轴承结构

图1为文中研究的波箔型动压径向气体轴承结构示意图，轴承座里包含以特殊材料制成的耐磨、耐高温的弹性波箔片以及平箔片，平箔片和波箔片的一端固定在轴承座上，另一端可自由伸缩。轴颈在旋转过程中，带动润滑气体由间隙大端向间隙小端运动，形成动压效应，润滑气膜因此产生承载力。轴承在工作时，弹性箔片之间的变形作用以及其之间的库仑摩擦作用，可以吸收不平衡涡动能量[21]，即使在受到不稳定涡动时，气体轴承也能保持较好的稳定性。

图1中，O为轴承的中心；Ob为轴颈的中心；e为轴颈的中心相对于轴承中心之间的偏心距；R为轴颈的半径；Ω为轴颈的角速度；hmax为最大楔形间隙处的气膜厚度；hmin为最小楔形间隙处的气膜厚度。

2 数学模型及求解

2.1 变密度变黏度湍流雷诺方程

传统的动压气体轴承多以空气为工质进行研究，空气的物性变化较小，但近临界区的S-CO2物性变化剧烈，微小的温度压力变化将会引起密度与黏度的剧烈变化。除此之外，由于S-CO2密度大、黏度小，导致润滑气膜中会出现湍流状态，在气体润滑模型中必须加以考虑。因此，一般形式的可压缩气体雷诺方程便不再适用于S-CO2动压气体轴承。文中采用文献[18]中由连续性方程和动量守恒方程导出的变密度变黏度层流雷诺方程，在此基础上，将湍流的影响以周向和轴向的湍流润滑系数加入变密度变黏度层流雷诺方程中，其中kx为周向湍流润滑系数，kz为轴向湍流润滑系数，这是与雷诺数有关的非线性函数[22]，可得如下形式的雷诺方程：

(1)

(2)

(3)

式中:h为润滑气膜实际厚度(mm)；C为轴承间隙(mm)；pa为径向轴承外部环境压力(Pa);μa为轴承外部环境下的S-CO2动力黏度(Pa·s)；ρa为轴承外部环境下S-CO2密度(kg/m3)；L为轴承长度(mm)。

对式(1)量纲一化后得到：

(4)

2.2 气膜厚度方程

箔片气体轴承在工作时，不同于刚性气体轴承，箔片的弹性变形会影响轴承的气膜厚度大小，数值分析时需考虑箔片的弹性变形方程。文中所采用的弹性箔片结构如图2所示，模型中不考虑平箔的刚度，忽略弹性箔片之间的相互摩擦作用以及波箔与轴承座的摩擦，则箔片柔度[23]为

(5)

式中：pa为轴承外部环境压力(Pa);s为波箔单元的长度(mm);l为半波箔长度(mm)；Eb为波箔材料弹性模量(Pa)；νb为波箔片的泊松系数;tb为波箔片厚度(mm)。

一定偏心率下轴颈静平衡位置处的气膜厚度由2部分组成，一部分为箔片变形前的几何间隙，另一部分为箔片变形的变形间隙。量纲一气膜厚度表达式为

(6)

2.3 变密度变黏度湍流雷诺方程的数值求解

采用有限差分数值方法，将轴承的计算域离散化为差分网格，用有限个网格节点代替连续的求解域，润滑气膜展开如图3所示的m×n的离散网格。图中Δθ、Δz分别为周向和轴向的网格间距，Δθ=2π/m，Δz=2/n，i、j代表周向和轴向的网格节点编号。

基于有限差分法，采用计算精度较高的中心差商，将控制方程中压力在圆周方向和轴向方向的压力偏导数用相邻节点函数值的中心差商近似表示。

(7)

式中：i+1/2和i-1/2分别为i后半点处和前半点处；j+1/2和j-1/2为j后半点处和前半点处。

根据式(7)对式(4)进行离散化。

离散式(4)左边第一项：

(8)

离散式(4)左边第二项：

(9)

离散式(4)右边项：

(10)

将式(8)、(9)、(10)代入量纲一化方程(4)整理得到如下矩阵方程组：

Pi,j=

(11)

式中：Ai,j、Bi,j、Ci,j、Di,j、Ei,j、Fi,j为系数矩阵，各系数矩阵计算式如公式(12)所示。

(12)

由于轴承的两端(z=1,z=n+1)以及弹性箔片的固定端与轴承外部是相通的，且气膜未展开之前，该部分(θ=1，θ=m+1)是重合的，所以这几个区域的压力均为径向气体轴承的工作环境压力，故量纲一求解边界条件为

(13)

为了加快迭代计算过程的收敛速度，文中采用松弛迭代法进行计算，迭代格式为

(14)

通常数值计算的迭代过程中，常采用一定的相对收敛准则来判断迭代是否达到精度要求，达到精度要求后，便可终止迭代。文中相对收敛准则为

(15)

2.4 计算对象及方法

2.4.1 计算对象

文中动压径向轴承的结构参数以及润滑气体参数见表1。

表1 动压径向气体轴承相关参数

2.4.2 密度黏度处理方法

对于S-CO2的物性处理，文中采用美国国家标准与技术研究院的NIST Refprop 软件获取S-CO2流体热物性，并通过软件编程调用其物性。

2.4.3 计算流程

整体迭代计算流程如图4所示。

2.5 静态性能计算

通过耦合变密度变黏度湍流雷诺方程与箔片变形方程，求解获得润滑气膜的气膜压力分布，在此基础上进行静态特性的计算，其中包括承载力、偏位角、摩擦力矩等。

文中轴颈外载荷方向假定垂直向下，因此气膜的合力方向垂直向上，与外载荷相互平衡。假定偏心方向上的力为Wn，垂直偏心方向的力为Wt，则其计算公式为

(16)

(17)

式中：负号方向代表力垂直于轴颈并指向里。

轴承的量纲一承载力为

(18)

有量纲承载力为

(19)

偏位角满足如下关系

(20)

(21)

作用在轴颈表面的摩擦力矩其量纲一化形式为

(22)

则有量纲摩擦力矩为

(23)

3 数值结果及分析

3.1 程序验证

根据图4所示流程，用MATLAB 编写数值计算程序。为验证程序的正确性，计算刚性表面轴承不同工况下的承载力、偏位角，并与文献[18]结果对比。对比结果如表2所示，其中L/D为轴承的长径比；Λ为轴承数；ε为轴颈偏心率；Re为初始雷诺数。从表2可以看出文中计算结果与文献[18]结果非常接近，最大误差为2.9%，因此验证了文中所用计算程序具有一定的准确性。

表2 文中计算结果与文献[18]结果对比

3.2 近临界区物性及轴承静态特性计算结果

图5和图6所示分别为S-CO2的密度与黏度在不同压力下随温度的变化规律。可以看出，S-CO2密度与黏度随温度的升高均呈现减小的趋势，且越接近临界点时，变化越显著。当压力为7.6 MPa，温度从305 K升高1 K时，密度减小约38%，黏度减小约39%。而与压力的关系则相反，当温度一定时，S-CO2密度与黏度随压力的升高而增大。当温度为306 K,压力从7.6 MPa升高至8.0 MPa时，密度增大约73.8%，黏度增大约80%。气体轴承的静态特性与工作介质的黏度、压力及温度都密切相关，因此，研究物性变化对轴承静态性能的影响至关重要。

图7和图8所示分别为转速48 000 r/min、偏心率0.6条件下的S-CO2润滑气膜的厚度分布与气膜压力分布。由图7可以看出，气膜厚度在轴承轴向上呈现出两端小，中间大的分布；在圆周方向上呈现为先下降后增大的趋势。轴向中间气膜厚度大的原因是箔片未变形时此处压力较高，弹性箔片发生变形，产生额外的间隙空间，气膜厚度增大。

图8中，量纲一最大气膜压力为1.01，约为7.6 MPa,量纲一最小气膜压力为0.987，约为7.42 MPa,气膜压力峰值出现在圆周方向150°左右，气膜压力在圆周角度呈现先增后减的趋势。这是由于转子带动黏性气体运转的过程中，间隙由大变小，再由小变大。间隙大的地方气膜厚度大，压力较低。间隙小的地方气膜厚度小，压力较高。

气膜的厚度和压力在轴承楔形间隙中的分布与轴承的性能息息相关，且温度、压力引起的物性变化也最直接体现在润滑气膜的变化上。图9和图10所示分别为7.515 MPa下，气膜厚度与气膜压力轴向中截面随温度变化的曲线。从图9可以看出，随着温度的升高，气膜压力峰值逐渐降低。温度为307 K时，量纲一气膜压力峰值约为1.015，而温度为313 K时，量纲一气膜压力峰值约为1.013，压力峰值降低约0.015 MPa。这是由于随着温度的升高，S-CO2黏性减小，动压效应减弱，压力有所降低。图10中，最小气膜厚度同样随着温度升高而减小。温度为307 K时，量纲一气膜厚度约为0.75，而温度为313K时，量纲一气膜厚度约为0.65，降低了约0.1。这是由于黏性减小导致动压效应减弱的缘故，压力有所降低，箔片变形减小，从而产生的额外间隙减小。

图11和图12所示分别为转速48 000 r/min、偏心率0.6条件下的近临界区承载力、摩擦力矩随温度和压力的变化关系。从图11可以看出，在305～331 K的近临界温度范围内，随着温度的升高，径向轴承承载力先急剧降低，然后逐渐减小，最终趋于平缓。越靠近临界点时，温度的变化对承载力影响越显著。当压力为7.6 MPa时，温度由305 K升至307 K，承载力减小了约22.8%;压力为8.0 MPa时，温度由305 K升至307 K,承载力减小了约19.3%。同样地，压力的变化对承载力的影响也非常显著，承载力随着压力的提高而增大。如图11所示，温度为307 K时，压力从7.6 MPa升高至8.0 MPa，承载力大幅提高约32%；温度为309 K时，压力从8.0 MPa升高至8.4 MPa，承载力提高约25%。相较于传统滚珠轴承，动压气体轴承本身承载力较小，若改变工作条件，径向轴承承载力降低会不利于压缩机转子的稳定运行。因此在设计S-CO2压缩机时，需设计合适的进出口温度及压力。

由图12可以看出，压力一定时，随着温度升高，径向轴承摩擦力矩先急剧降低，然后逐渐减小，最终趋于平缓。靠近临界点的区域摩擦力矩有大幅减小现象，这是由于黏度的大幅减小，润滑气膜与轴颈之间的黏性力减小导致的。压力为7.6 MPa时，当温度由305 K升高至307 K时，摩擦力矩减小约24%；压力为8.0 MPa时，当温度由305 K升高至309 K时，摩擦力矩减小约33.3%。同样地，压力的变化对摩擦力矩的影响也非常显著，摩擦力矩随着压力的提高而增大。307 K时，压力从7.6 MPa升高至8.0 MPa，摩擦力矩增大约28.8%。309 K时，压力从8.0 MPa升高至8.4 MPa时，摩擦力矩增大约32.5%。摩擦力矩的大幅提高一方面使得功率损耗增大,另一方面会使得轴颈高速运转过程中所产生的热量增多，以致于导致过高的温升使物性发生剧烈变化。

3.3 偏心率对径向轴承静态性能的影响

图13—15描述了T=314 K、p=7.515 MPa条件下，不同转速时承载力、偏位角、摩擦力矩随轴颈偏心率的变化。由图13可看出，同一转速下，径向气体轴承的承载力随着偏心率增大而增大；当转速提高时，承载力逐渐变大；在小偏心率时，各转速下的轴承承载力差距较小，随着偏心率的提高，差距逐渐变大。在转子运转过程中，当偏心率增大时，轴颈与弹性箔片之间的楔形间隙变小，在这样的情况下，轴颈挤压气膜的作用就显得愈加强烈，由此一来，楔形间隙中的气膜压力变高，径向轴承的承载力得以提高。

从图14可看出，随着偏心率的增大，偏位角逐渐减小；同一偏心率条件下，当转速增大时，偏位角呈减小的趋势。由于偏位角能表征系统的稳定性，偏位角越小，稳定性相对越好。因此转速的增大，在一定程度上有利于提高系统的稳定性。

从图15可看出，随着偏心率的增大，摩擦力矩呈现出加速增长的趋势；且转速越大，摩擦力矩越大，这是因为当转子的转速增大时，轴颈与S-CO2润滑气膜之间的摩擦加剧，摩擦功耗增大。另外从计算结果可以看出，相同偏心率变化时，转速越大时摩擦力矩相应变化也越大。

3.4 轴承数对静态性能的影响

轴承数是表征气体轴承工作状态的重要参数之一，由于超临界二氧化碳布雷顿循环压缩机的径向轴承工作环境压力高，相较于传统的空气动压轴承，轴承数的数值较小。图16—18描述的是T=314 K、p=7.515 MPa、长径比为1条件下，不同偏心率时承载力、偏位角以及摩擦力矩随轴承数的变化。从图16可看出，随着轴承数的增大，径向气体轴承的承载力呈现出减速增长的趋势；在小轴承数条件下，承载力较小，且不同偏心率时的承载力差值较小，而随着轴承数的增大，其差值不断增大，增大一定程度后趋于不变。

如图17所示，随着轴承数的增大，径向气体轴承的偏位角呈现出减速下降的趋势。由于偏位角能表征系统的稳定性，偏位角越小，稳定性相对越好。因此轴承数的增大，有利于保持系统的稳定。

由图18可以看出，随着轴承数的增大，径向气体轴承的摩擦力矩呈现出增长趋势，其中小偏心率下近乎线性增长，随着偏心率的增大，增长趋势变为减速增长。偏心率为0.6时，轴承数从0.001增大到0.01时，摩擦力矩增大约2.75倍。

3.5 长径比对静态性能的影响

轴承的长径比(L/D)是气体轴承的重要结构参数，文中研究长径比对径向轴承静态性能的影响。图19和图20所示分别为T=314 K、p=7.515 MPa、Λ=0.01条件下,不同偏心率时承载力和偏位角随长径比的变化情况。

如图19所示，随着长径比的增大，承载力呈现线性增长趋势；长径比为0.5时，即使在大偏心率情况下，承载力也较小，当偏心率为0.6时，长径比为2时的承载力为长径比为1条件的2.4倍，拥有较大的承载力。但长径比较大时轴承的计算应该考虑温度场的影响。为保证径向气体轴承能够对高速转子进行稳定支撑，在设计S-CO2压缩机径向气体轴承时应选择合适的长径比，以保证合适的承载力。

由图20可以看出，随着长径比的增大，偏位角呈现减速减小的趋势。在长径比从0.5增大到1时，偏位角急剧下降，如偏心率为0.6时，偏位角减小了约17.6°；然而长径比从1增大到2时，偏位角变化却很小，如偏心率为0.6情况下偏位角仅减小3.5°。

4 结论

考虑S-CO2实际物性，数值求解箔片动压轴承在不同温度、压力工作环境下的承载力和摩擦力矩等静态特性。主要结论如下：

(1)近临界区密度、黏度的剧烈变化会引起承载力和摩擦力矩的大幅变化，且越靠近临界点，变化程度越大；同一压力下，承载力和摩擦力矩随着温度的升高而减小，同一温度下，承载力和摩擦力矩随着压力的升高而增大。

(2)随着偏心率的增大，S-CO2轴承的承载力和摩擦力矩都有不同程度的增大，而偏位角减小；随轴承数的增大，S-CO2轴承的承载力和摩擦力矩增大，而偏位角减小；增大径向轴承的长径比，S-CO2轴承的承载力增大，而偏位角减小。