基于动压效应的柱塞泵滑靴副润滑特性仿真分析

王月鹏， 陈先勇， 曹华威， 宋 博， 李 嘉

(1.国网北京昌平供电公司， 北京 102200； 2.许继三铃专用汽车有限公司， 河南 许昌 461000；3.许继时代技术有限公司， 河南 许昌 461000；4.长安大学 公路养护装备国家工程实验室， 陕西 西安 710064)

引言

柱塞泵是重要的动力装置之一，是能够连续地提供不同流量和压力的介质。相比其他类型的油泵，柱塞泵转速范围大、增压能力强、可控性高[1]。因此，高性能柱塞泵的性能分析研究一直是被关注的重点。其中，柱塞泵中的几个主要摩擦副是柱塞泵性能和寿命的重要影响因素。而滑靴副长期工作在高速、重载的状态下，若其油膜太厚，密封性遭到破坏，泄漏增加，将大大降低泵的容积效率；若其油膜太薄或者无法形成，则滑靴副表面会发生磨损甚至烧坏，缩短泵的使用寿命。所以滑靴副的性能更直接制约柱塞泵高压化、高速化的技术发展。因此，要实现柱塞泵大流量、高压化、高速化、低噪声和长寿命等目标，掌握滑靴副的润滑特性，建立优异的润滑油膜是重要的技术手段之一[2]。然而，目前国内柱塞泵多为测绘仿制，且滑靴副大都按静压支承原理进行设计与分析。由于柱塞泵处于周期性交变重载的工作条件下，滑靴在斜盘上的运动轨迹是非规则的空间曲线，实际中滑靴副底面的楔形油膜因为相对运动将产生动压支承力，从而使滑靴副的压力场分布变得更加复杂[3]。所以，为精确反映滑靴副的实际润滑情况，有必要专门对柱塞泵滑靴副润滑特性进行研究。

1 国内外研究现状

目前，针对柱塞泵滑靴副润滑特性的相关研究受到了越来越多的重视。国外对于柱塞泵滑靴副的相关研究可追溯到20世纪70年代，主要代表是KOC E等[4-5]，确定了滑靴受到的离心力矩和球铰副摩擦力矩是导致滑靴发生倾覆现象的2个重要因素，会影响滑靴副油膜润滑状态。同时，MANRING N D等[6-7]在发展初期，在滑靴副的研究中也取得了一系列的研究成果，理论分析并建立了滑靴的动力学模型，根据受力情况对滑靴的倾覆趋势进行了预测，但忽略了油膜的挤压效应。由于柱塞泵工作过程中，滑靴侧倾可能与斜盘发生接触，并在高压下产生弹性变形。因此，BERGADA J M等[8-9]对滑靴副油膜可能发生的弹性混合润滑情况进行了研究。另一方面，随着数值分析方法和计算机仿真技术的兴起和发展，出现了采取数值计算和软件仿真相结合的方法进行柱塞泵关键摩擦副的研究工作及成果。 Purde大学基于C++设计了仿真软件，IVANTYSYNOVA M和WIECZ-OREK U等[10-11]应用该软件对轴向柱塞泵的配流副、柱塞副、滑靴副3处油膜进行了全面的润滑特性仿真研究。PELOSI M等[12]建立了滑靴副的热弹流动态仿真模型，对高压下的滑靴变形、油膜的温度场进行了数学建模和仿真研究，该模型涵盖了滑靴的动力学特性、压力场模型、滑靴和斜盘的变形以及滑靴副的热传递模型。

国内对于柱塞泵滑靴副的研究起步较晚，加上试验手段和条件匮乏，使得相关研究工作开展较少，主要集中在对柱塞泵滑靴副的润滑问题及影响规律方面[13-14]。郑家锦等[15]提出了一种多油室滑靴结构，并在滑靴表面开设梯形油槽，通过试验研究表明，这种结构可以减小滑靴的侧倾。胡新华[16]对柱塞副和滑靴副的运动受力进行了分析，计算了滑靴受到的倾覆力矩，并通过搭建试验台对油膜厚度进行了测量；此外，还提出了一种静压支承球铰副的设计方法[17]，并在考虑油膜挤压效应的情况下求解了滑靴油膜的静压支承力。孙营辉[18]针对滑靴副油膜厚度场、压力场以及滑靴自旋转速的测试，设计了试验台测试泵的机械结构，实现了厚度场与压力场的多点测量。于思淼[19]重点研究了滑靴副油膜的动压效应和挤压效应，应用MATLAB编制了滑靴副流体动压支撑分析软件。

