蒲偉,王家序,,周廣武,李俊陽(yáng),褚坤明

(1.四川大學(xué)空天科學(xué)與工程學(xué)院, 610065, 成都; 2.重慶大學(xué)機(jī)械傳動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 404100, 重慶)

卷吸速度方向與橢圓短軸成一夾角的彈流潤(rùn)滑漸近網(wǎng)格加密算法

蒲偉1,王家序1,2,周廣武1,李俊陽(yáng)2,褚坤明1

(1.四川大學(xué)空天科學(xué)與工程學(xué)院, 610065, 成都; 2.重慶大學(xué)機(jī)械傳動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 404100, 重慶)

針對(duì)數(shù)值求解卷吸速度方向與橢圓短軸成一夾角的重載薄膜彈流潤(rùn)滑問(wèn)題,提出將漸進(jìn)網(wǎng)格加密(progressive mesh densification,PMD)算法的思想應(yīng)用于此類問(wèn)題的求解,建立了此類問(wèn)題雷諾方程的離散差分格式,運(yùn)用快速傅里葉變換算法計(jì)算點(diǎn)接觸彈性變形,采用PMD網(wǎng)格跳轉(zhuǎn)方法加速迭代收斂,最后得到此類問(wèn)題的完全數(shù)值解和油膜形狀分布圖,并與前人的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比分析。結(jié)果表明:在重載薄膜潤(rùn)滑條件下,數(shù)值求解結(jié)果與前人的實(shí)驗(yàn)結(jié)果相一致;與直接在最終網(wǎng)格上進(jìn)行迭代相比,文中提出的方法可以在粗糙網(wǎng)格上迅速消除誤差,為高密度網(wǎng)格提供更好的初值,從而使計(jì)算收斂速度提高1倍以上,證明了PMD算法在求解此類問(wèn)題時(shí)具有良好的收斂性和有效性。

彈流潤(rùn)滑;漸進(jìn)網(wǎng)格加密算法;卷吸速度;橢圓接觸

弧齒錐齒輪和準(zhǔn)雙曲面齒輪等傳動(dòng)由于具有承載能力大、傳動(dòng)平穩(wěn)等優(yōu)點(diǎn),已廣泛應(yīng)用于航空、航天、車輛、工程機(jī)械、礦山機(jī)械、機(jī)床等工程領(lǐng)域。此類傳動(dòng)件與普通圓柱齒輪及滾動(dòng)軸承等傳動(dòng)件相比,潤(rùn)滑模型具有特殊性,其卷吸速度方向與接觸橢圓短軸并不重合,而是成一夾角。因此,研究卷吸速度方向與橢圓短軸成一夾角的彈流潤(rùn)滑的高效、穩(wěn)定數(shù)值解法,對(duì)揭示此類傳動(dòng)界面的潤(rùn)滑機(jī)理具有十分重要的意義。

