网格法和粒子法模拟齿轮搅油损失中的比较分析

宁 斌，王 斌，吕红明，王交龙，周雅杰

(1.江苏大学机械工程学院，江苏 镇江 212013)

(2.盐城工学院机械工程学院，江苏 盐城 224051)

(3.盐城工学院汽车学院，江苏 盐城 224051)

网格法是当前处理流体力学问题的常用方法，对于简单的流动问题能给出准确的结果。网格法指的是在仿真前用若干个子区域来代替原本的计算区域，且每个子区域内的物理量都是由该区域内的节点来决定，通过计算各个节点上的代数方程组来解决流体问题。目前，网格法已经成为研究流体领域问题的主要手段。如彭钱磊等[1]利用网格法对齿轮飞溅润滑性能进行数值仿真；王斌等[2]利用网格法对渐开线直齿轮在混合润滑状态下的效率进行了分析；尹雪梅等[3]在分析油膜轴承压力分布时采用了动网格技术；Blocken等[4]通过网格法对大气边界层进行了CFD模拟。

在处理自由面大变形问题时，网格法很难给出准确的结果。20世纪90年代，MPS法(moving particle semi-implicit method)作为一种无网格法开始受到人们关注[5]。MPS法以Lagrange方程为基础，用具有密度、速度、动量等物理量的流体粒子代替所要计算的流体域，通过研究流体粒子的情况来观察整个流域的变化情况。这种方法在解决自由表面发生破碎与翻卷等现象时有很大的优势。MPS法因为其显著的优点受到人们的关注，并在许多流体问题中得到应用，如：朱跃等[6]分析了MPS法在多相多组分界面流动的应用；Khayyer等[7]对提高MPS法的稳定性及准确性进行了研究；孙中国等[8]分析了MPS法在流体机械数值模拟中的应用；梁杨杨等[9]采用MPS法对液滴斜向撞击液面过程进行了数值模拟。

网格法和粒子法都可以用来模拟流体动力学问题，但网格法与MPS法在许多方面存在差异，如前期处理、算法、计算效率等。本文以齿轮搅油模型为例，分别采用网格法和MPS法进行齿轮搅油仿真并对比分析仿真结果与实验结果，为以后选用分析流体动力学问题的方法提供了参考。

1 数值方法

1.1 控制方程

1.1.1网格法控制方程

网格法中，液体流动问题的控制方程[10]为

(1)

式中：ρ为质量密度；φ为通用物理量；t为空间某时刻；u为质点速度；Γ为速度环量，定义为速度矢量沿封闭曲线的线积分；S为空间内质点或质点域的内能与势能总和；div为散度函数，反映物理量在各矢量方向上的分散程度；grad为梯度函数。此式为动量守恒方程以及能量守恒方程的通用公式。

1.1.2MPS法控制方程

在MPS法中，基本控制方程主要为Continum方程(质量守恒方程)和Navier-Stokes方程(动量守恒方程)，具体方程为：

(2)

(3)

在MPS法中，可以通过对粒子运动情况的分析，获得整个流体域的运动状况。在运动过程中，各个粒子都会对相邻粒子产生作用，对于粒子间相互作用，主要采用权函数[11]进行计算。

(4)

式中：W(r)表示粒子间的相互作用；r为粒子间的距离；re为粒子作用半径。

从式(4)可知：当r

1.2 算法的选用比较

1.2.1SIMPLE算法

本文选用网格计算方法中最为普遍的SIMPLE[12]算法。它是一种半隐式算法，通过压力修正理论对方程进行迭代求解，并使用上一次的迭代结果构建下一次的迭代，主要计算流程为：

1)设定压力值P*(P*为任意值)。

2)将P*与动量方程相结合，求解速度场，u*,v*,w*表示速度分量。

3)由于P*为空间某时刻的任意值，因此一般情况下通过动量方程获得的速度矢量u*也为任意值，即与周围空间点存在速度梯度误差，所以此时获得的速度场不满足连续方程，需要构建压力修正值P′，P′和P*相加后就可得到准确的压力值P,P=P′+P*。

4)通过P可以获得更加准确的速度场u,v,w(符合连续方程)。

5)将P设为新预测值P*，重新开始迭代计算。

1.2.2时间积分法

