胡 楊 孟永鋼

(1.上海大學(xué)機(jī)電工程與自動化學(xué)院 上海 200072；2.清華大學(xué)摩擦學(xué)國家重點(diǎn)實(shí)驗室 北京 100084)

可傾瓦徑向滑動軸承由于其良好的穩(wěn)定性，在大型高速渦輪機(jī)械中得到了廣泛的應(yīng)用。當(dāng)軸徑表面線速度較低時，可傾瓦滑動軸承運(yùn)行在層流狀態(tài)，可以采用DOWSON[1]提出的廣義雷諾方程來預(yù)測軸承性能。

然而，在大型渦輪機(jī)械中，由于轉(zhuǎn)子直徑大、轉(zhuǎn)速高，軸頸表面線速度可達(dá)數(shù)百米每秒，此時，可傾瓦滑動軸承運(yùn)行在湍流狀態(tài)。NG和PAN[2]提出了剪切驅(qū)動理論來模擬薄膜流動中的湍流現(xiàn)象。基于NG和PAN[2]的線性化湍流理論，TANIGUCHI等[3]對大型可傾瓦滑動軸承進(jìn)行了熱流體動力學(xué)分析，瓦面溫度、摩擦損失和偏心率的實(shí)驗結(jié)果與仿真結(jié)果吻合較好。

當(dāng)負(fù)載較大時，施加在瓦塊表面上的壓力使瓦塊產(chǎn)生彈性變形；另一方面，溫升也會使瓦塊產(chǎn)生較大的熱變形。軸瓦變形會引起油膜分布的變化，從而影響潤滑性能。因此，在重載條件下，有必要考慮軸瓦變形對潤滑性能的影響。ETTLES[4]在預(yù)測小型軸承的潤滑性能時，采用了一維梁彎曲理論來評估軸瓦變形。雖然一維[4]和二維瓦塊變形模型[5]具有較高的計算效率，但是當(dāng)可傾瓦滑動軸承的尺寸較大時，為了確保精度，采用三維瓦塊變形模型[6]是必要的。SANO等[7]通過有限元方法評估瓦塊變形，將TANIGUCHI等[3]的熱流體動力潤滑模型擴(kuò)展為熱彈流體動力潤滑模型，并基于所構(gòu)建模型，對直徑890 mm的直接潤滑的兩瓦可傾瓦滑動軸承進(jìn)行了研究。ARIHARA等[8]構(gòu)建了熱彈流體動力潤滑模型，對高速重載條件下的可傾瓦滑動軸承的靜態(tài)性能進(jìn)行了研究，并與實(shí)驗結(jié)果進(jìn)行了對比分析。MERMERAS等[9]描述了一種先進(jìn)的建模方法，并對新型900 mm三瓦可傾瓦徑向滑動軸承進(jìn)行優(yōu)化設(shè)計和驗證。HAGEMANN等[10]在不同供油流量條件下對定向潤滑五瓦可傾瓦滑動軸承進(jìn)行了實(shí)驗研究，通過與理論結(jié)果的對比發(fā)現(xiàn)速度對從富油到貧油的轉(zhuǎn)變影響很大。SAN ANDRéS等[11]評估了在潤滑油量過高和過低時可傾瓦滑動軸承的穩(wěn)態(tài)和動態(tài)性能。

上述可傾瓦滑動軸承的熱彈流體動力潤滑模型大多是采用有限差分法或有限元法進(jìn)行自主編程求解。由于編程代碼非開源和不易復(fù)制，限制了仿真工具在工業(yè)中的廣泛應(yīng)用。使用商業(yè)軟件可以使研究人員更多地關(guān)注物理現(xiàn)象，而不是復(fù)雜的數(shù)值求解方法。LOHNER等[12]通過結(jié)合COMSOL Multi-physics和MATLAB實(shí)現(xiàn)了HABCHI[13]的方法，提供了一種求解涂層圓柱線接觸的熱彈流體動力模型的新思路。受這些工作的啟發(fā)，本文作者采用了一種結(jié)合自主編程和商用軟件優(yōu)點(diǎn)的仿真方法。

