有限元彈性動力分析與優化在某氣門機構中的應用

孫先國，劉浩，倪計民

（1．同濟大學汽車學院，上海200092；2．舍弗勒集團發動機部，上海201804）

有限元彈性動力分析與優化在某氣門機構中的應用

孫先國1,2，劉浩2，倪計民1

（1．同濟大學汽車學院，上海200092；2．舍弗勒集團發動機部，上海201804）

對某發動機氣門機構建立了混合剛體-彈性體的參數化非線性有限元模型，進行了彈性動力學計算。通過回歸方程對模型進行了優化，降低了搖臂的質量，在保證搖臂強度的同時改善了凸輪-滾子的接觸性能，從而降低了磨損，提高了耐久性。

氣門機構有限元法彈性動力回歸優化

1 前言

以往的氣門機構分析，多是將氣門機構系統簡化成單質量或多質量集總參數模型[1-5]，然后進行運動學和動力學分析。這種方法的優點是建模方便，計算快速，理論分析方便。另一種方法是將部件看成沒有變形的剛體，并通過自由度耦合將各部件聯接起來，組成一個多剛體系統[6-7]。通過矢量力學或分析力學，建立起規則統一的動力學方程組并求解。

隨著氣門機構向輕量化和高速化方向的發展[8-10]，部件柔度和慣性力的加大，部件的彈性變形可能給氣門的運動輸出結果帶來誤差，因此有必要計入這種變形對計算精度的影響[11]。計入變形的影響后，還能夠在計算出動力學響應的同時得到部件的應力譜。

目前比較流行的包含彈性體多體系統的建模方法是相對坐標法。這種方法首先在彈性體上建立一個浮動坐標系，將彈性體的運動分解成浮動坐標系的牽連運動和其相對于浮動坐標系的變形運動的疊加。彈性變形則按無大范圍運動時理想邊界條件下的整體變形或模態函數進行離散，其中假設模態法在實際工程中得到了廣泛的應用，成為很多商業動力學軟件的理論基礎[12-14]。假設模態法將整個彈性體的變形場表示為一組模態函數和模態坐標的線性組合。因為計算之前需要先進行模態分析，所以只能考慮小變形的情況，而且模態函數的選取有一定的隨意性。

以Shabana[15]為代表的一派學者提倡將彈性體的大位移以及彈性變形都用相對總體坐標系的單元節點坐標描述，不區分物體的剛性運動和變形，對二者都按照連續介質力學的方式統一處理，在此基礎上發展出能處理部件在任意范圍內的大運動、大變形的非線性動力學模型，稱之為絕對坐標法。相對坐標法可以看作是多剛體動力學框架內的推廣。絕對坐標法實際上是一種有限元法，通過有限元網格離散的同時計算部件的剛體運動和變形（應力），而且可以直接計算彈性體之間的接觸－沖擊及滑動，這是假設模態法難以做到的。理論上非線性有限元的工具都有可能用來進行多體動力學分析。

有關用集總參數法或多剛體模型和假設模態法的剛－彈耦合模型的氣門機構動力學計算例子已經有很多，但應用非線性有限元法對氣門機構進行彈性動力分析尚少見相關的文獻。

ABAQUS是一款基于非線性連續介質理論的通用非線性有限元軟件，本文試圖在ABAQUS軟件的基礎上，建立起適合計算氣門機構的混合剛體－彈性體非線性有限元動力學模型。

2 模型處理

某發動機的氣門機構模型如圖1所示。有限元動力學模型的自由度比較多，如何正確反映模型的主要特點又盡可能地縮小模型的規模，是建模的首要考慮。比如剛度大的部件可以簡化成剛體，有限元模型盡量使用殼單元和梁單元，避免使用實體單元等。

圖1 總體幾何模型

2.1 凸輪

凸輪結構的剛度比較大，為降低模型規模，凸輪簡化成剛體，剛體的整體運動通過參考點控制（圖2）。參考點即為坐標原點（0，0，0）。由于假設凸輪繞參考點作等速轉動，不需要輸入轉動慣量。