综上可知，国内对滑靴副的研究由于起步晚相对较为局限，主要集中在滑靴副动力学特性以及润滑特性的理论研究。且上述研究对象基本都是平行柱塞的轴向柱塞泵滑靴副，而轴向倾斜式的球面斜盘轴向柱塞泵的滑靴副相关研究现阶段还处于对滑靴静压润滑特性的求解方面。对于动压效应、滑靴变形以及油膜温度场等关键问题还没有进行深入研究。为此，进行动压效应下的柱塞泵润滑特性仿真研究，对提高柱塞泵的性能具有一定的理论及应用价值。

因此，在某型带滑靴且柱塞倾斜的球面斜盘柱塞泵的基本参数基础上，进行滑靴副润滑特性的相关分析研究；其次，给出该型柱塞泵滑靴副油膜的润滑特性求解方法；进而，利用该求解方法进行动压效应的仿真分析，并基于支撑力计算验证方法的有效性；最后，对影响油膜压力分布的主要因素进行仿真分析研究，以确定影响规律。

2 柱塞泵参数

以某型带滑靴且柱塞倾斜的球面斜盘柱塞泵为对象，进行滑靴副润滑特性的相关分析研究。该型柱塞泵的基本元件有：转子、分油盘、柱塞(9个)、滑靴、斜盘和控制活塞等。其中，滑靴副主要由柱塞、滑靴、斜盘3个部件组成，为保证滑靴和斜盘间润滑为全膜润滑(即纯液体摩擦)，柱塞泵的滑靴副为带中心油池的阻尼孔型静压支承滑靴，斜盘为球面斜盘。因斜盘工作表面为凹球面，为保证良好贴合，滑靴底面也为近似球面。此外，由于斜盘表面的曲率半径远大于滑靴副的相关结构尺寸，为研究方便，在接下来的建模当中，假设滑靴底面和斜盘表面均为圆形平面。图1所示为该型柱塞泵滑靴副结构示意图，泵转子、柱塞、滑靴以及保持架等部分结构参数如表1所示。为了更好的展示其结构，进行了三维模型建模，模型示意图如图2所示。

图1 柱塞泵滑靴副结构示意图

图2 柱塞泵滑靴副三维模型示意图

该型柱塞泵燃油介质采用RP-3型航空煤油，其部分性能参数如下：

(1) 最大转速n： 4620 r·min-1

(2) 体积流量V： 400～10000 L·h-1

(3) 进口压力pin： 0.9 MPa

(4) 增压Δp： 20 MPa

(5) 功率P： ≤5.6 kW

此处由于本研究重点及篇幅限制，暂不给出各个参数的具体计算过程。

表1 部分结构参数

3 柱塞泵滑靴副动压效应

滑靴底面油膜产生的静压力不足以平衡柱塞腔的油压力，需将中心弹簧力、惯性力、离心力等其他力考虑进来。即静压支承滑靴稳定工作时的油膜支承力并不完全由静压支承力提供，滑靴副间还存在着动压效应，相应地存在动压支承力。滑靴在工作过程中受到离心力矩和摩擦力矩的作用，将使滑靴发生倾覆，从而在滑靴底面和斜盘表面形成一定夹角，产生楔形油膜，如图3所示。由流体力学动压润滑理论可知，滑靴与斜盘之间形成的楔形油膜由于相对运动会产生附加压力场，这是因为流体通过收敛间隙速度增大，产生额外的动压力，从而改变滑靴副油膜静压支承下的压力场分布，使整个滑靴副的润滑情况变得更为复杂。故燃油柱塞泵的静压支承滑靴实际上是靠动压和静压协调配合来实现有效润滑。

图3 柱塞泵滑靴副动压示意图

滑靴副工作过程中，存在一系列的润滑特性参数，滑靴副的实际润滑效果正是这些润滑特性参数共同作用的结果，主要包括承载力、抗倾覆力矩、泄漏量等，通常对相应的微观特性进行积分即可求得。