本文作者采用COMSOL和MATLAB相結(jié)合的方法，建立了考慮湍流的可傾瓦滑動軸承的高效熱彈流體動力潤滑模型。由于計算流體力學(xué)計算量大，且在惡劣工況下易發(fā)生失穩(wěn)，文中所建立的模型仍然是基于雷諾方程。模型中，采用COMSOL中的偏微分方程物理場求解考慮湍流的修正雷諾方程和能量守恒方程以及熱傳導(dǎo)方程，采用COMSOL中的熱應(yīng)力物理場求解壓力和溫度引起的瓦塊變形，利用MATLAB自編程序來實(shí)現(xiàn)子模型之間的耦合。在此基礎(chǔ)上，基于所構(gòu)建模型，研究了軸承瓦塊弧長、名義間隙、支承點(diǎn)偏移量和軸向?qū)挾葘蓛A瓦滑動軸承性能包括油膜厚度、壓力、瓦面溫度、摩擦功耗等的影響。

1 考慮湍流的可傾瓦滑動軸承熱彈流體動力潤滑模型

1.1 油膜厚度描述

文中以富油潤滑的四瓦可傾瓦滑動軸承作為研究對象，可傾瓦滑動軸承的示意圖如圖1(a)所示，坐標(biāo)系如圖1(b)所示。油膜厚度表達(dá)式如下：

圖1 可傾瓦滑動軸承示意與坐標(biāo)系

h(θ,z)=c+eXcosθ+eYsinθ-m·c·cos(θ-

θp)+δp·(R+tp)·sin(θ-θp)+δd(θ,z)

(1)

式中：c為徑向間隙；m為預(yù)緊；θp為支承點(diǎn)位置；δp為瓦塊傾斜角；R為轉(zhuǎn)子半徑；tp為瓦塊表面與支承點(diǎn)位置的距離；δd為瓦塊變形引起的油膜厚度變化；eX和eY分別為沿X和Y方向軸心位置分量。

1.2 考慮湍流的修正雷諾方程

為了考慮湍流效應(yīng)的影響，采用NG和PAN[2]推導(dǎo)的線性湍流潤滑理論。修正雷諾方程公式表達(dá)式如下：

(2)

其中，G1、G2和F1的表達(dá)式如下：

(3a)

(3b)

(3c)

(4)

式中：c為徑向間隙；N為轉(zhuǎn)速(r/min)；ψ為間隙比。量綱一化坐標(biāo)系(θ,η,ζ)定義如下：

θ=x/R,η=y/h,ζ=z/(L/2)

(5)

ξ1、ξ2、ξ3和ξ4定義如下：

(6a)

(6b)

(6c)

(6d)

式中：fc(y)和gc(y)為湍流函數(shù)，具體計算可參考文獻(xiàn)[3]。

1.3 三維能量方程

湍流狀態(tài)下，量綱一化的三維能量方程表示如下：

(7)

能量方程中量綱一化因子定義如下：

(8)

量綱一化的速度場表達(dá)如下：

(9a)

(9b)

(9c)

1.4 熱傳導(dǎo)方程

量綱一化柱坐標(biāo)下的瓦塊熱傳導(dǎo)方程表達(dá)式如下：

(10)

1.5 瓦塊變形

當(dāng)軸承負(fù)載很高時，需要考慮瓦塊的變形。瓦塊的變形包括兩部分，一部分是由于瓦面上油膜壓力造成，另一部分是由于瓦塊內(nèi)溫度分布不一致造成。瓦塊變形量ε的表達(dá)式[14]如下:

ε=D-1σ+αT(Tp-Tref)

(11)

式中：σ為應(yīng)力張量；D是各向同性材料假設(shè)下的彈性剛度矩陣；αT為熱膨脹張量；Tp為瓦塊溫度分布;Tref為固體材料熱膨脹的參考溫度。

1.6 邊界條件

熱邊界條件簡述如下：

(1)在流體和瓦塊的接觸面上，假設(shè)熱通量是連續(xù)的，并且流體和瓦塊接觸面的溫度是相同的。

(2)在流體與軸的接觸面上，假設(shè)軸表面溫度在圓周方向上是相等的，流體與軸之間的整體熱交換為0。

(3)在流體入口邊界處，溫度設(shè)定為入口混合溫度。根據(jù)MITSUI等[15]的研究，使用SUH和PALAZZOLO[16]提出的改進(jìn)入口混合溫度模型來計算瓦塊入口的潤滑油混合溫度，并假設(shè)熱油攜帶系數(shù)為0.8。

(4)在流體出口邊界和軸向兩端邊界，假設(shè)熱通量為0。

(5)在其他瓦塊表面，假設(shè)與周圍環(huán)境進(jìn)行熱對流。如文獻(xiàn)[3]所述，背面的熱對流系數(shù)為350 W/(m2·K)，其他表面的熱對流系數(shù)為115 W/(m2·K)。

