林起崟 魏正英 王宁 陈渭

(西安交通大学 机械制造系统工程国家重点实验室，陕西 西安 710049)

表面织构技术在摩擦学领域的应用越来越受到关注.研究发现应用表面织构技术可以有效降低轴承摩擦阻力，同时还可以提高其承载能力［1-2］.Etsion等［3］对表面织构对机械密封性能的影响进行了理论和实验研究.Siripuram、Fowell 和Dobrica 等［4-6］对表面织构的形状、深度等结构参数对轴承摩擦学性能的影响进行了分析，优化了相应的织构结构.

轴承表面往往含有成千上万个微织构，其摩擦学性能是由其表面所有微织构相互叠加、共同作用体现出的宏观效果.受限于计算机软硬件资源，传统的研究表面织构的方法无法建立完整的具有微织构的轴承模型，只能选择一个或多个微织构构建理论模型［2，7］，也就无法对大面积织构表面进行分析和优化，进而无法准确分析织构对轴承整体摩擦学性能的影响.研究已经证实，流固界面处的流体存在速度滑移，特别是存在表面织构时滑移现象更加显著［8-11］.织构表面界面处的速度滑移现象是表面上所有织构共同作用的综合效应.因此，理论分析时可以采用滑移边界条件代替织构表面进行研究.Salant等［12］的研究表明，表面织构类似滑移区域，滑移表面可以有效降低摩擦阻力，提高承载能力，还可以显著提高机械密封性能.目前主要采用滑移长度模型和极限剪切应力模型来表征滑移表面［12-15］，文中结合上述模型提出了一种双系数数值滑移边界，采用文献［16］的实验轴承进行分析，结果表明理论值与文献实验值吻合良好，验证了该双系数数值滑移边界的正确性.此外，文中采用一种基于压力变化的相变模型描述润滑膜空化现象，将该数值边界施加到有微织构的滑动轴承表面，研究了织构表面位置分布、面积大小等宏观参数对滑动轴承承载能力和空穴的影响规律.

1 理论方法

经典润滑理论所采用的Reynolds 方程是通过引入一系列假设性条件从完全Navier-Stokes 方程简化而来的，因此无法考虑润滑膜厚度方向的影响，还忽略了惯性力的影响，而文中轴承模型在高速运行时必须要考虑惯性力的影响，所以采用完全的Navier-Stokes 方程进行分析计算.描述连续性流体的质量守恒方程和Navier-Stokes 方程可用通用形式表示:

式中，ρ 是密度，t 是时间，U 是速度矢量.通用变量φ、耗散系数Γφ和通用源项Sφ表示为

式中，ui是i 方向的速度，μ 是黏度，p 是润滑液膜压力.

考虑到传统研究轴承空化现象的half-Sommerfeld 和Reynolds 边界方法无法分析润滑介质液相与气相之间的转变机制，文中采用一种基于压力变化的相变模型描述空化现象，Kunz 等［17］的研究认为该相变模型能够准确分析空化现象的气液转变机理.当润滑液膜压力小于其汽化压力即饱和蒸汽压pSat时相变发生，空穴产生，润滑介质由液相转变成气相，反之气相转变成液相.相变时气相、液相之间的质量传输率m 为

式中，上标+和-分别表征气相向液相转变和液相向气相转变，α 为液相体积分数，ρl和ρv分别为液相和气相密度，Cc、Cv、Σ∞和θ∞为经验常数.

液相体积分数α 的控制方程为

相变时控制方程中润滑介质物性参数采用其混合相的物性参数χ 进行表征:

式中，χl和χv分别表征为液相和气相物性参数.

滑移长度模型认为滑移速度与壁面处剪切率和滑移长度成正比，如式(5)所示.滑移长度b 是一个虚拟长度，定义为壁面处速度线性递减为0 时所在的位置到该壁面的距离，如图1 所示.极限剪切应力模型认为只有当壁面处流体剪切应力达到某一临界值时滑移才发生，且滑移时壁面处剪切应力等于该临界剪切应力值，如式(6)所示.

式中:Us是滑移速度;b 是滑移长度;u/y 是壁面剪切率;s是壁面剪切应力;0是壁面原始剪切应力，且 0 =μ(u/y);c是发生滑移的临界剪切应力.

图1 滑移长度示意图Fig.1 Schematic diagram of slip length

结合上述两种滑移模型，考虑到滑移速度Us与临近固液界面的流体速度Uin密切相关，文中通过引入两个滑移度系数f1和f2将壁面滑移速度与临近壁面流体速度Uin和壁面法向向量n 建立相应关系:

