悬臂止推箔片轴承三维流-固-热耦合数值分析

赵 琪, 强铭琛, 侯 予,2, 陈双涛,2, 赖天伟,2*

(1.西安交通大学 能源与动力工程学院, 陕西 西安 710049;2.西安交通大学 深低温技术与装备教育部重点实验室, 陕西 西安 710049)

气体箔片轴承具有结构简单、成本低、无污染和高速稳定等诸多优点[1-3]，其在高速旋转机械中应用十分广泛[4-6].为了满足更高功率密度旋转机械的需求，气体箔片轴承通常需要在高转速下工作.此时，由于气体黏性带来的高速剪切耗散热量会使得轴承工作温度迅速升高，热量若不能及时散失，润滑气膜温度过高时，可能会发生轴承热失效等故障[7-8]，因此进行箔片轴承的散热分析以及承载力特性分析在轴承性能研究中非常重要[9].

箔片轴承散热和承载力特性分析是1个比较复杂的流-固-热耦合问题，为了深入了解箔片轴承的承载特性，学者们已经进行了大量的理论和试验研究工作[10-12].Andrés和Kim[13]建立了波箔径向轴承的二维流-固-热耦合分析模型，模型中考虑到了物性变化和轴承部件的热膨胀.Lehn等[14]提出了1种考虑止推盘热膨胀与止推盘背面流体冷却效应的波箔止推轴承三维流-固-热耦合模型，该模型耦合了轴承不同部件，并提出了箔片之间以及箔片与轴承座间的接触热阻解析式，结果表明高转速下止推盘的受热弯曲变形会导致轴承承载能力下降.郭雨等[15]建立了悬臂轴承二维流-固-热分析模型，并开展了冷却气对轴承散热影响的试验研究，结果表明引入冷却空气后，箔片背面温度明显下降，而且冷却效果随着供应压力的增加而提高.张凯等[16]针对1种波箔-金属网混合型径向箔片轴承建立了流-固-热耦合模型，通过导热模型计算了气膜通过转子和壳体向外界的散热量.基于数学模型对不同载荷与转速下轴承温度分布进行了计算.Radil等[17]比较了3种箔片轴承热管理控制方法，分别通过冷却气直接垂直吹击转子内表面、冷却气沿轴向吹入空心转子产生对流以及冷却气吹入轴承支承部件的方式对轴承转子系统进行降温，试验结果表明在轴承支承结构引入冷却气流冷却效果最好.

箔片气体轴承中箔片变形会影响气膜中气体的流动，而气膜压力反过来又会作用于箔片，里面包含着复杂的流体和固体耦合作用关系，此外气膜中气体黏性剪切产生的热量又会影响气体和箔片的物理特性，因此建立三维的箔片轴承流-固-热耦合模型对轴承性能分析是非常重要的.本文中通过有限元商业软件COMSOL对带有冷却流道的悬臂型箔片轴承进行了三维流-固-热耦合数值分析，考虑了箔片周向和径向传热以及冷却气实际流动特性的影响.通过和试验结果对比验证了模型的准确性，分析了变轴承间隙和变转速条件下箔片轴承的承载力特性，得到了气膜压力场、温度场分布以及承载力等变化规律，研究了箔片的接触变形问题，此外分析了冷却气流速对轴承散热和承载力特性的影响.

1 物理与数值模型

1.1 物理模型

悬臂止推箔片轴承实物图和结构示意图分别如图1和图2所示.该轴承由1个焊接板与6片周向均匀分布的悬臂止推箔片构成.箔片外缘安装区域通过激光焊接在轴承座上，后缘自由搭接在相邻箔片上.每片箔片通过特殊的搭建的方式可以与止推盘之间形成楔形间隙，有利于形成动压润滑气膜.气膜力作用下，悬臂箔片搭接结构可以发生相对滑移，从而产生库伦阻尼耗散作用，可以有效的提高轴承转子系统的稳定性.转子转动方向与箔片延展方向相同.轴承座沿周向与径向开设有多条冷却气流槽道，可以加快箔片侧散热，从而降低箔片表面温度.悬臂止推箔片轴承的结构参数列于表1中.

表1 悬臂止推箔片轴承结构参数Table 1 Structure parameters of multi of thrust foil bearing

Fig.1 Physical diagram of multi-leaf thrust foil bearing图1 悬臂型止推箔片轴承实物图