2 可傾瓦滑動軸承熱彈流體動力潤滑模型的計算方法

對于某一瓦塊內(nèi)壓力場、溫度場和瓦塊變形場的計算通過COMSOL和MATLAB聯(lián)合仿真實(shí)現(xiàn)。利用COMSOL分別建立3個獨(dú)立的有限元模型，分別為壓力子模型、溫度子模型和變形子模型。壓力子模型采用COMSOL中偏微分方程物理場模塊進(jìn)行構(gòu)建，用于求解考慮湍流的修正雷諾方程，從而獲得油膜壓力。溫度子模型也采用COMSOL中偏微分方程物理場模塊進(jìn)行構(gòu)建，用于求解能量方程和瓦塊熱傳導(dǎo)方程，從而獲得流體和瓦塊溫度場分布。而變形子模型采用COMSOL中熱應(yīng)力物理場模塊進(jìn)行構(gòu)建，用于獲得瓦塊變形分布。COMSOL與MATLAB之間的數(shù)據(jù)傳遞和交換通過COMSOL的livelink for matlab函數(shù)實(shí)現(xiàn)[17-18]。在MATLAB中計算偏微分方程的系數(shù)，將其寫入文本文件，然后在COMSOL中通過插值函數(shù)讀取。利用COMSOL計算出的油膜壓力、瓦塊溫度和瓦塊變形量，通過mpheval、mphinterp、mphglobal等函數(shù)輸出到MATLAB程序中。此外，COMSOL子模型和MATLAB程序中的網(wǎng)格劃分是不同的。在MATLAB中，網(wǎng)格為16×16×10。而在COMSOL中，3種有限元模型有各自的網(wǎng)格劃分，如圖2所示。此外，壓力、溫度和瓦塊變形的收斂準(zhǔn)則采用相對誤差為10-3。而力矩平衡和負(fù)載平衡的相對誤差為10-2。

圖2 COMSOL中構(gòu)建的子模型的網(wǎng)格劃分

圖3所示為考慮湍流的可傾瓦滑動軸承熱流體動力潤滑模型總體計算流程。計算開始前，先假定軸心位置和各個瓦塊的傾斜角。各個瓦塊的力矩平衡通過牛頓-拉斐遜法實(shí)現(xiàn)，從而得到各個瓦塊的傾斜角。對于無負(fù)載的上瓦塊可能不存在力矩平衡點(diǎn)。如果沒有平衡位置，則給出一個臨時傾斜角度(該傾斜角可使力矩的絕對值降至最小)。在給定負(fù)載下，軸心位置同樣采用牛頓-拉斐遜法迭代獲得，直到流體承載力和負(fù)載平衡。表1所示為可傾瓦徑向滑動軸承結(jié)構(gòu)參數(shù)，表2所示為可傾瓦滑動軸承材料參數(shù)，表3所示為潤滑油參數(shù)。

圖3 考慮湍流的可傾瓦滑動軸承熱彈流體動力潤滑模型總體計算流程

表1 可傾瓦徑向滑動軸承結(jié)構(gòu)參數(shù)

表2 可傾瓦滑動軸承材料參數(shù)

表3 HP-8A潤滑油參數(shù)

3 結(jié)果與分析

為了驗證所建立的模型，在忽略瓦塊變形的情況下，將文中模型預(yù)測結(jié)果與TANIGUCHI等[3]的預(yù)測結(jié)果進(jìn)行比較。仿真中，轉(zhuǎn)速為3 000 r/min，載荷為180 kN，施加在2個瓦塊之間。采用直徑479 mm、長度300 mm的四瓦可傾瓦滑動軸承。瓦塊弧長80°，徑向間隙0.612 mm，瓦塊厚度121 mm。支承點(diǎn)偏移量50%，預(yù)緊為0。潤滑油型號為ISO VG32。圖4所示為文中模型和TANIGUCHI等[3]模型預(yù)測的油膜壓力、瓦塊表面溫度和油膜厚度結(jié)果比較?？梢杂^察到，文中的預(yù)測結(jié)果與TANIGUCHI等[3]的仿真結(jié)果吻合較好。

圖4 文中模型結(jié)果和 TANIGUCHI等[3]仿真結(jié)果的對比

