基于MatLab与C/C++混合编程的弹流润滑数值算法改进*

莫易敏 熊 钊 王 骏2 胡 杰 洪 叶

(1.武汉理工大学机电学院 湖北武汉 430070；2.上汽通用五菱汽车股份有限公司 广西柳州 545007)

弹流润滑的计算与工程实际联系紧密，存在于滚动轴承、滑动轴承、齿轮传动等多个领域[1]，其中主要的雷诺方程是二阶偏微分方程，以往依靠解析方法求解，必须经过许多简化处理才能获得近似解，因而理论计算具有很大的误差[2]。当今由于计算机技术的迅速发展，复杂的润滑问题有可能进行数值解算。目前，润滑计算作为一种理论依据在机械设计与润滑材料的选用分析中占有重要的地位[3]。

本文作者使用混合编程实现求解弹流润滑方程，不考虑压力迭代方法的影响下，探究了膜厚方程中膜厚常量变化量ΔH对迭代相对精度的影响，得出了较为合适的ΔH变化形式，并依据ΔH的大小对相对精度的影响程度，在整体压力收敛精度趋于稳定时，在极小相对精度处扩大ΔH，使相对进度减小，满足收敛要求，提前收敛，从而实现了弹流润滑的快速数值计算。

1 弹流润滑基本方程及其离散化

如图1所示，研究采用的是工程实际中最普遍的两接触体主平面重合椭圆接触弹流润滑模型[1]。

图1 主平面重合的椭圆接触模型

Reynolds方程为(只考虑x方向速度)[3]：

(1)

量纲一化并离散后：

式中：p为压力，不同位置压力不同，可以看作与x、y有关的二元函数；p(x,y)εi±1/2,j=0.5(εi,j+εi±1,j)，ε0=εi-1/2,j+εi+1/2,j+εi,j-1/2+εi,j+1/2。

不考虑温变时，等温情况下，润滑油黏度、密度、膜厚变化因素只需要考虑压力，黏度-压力方程[4]为

η=η0e(lnη0+9.67)×[(1+p/p0)z-1]

(2)

式中:z一般取0.68。

密度-压力方程为

(3)

式中：η0、ρ0分别为p=0时等温条件下润滑剂的黏度与密度。

膜厚方程[5]：

(4)

量纲一化后：

其中:v(x,y)为在油膜压力下接触表面产生的弹性形变方程[1]：

(5)

离散后：

式中：h0为膜厚常量，是定值；E是两接触表面的综合弹性模量；V(i,j)任意点的弹性形变。

载荷方程[6]：

(6)

量纲一化并离散后有：

需要注意的是：在对上述方程量纲一化的过程中，须将初始化的椭圆接触区长轴长作为量纲一化的基数，即量纲一化后参数a为长轴；进行数值计算中只有含有偏微分与积分的项，需要进行离散化[7]。

2 离散方程迭代求解方法的改进

基本方程中离散化后不考虑温度的情况下，需要求解的是2个独立的方程，即Reynolds方程和载荷平衡方程。其中Reynolds方程由于压力的迭代求解，载荷平衡方程由于膜厚常数调整的判定，文中选用的是较为简便的差分法进行压力迭代，其他的压力迭代法如直接迭代法和多重网格积分法等暂不做讨论。

现有的膜厚常数调整方法单一，均为膜厚常数的变化量ΔH固定周期逐步减小(常用的是每10次压力迭代，ΔH值减为原来的1/2)，方法固定，且对不同的参数没有较好的针对性，只是一个经验值。

本文作者将探究在不同的ΔH(膜厚常量变化量)的变化方式下相对精度的具体变化情况，寻找出压力迭代整体的相对精度较低的ΔH变化方式。在探究得到的ΔH与相对精度的规律上，通过改变ΔH的操作，以期能够降低原有的最小相对精度，从而实现算法的提前收敛，提高效率，使算法更具操作性。

3 MatLab与C/C++混合编程实现弹流润滑计算

在数值、科学和工程计算领域，FORTRAN、C/C++语言更多地用于底层函数开发、软件开发、单片机控制等，其运算速度很快，对大型数据的运算效率高，尤其对高性能计算和并行化计算方面，但是对数据进行画图、图像处理分析等方面，实现起来比较复杂[8]，并且FORTRAN的运算不够灵活。而MatLab是第四代的解释性编程语言，编程思维符合人的思维习惯，数值计算方面编程容易实现，并且有强大的数据、图像处理工具箱，但是其对于循环运算的速度很慢[9]。文中结合2种语言的编程优点，在MatLab平台下，将弹流润滑算法中耗时的主要循环部分弹性变形矩阵用C/C++语言实现，并将其编译成能被MatLab调用的动态链接库[10](.dll文件)。算法的具体实现如图2所示。

图2 混合编程弹流润滑数值计算流程

弹流润滑算例的参数取x、y方向上综合曲率半径分别为0.03、0.06 m，载荷w=1 000 N，初始黏度η0=0.02 Pa·s，两接触表面沿x方向平均速度us=3 m/s，综合弹性模量E=221 GPa。计算区域量纲一化X起始与终点坐标分别为-2.5、1.5，Y方向上则是-2、2。

