基于MatLab與C/C++混合編程的彈流潤滑數(shù)值算法改進(jìn)*

莫易敏 熊 釗 王 駿2 胡 杰 洪 葉

(1.武漢理工大學(xué)機(jī)電學(xué)院 湖北武漢 430070；2.上汽通用五菱汽車股份有限公司 廣西柳州 545007)

彈流潤滑的計(jì)算與工程實(shí)際聯(lián)系緊密，存在于滾動(dòng)軸承、滑動(dòng)軸承、齒輪傳動(dòng)等多個(gè)領(lǐng)域[1]，其中主要的雷諾方程是二階偏微分方程，以往依靠解析方法求解，必須經(jīng)過許多簡化處理才能獲得近似解，因而理論計(jì)算具有很大的誤差[2]。當(dāng)今由于計(jì)算機(jī)技術(shù)的迅速發(fā)展，復(fù)雜的潤滑問題有可能進(jìn)行數(shù)值解算。目前，潤滑計(jì)算作為一種理論依據(jù)在機(jī)械設(shè)計(jì)與潤滑材料的選用分析中占有重要的地位[3]。

本文作者使用混合編程實(shí)現(xiàn)求解彈流潤滑方程，不考慮壓力迭代方法的影響下，探究了膜厚方程中膜厚常量變化量ΔH對迭代相對精度的影響，得出了較為合適的ΔH變化形式，并依據(jù)ΔH的大小對相對精度的影響程度，在整體壓力收斂精度趨于穩(wěn)定時(shí)，在極小相對精度處擴(kuò)大ΔH，使相對進(jìn)度減小，滿足收斂要求，提前收斂，從而實(shí)現(xiàn)了彈流潤滑的快速數(shù)值計(jì)算。

1 彈流潤滑基本方程及其離散化

如圖1所示，研究采用的是工程實(shí)際中最普遍的兩接觸體主平面重合橢圓接觸彈流潤滑模型[1]。

圖1 主平面重合的橢圓接觸模型

Reynolds方程為(只考慮x方向速度)[3]：

(1)

量綱一化并離散后：

式中：p為壓力，不同位置壓力不同，可以看作與x、y有關(guān)的二元函數(shù)；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。

不考慮溫變時(shí)，等溫情況下，潤滑油黏度、密度、膜厚變化因素只需要考慮壓力，黏度-壓力方程[4]為

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

(2)

式中:z一般取0.68。

密度-壓力方程為

(3)

式中：η0、ρ0分別為p=0時(shí)等溫條件下潤滑劑的黏度與密度。

膜厚方程[5]：

(4)

量綱一化后：

其中:v(x,y)為在油膜壓力下接觸表面產(chǎn)生的彈性形變方程[1]：

(5)

離散后：

式中：h0為膜厚常量，是定值；E是兩接觸表面的綜合彈性模量；V(i,j)任意點(diǎn)的彈性形變。

載荷方程[6]：

(6)

量綱一化并離散后有：

需要注意的是：在對上述方程量綱一化的過程中，須將初始化的橢圓接觸區(qū)長軸長作為量綱一化的基數(shù)，即量綱一化后參數(shù)a為長軸；進(jìn)行數(shù)值計(jì)算中只有含有偏微分與積分的項(xiàng)，需要進(jìn)行離散化[7]。

2 離散方程迭代求解方法的改進(jìn)

基本方程中離散化后不考慮溫度的情況下，需要求解的是2個(gè)獨(dú)立的方程，即Reynolds方程和載荷平衡方程。其中Reynolds方程由于壓力的迭代求解，載荷平衡方程由于膜厚常數(shù)調(diào)整的判定，文中選用的是較為簡便的差分法進(jìn)行壓力迭代，其他的壓力迭代法如直接迭代法和多重網(wǎng)格積分法等暫不做討論。

現(xiàn)有的膜厚常數(shù)調(diào)整方法單一，均為膜厚常數(shù)的變化量ΔH固定周期逐步減小(常用的是每10次壓力迭代，ΔH值減為原來的1/2)，方法固定，且對不同的參數(shù)沒有較好的針對性，只是一個(gè)經(jīng)驗(yàn)值。

本文作者將探究在不同的ΔH(膜厚常量變化量)的變化方式下相對精度的具體變化情況，尋找出壓力迭代整體的相對精度較低的ΔH變化方式。在探究得到的ΔH與相對精度的規(guī)律上，通過改變ΔH的操作，以期能夠降低原有的最小相對精度，從而實(shí)現(xiàn)算法的提前收斂，提高效率，使算法更具操作性。

