滑片泵进出口压力对其内流特性影响的数值模拟

张群峰，闫盼盼，单建平，何万发

(1. 北京交通大学 土木建筑工程学院，北京100044;2. 中航工业哈尔滨东安发动机集团有限公司，黑龙江 哈尔滨150066)

0 引言

空化［1］是润滑油滑片泵工作过程中容易发生的一种物理现象。空化发生时，滑片泵的工作性能受到很大的影响，往往导致滑片泵工作寿命缩短和效率下降，产生振动和噪声。滑片泵内的空化主要为叶片表面空化和间隙空化［2］。叶片表面空化通常是由叶片背面的二次流引起的，间隙空化是润滑油经过叶片与定子壁面的径向间隙以及叶片轴向端面间隙引起局部流速增加，压力下降到润滑油的当地饱和蒸气压而形成的。目前，从公开发表的文献来看，研究人员［3-10］通过理论分析和经验公式对滑片泵工作性能和空化现象做了一定的研究。由于滑片泵叶片的工作面积是变化的，仅靠理论分析和经验公式难以对滑片泵的空化特性作进一步的研究。随着计算流体力学(CFD)的发展和计算机性能的提高，CFD 为进行滑片泵空化特性的研究提供了重要手段，已有研究人员对简化滑片泵模型(简化处理滑片头部和定子之间的间隙形状)进行了二维或三维的数值模拟［11-12］。

利用CFD 技术来模拟滑片泵的空化特性的关键问题是滑片泵在变容积过程中运动边界的处理以及空化区气液界面的处理。处理运动边界的方法通常有基于贴体网格的任意拉格朗日欧拉(ALE)方法［13］、重叠网格法［14］和浸入边界法［15］。由于浸入边界法存在边界上求解精度不高和存在体积误差问题，目前没有广泛使用，而ALE 方法和重叠网格法则得到了较为广泛的使用。空化区气液界面的处理主要涉及空化模型和气液界面的捕捉，目前在CFD软件中广泛采用的是流体体积(VOF)方法［16-18］和基 于Rayleigh-Plesset 方 程［1］的Rayleigh 空 化 模型［19-21］。本文通过自编程序与Star-cd 软件相结合，采用ALE 方法准确模拟滑片运动过程中的相对位置，精确捕捉滑片头部和定子之间的间隙形状。并利用该功能，结合Menter SST k-ω 湍流模型［22-23］和Rayleigh 空化模型对某一航空发动机滑片泵的内部流场进行数值模拟，通过计算结果和实验结果的比较确定合适的空化模型参数，分析滑片泵内的流动特性和空化产生、发展和溃灭的过程;研究滑片泵进、出口压力对滑片泵工作性能和空化特性的影响。

1 控制方程与数值方法

1.1 控制方程和空化模型

润滑油发生气化时，气液界面上的密度变化很大，而界面速度是连续的。为了保证数值计算的稳定，通过求解三维N-S 方程并采用VOF 方法和Rayleigh 空化模型可以模拟空化的产生和溃灭等物理过程。

液态润滑油和气化润滑油的连续方程分别为

式中:α 和β 分别为液态润滑油和气化润滑油的体积分数，满足α+β=1;Sβ为空化源项;ρl和ρv分别为液态润滑油和气化润滑油密度;t 为时间;vgj为运动网格速度;uj为流动速度分量;xj为节点坐标值;j=1，2，3.

液态润滑油和气化润滑油混合物的动量方程为

式中:ρ=ρlα+ρvβ;f 为液态润滑油和气化润滑油界面动量输运源项;为克罗内克算子，ui和uk为流动速度分量，xk为节点坐标值，i=1，2，3.

Sβ通过Rayleigh 空化模型来求解，该模型中，空化核的增长率利用Rayleigh-Plesset 方程来描述，该方程为

式中:R 为空化核半径;pv为当地饱和蒸气压;σ 为表面张力系数;p 为当地静压。令n 为单位液体体积内空化核数，对上述方程简化后，可得

由以上各式可以看出，影响空化产生和空化区大小的因素主要是液体的饱和蒸气压、当地静压、空化核数和空化核的初始半径。其中空化核数和气泡初始半径是需要人为给定的。研究表明，空化核的初始半径对空化区域的大小影响较小，而空化核数目的大小对空化区域的大小影响较大，计算中需要通过实验结果来确定。