圖2 凸輪剛體模型

2.2 滾子

滾子也即滾子軸承，整體剛度大，零件多，為了降低模型規模，滾子也簡化成剛體（圖3）。滾子的質量是16.8 g，對旋轉軸的轉動慣量是442 g·mm2。

圖3 滾子剛體模型

2.3 搖臂

搖臂是氣門機構系統中承受凸輪和彈簧動載荷的部件，剛度不高，應使用彈性體模型（圖4）。搖臂是鈑金件，適合用殼單元模擬，并可以方便地通過定義殼厚h將搖臂模型參數化，搖臂殼厚為3 mm。

圖4 搖臂有限元模型

搖臂使用通用縮減積分殼單元S4R，S4R單元可以根據定義的殼厚度自動選擇薄殼或厚殼公式。限于篇幅，S4R通用殼單元的有限元格式本文不再詳述，相關內容可以參考相關文獻[16]。

2.4 氣門與氣門彈簧

氣門雖然也有一定的變形，但本文氣門不是關注的重點，僅取與搖臂接觸的頂面做成剛體面，氣門簡化成參考點上的集中質量，如圖5所示。

圖5 氣門與氣門彈簧模型

除了氣門外，彈簧上座、氣門鎖夾的質量也都集中到氣門參考點上，三者質量之和為42.1 g。

彈簧對氣門機構的動力學響應有重要影響[17-18]，等距螺旋彈簧剛度可計算如下[19]，式（1）中各項參數見表1。

式中，

K——彈簧剛度，N/mm；

G——氣門彈簧剪切彈性模量，MPa；

d——鋼絲直徑，mm；

D——彈簧中徑，mm；

n——有效圈數。

表1 彈簧參數

根據式（1）計算得彈簧剛度K=39.5 N/mm，與搖臂類似，可以通過定義d將彈簧模型參數化。

為了模擬彈簧的顫振效果，將彈簧分成11段，質量點之間以彈簧－阻尼聯接，見圖5。

每點的質量可以計算如下[20]：

其中，ρ為彈簧密度（同表2），其余同式（1）。

取1.85Mi的彈簧質量集中到氣門參考點上（圖5），其余質量點均取Mi。

文獻[1]經試驗研究建議，等距螺旋彈簧阻尼比δ取0.01～0.011。本文取0.011，彈簧阻尼可計算如下：

K、M、c分別為彈簧的總剛度、總質量和總阻尼。由于各個質量點是串聯的關系，各質量點之間的彈簧阻尼要取總體剛度和阻尼的11倍。

3 模型設置

3.1 材料特性

氣門機構各部件均由鋼材組成。凸輪、滾子、彈簧、氣門分別簡化成剛體和集中質量，無需輸入彈性材料參數。搖臂材料是16MnCr5，材料的特性參數見表2。

表2 搖臂材料參數

3.2 接觸連接

凸輪－滾子（圖1）、搖臂－氣門（圖5）之間的接觸使用基于罰函數法的有限滑動、面－面接觸公式[21]，摩擦系數為0.01。液壓間隙調節器（HLA）不是關注的重點，可以簡化成固定的接地點。搖臂右側球面與接地點建立彈性連接（圖6），耦合全部自由度。彈性連接將一組節點與一個參考點相連，由參考點控制被連接節點的幾何運動。與剛性連接不同的是彈性連接允許被連接節點之間有相對位移，而且參考點所受載荷會根據被連接點的相對位置加權分配到各個節點上[21]。

圖6 搖臂與接地點的連接

由于搖臂簡化成殼單元，搖臂與滾子的接觸面（圖1）也退化成一圈節點（圖7），兩側分別與位于圓心的參考點建立彈性連接，耦合全部自由度。

圖7 搖臂與滾子的連接