如图2 所示，若临近壁面流体速度Uin在以壁面法向方向为负zn方向的坐标系xn-yn-zn下的3 个速度分量为(uin，vin，win)，滑移速度Us在该坐标系下的3 个速度分量为(us，vs，ws)，则该双系数滑移数值边界表示为:

图2 滑移速度与近壁面速度示意图Fig.2 Schematic diagram of slip velocity and velocity near wall

滑移度系数f1和f2的大小体现了织构滑移表面流体发生速度滑移的强弱程度，与织构滑移表面的特性密切相关.f1=f2表示壁面速度只有切向方向的速度分量，壁面法向方向的速度分量为0.当f1=f2=0 时，该双系数数值边界表征为无滑移静止壁面.当f1=1，f2=0 时，该双系数滑移数值边界与临界剪切应力值为0 时的极限剪切应力模型等效.在开源计算流体动力学程序包OpenFOAM 中编写该双系数数值滑移边界，植入到计算程序中，使之能施加到织构表面.

2 方法验证

为验证该双系数数值滑移边界的正确性，首先以文献［16］中的轴承为对象进行研究.该止推轴承含有12 个长17 mm、宽16.5 mm 的轴承块，轴承块上加工有2150(43 ×50)个直径100 μm、平均深度25 μm的圆形凹坑织构，凹坑中心距200 μm，如图3 所示.织构区域与轴承块宽度方向边距一侧为4.1 mm、另一侧为2.4 mm;织构区域与轴承块长度方向边距为8.4 mm，详细参数见文献［16］的正文描述.润滑油密度为848 kg/m3，黏度为0.022 Pa·s.采用提出的双系数数值滑移边界表征滑块织构表面，将该边界施加到织构所在的整个区域.轴承承载能力在实验时是作为恒定外力(150 N)加载的，而润滑液膜厚度是实验结果;在理论分析时，实验获得的润滑液膜厚度是作为已知量用于构建数值模型，而承载能力是理论分析的结果值.理论分析结果与文献的实验结果比较如表1 所示.300 r/min 时，f1=f2=0.206 5;400 r/min 时，f1=f2=0.237 5;500 r/min 时，f1=f2=0.2882.3 种转速下，摩擦扭矩理论与实验的最大相对误差为4.8%;理论分析得到的承载能力与实验时加载的150 N 恒定载荷最大偏差为0.05%.理论分析结果与实验值吻合良好，表明将该双系数数值滑移边界施加到织构表面表征织构表面的综合效果可以分析织构表面对轴承性能的影响.

图3 文献［16］的轴承模型Fig.3 Bearing model in reference［16］

表1 理论值与实验值比较Table 1 Comparison between theoretical and experimental results

3 滑动轴承模型

为给滑动轴承微织构的优化设计提供理论指导，文中采用上述双系数数值滑移边界对圆柱形滑动轴承轴瓦内壁面上的织构区域的位置分布、面积大小等宏观参数对轴承承载能力和空化的影响规律进行研究.滑动轴承模型如图4 所示，轴承宽度B =25 mm，转轴半径r=25 mm，轴承半径间隙c =50 μm，偏心率e=0.5.将润滑水膜周向展开，水膜最厚位于90°处，最薄位于270°处.M 表示织构区域在周向方向的起始位置，l 表示织构区域的周向长度，q 表示织构区域在轴向方向的起始位置，s 表示织构区域的轴向宽度.转速ω=1000rad/s(即9550r/min).润滑介质为水，密度ρ = 998.2 kg/m3，黏度μ =1.003 ×10-3Pa·s.无量纲压力P*=Pc2/(μrωB)，无量纲承载能力W*=Wc2/(μωBr3)，其中承载能力pdS，p 为压力，dS 为微元面积.假设只在轴瓦内表面的部分区域加工有微织构，轴瓦未加工织构的表面和轴颈表面为普通的无滑移表面，轴向两端面处润滑液膜压力等于环境压力.

图4 滑动轴承模型Fig.4 Journal bearing model

4 数值结果分析与讨论

为保证数值计算的精度和适当的计算时间，进行了网格独立解分析.30 ℃时纯水的汽化压力即饱和蒸汽压约为4053 Pa.选用无织构轴承模型构建了8 套网格系统.轴承承载能力随网格数变化如图5所示，当网格数超过2 ×105时随着网格数增加承载能力的变化很小，通过衡量计算精度和计算时间最终选用节点数目为2 ×105的网格系统作为后续模型的计算网格.

图5 网格大小对承载能力的影响Fig.5 Effects of node number on load-carrying capacity

