一維和三維聯合動力學仿真方法在活塞發動機曲柄連桿機構動力學分析中的應用*

藺紅鵬 杜發榮 徐 征 周 煜

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

引言

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

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

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

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

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

1.1 建立一維熱力循環及動力學分析模型

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

表1 發動機技術規格表

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

1.2 一維動力學計算結果

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

1）發動機氣缸內的工質壓力是整個動力裝置的原動力，氣壓力-曲軸轉角曲線如圖2所示。

圖1 發動機GT-POWER一維仿真模型

圖2 發動機氣壓力-曲軸轉角曲線

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

3）發動機各運動副的載荷及其分布

圖3 慣性力曲線

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

圖4 活塞銷載荷分布

圖5 連桿瓦載荷分布

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

2.1 運動彈性動力學理論與數學模型

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

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

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

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

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

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

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

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

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

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

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

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

表2 曲柄連桿機構材料表

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

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

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

2.3 ABAQUS動力學分析結果

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

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

曲軸在上止點后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的屈服極限，因此在最嚴苛工況下曲軸仍滿足強度要求。曲軸在工作中承受周期性動載荷，因此不僅要考察應力最大點，更應注意應力梯度較大的位置，這些位置容易產生振動和疲勞[23]，是曲軸最易破壞的區域。考察曲拐中面曲柄臂的應力梯度，如圖9所示。作一個循環內曲柄臂應力梯度曲線圖，如圖10所示，應力梯度與缸內壓力趨勢基本一致，在最高爆發壓力下應力梯度最大，為121.97 MPa。可以看出，一個循環內曲拐中面的應力梯度變化較大，容易在曲柄臂圓角處產生疲勞破壞。

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

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

圖10 曲拐中面應力梯度循環曲線圖

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

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

圖12 結構優化前應力分布

圖13 結構優化后應力分布

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

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

運動副載荷以1#缸為例，在5°CA、40°CA及75°CA轉角下曲柄銷的載荷大小分別為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 發動機GT-POWER與ABAQUS動力學分析結果的對比驗證

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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