模拟滑片泵的变容积工作过程时，网格发生变形或重构，为了避免控制方程在离散过程中产生人工质量源，应满足空间守恒定律(SCL)［24］:

式中:nf为网格单元的表面数目;Sk为网格单元的第k 个面;Sj为表面Sk的面积矢分量。

1.2 湍流模型

式中:k 为湍流动能;μt为湍流粘性系数，采用Menter SST k-ω 湍流模型来计算湍流粘性系数。SST k-ω湍流模型在近壁处采用Wilcox k-ω 模型［25］，在边界层边缘和自由剪切层采用k-ε 模式，这两个区域之间通过混合函数来过渡，因此SST k-ω 湍流模型兼有近壁面Wilcox k-ω 模型的稳定性及边界层外部k-ε模型高效性的优点。具体形式为

1.3 计算方法

采用有限体积法在同位网格中对控制方程进行离散化，动量方程在空间上离散采用基于MUSCL格式［26-27］的二阶精度差分格式。模拟空化和转子转动时，计算采用非稳态计算，用PISO 方法［28］求解离散方程。引入空化模型后，计算容易发散，为了提高计算的收敛速度和稳定性，计算初始时，不激活空化模型进行单相流计算，迭代到流场建立后，再激活空化模型进行计算。

2 网格生成与运动网格的实现

选用的滑片泵转轴直径为30.154 mm，叶片数为4 片，叶片转动90°为一周期，因受安装位置的限制，滑片泵的进口流道不是直接对着吸油腔，而是与吸油腔在转动方向上近似呈90°夹角，图1 为滑片泵的结构剖面图。采用ICEM-CFD 软件生成滑片泵进、出口流道、定子内流道的六面体计算网格，因滑片泵内部流道较为复杂，为了保证网格质量，分块生成计算网格，各块网格交接面非连续时，采用Couple 进行连接以保证各部分网格间的连通性，如图2 所示。通过自编程序，根据定子轮廓线的曲线函数和滑片的运动规律，构建滑片泵转动区域的内表面，采用代数网格生成法生成运动区域的网格并进行网格的光顺处理，网格总数为57 万左右。图2 局部放大图为端面间隙网格，这部分网格的分布随着阀片的转动而变化。图3为转动区域叶片附近的网格，自编网格生成程序可以实现叶片头部与定子壁面间隙和端面间隙的网格生成，同时能保证间隙附近网格有较好的正交性。

图1 结构剖面图Fig.1 Cross-sectional view of vane pump

图2 滑片泵计算网格Fig.2 Computational mesh of vane pump

计算时，根据物理时间确定滑片泵转轴的位置，生成转动区域的网格，并输出网格节点号和坐标值，求解器利用用户子程序功能，在进行当前时间步迭代之前读入网格节点和坐标值，实现了当前时间步的网格重构，这种方法克服了传统网格变形和重构时计算时间长和容易产生高畸变网格的缺点。

图3 叶片附近的网格Fig.3 Mesh near vane

3 边界条件和计算工况

实际运行中，滑片泵的进口和出口压力是脉动的，加载时变的压力值更为准确，但这些数据需要通过实验测得，而且也难以将动态压力与相应的阀片的转动位置对应起来，因此计算时滑片泵的进口和出口设为定值压力边界条件。滑片泵的所有壁面给定无滑移、无渗透边界条件，对于湍流问题，近壁区采用壁面函数。转子和定子网格交接面为滑移面，设定Attached 边界条件，保证转子和定子网格的联通。滑片泵的额定工况转速为3 268 r/min，滑片泵的进口绝对压力分别为47 kPa 、65 kPa、85 kPa 和101 kPa，滑片泵出口绝对压力分别为301 kPa、501 kPa和1 001 kPa，温度为140 ℃，动力粘度为0.002 36 Pa·s，密度为925.9 kg/m3. 考虑空化时，设油汽的密度为液态润滑油的1/10 000、粘度为液态润滑油的1/10，饱和蒸气压为100 Pa. 滑片泵叶片头部和定子壁面间隙设为0.05 mm，轴向端面间隙设为0.02 mm.

4 空化模型参数的确定

