S-CO₂介质止推箔片气体动压轴承特性研究

Investigations on Characteristics of Thrust Gas Foil Bearings Lubricated by S-CO₂

摘要

针对超临界二氧化碳（Supercritical Carbon Dioxide,S-CO₂）润滑波箔型止推箔片气体动压轴承，通过反向传播（Back Propagation,BP）神经网络算法提出S-CO₂的物性模型，并考虑轴承工作时的非理想气体效应，提出考虑湍流效应气体润滑模型、箔片结构力学模型和气膜平均温升计算方法，对止推箔片气体动压轴承的静动态特性进行研究，并分析不同结构参数对箔片气体轴承静动态特性的影响规律.结果表明，本文提出的物性模型准确度高，相关系数高达99.997%.以S-CO₂为润滑介质的止推箔片气体动压轴承具有更高的承载力，且在适当范围内减小最小初始气膜厚度或增加膜厚比可以提高轴承的承载力.以S-CO₂为介质的止推箔片气体动压轴承的动态刚度系数和动态阻尼系数均远高于常温常压空气介质下的止推箔片气体动压轴承.随着最小初始气膜厚度减小，轴承的动态刚度系数和动态阻尼系数均迅速增加.

Abstract

For thrust gas foil bearings integrating Supercritical Carbon Dioxide （S-CO₂） as the working medium， a Back Propagation （BP） neural network algorithm is employed to propose a physical property model for S-CO₂. To account for the non-idea gas behavior in the bearing， models encompassing gas lubrication with turbulence effects， foil structural mechanics， and the calculation of the average gas film temperature are introduced. Static and dynamic characteristics of the thrust gas foil bearing are studied and contrasted with various lubricant media. The impact of diverse structural parameters on the static and dynamic attributes of the gas foil bearing is analyzed. Results indicate that the physical property model presented in this paper attains high accuracy， boasting a correlation coefficient of 99.997%. The thrust gas foil bearing lubricated with S-CO₂ exhibits enhanced load-bearing capacity， with potential for further improvement by adjusting the minimum initial gas film thickness or increasing the film thickness ratio within a suitable range. The dynamic stiffness and damping coefficient of the thrust gas foil bearing utilizing S-CO₂ significantly surpass those employing air as the medium， underscoring its superior dynamic characteristics. Furthermore， a reduction in the minimum initial film thickness leads to a rapid increase in the dynamic stiffness coefficients and damping coefficients of the bearing.

Graphical abstract

关键词

止推箔片气体动压轴承 / 超临界二氧化碳（S-CO₂） / 反向传播神经网络 / 静态特性 / 动态特性

二氧化碳具有环保、无毒、物性稳定、成本低廉和传热特性好等优点^［1］，并且在超临界区域比热容较大、导热性能好、动力黏度低^［2-4］，因此，二氧化碳介质被广泛应用于化学工业、余热利用^［4］、太阳即热发电、核电气冷堆、动力机械等领域^［5］. 相较于传统的蒸汽动力循环系统，以S-CO₂为工作介质的布雷顿循环系统具有更高的效率和更紧凑的动力系统结构^［5-6］，因此，S-CO₂布雷顿循环系统具有更良好的经济性和更广阔的工程应用前景^［7-8］. 目前，以S-CO₂为介质的旋转机械正在向超高速、超高温和高能量密度等目标持续迈进^［9］，传统的油轴承和滚动轴承难以满足上述条件.箔片气体动压轴承因其具有无油、免维护、高转速、耐高温高压、结构简单紧凑等特点，被视为S-CO₂旋转机械转子系统的理想支承部件^［10］.

相较于传统的油轴承和滚动轴承，箔片气体动压轴承无需供油或供气装置，采用箔片气体动压轴承支承转子可以达到缩小旋转机械的体积、减小质量、降低成本等目的^［11］；箔片气体轴承具有更高的可靠性、稳定性以及耐高温高压等优点^［12-13］.由箔片气体轴承支承的转子具有超高转速、高效率、长寿命、免维护、强环境适应性等优点^［14］，被广泛应用于空气循环机、压缩机、微型燃气轮机等高速旋转机械中^［15-17］ .尽管箔片气体轴承已得到成熟应用，但针对特殊介质下的箔片气体轴承的理论研究仍然不够完善.

为研究基于S-CO₂介质的箔片气体轴承的工作特性， Xu等^［18］提出了一种考虑实际气体效应和湍流流体的静压止推箔片轴承三维热流体动力分析的高级分析工具，该工具可以计算高压二氧化碳气膜的压力分布、温度分布、承载能力、气膜厚度、功率损失等性能参数. Qin等^［19］开发了一种用于止推箔片气体动压轴承三维弹流动力学模拟的求解器和一种新的求解轴承内部结构变形的求解器，结果表明该求解器可以准确地模拟在密集气体中运行的止推箔片气体动压轴承工作情况. Kim^［20］提出了考虑实际气体效应和气膜内部湍流流动的箔片气体径向轴承三维热流体力学分析的高级设计/分析工具，通过该工具模拟了不同环境条件下轴承性能，但是模拟结果低于公开文献中的测量数据. Qin等^［21］建立了流体相的准三维流固模型，研究了S-CO₂循环作用下止推箔片气体动压轴承的流体力学性能，并对流体力学求解器Eilmer进行了扩展，发现离心惯性力对轴承性能具有重要影响，为轴承性能的潜在影响因素提供了新的见解. 温建全^［5］针对S-CO₂闭式布雷顿循环的支承系统，将S-CO₂视为非理想气体并采用线性插值拟合的方法获取其物理性质，并推导出箔片气体径向和止推轴承的变密度变黏度的湍流雷诺方程，求解了S-CO₂介质的轴承静动态特性. 朱鹏程等^［14］针对搭接顶箔式波箔型气体动压轴承，耦合湍流效应修正雷诺方程、力学分析模型和能量方程，并进行了与空气介质对比分析以及不同轴承结构参数对径向轴承静态特性的影响规律.