與以往自主編程實(shí)現(xiàn)的方法相比，采用COMSOL與MATLAB聯(lián)合仿真的方法可以直觀、方便地顯示仿真結(jié)果，這有助于理解可傾瓦徑向滑動軸承運(yùn)行中的物理現(xiàn)象。以下瓦塊2為例，圖5所示為可傾瓦滑動軸承下瓦塊2 的壓力場、溫度場和變形場的三維圖。其中，載荷為4 927 N，轉(zhuǎn)速為6 000 r/min。紅色箭頭所示為流體流動方向。如圖1所示，載荷施加在下瓦塊2和3之間,載荷由下瓦塊2和3承受。下瓦塊2表面的油膜壓力分布如圖5(a)所示。由于摩擦產(chǎn)生熱量，出口油溫高于進(jìn)口油溫，如圖5(b)和5(c)所示。靠近瓦塊中部的徑向變形增大了油膜厚度，而邊緣處的徑向變形對油膜厚度的影響相反，如圖5(d)所示。

圖5 可傾瓦滑動軸承瓦塊2 的壓力場、溫度場和變形場的三維圖

對于文中研究的四瓦可傾瓦滑動軸承，由于載荷施加在兩塊下瓦之間，載荷由兩塊下瓦承受，下瓦表面油膜壓力大，而上瓦表面的油膜壓力很小。并且，下瓦的瓦溫遠(yuǎn)高于上瓦的瓦溫。因此，下文只分析軸承參數(shù)變化對下瓦的影響。仿真中，負(fù)載為20 000 N，轉(zhuǎn)速為10 000 r/min。為了研究軸承結(jié)構(gòu)參數(shù)對可傾瓦滑動軸承性能的影響，將采用控制變量法進(jìn)行研究。軸承參數(shù)包括瓦塊弧長、徑向間隙、支承點(diǎn)偏移量和軸向?qū)挾取?/p>

3.1 瓦塊弧長對可傾瓦滑動軸承性能的影響

圖6所示為瓦塊弧長對可傾瓦滑動軸承性能的影響。隨著瓦塊弧長的增加，承載面積變大，面壓(F/(L×D))降低，從而油膜最大壓力會降低，如圖6(a)所示。同時，由于面壓降低，偏心率降低，油膜厚度增加；進(jìn)一步地，由于油膜增厚，油膜局部剪切程度降低，瓦面最大溫度降低，如圖6(b)所示。圖6(c)所示為油膜厚度隨瓦塊弧長的變化趨勢，可見油膜厚度隨瓦塊弧長的增加而增加。而由于偏心率降低，軸心位置會上浮，如圖6(d)所示。

圖6 瓦塊弧長對可傾瓦滑動軸承性能的影響

表4所示為瓦塊弧長對可傾瓦滑動軸承摩擦功耗的影響。

表4 瓦塊弧長對可傾瓦滑動軸承摩擦功耗的影響

由表4可以看出，隨著瓦塊弧長的增加，承載面積增加，雖然由于膜厚增加，油膜局部剪切變?nèi)酰蓛A瓦滑動軸承的摩擦功耗也會增加。此外，需要注意的是，由于載荷施加在2個下瓦之間，兩塊下瓦上的壓力、溫度、油膜厚度分布基本相同。

3.2 徑向間隙對可傾瓦滑動軸承性能的影響

圖7所示為徑向間隙對可傾瓦滑動軸承性能的影響。隨著徑向間隙的增大，由于負(fù)載不變，油膜壓力變化很小，如圖7(a)所示。此外，由于徑向間隙增大，導(dǎo)致偏心率增加，油膜厚度增加，油膜局部剪切變?nèi)?，瓦面最大溫度降低，如圖7(b)所示。油膜厚度隨徑向間隙的變化趨勢如圖7(c)所示，可見油膜厚度隨徑向間隙的增加而增加。由于偏心率的增加，軸心位置相應(yīng)下沉，如圖7(d)所示。表5所示為徑向間隙對可傾瓦滑動軸承摩擦功耗的影響。隨著軸向間隙的增加，油膜厚度增加，油膜剪切應(yīng)力降低，摩擦功耗會略有降低。

圖7 徑向間隙對可傾瓦滑動軸承性能的影響

表5 徑向間隙對可傾瓦滑動軸承摩擦功耗的影響

3.3 支承點(diǎn)偏移量對軸承性能的影響