在以往的點(diǎn)接觸彈流潤(rùn)滑研究中,絕大部分學(xué)者針對(duì)卷吸速度方向與接觸橢圓短軸重合的情況進(jìn)行了深入的研究,并取得了一系列的研究成果[1-4]。然而,針對(duì)卷吸速度方向與橢圓短軸成一夾角的彈流潤(rùn)滑問(wèn)題,研究進(jìn)展則較為緩慢,特別是此類重載彈流潤(rùn)滑問(wèn)題的穩(wěn)定數(shù)值解法一直是彈流領(lǐng)域的挑戰(zhàn)。Chittenden等曾采用直接迭代法分析了卷吸速度方向與橢圓短軸成一夾角的情況,但該方法只能針對(duì)輕載荷的情況進(jìn)行分析[5]。Wang等采用多重網(wǎng)格法對(duì)此類問(wèn)題進(jìn)行了分析,得到了此類問(wèn)題的完全數(shù)值解[6]。雖然多重網(wǎng)格法能很快消除迭代誤差,加速計(jì)算收斂,但針對(duì)重載薄膜潤(rùn)滑問(wèn)題時(shí),收斂性卻并不理想,而且多重網(wǎng)格法在迭代計(jì)算時(shí)通過(guò)在粗糙網(wǎng)格和精細(xì)網(wǎng)格之間反復(fù)轉(zhuǎn)換來(lái)達(dá)到收斂精度,但由于在重載薄膜潤(rùn)滑時(shí)油膜厚度對(duì)網(wǎng)格密度變得極為敏感,因此在每一層網(wǎng)格所求得的收斂值之間差距變大,導(dǎo)致在反復(fù)跳轉(zhuǎn)迭代時(shí)誤差并不能減小,迭代時(shí)容易導(dǎo)致溢出或死循環(huán)[7]。Stahl等近期發(fā)表了卷吸速度方向與橢圓主軸成一夾角的情況的實(shí)驗(yàn)解,得到了幾組油膜分布圖及中心油膜厚度,結(jié)果表明卷吸速度方向?qū)τ湍ば螤罴昂穸染酗@著影響[8]。漸近網(wǎng)格加密(progressive mesh densification,PMD)算法是朱東等在求解卷吸速度方向與橢圓短軸重合的潤(rùn)滑問(wèn)題時(shí)提出的,該方法能夠在保證數(shù)值精度的情況下顯著加快求解過(guò)程,特別是對(duì)于重載超薄膜潤(rùn)滑[7]。但是,目前該算法僅局限于對(duì)卷吸速度方向與橢圓短軸重合的潤(rùn)滑問(wèn)題的求解。

本文旨在將PMD算法的思想應(yīng)用于求解卷吸速度方向與橢圓短軸成一夾角的彈流潤(rùn)滑問(wèn)題,并將該算法的數(shù)值計(jì)算結(jié)果與Stahl等發(fā)表的實(shí)驗(yàn)數(shù)據(jù)[8]進(jìn)行對(duì)比分析,以驗(yàn)證該算法在求解重載條件下此類潤(rùn)滑問(wèn)題時(shí)的收斂性和正確性。

1 數(shù)學(xué)模型

1.1 接觸幾何模型

弧齒錐齒輪和準(zhǔn)雙曲面齒輪等傳動(dòng)的接觸潤(rùn)滑問(wèn)題可簡(jiǎn)化為如圖1所示的橢圓接觸模型,卷吸速度方向與橢圓的短軸成一夾角,a、b分別為赫茲接觸橢圓的長(zhǎng)、短軸,ue為綜合卷吸速度,u、v為ue在x、y方向的分量。

圖1 接觸潤(rùn)滑的橢圓接觸幾何模型

1.2潤(rùn)滑基本方程[8]

等溫條件下的點(diǎn)接觸非穩(wěn)態(tài)彈流潤(rùn)滑歸一化雷諾方程為

(1)

式中:u=uecosθ;v=uesinθ。

油膜的幾何方程為

(2)

彈性變形方程為

壓力黏度方程為

η=η0eαp

(3)

壓力密度方程為

(4)

載荷平衡方程為

(5)

1.3 方程的歸一化

將方程(1)～(5)中的量按下式進(jìn)行歸一化

歸一化后的雷諾方程為

(6)

式中

其他歸一化方程皆由此方法推出,不再贅述。

2 數(shù)值計(jì)算

PMD算法的思想是將離散的潤(rùn)滑方程組在較小密度的網(wǎng)格(A網(wǎng)格)上進(jìn)行迭代(見(jiàn)圖2),當(dāng)達(dá)到設(shè)定的收斂精度后,立即跳轉(zhuǎn)到較大密度的網(wǎng)格(B網(wǎng)格)上進(jìn)行迭代,在B網(wǎng)格根據(jù)A網(wǎng)格上迭代得到的油膜厚度和油膜壓力值,通過(guò)插值的方法確定其他未知節(jié)點(diǎn)的油膜壓力和油膜厚度,并以此作為B網(wǎng)格上迭代的初值,當(dāng)在B網(wǎng)格上迭代達(dá)到收斂精度后,再如上述方法進(jìn)行跳轉(zhuǎn),直到在預(yù)先設(shè)定的最終網(wǎng)格(最密網(wǎng)格)上得到收斂解。