油膜承载力由下式进行计算：

(1)

由于油膜动压效应的作用，滑靴底面将形成不均衡的压力分布，从而分别产生绕x轴和绕y轴的抗倾覆力矩，表达式如下：

(2)

燃油柱塞泵工作时，滑靴副始终保持着一定的油膜间隙，油液从中心油池经滑靴密封带泄漏到转子腔，泄漏量的计算公式为：

(3)

泄漏量是滑靴副关键的润滑特性参数，若泄漏量过大，则会使整个柱塞泵的容积效率降低；若要减小泄漏量，则需要使油膜厚度变小，但会影响滑靴副润滑效果。

4 基于动压效应的柱塞泵滑靴副润滑特性求解

由于动压效应的存在，滑靴副油膜的润滑状况变得错综复杂，其压力场和厚度场的求解也变得更加困难。此处给出考虑动压效应的滑靴副润滑特性数值求解方法，分别完成油膜压力场、油膜厚度场的求解。

4.1 假设条件

由于滑靴副间的流动是错综复杂的空间流动，基于滑靴副基本润滑性质分析，研究中对润滑油膜作如下假设：

(1) 忽略油膜受到的体积力和惯性力，如重力和离心力；

(2) 介质为Newton流体；

(3) 不考虑油膜中的湍流和涡流；

(4) 油液在滑靴副壁面无滑移；

(5) 在润滑油膜厚度方向上，流体黏性和压力保持不变；

(6) 忽略油膜曲率引起的速度方向的变化；

(7) 假设油膜温度场均匀分布。

由于滑靴副油膜的厚度非常小，设定滑靴副油膜为层流，流体黏性和压力在厚度方向保持不变等假设是合理的。在层流状态下，选取柱坐标系为润滑油膜参考坐标系，由于滑靴倾斜方向与运动速度方向并不一致，所以在滑靴径向和周向都有动压产生。根据N-S方程和连续性方程可推导得出柱坐标系下的适合油膜描述的控制方程为：

(4)

式中，vr,vθ为滑靴与斜盘相对滑动速度的径向和周向分量;ht,hb为滑靴顶部和斜盘底部壁面的厚度。方程式等号后前两项表示油膜在滑靴与斜盘楔形间隙下运动产生的动压效应；中间两项表示由滑靴表面的边界速度梯度引起的变化；后两项表示滑靴和斜盘的弹性变形引起的油液挤压效应。

4.2 滑靴副油膜压力场求解

选取有限体积法求解油膜压力场，对滑靴副油膜的控制方程式(4)进行偏微分方程离散求解。由简化N-S方程，即式(5)，积分得到油膜流速方程，沿膜厚方向对流速进行积分得到单位宽度截面流量，结合控制体积宽度得到整个体积流量，最后通过流量守恒方程进行求解。滑靴副油膜压力场求解流程如图4所示。

具体过程如下：

(5)

式(5)为不可压缩流体的简化N-S方程，将其转化成圆柱坐标形式并在平面上积分可得油膜的流速为：

(6)

图4 滑靴副油膜压力场求解流程

将式(6)沿油膜厚度方向积分，以确定单位宽度截面上通过的体积流量，为：

(7)

式中，qr为体积流量的径向分量;qθ为体积流量的周向分量。

对滑靴副求解域进行径向和周向网格划分，即将滑靴副润滑区域划分成有限个控制体积，并使每一个体积单元充分微小，以保证求解过程中的迭代精度。以每一个网格节点为中心，选取半网格范围空间区域作为离散控制体积。根据流量守恒原理，从离散控制体积各边界流入和流出的流量相等，整个体积内的流量增量为0，即：

(8)

利用式(7)可求出式(8)中的各项流量，分别如式(9)～式(12)所示：

(9)

(10)

(11)

(12)

将式(9)～式(12)带入式(8)中，对方程中压力p的偏微分项进行向后差分处理，可得到各节点p(i, j)的数值表达式，进而构成了一个非线性方程组，方程组中的系数是各节点待求的油膜压力的函数。采取数值迭代的方法求解这一非线性方程组，从而得到整个油膜的压力分布。求解中先给定初始的中心油膜厚度hc和滑靴最大倾斜角βmax，采用Gauss-Seidel超松弛迭代算法进行求解，并设定数值收敛条件为:

