计入柱塞套弹性变形的柱塞副摩擦与密封特性分析

2023-01-16闫康昊

中国机械工程 2023年1期

杭旸闫康昊黄丹

河海大学力学与材料学院，南京，211100

0 引言

高压油泵是共轨系统的重要液力部件。径向柱塞泵因结构紧凑、转速高、扭矩低等优点而被广泛应用于轻型车及乘用车的共轨系统中。在径向柱塞泵向高压化和高速化方向发展的同时，也对其寿命和可靠性提出了更高的要求。作为径向柱塞泵设计及优化的基础，其内部关键摩擦副的润滑与密封机理已成为国内外研究的热点。

高压油泵的工作条件较为恶劣，柱塞在缸体内往往处于倾斜状态。王慧敏等[1-2]通过相关参数对柱塞副微运动特性展开了研究，并获得了柱塞副瞬态泄漏特性。在计入微运动的基础上，QI等[3-4]对柱塞处于偏心状态下的油膜特性进行了分析，并从摩擦学的角度考虑了表面形貌效应对摩擦副润滑特性的影响，但未计入柱塞套弹性变形，无法准确分析摩擦功耗及泄漏量。目前，与上述有关的复杂条件下柱塞副摩擦学研究多见于轴向柱塞泵研究领域[5]。高彦军等[6]利用FLUENT软件对轴向柱塞泵运动特性进行了模拟，探究了油液特性对流量脉动的作用规律。NIE等[7]和王克龙等[8]均建立了柱塞副弹性流体润滑模型，但其计算模型均处于全液润滑状态。对于混合润滑，近年来国内外专家学者开展了大量研究。目前主流混合润滑模型分为统计学模型和确定性模型两大类[9]。统计学模型采用统计学参数表征微凸体对流动特性的影响，具有计算高效、求解方便的特点，应用较为广泛。ZHAO等[10]和SUN等[11]均基于平均雷诺方程对混合润滑状态下轴向柱塞泵润滑特性开展了研究，对柱塞副摩擦学机理分析作了重要补充。由于点、线接触问题计算耗时较长，所以目前确定性模型的主要应用领域为非共形表面问题，在柱塞副混合润滑研究中应用较少。上述研究为共轨高压油泵柱塞副摩擦与密封特性的分析提供了重要参考。鉴于径向柱塞泵与轴向柱塞泵的工作特性有很大区别，且计入表面形貌、柱塞套弹性变形的柱塞副摩擦与密封特性的研究仍较少见诸报道，因此对径向柱塞泵柱塞副开展摩擦学研究及相关参数影响规律的分析显得迫切而必要。

本文以某高压共轨径向柱塞泵为例，以计入柱塞套弹性变形的柱塞副混合润滑数值仿真模型(mixed-elastohydrodynamic，MEHD)为基础，重点研究在柱塞泵运行过程中柱塞与柱塞套间的摩擦与密封性能，并通过参数分析评价不同运行条件下柱塞套弹性变形及表面形貌对摩擦功耗泄漏量的影响，为高速高压径向柱塞泵的结构设计与优化提供参考。

1 数学力学模型

柱塞副作为高压共轨径向柱塞泵的关键摩擦副，对喷射系统的燃油喷射压力有很大影响。偏心轴是柱塞运动的重要驱轴，对柱塞副的运动特性、力学特性、可靠性影响显著。转速已定时，偏心轴做等速回转运动，柱塞的运动特性由偏心轴的轮廓曲线决定。柱塞偶件-偏心轴系统结构简图见图1，ω为凸轮运动的角速度。在凸轮的偏心运动下，柱塞在缸体内做往复运动。

图1 柱塞偶件-偏心轴系统Fig.1 Plunger coupled parts-eccentric shaft system

图2 柱塞副结构模型图Fig.2 Structural model diagram of the piston/ bushing interface

图3 柱塞副油膜模型Fig.3 Film model of the piston/bushing interface

考虑柱塞微运动及柱塞套弹性变形时，油膜厚度方程可表示为

(1)

式中，c为半径间隙；Δh为柱塞套弹性变形量。

为方便分析且不失一般性，假设密封流体为牛顿流体，沿膜厚方向的压力保持不变，流体流动为层流，忽略惯性效应的影响。在此条件下，可采用平均形式的Reynolds方程[12]来计算油膜压力：

(2)

除此之外，还需对Reynolds方程补充边界条件。边界条件的定义包含进/出口处的压力边界以及油膜展开处的连续性条件，具体如下：

对于油膜的入口、出口，分别有

(3)

在油膜展开处的左端和右端有