Fig.2 Schematic diagram of multi-leaf thrust foil bearing图2 悬臂型止推箔片轴承结构示意图

在轴承运行过程中，形成的周期性楔形润滑气膜和热阻分布如图3所示，气体通过楔形气膜时，动压效应的作用可以提供轴向负载能力.在润滑气膜压力作用下，箔片之间以及箔片与焊接板之间均会发生相应的接触.为了便于建模，本文中假设前一片悬臂箔片的后缘部分与后一片箔片的前缘部分完全接触.气膜中黏性耗散热主要通过止推盘侧与箔片侧进行散热，其中，R1为止推盘表面与空气的对流换热热阻，R2为止推盘热阻，R3和R5为箔片热阻，R7为轴承座热阻，R4和R6为摩擦副间的接触热阻，R8为轴承座表面与空气的对流换热热阻，R9为箔片与轴承座间的间隙热阻.

Fig.3 Diagram of lubrication gas film and thermal resistance distribution of multi-leaf foil bearing图3 悬臂止推箔片轴承润滑气膜和热阻分布示意图

1.2 数值模型

本文中选择了悬臂轴承试验过程中的几种典型工况进行数值模拟.悬臂箔片圆周方向周期性分布，因此只对1/6轴承结构进行建模分析.悬臂箔片轴承计算模型如图4所示，其主要由几部分组成，分别为进气流道、气膜、上层箔片、底层箔片、箔片下冷却通道、轴承座中气体冷却通道以及轴承座.轴承工作过程中，润滑气膜中的黏性剪切耗散是主要产热源，热量耗散主要通过箔片表面的对流换热以及热传导来实现.数值计算中，气膜上表面设置为无滑移边界条件，换热条件为强制对流，强制对流气体平均流速为气膜中径处线速度.轴承转速ω设置为100 000 r/min，环境温度T0为20 ℃，环境压力P0为1 atm.在流-固-热耦合计算中，上层箔片与润滑气膜交界面完全耦合，考虑了双向的力传递、位移变化以及热量传递过程.同时考虑了箔片间以及箔片与轴承座接触热阻的影响.冷却气从进气流道通入，然后分别进入润滑气膜、箔片下冷却通道和轴承座冷却通道，轴承座设置为刚性域，轴承座上下表面为自然对流换热条件，左右壁面设置为周期性边界条件.

Fig.4 Calculation model diagram of multi-leaf foil bearing图4 悬臂箔片轴承计算模型示意图

数值模型中采用稳态可压缩Navier-Stokes方程对止推轴承间隙内的气膜压力场进行求解，控制方程为

式中，ρ为空气的密度(kg/m3)， ∇为微分算子，是速度矢量(m/s)，为黏性应力，p为气膜压力(Pa)，T代表转置，μ为空气黏度(Pa·s).

通过三维能量方程[18]对止推轴承气膜温度场进行求解，能量方程为

式中，k为气体导热率[W/(m·K)]，∆T为温度梯度，Cp为气体比热容[J/(kg·K)]，是黏性耗散热量(W).

根据接触热阻物理模型[19]，不同接触面间的接触热阻如式(5)所示：

式中，hcont为间隙气体的热导率[W/(m2·K)]；keq为接触面的平均导热系数[W/(m2·K)]；meq和σeq分别为接触面等效平均斜度绝对值(μm)和标准差；Pf(N/m2)和Hm(GPa)分别为接触面压力和弹性材料的显微硬度.

考虑到箔片间以及箔片与轴承座间的接触问题的复杂性，采用罚函数的方式来求解接触变形问题，本文中箔片间摩擦力以及箔片与轴承座间摩擦力均设置为0.

采用固体力学有限元求解器来求解箔片变形问题，载荷、位移与刚度矩阵的平衡方程以及整体能量方程如下所示

式中，K，u和F分别为总刚度矩阵、测量点位移矢量和测量点载荷.kfoil为箔片的热导率[W/(m·K)].

得到润滑气膜表面的压力场后，轴承载荷可以通过下面的积分方程得到

式中，N代表箔片个数，在本文中为6，θ2和θ2为前缘和尾缘角度，Ri和Ro为箔片内径和外径，p0为环境压力(Pa).

1.3 模型验证

