基于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變大處理，能在判斷出壓力收斂精度趨于穩定時，在最低點處滿足收斂條件，完成數值計算，避免陷入死循環而提高計算效率。