發(fā)動機配氣機構運動學及動力學分析-重慶大學車輛工程專業(yè)本科學生畢業(yè)設計

1、<p>　　重慶大學本科學生畢業(yè)設計（論文）</p><p>　　發(fā)動機配氣機構運動學及動力學分析</p><p><b>　　學 生：黎明</b></p><p>　　學 號：20133790</p><p>　　指導教師：阮登芳（教授）</p><p>　　專 業(yè)：

2、車輛工程</p><p>　　重慶大學車輛工程學院</p><p><b>　　二零一七年五月</b></p><p>　　Graduation Design(Thesis) of Chongqing University</p><p>　　Kinematics and dynamics analysis for en

3、gine valve train</p><p>　　Undergraduate: Li Ming</p><p>　　Supervisor: Prof. Ruan Dengfang</p><p>　　Major: Vehicle Engineering</p><p>　　College of Vehicle Engineering<

4、;/p><p>　　Chongqing University</p><p><b>　　May 2017</b></p><p><b>　　摘 要</b></p><p>　　配氣機構是發(fā)動機的重要組成部分，其設計的合理與否直接影響到發(fā)動機的充氣效率以及換氣質量，因此對發(fā)動機的動力性、燃油經(jīng)濟

5、性、可靠性、有害物質排放、發(fā)動機噪聲和振動有較大的影響[]。而頂置凸輪軸式配氣機構由于能適應更高的轉速而在許多小型汽油機中廣泛使用。但是頂置凸輪軸由于搖臂傳動比是變值，所以其幾何關系要復雜很多[]。本文在已知凸輪對搖臂的運動規(guī)律的條件下，針對某125發(fā)動機的配氣機構，經(jīng)理論分析運動學規(guī)律，并用matlab計算出其氣門對轉角的理論升程、速度、加速度。在考慮氣門間隙及傳動機構變形的影響下，建立配氣機構運動的單自由度模型，得出運動二階微分方程

6、。利用matlab采用龍格——庫塔法計算出氣門的實際運動規(guī)律，對比氣門實際升程和理論升程，對該發(fā)動機配氣機構的 “飛脫”、“反跳”以及運轉的平穩(wěn)性進行動力學特性評價。從而完成了整個配氣機構的運動學及動力學計算。</p><p>　　關鍵詞：運動學，動力學，配氣機構，matlab，龍格庫塔法</p><p><b>　　ABSTRACT</b></p>&

7、lt;p>　　Valve train is an important part of the engine, which has directly affect to the engine's volumetric efficiency and the quality of ventilation, so there is also a greater influence to the engine power, fuel

8、 economy, reliability, emissions of harmful substances, engine noise and vibration. Because the overhead camshaft valve train can adapt to the higher speed, it is widely used in many small gasoline engine. But for the o

9、verhead camshaft, the drive ratio of the rocker is changed by the time, so</p><p>　　Key words： Valve train, Kinematics, Dynamics,Matlab</p><p><b>　　一、緒論</b></p><p>　　本課題

10、以某125型摩托車發(fā)動機的頂置凸輪式配氣機構為研究對象，分別對其進行了運動學分析、剛度計算、以及動力學分析，并由所得到的數(shù)據(jù)對該機構進行動力學評估，為該發(fā)動機配氣機構的合理設計奠定基礎。</p><p><b>　　1.1課題研究意義</b></p><p>　　配氣機構是發(fā)動機結構中最復雜、工作最繁重的部件之一，固然也是發(fā)動機的重要組成部分，配氣機構控制著進排氣門的

11、開啟時刻和開啟速度，根據(jù)每個氣缸的工作順序，控制進氣和排氣，實現(xiàn)換氣功能。</p><p>　　配氣機構可以視為兩部分組成，即氣門組和氣門傳動組。根據(jù)凸輪軸的位置不同，將配氣機構分為三類：下置凸輪軸式配氣機構、中置凸輪軸式配氣機構以及頂置凸輪軸式配氣機構。本課題便是以凸輪軸頂置式配氣機構為研究對象展開工作。配氣機構的運轉平穩(wěn)性一部分決定了發(fā)動機振顫程度的大小，已成為一項非常重要的設計參考性能指標[]。對配氣機構的

12、研究內容歸納起來主要有兩個方面，一方面是零部件的設計，另一方面是機構的動力學問題。零部件設計又可以簡單分成三類：氣門及氣門等零部件設，計凸輪型線優(yōu)化設計，搖臂機構運動學設計。因為凸輪作為整個機構的原動件，它直接控制整個機構的運動，所以凸輪型線的設計至關重要。而對于機構動力性能的研究，又主要表現(xiàn)在凸輪型線如何如何影響氣門的運動規(guī)律。配氣機構的設計最初為剛性設計，但沒有絕對剛性的零部件，在機械應力以及熱應力作用下總會存在變形，所以逐步發(fā)展為

13、彈性設計。</p><p>　　高速汽油機由于轉速較快，每分鐘轉速可達5000轉以上，為保證充氣效率，都采用頂置式氣門裝置。這種裝置都適合用凸輪軸的三種安裝形式，但是，如果采用下置式或者中置式的凸輪軸，由于氣門與凸輪軸的距離較遠，需要氣門挺桿和挺柱等輔助零件，這樣的結構勢必造成氣門傳動件較多，結構復雜，發(fā)動機體積也不可避免的變得很大，而且在高速運轉下還容易產生噪聲，運轉平穩(wěn)性較差。所以，現(xiàn)代轎車發(fā)動機一般都采用頂

14、置式凸輪軸，以期待能改變這種現(xiàn)象。將凸輪軸配置頂置靠近氣門，不僅縮短了凸輪軸與氣門之間的距離，而且省略了挺桿和挺柱，從而簡化了傳動機構，使得發(fā)動機的結構變得緊湊。更重要的是，這種凸輪軸頂置的安裝方式可以減少整個系統(tǒng)往復運動的質量。在國內的，不少小型汽油機已采用了凸輪軸頂置的配氣機構。頂置凸輪軸配氣機構的系統(tǒng)剛度雖然大，但實際上仍是彈性系統(tǒng)，可看作一個“彈簧——質量系統(tǒng)”[]。當氣門機構運轉工作的時候，機構會產生彈性變形，再加上熱應力的影

