溫小飛，周瑞平，袁 強(qiáng)，雷俊松

（1.武漢理工大學(xué)能源與動力工程學(xué)院，湖北武漢 430063；2.浙江海洋大學(xué)港航與交通運(yùn)輸工程學(xué)院，浙江舟山 316022）

船舶軸系軸承在運(yùn)行工況常因軸承負(fù)荷變化而導(dǎo)致漏油[1]、軸承發(fā)熱等問題，但一直以來相關(guān)船舶規(guī)范文件[2-3]中均對運(yùn)行工況的軸承負(fù)荷做出明確要求或限制。在相關(guān)研究方面，1992年徐龍祥等[4]考慮了軸承油膜的非線性彈性影響用傳遞矩陣建立了大型汽輪機(jī)發(fā)電機(jī)組軸承負(fù)荷分配計(jì)算方法，王剛志等[5]認(rèn)為精確的潤滑計(jì)算應(yīng)采用連續(xù)梁軸承載荷計(jì)算方法，周瑞[6]建立了艦船推進(jìn)軸系校中多目標(biāo)優(yōu)化模型對軸承負(fù)荷進(jìn)行計(jì)算，程鵬[7]提出以示功圖為輸入條件計(jì)算出不同的曲柄轉(zhuǎn)角下其連桿大端軸承和主軸承的負(fù)荷，彭澤軍[8]建立了大型汽輪發(fā)電機(jī)組軸承負(fù)荷分配計(jì)算的矩陣方程對軸承負(fù)荷進(jìn)行了計(jì)算與分析?？傮w而言，船舶軸系軸承負(fù)荷計(jì)算方法一般僅考慮油膜剛度（且為常數(shù)）的影響，而很少考慮軸承油膜動力學(xué)響應(yīng)（動態(tài)），主要原因在數(shù)學(xué)模型、計(jì)算方法等方面構(gòu)建的困難程度，因此本文將針對該問題提出考慮油膜動力學(xué)響應(yīng)的船舶軸系軸承穩(wěn)態(tài)負(fù)荷耦合計(jì)算方法，為穩(wěn)態(tài)軸承負(fù)荷計(jì)算提供方法。

1 船舶軸系軸承穩(wěn)態(tài)負(fù)荷耦合數(shù)學(xué)模型

在船舶軸系穩(wěn)態(tài)運(yùn)行狀態(tài)，船舶軸系-軸承系統(tǒng)以油膜為介質(zhì)實(shí)現(xiàn)了力學(xué)平衡，同時又通過各軸承軸頸位置關(guān)系達(dá)到了軸承負(fù)荷重新分配和再平衡，因此船舶軸系軸承穩(wěn)態(tài)負(fù)荷耦合數(shù)學(xué)模型可表示為式（1），其中包含了幾何耦合關(guān)系和力學(xué)耦合關(guān)系。

式中，R為軸系軸承負(fù)荷矩陣，Ceff為軸承影響系數(shù)矩陣，R0為軸系軸承的熱態(tài)負(fù)荷矩陣，可通過船舶軸系校中計(jì)算方法直接求解；F 為軸系軸承油膜力矩陣，Z 為軸系軸承穩(wěn)態(tài)位置矩陣，F(xiàn)e,i、Fφ，i、φi、Zi、ci、εi分別是第i個軸承的偏心方向油膜力、垂直方向油膜力、偏轉(zhuǎn)角、軸頸中心位置、軸承間隙、偏心率等，且均可應(yīng)用油膜動力學(xué)方程和邊界條件，通過迭代求解得到。

2 混合求解計(jì)算流程設(shè)計(jì)

根據(jù)式（1）所包含數(shù)學(xué)模型的物理含義和數(shù)學(xué)求解方法，將船舶軸系軸承穩(wěn)態(tài)負(fù)荷計(jì)算過程分解為三個子模塊，分別為軸系校中計(jì)算模塊、油膜動力學(xué)求解模塊和耦合求解模塊，其中軸系校中計(jì)算模塊和耦合求解模塊均為解析求解，而油膜動力學(xué)求解模塊為迭代求解。三個子模塊之間通過數(shù)據(jù)流進(jìn)行耦合，其中軸系校中計(jì)算模塊計(jì)算結(jié)果分別傳遞至油膜動力學(xué)模塊和耦合求解模塊，Ceff傳遞至耦合求解模塊；同時，耦合求解模塊除了耦合處理其他兩個模塊計(jì)算結(jié)果外還將結(jié)果反饋到油膜動力學(xué)模塊并進(jìn)行新的迭代求解，如圖1所示。

圖1 模塊間數(shù)據(jù)流示意圖Fig.1 Diagram of data flow among modules

3 子功能模塊設(shè)計(jì)與處理

3.1 船舶軸系校中計(jì)算模塊

船舶軸系可簡化為超不靜定梁問題，且假設(shè)船舶軸系兩端均為自由端，則船舶軸系校中計(jì)算數(shù)學(xué)模型可以通過式（2）表示。

