基于MPS方法的二级齿轮箱飞溅润滑特性研究

谢迟新 刘桓龙 贾瑞河 黎 强

1.先进驱动节能技术教育部工程研究中心，成都，6100312.西南交通大学机械工程学院，成都，610031 3.江苏新瑞齿轮系统有限公司，常州，213000

关键字：飞溅润滑；移动粒子半隐式法；流场特性；润滑特性；搅油功损

0 引言

随着城市轨道交通的快速发展，列车运转密度和速度不断提高，保障列车运营安全变得越来越重要[1]。传动齿轮箱是城市轨道车辆传动系统的重要组成部分，在轨道车的正常高效运行中发挥着重要作用；同时，齿轮箱内润滑油的润滑作用又对传动齿轮箱的正常工作起着至关重要的作用。在工作过程中，齿轮齿面之间的相互摩擦、齿轮与箱内油液及气流、轴承的摩擦等均会产生大量热量，导致润滑油油温上升。油温上升会导致润滑油黏度降低，齿轮附着油膜厚度减小，从而对齿轮和轴承的润滑造成很大影响[2]。润滑过程中，润滑油油位过低会导致润滑不足，造成齿轮和轴承等零件的快速磨损，而油位过高会造成过高的搅油功率损失，损失能量且产生额外热量，因此对齿轮箱飞溅润滑的润滑效果研究具有十分重要的意义。计算流体动力学(computational fluid dynamics, CFD)方法可以十分有效地分析齿轮箱内部的流场特性[3]。近年来，随着计算机运算能力的增强，国内外学者对齿轮箱润滑进行了相关数值仿真分析。姜义尧等[4]基于CFD方法建立了某直升机中减速器飞溅润滑数值仿真模型，对中减速器齿轮啮合处润滑油流量与齿轮转速、浸油深度的关系进行了研究。刘中令等[5]建立齿轮箱的二维VOF多相流模型，获得了不同速度工况下齿轮箱内部流场的相态、压力、速度分布图。陈黎卿等[6]运用高速摄影技术对透明分动器内部流场进行观察实验，得到CFD液面波动情况的数值模拟结果与实验结果的最大误差为12.2%。常帅[7]运用CFD方法对行星齿轮传动飞溅润滑的搅油功率损失进行研究，通过分析功率损耗来源建立了行星传动效率的数学模型。彭钱磊等[8]利用齿面移动法建立汽车驱动桥飞溅润滑模型，通过实验验证了驱动桥搅油功耗的仿真误差在8.73%以内。孙凯等[9]运用fluent动网格技术对中减速器螺旋锥齿轮三维模型进行了动态数值仿真分析。HU等[10]对直升机锥齿轮传动齿轮箱的搅油功率损失进行了研究，分析了齿轮转速、润滑油油位、温度、动力黏度、密度和直升机倾角对搅拌功率损失的影响，并通过实验对功率损失仿真结果进行验证，取得了良好一致性。LIU等[11-12]分别用有限体积法(finite volume method, FVM)和光滑粒子流体动力学(smooth particle hydrodynamics, SPH)方法建立了FZG齿轮试验机模型，对比了实验与仿真的搅油功率损失，并用高速摄像机对油液飞溅状态进行记录，结果表明FVM法对齿轮箱飞溅润滑的油液分布和搅油损失具有良好的预测性，SPH法在油液分布上与试验吻合良好，搅油损失与试验结果存在较大误差。JI等[13]用SPH法和粒子图像测速(particle image velocimetry, PIV)技术定性分析了由齿轮旋转而产生的气泡数量和大小，并对油面下的速度场和速度分布云图进行了综合分析。GROENENBOOM等[14]运用SPH与FVM耦合的方法SPH-FEM对涌浪、变速箱润滑和心血管流量进行研究，结果表明，SPH-FEM方法非常适合湍流和流固耦合模型。皮彪等[15]首次将移动粒子半隐式(moving particle semi-implicit, MPS)方法运用到一级传动减速器的润滑中，得到了润滑油液飞溅润滑的状态和速度分布场，为润滑系统分析提供了一种新方法。LI等[16]将MPS法应用到单齿轮搅油损失研究中，通过公式拟合了不同因素对搅油功率损失的影响，为齿轮搅油损失研究提供了新方法。