15、響會使氣門運動產生很大程度的變化，與理論分析脫離，使傳動鏈脫節(jié)，將出現(xiàn)一系列不利于發(fā)動機運轉的現(xiàn)象，諸如氣門的開閉不正常，振動和噪聲增大、零件易損壞等現(xiàn)象，從而限制了轉速的提高</p><p>　　1.2課題國內外研究狀況</p><p>　　1.2.1國外研究現(xiàn)狀</p><p>　　處在21世紀的今天，計算機能力已然十分強大，尤其是在計算分析領域。越來越多設計過

16、程中運用CAE仿真技術，仿真結果相比理論分析計算結果也能更加有效的預測配氣機構的性能。Seon-Yang Hwang等人發(fā)現(xiàn)氣門的沖擊噪聲主要是脈沖噪聲，他們分別對兩種不同挺柱的配氣機構進行噪聲分析，得出結論：既要保證發(fā)動機性能又同時達到降低噪聲的目的，可以通過優(yōu)化緩沖階段的凸輪型線來實現(xiàn)。Nitin P.Gokhale運用三維CFD分析氣門座溫度與潤滑液流速的關系，氣門座溫度升高，潤滑液流速降低，兩者成負相關。與此同時還對配氣機構的磨

17、損做了深入研究，發(fā)現(xiàn)磨損主要集中在搖臂頭氣門小端的磨損，凸輪-挺柱摩擦磨損和氣門-氣門座沖擊磨損，這三者占整個磨損的60%~70%。A.K.Jamkhande在凸輪型線設計研究中采用了B-樣條曲線，設計的凸輪型線不僅更加平滑而且有更低的越度，加速度曲線上不光順的地方明顯減少，在凸輪尖端具有更小負加速度。R SPrabakar和S PMangalaramanan將氣門彈簧視為彈性體分析，運用ADAMS發(fā)現(xiàn)氣門反跳一般不出現(xiàn)在低速發(fā)動機上，

18、并且對于高速汽油機，也只集中出現(xiàn)在極速轉</p><p>　　1.2.2國內研究現(xiàn)狀</p><p>　　由于我國近代科學技術起步較晚，導師我國對發(fā)動機配氣機構的研究也相對晚了許多，樊久銘等人針對裝有液壓調節(jié)器的CA488汽油配氣機構，提出了“4+N”多質量動力模型。劉忠民等研發(fā)了可用于配氣機構的動力學特性和耐久性研究，可以控制配氣機構溫度和潤滑條件的綜合試驗系統(tǒng)。李興蘭利用ADAMS軟件

19、構建氣門機構的動態(tài)模型，獲得振動和噪聲分析與預測的主要部件之間的相互作用，提供準確的數(shù)據(jù)和邊界條件。曹志芬將氣門機構的剛性以簡單的機械方式進行測量，使配氣機構彈性系統(tǒng)模型的動力學分析有了最基本的剛度參數(shù)?；贑AD / CAE綜合開發(fā)技術，朱文波等探討發(fā)動機氣門機構數(shù)值研究的思路和過程，為氣門機構設計，動力學分析和優(yōu)化提供參考，經(jīng)驗設計以定量形式出現(xiàn)，產品開發(fā)效率高。王曉輝等人采用反向優(yōu)化設計方法重新設計凸輪輪廓，首先用測量設備測量實際

20、輪廓，然后使用三維軟件來適應凸輪輪廓，進行擬合，然后進入AVL軟件優(yōu)化，既確保了配氣機構的基本性能，又提高了凸輪型線的連續(xù)性和加工性。白賀研究了凸輪軸扭轉振動對氣門機構的影響。 通過扭轉振動實驗提取凸輪軸的扭轉角，并測量凸輪軸的動態(tài)轉矩。結果表明，凸輪軸的扭轉振動隨著轉速的</p><p><b>　　1.3課題研究背景</b></p><p>　　發(fā)動機的充氣效率以

21、及換氣質量很大程度上受配氣機構的設計影響，而發(fā)動機的動力性、燃油經(jīng)濟性、有害物質排放等各方面性能業(yè)余配氣機構的設計有關。隨著國際環(huán)境對發(fā)動機性能和排放要求不斷提高，配氣機構的研究也越來越深入，尤其是在彈性系統(tǒng)下的動力學特性研究更加迫切。</p><p>　　隨著制造商對發(fā)動機配氣機構的進一步設計改進，目前4氣門發(fā)動機在市場上已經(jīng)越來越普遍。大功率低燃消耗從來都是車主最為在意的兩個點，并且也是設計思想中起主導作用的

22、兩大因素。近來可變氣門正時及升程控制技術在解決低油耗問題上，變現(xiàn)十分突出，解決了由單個氣門控制氣體進出燃油供給時，充氣效率低，燃油時而供給不足，時而供給過量，從而造成發(fā)動機功率不夠，排放嚴重等問題，使發(fā)動機在更加省油的同時還能在全段轉速輸出期間都更有力，動力更強勁，因此頗為受廣大車主的青睞。凸輪軸頂置使配氣機構凸輪軸在旋轉中的負荷相應減小，并且對于凸輪軸和氣門彈簧的要求也降到了最低。這是因為它在設計上沒有挺柱、搖臂和推桿，直接通過凸輪軸

23、上的凸輪來驅動氣門開閉，這使得該機構在結構上有了大大簡化。同時從維修角度來看，這也降低了成本。因此頂置凸輪軸式配氣機構從此進入設計者眼中最關注的對象，由于頂置凸輪軸式配氣機構優(yōu)越的性能，簡介的傳動鏈結構，使得發(fā)動機體型有所減小，市場上各類型的發(fā)動機對這種配氣機構的采用率也越來越高。另外，從物理性能的角度來看，凸輪軸架空設計思路進入和排氣通道轉小，氣流阻力小，使氣體出來更加平滑。因此隨著內燃機結構上的</p><p&g

24、t;<b>　　1.4課題研究內容</b></p><p>　　1．配氣機構運動學分析。在傳動鏈結構上進行發(fā)動機運動學理論分析，通過matlab編程分析計算，得到分析運動學的一種通用算法。研究配氣機構理論運動規(guī)律，編寫動力matlab學分析程序，獲取氣門理論升程曲線，速度曲線，加速度曲線，提供動力學分析原始數(shù)據(jù)。</p><p>　　2．配氣機構剛度計算。利用Cati