3 MatLab與C/C++混合編程實(shí)現(xiàn)彈流潤滑計(jì)算

在數(shù)值、科學(xué)和工程計(jì)算領(lǐng)域，F(xiàn)ORTRAN、C/C++語言更多地用于底層函數(shù)開發(fā)、軟件開發(fā)、單片機(jī)控制等，其運(yùn)算速度很快，對大型數(shù)據(jù)的運(yùn)算效率高，尤其對高性能計(jì)算和并行化計(jì)算方面，但是對數(shù)據(jù)進(jìn)行畫圖、圖像處理分析等方面，實(shí)現(xiàn)起來比較復(fù)雜[8]，并且FORTRAN的運(yùn)算不夠靈活。而MatLab是第四代的解釋性編程語言，編程思維符合人的思維習(xí)慣，數(shù)值計(jì)算方面編程容易實(shí)現(xiàn)，并且有強(qiáng)大的數(shù)據(jù)、圖像處理工具箱，但是其對于循環(huán)運(yùn)算的速度很慢[9]。文中結(jié)合2種語言的編程優(yōu)點(diǎn)，在MatLab平臺下，將彈流潤滑算法中耗時(shí)的主要循環(huán)部分彈性變形矩陣用C/C++語言實(shí)現(xiàn)，并將其編譯成能被MatLab調(diào)用的動(dòng)態(tài)鏈接庫[10](.dll文件)。算法的具體實(shí)現(xiàn)如圖2所示。

圖2 混合編程彈流潤滑數(shù)值計(jì)算流程

彈流潤滑算例的參數(shù)取x、y方向上綜合曲率半徑分別為0.03、0.06 m，載荷w=1 000 N，初始黏度η0=0.02 Pa·s，兩接觸表面沿x方向平均速度us=3 m/s，綜合彈性模量E=221 GPa。計(jì)算區(qū)域量綱一化X起始與終點(diǎn)坐標(biāo)分別為-2.5、1.5，Y方向上則是-2、2。

通過MatLab計(jì)時(shí)工具tic、toc指令，分別記錄MatLab與調(diào)用C/C++編譯成的dll文件實(shí)現(xiàn)計(jì)算變形矩陣VI每次迭代運(yùn)行時(shí)間，結(jié)果如表1所示。

變形矩陣VI具體是計(jì)算求解域內(nèi)65×65等距節(jié)點(diǎn)的壓力彈性變形，4層循環(huán)運(yùn)算次數(shù)達(dá)654次，而由運(yùn)行情況來說，10次迭代，2種方法結(jié)果一樣，但運(yùn)行時(shí)間相差巨大。因此，若有需要，將MatLab中耗時(shí)的改為C/C++程序，能顯著減少計(jì)算時(shí)間。

表1 不同編程方法計(jì)算時(shí)間對比

4 數(shù)值計(jì)算算法分析與改進(jìn)

4.1 算法收斂影響因素分析

彈流潤滑壓力、膜厚分布求解的可靠性，主要取決于結(jié)果是否滿足壓力迭代收斂條件與載荷平衡方程，而壓力迭代算法的選取與膜厚常數(shù)最終取值是算法的主要影響因素[7]。

求解彈流潤滑方面的問題，需聯(lián)立求解各類復(fù)雜的方程，涉及流體潤滑、載荷方程、固體表面彈性變形方程、潤滑劑黏度方程及密度方程，系統(tǒng)非線性強(qiáng)，不能得出解析解。目前，對該類問題的計(jì)算方法有逆解法、牛頓有限元法、復(fù)合直接迭代法、牛頓有限差分法、直接迭代法和多重網(wǎng)格積分法[11]。

文中采用牛頓有限差分法，而膜厚常數(shù)方面一般根據(jù)Hamrock-Dowson最小膜厚的經(jīng)驗(yàn)公式：

算出最小膜厚并量綱一化為Hmin，此時(shí)可以通過Hmin反求出膜厚常數(shù)H0的初值，之后根據(jù)載荷方程平衡調(diào)整膜厚常數(shù)：H0=H0±ΔH，此處一般取ΔH=0.1～0.05Hmin，并且每固定次迭代，ΔH的值就變?yōu)樵瓉淼?/2[9]。

壓力迭代收斂準(zhǔn)則采用的是相對精度判斷準(zhǔn)則：

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

