藺紅鵬 杜發(fā)榮 徐 征 周 煜

（北京航空航天大學能源與動力工程學院 北京 100191）

引言

二沖程航空重油發(fā)動機作為目前最具發(fā)展?jié)摿Φ男⌒突钊l(fā)動機之一，它繼承了活塞式航空發(fā)動機體積小、重量輕、功率密度大等優(yōu)點[1-2]，更適合機體空間有限的無人機[3-5]。航空重油是指餾分在航空煤油與柴油之間的航空油料[6]，與傳統(tǒng)航空汽油相比，重油具有粘度高，揮發(fā)性差等特點，從而提高了燃料使用和儲運的安全性，同時隨著燃料一體化的發(fā)展趨勢，重油的使用更有利于簡化后勤保障，在軍隊和通用航空領域應用前景廣闊[6-9]。

本研究中某型二沖程航空重油發(fā)動機，由于其將裝備于固定翼飛機或直升機作為動力，因此為保證飛機的飛行安全，對航空發(fā)動機提出更為嚴苛的要求[10]，如從動力輸出上，要有比地面動力更好的平順性和平衡性[11-12]，否則由此產(chǎn)生的不平衡力及力矩將轉(zhuǎn)變成翼載荷；同時要求發(fā)動機具有超高的可靠性，不能出現(xiàn)空中停車等重大事故[13]。曲柄連桿機構作為內(nèi)燃機的核心，歷來受到內(nèi)燃機研究人員的重視[14]，因此對曲柄連桿機構的動力學分析非常必要。

目前國內(nèi)外對曲柄連桿機構的動力學分析方法很多，且比較成熟和完善。曲柄連桿機構動力學分析方法有解析法、圖解法和數(shù)值方法[15]。解析法精度高，但對復雜機構工作量大，且實現(xiàn)起來有困難；圖解法形象直觀，但精度不夠，適用性較差。隨著計算機技術的不斷發(fā)展，用數(shù)值方法求解高度非線性的微分方程組成為了可能，使得復雜工程問題可得到更精確的解。前人在發(fā)動機動力學分析方面做了很多工作[16-21]，如張文燦[16]基于單缸機的一維動力學分析只基于一維模型進行，不能具體反映每個部件詳細的應力分布及變形；三維的數(shù)值仿真對模型做出很多假設，忽略部件彈性變形，如吳楠[19]將部分或全部部件視為剛體，或?qū)⑶B桿機構每個部件分開單獨進行彈性的有限元分析，雖然考慮了部件的彈性變形，但對于數(shù)值計算的邊界條件和非線性因素忽略較多，不能真實反映部件的受力和部件間的相互作用，在模型分析的完整性和精確度上都有不同程度的局限性。

本文提出基于GT-POWER與ABAQUS結合的動力學分析方法，更全面地對模型進行覆蓋360°全循環(huán)載荷分析，基于ABAQUS的三維分析采用運動彈性動力學方法，將所有部件視為彈性體，并考慮幾何與接觸等非線性因素，在基于GT-POWER的一維分析基礎上，以一維分析結果作為載荷邊界條件，對受載嚴重的沖程重點分析，可以減少計算點，更高效地得到更為精確的計算結果。分析中重點考慮基于二沖程和航空發(fā)動機的獨特要求，填補針對二沖程重油航空發(fā)動機動力學分析研究的空白，對重油航空發(fā)動機相關結構設計具有重要意義。

1 基于GT-POWER的發(fā)動機曲柄連桿機構一維動力學分析

1.1 建立一維熱力循環(huán)及動力學分析模型

GT-POWER涵蓋了發(fā)動機的熱力、機械、流體和結構等各個方面，內(nèi)燃機的模型庫豐富又完整，且計算速度快，功能強大。本文首先根據(jù)二沖程重油發(fā)動機的循環(huán)參數(shù)，應用GT-POWER建立發(fā)動機一維熱力循環(huán)模型，并進行參數(shù)優(yōu)化，建立一個符合性能要求的發(fā)動機熱力循環(huán)模型如圖1所示。本研究的對象為一臺直列兩缸的二沖程航空發(fā)動機，曲軸旋轉(zhuǎn)360°CA完成一個循環(huán)，兩缸發(fā)火順序相位差180°CA。發(fā)動機的基本參數(shù)如表1所示。