25、a為配氣機構建模，并導入Workbench中進行剛度計算，提供動力學分析分析所需數(shù)據(jù)。</p><p>　　3．配氣機構動力學計算。采用龍格庫塔法，利用Matlab編寫變搖臂比頂置凸輪軸配氣機構的動力學計算通用程序，并得出氣門實際升程曲線，速度曲線以及加速度曲線。</p><p>　　4．對所選發(fā)動機配氣機構的動力學特性進行評價，</p><p>　　二、氣門機構的

26、主要設計要求</p><p>　　發(fā)動機中的運行必須有著新氣不斷進入，廢棄不斷排除的過程。較早采用的控制其進氣、排氣的方式有很多，如滑閥控制（回流掃氣二沖程機）、旋轉閥控制（有過個別試驗機型）、氣門控制、混合控制（有排氣門的直流掃氣二沖程機）等等。對于當今市場主流四沖程發(fā)動機而言，氣門控制進排氣的方式已經(jīng)逐漸穩(wěn)定下來，因其可靠的密封性能，便拆裝的優(yōu)點，氣門已經(jīng)在配氣機構中占據(jù)著不可替代的地位?，F(xiàn)有產品幾乎全都采用

27、氣門控制進排氣。</p><p>　　氣門機構控制著換氣過程，其設計好壞直接影響著四沖程發(fā)動機方方面面的性能。</p><p>　　首先，由內燃機原理闡述，發(fā)動機的最大功率轉速nmax正比于進氣門通路面積Av除以活塞排量Vh，即nmax（Av/Vh），而升功率NlpeNnmax。拿同一類型的發(fā)動機作比較，因最大功率工況的peN差不多一樣大，可以說他們的升功率NmaxnmaxAv/Vh。&l

28、t;/p><p>　　其次，進氣門定時影響充氣系數(shù)隨轉速變化的情況，影響泵氣損失，影響換氣質量，因此對發(fā)動機的動力性，燃料經(jīng)濟性和有害排放有影響。此外，配氣機構的因運動副的摩擦而發(fā)生摩擦損失，與發(fā)動機的轉速有關，總的說來發(fā)動機轉速越低摩擦損失的比重越大，轉速越高，配氣機構摩擦損失所占比重越低。與此同時發(fā)動機的機械效率和油耗率也隨著配氣機構的機構摩擦損失與發(fā)動機中低轉速負荷工況下變化而變化。</p>&l

29、t;p>　　配氣機構各個工況轉速下所產生的噪聲以及工作不平穩(wěn)的振動很大程度影響整個發(fā)動機可靠性與否和噪聲水平高低。在高熱負荷下，由于材料性能的限制氣門將發(fā)生熱裂、凸輪在高速運轉下受到較高接觸應力的影響發(fā)生接觸疲勞、在各種應力綜合作用下氣門彈簧也可能發(fā)生疲勞失效。諸如此類都是重要的可靠性問題，在設計過程中不得不引起重視，對動力學特性做出做夠有效的分析評價。氣門機構是彈性系統(tǒng)，在工作過程中一旦發(fā)生振動，就有可能產生“飛脫”，“反跳”，

30、氣門“提前落座”一系列異常現(xiàn)象，并且產生這種異?，F(xiàn)象的概率會隨著轉速的提高而增加，這些等異?，F(xiàn)象的出現(xiàn)，將使得凸輪-從動件接觸面的磨損加劇，從而影響著發(fā)動機的壽命和運轉平穩(wěn)性以及工作可靠性。</p><p>　　根據(jù)這些情況，對氣門機構的設計理念和主要要求可總結為：</p><p>　　為使充氣量足夠大，在氣門口要保證足夠大的氣流通路面積；</p><p>　　機構

31、的動力學特應當滿足機構的“飛脫”和“反跳”不出現(xiàn)在任意轉速工況下；</p><p>　　為了抑制機構接觸表面的磨損現(xiàn)象，不應該有過高的熱應力和機械應力發(fā)生在重要零件表面，同時在材料方面也必須考慮到使用性能上的耐久性作；</p><p>　　便于制造和維修，成本較低。</p><p><b>　　三、運動學分析</b></p>&l

32、t;p>　　配氣機構運動學所研究的內容包括：①凸輪從動件（下置凸輪軸式氣門機構中的挺柱，頂置凸輪軸式的搖臂）運動規(guī)律與凸輪輪廓開關系，以及從動件與凸輪的相對滑動關系（本文不研究相對滑動）；②在不考慮機構彈性的條件下，凸輪從動件的運動規(guī)律與氣門運動規(guī)律的關系。</p><p>　　3.1凸輪廓線預處理</p><p>　　對于普通的下置凸輪式配氣機構的運動學計算，只要乘以或除以固定的搖

33、臂比即可，但是針對本課題研究的配氣機構采用了凸輪——搖臂——氣門如圖1所示這種頂置凸輪式配氣機構，這種配氣機構的搖臂比是隨著凸輪轉角變化而變化的，要通過配氣機構運動件的幾何關系及其運動規(guī)律利用三角函數(shù)關系換算得到。</p><p>　　對于經(jīng)搖臂驅動氣門的機構來說，無論搖臂是處于搖臂中間還是處于搖臂一端，運動學分析方法是一樣的。頂置凸輪軸的凸輪輪廓通常是按照它驅動平底挺柱時的挺柱升程規(guī)律給定、制造并用平底測頭檢驗

34、的。運動學分析要解決由已知的求不同凸輪轉角下?lián)u臂擺角，再進而求氣門理論運動規(guī)律的問題。</p><p>　　初始獲得的凸輪廓線數(shù)據(jù)是不光順的，顯然不符合設計理念要求。將數(shù)據(jù)導入matlab分析計算之前需要將給定的凸輪廓線數(shù)據(jù)在matlab中進行光順預處理，利用matlab分析軟件強大的數(shù)據(jù)插值，曲線擬合的功能，將凸輪廓線首先處理成成一條光滑的曲線，但是這時并沒有確切的升程與凸輪轉角的相應的函數(shù)關系式，只用采用樣條

35、插值的方法得到一組數(shù)據(jù)。導入excel制作成凸輪升程一覽表，如附錄表C，為matlab編程計算做準備。</p><p>　　3.2氣門理論運動規(guī)律與凸輪輪廓的關系</p><p>　　圖2給出了頂置凸輪軸及搖臂從動件的關系圖。圖中是搖臂柱面中心A在凸輪以頂點D與搖臂接觸時的位置。任一瞬時凸輪的轉角以與的夾角表示，搖臂的位置則以與的夾角W表示。圖中還給出了凸輪上同一點F與平底挺柱接觸時的升程