文献[4-11]主要基于FVM法在欧拉坐标系的基础上对二维和单级传动齿轮箱模型的飞溅润滑流场特性和搅油损失进行研究分析，结果表明该方法在齿轮箱飞溅润滑问题上具有较高的适用性和准确性。但是该方法存在三处弊端：一是在齿轮啮合处建立有限元模型的复杂性较高；二是对计算机硬件要求很高；三是齿轮啮合间隙非常小，容易造成网格扭曲畸变，导致计算无法进行。因此，FVM法在多级齿轮传动和复杂模型的研究上不占优势。文献[12-14]运用粒子法对齿轮箱飞溅润滑进行研究，这是一种基于拉格朗日坐标系的计算方法，对流体和固体物质本身进行离散，而不是对空间进行离散，它摆脱了网格的束缚，可以避免齿轮啮合处有限元建模复杂和间隙太小容易扭曲畸变的弊端，更适合复杂的空间几何运动。SPH法对流体压力采用显式求解，计算效率更高，但其准确性和稳定性低于MPS方法[17-18]。

移动粒子半隐式法(MPS)通过求解压力泊松方程得到流体的压力场，并通过压力梯度修正预测流体速度，计算更加准确。目前仅有皮彪等[15]、LI等[16]将MPS方法运用于一级传动齿轮箱润滑的研究中。本文以轨道车二级传动齿轮箱为研究对象，运用MPS法对齿轮箱内部流场进行数值仿真，分析了不同输入轴转速、初始润滑油油量和环境温度共12种工况下齿轮箱内部的润滑油流场特性；分析了齿轮关键啮合点的油液粒子数时域变化情况以及齿轮箱的搅油功率损失情况。针对FVM法和SPH法的缺点，本文利用MPS法较高的计算效率和准确性，最大程度地保留了几何模型的完整性，实现了复杂几何模型齿轮箱工况的多样化，同时保证了计算结果的准确性和有效性。

1 MPS法介绍

MPS法是由KOSHIZUKA等于1995年最先提出的一种完全的拉格朗日方法。在该方法中，流体通过有限数量的粒子表达，通过粒子将描述牛顿流体的基本力学Navier-Stokes方程离散化。该方法利用核函数描述粒子之间相互作用的强弱。为了实现控制方程的离散，建立了表达物理量的光滑模型、压力梯度模型和扩散模型[19]。MPS方法通过隐式方法求解压力方程，保证了流体的不可压缩性，并通过引入可接受较大时间步长的半隐式求解方式显著缩短了求解时间。由于MPS方法是一种无网格法，不受网格扭曲畸变的影响，因此在大变形水力学研究中得到了广泛应用。

1.1 控制方程

假设流体是不可压缩且具有黏性的，MPS方法的基础控制方程为拉格朗日形式的质量守恒方程和动量守恒方程[19]。

1.2 核函数

流体是由众多离散的粒子组成的，粒子的运动状态可以表征流体的运动状态，每个粒子包含一组属性值，如位置、质量、动量及能量等。在MPS方法中，粒子间的相互作用采用两粒子间的距离加权函数的方法进行计算，该核函数具有紧支集特性。

1.3 梯度模型和拉普拉斯模型

粒子间的物理量根据距离进行加权平均计算。梯度模型是粒子i和其作用域内所有邻域粒子j的梯度向量的加权平均值[19]。基于非正常扩散，利用拉普拉斯模型将粒子i的物理量根据邻域粒子的距离进行分配。

1.4 时间步长标准

仿真过程中，时间步长的选择十分重要，时间步长过大容易造成计算发散，太小则会增加不必要的计算时间。时间步长Δt的确定需满足以下条件：

(1)

式中，Δti为初始时间步长；C为克朗数；l0为粒子直径；umax为粒子最大速度；di为扩散系数；υ为流体的动力黏度，υmax为其最大值。

1.5 粒子数密度与自由表面

MPS法中，通过保证粒子的密度数恒定来确保流体的不可压缩性。粒子数密度ni指粒子i在核函数作用范围内该粒子和其相邻粒子的函数值的叠加[19]。MPS法将自由面粒子的压力设为0并参与压力泊松方程的迭代，因此，自由表面的判断对计算的准确性很重要。MPS法用下式来判别自由表面粒子[19]：