为了确定空化模型中需要人为设定的单位液体体内空化核的数目，首先对某一在滑片泵专用实验器上进行的实验工况来进行模拟，该工况的温度为140°，进口绝对压力为47 kPa，出口绝对压力为501 kPa，该工况的实验测得的平均流量为38.1 L/min.因空化核的初始半径对空化区域的大小影响较小，所以选取初始空化核半径为10-6m，分别设单位液体体积的空化核个数为109、1010、1011和1012. 表1 分别为不同空化核数时，滑片泵进口和出口的总平均流量。

表1 空化核数对计算流量的影响Tab.1 Influence of nuclei number on calculated flow rate

根据比较，单位液体体积的空化核个数为1010时，流量的计算值和实验值最接近，因此选择该单位液体体积的空化核个数作为空化模型的设置参数。

5 计算结果分析

图4 实验工况下空化体积分数的分布Fig.4 Distribution of cavitation volume fraction on a section of vane pump under the experimental condition

5.1 滑片泵的流动特性

首先对上述实验工况，且单位液体体积的空化核个数为1010的计算结果进行做进一步的分析。图4 为不同转角时，滑片泵某一截面上空化体积分数的分布，体积分数为1 表示空化后的润滑油蒸气，体积分数为0 表示液态润滑油，0 与1 之间的值为液态润滑和油润滑油蒸气的过渡区。图5 为滑片泵进口和出口瞬时流量与转子转角的关系曲线，进口瞬时流量随时间变化较小，而出口瞬时流量随时间变化较大。靠近吸油腔侧的进口流道内存在空化区，该区域的空化区域大小相对稳定，使得进口流量随转角变化不明显;在吸油腔内，叶片的背风侧存在空化区，这些空化随着叶片的转动而不断产生、发展和溃灭，并引起这些区域内质量的堆积(空化区减小)和释放(空化区增大)，从而引起出口流量的脉动。当叶片转动并打开封油区时，压油腔内的压力高于封油区内的压力，润滑油倒流回封油区内，此时，叶片背风区内的压力升高，空化区减小，出口流量为负值，呈现吸油状态。在周期时间内进口的平均流量为36.53 L/min，出口流量为37.66 L/min，二者的平均流量为37.10 L/min，进口和出口的平均流量基本一致，但由于引入空化模型后，数值上处理空化区域和液态润滑油区域的界面时会存在一定的数值误差，导致进、出口的平均流量有一定差别，但整体上计算是符合质量守恒定律的，计算值与实验值非常接近。

图5 滑油泵进、出口瞬时流量与转子转角的关系曲线Fig.5 Relation between inlet/outlet transient flow rates and rotor angle

图6 转动区域内空化体积分数的分布Fig.6 Distribution of cavitation volume fraction in rotating zone

图6为某时刻的叶片头部与定子壁面间隙和叶片轴向间隙的空化体积分数分布，与图4 相比可知，叶片间隙引起的空化区域的大小和强度比叶片表面空化小。叶片端面间隙引起较大的泄露速度，使得对叶片端部背风侧形成较大的空化区。

5.2 不同出口压力对滑片泵流动特性的影响

为了研究滑片泵不同出口压力对滑片泵流动特性的影响，在出口绝对压力为501 kPa 工况的基础上，增加了出口绝对压力为301 kPa 和1 001 kPa 的工况。图7 为不同出口压力下，滑片泵出口瞬时流量与转子转角的关系曲线。由图7 可知，出口压力增加，封油区打开后，从压油腔倒流回封油区的润滑油增加，使得出油口的回流流量增加，但与回流区相对应的转角范围减小。与出口压力301 kPa、501 kPa和1 001 kPa 相对应的平均流量分别为37.70 L/min、37.10 L/min 和36.55 L/min，可见出口压力增加，滑片泵的平均流量减小，但减小幅度不大。图8 为出口压力1001 kPa 的工况，封油区打开前后的空化体积分数的分布，与图4 相比，出口压力提高，叶片背风区内空化溃灭加速，出口瞬时流量与转子转角的关系曲线变得更陡，出口瞬时流量波动增大。

图7 出口压力对滑片泵出口瞬时流量的影响曲线Fig.7 Curves of outlet pressures influencing the outlet transient flow rate

图8 空化体积分数的分布Fig.8 Distribution of cavitation volume fraction

5.3 不同进口压力对滑片泵流动特性的影响