式（2）中 Mi、Pi、Zi、Ci分別為第 i節(jié)點(diǎn)的彎矩、集中力、軸承位置值和設(shè)定軸承位置值，Ii、Li、Ei、qi分別為第i單元的截面慣性矩、長度、材料的彈性模量和均布力。求解得到Mi值，通過式（3）即可計(jì)算得到軸系各軸承負(fù)荷Ri。

通過MATLAB編程技術(shù)，設(shè)計(jì)了船舶軸系校中計(jì)算模塊流程如圖2所示，其包含了幾何及物理模型數(shù)據(jù)庫文件bearing_data_model.xlsx，以 EXCEL 格式文件通過函數(shù)“xlsread”讀至計(jì)算程序，同時為了使程序具有更好的移植性，并結(jié)合實(shí)際計(jì)算需要編寫一系列自定義函數(shù)：幾何建模函數(shù)dimension_b（）、截面慣性矩計(jì)算函數(shù)I_cal（）、彎曲剛度計(jì)算函數(shù)EI_cal（）、均布力計(jì)算函數(shù)q_cal（）、系數(shù)矩陣A構(gòu)建函數(shù)A_matrix_b（）、系數(shù)矩陣C構(gòu)建函數(shù)C_matrix_b（）、軸承負(fù)荷計(jì)算函數(shù)bearing_F_cal（）等。在自定義函數(shù)基礎(chǔ)上，根據(jù)式（2-21）編制代碼 X_matrix=（A_matrix）^（-1）*C_matrix，計(jì)算得到軸系各截面的彎矩值和撓度值，最終得到R0和Ceff。

3.2 油膜力計(jì)算模塊

油膜力計(jì)算模塊中采用了無量綱化的滑動軸承油膜動力學(xué)數(shù)學(xué)模型即式（4）和Sommerfeld邊界條件即式（5）。通過數(shù)值計(jì)算可求解得到無量綱壓力P在θ和Y上的分布規(guī)律，在獲得求解域內(nèi)的壓力分布后，即可進(jìn)行油膜穩(wěn)態(tài)工況其他特征參數(shù)求解。

圖2 船舶軸系校中計(jì)算模塊流程圖Fig.2 Flow-process diagram on calculation module of ship shafting alignment

式（4）和式（5）中P、Y、?分別表示無量綱壓力、長度和偏心率，θ為偏轉(zhuǎn)角，φ周向角度，φ2為油膜破裂周向角。

軸承油膜力計(jì)算模塊以船舶軸系校中熱態(tài)軸承負(fù)荷、軸承幾何參數(shù)、潤滑油物性參數(shù)、軸轉(zhuǎn)速等作為輸入計(jì)算參數(shù)，以EXCEL文件形式保持，通過函數(shù)“xlsread”讀至計(jì)算程序，計(jì)算后得到油膜力、偏心率和軸頸位置數(shù)據(jù)；在模塊中包含了一系列自定義函數(shù):軸承參數(shù)構(gòu)建函數(shù)bearing_para_b（）、單軸承油膜動力學(xué)參數(shù)計(jì)算函數(shù)single_bf_cal（）、油膜力粗算函數(shù)F_cal（）、油膜壓力分布計(jì)算函數(shù)dispressure（）、作用力積分函數(shù)totalpressure（）、油膜力精算函數(shù)Fnext_cal（）等，具體流程如圖3所示。

由于油膜力計(jì)算精度與偏心率變化幅值數(shù)量級大小相關(guān)，若選擇不合適，將導(dǎo)致數(shù)據(jù)結(jié)果振蕩而無法收斂，為了解決這個問題，本文采用了逐級降低偏心率變化幅值數(shù)量級方法，實(shí)現(xiàn)了快速收斂計(jì)算并保證了計(jì)算精度。

3.3 耦合求解模塊

耦合求解模塊是以船舶軸系校中計(jì)算模塊和油膜力計(jì)算模塊的計(jì)算結(jié)果為基礎(chǔ)，根據(jù)油膜力和軸承負(fù)荷差值及其大小關(guān)系，確定偏心率微調(diào)策略，通過循環(huán)迭代最終達(dá)到軸系軸承負(fù)荷與油膜力的基本相等，以其誤差設(shè)定值（本文設(shè)定值為200 N）為收斂判據(jù)，其具體流程如圖4所示。在耦合求解模塊中除了調(diào)用其他兩個模塊數(shù)據(jù)外，還需反復(fù)調(diào)用兩個函數(shù)：油膜力計(jì)算調(diào)用函數(shù)F_e_cal（）和軸承負(fù)荷計(jì)算函數(shù)Rnext_cal（）。