但是，上述研究缺乏具体数学模型以准确地反映S-CO₂物性特征变化规律及动态分析模型，在一定程度上制约了箔片气体轴承在S-CO₂布雷顿循环系统中的应用.因此，针对止推箔片气体动压轴承，本文基于BP神经网络算法建立S-CO₂物性模型，并引入S-CO₂介质作为止推箔片气体动压轴承的工作介质，在计算过程中考虑S-CO₂介质的实际气体效应，提出考虑湍流效应的变黏度变密度气体润滑模型、箔片结构力学模型和平均温升计算方法，并建立了轴承静动态弹流耦合数值计算方法，分析了 S-CO₂介质下和常温常压空气介质下的止推箔片气体动压轴承特性对比以及不同轴承结构参数对轴承特性的影响规律.

1 理论模型

1.1 S-CO₂物性模型

与常温常压（0.1 MPa，300 K）下空气介质不同，S-CO₂在不同压强和温度下的密度、黏度呈强烈的非线性变化，因此无法使用理想气体状态方程来准确描述S-CO₂的物理性质. 而BP神经网络可以通过算法学习不同的数据模式和特征，处理数据间的非线性关系，从而可以作为S-CO₂物性模型的训练工具.

本文研究的止推箔片气体动压轴承温度为314~500 K，压力为7~8 MPa，建立的模型使用环境条件参数作为输入特征，训练数据集包含70 125个样本，其中70%用于训练网络，15%用于验证模型，15%用于测试模型.

图1为基于BP神经网络训练的S-CO₂物性模型结构图，该模型由输入层I、隐藏层H和输出层O组成. 其中，输入层由温度T和压力P组成，输出层由密度\rho、黏度\mu和比热容C_{\mathrm{P}}组成，隐藏层的层数Y由经验公式计算得到，如式（1）所示.

Y=\sqrt[]{m+n}+b（1）

式中：m为输入层的节点数量；n为输出层的节点数量；b为1~10中任意整数. 本文根据输入层节点和输出层节点数量，通过式（1）计算选取隐藏层的层数Y为10，形成2-10-3神经网络模型结构.

在BP神经网络训练过程中，温度和压力信息作为输入信息进入网络，依次传递到隐藏层和输出层. 每个神经元接受上一层所有神经元的输出作为输入信号进行加权求和，并加上一个偏置项，如式（2）所示.

f={\displaystyle \sum _{i=\mathrm{1}}^n}{\omega }_{ij}{x}_i+b_j（2）

式中：x_i为输入层的第i个神经元的输出；{\omega }_{ij}为连接输入层的第i个神经元和隐藏层的第j个神经元之间的连接权值；b_j为隐藏层的第j个神经元的偏置系数；f为激活函数.

在物性模型训练过程中，将训练数据集输入网络中，并计算输出结果与期望结果之间的误差，使用误差反向传播算法，将网络输出结果与真实结果进行比较，根据损失函数获得误差值，本文采用均方误差作为损失函数，如式（3）所示. 沿着网络反向传播误差，计算每个神经元对误差的贡献，根据每个神经元的贡献更新每个神经元之间的权重，使误差逐渐减小，达到期望的训练效果为止.

E=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{j=\mathrm{1}}^m}{\left(y_j-t_j\right)}^{\mathrm{2}}（3）

式中：E为损失函数；y_j为BP神经网络的预测值；t_j为真实值；m为神经元个数.

1.2 箔片结构力学模型

止推箔片气体动压轴承是一种具有柔性表面并以环境气体作为润滑介质的自作用式动压轴承，利用流体动压效应将周围环境气体吸入楔形间隙中，在推力盘表面和顶箔之间产生高压气膜以支承转子系统的轴向力.如图2所示，轴承主要由顶箔、波箔和轴承壳3部分组成.波箔的一端通过焊接或其他方式固定在轴承壳上，另一端自由，顶箔支承在波箔上，其中顶箔表面由倾斜表面和水平表面组成，与推力盘表面形成楔形收敛间隙.

本文采用Feng等^［12］提出的Link-Spring模型计算轴承箔片结构的等效刚度，该模型将波箔上的每个波形凸起单元简化为2个刚性连杆和水平放置的等效弹簧，如图3所示. 刚性连杆一端与相邻刚性连杆通过铰链连接，可以自由转动，另一端与弹簧连接，相邻弹簧通过刚性连杆连接.

在建立箔片结构的力学模型中，考虑了波箔的弹性、顶箔与波箔之间的摩擦、波箔与轴承壳之间的摩擦、波箔上波形凸起之间的相互作用等因素. 计算时，将顶箔和波箔划分成相同的网格分布，将每个波箔单元的等效刚度K_{\mathrm{V}}与相应的顶箔单元的等效刚度K_{\mathrm{t}\mathrm{o}\mathrm{p}}相加，即可得到整个箔片结构的等效刚度K_{\mathrm{f}}.

K_{\mathrm{f}}=K_{\mathrm{V}}+K_{\mathrm{t}\mathrm{o}\mathrm{p}}（4）

箔片结构的变形\mathrm{\Delta }H可通过直接刚度法进行计算，如式（5）所示.

F=K_{\mathrm{f}}\mathrm{\Delta }H（5）

式中：F表示由气膜压力产生的载荷.

1.3 静态分析模型

在止推箔片气体动压轴承的超临界工作条件下，S-CO₂介质气体由于具有较大的密度、较小的黏度和较高的工作转速（≥10 000 r/min），其气体流动雷诺数R{e}_{\mathrm{h}}较大（>10 000），气流处于湍流状态，因此基于层流假设推导的雷诺方程无法适用于描述S-CO₂介质下止推箔片气体动压轴承的流体特性.

本文在止推箔片气体动压轴承的流体润滑模型中引入湍流系数. 为获取流体湍流系数，Ng等^［22］将湍流剪切流动的“壁面定律”应用于湍流润滑的分析，假定湍流流体的各向同性特征，建立控制微分方程的线性化形式，获得流体湍流系数的倒数G_x和G_y曲线. 为简化雷诺方程的求解，本文通过曲线拟合获得流体湍流系数与雷诺数之间的关系.

K_{\theta }=\mathrm{12}+\mathrm{0.004}\mathrm{6}R{e}_{\mathrm{h}}^{\mathrm{0.97}}（6）

K_R=\mathrm{12}+\mathrm{0.002}\mathrm{5}R{e}_{\mathrm{h}}^{\mathrm{1.02}}（7）

