船舶可傾瓦推力軸承潤滑油膜的軸向動特性計算方法

張贛波，趙 耀，儲 煒，袁 華

（華中科技大學 船舶與海洋工程學院，武漢 430074）

船舶可傾瓦推力軸承潤滑油膜的軸向動特性計算方法

張贛波，趙 耀，儲 煒，袁 華

（華中科技大學 船舶與海洋工程學院，武漢 430074）

推力環和推力瓦之間的楔形潤滑油膜是實現螺旋槳推力傳遞的重要環節，其軸向動特性直接關乎船舶軸系轉子的縱向振動特性。文章分別論述了船舶可傾瓦推力軸承楔形潤滑油膜軸向動特性的一維流近似解析方法和二維流數值方法，在已求得油膜靜特性基礎上，分別結合偏導數法和小攝動法獲解了油膜動特性，推導了兩種方法計算油膜動特性的求解式，并給出了詳細計算過程。以某船舶可傾瓦推力軸承為算例，對兩種方法的計算結果作了比較分析，得到了最小油膜厚度、油膜承載力和油膜軸向剛度三者間的一般變化規律，討論了油膜動特性隨軸轉速的變化關系。計算結果為基于轉子動態模型研究軸系縱向振動的傳遞機理提供了油膜動特性數據。

可傾瓦；推力軸承；動壓潤滑；油膜動特性

Key words:tilting pad;thrust bearing;hydrodynamic lubrication;oil film dynamic characteristics

0 引 言

實船上，米歇爾式（Michell）和金斯伯雷式（Kingsbury）兩種類型的可傾瓦滑動式推力軸承由于具有結構緊湊、體積小、重量輕、摩擦系數小以及承載力大等優點獲得了廣泛應用。兩者的結構差異主要在于推力瓦的支撐結構形式，前者為單點球面支柱剛性支撐，后者為平衡塊交錯疊加剛性支撐[1-2]。推力軸承是船舶推進系統中的最重要傳力元件，其運行的安全穩定性直接關系到船舶的航行性能，故推力軸承設計要求有高度的可靠性，除必須與螺旋槳特性相匹配外，還需對其潤滑性能和動特性進行分析和準確預測[3]。

推力軸承還作為螺旋槳脈動推力誘導的主推進軸系縱向振動向船體傳遞的重要通道，其動力學參數與主推進軸系縱向振動特性密切相關[4]。從轉子動力學角度看，推力環與推力瓦之間的楔形潤滑油膜是軸系轉子連接軸承基座不可或缺的環節，但限于油膜自身的復雜性，以往對其軸向動特性的定量研究工作涉及較少。國內多數學者在研究船舶主推進軸系縱向振動特性時并未將油膜動特性納入軸系動力學模型中，關于潤滑油膜對軸系縱向振動特性的影響也始終未形成定論[5]。潤滑油膜作為推力軸承實現螺旋槳推力傳遞的重要組成部分，明確其動特性可完善對推力軸承的動力學建模。

Schwanecke[6]和Vassilopoulos[7]較早開展推力軸承潤滑油膜軸向剛度和阻尼的理論計算，其理論基礎是簡化的無限寬軸承潤滑理論，由于忽略油膜壓力沿徑向的二次分布，且未考慮油膜的熱效應，其分析方法還有待繼續深入。潤滑油膜動特性取決于推力瓦幾何結構參數、潤滑油物性參數和軸轉速等因素，并與這些因素是一種高度復雜的非線性關系，其求解過程較為繁復，需遵循先靜態后動態的計算順序[8-9]。

本文先簡要介紹了油膜軸向剛度的一維流近似解析方法，將扇形推力瓦近似為無限寬矩形滑塊，獲得了油膜承載力的解析式，求偏導得到了油膜軸向剛度；再基于二維流動壓潤滑理論論述了較為準確的油膜動特性計算方法，推導了油膜軸向剛度和阻尼的求解式，運用中差分法對動壓潤滑偏微分方程進行了離散數值求解，給出了詳細計算流程，討論了潤滑油膜動特性隨軸轉速的變化關系，為后續主推進軸系轉子系統軸向振動的傳遞機理分析作鋪墊。