表1 發(fā)動機技術規(guī)格表

在發(fā)動機熱力循環(huán)模型的基礎上建立GTPOWER的動力學分析模型。建模過程中一些參數(shù)需根據(jù)結構進行計算或?qū)嶋H情況確定，如曲柄臂轉(zhuǎn)動慣量、曲柄銷慣性矩以及軸承的徑向容差和潤滑油性能的選取等。建立的模型如圖1所示，在兩個模型聯(lián)合仿真的基礎上得到發(fā)動機整個循環(huán)內(nèi)的動力學參數(shù)，并與理論計算結果作對比驗證其正確性。

1.2 一維動力學計算結果

在動力分析中，不計自重和摩擦阻力，主要分析氣壓力、往復和旋轉(zhuǎn)慣性力等在曲柄連桿機構中的作用情況，分析活塞側壓力的大小和規(guī)律，連桿軸頸以及曲軸主軸頸等運動副的載荷分布規(guī)律。計算結果如下：

1）發(fā)動機氣缸內(nèi)的工質(zhì)壓力是整個動力裝置的原動力，氣壓力-曲軸轉(zhuǎn)角曲線如圖2所示。

圖1 發(fā)動機GT-POWER一維仿真模型

圖2 發(fā)動機氣壓力-曲軸轉(zhuǎn)角曲線

2）曲柄連桿機構各部件運動中產(chǎn)生慣性力，活塞往復運動的慣性力是導致發(fā)動機不平衡的主要原因，曲軸定速旋轉(zhuǎn)其旋轉(zhuǎn)慣性力為定值302N，如圖3所示。

3）發(fā)動機各運動副的載荷及其分布

圖3 慣性力曲線

作用在活塞頂上的氣壓力和運動件慣性力的最終作用表現(xiàn)在各個運動副所受的載荷及其分布規(guī)律。如圖4、圖5所示，其中活塞銷載荷與連桿小頭載荷互為反作用力；曲柄銷載荷與連桿瓦載荷互為反作用力。

圖4 活塞銷載荷分布

圖5 連桿瓦載荷分布

2 基于裝配體下的ABAQUS曲柄連桿機構動力學分析

2.1 運動彈性動力學理論與數(shù)學模型

彈性動力分析的任務是研究曲柄連桿機構在外部載荷和自身慣性力共同作用下的運動和受力情況。應用運動彈性動力學分析方法建立有限元模型，不僅考慮了系統(tǒng)的彈性變形，建立更精確的運動微分方程，而且具有有限元的運算模式統(tǒng)一，模型適應性廣的特點。本研究應用基于“瞬時結構假定”的運動彈性動力學分析方法，即曲柄連桿機構在運動循環(huán)中某一位置，將機構的形狀和作用于其上的負荷瞬時“凍結”，從而將機構利用結構的方法分析[22]?；静襟E是：

1）將曲柄連桿機構各構件劃分單元，在各單元指定點設置結點；

2）選擇位移模式，建立廣義坐標，從拉格朗日方程導出單元的運動微分方程；

3）將單元運動方程集合為系統(tǒng)的運動方程，其一般形式如下：

式中：[M]為系統(tǒng)質(zhì)量矩陣；[C]為系統(tǒng)阻尼矩陣；[K]為系統(tǒng)剛度矩陣；{F}為系統(tǒng)載荷矩陣；{x}為系統(tǒng)廣義坐標列陣。

式（1）的實質(zhì)是廣義的牛頓第二定律，機械振動學已經(jīng)給出了上式的解法。這個非齊次微分方程組的解包括兩部分，一是外載荷為零的自由振動齊次解，由于振動能量耗散，振動很快消失；另一個是系統(tǒng)在外界激勵下的穩(wěn)態(tài)響應。將上式解耦：

式（2）中{ψN}為正則振型矩陣，也即形函數(shù)矩陣，代入系統(tǒng)運動微分方程，并左乘{ψN}T，化簡后得：

方程解耦完成，令qNi=xqieiωt，代入上式：

式中：λi=ω/ωi，i=1，2，…，n。

求出系統(tǒng)運動微分方程可得到廣義坐標x，進而根據(jù)幾何方程、物理方程等求出曲柄連桿的應力應變及其真實的位移、速度、加速度以及零部件的變形等。