(4)

由于柱塞腔燃油压力较大，故它对黏度的影响不可忽略。本文采用Barus黏度-压力关系式[15]，即

μ=μ0exp(ap)

(5)

式中，μ0为环境黏度；a为黏压系数，通常取a=2×10-8Pa-1。

为了继续研究表面形貌对摩擦特性的影响，本文以高斯粗糙面来模拟表面形貌，建立粗糙面时采用指数形式的自相关函数：

(6)

参考陈辉等[16]使用AR模型的二维数字滤波技术生成高斯各项同性粗糙表面。局部膜厚hT则由式(1)与两粗糙面叠加所得。

以h=3σ为不同润滑状态(全液润滑和混合润滑)的分界线，并定义膜厚比H′=h/σ。采用GREENWOOD等[17]提出的粗糙表面接触模型，柱塞和柱塞套之间的接触压力可通过下式计算：

(7)

式中，η为粗糙表面的峰元密度；β为峰元曲率半径；E′为柱塞副综合弹性模量；F5/2(H′)可参考文献[18]计算得到。

流体平均剪应力可由下式计算：

(8)

式中，φfp、φfs、φf为剪切力因子[12-13]。

微凸体接触剪应力可由下式计算：

τc=fpcv

(9)

式中，f为经验摩擦因数，通常取0.08。

由此，总摩擦力包含流体的黏性摩擦力和粗糙度微凸峰之间的接触摩擦力，可计算为

(10)

则摩擦功耗为

Pf=Ffv

(11)

考虑表面粗糙度的泄漏量可表示为[19]

(12)

Reynolds方程求解收敛判断条件取为

(13)

式中，k为迭代次数；m、n为周向和轴向网格划分数目；pi，j为第(i，j)点的油膜压力。

柱塞受力分析如图4所示。柱塞在缸体内满足如下平衡条件。

图4 柱塞受力分析示意图Fig.4 Schematic diagram of piston force analysis

由∑Fy=0得

Fsk=(Fp+Fk+G+Fa+Ff)/cosα

(14)

由∑Fx=0得

Fsksinα=Fpx

(15)

由∑M=0得

Fsklsinα=Mz

(16)

(17)

式中，Fsk为凸轮对柱塞的作用力；Fp为高压油液对柱塞的作用力；Fk为回程辅助弹簧作用力；G为柱塞重力；Fa为轴向惯性力；α为凸轮压力角；l为柱塞长度。

2 联合求解与仿真方法

2.1 有限元模型

由于柱塞套结构较为简单，故采用Solid185单元进行模拟，建立的有限元分析模型如图5所示。由于柱塞套装配于缸体内部，故对该圆筒模型外表面及内端面施加固定约束，内表面留缸长度段施加油膜压力作为荷载，其余部分则施加柱塞腔工作压力。由于流体域采用有限差分法离散，而固体域则采用有限单元法，故计算时必须保证节点之间的一一对应。固体域与流体域节点的映射关系如图5所示。

图5 有限元模型固体域节点与流体域节点映射关系Fig.5 Node mapping between fluid domain and solid domain in finite element model

通过MATLAB与ANSYS互相调用结果文件构建联合仿真模型，实现FDM(finite difference method)与FEM(finite element method)的耦合求解。凸轮转速1500 r/min、入口压力120 MPa时的油膜压力分布及其对应的柱塞套弹性变形分布如图6所示。由图6可知，在油膜压力峰值或柱塞与铜套接触的位置，会产生较大的弹性变形，且整体变形分布与油膜压力分布一致，体现了本文方法的可靠性。

(a)油膜压力

(b)弹性变形图6 油膜压力及弹性变形分布图Fig.6 Distribution of film pressure and elastic deformation

2.2 数值求解

本文基于Reynolds边界条件，采用有限差分法求解润滑方程式(2)，应用超松弛迭代法求解油膜压力离散方程，运用辛普森积分法计算柱塞副摩擦与密封特性参数。根据柱塞微运动状况、表面形貌参数、柱塞结构参数，由式(1)计算各节点油膜厚度，以式(13)作为收敛判断条件求解润滑方程得到柱塞副各节点油膜压力分布，以式(15)及式(16)判断柱塞受力是否平衡，若不平衡则继续搜索满足条件的倾斜角及偏移量。最后根据式(7)～(12)计算各摩擦与密封特性参数。同时，为便于进行参数分析，本文从MEHD模型分离出仅考虑柱塞套弹性变形的EHD(electrohydrodynamic)模型及仅考虑表面形貌的MHD(mixed hydrodynamic)模型。