通过MatLab计时工具tic、toc指令，分别记录MatLab与调用C/C++编译成的dll文件实现计算变形矩阵VI每次迭代运行时间，结果如表1所示。

变形矩阵VI具体是计算求解域内65×65等距节点的压力弹性变形，4层循环运算次数达654次，而由运行情况来说，10次迭代，2种方法结果一样，但运行时间相差巨大。因此，若有需要，将MatLab中耗时的改为C/C++程序，能显著减少计算时间。

表1 不同编程方法计算时间对比

4 数值计算算法分析与改进

4.1 算法收敛影响因素分析

弹流润滑压力、膜厚分布求解的可靠性，主要取决于结果是否满足压力迭代收敛条件与载荷平衡方程，而压力迭代算法的选取与膜厚常数最终取值是算法的主要影响因素[7]。

求解弹流润滑方面的问题，需联立求解各类复杂的方程，涉及流体润滑、载荷方程、固体表面弹性变形方程、润滑剂黏度方程及密度方程，系统非线性强，不能得出解析解。目前，对该类问题的计算方法有逆解法、牛顿有限元法、复合直接迭代法、牛顿有限差分法、直接迭代法和多重网格积分法[11]。

文中采用牛顿有限差分法，而膜厚常数方面一般根据Hamrock-Dowson最小膜厚的经验公式：

算出最小膜厚并量纲一化为Hmin，此时可以通过Hmin反求出膜厚常数H0的初值，之后根据载荷方程平衡调整膜厚常数：H0=H0±ΔH，此处一般取ΔH=0.1～0.05Hmin，并且每固定次迭代，ΔH的值就变为原来的1/2[9]。

压力迭代收敛准则采用的是相对精度判断准则：

式中：e一般取0.01～10-6。

图3所示是ΔH=0.1Hmin，每3次迭代减半，迭代150次，取得的最优结果，为了显示清晰，只截取了X方向 -1.5、1.5。可知：压力分布出现了二次应力峰，但是由于材料刚度、尺寸较大，其应力峰的值并不大，因此出口处膜厚减小现象也不明显，但是计算结果是符合弹流润滑的一般特征的。

图3 弹流润滑压力与膜厚分布

图4示出了不同ΔH变化下相对精度随迭代次数变化规律。由不同情况下100次迭代曲线可以看出来，ΔH的不同变化形式是不会显著影响收敛相对精度的总体变化趋势，并且随着ΔH的不断减小，ΔH对压力迭代的影响可以忽略不计，此时主要影响因素归于迭代算法的特性。图4中在迭代次数大于30次后，相对精度在10-3～10-6范围内周期变化，变化周期大致为每22次迭代相对精度并未发散，这一特征符合压力迭代算法收敛的要求[8]。

图4 不同ΔH变化下相对精度随迭代次数变化曲线

4.2 膜厚变化对相对精度影响分析与改进

由4.1节可以知道当膜厚变化量ΔH较大时会对相对精度有较大影响，当迭代次数变大，ΔH减小到一定程度时影响可以忽略，故主要取迭代次数在8～20次的ΔH不同迭代次数与减小幅度两方面对收敛精度进行分析。

由图5(a)可知，在ΔH的变化周期不变的情况下，迭代后减小幅度越大，其整体收敛精度越小，并且出现两次局部最小值的迭代次数之差变小，但减小幅度过小后整体收敛精度变小不明显。由图5(b)可知，在迭代后ΔH减小幅度相同的情况下，变化周期越大，其整体相对精度越大，并且局部最小值出现的相隔次数也越大。因此，在选取ΔH的变化形式时选取变化周期较小(大于等于3次)，减小幅度较大时，可以得到较好的迭代收敛曲线。

由于得到较好的收敛曲线后，此时不一定能达到收敛条件，在已知ΔH较大能有效影响压力迭代的相对精度，故在判断出压力迭代的次数值时，将算法回滚，并增加ΔH的值。图6中是在每三次迭代ΔH变为1/2，迭代至43次时达到极小值，但是相对精度为1.630 7×10-5，没达到压力收敛要求；并且由图4可知，即使迭代至100次，仍然达不到相对精度要求。此时在保证膜厚分布合理的条件下，增加第43次迭代的ΔH值，计算验证可得当ΔH变为原来的64倍时，能在此次迭代达到收敛条件，并且压力分布(如图7所示)符合一般的弹流润滑规律；同时，将原本至少迭代100次以上的程序，变为迭代不到50次就可以满足要求，且不陷入死循环。

图5 相对精度随迭代次数变化曲线

图6 ΔH扩大后的相对精度

图7 ΔH扩大64倍后第43次迭代压力分布

5 结论

(1)混合编程能兼顾FORTRAN的计算速度以及MatLab的程序灵活、数据易于处理的优点。

(2)弹流润滑数值计算当相对精度在较低水平重复时，可以认为是得到收敛解，但现有的算法，并未对其进行判断，并且由于收敛条件设置的不合理，很容易导致程序陷入死循环。

(3)以实际的算法计算为依据，通过对膜厚常数的变化分析，得出合理设置ΔH变化形式，并辅以相对精度极小点处进行ΔH变大处理，能在判断出压力收敛精度趋于稳定时，在最低点处满足收敛条件，完成数值计算，避免陷入死循环而提高计算效率。