(2)

在仿真时，粒子处于自由表面，粒子周围的邻居较少，粒子数密度比初始粒子数密度小，所以β应取一个小于1的数。在MPS法理论中，β的建议取值范围为0.8≤β<1[20]。β值偏小时，容易将非自由表面粒子误判为自由表面粒子，引起压力振荡现象，给模拟的结果带来误差，这种情况在流体运动比较剧烈时尤为明显。本文的数值仿真流体运动剧烈，故β取0.97。

2 数值仿真

2.1 几何处理与MPS粒子生成

本文研究的模型为单速比两级平行轴斜齿轮减速器，采用飞溅润滑方式。齿轮箱润滑系统零部件众多，箱体结构复杂，而飞溅润滑主要关注箱体内表面、齿轮、挡油板等关键零部件对润滑油分布造成的影响，为了节省计算时间，对齿轮箱结构进行以下简化：

(1)去除对计算结果影响很小的零部件，如螺栓、螺母、螺钉等。

(2)由于油液运动时只有箱体内表面几何特征会对流体的分布产生影响，故去除箱体外表面无关特征。

(3)MPS法在处理齿轮对接触时容易造成发散，用齿面移动法对齿轮的齿面切除1%厚度可以有效保证仿真的顺利进行，即在保留所有轮齿和不改变齿轮安装位置的基础上，通过改变轮齿厚度来增大啮合区域间隙，以保证计算的正常进行[21]。目前，针对齿轮啮合处间隙狭小、给流体域网格划分和求解带来困难的问题，常用的齿轮建模处理方法包括无齿法[22]、单向切齿法[23]、双向切齿法[24]和分离法[25]。根据彭钱磊等[21]的研究，相对以上的齿轮处理方法，齿面移动法无需改变齿轮副的安装位置，且获得的润滑油飞溅效果符合实际情况，搅油功率与试验结果基本一致，最大误差在8%以内。齿面移动方式如图1所示，简化后的几何模型如图2所示。

图1 齿面移动法

图2 简化后几何模型

仿真前，运用MPS法对4种液位高度的液体进行粒子化处理。根据计算机性能，将粒子直径设置为2 mm。粒子化处理后，各液位的粒子数量见表1，其中，V表示润滑油体积，n表示粒子个数。任意时刻，齿轮箱内任意位置润滑油油液的体积V与粒子数目n的关系为

表1 不同油量下的粒子数量

V=kn

(3)

式中，k为常数，与粒子直径有关。

MPS油液粒子模型如图3所示。粒子化过程中，为避免粒子化失败问题，在几何模型处理过程中，经过仔细排查与处理，确保了齿轮运转时运动零部件的几何特征不发生干涉，同时保证了齿轮箱的几何封闭性。另外，通过CAD软件的几何处理功能，保证了不同油量工况下液位高度的准确性。

图3 油液粒子模型

2.2 参数设置

利用Particleworks软件对模型进行单相数值仿真。Courant number设置为0.2，时间步长为1.2×10-4s，参数设置满足CFL条件。压力与黏度模型采取隐式算法，液体表面张力采用potential模型，并采用双精度进行求解。考虑粒子受到的重力，重力加速度为9.8 m2/s。由于齿轮箱内润滑油流动情况复杂，采用LES湍流模型求解。

采用某品牌润滑油对齿轮箱进行润滑，不同温度下润滑油的属性见表2，表中，θ为润滑油温度，ρ为润滑油密度，γ为运动黏度。对不同转速、不同油量、不同温度共12种工况进行仿真分析，每种工况的求解物理时长设置为2 s。齿轮转速载荷的赋予采用斜坡方式，0～0.1 s为加速阶段，0.1～2 s为匀速阶段。基于产品具体参数保密的原因，不同工况下的转速和油量具体参数值隐去，由低至高4种转速分别用R1、R2、R3和R4代替，由少至多4种油量分别用V1、V2、V3和V4代替，不同工况参数见表3，表中，nI为输入轴转速。

表2 不同温度下润滑油的属性

表3 12种工况的初始油量和输入轴转速