R{e}_{\mathrm{h}}=\frac{\overline{\rho }\overline{H}}{\overline{\mu }}\frac{{\rho }_{\mathrm{0}}R\omega h_{\mathrm{2}}}{{\mu }_{\mathrm{0}}}（8）

式中：K_{\theta }和K_R分别为周向和径向的湍流修正系数；\overline{\mu }=\frac{\mu }{{\mu }_{\mathrm{0}}}；\overline{\rho }=\frac{\rho }{{\rho }_{\mathrm{0}}}；\overline{H}=\frac{H}{h_{\mathrm{2}}}；R为半径方向的有量纲长度；h_{\mathrm{2}}为最小初始气膜厚度；{\rho }_{\mathrm{0}}和{\mu }_{\mathrm{0}}分别为环境温度与压力下的密度和黏度；\omega为转子的角速度.

本文在流体润滑模型中考虑了湍流效应，在仅考虑层流的可压缩气体雷诺方程引入了湍流系数修正K_{\theta }和K_R，并进行无量纲化，得到变黏度变密度的湍流雷诺方程，如式（9）所示.

\frac{\mathrm{1}}{\overline{R}}\frac{\partial }{\partial \overline{R}}\left(\overline{R}\frac{\overline{\rho }{\overline{H}}^{\mathrm{3}}}{K_R}\frac{\mathrm{12}}{\overline{\mu }}\frac{\partial \overline{P}}{\partial \overline{R}}\right)+\frac{\mathrm{1}}{{\overline{R}}^{\mathrm{2}}}\frac{\partial }{\partial \theta }\left(\frac{\overline{\rho }{\overline{H}}^{\mathrm{3}}}{K_{\theta }}\frac{\mathrm{12}}{\overline{\mu }}\frac{\partial \overline{P}}{\partial \theta }\right)=\Lambda \frac{\partial }{\partial \theta }\left(\overline{\rho }\overline{H}\right)+\mathrm{2}\Lambda \frac{\partial }{\partial t}\left(\overline{\rho }\overline{H}\right)（9）

式中：\overline{P}=\frac{P}{p_{\mathrm{a}}}；\overline{R}=\frac{R}{R_{\mathrm{o}}}；\Lambda =\frac{\mathrm{6}{\mu }_{\mathrm{0}}\omega {R_{\mathrm{o}}}^{\mathrm{2}}}{p_{\mathrm{a}}{h_{\mathrm{2}}}^{\mathrm{2}}}；\theta为圆周方向的角度；p_{\mathrm{a}}为环境压力；R_{\mathrm{o}}为箔片外部半径；t为时间，在静态分析时，不考虑气膜随时间的变化，因此可以忽略含时间t的计算项.

由于S-CO₂的密度和黏度均为温度和压力的强非线性函数，因此在求解雷诺方程时无法将密度直接替换为压力，将黏度视为常数. 本文在求解雷诺方程时采用1.1节中基于BP神经网络算法训练的 S-CO₂物性模型获取密度和黏度.

通过式（10）和式（11）可求得止推轴承的轴向承载能力W和摩擦力矩M_{\mathrm{f}}.

W=p_{\mathrm{a}}{R}_{\mathrm{o}}^{\mathrm{2}}{\int }_{\mathrm{0}}^{\beta }{\int }_{R_{\mathrm{i}}/R_{\mathrm{o}}}^{\mathrm{1}}\left(\overline{P}-\mathrm{1}\right)\overline{R}\mathrm{d}\overline{R}\mathrm{d}\theta（10）

M_{\mathrm{f}}=p_{\mathrm{a}}{R}_{\mathrm{o}}^{\mathrm{2}}{h}_{\mathrm{2}}{\int }_{\mathrm{0}}^{\beta }{\int }_{R_{\mathrm{i}}/R_{\mathrm{o}}}^{\mathrm{1}}\left(\frac{\overline{R}\overline{H}}{\mathrm{2}}\frac{\partial \overline{P}}{\partial \theta }+\frac{\Lambda {\overline{R}}^{\mathrm{3}}}{\mathrm{6}\overline{H}}\right)\mathrm{d}\overline{R}\mathrm{d}\theta（11）

式中：\beta为箔片展角；R_{\mathrm{i}}为箔片内部半径.

1.4 温度计算模型

诸多学者采用气膜能量方程求解止推箔片气体动压轴承的温度分布，但是求解能量方程需要复杂的迭代计算，从而增加整个迭代过程的时间成本，也容易造成数据发散，因此，本文根据润滑膜的热平衡条件求解气膜的平均温升^［23］：

k{W}_{\mathrm{f}}=C_{\mathrm{P}}\rho \mathrm{\Delta }T{Q}_{\mathrm{o}\mathrm{u}\mathrm{t}}（12）

式中：k为半经验系数，取0.6；W_{\mathrm{f}}为气膜平均摩擦功率；C_{\mathrm{P}}为定压比热容；\mathrm{\Delta }T为气膜平均温升；Q_{\mathrm{o}\mathrm{u}\mathrm{t}}为气体端面泄流量.

黏性摩擦力f计算公式为：

f={\int }_{\mathrm{0}}^{\beta }{\int }_{R_{\mathrm{i}}}^{R_{\mathrm{o}}}\mu \left(\frac{\partial U}{\partial H}+\frac{\partial U}{\partial R}+\frac{\partial U}{\partial \theta }\right)\mathrm{d}\theta \mathrm{d}R（13）

式中：U为线速度.

由于在膜厚方向的速度梯度远大于周向和径向的速度梯度，因此黏性摩擦力计算式（13）可以简化为：

f={\int }_{\mathrm{0}}^{\beta }{\int }_{R_{\mathrm{i}}}^{R_{\mathrm{o}}}\left(\mu \frac{\partial U}{\partial H}\right)\mathrm{d}\theta \mathrm{d}R（14）

气膜平均摩擦功率W_{\mathrm{f}}的表达式为：

W_{\mathrm{f}}=fU（15）

止推箔片气体动压轴承在实际工作中压力分布会呈现中间高内外径两侧低的情况，因此扇形区域内润滑气体受到高压作用会向两侧流出，即端面泄流量.止推箔片气体动压轴承内部流体流动示意图如图4所示，止推轴承根据质量守恒定理，可得端面泄流量Q_{\mathrm{o}\mathrm{u}\mathrm{t}}的表达式为：