圖2 PMD算法示意圖

方程(1)～(5)中,彈性變形、潤(rùn)滑油黏度和密度都是壓力的函數(shù),將其與雷諾方程合并成一個(gè)微分積分方程組,采用復(fù)合迭代法[9]進(jìn)行求解,其中彈性變形的計(jì)算采用的是快速傅里葉變換(FFT)方法[10]。求解時(shí),在每一個(gè)未知節(jié)點(diǎn)壓力處,將該微分積分方程組離散成一個(gè)差分方程[11],表達(dá)式為

Ai,jPi-1,j+Bi,jPi,j+Ci,jPi+1,j=Fi,j

(7)

式中:Ai,j、Bi,j、Ci,j、Fi,j為已知項(xiàng),由已經(jīng)求解出的壓力分布算出;Pi-1,j、Pi,j、Pi+1,j是本次迭代需求解的未知量。

采用二階中心差分方法,將雷諾方程中左邊的壓力流項(xiàng)離散[11]

式中

將離散的壓力流拆分,按照式(3)的格式進(jìn)行改寫,并分散到各項(xiàng)系數(shù)中,結(jié)果如下

對(duì)雷諾方程右邊的2項(xiàng)都分別采用一階向后差分,然后按照式(7)的格式進(jìn)行改寫,并將其分散到各項(xiàng)系數(shù)中進(jìn)行疊加,方法與上述相同,此處不再贅述。

未知節(jié)點(diǎn)的PD采用拋物插值方法求出

式中:X、X0、X1、X2分別為PD、PA、PB、PC在高密度網(wǎng)格中的橫坐標(biāo)。

未知節(jié)點(diǎn)PE的計(jì)算方法與PD的相同。采用上述插值方法,可將高密度網(wǎng)格上各節(jié)點(diǎn)的油膜壓力和油膜厚度全部計(jì)算出來(lái),然后進(jìn)行迭代,直到達(dá)到設(shè)定的收斂精度。設(shè)定的收斂精度可取為

重復(fù)以上過(guò)程,直到達(dá)到設(shè)定的網(wǎng)格密度。由于網(wǎng)格密度每增加1倍,計(jì)算時(shí)間將會(huì)增加幾倍,綜合考慮計(jì)算時(shí)間和收斂精度,建議采用圖3所示的3種網(wǎng)格密度進(jìn)行迭代。

圖3 用于迭代計(jì)算的3種網(wǎng)格密度

在迭代過(guò)程中,壓力松弛因子ωp一般取0.005～0.3,在載荷大、膜厚小的情況下,ωp應(yīng)取相對(duì)小的值;H0的修正因子ωh一般取0.001～0.01。計(jì)算流程圖如圖4所示。

圖4 計(jì)算流程圖

3 結(jié)果與討論

為了驗(yàn)證本算法在求解重載薄膜潤(rùn)滑時(shí)的收斂性和正確性,將本算法的數(shù)值仿真結(jié)果與Stahl等的實(shí)驗(yàn)結(jié)果[8]進(jìn)行比較。數(shù)值仿真的相關(guān)輸入?yún)?shù)與文獻(xiàn)[8]中的實(shí)驗(yàn)參數(shù)保持一致,基本參數(shù)如下:θ={63.43°,73.90°,82.37°,90.00°},pH=1.03GPa,ue=0.6 m/s,Rx=4mm,Ry=12.7 mm,k=0.46;數(shù)值求解域?yàn)?2≤X≤2,-2≤Y≤2;采用3層網(wǎng)格,低密度網(wǎng)格參數(shù)為64×64,收斂精度為10-6;中密度網(wǎng)格參數(shù)為128×128,收斂精度為10-5;高密度網(wǎng)格參數(shù)為256×256,收斂精度為10-5;高密度網(wǎng)格有66 049個(gè)節(jié)點(diǎn)。數(shù)值仿真結(jié)果見(jiàn)表1～表4,對(duì)比結(jié)果如圖5所示。

