图1示出了波箔型气体箔片推力轴承结构示意图，轴承主要由轴承座、波箔、顶箔及推力盘组成，顶箔位于波箔上方，波箔为顶箔提供弹性支承，两箔片的一端共同固定于轴承座上，由此构成若干扇形瓦块。扇形瓦块结构参数如图2(a)所示，其中${R_1}$为轴承内径，${R_2}$为轴承外径，$\beta $为瓦块张角，b为节距比（$b = \gamma /\beta $）。图2(b)为扇形瓦块剖面图，推力盘转速较高时，推力盘将气体由进气口带入扇形瓦块区域，然后通过出气口排出，初始气膜间隙由${h_1}$减小为${h_2}$，顶箔和推力盘之间形成楔形润滑气膜，为推力盘提供支承力，实现气体动压润滑。

图  1  波箔型气体箔片推力轴承结构

Figure 1.  Structure of bump-type gas foil thrust bearings

图  2  扇形瓦块结构示意图

Figure 2.  Diagram of sector pad structure

波箔在工作中的变形如图3所示，气膜压力通过顶箔的传递，间接作用于波箔，使波箔发生弹性变形。由于波箔结构在轴承圆周方向上均匀分布，因此对单波进行的受力变形分析适用于整块波箔。

为便于计算分析，作如下假设：(1) 波箔的顶部和顶箔之间始终保持接触，单波所受载荷为作用于波峰的集中力；(2) 由于波箔和顶箔、波箔和轴承座之间的摩擦很小，故可忽略；(3) 波箔受力后的弹性变形量远小于原始尺寸，属于小变形；(4) 如图3(b)所示，由于单波周向变形量很小，故可认为波箔变形后仍保持为一个扇形，且内外径不变。

波箔的单波在工作过程中的变形如图4所示，单波横截面尺寸为${t_{\rm{B}}} \times L$（L为波箔内外半径之差），波拱半角为$\varphi $，波拱半长为l，其顶部受到竖直向下的载荷F，由于其横截面高度${t_{\rm{B}}}$远小于其轴线半径${R_{\rm{0}}}$，因此可利用卡氏定理计算其变形。相较于单波顶点的竖直位移$\Delta y$，其水平位移$\Delta x$对气膜厚度的影响非常小，故本文不予考虑。下面计算$\Delta y$，若规定使单波曲率增大的弯矩为正，则单波任意截面m-m上的弯矩M为：

M关于F的偏导数为：

对单波左半段进行积分，并将式(1)、(2)代入可得单波顶点的竖直位移为：

式中：E为箔片弹性模量。

单波横截面惯性矩可表示为${I_{\rm{B}}} = \dfrac{{L{t_B}^3}}{{12}}$，代入式(3)可得：

图  3  波箔变形示意图

Figure 3.  Diagram of bump foil deformation

图  4  单个波箔的变形示意图

Figure 4.  Diagram of single bump foil deformation

轴承在工作时，顶箔在气膜压力和波箔支承力的作用下将发生一定的弹性变形。为计算该弹性变形量，将顶箔视为二维梁进行分析。如图5所示，选取位于两单波顶点之间的某一段顶箔作为分析对象，由于其长度s（即波距）较小，故可将该区间内的气膜压力视为均布载荷q。

图  5  顶箔受力变形示意图

Figure 5.  Diagram of top foil deformation under force

变形方程可表示为：

式中：${w_1}$为顶箔竖直变形量，x为距左顶点的距离，${I_{\rm{T}}}$为顶箔横截面惯性矩。

对式(5)积分可得：

暂时仅考虑顶箔变形，假设波箔不变形，可得以下边界条件：

将式(7)代入式(6)，可求得常数项，得到：

箔片结构总体变形由顶箔和波箔共同决定，其值可表示为：

式中：${w_1}$为顶箔自身产生的挠度，可由式(8)求得，${w_2}$为顶箔在波箔变形影响下所产生的竖直位移。如图6所示，波箔两个顶点的竖直位移分别为${w_{{\rm{B1}}}}$和${w_{{\rm{B2}}}}$，水平位移分别为$\Delta {x_1}$和$\Delta {x_2}$，则由几何关系可得，距左顶点x处的位移${w_2}$为：

图  6  顶箔在波箔变形影响下所产生位移示意图

Figure 6.  Diagram of the displacement of top foil with the deformation of bump foil

联立可压缩气体的运动方程、连续性方程和能量守恒方程，并结合理想气体状态方程，可得用来描述顶箔与推力盘之间气体流动规律的等温定常Reynolds方程：

式中：r为径向坐标，$\theta $为周向坐标，p为气膜压力，h为气膜厚度，$\mu $为气体的动力黏度，$\omega $为转轴角速度。

将式(11)量纲为-化得：

式中：$\bar r = \dfrac{r}{{{R_2}}}$，$P = \dfrac{p}{{{p_0}}}$，$H = \dfrac{h}{{{\delta _h}}}$，$\Lambda = \dfrac{{6\mu \omega {R_2}^2}}{{{p_0}{\delta _h}^2}}$，其中${p_0}$为环境压力，${\delta _{\rm{h}}}$为进出口气膜厚度之差，$\Lambda $为轴承系数。