(13)

式中，m和n分别为径向和周向网格节点数；ε为从k次迭代到k+1次的相对误差，取[ε]=1×10-6。若迭代误差满足要求，油膜压力场求解完毕，通过积分可求得滑靴副油膜支承力、动压力矩、泄漏量等润滑特性参数。

4.3 滑靴副油膜厚度场求解

在滑靴底面与斜盘表面间形成一层楔形油膜，该滑靴副油膜厚度场是复杂的三维空间场，随着外部载荷、压力场分布以及滑靴姿态的改变而改变。滑靴副油膜厚度场求解流程如图5所示。

如图5可以看出，厚度场的求解过程中，需要结合压力场求解结果为初值，进而完成油膜各个力矩的计算，在判断状态后，进行循环迭代。迭代中，对中心油膜厚度和倾斜角进行更新。

为了更加准确地描述油膜的倾覆特性，通过滑靴倾斜角度和滑靴中心油膜厚度进行该楔形油膜的数学描述，示意图如图6所示。其中，需要对滑靴x轴、y轴2个方向的角位移进行求解，并计算出滑靴底面中心油膜厚度。由图6中的几何关系进行推导，可得厚度场的表达式为：

h=hc+r·sinβ·cos(Ω-φ)

(14)

式中，hc为中心油膜厚度;Ω为坐标系中任一点与x轴的夹角;φ为滑靴最大倾斜方位角。

图5 滑靴副油膜厚度场求解流程

图6 滑靴副油膜厚度计算示意图

滑靴受到的主要力矩为离心力矩、密封带油膜摩擦力矩以及油膜动压力矩。此外，还受到围绕柱塞球头相对转动受到球铰副油膜的摩擦力矩。

滑靴在x轴和y轴上受到的离心力矩为：

(15)

摩擦力矩为：

(16)

动压力矩为：

(17)

由于滑靴球窝与柱塞球头的相对运动幅度较小，假设球窝中心与球头中心重合，则球绞副摩擦力矩为：

MFx=-μbsFnRb

(18)

式中，μbs为摩擦系数，一般取常数。

对式(14)的厚度场进行求解，是完成滑靴倾斜角度和倾斜方位角，而滑靴倾斜角和倾斜方位角均由其绕x轴和y轴的转动角度φx和φy决定，满足如下关系：

(19)

假定球头和滑靴的接触点在球心和滑靴重心的连线上(z轴上)，则滑靴的运动速度由柱塞腔内柱塞的线性运动决定，为：

vz=vcosψ

(20)

式中，v为柱塞在柱塞腔内的相对运动速度。

对式(20)求导，得滑靴在z轴方向加速度，为：

(21)

根据受力平衡和力矩平衡条件，确定滑靴副动力学平衡模型，为：

(22)

式中，Jx为滑靴关于x轴的转动惯量;Jy为滑靴关于y轴的转动惯量。

从式(22)可知，滑靴副的动力学平衡模型是一个非线性方程组，运用Newton-Raphson方法进行迭代求解得出φx和φy，可得到整个滑靴副的厚度场。

5 基于动压效应的柱塞泵滑靴副润滑特性仿真

对柱塞泵滑靴副油膜的润滑特性进行数值模拟，分析其动压效应下的润滑特性，并给出了验证；同时，分析中心油膜厚度、滑靴最大倾斜角、转子转速以及进口压力对滑靴副油膜压力分布的影响。

5.1 动压效应仿真结果及验证

以上述求解方法进行考虑动压效应下的仿真求解。图7和图8所示为滑靴副1个周期内的中心油膜厚度和最大倾斜角计算结果。

图7 滑靴副中心油膜厚度

图8 滑靴最大倾斜角

图中，最大倾斜角是绝对值，没有考虑其与坐标系位置的关系。图7中，在考虑油膜动压效应的情况下，油膜厚度在α=180°状态(吸排区切换状态)时，油膜厚度不再发生跳变，这正是动压效应的作用。图8中，滑靴在α=180°状态(吸排区切换状态)时发生最大倾斜，滑靴在该点将产生最大动压支承力。