表1θ=63.43°時(shí)每一層網(wǎng)格上迭代收斂后的油膜厚度值

網(wǎng)格參數(shù)中心油膜厚度/nm最小油膜厚度/nm64×6420438128×12822860256×25623465

表2θ=73.90°時(shí)每一層網(wǎng)格上迭代收斂后的油膜厚度值

網(wǎng)格參數(shù)中心油膜厚度/nm最小油膜厚度/nm64×6421535128×12824354256×25625159

表3θ=82.37°時(shí)每一層網(wǎng)格上迭代收斂后的油膜厚度值

網(wǎng)格參數(shù)中心油膜厚度/nm最小油膜厚度/nm64×6422631128×12826753256×25627758

表4θ=90.00°時(shí)每一層網(wǎng)格上迭代收斂后的油膜厚度值

網(wǎng)格參數(shù)中心油膜厚度/nm最小油膜厚度/nm64×6423147128×12827955256×25629263

(a)油膜形狀對(duì)比

(b)中心油膜厚度對(duì)比

(c)最小油膜厚度對(duì)比

從表1～表4中可以看出,高密度網(wǎng)格層上迭代得到的油膜厚度值大于低密度網(wǎng)格層上迭代所得的油膜厚度值,且隨著網(wǎng)格的逐漸加密,油膜厚度的增量逐漸減小,這與朱東等在研究卷吸速度方向與橢圓短軸重合時(shí)得出的結(jié)論[7]是一致的。綜合考慮計(jì)算精度和計(jì)算時(shí)間,網(wǎng)格密度沒(méi)有再進(jìn)一步增大。

圖5a為油膜形狀對(duì)比,數(shù)值仿真圖形采用了顏色漸變的方式來(lái)繪制。文獻(xiàn)[8]中采用的橢球試件的表面粗糙度Ra=0.03μm。根據(jù)文獻(xiàn)[12]可知,表面粗糙度對(duì)油膜厚度有一定影響,影響程度與表面的紋理和粗糙度大小有關(guān)。文獻(xiàn)[8]的實(shí)驗(yàn)所用橢球試件的粗糙度較小,故對(duì)油膜厚度的影響程度也較小,可以用于驗(yàn)證彈流潤(rùn)滑模型的光滑解。從圖5中可以看出,油膜形狀隨著卷吸速度方向的變化而變化,沿著潤(rùn)滑油流動(dòng)的方向,彈流潤(rùn)滑出口區(qū)的間隙呈發(fā)散狀,導(dǎo)致接觸區(qū)的兩側(cè)面壓力急劇下降而低于赫茲壓力,使得固體表面產(chǎn)生突起,隨著潤(rùn)滑油流動(dòng)方向的改變,固體表面突起的位置不斷變化,所以形成的油膜形狀也隨著流動(dòng)方向發(fā)生變化。數(shù)值仿真得到的油膜形狀及油膜厚度值都與文獻(xiàn)[8]的實(shí)驗(yàn)結(jié)果相一致。算例的最大赫茲接觸壓力pH達(dá)到了1.03GPa,油膜厚度最小為40nm,收斂精度達(dá)到了10-5。