Q_{\mathrm{o}\mathrm{u}\mathrm{t}}=Q_{\mathrm{s}\mathrm{u}\mathrm{c}}-Q_{\mathrm{r}\mathrm{e}\mathrm{c}}（16）

式中：Q_{\mathrm{s}\mathrm{u}\mathrm{c}}为从周围环境吸入的流量；Q_{\mathrm{r}\mathrm{e}\mathrm{c}}为气膜内部循环的流量. Q_{\mathrm{s}\mathrm{u}\mathrm{c}}和Q_{\mathrm{r}\mathrm{e}\mathrm{c}}可根据式（17）和式（18）计算求得^［24-25］.

\begin{array}{l}{Q}_{\mathrm{s}\mathrm{u}\mathrm{c}}=\omega R_{\mathrm{o}}^{\mathrm{2}}{h}_{\mathrm{2}}{\int }_{\mathrm{0}}^{\beta }{\int }_{\mathrm{0}}^{\overline{H}}\overline{V_{\mathrm{r}}}\left|{}_{\overline{R}=\overline{R_{\mathrm{i}}}}\right.\overline{R_{\mathrm{i}}}\overline{H}\mathrm{d}\theta \mathrm{d}\overline{H}+\\ \omega R_{\mathrm{o}}^{\mathrm{2}}{h}_{\mathrm{2}}{\int }_{\mathrm{0}}^{\beta }{\int }_{\mathrm{0}}^{\overline{H}}\overline{V_{\mathrm{r}}}\left|{}_{\overline{R}=\overline{R_{\mathrm{o}}}}\right.\overline{R_{\mathrm{o}}}\overline{H}\mathrm{d}\theta \mathrm{d}\overline{H}\end{array}（17）

Q_{\mathrm{r}\mathrm{e}\mathrm{c}}=\omega R_{\mathrm{o}}^{\mathrm{2}}{h}_{\mathrm{2}}{\int }_{R_{\mathrm{i}}/R_{\mathrm{o}}}^{\mathrm{1}}{\int }_{\mathrm{0}}^{\overline{H}}\overline{V_{\theta }}\left|{}_{\theta =\beta }\right.\overline{H}\mathrm{d}\overline{R}\mathrm{d}\overline{H}（18）

式中：\overline{V_{\mathrm{r}}}为润滑气体无量纲径向流速；\overline{V_{\theta }}为润滑气体无量纲周向流速；\overline{R_{\mathrm{i}}}为无量纲径向位置尺寸.

因此，平均温升\mathrm{\Delta }T的计算式如下：

\mathrm{\Delta }T=\frac{k{W}_{\mathrm{f}}}{C_{\mathrm{P}}\rho Q_{\mathrm{o}\mathrm{u}\mathrm{t}}}（19）

1.5 动态分析模型

为获得超临界二氧化碳介质下的止推箔片气体动压轴承的动态特性，使用扰动法对非定常状态下的雷诺方程进行处理. 当止推箔片气体动压轴承在正常工况下工作时，会达到平衡位置，该平衡位置即为式（9）求解得到的静平衡位置，此时对推力盘施加一个轴向方向的位移扰动量\mathrm{\Delta }Z和速度扰动量\mathrm{\Delta }\dot{Z}，将无量纲气膜压力\overline{P}、气膜厚度\overline{H}和箔片变形量\mathrm{\Delta }\overline{H}在静平衡位置进行一阶泰勒展开：

\begin{array}{l}\overline{P}={\overline{P}}_{\mathrm{0}}+{\overline{P}}_Z\mathrm{\Delta }\overline{Z}+{\overline{P}}_{\dot{Z}}\mathrm{\Delta }\overline{\dot{Z}}\\ \overline{H}={\overline{H}}_{\mathrm{0}}+{\overline{H}}_Z\mathrm{\Delta }\overline{Z}+{\overline{H}}_{\dot{Z}}\mathrm{\Delta }\overline{\dot{Z}}\\ \mathrm{\Delta }\overline{H}=\mathrm{\Delta }{\overline{H}}_{\mathrm{0}}+\mathrm{\Delta }{\overline{H}}_Z\mathrm{\Delta }\overline{Z}+\mathrm{\Delta }{\overline{H}}_{\dot{Z}}\mathrm{\Delta }\overline{\dot{Z}}\end{array}（20）

将式（20）代入式（9），略去高阶项，可获得静平衡位置的稳态扰动方程、位移扰动方程和速度扰动方程.

稳态扰动方程为：

\begin{array}{l}\frac{\mathrm{1}}{\overline{R}}\frac{\partial }{\partial \overline{R}}\left(\overline{R}\frac{\overline{{\rho }_{\mathrm{0}}}{\overline{H_{\mathrm{0}}}}^{\mathrm{3}}}{K_{R\mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_{\mathrm{0}}}}{\partial \overline{R}}\right)+\\ \frac{\mathrm{1}}{{\overline{R}}^{\mathrm{2}}}\frac{\partial }{\partial \theta }\left(\frac{\overline{{\rho }_{\mathrm{0}}}{\overline{H_{\mathrm{0}}}}^{\mathrm{3}}}{K_{\theta \mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_{\mathrm{0}}}}{\partial \theta }\right)=\Lambda \frac{\partial }{\partial \theta }\left(\overline{{\rho }_{\mathrm{0}}}\overline{H_{\mathrm{0}}}\right)\end{array}（21）

位移扰动方程为：

\frac{\mathrm{1}}{\overline{R}}\frac{\partial }{\partial \overline{R}}\left(\overline{R}\frac{\overline{{\rho }_{\mathrm{0}}}{\overline{H_{\mathrm{0}}}}^{\mathrm{3}}}{K_{R\mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_Z}}{\partial \overline{R}}+\mathrm{3}\overline{R}\overline{H_Z}\frac{\overline{{\rho }_{\mathrm{0}}}{\overline{H_{\mathrm{0}}}}^{\mathrm{2}}}{K_{R\mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_{\mathrm{0}}}}{\partial \overline{R}}\right)+