根据图2的几何关系，可得推力盘无倾斜时的气膜厚度表达式：

式中：${w_t}(r,\theta )$为箔片结构的总体变形量，$g(\theta )$可表示为：

故量纲为-气膜厚度表达式为：

式(12)所示的Reynolds方程为非线性偏微分方程，本文采用有限差分法及Newton-Raphson迭代法进行求解^[13]。当轴承稳定工作时，其每个扇形瓦块内的气膜压力分布及气膜厚度分布相同，故取任一个瓦块作为计算域即可。如图7所示，轴承气膜承载区以波箔固定端和自由端为周向边界，以轴承内外径处为径向边界。将其离散为$m \times n$的网格，压力值${P_{i,j}}$表示各节点上气膜压力值，每格步长为：$\Delta \theta = \dfrac{\beta }{m}$，$\Delta \bar r = \dfrac{{{R_2} - {R_1}}}{{n{R_2}}}$。

为提高计算精度，采用中差商表示偏微分项：

将式(16)~(18)代入式(12)，化简得：

式中：${A_{i,j}} = \dfrac{{{{\left( {P{H^3}\bar r} \right)}_{i,j + 1/2}}}}{{{{\bar r}_{i,j}}\Delta {{\bar r}^2}}}$，${B_{i,j}} = \dfrac{{{{\left( {P{H^3}\bar r} \right)}_{i,j - 1/2}}}}{{{{\bar r}_{i,j}}\Delta {{\bar r}^2}}}$，${C_{i,j}} = \dfrac{{{{\left( {P{H^3}} \right)}_{i + 1/2,j}}}}{{{{\bar r}_{i,j}}^2\Delta {\theta ^2}}}$，${D_{i,j}} = \dfrac{{{{\left( {P{H^3}} \right)}_{i - 1/2,j}}}}{{{{\bar r}_{i,j}}^2\Delta {\theta ^2}}}$，${E_{i,j}} = {A_{i,j}} + {B_{i,j}} + {C_{i,j}} + {D_{i,j}}$，${F_{i,j}} = \Lambda \dfrac{{{{\left( {PH} \right)}_{i + 1/2,j}} - {{\left( {PH} \right)}_{i - 1/2,j}}}}{{\Delta \theta }}$

求解过程中引入边界条件，可压缩气体润滑中不存在气膜破裂现象，轴承承载区域各边界均与大气环境连通，故其压力与环境压力相同，无量纲形式的边界条件为：$P\left( {0,\bar r} \right) = 1{\rm{ }},P\left( {\beta ,\bar r} \right) = 1{\rm{ }},P\left( {\theta ,\dfrac{{{R_1}}}{{{R_2}}}} \right) = 1{\rm{ }}, P\left( {\theta ,1} \right) = 1$。

引入下列相对收敛准则：

图8示出了本文流固耦合的数值计算过程。

图  7  网格划分及差商示意图

Figure 7.  Mesh generation and difference quotient

图  8  计算流程图

Figure 8.  Calculation flow chart

为验证上述数值解法的正确性，选用与文献[4]相同的刚性气体推力轴承参数，按上述算法借助MATLAB编程计算其单个瓦块的承载力，计算结果对比如表1所示。由对比结果可知，对于不同气膜间隙比和轴承数下的单个瓦块承载力，本文计算所得结果与参考文献结果的误差均小于3%，具有较高的吻合度。

$\frac{{{h_1}}}{{{h_2}}}$	$\Lambda $	Dimensionless bearing capacity $\bar W \times {10^2}$		Deviation/%
		Results in this paper	Results from literature
2	0.1	0.0137	0.014	−2.14
	1	0.1402	0.144	−2.64
	10	1.6000	1.612	−0.74
	20	3.3276	3.33	−0.07
	40	6.2657	6.246	0.32
5	0.1	0.1866	0.192	−0.03
	1	2.3580	2.343	0.64
	10	22.7503	22.2	2.48
	20	35.5801	34.66	2.65
	40	47.7257	46.4	2.86

表 1  计算结果对比

Table 1.  Comparison of calculation results

在刚性气体推力轴承的基础上，以某一个具体波箔型气体箔片推力轴承为算例进行研究,计算时取$m = 60$，$n = 50$。轴承瓦块数目为6，各瓦块波纹数目均为12，其他轴承结构及运行参数如表2所示^[14]。图9至图11分别示出了该波箔型气体箔片推力轴承的量纲为-气膜压力分布图、量纲为-气膜厚度分布图及箔片结构变形量分布图。

图12示出了箔片推力轴承和刚性推力轴承在扇形瓦块径向中点处（即$\overline r = 0.5$）的量纲为-气膜厚度对比图及其局部放大图。由图12(a)可知，两者气膜厚度的大小相近，且沿周向的变化趋势相似；由图12(b)可知，箔片推力轴承的气膜厚度略大于刚性推力轴承，且有小幅波动。图13示出了不同节距比下箔片推力轴承和刚性推力轴承的承载力及摩擦力矩对比图。由图13(a)可知，箔片推力轴承的承载力略小于刚性推力轴承，且当节距比b在0.6附近时，两者承载力均达到最大值。图13(b)可知，箔片推力轴承的摩擦力矩略小于刚性推力轴承，两者均随节距比b的增大而减小，且近似于线性变化。