ABAQUS的CONNECTOR單元提供了各種模擬2個點之間相對運動關系的方法；其中Hinge連接可以模擬2個點之間僅有一個自由度的相對轉動[21]，符合滾子與搖臂之間的運動關系。

3.3 載荷工況

求解分以下幾個步驟：

1）沿彈簧作用方向移動彈簧參考點（圖5），使彈簧產生一個初始預緊力，同時凸輪－滾子、搖臂－氣門建立接觸關系。

2）轉動凸輪參考點，帶動凸輪以及滾子搖臂工作，假設凸輪軸轉速是3 000 r/min（對應曲軸轉速6 000 r/min）等速旋轉。

3）凸輪工作時，固定凸輪參考點和搖臂接地點除轉動方向的所有自由度，固定彈簧參考點。

4 模型計算

4.1 計算方法

考慮無阻尼結構動力學方程

式中，

M——質量矩陣；

K——為剛度矩陣；

F——是外力向量；

U——是位移向量。

方程（4）的初始條件

對t+Δt時刻，將方程（4）、（5）表示為如下形式

以上是HHT隱式積分格式[22]，其中β=1/4 (1-α2)，γ=1/2-α，-1/2≤α≤0。HHT法是Newmark法的推廣，當α=0等同于Newmark法，二階精度，無條件穩定。HHT法的優點是通過參數α自動在方程（4）中引入一個數值阻尼，并可以通過α方便地控制其大小。

由于網格離散形成的非光滑表面，有限元動力學模型在計算接觸－沖擊的時候會帶來沒有意義的高頻響應。HHT法通過微小的數值阻尼就可以抑制高頻噪聲，而對低頻的求解沒什么影響，因此很適合用來求解有限元動力學問題。

4.2 氣門的動態響應

取α=-0.414 21，每周計算360步，時間步長為5.5×10-5s。氣門參考點的升程、速度、加速度曲線如圖8所示，凸輪-滾子和氣門-搖臂的接觸力如圖9所示。

圖8 氣門響應特性

由圖8可以看出，由于彈簧阻尼的作用，隨著凸輪升程的結束，速度、加速度都會迅速衰減至穩定狀態，符合實際的工作情形。由圖9可以看出，接觸力都有明顯的波動效應，升程結束后迅速衰減。如果接觸力小于0，氣門和搖臂之間會出現飛脫，凸輪在旋轉的同時氣門升程卻沒有變化，氣門不能正常工作，產生接觸損失。本文中接觸力在一個凸輪周期內都大于0。

圖9 凸輪－滾子和氣門－搖臂的接觸力

4.3 搖臂的強度計算

有限元彈性動力分析能在求解動力學積分的同時對包含單元的彈性體直接計算出每一步的應力，最大凸輪接觸力時刻應力分布如圖10所示。A是最大的應力集中點，在工作周期內A點的應力譜如圖11所示。

可見A點應力譜和凸輪接觸力非常相似，說明凸輪接觸力是影響搖臂強度最重要的因素。厚為3 mm的搖臂模型A點應力最大值為103 MPa。

圖10 搖臂應力分布

圖11 A點應力譜

凸輪－滾子的接觸應力Scon是氣門機構性能最重要的評價指標之一。2個外圓柱之間的線接觸可以按赫茲接觸公式計算[23]。

式中：

F——接觸力；

W——接觸線寬度；

ρ1，ρ2——接觸點處凸輪與滾子的曲率半徑；

μ1，μ2——凸輪與滾子的泊松比；

E1，E2——凸輪與滾子的彈性模量。

表3中ρ1是最大接觸力時刻凸輪的曲率半徑，可以通過三點定圓法計算得到。凸輪－滾子最大接觸力為1 514 N，根據式（9）相應的接觸應力為921.9 MPa。

根據設計規范（表4），搖臂強度和接觸應力滿足要求，搖臂強度良好。但接觸應力偏高，有必要對模型進行進一步優化。

表3 接觸輸入參數