\frac{\mathrm{1}}{{\overline{R}}^{\mathrm{2}}}\frac{\partial }{\partial \theta }\left(\frac{\overline{{\rho }_{\mathrm{0}}}{\overline{H_{\mathrm{0}}}}^{\mathrm{3}}}{K_{\theta \mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_Z}}{\partial \theta }+\mathrm{3}\overline{H_Z}\frac{\overline{{\rho }_{\mathrm{0}}}{\overline{H_{\mathrm{0}}}}^{\mathrm{2}}}{K_{\theta \mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_{\mathrm{0}}}}{\partial \theta }\right)=

\Lambda \frac{\partial }{\partial \theta }\left(\overline{{\rho }_{\mathrm{0}}}\overline{H_Z}\right)-\mathrm{2}\Lambda \gamma \overline{{\rho }_{\mathrm{0}}}\overline{H_Z}（22）

速度扰动方程为：

\begin{array}{l}\frac{\mathrm{1}}{\overline{R}}\frac{\partial }{\partial \overline{R}}\left(\overline{R}\frac{\overline{{\rho }_{\mathrm{0}}}{{\overline{H}}_{\mathrm{0}}}^{\mathrm{3}}}{K_{R\mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_{\dot{Z}}}}{\partial \overline{R}}+\mathrm{3}\overline{R}\overline{H_{\dot{Z}}}\frac{\overline{{\rho }_{\mathrm{0}}}{{\overline{H}}_{\mathrm{0}}}^{\mathrm{2}}}{K_{R\mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_{\mathrm{0}}}}{\partial \overline{R}}\right)+\\ \frac{\mathrm{1}}{{\overline{R}}^{\mathrm{2}}}\frac{\partial }{\partial \theta }\left(\frac{\overline{{\rho }_{\mathrm{0}}}{{\overline{H}}_{\mathrm{0}}}^{\mathrm{3}}}{K_{\theta \mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_{\dot{Z}}}}{\partial \theta }+\mathrm{3}\overline{H_{\dot{Z}}}\frac{\overline{{\rho }_{\mathrm{0}}}{{\overline{H}}_{\mathrm{0}}}^{\mathrm{2}}}{K_{\theta \mathrm{0}}}\frac{\mathrm{12}}{\overline{{\mu }_{\mathrm{0}}}}\frac{\partial \overline{P_{\mathrm{0}}}}{\partial \theta }\right)=\end{array}\Lambda \frac{\partial }{\partial \theta }\left(\overline{{\rho }_{\mathrm{0}}}\overline{H_{\dot{Z}}}\right)-\mathrm{2}\Lambda \gamma \overline{{\rho }_{\mathrm{0}}}\overline{H_{\dot{Z}}}（23）

式中：下标0表示静平衡位置；\overline{P_Z}和\overline{P_{\dot{Z}}}分别为位移扰动和速度扰动下的无量纲气膜压力；\overline{H_Z}和\overline{H_{\dot{Z}}}分别为位移扰动和速度扰动下的无量纲气膜厚度；\gamma为频率比，由激振频率\nu和转子角速度\omega相除得到.

\gamma =\frac{\nu }{\omega }（24）

依次对式（21）、式（22）和式（23）进行求解，得到\overline{P_Z}和\overline{P_{\dot{Z}}}，在求解域内进行积分，可以得到止推箔片气体动压轴承的轴向动态刚度系数K_Z和动态阻尼系数D_Z，计算公式如下：

K_Z=-\frac{p_{\mathrm{a}}{R}_{\mathrm{o}}^{\mathrm{2}}}{h_{\mathrm{2}}}{\int }_{\mathrm{0}}^{\beta }{\int }_{R_{\mathrm{i}}/R_{\mathrm{o}}}^{\mathrm{1}}\overline{P_Z}\overline{R}\mathrm{d}\overline{R}\mathrm{d}\theta（25）

D_Z=-\frac{p_{\mathrm{a}}{R}_{\mathrm{o}}^{\mathrm{2}}}{h_{\mathrm{2}}\nu }{\int }_{\mathrm{0}}^{\beta }{\int }_{R_{\mathrm{i}}/R_{\mathrm{o}}}^{\mathrm{1}}\overline{P_{\dot{Z}}}\overline{R}\mathrm{d}\overline{R}\mathrm{d}\theta（26）

2 数值计算求解方法

止推箔片气体动压轴承的静动态特性数值求解流程如图5所示. 首先基于变密度变黏度的湍流雷诺方程求解气膜压力分布. 其次基于Link-Spring模型求解轴承刚度，然后基于刚度和气膜压力可求出箔片变形量，并与初始设置的气膜厚度进行运算，可求解轴承的气膜厚度变化.随后根据热平衡条件计算气膜平均温升，将上述3个模型和公式在MATLAB中耦合. 每次迭代计算时，都需要基于变密度变黏度的湍流雷诺方程求解气膜压力分布和基于Link-Spring模型^［12］求解箔片变形，基于温度计算方法求解气膜平均温升，基于BP神经网络算法训练的物性模型计算出每一个节点新的密度和黏度，用于下一步的迭代计算中. 满足收敛条件1后，获得止推箔片气体动压轴承静态特性参数，此时达到了动态分析中的静平衡位置. 耦合静平衡位置雷诺方程、位移扰动方程和速度扰动方程，直至满足收敛判断条件2，计算止推箔片气体动压轴承的动态轴向刚度系数和阻尼系数. 其中，收敛判断条件1和收敛判断条件2分别为：

\frac{P-P_{\mathrm{o}\mathrm{l}\mathrm{d}}}{P_{\mathrm{o}\mathrm{l}\mathrm{d}}}\le {\mathrm{10}}^{-\mathrm{3}}（27）

\frac{P_Z-P_{Z\mathrm{o}\mathrm{l}\mathrm{d}}}{P_{Z\mathrm{o}\mathrm{l}\mathrm{d}}}\le {\mathrm{10}}^{-\mathrm{3}},\frac{P_{\dot{Z}}-P_{\dot{Z}\mathrm{o}\mathrm{l}\mathrm{d}}}{P_{\dot{Z}\mathrm{o}\mathrm{l}\mathrm{d}}}\le {\mathrm{10}}^{-\mathrm{3}}（28）

式中：P_{\mathrm{o}\mathrm{l}\mathrm{d}}为前一次循环过程P的计算结果；P_{Z\mathrm{o}\mathrm{l}\mathrm{d}}为前一次循环过程P_Z 的计算结果；P_{\dot{Z}\mathrm{o}\mathrm{l}\mathrm{d}}为前一次循环过程P_{\dot{Z}}的计算结果.