1 一維流近似解析方法

圖1 可傾瓦推力軸承結構示意圖Fig.1 Schematic diagram of tilting pad thrust bearing

可傾瓦推力軸承內部結構示意圖如圖1所示，推力環前后各有一排推力瓦，分別用于承受螺旋槳正倒車推力。圖中各符號意義為：O為推力瓦幾何中心，ω為軸轉動角速度，φ為推力瓦圓心角，ri、ro分別為推力瓦的內外半徑，b=ro-ri為推力瓦徑向寬度，rp、θp分別為推力瓦支撐中心Op的徑向和周向坐標，h2，h1分別為進出油口邊油膜平均厚度。考慮推力瓦沿推力環分布的圓周對稱性，各推力瓦與推力環之間楔形間隙內的油膜具有完全相同的幾何形狀，故只需討論一個楔形油膜的潤滑特性，即可獲得推力軸承潤滑油膜的總特性。

早期限于求解技術的發展，潤滑油膜動特性（主要為剛度）的計算多基于簡化的一維流動壓潤滑理論。一維流理論將實際扇形推力瓦近似簡化為平面矩形滑塊，只考慮油膜壓力沿周向的二次分布，實際上相當于徑向無限寬推力瓦，這樣可將二維雷諾方程降階為一維雷諾方程。平面矩形滑塊油膜承載力的求解式為[10]：

式中：μ為軸承熱平衡后的潤滑油動力粘度；l=（ri+ro）φ/2為推力瓦平均半徑處的周向長度；v=（ri+ro）ω/ 2為推力瓦平均半徑處的線速度；a=h2/h1為推力瓦傾斜度。

由力矩平衡條件，可推導出推力瓦傾斜度a和推力瓦支撐中心與進油口邊的距離L的關系式為：

船舶可傾瓦推力軸承推力瓦周向偏心距e的取值一般在0.05～0.1l范圍內[11]，通常取為e=0.08l，則L=0.58l，代入（2）式可倒推出推力瓦傾斜度a。結合油膜承載力與螺旋槳靜推力平衡條件，將已求解的推力瓦傾斜度a代入（1）式可計算出油口邊油膜厚度h1。由此可見，一旦推力瓦支撐中心取定后，隨著螺旋槳推力的變化，油膜厚度h1和推力瓦傾斜度a會自行調整。這正是可傾瓦推力軸承的結構特點。由于（2）式是超越方程，無法得到推力瓦傾斜度a的解析解，可通過一組a和L的離散值由Lagrange插值公式間接獲解。a和L的對應值列于表1，當取L=0.58l時，插值可得推力瓦傾斜度a= 2.247 8。

為確定熱平衡后的潤滑油動力粘度，還需求解潤滑油膜摩擦力和流量，兩者求解式分別為[10]：

由絕熱流動熱平衡方程可得潤滑油溫升為：

式中：c、ρ分別為潤滑油比熱容和密度。

設潤滑油進油口邊溫度為t0，則熱平衡后的潤滑油平均溫度為：

式中：α為考慮實際軸承傳導散熱的影響予以修正的系數。船舶推力軸承屬于中低速軸承，可取α= 0.8。

獲解潤滑油熱平衡溫度后，可由粘溫關系獲得對應溫度的潤滑油動力粘度。這里借用美國材料試驗協會（ASTM）推薦的粘溫指數公式：

式中：系數A和B可由潤滑油兩個已知溫度下的粘度值確定。

將油膜承載力對油膜厚度求偏導，并考慮推力瓦支撐偏心的影響，得到油膜軸向剛度為：

式中：Z為推力瓦數，μ和h1可結合（1）～（6）式獲解。