3 齿轮箱流场特性分析

通过后处理，得到不同工况下粒子的运动形态和物理特性，以便更好地分析齿轮箱润滑油飞溅润滑的流场特性。本节通过齿轮箱内部粒子的速度场、压力场和迹线分布的可视化云图对齿轮箱流场特性进行定性分析。齿轮箱内部相关结构如图4所示。

图4 齿轮箱体内部结构

3.1 速度场分析

图5为不同工况下粒子的速度场分布云图。从工况1～4可看到，随着输入轴转速的增大，齿轮箱内油液粒子的速度增大，更多的滑油被齿轮搅起达到箱体顶部，在挡油板作用下落入集油槽，对传动系统润滑起到积极作用。从工况5～8可看到，随着滑油初始油量的增加，齿轮的搅油量增大，齿轮箱的传动系统得到更大油量，润滑效果更好。从工况9～12可看到，相对于环境温度40 ℃，当温度为-40 ℃和-25 ℃时，滑油的运动黏度大幅增大，流动性降低，导致滑油大股聚集在一起，其质量增大，粒子速度降低，同时滑油粒子分布的均匀性减弱，集油槽油量减少，润滑效果较差。而当环境温度为60 ℃时，相对于40 ℃，滑油的运动黏度减小约一半，黏附在轮齿上的作用力减小，齿轮的搅油量减少，润滑效果相对差。综上，在4种环境温度工况下，温度为40 ℃时润滑效果最好。

图5 不同工况下速度场分布云图

3.2 压力场分析

图6显示了轮系的旋转方向。图7为不同工况下粒子的压力场分布云图。从图7中可以看到，粒子压力值最大的位置为区域Ⅰ，主要由于该区域粒子受到左侧回流油液粒子的冲击，回流由未到达箱体顶部的油液在重力作用下沿着箱体壁面流下形成。粒子压力值相对较小的位置为区域Ⅱ，即输出轴齿轮正下方的齿轮箱底部区域，主要由于齿轮箱内油液粒子的振荡冲击力和所受到的压力。其他粒子压力值较大的区域主要位于油液粒子集中处，如集油槽、导油板。

图7 不同工况下压力场分布云图

3.3 迹线分析

图8为不同工况下粒子的迹线分布云图。从工况1～4可看到，随着输入轴转速的增大，迹线的速度逐渐增大，迹线的分布区域也逐渐增大。当转速为R4时，齿轮箱内迹线的分布区域最广，整体的平均速度最大。从工况5～8可看到，随着滑油初始油量的增加，齿轮箱内迹线的分布区域逐渐增大，尤其当初始油量为V3和V4时，齿轮箱内迹线的密度明显增大，可认为润滑效果得到了明显增强。从工况9～12可看到，当温度为-40 ℃和-25 ℃时，迹线的分布主要集中在各齿轮的齿面上，齿轮箱顶部和集油槽的迹线分布很少，可能存在润滑不良风险。当温度为40 ℃和60 ℃时，可看到温度为40 ℃时齿轮箱迹线的分布区域和分布密度更大，润滑效果更好。

图8 不同工况下迹线分布云图

4 啮合点润滑情况分析

为进一步定量分析齿轮传动系统的润滑情况，选取图9中啮合处Ⅱ，提取各工况下仿真结果数据，得到啮合处Ⅱ粒子数量时域变化曲线，如图10～图12所示。因为在工程运行情况下，啮合处Ⅱ更容易发生润滑不良风险，故本节仅对啮合处Ⅱ的润滑情况进行分析。

图9 齿轮啮合区域

图10 不同转速工况下啮合处粒子数变化曲线

图11 不同油量工况啮合处粒子数变化曲线

图12 不同温度工况下啮合处粒子数变化曲线

从图10可以发现，环境温度为40 ℃、初始滑油量为V3时，在4种不同转速工况下，启动阶段，啮合处均无粒子经过，随着时间推移，粒子数振荡上升，最终在一个小幅度范围内波动，逐渐达到稳态，且最终的稳态粒子数平均值相差很小。说明在所给温度和油量条件下，啮合处的稳态滑油量在输入轴转速[R1,R4]范围内基本不变。