MPS法通过半隐式的时间积分法来解决问题，该算法与SIMPLE算法一样，主要用来解决不可压缩流体的问题。时间积分法求解过程包含两步：显式修正和隐式修正。具体计算流程如下：

(5)

(6)

2)通过计算Poisson方程，得到下一时刻的压力；

(7)

3)由所获得的压力，获得该时刻的粒子速度及位置：

(8)

(9)

1.3 边界处理方式比较

1.3.1边界层理论

网格法中，流体的控制方程为非线性方程，这使得边界问题很难处理。本文以齿轮搅油为模型，在搅油过程中，油液在箱体壁面以及齿轮表面产生湍流。湍流区可以分为与箱体壁相互作用的近壁区和远离箱体壁的完全湍流区。由于受到箱体壁的影响，油液在近壁区的流动与在完全湍流区的流动有较大差异。根据与箱体壁的垂直距离，将近壁区划分为3个子层，依次为黏性底层、过渡层以及对数率层，各子层高度和速度如下：

(10)

(11)

1.3.2MPS法边界处理方法

在MPS法中，流体粒子在边界处易发生粒子缺失以及粒子穿透的现象，这样将无法保证仿真结果的准确性。本文通过在边界处以及边界外布置多层边界粒子[13]来防止这类现象的发生。由图1可知，边界粒子分为2类：第一类为与流体接触的壁面粒子，这一类粒子主要分布在固壁上的壁面粒子；第二类为分布在边界外的壁外粒子。在计算过程中，采用权函数(式(4))计算壁面粒子与流体粒子。若计算时不考虑壁外粒子，一方面会使得壁面粒子被认定为流体粒子，另一方面会使计算所得的粒子密度小于实际值，所以需要在边界外分布多层粒子，壁外粒子的层数由权函数中re的大小决定。

图1 壁面边界的粒子分布

在求解过程中，为了防止粒子在与壁面接触后穿过边界，需要计算壁面粒子的压力，这样能够确保液体粒子撞到壁面后能够回弹。

2 仿真过程与结果比较分析

本文分别采用网格法与MPS法进行齿轮搅油仿真，并对仿真结果进行比较分析。首先齿轮搅油模型是一个流体领域的典型案例，其次齿轮搅油损失是测量电动汽车减速箱功耗的重要参考指标，具有实际意义。不仅如此，在齿轮啮合转动过程中，油液面会发生较大的变形，尤其是当油液撞击到箱体壁后，液面会发生翻卷、破碎等自由面流动现象。因此分别用网格法和MPS法对齿轮搅油模型进行数值仿真并比较获得的仿真结果,具有一定的参考价值与研究意义。

2.1 建立搅油模型

建立NX搅油齿轮箱模型，如图2所示。一对啮合齿轮安装在箱体中间，将油液倒入箱体内，油液宽度L=360 mm，高度H=145 mm，油液面与齿轮中心在同一水平面，箱体容积为360 mm×290 mm×100 mm。表1为齿轮副基本参数，其中主动轮为G1、从动轮为G2。

表1 齿轮副基本参数

图2 搅油模型几何尺寸

2.2 前期处理比较

2.2.1网格法前期处理

在使用网格法进行仿真前要对模型进行网格划分，由于齿槽间的间隙过小，使得在啮合区域所划分的网格质量较差，会对仿真结果产生影响，因此本文将两个齿轮分别向两侧移动1 mm，增加啮合区域间隙,从而确保网格的质量。

在求解过程中，模型的网格数越少，则求解速度越快。本文通过局部加密法来划分网格，对轮齿周围区域进行加密，并在齿槽区域划分双层网格，这样能够减少求解时间、提高仿真准确性，局部网格加密图如图3所示[14]。

图3 局部网格加密图

随后将划分好的网格模型导入仿真软件内，采用压力基求解，如图4所示。由于齿轮在转动时，流体会产生剧烈飞溅，油液面无法确定，因此由流体单元网格与总体积比函数F来确定自由面。

图4 VOF两相流设置

2.2.2MPS法前期处理

分别定义箱体、油液和齿轮3个部分的物理属性，如图5所示，将油液的物理属性定义为Fluid(流体)，齿轮组与箱体的物理属性定义为Polygon(多面体)。Fluid为MPS法仿真软件中最为重要的元素，由粒子充当流体的角色，可以定义流体密度、运动黏度等物理量。而Polygon为盛放液体的容器，以自身轮廓为边界条件，表达液体和粒子的接触。