一維流理論未考慮油膜壓力沿徑向的變化，計算油膜壓力呈半圓柱面分布，在同樣油膜厚度情形下，一維流計算油膜承載力必定大于實際軸承外荷載。當推力瓦寬長比滿足b/l＞3時，一維流計算油膜承載力同實際值較為接近[12]，而實際船舶可傾瓦推力軸承推力瓦的寬長比一般為0.8～1，故按（7）式計算潤滑油膜軸向剛度只能是近似的。

2 二維流數值方法

鑒于一維流理論的不足，且無法獲解油膜阻尼，還需借助于二維流動壓潤滑理論，以實際有限寬扇形推力瓦為潤滑油膜載體。

2.1 潤滑油膜控制方程

（1）雷諾方程

雷諾方程可看作Navier-Stokes方程的特殊形式。對于可傾瓦推力軸承，建立如下柱坐標形式的穩態雷諾方程[10]：

（2）能量方程

潤滑油粘性使其在運動過程中摩擦生熱，導致進入楔形間隙內的油膜溫度升高。流體動壓潤滑以對流散熱為主要熱傳導形式，可近似為絕熱流動，即摩擦熱能全部由潤滑油流動帶走。由功、能守恒原理可推導出潤滑油膜的絕熱流動能量方程為[10]：

（3）膜厚方程

軸系穩定旋轉過程中，推力瓦繞支撐點發生微小傾斜，可近似為線支撐，如圖2所示。

圖2 推力瓦傾斜示意圖Fig.2 Tilting sketch of pad

由幾何關系可推導出油膜厚度方程為：

式中：hp為推力瓦支撐中心位置對應的油膜厚度，γp為推力瓦繞支撐線OP的傾斜角。一般地，γp很小，有tanγp≈γp。

設最小油膜厚度位置坐標為 （rm，θm），代入（10）式可得最小油膜厚度為：

對（11）式移項整理，代入（10）式，可得用最小油膜厚度表示的油膜厚度方程為：

雷諾方程和能量方程都是偏微分方程，且兩者相互耦合，難以獲得解析解，只能借助于數值計算方法。為改善數值計算的穩定性，并使分析結果具有普遍性，有必要對上述控制方程進行無量綱化處理。定義如下無量綱量：

式中：μ0為對應進油口邊溫度t0的潤滑油動力粘度。

結合上述無量綱量，直接推導出潤滑油膜無量綱控制方程分別為：

潤滑油膜靜特性包括油膜承載力和力矩等，其中，油膜力矩是油膜分布壓力對推力瓦支撐中心的合力矩。兩者將作為油膜壓力場數值計算的控制條件，表達式分別為：

2.2 邊界條件和控制條件

潤滑油膜的邊界條件包括壓力、溫度和粘度。其中，壓力邊界條件針對推力瓦四周邊，溫度和粘度邊界條件只針對進油口邊。四周邊壓力邊界條件為：

進油口邊溫度和粘度邊界條件為：

潤滑油膜的控制條件是油膜承載力等于螺旋槳靜推力以及油膜力矩等于零，即：

式中：KT為推力系數，一般由螺旋槳敞水性能試驗獲得，ρw為水密度，n為螺旋槳轉速，D為螺旋槳直徑。

2.3 數值計算

以單個推力瓦與推力環之間的潤滑油膜為求解域，將潤滑油膜劃分為有限網格，如圖3所示。將油膜沿周向劃分為m格，分別編號1，L，m+1，周向步長為Δθ=φ/m；沿徑向劃分為n格，分別編號1，L，n+ 1，徑向步長為ΔR=1/n。

圖3 潤滑油膜有限差分網格Fig.3 Finite difference grid of oil film

這里直接給出無量綱雷諾方程和能量方程的中心差分格式：

式中：各系數分別為：