为了研究滑片泵不同进口压力对流动特性的影响，在进口绝对压力为47 kPa 工况的基础上，增加了进口绝对压力为65 kPa、85 kPa、101 kPa 的工况。图9 为不同进口压力下，滑片泵出口瞬时流量与转子转角的关系曲线。由图9 可知，进口压力增加，封油区打开，从压油腔倒流回封油区的润滑油减少，出油口的回流减小，与回流相对应的转角范围减小，进口压力101 kPa 时，出口回流基本消失。与进口压力47 kPa、65 kPa、85 kPa 和101 kPa 的工况相对应的平均流量分别为37.10 min、42.33 L/min、48.07 L/min 和51.89 L/min，可见随着进口压力的提高，滑片泵平均流量增加，且增幅较大。图10 为进口压力101 kPa的工况，封油区打开前后的空化体积分数的分布，与图4 相比，叶片背风区内空化区域减小，但由于出口压力不变，因此叶片背风区内空化溃灭速率变化不大，使得不同进口压力时，出口瞬时流量与转子转角的关系曲线的斜率基本一致，出口瞬时流量波动减小。

图9 进口压力对滑片泵出口瞬时流量的影响曲线Fig.9 Curves of inlet pressures influencing the outlet transient flow rate

图10 空化体积分数的分布Fig.10 Distribution of cavitation volume fraction

6 结论

利用Star-cd 用户子程序功能调用自编程序，实现滑片泵运动区域的网格重构和叶片与定子间隙的网格生成，以单位液体体积的空化核个数1010为空化模型的设置参数，模拟了滑片泵在不同进出口压力状态下的流动特性，得到如下结论:

1)在进、出口设定压力边界条件时，滑片泵工作过程中，由于进口流道内有一相对稳定的空化区，瞬态进口流量波动较小，由于在吸油区和封油区内出现动态变化的空化区，而使得出口流量有较大的波动。

2)滑片泵进口压力增加，叶片背风区内空化区域减小，但空化溃灭速率变化不大，使得出口瞬时流量波动减小;出口压力增加，叶片背风区内空化区域变化不大，但空化溃灭加速，使得出口瞬时流量波动增大。

3)滑片泵进口压力增加，出油口的回流减小，与回流相对应的转角范围减小，使得滑片泵的平均流量呈较大幅度的增加;滑油泵出口压力增加，出油口的回流流量增加，与回流区相对应的转角范围减小，使得滑片泵的平均流量呈小幅减小。

References)

［1］Brennen C E. Cavitation and bubble dynamics［M］. New York:Cambridge University Press，2014.

［2］吴玉林. 流体机械及工程［M］. 北京:中国环境科学出版社，2005.WU Yu-lin. Fluid machinery engineering［M］. Beijing:China Environmental Science Press，2005. (in Chinese)

［3］侯训波，孙玉清，宋昌平，等.叶片泵定子曲线及叶片对流量脉动的影响［J］. 液压与气动，2008，11(8):57 -59.HOU Xun-bo，SUN Yu-qing，SONG Chang-ping，et al. Analysis on cam ring and vane arousing flow fluctuation for hydraulic vane pump［J］. Chinese Hydraulics ＆ Pneumatics，2008，11(8):57 -59. (in Chinese)

［4］张海竹，卢勇，张薇，等.高压叶片泵流体泄漏研究［J］.流体机械，2009，37(3):1 -5.ZHANG Hai-zhu，LU Yong，ZHANG Wei，et al. Analysis on the flow leakage in high pressure vane pump［J］. Fluid Machinery，2009，37(3):1 -5.(in Chinese)

［5］王力，权龙. 单作用叶片泵瞬时流量的分析与计算［J］. 振动、测试与诊断，2006，26(3):188 -191.WANG Li，QUAN Long. Transient flux analysis and calculation on single-acting vane pump［J］. Journal of Vibration，Measurement＆ Diagnosis，2006，26(3):188 -191. (in Chinese)

［6］李少年，魏列江，王峥嵘，等.油液中气泡含量对高压叶片泵工作腔油液压力和叶片受力的影响［J］. 液压与气动，2010，13(9):91 -94.LI Shao-nian，WEI Lie-jiang，WANG Zheng-rong，et al. Influence of the bubble content in oil on working cavity pressure and force of vane in high pressure intra-vane type pump［J］，Chinese Hydraulics ＆ Pneumatics，2010，13(9):91 -94.(in Chinese)