为实现FDM与FEM的耦合计算，本文在求解雷诺方程迭代过程的每一步，存储油膜压力为.txt文件，供FEM方法调用。FEM计算得到径向节点位移后自动写入创建好的.txt文件，MATLAB读取并重组后对应叠加到油膜厚度上，更新的油膜厚度分布参与下一次迭代，如此循环，直至满足所有判定条件。数值求解流程示意图见图7。

图7 数值计算流程图Fig.7 Flow chart of numerical calculation

3 模型验证与算例分析

取典型柱塞副，应用本文提出的计入表面形貌和弹性变形的MEHD模型和联合仿真方法对柱塞副摩擦与密封特性进行求解与分析。模型计算参数如表1所示。若不进行特殊说明，本文算例均取表面纹理为各向同性、两表面粗糙度均为0.8 μm且凸轮转角为90°(即凸轮压力角最大，此时处于排油阶段)。

表1 MEHD模型计算参数Tab.1 Calculation parameters for MEHD model

3.1 网格无关性分析

由于油膜各个方向上都包含流动信息且在厚度方向上尺度小至几微米，而油膜周向和轴向的尺度为数十毫米，因此，仿真计算结果对网格质量较为敏感。为保证计算精度，同时又提高计算收敛速度、节省计算资源，合理的网格划分显得尤为重要。本文采用图8所示的4种网格划分方案，选取展开角为π时对应的轴向油膜压力分布进行网格无关性分析。如图8所示，油膜轴向压力分布在网格划分180×50以后已趋于稳定，且与其他相关文献相比这样的划分已相当严格，故本文后续分析中采取此方案进行网格划分。

图8 网格无关性分析Fig.8 Grid independence analysis

3.2 模型验证

(a)吸油状态油膜压力分布对比图

(b)排油状态油膜压力分布对比图图9 模型验证Fig.9 Model validation

3.3 工作条件对最大微凸体接触压力的影响

图10为凸轮转速为1500 r/min，入口压力从80 MPa增加至160 MPa时油膜厚度分布及膜厚比等值线分布图，图11为对应的油膜压力分布图。由图10可知，柱塞副间隙油膜在80 MPa时膜厚波动幅度为2.097～8.018 μm，在160 MPa时波动幅度为2.740～8.327 μm。入口压力增大时，柱塞副两端压力差增大，油膜动压压力增加，更易满足柱塞副力平衡条件。因此，当工作转速一定时，随着入口压力的增大，柱塞副油膜厚度波动幅度减小，偏心状态减弱。

(a)pin=80 MPa

(b)pin=120 MPa

(c)pin=160 MPa图10 不同入口压力下油膜厚度分布及膜厚比等值线分布Fig.10 Oil film thickness distribution and film thickness ratio isoline distribution under different inlet pressures

(a)pin=80 MPa

(b)pin=120 MPa

(c)pin=160 MPa图11 不同入口压力下油膜压力分布Fig.11 Oil film pressure distribution under different inlet pressures

由于柱塞内外端面圆心处于同一直线上，故两端的接触程度变化一致，只需通过最大微凸体接触压力即可评价其接触、磨损程度。图12所示为凸轮转速1500 r/min时不同计算模型下最大微凸体接触压力随入口压力的变化关系。由图12可知，最大微凸体接触压力随着入口压力的增大而减小，且MEHD结果始终小于MHD结果。由此可见，不考虑柱塞套弹性变形将高估微凸体接触压力，从而高估其磨损程度。同时，由两种模型计算结果的比值可以看出，随着入口压力的增大，弹性变形对微凸体接触压力的计算影响越来越大。此外，低压时增大入口压力对减少接触、改善润滑的效果要优于高压时的效果。

图12 最大微凸体接触压力随入口压力的变化关系Fig.12 Variation of maximum contact pressure of micro convex body under different inlet pressures

图13为入口压力为120 MPa，凸轮转速从1000 r/min增加至3000 r/min时油膜厚度分布及膜厚比等值线分布图，图14为对应的油膜压力分布图。凸轮转速为1000 r/min时，膜厚波动幅度为2.915～7.774 μm，而在3000 r/min时则达到了1.740～8.572 μm。转速的增加使得柱塞副需要产生更大的动压压力以达到力平衡，因此油膜形态变化更为剧烈，导致了易接触区域更小的油膜厚度，柱塞副的润滑情况也愈加恶劣。

(a)n=1000 r/min

(b)n=2000 r/min

(c)n=3000 r/min图13 不同凸轮转速下油膜厚度分布及膜厚比等值线分布Fig.13 Oil film thickness distribution and film thickness ratio isoline distribution under different cam speeds