二維流變粘度潤滑分析需聯立求解雷諾方程、能量方程、粘溫方程和膜厚方程，其數值計算流程為：首先以潤滑油入口溫度所確定的粘度值作為潤滑油粘度的初值，求解雷諾方程獲得等粘度條件下的油膜壓力場，由所求油膜壓力場求解能量方程獲得油膜溫度場，檢查壓力場和溫度場是否滿足給定收斂條件，若不滿足，按粘溫方程獲解新的油膜粘度，再次計算雷諾方程和能量方程，直至壓力場和溫度場滿足收斂條件；得到收斂的壓力場后，繼續檢查是否滿足控制條件，如滿足，結束計算，否則修改相應變量，重新計算。數值計算流程框圖如圖4所示。內層迭代循環是油膜壓力和溫度，外層迭代循環是油膜承載力和力矩，修改變量是最小油膜厚度和推力瓦傾角。

編制計算程序時，為避免迭代發散或者陷入無限循環，必須采取一些措施[12]，如設置無量綱計算變量（P、T、H）的上下限(需根據實際推力軸承的潤滑性能選取)；引入松弛因子，分別運用超松弛迭代法和低松弛迭代法計算油膜壓力場和溫度場；設置內外層循環的最大迭代次數等。

2.4 潤滑油膜動特性

完成靜平衡位置的油膜壓力場和溫度場的求解是分析油膜動特性的前提。油膜軸向動特性實際上反映油膜承載力的相應變化，故必須以非定常雷諾方程作為分析基礎。在無量綱穩態雷諾方程中加入擠壓油膜項，表達式為：

圖4 二維流動壓潤滑數值計算流程框圖Fig.4 Numerical calculation flow diagram of 2-D hydrodynamic lubrication

推力環在螺旋槳脈動推力作用下在其靜平衡位置附近作微幅軸向振動擠壓油膜，引起油膜厚度的擾動。油膜厚度可表示為：

式中：H0是靜平衡位置油膜厚度，Xk和ξk分別是各簡諧次軸向振動位移幅值和相位。

由（20）式可推導如下關系式成立：

油膜厚度的擾動引起油膜壓力的擾動。小擾動條件下，油膜壓力可近似表示為推力環偏離靜平衡位置的軸向位移和速度的線性函數，此時用一個剛度和阻尼系數表示油膜的動特性。將油膜動壓力在靜平衡位置展開為軸向位移ΔH和速度Δ的一階Taylor級數為：

式中：P0是靜平衡位置油膜壓力分布值，下標0表示在靜平衡位置取導數。

將油膜動壓力沿推力瓦面積分求得油膜動承載力為：

式中：Ft0是靜平衡位置油膜承載力，對應于螺旋槳靜推力。

由（23）式可直接得到油膜無量綱軸向剛度和阻尼求解式為：

由于無法得到油膜壓力的顯式表達式，由（24）式直接計算油膜動特性幾乎不可能。為此，可將（20）式和（22）式代入（19）式，結合（21）式，比較方程兩端對應項，同時略去高階小量，整理得：

（25）式左端項與穩態雷諾方程式（13）形式完全相同，只是右端項稍有不同。可運用同計算雷諾方程的有限差分法對方程進行超松弛迭代求解得到，再代入（24）式，由Simpson數值積分法獲得油膜的無量綱軸向剛度Ko和阻尼Co，量綱化后即可得實際潤滑油膜的動特性。結合已定義的無量綱量，容易推得油膜量綱動特性與無量綱動特性的轉化關系式為：

從（26）式看出，推力瓦徑向寬度、油膜初始動力粘度、最小油膜厚度和軸轉速是決定油膜動特性的關鍵因素。

3 算例分析

某可傾瓦動壓潤滑滑動推力軸承的推力瓦結構參數如表2所示。所用潤滑油牌號為ISO VG46，密度為890 kg/m3，比熱容為1 922.8 J/kg°C，在20°C和40°C時的運動粘度分別為105 mm2/s、46 mm2/ s，進油口溫度為35°C，對應的動力粘度為0.051 Pa·s。