表4 設計規范（MPa）

5 模型優化

5.1 設計方案

使用殼單元搖臂模型和多點質量彈簧模型的優點之一就是可以方便地通過定義相關幾何尺寸對模型參數化以便優化分析。選取搖臂厚度h和鋼絲直徑d作為設計變量，設計空間分別為h[2.1，3.0]和d[3.47，3.7]，厚度h取7個點，直徑d取6個點，在二維空間上均勻布置，共42個方案，任一方案均是參數（h，d）的組合。

根據前述有限元模型和參數設計方案，搖臂A處最大應力和凸輪－滾子接觸應力如表5和表6所示；搖臂A點應力和接觸應力的響應面如圖12所示。

表5 搖臂A點應力（MPa）

表6 接觸應力（MPa）

圖12 凸輪－滾子接觸應力

從圖表中可以看出無論是搖臂強度還是凸輪－滾子接觸強度的響應面都接近平面，沒有明顯的凸面或凹面，并且單調。

5.2 回歸優化分析

響應面法是基于回歸理論的優化方法。根據一定數量的測試結果，建立輸出參數和輸入變量的經驗方程，通過求解回歸方程達到優化設計的目的[24]。多項式響應面模型是應用比較廣泛的一種近似模型，通常使用二元二次模型[25]

實際應用中如果沒有明顯的凸面或凹面（圖12），可以略去式（10）中的高階項，簡化為平面模型

使用Excel軟件自帶的多元線性回歸工具分別對搖臂強度和接觸強度結果進行回歸分析[26]，得到搖臂強度Slever和接觸強度Scon的回歸方程

氣門機構優化的目的是降低凸輪－滾子之間的接觸強度使其進入“良好”區域，同時保證搖臂的強度在可接受區域內，也就是Slever≤125，Scon≤900。

將Slever=125，Scon=900代入式（12,13），可以得到搖臂厚度h=2.364，鋼絲直徑d=3.635。該計算結果只是由回歸方程得到的理論值，考慮到搖臂實際的加工制造和彈簧元件選型的限制，取h為2.4 mm，d為3.6 mm，代入式（12,13），得Slever=122.1 MPa，Scon=895.7 MPa，符合設計要求。

將優化后的參數（2.4,3.6）重新代入模型，計算得到A點應力為119.8 MPa，凸輪－滾子接觸應力為891.1 MPa，與前述回歸方程計算結果的誤差分別為1.88%和0.51%。凸輪－滾子接觸力、氣門加速度和搖臂－氣門接觸力對比如圖13至圖15所示。

圖13 滾子接觸力對比

圖14 氣門加速度

6 結論

針對某發動機氣門機構建立了非線性參數化有限元模型，進行了隱式動力學分析，得到了搖臂和氣門的強度及動力學響應，并對結果進行了優化分析。

圖15 氣門接觸力

經回歸優化設計后的氣門機構搖臂厚度由3 mm降到2.4 mm，減重達20%；凸輪－滾子間的接觸應力由921.9 MPa降到了891.1MPa，同時搖臂應力也滿足設計規范。接觸應力的降低有利于提高部件的疲勞耐久性和減少磨損。

由圖13及式（12）和式（13）可以看到，搖臂強度和接觸強度都是彈簧鋼絲直徑d的單調函數，d減少，搖臂強度和接觸強度也都會降低。在前述優化結果的基礎上，模型還有進一步優化的潛力，但前提是保證氣門與搖臂的接觸力大于0。

本例只是選取了搖臂在凸輪周期內的峰值應力校核搖臂強度。如果對搖臂應力譜用雨流計數法[27]之類的載荷統計工具，可以計算出搖臂在一個工作周期內所有的應力循環，結合疲勞分析理論得到更精確的強度壽命，這也是今后進一步開展工作的方向。

1尚漢冀．內燃機配氣凸輪機構：設計與計算[M].上海：復旦大學出版社，1988．