圖8所示為支承點(diǎn)偏移量對可傾瓦滑動軸承性能的影響。從圖8(a)可以看出，隨著支承點(diǎn)偏移量的增加，瓦塊2和瓦塊3的油膜壓力分布差異性越來越大。隨著支承點(diǎn)偏移量的增加，瓦塊2的油膜壓力增大，而瓦塊3的油膜壓力降低，如圖8(a)所示。與此同時，瓦塊2的油膜厚度減小，油膜局部剪切變強(qiáng)，瓦塊表面最大溫度增大，相反，瓦塊3的油膜厚度減小，油膜局部剪切變?nèi)?，瓦塊表面最大溫度降低，進(jìn)而瓦塊2的表面溫度大于瓦塊3的表面溫度，如圖8(b)所示。隨著支承點(diǎn)偏移量的增大，瓦塊2 和瓦塊3的油膜厚度的變化趨勢相反，如圖8(c)所示。而對于軸心位置有向右上方移動的趨勢，如圖8(d)所示。表6所示為支承點(diǎn)偏移量對可傾瓦滑動軸承摩擦功耗的影響?？梢钥闯?，隨著支承點(diǎn)偏移量的增加，承載能力會有所提高，可傾瓦滑動軸承的摩擦功耗略微增加。

圖8 支承點(diǎn)偏移量對可傾瓦滑動軸承性能的影響

表6 支承點(diǎn)偏移量對可傾瓦滑動軸承摩擦功耗的影響

3.4 軸向?qū)挾葘S承性能的影響

圖9所示為軸向?qū)挾葘蓛A瓦滑動軸承性能的影響。可以看出，隨著瓦塊軸向?qū)挾鹊脑黾?，承載面積增加，承載能力增大,而面壓降低，油膜最大壓力降低，如圖9(a)所示。面壓降低會導(dǎo)致偏心率減小，從而油膜變厚。由于油膜增厚，油膜局部剪切程度降低，瓦面最大溫度會相應(yīng)降低，如圖9(b)所示。油膜厚度隨著軸向?qū)挾鹊淖兓厔萑鐖D9(c)所示，可見油膜厚度隨著軸向?qū)挾鹊脑龃蠖鴾p小。隨著偏心率降低，軸心位置有上浮的趨勢，如圖9(d)所示。表7所示為軸向?qū)挾葘蓛A瓦滑動軸承摩擦功耗的影響。一方面，隨著軸向?qū)挾鹊脑黾?，接觸面積增加；另一方面，由于膜厚增加，油膜局部剪切變?nèi)?。相比之下，接觸面積增大的影響要比油膜厚度局部剪切變?nèi)醯挠绊戯@著，因此，可傾瓦滑動軸承的摩擦功耗會大幅增加。

圖9 軸向?qū)挾葘蓛A瓦滑動軸承性能的影響

表7 軸向?qū)挾葘蓛A瓦滑動軸承摩擦功耗的影響

4 結(jié)論

結(jié)合商用軟件和自主編程的優(yōu)點(diǎn)，采用COMSOL和MATLAB聯(lián)合仿真的方法構(gòu)建了考慮湍流的可傾瓦滑動軸承熱彈流體動力潤滑模型。通過與文獻(xiàn)對比，驗證了模型的準(zhǔn)確性?；谒鶚?gòu)建模型，研究了軸承參數(shù)包括瓦塊弧長、徑向間隙、支承點(diǎn)偏移量和寬度對可傾瓦滑動軸承性能的影響。主要結(jié)論如下：

(1)隨著瓦塊弧長的增加，油膜壓力降低，油膜厚度增加，瓦面溫度降低，而軸心位置略有上浮，可傾瓦滑動軸承的摩擦功耗也會增加。

(2)隨著徑向間隙的增大，油膜壓力變化很小，瓦面溫度降低，油膜厚度增加，軸心位置相應(yīng)下沉，摩擦功耗會略有降低。

(3)隨著支承點(diǎn)偏移量的增加，瓦塊2的油膜壓力、溫度增大，而油膜厚度降低；而瓦塊3的油膜壓力，溫度降低，而油膜厚度增加。軸心位置有向右上方移動的趨勢，承載能力、摩擦功耗略微增加。

(4)隨著瓦塊軸向?qū)挾鹊脑黾?，油膜壓力降低，瓦面溫度降低，油膜厚度增加，軸心位置有上浮的趨勢，摩擦功耗大幅增加。

對于可傾瓦滑動軸承，瓦面最大溫度和摩擦功耗是2個重要的指標(biāo)。瓦面溫度關(guān)系到軸承安全運(yùn)行，而摩擦功耗會影響節(jié)能減排。因此，需要對軸承參數(shù)進(jìn)行優(yōu)化選擇。文中的研究為可傾瓦滑動軸承的優(yōu)化設(shè)計提供有價值的參考。