［7］那焱青，王峥嵘，李少年，等.双作用子母叶片泵瞬时流量的分析［J］.兰州理工大学学报，2004，30(6):58 -60.NA Yan-qing，WANG Zheng-rong，LI Shao-nian，et al. Investigation of transient flow in double-acting pump with composite vanes［J］，Journal of Lanzhou University of Technology，2004，30(6):58 -60.(in Chinese)

［8］Osama A H. Theoretical modeling of sliding vane compressor with leakage［J］. International Journal of Refrigeration，2009，32(2):1555 -1562.

［9］Giuffrida A，Lanzafame R. Cam shape and theoretical flow rate in balanced vane pumps［J］. Mechanism and Machine Theory，2005，40(3):353 -369.

［10］Lu Y，Zhang W，Zhao Y Y，et al. Studies on several key problems of water hydraulic vane pump［J］. Industrial Lubrication and Tribology，2011，63(3):134 -141.

［11］Danardono D，Kim K S，Roziboyev E，et al. Design and optimization of an LPG roller vane pump for suppressing cavitation［J］.International Journal of Automobile Technology，2010，11(3):323 -330.

［12］Wang D M，Ding H，Jiang Y，et al. Numerical modeling of vane oil pump with variable displacement，2012-01-0637［R］. US:SAE，2012.

［13］Erzincanli B，Sahin M. An arbitary Lagrangian-Eulerian formation for solving moving boundary problems with large displacements and rotations［J］. Journal of Computational Physics，2013，255(15):660 -679.

［14］Liu J X，Akay H，Ecer A，et al. Flows around moving bodies using a dynamic unstructured overset-grid method［J］. International Journal of Computational Fluid Dynamics，2010，24(6):187-200.

［15］ildirim B，Lin S，Mathur S，et al. A parallel implementation of fluid-solid interaction solver using an immersed boundary method［J］. Computers ＆ Fluids，2013，86(5):251 -274.

［16］Hirt C W，Nichols B D. Volume of fluid(VOF)method for the dynamics of free boundary［J］. Journal of Computational Physics，1981，39:201 -225.

［17］Ubbink O，Issa R I. A method for capturing sharp fluid interfaces on arbitrary meshes［J］. Journal of Computational Physics，1999，153:26 -50.

［18］Jofre L，Lehmkuhl O，Castro J，et al. A 3-D volume-of-fluid advection method based on cell-vertex velocities for unstructured meshes［J］. Computers ＆ Fluids，2014，94(1):14 -29.

［19］Bakir F，Rey R，Gerber A G. Numerical and experimental investigations of the cavitating behavior of an inducer［J］. International Journal of Rotating Machinery，2004，10(1):15 -25.

［20］Ahuja V，Hosangadi A，Arunajatesan S. Simulation of cavitating flows using hybrid unstructured meshes［J］. Journal of Fluids Engineering，2001，123(2):331 -340.

［21］Fard M B，Nikseresht A H. Numerical simulation of unsteady 3D cavitating flows over axisymmetric cavitators［J］. Scientia Iranica B，2012，19(5):1258 -1264.

［22］Menter F R. Zonal two-equation k-w turbulence model for aerodynamic flows，1993-2906［R］. US:AIAA，1993.

［23］Langtry R B Menter F R. Transition modeling for general CFD application in aeronautics，2005-522［R］. US:AIAA，2005.

［24］Demirdzic I，Peric M. Space conservation law in finite volume calculations of fluid flow［J］. International Journal for Numerical Methods in Fluids，1988，8:1037 -1050.

［25］Wilcox D C. Turbulence modeling for CFD ［M］. California:DCW Industries Inc，2006.

［26］Van L B. Towards the ultimate conservative difference scheme V.a second order sequel to Godunov's method［J］. Journal of Computational Physics，1979，32:101 -136.

［27］Buffard T，Clain S. Monoslope and multislope MUSCL methods for unstructured meshes［J］. Journal of Computational Physics，2010，229(10):3745 -3776.

［28］Oliveira P J，Issa R I. An improved PISO algorithm for the computation of buoyancy-driven flows［J］. Numerical Heat Transfer，2001，40:473 -493.