表2 推力瓦結構參數Tab.2 Parameters of pad

可傾瓦推力軸承的結構優勢在于軸旋轉過程中會自行形成一個收斂的楔形油膜和足夠大的最小油膜厚度。由于螺旋槳推力是軸轉速的函數，故不同軸轉速對應的油膜動特性必然存在差異。表3列出典型轉速下油膜軸向剛度的一維流和二維流計算結果。在分析過程中，螺旋槳靜推力取為軸轉速的二次方關系。比較看出，油膜軸向剛度一維流計算值均小于二維流計算值，分析有兩方面原因：二維流考慮潤滑油膜的徑向端泄，計算油膜壓力呈拋物面分布，與實際油膜壓力場較為接近，在相同螺旋槳推力作用下，二維流油膜厚度必定小于一維流；一維流計算油流量偏小，油膜溫升大，潤滑油粘度降低。潤滑油膜一維流模型相對二維流模型較為粗糙，但計算量顯著減小。從算例結果看，油膜軸向剛度二維流計算值約為一維流計算值的2倍，但更準確的一維流模型修正值還需要更多算例分析予以獲得。

圖5給出二維流變粘度計算的油膜軸向剛度和阻尼隨軸轉速的變化關系。從圖5看出，油膜無量綱軸向剛度和阻尼都在一個平均值附近波動，且波動幅值很微小，表明油膜無量綱特性受軸轉速影響較小，這一規律有助于對油膜動特性的預估.實際油膜軸向剛度和阻尼都隨軸轉速增加遞增，且兩者數值均較大。對于獨立安裝的船用推力軸承，其縱向剛度范圍為1×109～5×109N/m，比較可知，在軸中高轉速行程內，油膜軸向剛度已超過推力軸承金屬實體剛度。

表3 油膜剛度一維流和二維流計算值比較Tab.3 Comparison of oil film stiffness with 1-D and 2-D

圖5 油膜軸向剛度和阻尼隨軸轉速的變化關系Fig.5 Oil film stiffness and damping versus shafting rotating speed

油膜承載力決定于最小油膜厚度，而油膜軸向剛度又與最小油膜厚度相關，三者間的相互關系可用如圖6所示的推力軸承動壓特性曲線表示。從圖6可見，最小油膜厚度與油膜承載力和油膜剛度呈負增長關系。隨軸轉速增加，油膜厚度迅速減薄，相應的油膜承載力和剛度迅速增加以適應增大的螺旋槳推力。

圖6 推力軸承動壓特性曲線Fig.6 Dynamic characteristic curve of thrust bearing

4 結 論

本文探討了可傾瓦推力軸承楔形潤滑油膜的軸向動特性，詳細給出了兩種方法的計算流程，結合某型推力軸承算例，得出了以下一些結論：

（1）油膜動特性分析需依據流體動壓潤滑理論，按對推力瓦模型的簡化處理，分為一維流和二維流方法。一維流理論是近似方法，對油膜剛度的計算值相對偏低，二維流計算方法較為準確，但計算量顯著增加；

（2）最小油膜厚度和軸轉速是決定油膜動特性的最關鍵因素，油膜承載力和油膜軸向剛度與最小油膜厚度呈此消彼長關系；

（3）受螺旋槳脈動推力對潤滑油膜的擠壓效應，隨軸轉速增加，油膜軸向剛度和阻尼都逐漸增大，且兩者量級均很高。在一定軸轉速下，油膜剛度值超過推力軸承金屬實體剛度。

[1]CB3103-81.船舶推進軸系滑動推力軸承[S].1982.

[2]胡榮華.船用滑動推力軸承結構設計研究[J].船舶工程,2007,29(5):60-64. Hu Ronghua.Structural design research of the marine thrust bearing[J].Marine Engineering,2007,29(5):60-64.

[3]Sternlicht D B,Reid J C,Arwas E B.Review of propeller shaft thrust bearings[J].A.S.N.E.Journal,1959,(5):277-289.