图9给出了在α=180°状态(吸排区切换状态)时，考虑动压效应的油膜压力分布。

图9 考虑动压效应的油膜压力分布

在该条件下，油膜动静压混合支承力可由式(23)对压力积分求得：

(23)

由静压润滑理论可知油膜的静压支承力为：

(24)

因此，油膜的动力支承力为：

Fd=Fhd-Fj

(25)

可求得该状态油膜动静压混合支承力，即：全部承载力为2509.5 N，静压支承力为1790 N，动压支承力为719.5 N。而滑靴副受力分析结果得到α=180°时滑靴所受合压紧力大约2532.5 N(此处受力理论分析不详细阐述)，因此，数值计算结果和理论分析结果所得到的动静压混合支承力几乎吻合。同时，证明了柱塞泵滑靴副油膜的承载力主要由静压支承力提供，但动压支承力体现的动压效应不可忽略。值得注意的是，有必要利用相关实验进行仿真的有效性验证，作为后续研究重点。

此外，如图10和图11所示为滑靴倾斜角β和滑靴旋转速度ω为0时的油膜压力分布 。

图10 β=0时的油膜压力分布

图11 ω=0时的油膜压力分布

图10和图11中，没有出现明显的动压峰值，进而计算求得动压支承力，β=0时为42.8 N，ω=0时为71.1 N，表明两种情况下滑靴副油膜中依旧存在少量的动压力，这是由于滑靴副流道中的介质本身流动，流体只要具有速度就会产生相应的动压力，通过间隙时介质流动速度增加，使动压力增大。

通过上述3种情况下的仿真结果分析可以得到，考虑动压效应时，数值计算结果和理论分析结果所得到的动静压混合支承力几乎吻合。同时，滑靴副油膜的动压力和静压力相对独立，相互影响较小，且在滑靴工作过程中静压作用更为稳定，在油膜承载力中占主导地位。而动压作用虽然更小，但动压支承力体现的动压效应不可忽略。

5.2 压力分布影响因素分析

滑靴副油膜形状主要由中心油膜厚度hc和滑靴最大倾斜角β确定，不同的油膜形状将使滑靴油膜压力分布产生巨大差异，影响其润滑特性。因此，逐一对上述两种影响因素进行仿真分析。

图12 不同中心油膜厚度下的压力分布仿真结果

首先，分析中心油膜厚度对油膜压力分布的影响规律。选择在滑靴最大倾斜角不变时，对不同中心油膜厚度下的滑靴副油膜压力进行仿真分析。如图12所示为β取0.00075 rad，hc分别取10，11，12 μm的油膜压力分布仿真结果。

由图12可知，当滑靴最大倾斜角保持不变时，不同的中心油膜厚度下，滑靴油膜的静压作用区域压力分布几乎不发生改变，动压力峰值随着中心油膜厚度的增加而减小。油膜整体的静压支承力保持不变，动压支承力减小，但静压力分布随着中心油膜厚度的增加其压力梯度向边缘发散。此外，油膜的动压效应对于中心油膜厚度的变化十分敏感，在最大倾斜半径上，10 μm中心膜厚下的最大压力可达40 MPa以上，1 μm的中心膜厚变化可能造成不少于20 MPa的动压力差值。同时，最大倾斜角固定时，中心油膜厚度的增加会使楔形油膜的收敛性减弱，从而抑制动压力的产生，使收敛间隙的动压效应减弱，相应的动压力峰值也随之减小。

其次，分析滑靴最大倾斜角对油膜压力分布的影响规律。选择在中心油膜厚度不变时，对不同滑靴最大倾斜角下的滑靴副油膜压力分布进行仿真分析。如图13所示为hc取11 μm，β分别取0.0006, 0.00075,0.0009 rad的油膜压力分布仿真结果。

图13 不同最大倾斜角下的压力分布仿真结果

由图13可知，滑靴中心油膜厚度固定时，随着滑靴最大倾斜角的增加，楔形油膜的收敛性增强，从而促进动压力的产生，使收敛间隙的动压效应更为显著。通过计算动压支承力可知，相比与中心油膜厚度，动压效应对于最大倾斜角的变化更为敏感。同时，在这一过程中油膜的静压作用区域压力分布几乎没有发生变化。