为重点研究织构区域的分布位置和面积对轴承摩擦学性能的影响规律，后续分析取滑移度系数f1=f2=1.在分析织构区域在周向方向的分布位置和面积大小的影响时，织构区域在轴向方向上的分布保持不变，其轴向宽度等于轴承宽度(s =B)，汽化压力为101325 Pa.研究结果表明，当织构区域位于空穴区时，表面织构对润滑液膜压力和轴承承载能力没有影响或者影响很小.不同织构区域工况下，润滑液膜中心截线压力分布如图6 所示.由图6(a)可见，无织构时，流体流经收敛间隙产生的流体动压力的最大压力峰位于237°处，空穴主要分布在0°～198°和280°～360°两个区域.当织构区域位于280°～300°时，其压力和空穴分布与无织构时的分布一样;两种模型的承载能力也相差很小(见表2)，在数值误差允许的范围内可以认为承载能力也相同.当润滑液膜压力低于汽化压力时，润滑介质由液相变为气相，空穴产生，此时空穴区压力恒等于汽化压力，因此当织构位于空穴区域时不会对液膜压力产生影响.

由图6(a)、(c)可见，织构区域位于200°～220°和位于180°～200°的最大压力峰也位于237°处，与无织构时一样，这个最大的压力峰是流体流经收敛间隙时产生的流体动压力.织构区域位于180°～200°和位于200°～220°时，还存在另外一个较小的压力峰，分别位于200°和220°处.第2 个较小的压力峰是由表面织构导致的流体动压，其位置正好位于表面织构区域末端(即沿流体运动方向的下游端).此外，在表面织构区域前端(亦即上游端)会产生一个低压区.如果表面织构位于流体压力上升区域(p/θ ＞0，θ 为周向方向的角度)，收敛间隙和表面织构产生的两种流体动压相互叠加和互相增强，提高了液膜最大压力和轴承承载能力;而如果表面织构位于流体压力下降区域(p/θ ＜0)，织构区域上游端的低压效应(或称为负动压效应)会将收敛间隙导致的动压区域截断，润滑液膜破裂，加剧空化现象，此时虽然仍存在两个压力峰，但是表面织构的影响是负面的，将降低轴承承载能力.由图6(b)可见，当织构区域位于240°～260°时，两个压力峰分别位于225°和260°;当织构区域位于260°～270°时，两个压力峰分别位于234°和270°;这两种情况的承载能力均小于无织构时的承载能力(见表2).

图6 不同织构区域工况下的压力分布Fig.6 Pressure distribution at different texture regions

如果表面织构区域末端位于流体压力上升区域前端的油膜破裂位置附近，织构区域末端的动压效应能起作用，压力分布曲线在织构区域末端存在一个较小的压力峰，同时又增强了流体流经收敛间隙产生的动压效应，承载能力增加.由图6(c)可见，当织构区域位于180°～200°时，在织构区域末端即200°处存在一个较小的压力峰;此时仍有部分织构位于流体压力上升区域，所以其承载能力大于无织构时的承载能力(见表2).如果表面织构区域前端位于流体压力上升区域前端的油膜破裂位置附近，织构区域前端的低压效应起作用，空穴面积增加，同时织构区域末端的动压效应也起作用，压力分布曲线在织构区域末端存在一个较小的压力峰，同时又增强了流经收敛间隙产生的动压效应，承载能力也增加，大于无织构时的承载能力.由图6(a)可见，当织构区域位于200°～220°时，空穴面积增加了，油膜破裂位置移到了210°附近;同时在织构区域末端即220°处存在一个较小的压力峰，其承载能力大于无织构时的承载能力(见表2).虽然织构区域位于200°～220°时的液膜最高压力大于织构区域位于180°～200°时的最高压力，但是由于织构区域位于200°～220°时空穴面积增加了，导致其承载能力小于180°～200°时的承载能力.如果表面织构正好位于流体压力下降区域的末端，因为靠近空穴区，织构区域前端的低压效应能起作用，空穴面积增加，而织构区域末端的动压效应将不起作用，整体压力分布曲线只有一个压力峰，承载能力降低.由图6(b)可见，当织构区域位于270°～280°时，压力分布曲线只有一个压力峰，空穴分别位于0°～198°和270°～360°两个区域，相比无织构工况空穴面积变大，承载能力变小(见表2).

表2 不同织构区域工况下的轴承承载能力Table 2 Load-carrying capacity of bearing at different texture regions

由图6(a)、(c)可见，如果收敛间隙所导致的流体动压的压力峰正好位于织构区域内，两个压力峰将会合并成一个.当织构区域位于220°～240°和180°～240°时，单一的压力峰均位于240°;当织构区域位于180°～260°时，单一的压力峰位于260°.由表2 可见，保持织构区域在周向方向上的起始位置不变(位于180°)，当织构区域的末端从200°延伸到240°时承载能力增加了18.4%，表明承载能力随着位于流体压力上升区域的织构表面的面积的增大而增大;表面织构位于180°～260°时的承载能力小于表面织构位于180°～240°时的承载能力，但仍大于无织构时的承载能力，表明如果位于流体压力上升区域的织构表面的面积大于位于流体压力下降区域的织构表面的面积，轴承的承载能力仍将大于无织构时的承载能力.