本文对止推箔片气体动压轴承进行静动态分析的初始参数如表1所示.

3 计算结果与分析

3.1 物性模型训练结果与对比分析

本文基于BP神经网络建立了S-CO₂的物性模型，该模型可以准确地描绘温度、压力与密度、黏度、比热容之间的映射关系，训练迭代次数达970次，均方差达0.050 7，如图6所示，可以看出误差均方差随着训练次数增加逐渐减小，直至所规定的期望误差停止训练. 图7为基于BP神经网络的物性模型训练结果，其中圆圈表示数据点，直线表示实际黏度值、密度值与预测密度值、黏度值之间的最佳拟合曲线. 其训练过程、验证过程和测试过程的回归系数均高达99.997%，表明该模型具有很高的可靠性和可信度. 物性模型与输出结果对比如图8所示，由图8可知，物性模型实际输出的结果与NIST数据库中得到的数据具有很好的吻合性，说明基于BP神经网络训练的物性模型能很好地反映出压力、温度与密度、黏度等物性之间的映射关系.

3.2 静态特性计算结果及验证

为证明本文所提出的模型具有可行性与有效性，将本文的计算结果与文献［19］中数据进行验证分析，取半径中间处不同周向角度下的无量纲压力值进行对比验证，模型验证结果如表2所示.由表2可知，本文所提出的流体润滑模型的计算结果与文献［19］中的数据具有较高的吻合度，因此本文所提出的理论模型具有较高的可行性与有效性.

图9（a）和图9（b）分别为2种介质下无量纲气膜压力分布情况，图9（c）为轴承沿径向展开的无量纲气膜压力分布对比. 从图9（c）中可以看出，在相同轴承参数下，2种介质的止推箔片气体动压轴承的无量纲气膜压力的变化趋势大致相同，在高压条件下， S-CO₂介质的止推箔片气体动压轴承的无量纲气膜压力相对较小，这是因为在高压条件下S-CO₂的密度比常温常压下的空气介质要大得多，在雷诺方程中，密度作为分数的分子，其数值的增大会导致雷诺方程中无量纲压力减小.

图9（d）是在50×10³ r/min的转速和相同轴承参数情况下，止推箔片气体动压轴承在S-CO₂介质和常温常压空气介质下的无量纲气膜厚度变化量分布情况，其中\mathrm{\Delta }{\overline{H}}_{\mathrm{空气}}为常温常压空气介质下的无量纲气膜厚度变化量，\mathrm{\Delta }{\overline{H}}_{\mathrm{S}-\mathrm{C}{\mathrm{O}}_{\mathrm{2}}}为S-CO₂介质下的无量纲气膜厚度变化量.从图9（d）中可以看出，在相同轴承参数下，二者的无量纲气膜厚度变化量分布呈现相同的变化趋势，在相同转速下，S-CO₂介质的气膜厚度变化量比常温常压下空气介质的变化量要大，这是因为前者的气膜不需要达到相对较高的刚度就足以支承高转速下转子的轴向载荷，还可能因为S-CO₂介质的黏度比空气介质的黏度大，在激发动压效应时，会有更多的S-CO₂气体被带入轴承间隙中，进而增大轴承气膜厚度.

图10是基于S-CO₂介质的止推箔片气体动压轴承和基于常温常压下空气介质的止推箔片气体动压轴承在不同转速下的承载能力对比，其中W_{\mathrm{空气}}为常温常压空气介质的承载力，W_{\mathrm{S}-\mathrm{C}{\mathrm{O}}_{\mathrm{2}}}为S-CO₂介质下的承载力. 从图10中可以看出，随着转速的上升，2种介质的止推箔片气体动压轴承的承载力均显著提高. 相比于常温常压下空气介质的轴承，以S-CO₂为介质的轴承承载能力有较大提升，并且在转速升高时承载能力提升得更为迅速. 这是由于S-CO₂的密度远大于空气介质，并且环境压力也远大于空气介质，根据式（10）可得S-CO₂介质下的承载力远大于常温常压空气介质下的承载力.此外，润滑模型中湍流系数大于12，根据流体润滑模型可知，湍流系数增大类似于增大动态黏度或轴承数，进而提升轴承的承载能力.

图11是在不同转速下，止推轴承正常运转时气膜的平均温升变化. 从图11中可以看出，随着转速上升，气膜平均温升呈递增趋势.根据式（15），随着转速上升，气膜的平均摩擦功率相应提高. 根据 式（19），当摩擦功率增大时，气膜平均温升随之增大；随着转速上升，止推轴承的气膜压力增大，进一步增大了压力梯度，同时，S-CO₂的黏度随着压力的增加而增加，根据式（13）可得，随着压力的增大，气膜的黏性摩擦力随之增大，导致气膜的摩擦功率增大，从而使气膜温升加剧.

3.3 气膜厚度对轴承静态特性的影响

图12是基于S-CO₂介质的止推箔片气体动压轴承在不同转速、不同最小初始气膜厚度及相同膜厚比下的静态承载力变化规律. 从图12中可以看出，在增大初始气膜厚度的情况下，止推轴承的承载力显著下降，并且在高转速下下降更为明显，最小气膜间隙从0.04 mm增大至0.06 mm，轴承承载能力下降为原来的40%. 这是因为当最小初始气膜厚度增大时，气膜厚度H作为雷诺方程中的分子增大，使得同为分子的P减小，动压效应变弱，较少的S-CO₂气体被带入轴承中，导致气膜的刚度会有所下降，从而影响止推轴承的静态承载力.

图13是基于S-CO₂介质的止推箔片气体动压轴承在不同转速、不同膜厚比下的静态承载力变化规律. 从图13中可以看出，随着膜厚比的增加，轴承的静态承载力也随之增加. 这是因为止推箔片气体动压轴承的顶箔是由一个倾斜平面和一个水平平面构成的，而膜厚比越大，其倾斜平面的倾斜角度越大，从而使得动压效应越强，越多的S-CO₂气体被拽入至轴承中，从而形成更大的气膜压力，进而导致止推轴承的承载力上升，这与文献［26］的分析结果基本一致.