當(dāng)卷吸速度方向與橢圓短軸成一夾角時(shí),從前面的差分方程中可以看出,與卷吸速度方向重合于橢圓主軸時(shí)相比,離散方程(7)中多了雷諾方程中沿y方向的剪切流,從而加大了迭代系數(shù)矩陣的擾動(dòng),使得此聯(lián)立方程組具有極強(qiáng)的非線性,而以往的數(shù)值方法在分析重載工況時(shí)收斂性并不理想。上述驗(yàn)證計(jì)算表明,PMD算法有效解決了這一難題,在粗糙的網(wǎng)格上能迅速消除迭代誤差,為高密度網(wǎng)格提供更好的初值,從而可加速收斂。因此,針對(duì)卷吸速度方向與橢圓短軸成一夾角的彈流潤(rùn)滑,PMD算法在求解重載薄膜潤(rùn)滑時(shí)具有良好的加速收斂特性。

此外,為了進(jìn)一步驗(yàn)證本算法的加速收斂性能,對(duì)2組算例進(jìn)行了對(duì)比分析:一組算例使用PMD算法,采用3層網(wǎng)格,網(wǎng)格參數(shù)與前述例子相同;另一組算例直接在256×256的網(wǎng)格上進(jìn)行迭代?；据斎?yún)?shù)如下:θ=50°,pH=2.02GPa,ue=0.5～15m/s,Rx=12.7 mm,Ry=36.1 mm,k=2。在相同收斂精度(10-5)下,計(jì)算結(jié)果如表5、表6、圖6所示,其中t1是在256×256的網(wǎng)格上直接迭代的計(jì)算時(shí)間,t2是采用PMD方法的計(jì)算時(shí)間。

表5 PMD算法的計(jì)算結(jié)果

表6 相對(duì)計(jì)算時(shí)間

(a)ue=15m/s (b)ue=5m/s (c)ue=0.5m/s

從表6中可以看出,與直接迭代相比,采用PMD方法可使計(jì)算速度提高1倍以上,并能夠顯著加快收斂速度。此外,這2組算例的最大赫茲接觸壓力達(dá)到了2.02GPa,最小油膜厚度低至8 nm,進(jìn)一步驗(yàn)證了PMD算法在求解重載薄膜潤(rùn)滑時(shí)具有良好的收斂性和有效性。

4 結(jié) 論

(1)本文給出了針對(duì)卷吸速度方向與橢圓短軸成一夾角時(shí)彈流潤(rùn)滑方程組的差分格式,得出了相應(yīng)的PMD數(shù)值求解算法,包括插值、迭代跳轉(zhuǎn)等具體步驟。

(2)將PMD算法求解的數(shù)值結(jié)果與Stahl等的實(shí)驗(yàn)研究結(jié)果[8]進(jìn)行了比較,結(jié)果表明二者相當(dāng)一致,證明對(duì)于重載薄膜潤(rùn)滑,PMD算法在求解卷吸速度方向與橢圓短軸成一夾角的彈流潤(rùn)滑問(wèn)題時(shí)具有良好的收斂性和有效性。

(3)與直接在最終網(wǎng)格上進(jìn)行迭代相比,針對(duì)卷吸速度方向與橢圓短軸成一夾角的彈流潤(rùn)滑問(wèn)題,采用PMD算法可以在粗糙網(wǎng)格上迅速消除誤差,為高密度網(wǎng)格提供更好的初值,從而顯著提高計(jì)算收斂速度。

[1] ZHU Dong, WANG Jane Q. Elastohydrodynamic lubrication: a gateway to interfacial mechanics: review and prospect [J]. Journal of Tribology, 2011, 133(4): 1-14.

[2] HAMROCK B J, DOWSON D. Isothermal elastohydro-dynamic lubrication of point contacts: part I Theoretical formulation [J]. Journal of Lubrication Technology, 1976, 98(2): 223-228.

[3] HAMROCK B J, DOWSON D. Isothermal elastohydro-dynamic lubrication of point contacts: part II Ellipticity parameter results [J]. Journal of Lubrication Technology, 1976, 98(3): 375-381.