本文中通过与文献[15]数据对比来验证悬臂轴承三维流-固-热耦合模型计算的准确性，不同轴承间隙下承载力与文献结果对比如图5所示，数值计算承载力与试验测试承载力比较接近，表明了该数学模型的有效性，轴承承载力随轴承间隙的减小而增大.在相同的转速下，保持较小的轴承间隙可以获得更大的气膜压力，从而获得更高的轴承承载力，文献中温度参考测量点的位置示意图如图6所示.悬臂型止推箔片轴承稳态运行时测量点温度的数值计算值与文献中的试验测试值随轴承承载力的变化趋势对比如图7所示，并列于表2中.结果表明，测量点温度的数值计算值与文献中试验测试值随承载力变化均呈现近似线性变化，由于接触热阻影响和箔片的热传导作用使得箔片背面内外径处温度趋于接近，数值计算值和文献中试验测试值吻合结果较好，证明了数值模型传热特性计算的准确性.本文中考虑了箔片沿径向与周向的热传导，高压区温度测量点(T1~T3)和低压区温度测量点(T4~T6)一致性较好.承载力较大时，轴承高压区测量点(T1~T3)温度数值计算值略高于文献试验值.T2测点在承载力为23.5 N时，文献试验与理论结果间温度误差较大，约为6.4 ℃，承载力为10 N时，温度误差较小，约为0.5 ℃.多种因素可能导致数值计算值与文献试验测试值的差异，如测量误差(实际测量时热电偶部分暴露在冷却气中)、热电偶安装位置误差以及数学模型假设(本文中未考虑转子热变形以及湍流等因素的影响，其可能会对气体润滑方程的求解产生影响)等.另外由于实际运行过程中箔片表面涂层会发生磨损，这会影响实际的接触热阻，进而造成理论与文献试验结果之间的误差.

表2 箔片部分测点温度计算结果与文献[15]结果对比Table 2 Comparison of foil temperatures at measurement points in the foil with literature results

Fig.5 Comparison of bearing load capacity at different bearing clearances with literature results图5 不同轴承间隙下承载力与文献[15]结果对比

Fig.6 Schematic diagram of temperature measurement points[15]图6 箔片温度测点示意图[15]

Fig.7 Comparison of foil temperature at different bearing loads with literature results图7 不同承载力下箔片温度与文献[15]结果对比

2 结果与讨论

2.1 变轴承气膜间隙

不同轴承间隙下，悬臂型止推箔片轴承的箔片变形(W)、气膜无量纲压力(P)分布和气膜温度场(T)分别如图8、图9和图10所示.由于箔片内外缘的弦长差异，箔片内外缘区域具有一定的刚度差异，导致箔片形变量沿径向持续增大.由于箔片下方间隙分布沿周向的差异，箔片形变量沿着周向亦具有一定差异，箔片形变主要集中在气膜前缘位置.由气膜前缘开始，随着气膜厚度逐渐减小，润滑介质逐渐被压缩从而使气膜压力逐渐增大.气膜峰值压力位于气膜中部位置，而气膜高压区主要位于气膜后缘位置，因为该区域具有较小的气膜厚度分布.气膜压力在度过峰值后逐渐降低，并在气膜出口处迅速降低至环境压力.轴承后缘的外侧区域压力相对较小，是因为其(因箔片弦长较大而刚度较小)具有相对较大箔片局部变形而导致较大的气膜厚度分布所致.在小间隙下，气膜高压区具有较高且更为平缓的压力分布，且轴承间隙越小，气膜高压区最大压力位置越靠近箔片后缘.在较大间隙下，动压效应减弱，高压区面积和无量纲气膜压力逐渐减小.轴承间隙从14 μm减小到8 μm时，无量纲压力最大值从1.22增加到了1.81，提升了约50%.气膜温度沿着前缘至后缘方向逐渐升高，且由于轴承外侧区域较高的线速度，使得轴承外侧区域温度较高.气膜前缘位置处由于较大的气膜厚度与较小的气膜压力梯度，气膜温度相对较小.轴承间隙从14 μm减小到8 μm时，气膜温度最大值从105增加到了206 ℃，提升了约100%.

Fig.8 Foil deformation at different bearing clearances图8 不同轴承间隙下箔片变形量

Fig.9 Dimensionless gas film pressure at different bearing clearances图9 不同轴承间隙下无量纲气膜压力分布