36、y和轉角α，以及半徑為的球面挺柱接觸時的升程和凸輪轉角。本次研究的嘉陵125發(fā)動機配氣機構，是一個朝樞軸轉的凸輪，規(guī)定以氣門全開時的轉角為00，在凸輪上升段轉角為負，下降段轉角為正。</p><p>　　在得到凸輪升程規(guī)律的基礎上，計算氣門運動規(guī)律是分析的基礎。首先從每一個（y，α）數(shù)據(jù)求對應的（W，）數(shù)據(jù)。在△中，利用三角函數(shù)關系可得：</p><p><b>　?。?.1）&

37、lt;/b></p><p>　　式中、和R的含義見圖2，針對本課題所分析機構，通過搖臂尺寸相關數(shù)據(jù)計算得知,，而且由△知:</p><p><b>　?。?.2）</b></p><p>　　其中凸輪半徑，。代表著升程對轉角的一階導數(shù)，，單位（mm/弧度）。</p><p>　　凸輪轉角與，α的關系為：</

38、p><p><b>　?。?.3）</b></p><p>　　ε、K和的計算式分如下：</p><p><b>　?。?.4）</b></p><p><b>　?。?.5）</b></p><p><b>　　（3.6）</b>&l

39、t;/p><p><b>　　式中：</b></p><p><b>　?。?.7）</b></p><p>　　接著再由（W，）求對應的氣門升程h。由圖3可見，當搖臂與凸輪基圓接觸時，搖臂長臂與氣門桿端面間有間隙Δ=0.05mm,而氣門桿端面比搖臂中心高出z=1.65mm，這時搖臂的位置角為：</p><

40、p><b>　?。?.8）</b></p><p><b>　?。?.9）</b></p><p>　　其中調節(jié)螺釘球半徑，經(jīng)過計算得，；與凸輪轉角對應的β和氣門理論升程h為：</p><p><b>　?。?.10）</b></p><p><b>　?。?.

41、11）</b></p><p>　　至此，由每一組（y，α）計算對應的（W，）再計算對應的（h，）的過程就完成了。但是，等步長遞增的α所對應的是不按等步長遞增的。在進行動力學分析時，必須要有等步長遞增的對應氣門理論升程h數(shù)據(jù)。因此，對于理論分析運算結果，還必需進行擬合與插值處理才能得到所要求的升程—轉角數(shù)據(jù)[]。</p><p>　　3.3運動學理論分析后的計算結果</p

42、><p>　　將公式（3.1）到（3.11）用matlab編程計算得到如圖4理論氣門升程曲線，matlab程序參見附錄1。</p><p>　　圖4.氣門理論升程曲線</p><p>　　并利用matlab樣條插值功能得到等步長遞增的所對應的氣門升程。擬合，插值處理后得到氣門的速度，加速度如圖5，其中速度與加速度值以轉速7000r/min換算得出。</p>

43、<p>　　圖5氣門速度以及加速度曲線</p><p>　　至此，配氣機構的運動學分析就已全部結束，氣門的最大升程當，hmax=7.815mm；氣門最大正速度，y’max=7.855m/s；氣門最大負速度，y’max=-14.47m/s。氣門加速度情況略微復雜，上升段最大正加速度y”max=1713m/s2，下降段最大正加速度y”max=11500m/s2，最大負加速度y”max=-11600m/s2

44、。加速度的劇烈波動是由速度曲線上的一些尖點造成，消除速度波動尖點，光順速度曲線逆向優(yōu)化設計凸輪廓線可以解決這一問題。</p><p>　　若用高階多項式擬合所得到的氣門理論升程曲線，并再求導，得出氣門理論速度以及加速度，將有如下結果：</p><p>　　圖6光順處理后的速度及加速度曲線</p><p>　　通過擬合前后對比，可以清晰看出，若凸輪型線數(shù)據(jù)給的足夠光滑

45、，我們所得結果也將會是光滑的，并且從加速度曲線上面可以看出，加速度最值大大降低，并且加速度沒有突變，氣門的運動不發(fā)生剛性沖擊。所以在配氣機構設計中，需要我們逆向設計優(yōu)化凸輪型線。</p><p><b>　　四、動力學分析</b></p><p>　　4.1動力學理論分析</p><p>　　實際氣門機構是一個彈性系統(tǒng)，配氣機構動力學正是要研究

46、考慮彈性的情況下的氣門的運動規(guī)律。運動學分析中得出的氣門理論運動規(guī)律，這為進一步做動力學分析提供了可能，理論氣門升程像是一個激振力函數(shù)，而真實情況下氣門的運動歸路則是對應這種激勵下的響應。由于配氣機構是有彈性的，相當于彈簧——質量系統(tǒng)，氣門的實際運動規(guī)律會偏離運動學分析所得的理論運動規(guī)律[]。由于氣門存在彈簧預緊力PS0，在剛消除氣門間隙的瞬時氣門還不會打開，要等氣門機構發(fā)生一定變形而變形力PS0時才能頂開氣門。因此氣門的實際運動要比氣

47、門的理論運動始點晚一些。</p><p>　　發(fā)動機在高速運行時，因彈性振動的影響，氣門的實際運動情況與理論規(guī)律差別較大，甚至有時氣門實際升程大于理論升程，這時由凸輪至氣門的整個傳動鏈將出現(xiàn)脫節(jié)之處，稱作“飛脫”；另一方面，氣門落座時刻早于理論落座點，以高速沖擊氣門座并可能“反跳”。正確設計的配氣機構不容許在發(fā)動機工作轉速范圍內出現(xiàn)飛脫和反跳。為此配氣機構的動力學特性需要進行專業(yè)的校核。</p>&

48、lt;p>　　為分析實際氣門運動規(guī)律及相應的各構件受力情況，需要建立動力學分析模型，寫出描述其運動情況的微分方程式然后求解。求解結果的可信度與模型本身有關，即取決于該模型能在多大程度上反映機構的特點。</p><p>　　已提出的氣門機構動力學計算模型有許多種，如單自由度模型、多自由度模型、有限元模型等等，其中單自由度模型是最基本而又最簡單的一種模型。</p><p>　　圖9中給出

