于存銀 國 杰 趙 建 張 恒 王仁智

（1.中國船舶及海洋工程設(shè)計研究院 上海 200011；2.哈爾濱工程大學(xué) 動力與能源工程學(xué)院 哈爾濱 150001）

0 引 言

內(nèi)燃機的活塞-連桿-曲軸-軸承和各個部件連接處的油膜是一個復(fù)雜的系統(tǒng)，燃?xì)鈮毫瓦\動件的慣性力通過連桿作用在曲柄銷上，最后力通過曲軸—軸承系統(tǒng)傳遞到整機。整個系統(tǒng)中，主軸承和曲軸之間通過油膜連接，油膜作為曲軸運動的動態(tài)邊界條件，一方面避免其與主軸承直接產(chǎn)生干摩擦，另一方面油膜充當(dāng)非線性的彈簧阻尼可以減小曲軸的振動。此外，曲軸的運動也會直接影響主軸承油膜的潤滑參數(shù)，決定潤滑質(zhì)量。所以內(nèi)燃機的曲軸-油膜-軸承系統(tǒng)是一個強耦合的系統(tǒng)，兩者相互作用、相互影響。

在考慮油膜潤滑求解軸心軌跡中，最早采用了3 種靜力學(xué)方法：Hahn 方法、Holland 方法以及遷移率法，然而這3 種方法并未考慮軸頸的動力學(xué)效應(yīng)。1982 年，陳伯賢和裘祖干將軸頸的慣性項計入軸心軌跡的動平衡方程中，進(jìn)一步優(yōu)化了動載滑動軸承軸心軌跡的算法。2003年，戴旭東等通過ADAMS 和自編Fortran 程序的聯(lián)合分析了曲軸剛體動力學(xué)和主軸承潤滑耦合的特性。之后，合肥工業(yè)大學(xué)桂長林等又在此基礎(chǔ)上結(jié)合曲軸的柔性體動力學(xué)、曲軸軸頸傾斜等因素研究了曲軸的軸心軌跡和主軸承潤滑特性，并采用柔度矩陣法進(jìn)行了計入軸瓦彈性變形對軸承潤滑的影響。方斌針對實船水潤滑軸承總結(jié)出低黏彈流動壓潤滑剛度系數(shù)，并考慮了軸承支撐剛度的影響。2017年，大連海事大學(xué)的魏立隊等利用商業(yè)軟件建立了機體-曲軸的彈流潤滑（Elastic Hydrodynamic Lubrication，EHD）耦合柔性體模型，結(jié)合有限體積法和有限元方法研究了船舶柴油機軸系振動特性。為了適應(yīng)內(nèi)燃機設(shè)計的快速化、簡便化和低成本要求，實現(xiàn)建模方法的通用化和計算方法的模塊化是進(jìn)行內(nèi)燃機高性能設(shè)計的必然要求。2000 年，MOURELATOS Z P采用子結(jié)構(gòu)法對曲軸和機體進(jìn)行有限元建模，然后利用NASTRAN 進(jìn)行縮減模態(tài)自由度，并將提取的質(zhì)量、剛度和阻尼矩陣數(shù)據(jù)輸入到Fortran 里計算，大大減少了計算量并提高了計算速度。2003年，密西根大學(xué)的MA Zhengdong和PERKINS N C針對某整機的簡化計算模型，提出一種快速計算多體動力學(xué)的方法，但計算局限于穩(wěn)態(tài)條件下。2006 年，我國的雷宣揚等基于有限梁單元構(gòu)造了整體曲軸模型，從而實現(xiàn)快速建模和仿真。

基于之前學(xué)者的理論和研究成果，為了實現(xiàn)建模計算的快速化和能夠方便結(jié)構(gòu)動力學(xué)的修改，本文提出了一種新的曲軸-軸承系統(tǒng)的建模計算方法，利用Fortran 語言編寫計算程序。在之前學(xué)者采用基于靜力變形的柔度矩陣法研究的彈流潤滑基礎(chǔ)上，借助于APDL 考慮了軸承組的瞬態(tài)動力學(xué)效應(yīng)，使彈流潤滑仿真計算更貼近實際物理工況，并且和不考慮軸承組變形的流體動力潤滑（Hydrodynamic Lubrication，HD）計算結(jié)果進(jìn)行了對比研究。