Fig.10 Temperature distribution of gas film under different bearing clearance图10 不同轴承间隙下气膜温度场分布

2.2 变轴承转速

不同转速下，悬臂止推箔片轴承的箔片变形量、气膜无量纲压力分布和气膜温度场分别如图11、图12和图13所示.由于箔片不同位置处的刚度差异，箔片形变主要集中在气膜前缘位置.轴承转速越大，动压效应增强，箔片整体的变形量也逐渐增大，轴承转速为140 000 r/min时，箔片最大变形量约为4.1 μm.随着转速的增加，轴承动压效应增强，承载区气膜压力显著提高.轴承转速越高，气膜高压区最大压力位置越靠近箔片后缘.轴承间隙从80 000 r/min增加到140 000 r/min时，气膜中径处无量纲压力最大值从1.21增加到了1.73，提升了约43%.轴承转速为140 000 r/min时，无量纲气膜压力最大值约为1.73.转速较大时，气体黏性剪切速度变大，黏性耗散热量增多，气膜温度逐渐升高，轴承转速从80 000 r/min增加到140 000 r/min时，气膜温度最大值从85 ℃增加到了345 ℃，提升了约305%.这说明转速升高对箔片温度场影响很大，气膜温度与转速呈现非线性增长的趋势.

Fig.11 Foil deformation at different rotation speeds图11 不同转速下箔片变形量

Fig.12 Dimensionless gas film pressure distribution at different rotation speeds图12 不同转速下无量纲气膜压力分布

Fig.13 Gas film temperature distribution at different rotation speed图13 不同转速下气膜温度场分布

2.3 变冷却气流速

本文中通过改变气体流速来研究冷却气对箔片轴承散热特性的影响.不同冷却气流速下，气膜表面沿径向方向温度如图14(a)所示，气膜温度最大值以及轴承承载力如图14(b)所示.由于靠近箔片外缘位置处线速度增大，气体黏性剪切耗散热量增多，沿着径向方向，气膜表面温度逐渐升高，近似呈现线性变化关系，由于内外径气体剪切速度存在较大差异，其温差非常大.随着冷却气流速的增大，气膜表面温度逐渐降低，流速从1 m/s增加到40 m/s时，气膜表面内缘位置处温度降低了约35 ℃左右，外缘位置降低了约21 ℃，气膜表面温度最大值从258 ℃降低至240 ℃.冷却气流速增大到一定值时，气膜表面降温效果不再发生明显变化，这说明应该控制冷却气流速，使得在保证轴承散热效果的同时仍具有一定的经济性.气膜温度降低时，气体黏性减小，随着冷却气流速增大，轴承承载力呈现下降趋势.流速从1 m/s增加至40 m/s时，轴承承载力减小了约4%.

Fig.14 (a) Gas film temperature, (b) bearing load and maximum temperature of gas film at different cooling gas flow rates(ω=140000 r/min, meq=0.4, σeq=0.4)图14 不同冷却气流速下(a)气膜温度以及(b)轴承承载力和气膜表面温度最大值(ω=140000 r/min，meq=0.4，σeq=0.4)

3 结论

为了深入分析箔片轴承的散热特性和承载特性，本文中建立了悬臂止推箔片轴承三维流-固-热耦合计算模型，研究了轴承间隙、转速以及冷却气流速变化对悬臂箔片轴承气膜温度场和承载力等的影响规律，得到的结论如下：

a.本文中建立了悬臂止推箔片轴承三维流-固-热耦合计算模型，考虑了箔片的径向和周向传热以及冷却气实际流动特性，另外，将数值计算结果与文献试验数据进行了比较，验证了模型的准确性.

b.随着轴承转速增大或者轴承间隙减小，轴承动压效应增强，气膜高压区面积和压力峰值都逐渐增大，气膜高压区主要位于气膜后缘位置.相比于轴承间隙变化，轴承转速升高对润滑气膜温度场影响更大，气膜温度与转速呈现非线性增长的趋势.转速从80 000 r/min增加至14000 r/min时，气膜温度最大值从85 ℃增大至345 ℃，提升了约305%.

c.随着冷却气流速的增大，气膜表面温度和轴承承载力均逐渐减小.冷却气流速增大到一定值时，降温效果不再发生明显变化.可以适当控制冷却气流速，使得可以在保证轴承散热效果的同时仍具有一定的经济性.