49、了單自由度模型。其中是氣門組建和彈簧的當量質量，是搖臂、推桿、挺柱等驅動件的當量質量；是氣門彈簧剛度，k是從搖臂長臂端到凸輪軸為止的整個氣門機構驅動部分的剛度，是氣門座的剛度；b是驅動部分的當量內阻尼系數(shù)，是氣門座阻尼系數(shù)，是反映氣門組件粘性外阻尼及氣門彈簧內阻尼的當量阻尼系數(shù)，是驅動部分的當量外阻尼系數(shù)。</p><p>　　在模型中當量凸輪繞旋轉，驅使從動件運動，我們認為轉軸軸心是剛性的，從動件質量為0。從動

50、件的位移Y在動力學分析之前已經(jīng)從運動學分析中得到力理論運動規(guī)律，是由實際凸輪型線根據(jù)氣門機構傳動鏈的幾何關系換算到氣門理論升程[]。計算y的起點（氣門的初始位置）定在氣門座被壓縮了之處。質量的位移y相當于氣門的實際升程。</p><p>　　模型中的運動質量通常取為：</p><p><b>　　（4.1）</b></p><p><b&

51、gt;　?。?.2）</b></p><p>　　為搖臂的轉動慣量，用solidworks分析得到（kg·m^2），對于變搖臂比的氣門機構就取其變化范圍內的中值。經(jīng)計算i取i=2.1。是氣門、氣門彈簧座和氣門鎖夾這些氣門組件的當量質量，代表氣門彈簧的質量，其中，，取為：</p><p>　　其中，是氣門全開時的搖臂位置角。達倫伯原理指出，作用于運動質量上的力與慣性力符

52、合力平衡關系，故對于消除氣門座變形階段的運動情況應有：</p><p><b>　　整理后得到:</b></p><p><b>　?。?.3）</b></p><p>　　式中M表示此階段參與運動的質量，考慮到，通常都用位移對轉角的倒數(shù),來表示，即,,,所以上式可寫成：</p><p><b

53、>　　（4.4）</b></p><p>　　式中： </p><p><b>　?。?.5）</b></p><p>　　相當于激振力，是“激振力函數(shù)”。它只決定于，即只決定于凸輪輪廓和氣門機構傳動鏈的設計。其中是氣門彈簧預緊力，其值為，是凸輪軸轉速，本次研究針對2個轉速的情況作了運動學分析，分別是

54、本發(fā)動機的額定轉速7000r/min，高速轉速為12000r/min.</p><p>　　上述的微分方程可以通過matlab，龍格庫塔法迭代求解。</p><p><b>　　4.2搖臂比i</b></p><p>　　此搖臂比對前述運動學分析毫無意義，但在分析搖臂受力時需了解氣門對搖臂的作用力PvR與凸輪對搖臂的作用力PcR之比。PvR作用

55、于搖臂和氣門桿端面的接觸點C并指向B，PcR作用于搖臂和凸輪的接觸點并指向A（見圖2，圖3）。若忽略搖臂的慣性力矩不計，可得搖臂比公式如下：</p><p><b>　?。?.6）</b></p><p>　　其中J就是圖2中與的夾角:</p><p><b>　?。?.7）</b></p><p>

56、;　　若R，和是對應等差遞增求出來的，則所得的J和i就對應不等差遞增的。在動力學分析時，對于變搖臂比的氣門機構就取其變化范圍內的中值，</p><p><b>　　圖8搖臂比計算結果</b></p><p>　　經(jīng)matlab計算，取i=1.6。</p><p><b>　　4.3搖臂剛度計算</b></p>

57、<p>　　一個機構的剛度是指彈性體抵抗變形（彎曲、拉伸、壓縮等）的能力。頂置凸輪軸式結構，結構緊湊質量輕，體積小，但整個系統(tǒng)卻具有較高的剛度，因此頂置凸輪軸式配氣機構在很高轉速下工況下同樣能正常工作。因此對于頂置凸輪軸式配氣機構，其很重要的一個優(yōu)點就是機構剛度大，因此機構剛度的計算對本課題的研究又很重要的意義。圖5所示為機構剛度測量的理論簡圖。</p><p>　　由于氣門與氣門座剛度對機構剛度的影

58、響很小，對動力學分析的影響也很小，因此本次分析假設氣門與氣門座為剛體，只考慮了彈簧和搖臂的剛度。而彈簧的剛度已知，因此對機構搖臂剛度的計算就成為了該配氣機構剛度計算的主要任務。</p><p>　　利用catia完成搖臂3D建模: </p><p>　　圖10 Catia3D搖臂模型</p><p>　　將其導入workbench分析軟件中，搖臂材料選用合金鑄鐵，如

59、圖5施加約束，劃分75209個單元網(wǎng)格：</p><p>　　圖11.劃分單元網(wǎng)格</p><p>　　并施加P=1000N作用力，在這個作用力下測得搖臂總變形最大為0.085mm。</p><p><b>　　圖12.搖臂應變</b></p><p><b>　　所以搖臂的剛度：</b></

60、p><p>　　k=1000/0.085=1.17×104N/mm</p><p>　　4.4解動力學微分方程</p><p>　　下面的式子即為需要求解的方程：</p><p><b>　　（4.8）</b></p><p>　　即二階常微分方程。因為在求解區(qū)間的起點，即凸輪轉角為而氣門間

61、隙剛剛消除(Y=0)的時刻，y和是已知的，即</p><p>　　采用龍格庫塔法從這些初值開始逐步迭代，可得出該微分方程式在一系列計算點上的數(shù)值解。</p><p>　　設原方程為,并已知某一轉角為時的升程和速度各為，。取轉角步長為Δα=0.5°（一般取為0.5o~1.0o），則下一個計算點上的升程和速度可由和用下式算出，即</p><p><b&g

62、t;　　（4.9）</b></p><p><b>　?。?.9）</b></p><p>　　式中的 ——即 （4.10）</p><p><b>　　（4.11）</b></p><p><b>　?。?.12）</b>&

63、lt;/p><p><b>　?。?.13）</b></p><p>　　這樣從n=0初始點位置開始逐點迭代計算。每一下的升程和速度得出后，該點的加速度也就得出了。</p><p>　　微分方程中阻尼系數(shù)b、、和，在計算前必須獲得，但從以前的相關研究來看這些系數(shù)目前還無法分別測定。但是阻尼系數(shù)是可以推算的。已經(jīng)有實踐表明，阻尼值對整個動力學計算結果