${R_1}$/m	${R_2}$/m	$\beta $/(°)	b	${t_{\rm{B}}}$/m	${t_{\rm{T} }}$/m	${R_0}$/m	$\varphi $/rad
$2.54 \times {10^{ - 2}}$	$5.08 \times {10^{ - 2}}$	60	0.5	$1.016 \times {10^{ - 4}}$	$1.524 \times {10^{ - 4}}$	${\rm{1}}{\rm{.2125}} \times {10^{ - 3}}$	0.8364
l/m	E/Pa	μ/(Pa·s)	p0/Pa	${h_1}$/m	${h_2}$/m	n/$(r \cdot {\min ^{ - 1}})$
${\rm{9}} \times {10^{ - {\rm{4}}}}$	$2.14 \times {10^{11}}$	$1.932 \times {10^{ - 5}}$	$1.01325 \times {10^5}$	$5.2 \times {10^{ - 5}}$	$2.6 \times {10^{ - 5}}$	30000

表 2  轴承结构及运行参数

Table 2.  Parameters of bearing structure and operation

图  9  量纲为-气膜压力分布

Figure 9.  Distribution of dimensionless gas film pressure

图  10  量纲为-气膜厚度分布

Figure 10.  Distribution of dimensionless gas film thickness

图  11  箔片结构变形量分布

Figure 11.  Distribution of foil structure deformation

图  12  量纲为-气膜厚度对比及其局部放大图

Figure 12.  Comparison of dimensionless gas film thickness and its partial enlargement drawing

图  13  箔片轴承与刚性轴承润滑性能对比

Figure 13.  Comparison of lubrication performance of foil bearing and that of rigid bearing

图14至图17为其他参数均不变时，各节距比b下波箔厚度${t_{\rm{B}} }$、波拱半长l、波拱半径${R_0}$和顶箔厚度${t_{\rm{T}}}$对轴承润滑性能的影响关系图。由图可知：轴承的承载力和摩擦力矩先随波箔厚度增大而快速增大，当波箔厚度达到一定值后，两者增大速度均趋于平缓；承载力和摩擦力矩均随波拱半长增大而缓慢减小；波拱半径较小时，承载力和摩擦力矩受其变化影响较小，波拱半径较大时，两者均随其增大而减小；随着顶箔厚度的增大，承载力和摩擦力矩先快速增大，之后几乎不再变化。因此，可通过增大波箔厚度${t_{\rm{B}}}$、顶箔厚度${t_{\rm{T}}}$或减小波拱半长l、波拱半径${R_0}$，来提高波箔型气体箔片推力轴承的承载力，反之可降低轴承摩擦力矩。

图  14  不同节距比下波箔厚度对润滑性能的影响

Figure 14.  Influences of bump foil thickness on lubrication performance with different pitch ratios

图  15  不同节距比下波拱半长对润滑性能的影响

Figure 15.  Influences of half bump length on lubrication performance with different pitch ratios

图  16  不同节距比下波拱半径对润滑性能的影响

Figure 16.  Influences of bump radius on lubrication performance with different pitch ratios

图  17  不同节距比下顶箔厚度对润滑性能的影响

Figure 17.  Influences of top foil thickness on lubrication performance with different pitch ratios

(1) 以波箔型气体箔片推力轴承为研究对象，运用弹性力学理论构建了波箔和顶箔的受力变形模型，得到了考虑箔片弹性变形的气膜厚度方程，结合可压缩气体Reynold方程，构建了一种流固耦合算法，并利用有限差分法和Newton-Raphson迭代法进行求解，借助MATLAB编程实现了对轴承润滑性能的精确数值计算。结合具体文献对所建立的计算模型进行了对比验证，所得结果与文献结果有良好的吻合度，证明了本文算法的可靠性。

(2) 通过对具体算例的研究，得到了轴承的量纲为-气膜压力分布、量纲为-气膜厚度分布与箔片结构变形量分布。在此基础上对比研究了箔片推力轴承和刚性推力轴承的润滑性能，结果表明，箔片轴承的气膜厚度较大，承载力和摩擦力矩较小。结合算例进一步研究了波箔厚度${t_{\rm{B}}}$、波拱半长l、波拱半径${R_0}$和顶箔厚度${t_{\rm{T}}}$对轴承润滑性能的影响规律，得出具体的影响规律曲线图。结果表明：轴承的承载力和摩擦力矩随${t_{\rm{B}}}$或${t_{\rm{T}}}$的增大而增大，随l或${R_0}$的增大而减小。故较大的${t_{\rm{B}}}$、${t_{\rm{T}}}$或较小的l、${R_0}$可提高轴承的承载力，反之可减小轴承的摩擦力矩，有利于降低摩擦功耗，减少发热量。