5 结语

文中提出双系数数值滑移边界表征表面织构宏微观相互作用的综合效果并施加到有微织构的轴承表面，采用基于压力变化的相变模型描述空化现象，分析了织构区域位置分布、面积大小等参数对滑动轴承承载能力和空穴的影响规律.首先采用文献实验轴承进行分析，理论结果与文献结果吻合良好，验证了该双系数数值滑移边界和分析方法的正确性.研究表明，织构位于空穴区域时不会影响润滑液膜压力和轴承承载能力;织构区域前端会存在一个低压区(低压效应)，末端存在一个高压区(动压效应);当织构位于压力上升区域时，收敛间隙所产生的动压效应与织构产生的动压效应相互叠加，互相增强，可以提高承载能力，且该区域织构表面面积越大，承载能力也越大;但当织构位于压力下降区域时，织构前端低压效应将截断收敛间隙所产生的动压效应，导致润滑液膜破裂加剧空化现象并降低承载能力.合理的表面织构设计可以显著提高轴承摩擦学性能.

此外，文中如此处理织构区域仅是数值上的织构表面，在实际进行织构设计时为使轴承摩擦学性能最优，相关几何参数可以参考他人重点研究织构几何尺寸影响的相关成果，以使文中的结论与相关成果互相补充.

［1］Tønder K.Inlet roughness tribodevices:dynamic coefficients and leakage［J］.Tribology International，2001，34(12):847-852.

［2］Cupillard S，Cervantes M J，Glavatskih S.Pressure buildup mechanism in a textured inlet of a hydrodynamic contact［J］.Journal of Tribology，2008，130(2):021701/1-10.

［3］Etsion I，Burstein L.A model for mechanical seals with regular microsurface structure ［J］.Tribology Transactions，1996，39(3):677-683.

［4］Siripuram R B，Stephens L S.Effect of deterministic asperity geometry on hydrodynamic lubrication［J］.Journal of Tribology，2004，126(3):527-534.

［5］Fowell M，Olver A V，Gosman A D，et al.Entrainment and inlet suction:two mechanisms of hydrodynamic lubrication in textured bearings［J］.Journal of Tribology，2007，129(2):336-347.

［6］Dobrica M B，Fillon M，Pascovici M D，et al.Optimizing surface texture for hydrodynamic lubricated contacts using a mass-conserving numerical approach ［J］.Proceedings of the Institution of Mechanical Engineers，Part J:Journal of Engineering Tribology，2010，224(8):737-750.

［7］Sahlin F，Glavatskih S B，Almqvist T，et al.Two-dimensional CFD-analysis of micro-patterned surfaces in hydrodynamic lubrication［J］.Journal of Tribology，2005，127(1):96-102.

［8］Choi C H，Westin K，Breuer K S.Apparent slip flows in hydrophilic and hydrophobic microchannels［J］.Physics of Fluids，2003，15(10):2897-2902.

［9］Tretheway D C，Meinhart C D.Apparent fluid slip at hydrophobic microchannel walls ［J］.Physics of Fluids，2002，14(3):L9-L12.

［10］Zhu Y，Granick S.Limits of the hydrodynamic no-slip boundary condition［J］.Physical Review Letters，2002，88(10):106102/1-4.

［11］Pit R，Hervet H，Leger L.Direct experimental evidence of slip in hexadecane:solid interfaces［J］.Physical Review Letters，2000，85(5):980-983.

［12］Salant R F，Fortier A E.Numerical analysis of a slider bearing with a heterogeneous slip/no-slip surface［J］.Tribology Transactions，2004，47(3):328-334.

［13］Spikes H A.The half-wetted bearing.part 1:extended Reynolds equation［J］.Proceedings of the Institution of Mechanical Engineers，Part J:Journal of Engineering Tribology，2003，217(1):1-14.

［14］Rao T.Analysis of single-grooved slider and journal bearing with partial slip surface［J］.Journal of Tribology，2010，132(1):014501/1-7.

［15］Fatu A，Maspeyrot P，Hajjam M.Wall slip effects in(elasto)hydrodynamic journal bearings［J］.Tribology International，2011，44(7/8):868-877.

［16］Marian V G，Gabriel D，Knoll G，et al.Theoretical and experimental analysis of a laser textured thrust bearing［J］.Tribology Letters，2011，44(3):335-343.

［17］Kunz R F，Boger D A，Stinebring D R，et al.A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction［J］.Computers ＆ Fluids，2000，29(8):849-875.