2劉忠民，俞小莉，沈瑜銘．配氣機構動力學模型的比較研究[J]．浙江大學學報，2005：39卷（12）．

3楊明軒.高速凸輪機構動力學模型研究[J]．機械，2007：（7）．

4 Husselman M．Modeling and Veriflication of Valve Train Dynamics in Engines[D]．Stellenbosch University，2005．

5 Desai H D，Patel V K．Computer Aided Kinematic and Dynamic Analysis of Cam and Follower[C]．London：Proceedings of the World Congress on Engineering，WCE 2010，2010：June 30-July 2．

6王曉，張保成．基于多體系統仿真的內燃機配氣機構動力學分析[J]．內燃機，2008（1）．

7 Cofaru C，Sandu G，Jelenschi L，Multibodyanalysis of the finger follower valvetrain system[J]．Transilvania University of Brasov,2010．

8倪計民．內燃機原理[M]．上海：同濟大學出版社，1997．

9張正有．發動機配氣機構發展綜述[N]．重慶工學院汽車學院，2002．

10路瓊瓊，李智，雷晶．內燃機配氣機構技術現狀及發展[J]．機械，2009年，36（4）．

11張策．機械動力學[M]．北京：高等教育出版社，2000．

12赫永剛．基于多體系統的配氣機構動力學特性研究[D].浙江大學機械與能源工程學院，2006．

13楊輝．剛－柔耦合動力學系統的建模理論與實驗研究[D]．上海交通大學，2002．

14 Dalpiaz G，Rivola A．A model for the elasto-dynamic analysis of a desmodromic valve train,tenth world congress on the therory of machines and mechanics[J]．Finland：Oulu,1999：June 20-24．

15 Shabana A A.Computational continuummechanics[M]. Cambridge University Press，2008．

16 Abaqus Theory Manual.Dassault Système,2011．

17傅紅良，林運．氣門彈簧特性對配氣機構的影響[J]．柴油機設計與制造，2007（3）.

18徐小明，蔣炎坤，劉志恩．基于分離變量法的發動機配氣機構彈簧動力學模擬[J]．柴油機，2007：（29）4．

19季文美等．機械振動[M]．北京：科學出版社，1985．

20 Seidlitz S．Valve train dynamics-A computer study[C]．SAE 890620．

21 Abaqus Analysis Manual.Dassault Système,2011．

22 Hilber H M，Hughes T J R，Taylor R L．Improved Numerical Dissipation for Time Integration Algorithms in Structural Dynamics[J]．Earthquake engineering and structural dynamics，Vol.5,283-292,1977．

23Johnson K L著．接觸力學[M]，徐秉業等譯．北京：高等教育出版社，1992．

24貴勁松，康海貴．結構可靠度分析的響應面法及其MATLAB實現[J]．計算力學學報，2004，21（6）．

25付凌暉，王惠文．多項式回歸的建模方法比較研究[J]．數理統計與管理，2004，23（1）．

26艾自勝等．Excel在曲線擬合中的應用[J]．蘇州大學學報（醫學版），2008，28（5）．

27徐灝．疲勞強度[M]．北京：高等教育出版社，1988．

The Application of FEM Elastic-dynamic Analysis and Optimization in a Valve Train

Sun Xianguo1,2,Liu Hao2,Ni Jimin1
（1 Automobile Academy,Tongji University;2 Engine Dept.,Schaeffler Group,China）

An elastic-dynamic analysis of a valve train is carried with parametric nonlinear rigid-flexible FEM mode,and then the mode is optimized by regression equation.As a result,the mass of the lever is reduced,and the contact performance between cam and roller is improved when the strength of the lever is safe at the same time,so the fatigue durability is improved and the wear of cam and roller decreased.

valve train,FEM,elastic-dynamic,regression optimization

10.3969/j.issn.1671-0614.2012.04.002

來稿日期：2012-08-07

孫先國（1981-），男，高級工程師，主要研究方向為機械結構有限元分析。