64、不敏感。因此在此次計算中取。氣門座剛度很大，為了簡化計算，計算時不考慮氣門座的變形。將公式（4.4）、（4.5）（4.9）~(4.13)用matlab編程。</p><p>　　龍各庫塔法要求迭代步長為0.5°，按照k1、k2、k3、k4的計算要求，需要我們找出轉角增長0.25°倍數(shù)的轉角所對應的左右理論升程Y的值，在matlab中樣條插值每隔0.25°計算一次，便可得到所需Y值。&

65、lt;/p><p>　　圖13.計算按步長增長的Y值</p><p>　　詳細matlab程序過程參見附錄B。</p><p>　　4.5動力學分析結果</p><p>　　龍哥庫塔法matlab分析程序見附錄B，7000r/min額定轉速下動力學分析計算結果：</p><p>　　圖 14. 7000r/min轉速下實際

66、升程與理論升程對比</p><p>　　圖15.7000r/min轉速下速度曲線</p><p>　　圖16.7000r/min轉速下加速度曲線</p><p>　　最大升程為ymax=7.785mm，最大正速度Vmax=7.893m/s，最大負速度Vmax=-18.51m/s。加速度局部有較大波動，由速度在該區(qū)域波動引起，若解決該波動，加速度峰值勢必會極大下降，所

67、以本次分析中，加速度的極值具有很少的參考性。</p><p>　　12000r/min高速轉速下實際生程與理論升程對比：</p><p>　　圖17.12000r/min轉速下理論升程與實際升程對比</p><p>　　圖18.12000r/min轉速下速度曲線</p><p>　　圖19.12000r/min轉速下加速度曲線</p&g

68、t;<p>　　最大升程為ymax=7.774mm，最大正速度Vmax=7.886m/s，最大負速度Vmax=-20.85m/s。加速度局部有較大波動，由速度在該區(qū)域波動引起，若解決該波動，加速度峰值勢必會極大下降，所以本次分析中，加速度的極值具有很少的參考性。</p><p><b>　　五、動力學特性評價</b></p><p>　　5.1“飛脫”和

69、“反跳”</p><p>　　通過上述計算得到了氣門的實際升程曲線之后，就可以對該配氣機構進行動力學性能評價了，即通過對比由運動學分析得到的氣門理論升程曲線與由動力學分析得到的氣門實際升程曲線，檢驗該配氣機構在發(fā)動機工作轉速范圍內是否會出現(xiàn)“飛脫”和“反跳”。</p><p>　　本文分別對發(fā)動機在額定轉速、高速2種情況下的動力學進行了分析，本發(fā)動機額定轉速為7000轉/分，所取的高速轉速

70、為12000轉/分，由圖中轉速與升程的關系曲線可以得知：</p><p>　　上面結果圖中顯示氣門實際升程曲線和理論升程曲線的對比，可以得出以下結論：</p><p>　　1.該發(fā)動機在額定轉速（7000r/min）和常用高速工作轉速（12000r/min）下的氣門實際升程曲線和理論升程曲線基本重合，說明該發(fā)動機配氣機構的動力性能良好，即凸輪型線的設計符合要求。</p>&l

71、t;p>　　2.在該發(fā)動機的高速轉速（12000r/min）的情況下，在一些上升工作段上出現(xiàn)了氣門實際升程大于理論升程，即出現(xiàn)了傳動鏈的脫節(jié)之處——“飛脫”，不過升程之差卻很小；在氣門落座階段也并未出現(xiàn)氣門“反跳”現(xiàn)象。因此，此發(fā)動機在高速轉速的情況下工作狀態(tài)也比較穩(wěn)定。</p><p>　　3.氣門實際開啟階段會略晚于理論開啟時間，同時實際氣門的最大升程也略小于理論氣門升程，這是由于氣門的彈性所引起； 在

72、氣門升程下降段出現(xiàn)了幅度較大的波動，可能是由于配氣機構傳動鏈的剛度不夠。</p><p>　　通過兩種狀態(tài)轉速下的動力學特性分析，該發(fā)動機并不存在較大的”飛脫”，“反跳”現(xiàn)象，說明該發(fā)動機的動力學性能良好。</p><p>　　5.2各參數(shù)對配氣系統(tǒng)的影響</p><p>　　從氣門運動微分方程可以看出分析，影響氣門的運動規(guī)律因素有很多，比如發(fā)動機轉速、傳動鏈剛度、

73、系統(tǒng)質量、凸輪型線等 []。</p><p>　　1. 配氣機構與凸輪轉速</p><p>　　氣門加速度與其所受慣性力與凸輪軸的轉速的平方成正比，在不改變配氣機構其它因素條件下，發(fā)動機轉速增大，從動力學分析結果來看氣門運動的加速度急劇增大。</p><p>　　2. 配氣機構與傳動鏈剛度</p><p>　　傳動鏈末端是氣門，氣門和氣門彈簧

74、成對組合使用。氣門與傳動件保持接觸，準備隨凸輪運動。隨著彈簧剛度的增加，氣門發(fā)生“飛脫”可能性降低，趨勢減弱?！帮w脫”的幅值和轉角減少但彈簧剛度增加另一方面將使得配氣機構受力增大，受力增大必然導師變形增加，磨損加劇。一個傳動鏈整體變形取決個各個環(huán)節(jié)的變形量，傳動鏈的整體剛度影響著氣門運動的平順性。</p><p>　　3.質量對配氣機構升程、速度、加速度影響</p><p>　　,隨著系統(tǒng)

75、質量的增大，系統(tǒng)所受慣性力成正比增大。質量減少，受力減小，所帶來的變形量對氣門實際升程曲線的影響也越小，“飛脫”可能性減少，，速度，加速度也隨之減小，運轉更加平穩(wěn)，可靠性增加。</p><p>　　4. 配氣機構與凸輪型線</p><p>　　凸輪型線決定氣門的理論升程規(guī)律，由于氣門機構實際為一個彈性系統(tǒng)，配氣機構真實的升程、速度、加速度是對應一個具有固有頻率的疊加[]。</p>

76、;<p><b>　　六、結 論</b></p><p>　　本文以某125型摩托車發(fā)動機頂置凸輪軸式配氣機構為研究對象，完整的完成了對該機構從理論分析，到模型的建立，到編程計算，到結論分析的全過程，并著重分析了該汽油機配氣機構運動學和動力學特性。具體實施為根據(jù)該高速汽油機的變搖臂比頂置凸輪軸配氣機構的幾何參數(shù)和凸輪的凸輪升程表在Matlab軟件中算出該配氣機構的氣門理論運動