[4]張贛波,趙 耀.船舶主推進軸系縱向振動的彈性波解析研究[J].中國造船,2012,53(3):140-150. Zhang Ganbo,Zhao Yao.Elastic wave analysis on longitudinal vibration of marine propulsion shafting[J].Shipbuilding of China,2012,53(3):140-150.

[5]趙 耀,張贛波,李良偉.船舶推進軸系縱向振動及其控制技術研究進展[J].中國造船,2011,52(198):259-269. Zhao Yao,Zhang Ganbo,Li Liangwei.Review of advances on longitudinal vibration of ship propulsion shafting and its control technology[J].Shipbuilding of China,2011,52(198):259-269.

[6]Schwanecke H.Investigations on the hydrodynamic stiffness and damping of thrust bearing in ship[J].Transactions of the Institute of Marine Engineers,1979,(91):68-77.

[7]Vassilopoulos M L.Methods for computing stiffness and damping properties of main propulsion thrust bearing[J].International Shipbuilding Progress,1982,329(29):13-31.

[8]陳 渭.流體動壓潤滑推力軸承動特性及其對轉子橫向振動狀態影響的研究[D].西安:西安交通大學,1991. Chen wei.The dynamic performances of hydrodynamic lubricated thrust bearing and the influence on crosswise vibration state of rotor[D].Xi’an:Xi’an Jiaotong University,1991.

[9]李 忠,袁小陽,朱 均.可傾瓦推力軸承的線性和非線性動特性研究[J].中國機械工程,2000,11(5):560-563. Li Zhong,Yuan Xiaoyang,Zhu Jun.A study of linear and non-linear dynamic performance for tilting-pad thrust bearing [J].China Mechanical Engineering,2000,11(5):560-563.

[10]溫詩鑄,黃 平.摩擦學原理[M].(第二版)北京:清華大學出版社,2002.

[11]商圣義.民用船舶動力裝置[M].北京:人民交通出版社,1996.

[12]楊沛然.流體潤滑數值分析[M].北京:國防工業出版社,1998.

Calculation method for axial dynamic characteristics of lubricant oil film in marine tilting pad thrust bearing

ZHANG Gan-bo,ZHAO Yao,CHU Wei,YUAN Hua
(School of Naval Architecture&Ocean Engineering,Huazhong University of Science&Technology,Wuhan 430074,China)

The oil film between thrust collar and tilting pads is one of significant transmitters for propeller thrust transmission.The axial dynamic characteristics of oil film make great influence on longitudinal vibration of propulsion shafting.In this paper,two methods are presented to evaluate the longitudinal stiffness and damping of oil film,namely one-dimension approximate analytical method and two-dimension numerical method.For the first method,the pad is interpreted as rectangle slider.While for the second method,the actual sector pad is directly used as objective.As the premise of dynamic analysis,the static characteristics of oil film with action of propeller steady thrust are firstly carried out based on dynamic lubrication principles.And then the calculations of oil film dynamic characteristics are derived in great detail by the partial derivative method and the small perturbation method.Taking the example of a marine thrust bearing, the comparison of results with two proposed methods is conducted.Further,the relationship among minimum film thickness,film carrying capacity and film stiffness is analyzed,and so also is for the rule between oil film dynamic characteristics and shafting rotating speed.The values of oil film dynamic characteristics can be applied for transmission research of propulsion shafting longitudinal vibration.

U664.3

：Adoi:10.3969/j.issn.1007-7294.2017.05.011

1007-7294（2017）05-0603-10

2016-12-07

教育部高等學校博士學科點專項科研基金項目（20130142110014）；國家自然科學基金委員會資助項目（51479078）

張贛波（1987-），男，博士研究生，E-mail：hustzgb@126.com；趙 耀（1958-），男，博士，教授，博士生導師，E-mail：yzhaozzz@hust.edu.cn。