耦合求解模塊關(guān)鍵問題是步長選擇與收斂判據(jù)的問題，本文首先對軸系軸承負(fù)荷與油膜力差值的絕對值進(jìn)行大小排序，借鑒船舶軸系合理校中計(jì)算的軸承變位調(diào)整思路，選擇其中差值最大的兩個軸承進(jìn)行軸頸偏心率微調(diào)量作為求解步長，以前后迭代值之差作為收斂判據(jù)。由于實(shí)際差值有可能為“正”或“負(fù)”，因此對應(yīng)偏心率調(diào)整也具有方向性；若軸承負(fù)荷差值為正，則其對應(yīng)軸頸偏心率將增加一個偏心率增幅即“delt_e”，反之則減小一個偏心率增幅。同時為了提高收斂速度和計(jì)算精度，軸頸偏心率增幅即求解步長，根據(jù)前后兩次迭代值之差進(jìn)行數(shù)量級調(diào)整即當(dāng)“delta_FR”大于2 000 N時求解步長為0.001，隨著“delta_FR”下降至2 000 N時求解步長為0.000 1，進(jìn)一步“delta_FR”小于1 000 N時求解步長調(diào)整為0.000 01。

圖3 油膜力計(jì)算模塊流程圖Fig.3 Flow-process diagram on calculation module of oil film force

圖4 耦合求解模塊流程圖Fig.4 Flow-process diagram of the coupling calculation module

4 計(jì)算結(jié)果分析與驗(yàn)證

以某64 000 DWT散貨船為例對船舶軸系軸承穩(wěn)態(tài)負(fù)荷耦合計(jì)算方法進(jìn)行分析，船舶軸系軸承穩(wěn)態(tài)負(fù)荷計(jì)算過程涉及到軸承、軸頸幾何尺寸及長細(xì)比等參數(shù)，見表1，表中包含了8個軸承及其相應(yīng)軸頸的基本幾何尺寸，依據(jù)長細(xì)比分類原則，除了尾軸承為長軸承外，均為有限長軸承，“*”表示等效直徑。船舶軸系軸承穩(wěn)態(tài)負(fù)荷計(jì)算除與軸頸、軸承、軸承間隙等幾何尺寸相關(guān)外，還與軸承潤滑油物性參數(shù)、工作條件等因素關(guān)系密切，根據(jù)船廠提供的滑油參數(shù)可知正常運(yùn)行溫度為40℃，滑油粘度為0.092 4 Pa.s。

表1 船舶軸系軸承的幾何參數(shù)表Tab.1 Geometric parameters on bearings of the marine shafting

通過MATLAB計(jì)算得到軸系軸承穩(wěn)態(tài)負(fù)荷變化曲線，如圖5所示，圖中8個軸承穩(wěn)態(tài)負(fù)荷隨著轉(zhuǎn)速變化而變化，其中1號軸承（尾軸承）穩(wěn)態(tài)負(fù)荷隨著轉(zhuǎn)速升高略有下降、中間軸承穩(wěn)態(tài)負(fù)荷隨著轉(zhuǎn)速升高逐漸增大，而主機(jī)主軸承穩(wěn)態(tài)負(fù)荷有部分增加、有部分減小。這些軸承負(fù)荷變化曲線與傳統(tǒng)軸系校中計(jì)算得到軸承負(fù)荷為定值在變化規(guī)律上體現(xiàn)已是完全不同，體現(xiàn)了船舶軸系軸承負(fù)荷在軸承油膜動力學(xué)響應(yīng)影響下軸系軸承重新分配的特征

為了驗(yàn)證數(shù)學(xué)模型和計(jì)算程序，進(jìn)行了實(shí)船試驗(yàn)并整理得到2號軸承橫向間隙變化的實(shí)測數(shù)據(jù)，并與理論計(jì)算值變化曲線進(jìn)行了比較，如圖6所示。通過比較分析可得：中間軸承橫向間隙的測量值與應(yīng)用船舶軸系軸承穩(wěn)態(tài)負(fù)荷耦合數(shù)學(xué)模型計(jì)算后得到的理論值基本吻合，驗(yàn)證了本文采用的數(shù)學(xué)模型、編程思想和計(jì)算程序。

圖5 軸系軸承穩(wěn)態(tài)負(fù)荷值變化曲線Fig.5 Curves on steady load of shafting bearings

5 結(jié)語

船舶軸系軸承穩(wěn)態(tài)負(fù)荷受到軸系幾何結(jié)構(gòu)、船舶運(yùn)行狀態(tài)、滑油參數(shù)等多重因素的影響，計(jì)算過程非常復(fù)雜，通過MATLAB編程可根據(jù)數(shù)學(xué)模型及耦合關(guān)系較容易地實(shí)現(xiàn)了計(jì)算和得到了預(yù)期結(jié)果，并大大節(jié)省了計(jì)算周期。本計(jì)算程序具有通用性，可直接用于船舶軸系校中計(jì)算和滑動軸承油膜動力學(xué)特性分析，也可延伸到船舶軸系軸承動態(tài)負(fù)荷計(jì)算。

圖6 中間軸承橫向間隙理論值與測量值比較圖Fig.6 Comparison on transverse clearance between theoretical value and measured value