图5 定义各组件物理属性

2.3 仿真结果比较

本文分别采用ANSYS[14]和Particleworks两种软件对齿轮搅油进行模拟仿真。如图6所示，初始时，齿轮副未开始转动，油液处于静平衡状态，在t=0时，主动轮G1开始转动，通过齿轮副的转动带动油液的运动。在计算前，将压力分布设置为静压，油液密度ρ0=900 kg/m3,重力加速度g=9.8 m/s2,运动黏度系数υ=1.82×10-4m2/s。

图6 初始油液位置

图7与图8显示的是齿轮副刚开始转动和油液开始翻卷的典型时刻的仿真图像。通过油液的分布状况可以看出，尽管流场变化比较剧烈，液面发生了翻卷，但网格法与MPS法均能很好地捕捉到这些特征。

图7 齿轮转动0.01 s时流域仿真图

图8 齿轮转动0.1 s时流域仿真图

从细节上看，网格法与MPS法存在着一些差别，如图7所示，在0.01 s仿真刚开始时，在网格法中自由面变化并不是很剧烈，但在粒子法仿真中能够清楚地看出一些油液粒子运动较为剧烈，这种现象在齿轮转动0.1 s后更为明显(图8)。这主要是由于MPS采用了式(4)的权函数，当相邻的粒子在靠近时，粒子数的密度快速增加，使得粒子飞溅。

2.4 数值结果比较

为了比较网格法与MPS法对于齿轮搅油的准确性，用上述搅油模型进行实验，并将实验数据与数值仿真获得的数据进行比较，图9、图10为在不同转速与转向情况下，齿轮副搅油损失变化曲线和数据点。由图可知，网格法与粒子法的计算结果逐渐收敛于相近的结果。由于实验测量时存在着装配误差、传感器测量误差以及其他一些不可控因素，而仿真数据是理想状态下的齿轮搅油功耗，因此两者之间存在一些偏差，但误差大小仍在允许范围内。

图10 搅油功率损失(主动轮顺时针转动)

图9 搅油功率损失(主动轮逆时针转动)

2.5 计算效率及精度的比较

因为网格法和MPS法采用不同的算法，所以在计算效率方面同样存在一定的差异。在网格法中，大量的时间会花费在网格划分中，且网格的质量在很大程度上影响计算的效率和精度。在MPS法中粒子的大小也会影响计算的效率和精度。为了能够准确地比较两类方法的计算效率，本文通过对网格的大小与粒子的大小进行控制，使网格的体积与粒子的体积接近，由此对比在网格体积与例子体积接近时两类方法的计算效率。本文中最大网格体体积约为7.78×10-8m3，最小网格体体积约为3.93×10-12m3。因此本文设粒子的直径为1 mm，则每个粒子体积约为4.19×10-9m3，处于最大网格体积与最小网格体积之间。计算时使用相同的时间步长1×10-4s。

经比较，MPS法的计算时间约为网格法的三分之一。由此可得，MPS法在计算效率上有较大优势。

在计算精度方面，由于网格法软件的限制，本文将齿轮啮合间隙有意扩大了2 mm，加大了齿轮啮合时的空隙和啮合齿轮之间的泵油量，使得齿轮间泵油损失增加，在一定程度上降低了计算的准确性。由此可知，MPS法在计算精度上具有一定的优势。

3 结论

本文通过网格法与MPS法分别对齿轮搅油进行了数值模拟和比较分析，结论如下：

1)在仿真结果方面，网格法和MPS法均能捕捉到齿轮转动时的自由液面的变化状况，说明两种方法都能准确地解决流体动力学方面的问题。

2)在计算搅油损失结果方面，两种方法所计算的结果与实验结果相吻合。

3)在计算效率方面，本文通过控制网格和粒子的体积大小来比较两者的计算效率，使得MPS法的计算效率明显高于网格法。

4)在计算精度方面，由于网格法在进行网格划分前需扩大齿轮间隙，降低了困油区的密闭性，使得搅油时的泵油量高于实际值；而MPS法则完全模拟真实搅油状态，相对计算精度较高。