3.4 动态特性计算结果及对比分析

为验证本文所提出的动态分析模型的准确性，采用与文献［9］的相同形状尺寸与环境参数的模型进行动态分析，并将本文的计算结果与文献［9］中结果进行对比. 对比结果如图14所示，由图14可知，本文所提出的动态分析模型的计算结果与文献［9］中结果具有较高的吻合度.

图15是在不同转速情况下，S-CO₂介质和常温常压空气介质下的止推箔片气体动压轴承的动态刚度系数和动态阻尼系数对比结果，其中K_{Z\mathrm{空气}}和K_{Z\mathrm{S}-\mathrm{C}{\mathrm{O}}_{\mathrm{2}}}分别为常温常压空气介质下和S-CO₂介质下的动态刚度系数，D_{Z\mathrm{空气}}和D_{Z\mathrm{S}-\mathrm{C}{\mathrm{O}}_{\mathrm{2}}}分别为常温常压空气介质下和S-CO₂介质下的动态阻尼系数. 从15图中可以看出，随着转速上升，两种介质下的动态刚度系数均呈现增大的趋势. 这是因为止推箔片气体动压轴承的动态特性主要受润滑气膜刚度和箔片结构的影响，随着转速升高，气膜刚度显著增大，导致轴承动态刚度系数增大，同时，在转子工作过程中，润滑气膜减缓了气膜能量的耗散，从而降低了动态阻尼系数. 相较于常温常压空气介质，以S-CO₂为介质的止推箔片气体动压轴承具有更大的动态刚度系数和动态阻尼系数，这是因为S-CO₂介质下的轴承工作时产生湍流效应，其湍流系数大于12，进而使得扰动方程中的\overline{P_Z}和\overline{P_{\dot{Z}}}进一步增大，导致S-CO₂介质下的动态刚度系数和动态阻尼系数更大，这意味着S-CO₂介质的止推箔片气体动压轴承具有更好的抵抗变形能力和良好的抗振性.

图16是基于S-CO₂介质的止推箔片气体动压轴承在不同转速、不同最小初始气膜厚度下，相同膜厚比的动态特性变化规律. 从图16中可以看出，最小初始气膜厚度对轴承的动态刚度系数和动态阻尼系数均有显著影响，随着最小初始气膜厚度的增大，轴承的动态刚度系数会减小，动态阻尼系数也呈现相同的变化趋势，最小初始气膜厚度从0.04 mm增至0.06 mm，轴承的动态刚度系数减小为原来的60%，动态阻尼系数减小为原来的70%. 这是因为最小初始气膜厚度减小，意味着轴承的动压效应增大，更多的S-CO₂气体被带入轴承间隙中，进而导致轴承受到更大的载荷，轴承在受到更大的载荷时，箔片结构发生了更大的变形、气膜被进一步压缩导致气膜厚度减小，气膜刚度增加，进而增大了轴承的动态刚度系数和动态阻尼系数.

4 结 论

相较于传统蒸汽动力循环系统，S-CO₂布雷顿循环系统具有更高的效率、更紧凑的结构以及良好的工程应用前景，而具有无油、高转速、耐高温高压、系统结构简单紧凑等特点的箔片气体轴承成为S-CO₂旋转机械转子系统的理想支承部件.

本文基于BP神经网络算法建立了S-CO₂的物性模型，提出了考虑湍流效应的箔片气体动压轴承综合力学模型，建立了轴承静动态弹流耦合数值计算方法，分析了以S-CO₂为介质的轴承特性以及不同轴承结构参数对轴承特性的影响规律.本文所得到的结论如下：

1）本文基于BP神经网络建立的物性模型可靠性高.通过比较基于BP神经网络训练的S-CO₂模型和NIST数据库中在不同温度、压力条件下输出的黏度、密度等数据，结果表明，训练过程相关系数高达99.997%，基于BP神经网络训练的S-CO₂物性模型的输出数据与NIST数据库的数据具有较好的吻合性.

2）增大膜厚比或减小最小初始气膜厚度可以提升其静态性能.相对于常温常压下空气介质，S-CO₂介质箔片气体动压轴承不仅具有更高的承载力，而且随着转速的升高，其提升程度也更为显著.此外，减小最小初始气膜厚度或增加膜厚比均能显著提高S-CO₂介质下的承载力，当转速不变时，将最小初始气膜厚度从0.04 mm增至0.06 mm，轴承承载能力下降为原来的40%.

3）S-CO₂介质箔片气体轴承与常温常压下空气介质相比具有更好的动态特性， 随着最小初始气膜厚度的减小，动态刚度系数和动态阻尼系数也随之增加；当转速不变，最小初始气膜厚度从0.04 mm增至0.06 mm时，轴承动态刚度系数减小为原来的60%，动态阻尼系数减小为原来的70%.

参考文献

[1] 浮翔. 新型跨临界CO₂冷热电联合循环性能研究与分析［D］. 重庆：重庆大学， 2014.

[2] FU X. Performance and exergy analyses of novel combined cooling， heating and power with transcritical CO₂ cycles［D］. Chongqing：Chongqing University，2014. （in Chinese）

[3] ZHANG X R， YAMAGUCHI H， FUJIMA K，et al ．Theoretical analysis of a thermodynamic cycle for power and heat production using supercritical carbon dioxide［J］．Energy，2007，32（4）：591-599．

[4] DANG C B， IINO K， FUKUOKA K，et al ．Effect of lubricating oil on cooling heat transfer of supercritical carbon dioxide［J］．International Journal of Refrigeration，2007，30（4）：724-731．

[5] BROUGH D， JOUHARA H. The aluminium industry： a review on state-of-the-art technologies，environmental impacts and possibilities for waste heat recovery［J］. International Journal of Thermofluids， 2020（1）： 100007．

[6] 温建全 ．超临界二氧化碳介质箔片轴承弹流耦合研究［D］．哈尔滨： 哈尔滨工业大学， 2017．

[7] WEN J Q ．Theoretical study on characteristics of compliant foil bearings lubricated with supercritical carbon dioxide［D］．Harbin：Harbin Institute of Technology，2017．（in Chinese）

[8] DOSTAL V， HEJZLAR P， DRISCOLL M J. High-performance supercritical carbon dioxide cycle for next-generation nuclear reactors［J］. Nuclear Technology， 2006， 154（3）： 265-282．