为了更为清晰的反应2个因素对油膜压力的影响规律，根据上述仿真结果，对2个因素的油膜动静压支承力进行计算。表2为不同中心油膜厚度的油膜动静压支承力结算结果，表3为不同最大倾斜角的油膜动静压支承力结算结果。

表3 不同最大倾斜角下的油膜动静压支承力

从表2及表3的计算结果来看，中心油膜厚度和最大倾斜角的变化能够影响油膜的动静压支承力，随着中心厚度增大，静压支撑力稍微增强，影响较小，而动压支撑力影响较大。随着最大倾斜角增大，静压支撑力稍微增强，而动压支撑力影响较大。总的来说，中心油膜厚度、滑靴最大倾斜角主要对动压效应产生影响，均是滑靴副油膜的压力分布的影响因素。

5.3 主要结构参数影响分析

通过改变滑靴底面工作半径，对其影响滑靴的润滑特性等规律进行分析。滑靴底面工作半径r2分别取9, 10, 11 mm，其他仿真参数保持一致，取hc=11 μm，β=0.00075 rad，进口压力为10.29 MPa，出口压力为0.3 MPa，转子转速为4620 r/min,中心油池半径固定为5 mm，得到滑靴底面工作半径对油膜压力分布的影响仿真结果如图14所示。

图14 滑靴底面工作半径对油膜压力分布的影响

由图14分析可知，滑靴底面工作半径对滑靴副油膜的压力分布影响很大，随着r2的增大，油膜的静压作用区和动压效应区的压力梯度都随之增大，整个油膜的承载力有了很大的提升。因此，在满足滑靴尺寸设计要求的前提下，通过增大滑靴底面工作半径的方式对滑靴结构进行优化，不仅可以提高滑靴副油膜整体的承载能力，还可以增强油膜的动压效应，从而达到提升滑靴抗倾覆能力的目的。

5.4 油膜温度分布影响分析

在求解中首先对雷诺方程进行求解获得压力分布，再带入积分方程进行求解得到温度的分布。图15为不同偏心率(e=0.1,0.4,0.6时)滑靴副的温度分布仿真结果。

可以看出，温度峰值主要有4个。分别位于轴承中心位置，最小油膜厚度的端部以及承载部位的前部。主要是因为温度分布与油膜厚度及压力梯度相关。从能量方程的关系来看，温度与膜厚成反比，所以在最大膜厚处温度最低，而在最小膜厚处温度最高。此外，与压力在周向方向的梯度的平方成正比处温度最高。温度还与压力在轴向方向的梯度平方成正比，因此温度会在宽度方向上呈现出中间低、两边高的分布。

图15 不同偏心率滑靴副的温度分布仿真结果

6 结论

为了更好的分析柱塞泵中滑靴副的润滑特性，提出了一种考虑动压效应的滑靴副润滑特性仿真求解方法。并基于该方法，对某型柱塞泵进行了润滑特性仿真分析研究。主要得到了以下结论：

(1) 考虑柱塞泵滑靴副中的动压效应，提出的滑靴副润滑特性的仿真求解方法，能够实现其压力场和厚度场的求解；仿真分析结果表明，数值计算结果和理论分析结果所得到的油膜动静压支承力几乎吻合，证明了该方法的有效性；

(2) 通过动压效应的分析结果表明，滑靴副油膜承受了动压、静压两种形式的混合支承力，动压力和静压力相对独立，相互影响较小，且在滑靴工作过程中静压作用更为稳定，在油膜承载力中占主导地位。而动压作用虽然更小，但动压支承力体现的动压效应不可忽略；

(3) 分析了滑靴副中心油膜厚度和最大倾斜角对油膜压力分布的影响规律。结果表明，中心油膜厚度、滑靴最大倾斜角是动压效应的主要影响因素。随着中心厚度增大，静压支撑力稍微增强，影响较小，而动压支撑力影响较大。随着最大倾斜角增大，静压支撑力稍微增强，而动压支撑力影响较大。此外，通过改变滑靴底面工作半径可以实现滑靴结构的优化。

所得出的结论对实现高性能柱塞泵的研制具有一定的工程实践意义。