1 模型建立和理論基礎(chǔ)

1.1 曲軸-軸承系統(tǒng)的激勵源計算

圖1 為4102 柴油機的活塞-連桿-曲軸模型。

圖1 4102 柴油機活塞-連桿-曲軸模型

如該圖所示，曲軸-軸承系統(tǒng)的主要激勵源是作用于活塞上的燃?xì)獗l(fā)壓力和活塞-連桿等運動部件的慣性力，其動力學(xué)計算簡化示意圖如下頁圖2所示。

圖2 簡化示意圖

曲柄銷處作用力，連桿推力沿著連桿傳遞到曲柄銷處，同時分解為曲柄銷中心處的切向力和法向力，由此得到作用在曲軸上的的載荷。

1.2 曲軸的有限元模型

1.2.1 梁單元

如圖3 所示，空間桿結(jié)構(gòu)的梁單元，每個梁單元有2 個節(jié)點，每個節(jié)點擁有6 個自由度，即3 個平動自由度和3 個轉(zhuǎn)動自由度。當(dāng)桿受到外部載荷時，桿會發(fā)生3 個方向的位移變形以及3 個平面的彎曲轉(zhuǎn)動變形。平動自由度的正方向是坐標(biāo)軸正方向，轉(zhuǎn)動自由度的正方向是以繞坐標(biāo)軸方向右手螺旋定理確定的方向為正。

圖3 空間桿結(jié)構(gòu)的梁單元

1.2.2 曲軸有限元模型的構(gòu)造

本文基于梁單元構(gòu)造了曲軸的有限元模型。我們將曲軸分為4 個曲拐，因為4 個曲拐的結(jié)構(gòu)以及材料參數(shù)相同，這樣我們可以采用子結(jié)構(gòu)法的思想，將曲軸分為4 個曲拐。只需要構(gòu)造出其中1 個曲拐的有限元模型，便可計算整個曲軸的有限元模型，因此我們先構(gòu)造出1 個曲拐的有限元模型。

我們將單個曲拐模型分為7 個梁單元，單個曲拐簡化之后如圖4 所示。

圖4 單曲拐的有限梁單元模型

由于不同曲拐之間只是相位差不同，將1 個單曲拐模型簡化成梁單元的超級單元之后，只需要整體調(diào)用梁單元集成的單曲拐有限元模型，從而實現(xiàn)整個曲軸有限梁單元的建立，其集成流程見圖5。

圖5 建模流程圖

經(jīng)過上述方法構(gòu)造出如下頁圖6 所示的曲軸整體有限梁單元模型。該模型的位移邊界條件是以油膜作支撐，載荷邊界條件是作用于曲柄銷上的激勵源，整體支撐的機體認(rèn)為是剛體。

圖6 曲軸有限梁單元模型

1.3 軸承潤滑理論

基于等溫、不可壓縮、等黏度的牛頓層流流體動載滑動軸承的Reynolds 方程為：

計入軸頸在軸承中發(fā)生傾斜時的油膜厚度方程：

圖7 軸頸傾斜的主軸承示意圖

軸承油膜反力和反力矩為：

1.4 軸承座模型

軸承座的材料是灰鑄鐵（HT250），灰鑄鐵的彈性模量是157 GPa，泊松比是0.27，密度為 7 250 kg/m，采用的單元類型是Solid185，主軸承座的位移邊界條件是在主軸承座的上部與機體連接的部分全部約束，認(rèn)為剛度無限大，有限元模型和約束位置見圖8。軸瓦表面的載荷邊界條件是每一載荷步下的油膜力施加到相應(yīng)的軸瓦表面節(jié)點上，軸瓦表面的網(wǎng)格節(jié)點和油膜的節(jié)點一一對應(yīng)。

圖8 軸承座有限元模型

2 數(shù)值計算方法

利用Fortran 語言自編的基于有限元方法計算曲軸動力學(xué)的程序和基于有限差分方法的求解軸承潤滑的程序兩者進(jìn)行耦合計算，求解步驟如下，計算流程見圖9。

圖9 計算流程圖

（1）設(shè)置主程序的計算仿真時間、計算步長、收斂精度和輸出結(jié)果的控制。