77、規(guī)律，再運用動力學分析的方法建立配氣機構的單質量模型及微分方程式、在Matlab軟件編程計算得到氣門的實際運動規(guī)律。本文通過對比由計算結果繪制的氣門理論升程曲線和實際升程曲線，即該配氣機構的動力學性能分析，對該發(fā)動機配氣機構節(jié)能型了評估。</p><p>　　本課題研究中，分析結果表明，隨著轉速的提高，氣門漸漸會出現(xiàn)“飛脫”和“反跳”，但在發(fā)動機額定轉速以及較高轉速的情況下，該發(fā)動機沒有出現(xiàn)上述問題。因此，該配氣

78、機構總體性能良好，設計符合要求。如果需要對其他頂置凸輪式配氣機構做動力學分析，由于分析方法大致相同，可直接使用本課題中的動力學通用程序，以及單自由度模型，得到分析結果。因此本課題的研究擁有比較好的通用性。</p><p>　　配氣機構的動力學分析的一個有效的途徑便是單自由度動力學模型，該模型已經(jīng)能夠成熟的被龍格庫塔法的迭代算法求出。通過對該機構建立單自由度系統(tǒng)的動力學模型，由已知氣門理論運動規(guī)律求氣門實際運動規(guī)律

79、的動力學計算，得到了氣門實際運動規(guī)律曲線，經(jīng)過對曲線進行分析后，對該發(fā)動機配氣機構進行評估，從而保證所選擇的設計參數(shù)和凸輪型線能使配氣機構平穩(wěn)正常工作。因此這個模型還兼具優(yōu)化設計的作用。但是頂置凸輪軸式配氣機構最大的特點便是搖臂比是變化的，并且變化的幅度比較大，在不同的轉角下?lián)u臂和系統(tǒng)的剛度都是隨著搖臂比的變化而變化的，將系統(tǒng)的剛度當作一個定量，影響了分析的精確性。實踐表明單自由度模型對氣門機構的運動情況做出的分析基本上能達到工程精度要

80、求，但更詳細的情況，諸如傳動鏈的飛脫究竟發(fā)生在哪一環(huán)節(jié)上，各接觸副和氣門的動載荷究竟有多大以及氣門彈簧的震顫問題等等，是無法在單自由度模型中得到，因此往往在配氣機構基本方案確定后，進一步用多自由度模型進行一些補充驗算。</p><p><b>　　七、展 望</b></p><p>　　配氣機構在發(fā)動機中控制著換氣過程，其運動的設計直接影響內燃機的可靠性及性能。發(fā)動

81、機配氣機構應保證充氣系數(shù)盡可能高，開啟和關閉進排氣門的時刻精準，各氣缸換氣良好，合理。四沖程發(fā)動機大多采用凸輪——氣門式配氣機構，工作可靠，運轉穩(wěn)定。本次課題針對頂置凸輪式配氣機構做了運動學分析計算，建立了單自由度模型進行動力學分析，計算工作量不大。單自由度模型是經(jīng)常使用在氣門機構設計方案初選階段，諸如傳動鏈的飛脫究竟發(fā)生在哪一環(huán)節(jié)上，彈簧的顫振是否會導致過大的應力或某些簧圈相碰，各接觸副和氣門桿的動載荷有多大等等，就無法由單自由度模型

82、分析中得到[]。因此往往在氣門機構基本方案確定后，對于一些特殊工況下尤其是在接近極限情況下，機構的響應究竟具有何種特性，還需要我們進一步分析驗證，這就需要我們用到多自由度模型進行一些補充驗算，雖然多自由度模型分析得到的結果更加準確，更具有參考性，但是由于自由度數(shù)多，從模型分析所得的運動微分方程就越為復雜，求解復雜的微分方程是一項十分繁瑣的工作，目前并沒有哪一款軟件或哪一種算法可以被通用在復雜的微分方程求解問題上，所以對配氣機構多自由度模

83、型的研究還有待進行。</p><p>　　對于配氣機構，相對于傳統(tǒng)的凸輪控制式配氣機構，現(xiàn)在有一種比較新型的設計——無凸輪發(fā)動機配氣機構（camless engine valve trains）。正如其名，該配氣機構中沒有凸輪，與目前比較偏重機械方面的設計來說，該機構利用了電磁控制，結構更小巧、簡單，并且控制精確，在控制氣門正時和氣門升程方面有著更好的效果。該發(fā)動機缸頭設計減少了20％的燃料消耗和污染物排放量。該

84、機構由智能氣門驅動系統(tǒng)（SVA）代替了由凸輪帶、凸輪軸和液壓凸輪配件組成的傳統(tǒng)的機械發(fā)動機氣門機構。在無凸輪發(fā)動機中，每個發(fā)動機氣門由一對置于缸蓋上表面上的氣門驅動器直接控制。</p><p>　　發(fā)動機基本排量不變的前提下進行配氣機構的優(yōu)化設計師當今發(fā)動機配氣機構設計的發(fā)展趨勢所在，提高各方面性能指標。為了更大范圍內提升內燃機動力指標、經(jīng)濟指標和生態(tài)指標，早先有多氣門配氣機構，隨后又發(fā)展有可變氣門運動配氣機構，

85、而傳統(tǒng)配氣機構相對于可變氣門運動配氣機構提升內燃機動力指標、經(jīng)濟性的能力這方面的能力就要遜色許多[]。對配氣機構的研究還有很大的發(fā)展空間，只有不懈的創(chuàng)新、追求和探索，才能使之不斷改進，不斷完善，更加廣泛的被使用。</p><p><b>　　致 謝</b></p><p>　　光陰似箭，歲月如梭，匆匆大學四年即將走到盡頭，是結束亦是開始。這里有我熟悉的教學樓、實驗室

86、，有我熱愛魅力校園風光。不經(jīng)意間，在這熟悉的校園中，我度過了人生中最為珍惜的一段時光。大學四年期間，有勤奮時的朝暮勤勉，有懈怠時的慵懶無為。但始終感受到整個時代的進步與科學技術發(fā)展的迅速，內心時而會不斷激勵自己，將自身埋在圖書館浩瀚的書海里，以提升個人素質，思想境界，學習上也不斷勤勉自己，努力專研專業(yè)學術知識，無奈大學四年匆匆過去，心中不免深深嘆息韶光易逝。 </p><p>　　畢業(yè)前最后一次鄭重的檢驗便是畢