2.2 有限元模型的建立與邊界條件

發(fā)動機曲柄連桿機構系統(tǒng)主要包括活塞、活塞銷、連桿以及曲軸等，不同部件采用不同材料，具體如表2所示。由于活塞為非對稱結構，為提高仿真精度，考慮曲軸潤滑油道及曲柄銷減重孔，因此采用全模型。結合發(fā)動機的基本參數(shù)和制造工藝，在CATIA中建立基于特征控制的發(fā)動機三維模型，在Hypermesh中進行結構離散。其中缸套和活塞銷為六面體一階單元C3D8R，其余部件采用C3D4四面體單元，使網(wǎng)格相對于模型具有較高的幾何相符率。對預計應力較大區(qū)域進行網(wǎng)格加密，如活塞銷孔內(nèi)沿及支撐肋等。共生成664 028個單元，185 388個節(jié)點，如圖6所示。

表2 曲柄連桿機構材料表

圖6 曲柄連桿機構有限元模型

發(fā)動機循環(huán)的動力學分析采用基于“瞬時結構假定”的運動彈性動力學分析方法，每隔10°CA選取一個計算點，對不同計算點的模型采用不同的邊界條件。

為消除模型剛體位移，曲軸主軸頸6個自由度全約束，保證結構穩(wěn)定。載荷邊界條件為兩個活塞頂面承受來自燃燒室的氣壓力，可通過發(fā)動機熱力循環(huán)軟件GT-POWER仿真得到?？紤]運動副存在間隙或接觸的狀態(tài)非線性，在氣缸套和活塞裙部、連桿大頭和曲柄銷等接觸面設置接觸對，接觸對主從面合理的選擇不僅可以保證數(shù)值計算的順利進行，還可減少接觸搜尋算法的計算成本。

2.3 ABAQUS動力學分析結果

應用ABAQUS/Standard求解器通過運動彈性動力學方法求解，獲得曲柄連桿機構變形及應力結果。由于在氣缸達到最高爆發(fā)壓力11 MPa時（上止點后5°CA）機構受載最惡劣，因此取上止點后5°CA、40°CA及75°CA對其應力應變、變形及載荷作分析。如圖7所示。

從圖7a看出曲軸轉(zhuǎn)角為5°CA時機構受載最惡劣，部分應力超出了材料的許用值，但經(jīng)過結構優(yōu)化配合熱處理工藝可消除危險截面。如連桿小頭出現(xiàn)基于赫茲理論的邊沿效應，應力達到1 006 MPa，本文對此結構優(yōu)化后重新計算，應力值降為407.43 MPa，如圖13所示，再配以熱處理及噴丸等工藝可以進一步提高連桿小頭承載能力。優(yōu)化后的機構在最高爆發(fā)壓力下應力最大且安全，隨著曲軸轉(zhuǎn)角的繼續(xù)增大，應力也相應減小，如圖18所示，因此在其它轉(zhuǎn)角下機構都是安全的。

曲軸在上止點后5°CA、40°CA及75°CA的應力分布如圖8所示。

圖7 曲柄連桿機構在5°CA、40°CA及75°CA時應力分布

圖8 曲軸部件在5°CA、40°CA及75°CA時應力分布

從圖8a知，曲柄銷減重孔、曲柄銷和曲柄臂連接處及主軸頸支撐點應力較大，最大為287.4 MPa，遠小于曲軸材料Cr40的屈服極限，因此在最嚴苛工況下曲軸仍滿足強度要求。曲軸在工作中承受周期性動載荷，因此不僅要考察應力最大點，更應注意應力梯度較大的位置，這些位置容易產(chǎn)生振動和疲勞[23]，是曲軸最易破壞的區(qū)域?？疾烨罩忻媲鄣膽μ荻?，如圖9所示。作一個循環(huán)內(nèi)曲柄臂應力梯度曲線圖，如圖10所示，應力梯度與缸內(nèi)壓力趨勢基本一致，在最高爆發(fā)壓力下應力梯度最大，為121.97 MPa。可以看出，一個循環(huán)內(nèi)曲拐中面的應力梯度變化較大，容易在曲柄臂圓角處產(chǎn)生疲勞破壞。