（2）調(diào)用曲軸動力學(xué)的程序，輸入曲軸的結(jié)構(gòu)參數(shù)和尺寸，設(shè)置曲軸的初始條件，包括主軸頸的軸心的位置、速度以及主軸承處的軸頸的傾斜角。

（3）調(diào)用主軸承油膜潤滑的程序，輸入主軸承的結(jié)構(gòu)參數(shù)和材料參數(shù)，設(shè)置計算的收斂精度。

（4）仿真開始，曲軸動力學(xué)輸出主軸頸的軸心位移、速度和主軸頸在主軸承處發(fā)生的傾斜程度，通過瞬態(tài)動力學(xué)求解軸瓦表面的受載的變形量修正油膜厚度，油膜潤滑程序得到曲軸的運行參數(shù)計算出主軸承的油膜壓力，通過積分計算出作用在曲軸上的油膜反力和力矩。

（5）仿真時間達(dá)到，計算結(jié)果穩(wěn)定，輸出計算結(jié)果，否則重新設(shè)定仿真時間，重新計算。

3 模型計算和結(jié)果分析

以4102 柴油機為研究對象，驗證上述計算方法。4102 柴油機的輸入計算參數(shù)見下頁表1，下頁圖10是氣缸壓力曲線，其計算的發(fā)火順序是1-3-4-2。

圖10 氣缸壓力曲線

表1 柴油機計算參數(shù)

3.1 曲軸主軸承油膜反力

圖11 和圖12 分別是通過EHD 和HD 計算的各個主軸承的油膜反力，由圖可見，總體趨勢一致。第2、4 主軸承的油膜反力峰值最大，第1、5 軸承次之，中間第3 主軸承的油膜力峰值最小。EHD 計算的主軸承反力峰值依次是26 547 N、34 811 N、22 432 N、33 512 N、27 654 N，HD 計算的主軸承反力峰值依次是25 986 N、36 075 N、22 987 N、36 780 N、26 883 N。由此可見：考慮軸承的瞬態(tài)彈性變形，整體的軸承反力是降低的，但是曲軸兩端主軸承的油膜反力整體提高了，這體現(xiàn)了曲軸動力學(xué)和軸承座動力學(xué)之間的耦合關(guān)系。

圖11 彈流潤滑的主軸承油膜反力

圖12 流體動力潤滑的主軸承油膜反力

3.2 曲軸主軸頸動力學(xué)特性

圖13 和圖14 所示分別是EHD 和HD 計算的主軸頸的受載傾斜程度，比較計算結(jié)果見表2。主軸頸在主軸承處由于受載變形產(chǎn)生的傾斜，從圖中可以發(fā)現(xiàn)主軸頸1、5 的傾斜程度最嚴(yán)重。

圖13 彈流潤滑的主軸頸的傾斜程度

圖14 流體動力潤滑的主軸頸的傾斜程度

表2 5 個主軸頸的動力學(xué)特性比較

這種情況主要是因為主軸頸1、5 處于曲軸兩端，只有一側(cè)因受外部載荷而彎曲變形，從而導(dǎo)致傾斜嚴(yán)重；而主軸頸2、3、4 兩邊同時受載，所以在一定程度上可實現(xiàn)兩側(cè)力和力矩的平衡抵消，尤其是主軸頸3 的傾斜程度最小。在EHD計算下，因油膜壓力的偏置使軸承發(fā)生彈性變形、油膜厚度增加，油膜壓力分布趨于減小，使得作用于主軸頸的傾斜力矩減小，所以整體傾斜程度減小。