87、業(yè)論文，這是對自我所學知識，科研能力的此全面檢驗。從課題的開展到論文完成的整個過程中，最重要的是給予我指導的導師阮登芳教授，另外本專業(yè)的許多同學也給予我不少軟件分析上的幫助，在這里我要由衷的感謝他們，謝謝你們!首先，我要感謝學校和學院為我們建設了這樣一個美麗的校園學習環(huán)境，和專業(yè)的教學實驗室實驗器材，讓我們可以理論聯(lián)系實際，學習更加深入。同時我還要感謝本學院所有的老師們，感謝他們在這四年來對我的悉心教導。師之所存，道之所存也。正是有了他

88、們這一群教學嚴謹細致、工作一絲不茍的教育工作者，科學知識，人生哲理才得以最大化的延續(xù)。他們的直言不諱，使我在今后的人生中收獲甚多。 </p><p>　　其次，對于我的指導老師阮登芳教授我要特別感謝一次，開題，譯文翻譯，國內外文獻查閱，軟件分析，解題整個論文的撰寫工作無不透露著老師的心血，整個過程都有老師的陪伴，這份良苦用心讓我感動不已。在老師的悉心指導下，艱苦奮斗，查閱各種文獻方法這篇論文才算順利完成。在這個

89、過程中，導師指點我查找分析文獻資料，與我共同分析課題分析結果的合理性，一同探究期間計算程序出現(xiàn)的問題。每一次的幫助我都十分受用，在此對老師致以深深的謝意。</p><p>　　最后，本專業(yè)同學以及其他專業(yè)同學在我整個畢業(yè)設計過程中都有或多或提供指導和幫助，這讓我深深感受到同學間的那種情誼，在這論文完成的時刻，我要謝謝你們，謝謝你們的精神鼓勵和技術支持，風里雨里我們都曾一起走過，讓我在整個過程中感到無比幸福！感謝所

90、有關心和幫助過我的老師、同學、朋友，謝謝你們！在論文完成之際，我要向我的指導老師致以最衷心的感謝和深深的敬意!</p><p>　　最后的最后，衷心地感謝在百忙之中抽空評閱本文和參加答辯的各位專家、教授!</p><p>　　附錄A：matlab運動學分析程序</p><p>　　已知凸輪升程--轉角（y，α）數(shù)據(jù)求氣門升程—轉角（h，αr）matlab程序實現(xiàn)&l

91、t;/p><p><b>　　la=39.61;</b></p><p><b>　　lc=30.28;</b></p><p><b>　　ro=12.5;</b></p><p><b>　　rs=30;</b></p><p>　

92、　ymax=6.532;</p><p><b>　　p=10;</b></p><p><b>　　o=0.05;</b></p><p><b>　　lb=33.20;</b></p><p><b>　　z=1.65;</b></p>

93、<p>　　yys=gradient(ys)./gradient(as);</p><p>　　for i=1:length(ys)</p><p><b>　　y=ys(i);</b></p><p><b>　　a=as(i);</b></p><p>　　yy=yys(i); &

94、lt;/p><p>　　R=(yy.^2)+(ro+rs+y).^2;</p><p>　　w=acosd((la.^2+lc.^2-R)/(2*la*lc));</p><p>　　rr=R.^0.5;</p><p>　　e=atand(yy/(ro+rs+y));</p><p>　　k=acosd((lc.^2+R

95、-la.^2)/(2*lc*rr)); </p><p>　　Rd=ro+rs+ymax;</p><p>　　kd=acosd((lc.^2+Rd.^2-la.^2)/(2*lc*Rd));</p><p>　　ar1=a+e+(k-kd);</p><p>　　ar2=a+e-(k-kd);</p><p>　　w

96、o=acosd((la.^2+lc.^2-(ro+rs).^2)/(2*la*lc));</p><p>　　bo=asind((p+z+o)/(lb));</p><p>　　b=bo-(w-wo);</p><p>　　h=lb*(sind(bo)-sind(b))-o;</p><p>　　shuchu(1,i)=w;</p>

97、;<p>　　shuchu(2,i)=ar1;</p><p>　　shuchu(3,i)=ar2;</p><p>　　shuchu(4,i)=h;</p><p><b>　　end</b></p><p>　　[FileName,PathName]=uiputfile({'*.xls'

98、,'EXCEL Files(*.xls)';'*.*','All Files(*.*)'}liming.xls');</p><p>　　xlswrite([PathName,FileName],shuchu); liming.xls</p><p>　　附錄B:動力學分析計算基本程序</p><p><

99、b>　　k=11700;</b></p><p><b>　　nc=12;</b></p><p>　　ps0=205.8;</p><p>　　ks=41.356;</p><p>　　beta1=17.73;</p><p><b>　　M=0.695;</b

100、></p><p><b>　　det=0.5;</b></p><p><b>　　y=0;</b></p><p><b>　　yy=0;</b></p><p>　　A=zeros(1,386);</p><p>　　B=zeros(1,38

101、6);</p><p>　　C=zeros(1,386);</p><p>　　for n=1:1:386</p><p>　　k1=(k/(36*nc^2*M))*Y(2*n-1)-ps0/(36*nc^2*M)-(k+ks)/(36*nc^2*M)*y-beta1/(6*nc*M)*yy;</p><p><b>　　C(n)=

102、k1;</b></p><p>　　k2=(k/(36*nc^2*M))*Y(2*n)-ps0/(36*nc^2*M)-(k+ks)/(36*nc^2*M)*(y+det/2*yy)-beta1/(6*nc*M)*(yy+det/2*k1);</p><p>　　k3=(k/(36*nc^2*M))*Y(2*n)-ps0/(36*nc^2*M)-(k+ks)/(36*nc^2*

103、M)*(y+det/2*yy+det^2/4*k1)-beta1/(6*nc*M)*(yy+det/2*k2);</p><p>　　k4=(k/(36*nc^2*M))*Y(2*n+1)-ps0/(36*nc^2*M)-(k+ks)/(36*nc^2*M)*(y+det*yy+det^2/2*k2)-beta1/(6*nc*M)*(yy+det*k3);</p><p>　　y=y+de

104、t*yy+det^2/6*(k1+k2+k3);</p><p><b>　　A(n)=y;</b></p><p>　　yy=yy+det/6*(k1+2*k2+2*k3+k4);</p><p><b>　　B(n)=yy;</b></p><p><b>　　end</b>