[4] 朱東, 溫詩(shī)鑄, 鄭林慶. 熱彈性流體動(dòng)力潤(rùn)滑橢圓接觸問(wèn)題的完全數(shù)值解 [J]. 機(jī)械工程學(xué)報(bào), 1984, 20(1): 1-15.

ZHU Dong, WEN Shizhu, ZHENG Linqing. A full numerical solution for the thermoelastohydynamic problem in elliptical contacts [J]. Chinese Journal of Mechanical Engineering, 1984, 20(1): 1-15.

[5] CHITTENDEN R J, DOWAON D, DUNN J F. A theoretical analysis of the isothermal elastohydrodynamic lubrication of concentrated contacts: II General case, with lubricant entrainment along either principal axis of the Hertzian contact ellipse or at some intermediate angle [J]. Proceedings of the Royal Society of London: A Mathematical and Physical Sciences, 1985, 397(1813): 271-294.

[6] WANG Jing, QU Shiyue, YANG Peiran. Simplified multigrid technique for the numerical solution to the steady-state and transient EHL line contacts and the arbitrary entrainment EHL point contacts [J]. Tribology International, 2001, 34(3): 191-202.

[7] ZHU Dong. On some aspects in numerical solution of thin-film and mixed EHL [J]. Proc Inst Mech Eng: Part J J Eng Tribol, 2007, 221(5): 561-579.

[8] STAHL K, MICHAELIS K, MAYER J. Theoretical and experimental investigations on EHL point contacts with different entrainment velocity directions [J]. Tribology Transactions, 2013, 56(5): 728-738.

[9] AI Xiaolan. Numerical analysis of elastohydrodynamically lubricated line and point contacts with rough surfaces by using semi-system and multigrid methods [D]. Evanston, IL, USA: Northwestern University, 1993.

[10]LIU Shuangbiao, WANG Qian, LIU Geng. A versatile method of discrete convolution and FFT (DC-FFT) for contact analysis [J]. Wear, 2000, 243(1): 101-110.

[11]REN Ning. Advanced modeling of mixed lubrication and its mechanical and biomedical applications [D]. Evanston, IL, USA: Northwestern University, 2009.

[12]ZHU Dong, HU Yuan-Zhong. Effects of rough surface topography and orientation on the characteristics of EHD and mixed lubrication in both circular and elliptical contacts [J]. Tribology Transactions, 2001, 44(3): 391-398.

(編輯 葛趙青)

ProgressiveMeshDensification(PMD)MethodinElastohydrodynamicLubricationofEllipticalContactswithArbitraryEntrainment

PU Wei1,WANG Jiaxu1,2,ZHOU Guangwu1,LI Junyang2,CHU Kunming1

To numerically solve the elastohydrodynamic lubrication problem that entraining velocity vector has an oblique angle with the minor axis of Hertzian contact ellipse under heavy loading and thin film lubrication conditions, the progressive mesh densification (PMD) method is used, and a discrete differential scheme for Reynolds equation is set up. The fast Fourier transform is adopted to calculate the elastic deformation for point contact, and the PMD method is used to accelerate the iterative convergence. Moreover, a complete numerical solution and the film thickness distribution are obtained and compared with the experimental results from literatures to validate the model. The results show that, under heavy loading and thin film lubrication conditions, this numerical solution is in accordance with the experimental results, proving that the PMD method is effective and has fine convergent behavior for solving this kind of problem.

elastohydrodynamic lubrication; progressive mesh densification method; entraining velocity; elliptical contact

2014-02-08。

蒲偉(1989—),男,博士生;王家序(通信作者),男,教授,博士生導(dǎo)師。

國(guó)家自然科學(xué)基金資助項(xiàng)目(51375506);教育部高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金資助項(xiàng)目(2012018113001);中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目(2014SCU11009)。

時(shí)間:2014-06-13

10.7652/xjtuxb201409016

TH132

:A

:0253-987X(2014)09-0095-06

網(wǎng)絡(luò)出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20140613.1457.003.html