(a)n=1000 r/min

(b)n=2000 r/min

(c)n=3000 r/min图14 不同凸轮转速下油膜压力分布Fig.14 Oil film pressure distribution under different cam speeds

图15所示为入口压力120 MPa时不同计算模型下最大微凸体接触压力随凸轮转速的变化关系。低速状态下，弹性变形对微凸体接触压力的影响较小。随着转速提高至3000 r/min，两种模型计算结果的比值达到了0.66，此时柱塞套弹性变形对接触状态的评估起着极为重要的作用。此外，由图15可知，MEHD模型计算得到的微凸体接触压力增长趋势更为平缓。

图15 最大微凸体接触压力随凸轮转速的变化关系Fig.15 Variation of maximum contact pressure of micro convex body under different cam speeds

3.4 工作条件对摩擦功耗的影响

摩擦功耗反映了能量的耗散，是评价系统机械效率的重要指标。图16所示为凸轮转速1500 r/min时不同计算模型下摩擦功耗随入口压力的变化关系。由式(8)可知，柱塞所受摩擦力是由压差流及剪切流两者共同作用产生的。当柱塞为排油状态时，压差流产生的摩擦力与压差方向相同，而剪切流产生的摩擦力则与柱塞运动方向相反。随着工作压力的增大，柱塞副油膜微元体间切应力的增加导致黏性摩擦功率损失的增大。因此，各模型计算所得的摩擦功耗均呈增长趋势。对比MEHD结果及EHD结果可知，微凸体接触增大了摩擦力，计入表面形貌使得润滑环境更为恶劣。对比MEHD结果及MHD结果可知，在相同工况下，考虑柱塞套弹性变形时，将得到更大的摩擦功耗。对比三种模型的计算结果可知，在计算摩擦功耗时弹性变形比表面形貌更不可忽视。

图16 摩擦功耗随入口压力的变化关系Fig.16 Variation of friction power loss under different inlet pressures

图17所示为入口压力120 MPa时不同计算模型下摩擦功耗随凸轮转速的变化关系。随着凸轮转速的增大，各模型计算所得的摩擦功耗也随之增大。其中，MEHD模型结果最大，其次为EHD模型，MHD模型结果最小。同样的，EHD模型计算所得结果与MEHD模型更接近，相差不超过5%。

图17 摩擦功耗随凸轮转速的变化关系Fig.17 Variation of friction power loss under different cam speeds

3.5 工作条件对泄漏量的影响

泄漏量是衡量柱塞副容积效率的重要指标。图18所示为凸轮转速1500 r/min时不同计算模型下泄漏量随入口压力的变化关系。在排油阶段，泄漏量由剪切流与压差流的流量相减而得。定义内泄漏为正，外泄漏为负。不考虑弹性变形时，随着入口压力的不断增大，向外的压差流量不断增大，虽柱塞倾斜程度亦发生变化，但对向内的剪切流量影响较小，故此时泄漏量不断减小，并表现为内泄漏。考虑弹性变形时，间隙的增大对剪切流量与压差流量均有不同程度的增大影响，并且随着入口压力增大，弹性变形对压差流量的增大影响更为显著，因此MEHD模型所得结果要小于MHD模型所得结果，且两者的差值越来越大。此外，还可发现MHD模型计算所得结果变化较缓，入口压力从80 MPa增长至160 MPa的过程中泄漏量仅减少了8.86%。EHD模型计算所得结果与MEHD模型更为接近，但仅考虑弹性变形将高估系统的泄漏量，带来不可忽视的误差。

图18 泄漏量随入口压力的变化关系Fig.18 Variation of leakage under different inlet pressures

图19所示为入口压力120MPa时不同计算模型下泄漏量随凸轮转速的变化关系。各模型计算所得的泄漏量随转速的提高而增大，表现为越来越大的内泄漏。随着凸轮转速的提高，向内的剪切流量增大，泄漏量显著增大。与MEHD模型相比，EHD模型所得结果略微偏大，这也表明表面粗糙度对压差流量的增大影响较为显著，因此MEHD模型得到较小的内泄漏值。同时还可发现存在临界转速vcr≈2300 r/min，当转速小于vcr时，MEHD模型结果小于MHD模型结果，转速大于vcr则相反。造成该现象主要是由于转速的提高使得弹性变形对剪切流的影响越来越大，计入弹性变形时外泄漏的增加愈加难以平衡内泄漏的增加，直至2500 r/min时，弹性变形对剪切流量的增大影响已大于压差流量，故此时得到与前面相反的结果。