連桿在上止點后5°CA、40°CA及75°CA的應力分布如圖11所示。

圖9 曲拐中面應力梯度云圖

圖10 曲拐中面應力梯度循環(huán)曲線圖

圖11 連桿在5°CA、40°CA及75°CA時的應力分布

由圖11可知，連桿小頭承壓部位由于活塞銷變形引起邊緣效應，邊緣超出許用應力，優(yōu)化措施是在小頭內(nèi)孔邊沿作1×45°倒角，并采用局部硬化的熱處理，減小應力集中。對裝配了優(yōu)化后的連桿機構進行計算，連桿小頭應力分布對比如圖12、圖13所示。

圖12 結構優(yōu)化前應力分布

圖13 結構優(yōu)化后應力分布

由圖13可知結構優(yōu)化后連桿應力分布正常，小頭孔下邊沿最大應力為407.43 MPa，小于42CrMoV材料588 MPa的屈服極限，實際連桿還需熱處理噴丸等強化工藝。連桿是二力桿，對變形較為敏感，長度方向變形量不能超過20μm，否則影響發(fā)動機裝配和正常運作。如縱向變形過大不僅引起襯套咬合甚至可能發(fā)生活塞碰氣門等故障。任意節(jié)點實際位移=剛體位移+變形位移，因此考慮兩點的變形時須找一中性基準點，此點只有剛體位移而沒有變形位移，選擇所關心點的絕對位移與中性點絕對位移的差便是該點的變形量。以5°CA曲軸轉(zhuǎn)角為例，選取連桿上具有代表性的小頭節(jié)點、大頭上端節(jié)點的變形數(shù)據(jù)，處理得到兩點基于中性點的變形量，如圖14所示。

由圖14曲線知，連桿變形量基本符合要求，不影響發(fā)動機的裝配和正常工作，同時又做到了連桿的輕量化。

運動副載荷以1#缸為例，在5°CA、40°CA及75°CA轉(zhuǎn)角下曲柄銷的載荷大小分別為31 705.19 N、9 456.52 N和2 514.7 N，活塞側推力大小分別為763.56 N、1 554.91 N和724.04 N，與理論計算結果符合較好。

圖14 連桿縱向變形量曲線

3 結果對比分析與討論

3.1 發(fā)動機GT-POWER與ABAQUS動力學分析結果的對比驗證

圖15 活塞側推力的一三維載荷曲線

二沖程為平面曲軸，兩缸載荷的分析結果大小方向一致，僅相位差180°CA，因此僅以1#缸作分析。

提取GT-POWER動力學分析結果及對ABAQUS分析結果后處理，得到1#缸活塞側推力及曲柄銷載荷在一個循環(huán)內(nèi)的一三維結果對比，如圖15、圖16中曲線所示。理論計算與基于ABAQUS的分析在大小和方向上基本一致，GT-POWER的一維分析與三維的ABAQUS在趨勢上符合較好，但與理論結果差異較大，如圖17所示，最大差異在3%左右。其原因主要是兩種計算方法的不同，基于有限元的ABAQUS動力學分析相對假設條件較少，且網(wǎng)格劃分合理時其誤差主要來源于截斷誤差；一維動力學分析需輸入較多設定參數(shù)，其中一些參數(shù)根據(jù)同級別發(fā)動機采用經(jīng)驗值，或在計算零部件慣量等參數(shù)時對幾何模型做了簡化。因此由于條件限制GTPOWER采用了較多的假設及簡化，使計算結果出現(xiàn)較大誤差。但相對理論或三維結果，其平均偏差保持在3%以內(nèi)，對精度要求不高的分析仍具有一定參考意義。

需注意的是1#缸活塞側推力在一個循環(huán)內(nèi)有兩個峰值，這對發(fā)動機平衡性和動力輸出特性都是有害的，可加裝飛輪甚至平衡軸來改善。GT-POWER與ABAQUS關于1#缸曲柄銷載荷的分析結果相符較好。通過ABAQUS對發(fā)動機動力學分析，可得到曲柄銷整環(huán)的應力分布，應力最小位置適宜布置潤滑油孔，從而改善曲柄銷的潤滑條件，使油膜更容易形成，提高發(fā)動機的使用壽命。