圖3所示是ΔH=0.1Hmin，每3次迭代減半，迭代150次，取得的最優(yōu)結(jié)果，為了顯示清晰，只截取了X方向 -1.5、1.5?？芍簤毫Ψ植汲霈F(xiàn)了二次應(yīng)力峰，但是由于材料剛度、尺寸較大，其應(yīng)力峰的值并不大，因此出口處膜厚減小現(xiàn)象也不明顯，但是計(jì)算結(jié)果是符合彈流潤滑的一般特征的。

圖3 彈流潤滑壓力與膜厚分布

圖4示出了不同ΔH變化下相對精度隨迭代次數(shù)變化規(guī)律。由不同情況下100次迭代曲線可以看出來，ΔH的不同變化形式是不會(huì)顯著影響收斂相對精度的總體變化趨勢，并且隨著ΔH的不斷減小，ΔH對壓力迭代的影響可以忽略不計(jì)，此時(shí)主要影響因素歸于迭代算法的特性。圖4中在迭代次數(shù)大于30次后，相對精度在10-3～10-6范圍內(nèi)周期變化，變化周期大致為每22次迭代相對精度并未發(fā)散，這一特征符合壓力迭代算法收斂的要求[8]。

圖4 不同ΔH變化下相對精度隨迭代次數(shù)變化曲線

4.2 膜厚變化對相對精度影響分析與改進(jìn)

由4.1節(jié)可以知道當(dāng)膜厚變化量ΔH較大時(shí)會(huì)對相對精度有較大影響，當(dāng)?shù)螖?shù)變大，ΔH減小到一定程度時(shí)影響可以忽略，故主要取迭代次數(shù)在8～20次的ΔH不同迭代次數(shù)與減小幅度兩方面對收斂精度進(jìn)行分析。

由圖5(a)可知，在ΔH的變化周期不變的情況下，迭代后減小幅度越大，其整體收斂精度越小，并且出現(xiàn)兩次局部最小值的迭代次數(shù)之差變小，但減小幅度過小后整體收斂精度變小不明顯。由圖5(b)可知，在迭代后ΔH減小幅度相同的情況下，變化周期越大，其整體相對精度越大，并且局部最小值出現(xiàn)的相隔次數(shù)也越大。因此，在選取ΔH的變化形式時(shí)選取變化周期較小(大于等于3次)，減小幅度較大時(shí)，可以得到較好的迭代收斂曲線。

由于得到較好的收斂曲線后，此時(shí)不一定能達(dá)到收斂條件，在已知ΔH較大能有效影響壓力迭代的相對精度，故在判斷出壓力迭代的次數(shù)值時(shí)，將算法回滾，并增加ΔH的值。圖6中是在每三次迭代ΔH變?yōu)?/2，迭代至43次時(shí)達(dá)到極小值，但是相對精度為1.630 7×10-5，沒達(dá)到壓力收斂要求；并且由圖4可知，即使迭代至100次，仍然達(dá)不到相對精度要求。此時(shí)在保證膜厚分布合理的條件下，增加第43次迭代的ΔH值，計(jì)算驗(yàn)證可得當(dāng)ΔH變?yōu)樵瓉淼?4倍時(shí)，能在此次迭代達(dá)到收斂條件，并且壓力分布(如圖7所示)符合一般的彈流潤滑規(guī)律；同時(shí)，將原本至少迭代100次以上的程序，變?yōu)榈坏?0次就可以滿足要求，且不陷入死循環(huán)。

圖5 相對精度隨迭代次數(shù)變化曲線

圖6 ΔH擴(kuò)大后的相對精度

圖7 ΔH擴(kuò)大64倍后第43次迭代壓力分布

5 結(jié)論

(1)混合編程能兼顧FORTRAN的計(jì)算速度以及MatLab的程序靈活、數(shù)據(jù)易于處理的優(yōu)點(diǎn)。

(2)彈流潤滑數(shù)值計(jì)算當(dāng)相對精度在較低水平重復(fù)時(shí)，可以認(rèn)為是得到收斂解，但現(xiàn)有的算法，并未對其進(jìn)行判斷，并且由于收斂條件設(shè)置的不合理，很容易導(dǎo)致程序陷入死循環(huán)。

(3)以實(shí)際的算法計(jì)算為依據(jù)，通過對膜厚常數(shù)的變化分析，得出合理設(shè)置ΔH變化形式，并輔以相對精度極小點(diǎn)處進(jìn)行ΔH變大處理，能在判斷出壓力收斂精度趨于穩(wěn)定時(shí)，在最低點(diǎn)處滿足收斂條件，完成數(shù)值計(jì)算，避免陷入死循環(huán)而提高計(jì)算效率。