[9] LI M J， ZHU H H， GUO J Q，et al ．The development technology and applications of supercritical CO₂ power cycle in nuclear energy，solar energy and other energy industries［J］. Applied Thermal Engineering，2017，126：255-275．

[10] 谢永慧， 王雨琦， 张荻， 等. 超临界二氧化碳布雷顿循环系统及透平机械研究进展［J］. 中国电机工程学报， 2018， 38（24）：7276-7286．

[11] XIE Y H， WANG Y Q， ZHANG D， et al. Review on research of supercritical carbon dioxide brayton cycle and turbomachinery［J］．Proceedings of the CSEE，2018，38（24）：7276-7286．（in Chinese）

[12] 李卓聪，李志刚，李军 ．超临界二氧化碳（S-CO₂）气体箔片止推轴承动力特性研究［J］. 工程热物理学报，2020，41（12）：2964-2969．

[13] LI Z C， LI Z G， LI J ．Investigation on the dynamic characteristics of gas foil thrust bearings for supercritical carbon dioxide（S-CO₂） cycles［J］．Journal of Engineering Thermophysics，2020，41（12）：2964-2969．（in Chinese）

[14] DELLACORTE C， VALCO M J ．Load capacity estimation of foil air journal bearings for oil-free turbomachinery applications［J］．Tribology Transactions，2000，43（4）：795-801．

[15] 张凯， 余可， 郭志阳， 等. 嵌入弹簧式气体箔片轴承动态特性实验研究［J］. 湖南大学学报（自然科学版）， 2020， 47（12）：38-45．

[16] ZHANG K， YU K， GUO Z Y，et al ．Experimental investigation on dynamic characteristics of nested compression spring gas foil bearing［J］. Journal of Hunan University （Natural Sciences），2020， 47（12）： 38-45．（in Chinese）

[17] FENG K， KANEKO S ．Analytical model of bump-type foil bearings using a link-spring structure and a finite-element shell model［J］．Journal of Tribology， 2010， 132（2）： 021706．

[18] 邱烁现， 徐方程. 平箔片凹陷对止推箔片轴承性能影响试验研究［J］. 机械工程学报， 2022， 58（15）： 226-232．

[19] QIU S X， XU F C ．Experimental study on the effect of top foil sagging on the performance of thrust foil bearing［J］．Journal of Mechanical Engineering，2022，58（15）：226-232．（in Chinese）

[20] 朱鹏程， 管玉坤， 门日秀， 等. 超临界二氧化碳系统用箔片气体动压轴承静特性研究［J］. 振动工程学报， 2024， 37（5）：875-884．

[21] ZHU P C， GUAN Y K， MEN R X， et al. Static characteristics of gas foil bearing for supercritical carbon dioxide system［J］．Journal of Vibration Engineering，2024，37（5）：875-884．

[22] DELLACORTE C， BRUCKNER R J ．Remaining technical challenges and future plans for oil-free turbomachinery［J］．Journal of Engineering for Gas Turbines and Power，2011，133（4）：042502．

[23] 冯凯，黄明，李成勤，等.径向气体箔片轴承高速重载测试及实验台搭建［J］.湖南大学学报（自然科学版），2020，47（2）：1-5.

[24] FENG K， HUANG M， LI C Q， et al. Test of high speed and heavy load gas foil bearing and set up of test rig［J］.Journal of Hunan University（Natural Sciences），2020，47（2）：1-5. （in Chinese）

[25] 刘良军. 弹性箔片气体动压推力轴承的静动态特性和热特性研究［D］.长沙： 湖南大学， 2015.

[26] LIU L J. Theoretical considerations of hydrodynamic and thermal characteristics of gas foil thrust bearings［D］. Changsha：Hunan Universuty，2015. （in Chinese）

[27] XU F C， KIM D. Three-dimensional turbulent thermo-elastohydrodynamic analyses of hybrid thrust foil bearings using real gas model［C］//Proceedings of ASME Turbo Expo 2016：Turbomachinery Technical Conference and Exposition，June 13-17， 2016. Seoul： American Society of Mechanical Engineers， 2016.

[28] QIN K， JAHN I， GOLLAN R，et al ．Development of a computational tool to simulate foil bearings for supercritical CO₂ cycles［J］．Journal of Engineering for Gas Turbines and Power，2016，138（9）：092503．

[29] KIM D ．Design space of foil bearings for closed-loop supercritical CO₂ power cycles based on three-dimensional thermohydrodynamic analyses［J］．Journal of Engineering for Gas Turbines and Power， 2016， 138（3）： 032504．

[30] QIN K， JAHN I H， JACOBS P A. Effect of operating conditions on the elastohydrodynamic performance of foil thrust bearings for supercritical CO₂ cycles［J］. Journal of Engineering for Gas Turbines and Power，2017，139（4）：042505．

[31] NG C W， PAN C H T. A linearized turbulent lubrication theory［J］．Journal of Basic Engineering，1965，87（3）：675-682．

[32] 温诗铸， 黄平， 田煜， 等. 摩擦学原理［M］. 5版. 北京： 清华大学出版社， 2018： 19-108．

[33] WEN S Z， HUANG P， TIAN Y， et al. Principles of tribology［M］．5th ed．Beijing：Tsinghua University Press，2018：19-108．（in Chinese）

[34] 赵雪源 ．新型高阻尼气体箔片轴承及其转子系统的理论和实验研究［D］．长沙： 湖南大学，2017．

[35] ZHAO X Y ．Theoretical analysis and experimental investigation on the performance of a novel gas foil bearing with high structural damping and its rotordynamic response［D］. Changsha：Hunan University，2017．（in Chinese）

[36] ZHANG K， ZHAO X Y， FENG K，et al ．Thermohydrodynamic analysis and thermal management of hybrid bump-metal mesh foil bearings：experimental tests and theoretical predictions［J］．International Journal of Thermal Sciences，2018，127：91-104．

[37] HESHMAT C A， XU D S， HESHMAT H ．Analysis of gas lubricated foil thrust bearings using coupled finite element and finite difference methods［J］．Journal of Tribology，2000，122（1）：199-204．

基金资助

国家重点研发计划资助项目（2021YFF0603000）