以EHD 和HD 計算的曲軸徑向振動位移與比較計算結(jié)果見圖15、圖16 和表2。由于在外載荷的作用下，曲軸軸頸的運動規(guī)律是周期性的，曲軸的外載荷主要源于燃?xì)鈮毫屯鶑?fù)慣性力以及離心慣性力，但是振動規(guī)律主要還是取決于氣缸壓力和往復(fù)慣性力的強迫作用。由于主軸頸1 的相鄰氣缸是氣缸1，故其運動規(guī)律主要受到氣缸1的載荷影響，四缸機的發(fā)火順序是1-3-4-2，所以主軸頸1 在氣缸1 發(fā)火時的徑向振動最劇烈，受到其他氣缸的發(fā)火作用的效果由近及遠(yuǎn)；主軸頸2 同時受到氣缸1 和氣缸2 的作用，所以其徑向振動的峰值主要在180°和360°附近；受到柴油機對稱結(jié)構(gòu)的影響，主軸頸3 位于中間位置，所以其同時受到相鄰氣缸2 和氣缸3 的影響，180°和540°附近的徑向振動最劇烈，其次受到較遠(yuǎn)氣缸的影響，所以其振動規(guī)律具有極強的對稱性和周期性；因為主軸頸5、4 是和主軸頸1、2 以主軸頸3 為對稱，所以其振動規(guī)律是相似的，只是受到發(fā)火順序的影響，存在相位差，但幅值和趨勢是一致的，總體上的峰值和趨勢變化一致。

圖15 彈流潤滑的曲軸徑向振動位移

圖16 流體動力潤滑的曲軸徑向振動位移

EHD 計算的主軸頸徑向振動位移整體偏小，主要原因是軸承的瞬態(tài)彈性變形吸收了一部分能量，導(dǎo)致主軸承潤滑參數(shù)發(fā)生改變。主軸頸5 的振動位移變化和其他主軸頸的振動位移變化趨勢不一致，主要原因在于主軸頸5 發(fā)生了嚴(yán)重傾斜。

3.3 主軸承潤滑特性

下頁圖17 至圖20 所示分別是以EHD 和HD計算的主軸承的最小油膜厚度以及最大油膜壓力，具體計算結(jié)果差異比較參見表3。主軸承的最小油膜厚度和最大油膜壓力在一個工作循環(huán)內(nèi)的變化規(guī)律基本是負(fù)相關(guān)的。最小油膜厚度越小，其最大油膜壓力越大，變化規(guī)律基本符合油膜潤滑的動壓效應(yīng)，但是由于油膜壓力同時受到軸頸的擠壓效應(yīng)影響，所以從圖中也可以發(fā)現(xiàn)，即使油膜厚度相同，但是最大油膜壓力卻相差很大，這說明油膜擠壓效應(yīng)在流體潤滑計算中是不可忽略的。EHD 的計算結(jié)果顯示主軸承的最大油膜壓力都降低了，尤其是主軸承2 和 4 變化最明顯。主軸承2 的最大油膜壓力降低最大，但是最小油膜厚度并沒有增加。通過前頁的表2 也可判斷，有可能是軸頸2 發(fā)生嚴(yán)重傾斜，由此也證明了考慮軸頸傾斜的必要性。

表3 5個主軸承的潤滑特性比較

圖17 彈流潤滑的主軸承最小油膜厚度

圖20 流體動力潤滑的主軸承最大油膜壓力

圖18 流體動力潤滑的主軸承最小油膜厚度

圖19 彈流潤滑的主軸承最大油膜壓力

4 結(jié) 論

（1）通過有限元和有限差分的方法實現(xiàn)了活塞-連桿-曲軸-油膜-軸承系統(tǒng)的建模，實現(xiàn)了參數(shù)化和模塊化，能夠進(jìn)行系統(tǒng)的快速數(shù)值仿真，具有通用性，便于結(jié)構(gòu)動力學(xué)修改。

（2）整體曲軸-軸承系統(tǒng)的動力學(xué)特性和潤滑特性的變化都呈現(xiàn)對稱性和周期性的規(guī)律，這種現(xiàn)象和柴油機的結(jié)構(gòu)相符合。

（3）將彈流潤滑和流體動力潤滑的計算結(jié)果進(jìn)行了對比，發(fā)現(xiàn)考慮軸承的瞬態(tài)動力學(xué)計算的軸承變形對曲軸動力學(xué)以及軸承潤滑特性影響很大，主軸頸的徑向振動位移和傾斜程度都有減小趨勢，尤其第3 主軸承的軸頸傾斜程度減小23%；所有主軸承的最大油膜壓力都減小，尤其第2 和第4 主軸承的最大油膜壓力分別減小29.3%和23.9%。

（4）目前只是研究了考慮軸頸傾斜影響的曲軸動力學(xué)和軸承潤滑特性，而沒有考慮到在油膜載荷下的軸承座動力學(xué)特性，下一步將對軸承座的動變形進(jìn)行評價。