从图11可以看到，粒子数最终达到动态平衡。在其他边界条件一定的情况下，随着初始油量的增加，啮合处粒子数也随之增加。

从图12可以看到，环境温度为-40 ℃和-25 ℃时，啮合处更早出现油液粒子，主要由于低温下滑油的黏度增大，油液粒子在齿面的黏附力更大，不容易被齿轮旋转时的离心力甩至空中。温度为40 ℃时，啮合处粒子数约是60 ℃时粒子数的2.5倍，主要由于60 ℃时滑油运动黏度约为40 ℃时滑油运动黏度的一半，在齿面的黏附力更小。说明在4种环境温度下，40 ℃时润滑效果最好，-40 ℃时润滑稳定性较差。

5 齿轮搅油功率损失分析

飞溅润滑齿轮箱在工作时会产生啮合摩擦功率损失、风阻功率损失和搅油功率损失，搅油功率损失约占总功率损失的30%[25]。降低搅油功率损失可有效降低齿轮箱功率损失和发热量，显著提高齿轮箱的传动效率和使用寿命。从仿真数据中提取出各工况的搅油力矩损失时域变化曲线，如图13～图15所示。

图13 不同转速工况下力矩损失变化曲线

图14 不同油量工况下力矩损失变化曲线

图15 不同温度工况下力矩损失变化曲线

从图13可以看到，在前期力矩损失很小，随着时间推移呈现振荡上升，经过一段时间大幅波动，最后在小幅度范围内振荡，达到动态稳定状态。主要由于齿轮加速采用斜坡方式，初始阶段转速低，力矩损失小，然后流场进入剧烈变化阶段，力矩损失迅速增大，经过一段时间的波动流场逐步稳定，力矩损失逐渐趋于稳定。另外可以发现，随着输入轴转速增大，齿轮的搅油力矩损失随之增大，且当转速为N3和N4时，搅油力矩损失大幅提升，力矩损失与转速成非线性关系。

从图14可以看到，力矩损失变化趋势与图13相似，且初始滑油油量越多，力矩损失增大的响应速度越快。同时，滑油量越多，齿轮搅油力矩损失也随之增大。

从图15可以看到，当环境温度分别为-40 ℃和-25 ℃时，齿轮的搅油力矩损失很大，主要由于温度对润滑油黏度的影响很大，-40 ℃时滑油运动黏度约为40 ℃时的滑油运动黏度1000倍，增大了齿轮的运转阻力，说明齿轮箱在-40 ℃的低温下会存在工作困难问题。齿轮箱搅油功率损失公式为

(4)

式中，Ti为单个齿轮的搅油力矩损失；ni为单个齿轮的转速；P为总搅油功率损失。

根据上式得到不同工况下齿轮箱搅油功率损失，如图16所示，可以看到，齿轮箱的搅油功率损失随着输入轴转速和初始滑油油量的增大而增加，随着环境温度的升高而降低，且均为非线性变化。值得注意的是，当输入轴转速为高转速时，相对中低转速，搅油力矩损失会大幅增加。另外，环境温度为-40 ℃时，齿轮箱可能存在工作困难的危险。

图16 不同工况下功率损失对比

6 结论

(1)将移动粒子丰隐式(MPS)法成功运用到二级斜齿轮传动齿轮箱复杂模型的飞溅润滑可视化仿真中，为复杂模型的润滑研究提供了新的研究方法。

(2)提取了不同工况下的速度场、压力场和迹线分布的可视化云图，通过可视化云图对齿轮箱飞溅润滑流场特性进行了定性分析。

(3)提取了不同工况下齿轮啮合点的油液粒子数时域变化情况，结果表明：啮合处滑油粒子数在输入轴转速[R1,R4]范围内基本不变；初始油量越大，啮合处粒子数也随之增大；在-40 ℃、-25 ℃、40 ℃、60 ℃四种环境温度下，40 ℃时啮合处滑油粒子数最多，润滑效果最好，-40 ℃条件下润滑稳定性较差。

(4)提取了各工况齿轮的搅油力矩损失时域变化曲线，计算分析了各工况的齿轮箱搅油功率损失情况。结果表明：搅油功率损失与输入轴转速和初始滑油油量成正相关，与环境温度的提高成负相关，且均为非线性变化。