圖16 曲柄銷載荷的一三維載荷曲線

圖17 不同載荷形式一三維結果差異曲線

3.2 曲柄連桿機構最大應力分布曲線

曲柄連桿機構由不同材料部件組成，不同材料屈服極限也有較大差異，應力大小對于材料的強度剛度是與材料緊密相關的。研究40Cr的曲軸，分析其應力大小及其位置，提取曲軸一個循環(huán)內(nèi)每個計算點的最大應力，判斷其每一曲軸轉(zhuǎn)角下的可靠性?？傻萌鐖D18所示的最大應力與曲軸轉(zhuǎn)角曲線。

圖18 40Cr曲軸的最大應力曲線

由于采用的是基于“瞬時結構假定”的運動彈性動力學分析方法，因此圖18曲線是由分布密集的有限個計算點擬合而成的樣條曲線。通過后期增加的計算點的驗證，其結果符合較好，通過該曲線可以預測一個循環(huán)內(nèi)其他點的最大應力位置。

3.3 針對二沖程航空發(fā)動機的動力學結構設計改良方案

曲拐布置及點火正時、配氣與供油系統(tǒng)的匹配，決定了二沖程發(fā)動機結構型式、運行方式與常見的四沖程發(fā)動機有較大差異。二沖程發(fā)動機一個循環(huán)曲軸旋轉(zhuǎn)360°CA，而四沖程為720°CA，因此二沖程發(fā)動機運動副等所受交變載荷的頻率是四沖程的兩倍。同時造成熱負荷比四沖程發(fā)動機高，總之二沖程發(fā)動機的運行工況更為惡劣。因此針對二沖程航空發(fā)動機的動力學結構設計具有更嚴苛的標準，不僅要求良好的配氣、掃氣及燃燒組織，同時針對動力學的分析應盡可能接近發(fā)動機真實的動力學過程，減少假設和簡化，對重點區(qū)域詳細分析。如轉(zhuǎn)動副載荷的分析指導滑動軸承的結構設計，保證良好的壓力潤滑，提高發(fā)動機運行品質(zhì)和壽命；側壓力的分析指導活塞裙部型線的優(yōu)化，使活塞在正常運行中形成良好的楔形油膜，從而延長缸套壽命和提高發(fā)動機的效率。

4 結論

1）在熱力循環(huán)的基礎上，基于GT-POWER的動力學分析，更側重部件間運動副載荷大小隨曲軸轉(zhuǎn)角的變化規(guī)律，對發(fā)動機在整個循環(huán)內(nèi)的載荷分布有直接的反映；基于有限元理論，應用運動彈性動力學理論的ABAQUS動力學分析，更注重某個瞬態(tài)的應力分布和部件的變形情況，對部件每一部分受力有更詳細的分析。GT-POWER與ABAQUS結果誤差保持在3%以內(nèi)。

2）GT-POWER的示功圖等輸出參數(shù)可作為ABAQUS分析的載荷邊界條件。GT-POWER與ABAQUS的聯(lián)合仿真，結合一維分析計算速度快及三維分析計算精度高的優(yōu)勢，對特定動力學過程進行重點分析，不僅可以得到三維非線性分析的精確結果，同時可減少三維計算量，節(jié)約計算成本。

3）一維和三維結合的動力學分析方法雖在滿足精度的同時減少了計算量，但目前無試驗數(shù)據(jù)與之對比驗證。對于聯(lián)合仿真結果的試驗數(shù)據(jù)驗證以后可以進一步作深入研究。

4）聯(lián)合仿真結果可指導結構設計及優(yōu)化。如在連桿小頭內(nèi)孔邊沿作1×45°倒角后，最大應力由邊沿效應引起的超過1 000 MPa降為407.43 MPa，連桿縱向變形量小于20μm，保證了發(fā)動機的正常運轉(zhuǎn)；最高爆發(fā)壓力11 MPa下，曲柄銷載荷最大為31 705.19 N，活塞側推力為763.56 N，根據(jù)計算結果，可指導曲柄銷壓力潤滑油孔位置確定，保證發(fā)動機良好潤滑，從而提高了壽命，同時根據(jù)側推力可為發(fā)動機傾覆力矩的平衡提供參考，優(yōu)化發(fā)動機的動力學輸出特性。