胡 楊 孟永鋼

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

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

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

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

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

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

1 考慮湍流的可傾瓦滑動(dòng)軸承熱彈流體動(dòng)力潤(rùn)滑模型

1.1 油膜厚度描述

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

圖1 可傾瓦滑動(dòng)軸承示意與坐標(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)的線性湍流潤(rùn)滑理論。修正雷諾方程公式表達(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ù)，具體計(jì)算可參考文獻(xiàn)[3]。

1.3 三維能量方程

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

(7)

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

(8)

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

(9a)

(9b)

(9c)

1.4 熱傳導(dǎo)方程

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

(10)

1.5 瓦塊變形

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

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

(11)

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

1.6 邊界條件

熱邊界條件簡(jiǎn)述如下：

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

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

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

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

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

2 可傾瓦滑動(dòng)軸承熱彈流體動(dòng)力潤(rùn)滑模型的計(jì)算方法

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

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

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

圖3 考慮湍流的可傾瓦滑動(dòng)軸承熱彈流體動(dòng)力潤(rùn)滑模型總體計(jì)算流程

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

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

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

3 結(jié)果與分析

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

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

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

圖5 可傾瓦滑動(dòng)軸承瓦塊2 的壓力場(chǎng)、溫度場(chǎng)和變形場(chǎng)的三維圖

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

3.1 瓦塊弧長(zhǎng)對(duì)可傾瓦滑動(dòng)軸承性能的影響

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

圖6 瓦塊弧長(zhǎng)對(duì)可傾瓦滑動(dòng)軸承性能的影響

表4所示為瓦塊弧長(zhǎng)對(duì)可傾瓦滑動(dòng)軸承摩擦功耗的影響。

表4 瓦塊弧長(zhǎng)對(duì)可傾瓦滑動(dòng)軸承摩擦功耗的影響

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

3.2 徑向間隙對(duì)可傾瓦滑動(dòng)軸承性能的影響

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

圖7 徑向間隙對(duì)可傾瓦滑動(dòng)軸承性能的影響

表5 徑向間隙對(duì)可傾瓦滑動(dòng)軸承摩擦功耗的影響

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

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

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

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

3.4 軸向?qū)挾葘?duì)軸承性能的影響

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

圖9 軸向?qū)挾葘?duì)可傾瓦滑動(dòng)軸承性能的影響

表7 軸向?qū)挾葘?duì)可傾瓦滑動(dòng)軸承摩擦功耗的影響

4 結(jié)論

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

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

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

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

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

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