发动机配气机构运动学及动力学分析-重庆大学车辆工程专业本科学生毕业设计

发动机配气机构运动学及动力学分析-重庆大学车辆工程专业本科学生毕业设计

重庆大学本科学生毕业设计（论文）重庆大学本科学生毕业设计（论文）发动机配气机构运动学及动力学分析学生：黎明学号：指导教师：阮登芳（教授）专业：车辆工程重庆大学车辆工程学院二零一七年五月\n重庆大学本科学生毕业设计（论文）GraduationDesign(Thesis)ofChongqingUniversityKinematicsanddynamicsanalysisforenginevalvetrainUndergraduate:LiMingSupervisor:Prof.RuanDengfangMajor:VehicleEngineeringCollegeofVehicleEngineeringChongqingUniversityMay2017\n重庆大学本科学生毕业设计（论文）摘要摘要配气机构是发动机的重要组成部分，其设计的合理与否直接影响到发动机的充气效率以及换气质量，因此对发动机的动力性、燃油经济性、可靠性、有害物质排放、发动机噪声和振动有较大的影响[[]秦元伟.摩托车发动机配气机构的运动学和动力学分析[D].重庆大学,2008.]。而顶置凸轮轴式配气机构由于能适应更高的转速而在许多小型汽油机中广泛使用。但是顶置凸轮轴由于摇臂传动比是变值，所以其几何关系要复杂很多[[]廖祥兵,王军,张立军.顶置凸轮轴凸轮型线的优化设计[J].内燃机学报,2001,19(6):588-592.]。本文在已知凸轮对摇臂的运动规律的条件下，针对某125发动机的配气机构，经理论分析运动学规律，并用matlab计算出其气门对转角的理论升程、速度、加速度。在考虑气门间隙及传动机构变形的影响下，建立配气机构运动的单自由度模型，得出运动二阶微分方程。利用matlab采用龙格——库塔法计算出气门的实际运动规律，对比气门实际升程和理论升程，对该发动机配气机构的“飞脱”、“反跳”以及运转的平稳性进行动力学特性评价。从而完成了整个配气机构的运动学及动力学计算。关键词：运动学，动力学，配气机构，matlab，龙格库塔法IV\n重庆大学本科学生毕业设计（论文）ABSTRACTABSTRACTValvetrainisanimportantpartoftheengine,whichhasdirectlyaffecttotheengine'svolumetricefficiencyandthequalityofventilation,sothereisalsoagreaterinfluencetotheenginepower,fueleconomy,reliability,emissionsofharmfulsubstances,enginenoiseandvibration.Becausetheoverheadcamshaftvalvetraincanadapttothehigherspeed,itiswidelyusedinmanysmallgasolineengine.Butfortheoverheadcamshaft,thedriveratiooftherockerischangedbythetime,soithasamorecomplexgeometryrealationship.Withknowingthelawofmotionofcamontherocker'scondition,inthearticle,thedisplacementofthevalveiscalculated.Inconsideringthevalveclearanceandthedrivemechanismundertheinfluenceofdeformation,theactualvalvemovementruleiscalculatedbyusingtheRunge-Kuttamethod,andtherunningspeediscalculatedwiththeconditionsthatthetransmissionchainisflyingoffandreboundwhicharenotinthenormalconditions.Thenthekinematicsanddynamicscalculationsofthevalvetrainarecompleted.Andonthisbasis,withjoiningthemodalanalysisofthevalve,thetheoreticalbasisforthevalvetraindesignareprovided.Avalvetrainofa125motorcycleengineischosenfortheobjectofstudyinthissubject.Keywords：Valvetrain,Kinematics,Dynamics,MatlabIV\n重庆大学本科学生毕业设计（论文）目录目录摘要IABSTRACTII一、绪论11.1课题研究意义11.2课题国内外研究状况21.2.1国外研究现状21.2.2国内研究现状31.3课题研究背景31.4课题研究内容4二、气门机构的主要设计要求6三、运动学分析83.1凸轮廓线预处理83.2气门理论运动规律与凸轮轮廓的关系93.3运动学理论分析后的计算结果11四、动力学分析134.1动力学理论分析134.2摇臂比i154.3摇臂刚度计算164.4解动力学微分方程174.5动力学分析结果19五、动力学特性评价245.1“飞脱”和“反跳”245.2各参数对配气系统的影响24六、结论26七、展望27致谢28附录A：matlab运动学分析程序29IV\n重庆大学本科学生毕业设计（论文）目录附录B:动力学分析计算基本程序30参考文献31IV\n重庆大学本科学生毕业设计（论文）绪论一、绪论本课题以某125型摩托车发动机的顶置凸轮式配气机构为研究对象，分别对其进行了运动学分析、刚度计算、以及动力学分析，并由所得到的数据对该机构进行动力学评估，为该发动机配气机构的合理设计奠定基础。1.1课题研究意义配气机构是发动机结构中最复杂、工作最繁重的部件之一，固然也是发动机的重要组成部分，配气机构控制着进排气门的开启时刻和开启速度，根据每个气缸的工作顺序，控制进气和排气，实现换气功能。配气机构可以视为两部分组成，即气门组和气门传动组。根据凸轮轴的位置不同，将配气机构分为三类：下置凸轮轴式配气机构、中置凸轮轴式配气机构以及顶置凸轮轴式配气机构。本课题便是以凸轮轴顶置式配气机构为研究对象展开工作。配气机构的运转平稳性一部分决定了发动机振颤程度的大小，已成为一项非常重要的设计参考性能指标[[]徐小明.车用发动机配气机构动力学及其振动特性分析[D].华中科技大学,2007.]。对配气机构的研究内容归纳起来主要有两个方面，一方面是零部件的设计，另一方面是机构的动力学问题。零部件设计又可以简单分成三类：气门及气门等零部件设，计凸轮型线优化设计，摇臂机构运动学设计。因为凸轮作为整个机构的原动件，它直接控制整个机构的运动，所以凸轮型线的设计至关重要。而对于机构动力性能的研究，又主要表现在凸轮型线如何如何影响气门的运动规律。配气机构的设计最初为刚性设计，但没有绝对刚性的零部件，在机械应力以及热应力作用下总会存在变形，所以逐步发展为弹性设计。高速汽油机由于转速较快，每分钟转速可达5000转以上，为保证充气效率，都采用顶置式气门装置。这种装置都适合用凸轮轴的三种安装形式，但是，如果采用下置式或者中置式的凸轮轴，由于气门与凸轮轴的距离较远，需要气门挺杆和挺柱等辅助零件，这样的结构势必造成气门传动件较多，结构复杂，发动机体积也不可避免的变得很大，而且在高速运转下还容易产生噪声，运转平稳性较差。所以，现代轿车发动机一般都采用顶置式凸轮轴，以期待能改变这种现象。将凸轮轴配置顶置靠近气门，不仅缩短了凸轮轴与气门之间的距离，而且省略了挺杆和挺柱，从而简化了传动机构，使得发动机的结构变得紧凑。更重要的是，这种凸轮轴顶置的安装方式可以减少整个系统往复运动的质量。在国内的，不少小型汽油机已采用了凸轮轴顶置的配气机构。顶置凸轮轴配气机构的系统刚度虽然大，5\n重庆大学本科学生毕业设计（论文）绪论但实际上仍是弹性系统，可看作一个“弹簧——质量系统”[[]秦培军.一种新型可变气门正时机构的研究[D].同济大学,2006.]。当气门机构运转工作的时候，机构会产生弹性变形，再加上热应力的影响会使气门运动产生很大程度的变化，与理论分析脱离，使传动链脱节，将出现一系列不利于发动机运转的现象，诸如气门的开闭不正常，振动和噪声增大、零件易损坏等现象，从而限制了转速的提高[[]李香梅,张文平,黄刚.柴油机配气机构动力学计算模拟[C]//全国振动工程及应用学术会议.2006.]。随着发动机的转速增加，机构运动件的振动也将增大。所以在发动机高速运转时，情况将更严重，甚至出现传动链脱节——“飞脱现象”，及过高的落座速度——“反跳现象”，甚至破坏配气机构的正常工作。因此，配气机构的运动学以及动力学分析工作显得尤为重要，并且值得设计工作者对其进行不断的优化设计，以提高发动机综合性能。1.2课题国内外研究状况1.2.1国外研究现状处在21世纪的今天，计算机能力已然十分强大，尤其是在计算分析领域。越来越多设计过程中运用CAE仿真技术，仿真结果相比理论分析计算结果也能更加有效的预测配气机构的性能。Seon-YangHwang等人发现气门的冲击噪声主要是脉冲噪声，他们分别对两种不同挺柱的配气机构进行噪声分析，得出结论：既要保证发动机性能又同时达到降低噪声的目的，可以通过优化缓冲阶段的凸轮型线来实现。NitinP.Gokhale运用三维CFD分析气门座温度与润滑液流速的关系，气门座温度升高，润滑液流速降低，两者成负相关。与此同时还对配气机构的磨损做了深入研究，发现磨损主要集中在摇臂头气门小端的磨损，凸轮-挺柱摩擦磨损和气门-气门座冲击磨损，这三者占整个磨损的60%~70%。A.K.Jamkhande在凸轮型线设计研究中采用了B-样条曲线，设计的凸轮型线不仅更加平滑而且有更低的越度，加速度曲线上不光顺的地方明显减少，在凸轮尖端具有更小负加速度。RSPrabakar和SPMangalaramanan将气门弹簧视为弹性体分析，运用ADAMS发现气门反跳一般不出现在低速发动机上，并且对于高速汽油机，也只集中出现在极速转速附近；并且发现双气门式弹簧的刚度是一个变化的值，随着转速的提高，刚度是会有所减小的，刚度的波动与变化的摇臂比和线圈的不确定运动有关。YeongchingLin等在分析单个气门模型时，用动力学分析软件得到配气机构零件的载荷激励，在有限元模型中去评估每个零件的许用应力和最大应力，大大的节约了配气机构的开发时间与成本。这一流程不仅为丰田汽车技术中心设计配气结构提供了指南，也为同类型分析建立了一个基本方法模型，但实则也是计算机软件技术进步带来的设计上的便捷。5\n重庆大学本科学生毕业设计（论文）绪论1.2.2国内研究现状由于我国近代科学技术起步较晚，导师我国对发动机配气机构的研究也相对晚了许多，樊久铭等人针对装有液压调节器的CA488汽油配气机构，提出了“4+N”多质量动力模型。刘忠民等研发了可用于配气机构的动力学特性和耐久性研究，可以控制配气机构温度和润滑条件的综合试验系统。李兴兰利用ADAMS软件构建气门机构的动态模型，获得振动和噪声分析与预测的主要部件之间的相互作用，提供准确的数据和边界条件。曹志芬将气门机构的刚性以简单的机械方式进行测量，使配气机构弹性系统模型的动力学分析有了最基本的刚度参数。基于CAD/CAE综合开发技术，朱文波等探讨发动机气门机构数值研究的思路和过程，为气门机构设计，动力学分析和优化提供参考，经验设计以定量形式出现，产品开发效率高。王晓辉等人采用反向优化设计方法重新设计凸轮轮廓，首先用测量设备测量实际轮廓，然后使用三维软件来适应凸轮轮廓，进行拟合，然后进入AVL软件优化，既确保了配气机构的基本性能，又提高了凸轮型线的连续性和加工性。白贺研究了凸轮轴扭转振动对气门机构的影响。通过扭转振动实验提取凸轮轴的扭转角，并测量凸轮轴的动态转矩。结果表明，凸轮轴的扭转振动随着转速的增加和凸轮轴的扭转振动的影响而增大，凸轮轴的扭转振动也增加了凸轮轴的扭转振动。1.3课题研究背景发动机的充气效率以及换气质量很大程度上受配气机构的设计影响，而发动机的动力性、燃油经济性、有害物质排放等各方面性能业余配气机构的设计有关。随着国际环境对发动机性能和排放要求不断提高，配气机构的研究也越来越深入，尤其是在弹性系统下的动力学特性研究更加迫切。随着制造商对发动机配气机构的进一步设计改进，目前4气门发动机在市场上已经越来越普遍。大功率低燃消耗从来都是车主最为在意的两个点，并且也是设计思想中起主导作用的两大因素。近来可变气门正时及升程控制技术在解决低油耗问题上，变现十分突出，解决了由单个气门控制气体进出燃油供给时，充气效率低，燃油时而供给不足，时而供给过量，从而造成发动机功率不够，排放严重等问题，使发动机在更加省油的同时还能在全段转速输出期间都更有力，动力更强劲，因此颇为受广大车主的青睐。凸轮轴顶置使配气机构凸轮轴在旋转中的负荷相应减小，并且对于凸轮轴和气门弹簧的要求也降到了最低。这是因为它在设计上没有挺柱、摇臂和推杆，直接通过凸轮轴上的凸轮来驱动气门开闭，这使5\n重庆大学本科学生毕业设计（论文）绪论得该机构在结构上有了大大简化。同时从维修角度来看，这也降低了成本。因此顶置凸轮轴式配气机构从此进入设计者眼中最关注的对象，由于顶置凸轮轴式配气机构优越的性能，简介的传动链结构，使得发动机体型有所减小，市场上各类型的发动机对这种配气机构的采用率也越来越高。另外，从物理性能的角度来看，凸轮轴架空设计思路进入和排气通道转小，气流阻力小，使气体出来更加平滑。因此随着内燃机结构上的不断强化以及性能上面的优化，能够适应更高转速的顶置凸轮轴配气机构，更能适用于大多数新型汽油机。对于现代的中小型车用发动机，都有着较高的额定转速，因此顶置凸轮轴式配气机构也越来越应用广泛，它不仅减少了从凸轮到气门之间的传动零件，更重要的是减轻了配气机构的运动质量，提高了系统的刚度，使配气系统可在很高转速下正常工作[[]余志敏.柴油机配气凸轮型线优化设计及其配气相位优化[D].武汉理工大学,2009.]。每个气门在发动机的一个工作循环中只启闭一次，设由此而引起的振动，在气门关闭期间能完全衰减。这样配气机构每个循环的条件都相同，因此，只需研究气门的一次启、闭过程即可[[]施玉春.DA1发动机配气机构可靠性研究[D].哈尔滨工程大学,2006.]。有必要在气门机构的弹性变形下计算出气门运动的模型，但有必要计算出单一质量模型最为简洁的模型。但顶置凸轮轴的摇臂比发生变化，几何关系比多自由度模型更复杂，需要经过一系列转换才能从平底挺柱运动中知道凸轮从动位移，速度，加速度规律，计算阀门的位移，速度和加速度的规律。为了得出更为理想的气门运动规律，而且从所需的阀门运动规律来看，反向找到凸轮轮廓，但这个计算比较复杂，并且所得的凸轮轮廓值对应的凸轮角度不均匀，还必须使用插值方法来计算凸轮的均匀角度对应凸轮廓线值。所以一般使用选择的凸轮轮廓，找到阀门运动规律，然后计算动态性能。实际气门落座速度在允许范围内。基本可以认为凸轮轮廓设计满足要求。1.4课题研究内容1．配气机构运动学分析。在传动链结构上进行发动机运动学理论分析，通过matlab编程分析计算，得到分析运动学的一种通用算法。研究配气机构理论运动规律，编写动力matlab学分析程序，获取气门理论升程曲线，速度曲线，加速度曲线，提供动力学分析原始数据。2．配气机构刚度计算。利用Catia为配气机构建模，并导入Workbench中进行刚度计算，提供动力学分析分析所需数据。3．配气机构动力学计算。采用龙格库塔法，利用Matlab编写变摇臂比顶置凸轮轴配气机构的动力学计算通用程序，并得出气门实际升程曲线，速度曲线以及加速度曲线。5\n重庆大学本科学生毕业设计（论文）绪论4．对所选发动机配气机构的动力学特性进行评价，5\n重庆大学本科学生毕业设计（论文）气门机构的主要设计二、气门机构的主要设计要求发动机中的运行必须有着新气不断进入，废弃不断排除的过程。较早采用的控制其进气、排气的方式有很多，如滑阀控制（回流扫气二冲程机）、旋转阀控制（有过个别试验机型）、气门控制、混合控制（有排气门的直流扫气二冲程机）等等。对于当今市场主流四冲程发动机而言，气门控制进排气的方式已经逐渐稳定下来，因其可靠的密封性能，便拆装的优点，气门已经在配气机构中占据着不可替代的地位。现有产品几乎全都采用气门控制进排气。气门机构控制着换气过程，其设计好坏直接影响着四冲程发动机方方面面的性能。首先，由内燃机原理阐述，发动机的最大功率转速nmax正比于进气门通路面积Av除以活塞排量Vh，即nmax∝（Av/Vh），而升功率Nl∝peNnmax。拿同一类型的发动机作比较，因最大功率工况的peN差不多一样大，可以说他们的升功率Nmax∝nmax∝Av/Vh。其次，进气门定时影响充气系数随转速变化的情况，影响泵气损失，影响换气质量，因此对发动机的动力性，燃料经济性和有害排放有影响。此外，配气机构的因运动副的摩擦而发生摩擦损失，与发动机的转速有关，总的说来发动机转速越低摩擦损失的比重越大，转速越高，配气机构摩擦损失所占比重越低。与此同时发动机的机械效率和油耗率也随着配气机构的机构摩擦损失与发动机中低转速负荷工况下变化而变化。配气机构各个工况转速下所产生的噪声以及工作不平稳的振动很大程度影响整个发动机可靠性与否和噪声水平高低。在高热负荷下，由于材料性能的限制气门将发生热裂、凸轮在高速运转下受到较高接触应力的影响发生接触疲劳、在各种应力综合作用下气门弹簧也可能发生疲劳失效。诸如此类都是重要的可靠性问题，在设计过程中不得不引起重视，对动力学特性做出做够有效的分析评价。气门机构是弹性系统，在工作过程中一旦发生振动，就有可能产生“飞脱”，“反跳”，气门“提前落座”一系列异常现象，并且产生这种异常现象的概率会随着转速的提高而增加，这些等异常现象的出现，将使得凸轮-从动件接触面的磨损加剧，从而影响着发动机的寿命和运转平稳性以及工作可靠性。根据这些情况，对气门机构的设计理念和主要要求可总结为：a)为使充气量足够大，在气门口要保证足够大的气流通路面积；b)机构的动力学特应当满足机构的“飞脱”和“反跳”不出现在任意转速工31\n重庆大学本科学生毕业设计（论文）气门机构的主要设计况下；a)为了抑制机构接触表面的磨损现象，不应该有过高的热应力和机械应力发生在重要零件表面，同时在材料方面也必须考虑到使用性能上的耐久性作；b)便于制造和维修，成本较低。31\n重庆大学本科学生毕业设计（论文）运动学分析三、运动学分析配气机构运动学所研究的内容包括：①凸轮从动件（下置凸轮轴式气门机构中的挺柱，顶置凸轮轴式的摇臂）运动规律与凸轮轮廓开关系，以及从动件与凸轮的相对滑动关系（本文不研究相对滑动）；②在不考虑机构弹性的条件下，凸轮从动件的运动规律与气门运动规律的关系。3.1凸轮廓线预处理对于普通的下置凸轮式配气机构的运动学计算，只要乘以或除以固定的摇臂比即可，但是针对本课题研究的配气机构采用了凸轮——摇臂——气门如图1所示这种顶置凸轮式配气机构，这种配气机构的摇臂比是随着凸轮转角变化而变化的，要通过配气机构运动件的几何关系及其运动规律利用三角函数关系换算得到。对于经摇臂驱动气门的机构来说，无论摇臂是处于摇臂中间还是处于摇臂一端，运动学分析方法是一样的。顶置凸轮轴的凸轮轮廓通常是按照它驱动平底挺柱时的挺柱升程规律y=f(α)给定、制造并用平底测头检验的。运动学分析要解决由已知的y=f(α)求不同凸轮转角下摇臂摆角，再进而求气门理论运动规律的问题。初始获得的凸轮廓线数据是不光顺的，显然不符合设计理念要求。将数据导入matlab分析计算之前需要将给定的凸轮廓线数据在matlab中进行光顺预处理，利用matlab分析软件强大的数据插值，曲线拟合的功能，将凸轮廓线首先处理成成一条光滑的曲线，但是这时并没有确切的升程与凸轮转角的相应的函数关系式，只用采用样条插值的方法得到一组数据。导入excel制作成凸轮升程一览表，如附录表C，为matlab编程计算做准备。31\n重庆大学本科学生毕业设计（论文）运动学分析3.2气门理论运动规律与凸轮轮廓的关系图2给出了顶置凸轮轴及摇臂从动件的关系图。图中Ad是摇臂柱面中心A在凸轮以顶点D与摇臂接触时的位置。任一瞬时凸轮的转角以OcAd与OcD的夹角αr表示，摇臂的位置则以OA与OOc的夹角W表示。图中还给出了凸轮上同一点F与平底挺柱接触时的升程y和转角α，以及半径为rs的球面挺柱接触时的升程ys和凸轮转角αs。本次研究的嘉陵125发动机配气机构，是一个朝枢轴转的凸轮，规定以气门全开时的转角为00，在凸轮上升段转角为负，下降段转角为正。在得到y=f(α)凸轮升程规律的基础上，计算气门运动规律是分析的基础。首先从每一个（y，α）数据求对应的（W，αr）数据。在△OAOc中，利用三角函数关系可得：图2顶置凸轮轴与摇臂的运动关系图W=arccosla2+lc2-R22lalc（3.1）式中la、lc和R的含义见图2，针对本课题所分析机构，通过摇臂尺寸相关数据计算得知la=39.61mm,lc=30.28mm，而且由△OcAE知:R=OcE2+AE2=(y')2+(r0+rs+y)2（3.2）其中凸轮半径r0=12.5mm，rs=30mm。y'代表着升程对转角的一阶导数，31\n重庆大学本科学生毕业设计（论文）运动学分析y'=dydα×180°π，单位（mm/弧度）。凸轮转角αr与αs，α的关系为：αr=αs-K-Kd=α+ε-(K-Kd)（3.3）ε、K和Kd的计算式分如下：ε=arctanOcEAE=arctan⁡y'r0+rs+y（3.4）K=arccos⁡lc2+R2-la22lcR（3.5）Kd=arccos⁡lc2+Rd2-la22lcRd（3.6）式中：Rd=R│α=0=r0+rs+ymax（3.7）接着再由（W，αr）求对应的气门升程h。由图3可见，当摇臂与凸轮基圆接触时，摇臂长臂与气门杆端面间有间隙Δ=0.05mm,而气门杆端面比摇臂中心高出z=1.65mm，这时摇臂的位置角为：W0=arccosla2+lc2-（r0+rs）22lalc（3.8）31\n重庆大学本科学生毕业设计（论文）运动学分析β0=arcsin⁡ρ+z+Δlb（3.9）其中调节螺钉球半径ρ=10mm，经过计算得，lb=33.28mm；与凸轮转角αr对应的β和气门理论升程h为：β=β0-ΔW=β0-(W-W0)（3.10）h=lbsinβ0-sinβ-Δ（3.11）至此，由每一组（y，α）计算对应的（W，αr）再计算对应的（h，αr）的过程就完成了。但是，等步长递增的α所对应的αr是不按等步长递增的。在进行动力学分析时，必须要有等步长递增的αr对应气门理论升程h数据。因此，对于理论分析运算结果，还必需进行拟合与插值处理才能得到所要求的升程—转角数据[[]谢宗法,LULi,孔超,等.内燃机顶置配气凸轮型线的设计[J].山东大学学报工学版,2008,38(4):84-88.]。3.3运动学理论分析后的计算结果将公式（3.1）到（3.11）用matlab编程计算得到如图4理论气门升程曲线，matlab程序参见附录1。图4.气门理论升程曲线图4气门理论升程曲线并利用matlab样条插值功能得到等步长递增的αr所对应的气门升程。拟合，插值处理后得到气门的速度，加速度如图5，其中速度与加速度值以转速7000r/min换算得出。31\n重庆大学本科学生毕业设计（论文）运动学分析图5气门速度以及加速度曲线至此，配气机构的运动学分析就已全部结束，气门的最大升程当αr=-6°，hmax=7.815mm；气门最大正速度αr=-40°，y’max=7.855m/s；气门最大负速度αr=12°，y’max=-14.47m/s。气门加速度情况略微复杂，上升段最大正加速度y”max=1713m/s2，下降段最大正加速度y”max=11500m/s2，最大负加速度y”max=-11600m/s2。加速度的剧烈波动是由速度曲线上的一些尖点造成，消除速度波动尖点，光顺速度曲线逆向优化设计凸轮廓线可以解决这一问题。若用高阶多项式拟合所得到的气门理论升程曲线，并再求导，得出气门理论速度以及加速度，将有如下结果：图6光顺处理后的速度及加速度曲线通过拟合前后对比，可以清晰看出，若凸轮型线数据给的足够光滑，我们所得结果也将会是光滑的，并且从加速度曲线上面可以看出，加速度最值大大降低，并且加速度没有突变，气门的运动不发生刚性冲击。所以在配气机构设计中，需要我们逆向设计优化凸轮型线。31\n重庆大学本科学生毕业设计（论文）动力学分析四、动力学分析4.1动力学理论分析实际气门机构是一个弹性系统，配气机构动力学正是要研究考虑弹性的情况下的气门的运动规律。运动学分析中得出的气门理论运动规律，这为进一步做动力学分析提供了可能，理论气门升程像是一个激振力函数，而真实情况下气门的运动归路则是对应这种激励下的响应。由于配气机构是有弹性的，相当于弹簧——质量系统，气门的实际运动规律会偏离运动学分析所得的理论运动规律[[]高巍,赵慧超.CA6110发动机配气机构采用滚子挺柱分析[C]//中国汽车工程学会发动机分会、中国内燃机学会汽油机煤气分会2002年度联合学术年会.2002.]。由于气门存在弹簧预紧力PS0，在刚消除气门间隙的瞬时气门还不会打开，要等气门机构发生一定变形而变形力≥PS0时才能顶开气门。因此气门的实际运动要比气门的理论运动始点晚一些。图7单自由度模型发动机在高速运行时，因弹性振动的影响，气门的实际运动情况与理论规律差别较大，甚至有时气门实际升程大于理论升程，这时由凸轮至气门的整个传动链将出现脱节之处，称作“飞脱”；另一方面，气门落座时刻早于理论落座点，以高速冲击气门座并可能“反跳”。正确设计的配气机构不容许在发动机工作转速范围内出现飞脱和反跳。为此配气机构的动力学特性需要进行专业的校核。为分析实际气门运动规律及相应的各构件受力情况，需要建立动力学分析模型，写出描述其运动情况的微分方程式然后求解。求解结果的可信度与模型本身有关，即取决于该模型能在多大程度上反映机构的特点。已提出的气门机构动力学计算模型有许多种，如单自由度模型、多自由度模型、有限元模型等等，其中单自由度模型是最基本而又最简单的一种模型。图9中给出了单自由度模型。其中M1是气门组建和弹簧的当量质量，M2是31\n重庆大学本科学生毕业设计（论文）动力学分析摇臂、推杆、挺柱等驱动件的当量质量；ks是气门弹簧刚度，k是从摇臂长臂端到凸轮轴为止的整个气门机构驱动部分的刚度，kz是气门座的刚度；b是驱动部分的当量内阻尼系数，bz是气门座阻尼系数，β1是反映气门组件粘性外阻尼及气门弹簧内阻尼的当量阻尼系数，β2是驱动部分的当量外阻尼系数。在模型中当量凸轮绕旋转，驱使从动件运动，我们认为转轴轴心是刚性的，从动件质量为0。从动件的位移Y在动力学分析之前已经从运动学分析中得到力理论运动规律，是由实际凸轮型线根据气门机构传动链的几何关系换算到气门理论升程[[]徐波.摩托车发动机顶置凸轮机构的分析与研究[D].华中科技大学,2002.]。计算y的起点（气门的初始位置）定在气门座被压缩了δ0=Ps0/ks之处。质量M1的位移y相当于气门的实际升程。模型中的运动质量通常取为：M1=Mvz+13ms（4.1）M2=IRlv2（4.2）IR为摇臂的转动惯量，用solidworks分析得到IR=81.2（kg·m^2），对于变摇臂比的气门机构就取其变化范围内的中值。经计算i取i=2.1。Mvz是气门、气门弹簧座和气门锁夹这些气门组件的当量质量，ms代表气门弹簧的质量，其中Mvz=0.033kg，ms=0.6kg，lv取为：12lbcosβv0+cosβd其中βv0=arcsinρ+zlb=20.2°，βd是气门全开时的摇臂位置角。达伦伯原理指出，作用于运动质量上的力与惯性力符合力平衡关系，故对于消除气门座变形阶段的运动情况应有：kY-y+kzδ0-y+bY-y-bzy-My-β1+β2y-ksy-Ps0=0整理后得到:My+b+bz+β1+β2y+k+kz+ksy=bY+kY+kzδ0-Ps0（4.3）式中M表示此阶段参与运动的质量，考虑到Y，y通常都用位移对转角的倒数Y',y'来表示，即Y=6ncY',y=6ncy',y=36nc2y'',所以上式可写成：y''+b+bz+β1+β26ncMy'+k+ks+kz36nc2My=Fyα（4.4）31\n重庆大学本科学生毕业设计（论文）动力学分析式中：Fyα=b6ncMY'+k36nc2MY+kzδ0-Ps036nc2M（4.5）相当于激振力，是“激振力函数”。它只决定于Y=fα，即只决定于凸轮轮廓和气门机构传动链的设计。其中Ps0是气门弹簧预紧力，其值为205.8N，nc是凸轮轴转速，本次研究针对2个转速的情况作了运动学分析，分别是本发动机的额定转速7000r/min，高速转速为12000r/min.上述的微分方程可以通过matlab，龙格库塔法迭代求解。4.2摇臂比i此摇臂比对前述运动学分析毫无意义，但在分析摇臂受力时需了解气门对摇臂的作用力PvR与凸轮对摇臂的作用力PcR之比。PvR作用于摇臂和气门杆端面的接触点C并指向B，PcR作用于摇臂和凸轮的接触点并指向A（见图2，图3）。若忽略摇臂的惯性力矩不计，可得摇臂比公式如下：i=lblc×cosβsinJ（4.6）其中J就是图2中OA与AF的夹角:J=arccosla2+R2-lc22laR+ε（4.7）若R，ε和β是对应等差递增α求出来的，则所得的J和i就对应不等差递增的αr。在动力学分析时，对于变摇臂比的气门机构就取其变化范围内的中值，图8摇臂比计算结果经matlab计算，取i=1.6。31\n重庆大学本科学生毕业设计（论文）动力学分析4.3摇臂刚度计算图9机构刚度测量简图一个机构的刚度是指弹性体抵抗变形（弯曲、拉伸、压缩等）的能力。顶置凸轮轴式结构，结构紧凑质量轻，体积小，但整个系统却具有较高的刚度，因此顶置凸轮轴式配气机构在很高转速下工况下同样能正常工作。因此对于顶置凸轮轴式配气机构，其很重要的一个优点就是机构刚度大，因此机构刚度的计算对本课题的研究又很重要的意义。图5所示为机构刚度测量的理论简图。由于气门与气门座刚度对机构刚度的影响很小，对动力学分析的影响也很小，因此本次分析假设气门与气门座为刚体，只考虑了弹簧和摇臂的刚度。而弹簧的刚度已知，因此对机构摇臂刚度的计算就成为了该配气机构刚度计算的主要任务。利用catia完成摇臂3D建模:图10Catia3D摇臂模型将其导入workbench分析软件中，摇臂材料选用合金铸铁，如图5施加约束，划分75209个单元网格：31\n重庆大学本科学生毕业设计（论文）动力学分析图11.划分单元网格并施加P=1000N作用力，在这个作用力下测得摇臂总变形最大为0.085mm。图12.摇臂应变所以摇臂的刚度：k=1000/0.085=1.17×104N/mm4.4解动力学微分方程下面的式子即为需要求解的方程：y''=Fyα-k+ks+kz36nc2My-b+bz+β1+β26ncMy'（4.8）31\n重庆大学本科学生毕业设计（论文）动力学分析即二阶常微分方程y''=fα，y，y'。因为在求解区间的起点，即凸轮转角为α0而气门间隙刚刚消除(Y=0)的时刻，y和y'是已知的，即yα0=0y'α0=0采用龙格库塔法从这些初值开始逐步迭代，可得出该微分方程式在一系列计算点αn(n=0，1，2，┅N)上的数值解。设原方程为y''=fα，y，y',并已知某一转角为αn时的升程和速度各为yαn=yn，y'αn=yn'。取转角步长为Δα=0.5°（一般取为0.5o~1.0o），则下一个计算点αn+1=αn+Δα上的升程yn+1和速度yn+1'可由yn和yn'用下式算出，即yn+1=yn+Δα·yn'+(Δα)26(k1+k2+k3)（4.9）yn+1'=yn'+Δα6(k1+2k2+2k3+k4)（4.9）式中的k1=fαn，yn，yn'——即y''（4.10）k2=fαn+Δα2，yn+Δα2·yn'，yn'+Δα2·k1（4.11）k3=fαn+Δα2，yn+Δα2·yn'+(Δα)24·k1，yn'+Δα2·k2（4.12）k4=fαn+Δα，yn+Δα·yn'+(Δα)22·k2，yn'+Δα·k3（4.13）这样从n=0初始点位置开始逐点迭代计算。每一αn下的升程yn和速度yn'得出后，该点的加速度y''=fα，y，y'也就得出了。微分方程中阻尼系数b、bz、β1和β2，在计算前必须获得，但从以前的相关研究来看这些系数目前还无法分别测定。但是阻尼系数b+β1+β2是可以推算的。已经有实践表明，阻尼值对整个动力学计算结果不敏感。因此在此次计算中取b=0，bz=0，β2=0，β1=0.2·k(M1+M2)。气门座刚度很大，为了简化计算，计算时不考虑气门座的变形。将公式（4.4）、（4.5）（4.9）~(4.13)用matlab编程。龙各库塔法要求迭代步长为0.5°，按照k1、k2、k3、k4的计算要求，需要我们找出转角增长0.25°倍数的转角所对应的左右理论升程Y的值，在matlab中样条插值每隔0.25°计算一次，便可得到所需Y值。31\n重庆大学本科学生毕业设计（论文）动力学分析图13.计算按步长增长的Y值详细matlab程序过程参见附录B。4.5动力学分析结果龙哥库塔法matlab分析程序见附录B，7000r/min额定转速下动力学分析计算结果：图14.7000r/min转速下实际升程与理论升程对比31\n重庆大学本科学生毕业设计（论文）动力学分析图15.7000r/min转速下速度曲线图16.7000r/min转速下加速度曲线最大升程为ymax=7.785mm，最大正速度Vmax=7.893m/s，最大负速度31\n重庆大学本科学生毕业设计（论文）动力学分析Vmax=-18.51m/s。加速度局部有较大波动，由速度在该区域波动引起，若解决该波动，加速度峰值势必会极大下降，所以本次分析中，加速度的极值具有很少的参考性。12000r/min高速转速下实际生程与理论升程对比：图17.12000r/min转速下理论升程与实际升程对比31\n重庆大学本科学生毕业设计（论文）动力学分析图18.12000r/min转速下速度曲线31\n重庆大学本科学生毕业设计（论文）动力学分析图19.12000r/min转速下加速度曲线最大升程为ymax=7.774mm，最大正速度Vmax=7.886m/s，最大负速度Vmax=-20.85m/s。加速度局部有较大波动，由速度在该区域波动引起，若解决该波动，加速度峰值势必会极大下降，所以本次分析中，加速度的极值具有很少的参考性。31\n重庆大学本科学生毕业设计（论文）运动学特性分析五、动力学特性评价5.1“飞脱”和“反跳”通过上述计算得到了气门的实际升程曲线之后，就可以对该配气机构进行动力学性能评价了，即通过对比由运动学分析得到的气门理论升程曲线与由动力学分析得到的气门实际升程曲线，检验该配气机构在发动机工作转速范围内是否会出现“飞脱”和“反跳”。本文分别对发动机在额定转速、高速2种情况下的动力学进行了分析，本发动机额定转速为7000转/分，所取的高速转速为12000转/分，由图中转速与升程的关系曲线可以得知：上面结果图中显示气门实际升程曲线和理论升程曲线的对比，可以得出以下结论：1.该发动机在额定转速（7000r/min）和常用高速工作转速（12000r/min）下的气门实际升程曲线和理论升程曲线基本重合，说明该发动机配气机构的动力性能良好，即凸轮型线的设计符合要求。2.在该发动机的高速转速（12000r/min）的情况下，在一些上升工作段上出现了气门实际升程大于理论升程，即出现了传动链的脱节之处——“飞脱”，不过升程之差却很小；在气门落座阶段也并未出现气门“反跳”现象。因此，此发动机在高速转速的情况下工作状态也比较稳定。3.气门实际开启阶段会略晚于理论开启时间，同时实际气门的最大升程也略小于理论气门升程，这是由于气门的弹性所引起；在气门升程下降段出现了幅度较大的波动，可能是由于配气机构传动链的刚度不够。通过两种状态转速下的动力学特性分析，该发动机并不存在较大的”飞脱”，“反跳”现象，说明该发动机的动力学性能良好。5.2各参数对配气系统的影响从气门运动微分方程可以看出分析，影响气门的运动规律因素有很多，比如发动机转速、传动链刚度、系统质量、凸轮型线等[[]张彤,焦永和,郑健健,等.175FM发动机配气机构的可靠性仿真分析[J].燕山大学学报,2005,29(3):237-241.]。1.配气机构与凸轮转速气门加速度与其所受惯性力与凸轮轴的转速的平方成正比，在不改变配气机31\n重庆大学本科学生毕业设计（论文）运动学特性分析构其它因素条件下，发动机转速增大，从动力学分析结果来看气门运动的加速度急剧增大。2.配气机构与传动链刚度传动链末端是气门，气门和气门弹簧成对组合使用。气门与传动件保持接触，准备随凸轮运动。随着弹簧刚度的增加，气门发生“飞脱”可能性降低，趋势减弱。“飞脱”的幅值和转角减少但弹簧刚度增加另一方面将使得配气机构受力增大，受力增大必然导师变形增加，磨损加剧。一个传动链整体变形取决个各个环节的变形量，传动链的整体刚度影响着气门运动的平顺性。3.质量对配气机构升程、速度、加速度影响F=d2Md2t,随着系统质量的增大，系统所受惯性力成正比增大。质量减少，受力减小，所带来的变形量对气门实际升程曲线的影响也越小，“飞脱”可能性减少，，速度，加速度也随之减小，运转更加平稳，可靠性增加。4.配气机构与凸轮型线凸轮型线决定气门的理论升程规律，由于气门机构实际为一个弹性系统，配气机构真实的升程、速度、加速度是对应一个具有固有频率的叠加[[]江孝礼.四行程发动机配气机构可靠性研究[D].武汉汽车工业大学武汉理工大学,2001.]。31\n重庆大学本科学生毕业设计（论文）结论六、结论本文以某125型摩托车发动机顶置凸轮轴式配气机构为研究对象，完整的完成了对该机构从理论分析，到模型的建立，到编程计算，到结论分析的全过程，并着重分析了该汽油机配气机构运动学和动力学特性。具体实施为根据该高速汽油机的变摇臂比顶置凸轮轴配气机构的几何参数和凸轮的凸轮升程表在Matlab软件中算出该配气机构的气门理论运动规律，再运用动力学分析的方法建立配气机构的单质量模型及微分方程式、在Matlab软件编程计算得到气门的实际运动规律。本文通过对比由计算结果绘制的气门理论升程曲线和实际升程曲线，即该配气机构的动力学性能分析，对该发动机配气机构节能型了评估。本课题研究中，分析结果表明，随着转速的提高，气门渐渐会出现“飞脱”和“反跳”，但在发动机额定转速以及较高转速的情况下，该发动机没有出现上述问题。因此，该配气机构总体性能良好，设计符合要求。如果需要对其他顶置凸轮式配气机构做动力学分析，由于分析方法大致相同，可直接使用本课题中的动力学通用程序，以及单自由度模型，得到分析结果。因此本课题的研究拥有比较好的通用性。配气机构的动力学分析的一个有效的途径便是单自由度动力学模型，该模型已经能够成熟的被龙格库塔法的迭代算法求出。通过对该机构建立单自由度系统的动力学模型，由已知气门理论运动规律求气门实际运动规律的动力学计算，得到了气门实际运动规律曲线，经过对曲线进行分析后，对该发动机配气机构进行评估，从而保证所选择的设计参数和凸轮型线能使配气机构平稳正常工作。因此这个模型还兼具优化设计的作用。但是顶置凸轮轴式配气机构最大的特点便是摇臂比是变化的，并且变化的幅度比较大，在不同的转角下摇臂和系统的刚度都是随着摇臂比的变化而变化的，将系统的刚度当作一个定量，影响了分析的精确性。实践表明单自由度模型对气门机构的运动情况做出的分析基本上能达到工程精度要求，但更详细的情况，诸如传动链的飞脱究竟发生在哪一环节上，各接触副和气门的动载荷究竟有多大以及气门弹簧的震颤问题等等，是无法在单自由度模型中得到，因此往往在配气机构基本方案确定后，进一步用多自由度模型进行一些补充验算。31\n重庆大学本科学生毕业设计（论文）展望七、展望配气机构在发动机中控制着换气过程，其运动的设计直接影响内燃机的可靠性及性能。发动机配气机构应保证充气系数尽可能高，开启和关闭进排气门的时刻精准，各气缸换气良好，合理。四冲程发动机大多采用凸轮——气门式配气机构，工作可靠，运转稳定。本次课题针对顶置凸轮式配气机构做了运动学分析计算，建立了单自由度模型进行动力学分析，计算工作量不大。单自由度模型是经常使用在气门机构设计方案初选阶段，诸如传动链的飞脱究竟发生在哪一环节上，弹簧的颤振是否会导致过大的应力或某些簧圈相碰，各接触副和气门杆的动载荷有多大等等，就无法由单自由度模型分析中得到[[]DavidJW.OptimalDesignofHigh-speedValveTrainSystem.SAE94502,1994]。因此往往在气门机构基本方案确定后，对于一些特殊工况下尤其是在接近极限情况下，机构的响应究竟具有何种特性，还需要我们进一步分析验证，这就需要我们用到多自由度模型进行一些补充验算，虽然多自由度模型分析得到的结果更加准确，更具有参考性，但是由于自由度数多，从模型分析所得的运动微分方程就越为复杂，求解复杂的微分方程是一项十分繁琐的工作，目前并没有哪一款软件或哪一种算法可以被通用在复杂的微分方程求解问题上，所以对配气机构多自由度模型的研究还有待进行。对于配气机构，相对于传统的凸轮控制式配气机构，现在有一种比较新型的设计——无凸轮发动机配气机构（camlessenginevalvetrains）。正如其名，该配气机构中没有凸轮，与目前比较偏重机械方面的设计来说，该机构利用了电磁控制，结构更小巧、简单，并且控制精确，在控制气门正时和气门升程方面有着更好的效果。该发动机缸头设计减少了20％的燃料消耗和污染物排放量。该机构由智能气门驱动系统（SVA）代替了由凸轮带、凸轮轴和液压凸轮配件组成的传统的机械发动机气门机构。在无凸轮发动机中，每个发动机气门由一对置于缸盖上表面上的气门驱动器直接控制。发动机基本排量不变的前提下进行配气机构的优化设计师当今发动机配气机构设计的发展趋势所在，提高各方面性能指标。为了更大范围内提升内燃机动力指标、经济指标和生态指标，早先有多气门配气机构，随后又发展有可变气门运动配气机构，而传统配气机构相对于可变气门运动配气机构提升内燃机动力指标、经济性的能力这方面的能力就要逊色许多[[]钱向明,熊永森.浅析柴油机配气机构的发展现状[J].农业装备与车辆工程,2009(4):44-46.[15]刘臭陈国良，陈鹰，顶置凸轮轴式配气机构凸轮轮廓的数值计算，内燃机学报,2002,(3):283-286[16]吕林王勇波，车用发动机配气机构运动学和动力学分析，武汉理工大学学报，2006.12[17]吴兆汗车辆发动机设计，国防工业出版社，1982[18]尚汉翼内燃机配气凸轮机构设计与计算，上海：复旦大学出版社，1998]。对配气机构的研究还有很大的发展空间，只有不懈的创新、追求和探索，才能使之不断改进，不断完善，更加广泛的被使用。31\n重庆大学本科学生毕业设计（论文）参考文献致谢光阴似箭，岁月如梭，匆匆大学四年即将走到尽头，是结束亦是开始。这里有我熟悉的教学楼、实验室，有我热爱魅力校园风光。不经意间，在这熟悉的校园中，我度过了人生中最为珍惜的一段时光。大学四年期间，有勤奋时的朝暮勤勉，有懈怠时的慵懒无为。但始终感受到整个时代的进步与科学技术发展的迅速，内心时而会不断激励自己，将自身埋在图书馆浩瀚的书海里，以提升个人素质，思想境界，学习上也不断勤勉自己，努力专研专业学术知识，无奈大学四年匆匆过去，心中不免深深叹息韶光易逝。毕业前最后一次郑重的检验便是毕业论文，这是对自我所学知识，科研能力的此全面检验。从课题的开展到论文完成的整个过程中，最重要的是给予我指导的导师阮登芳教授，另外本专业的许多同学也给予我不少软件分析上的帮助，在这里我要由衷的感谢他们，谢谢你们!首先，我要感谢学校和学院为我们建设了这样一个美丽的校园学习环境，和专业的教学实验室实验器材，让我们可以理论联系实际，学习更加深入。同时我还要感谢本学院所有的老师们，感谢他们在这四年来对我的悉心教导。师之所存，道之所存也。正是有了他们这一群教学严谨细致、工作一丝不苟的教育工作者，科学知识，人生哲理才得以最大化的延续。他们的直言不讳，使我在今后的人生中收获甚多。其次，对于我的指导老师阮登芳教授我要特别感谢一次，开题，译文翻译，国内外文献查阅，软件分析，解题整个论文的撰写工作无不透露着老师的心血，整个过程都有老师的陪伴，这份良苦用心让我感动不已。在老师的悉心指导下，艰苦奋斗，查阅各种文献方法这篇论文才算顺利完成。在这个过程中，导师指点我查找分析文献资料，与我共同分析课题分析结果的合理性，一同探究期间计算程序出现的问题。每一次的帮助我都十分受用，在此对老师致以深深的谢意。最后，本专业同学以及其他专业同学在我整个毕业设计过程中都有或多或提供指导和帮助，这让我深深感受到同学间的那种情谊，在这论文完成的时刻，我要谢谢你们，谢谢你们的精神鼓励和技术支持，风里雨里我们都曾一起走过，让我在整个过程中感到无比幸福！感谢所有关心和帮助过我的老师、同学、朋友，谢谢你们！在论文完成之际，我要向我的指导老师致以最衷心的感谢和深深的敬意!最后的最后，衷心地感谢在百忙之中抽空评阅本文和参加答辩的各位专家、教授!31\n重庆大学本科学生毕业设计（论文）附录附录A：matlab运动学分析程序已知凸轮升程--转角（y，α）数据求气门升程—转角（h，αr）matlab程序实现la=39.61;lc=30.28;ro=12.5;rs=30;ymax=6.532;p=10;o=0.05;lb=33.20;z=1.65;yys=gradient(ys)./gradient(as);fori=1:length(ys)y=ys(i);a=as(i);yy=yys(i);R=(yy.^2)+(ro+rs+y).^2;w=acosd((la.^2+lc.^2-R)/(2*la*lc));rr=R.^0.5;e=atand(yy/(ro+rs+y));k=acosd((lc.^2+R-la.^2)/(2*lc*rr));Rd=ro+rs+ymax;kd=acosd((lc.^2+Rd.^2-la.^2)/(2*lc*Rd));ar1=a+e+(k-kd);ar2=a+e-(k-kd);wo=acosd((la.^2+lc.^2-(ro+rs).^2)/(2*la*lc));bo=asind((p+z+o)/(lb));b=bo-(w-wo);h=lb*(sind(bo)-sind(b))-o;shuchu(1,i)=w;shuchu(2,i)=ar1;shuchu(3,i)=ar2;shuchu(4,i)=h;end[FileName,PathName]=uiputfile({'*.xls','EXCELFiles(*.xls)';'*.*','AllFiles(*.*)'}liming.xls');xlswrite([PathName,FileName],shuchu);liming.xls31\n重庆大学本科学生毕业设计（论文）附录附录B:动力学分析计算基本程序k=11700;nc=12;ps0=205.8;ks=41.356;beta1=17.73;M=0.695;det=0.5;y=0;yy=0;A=zeros(1,386);B=zeros(1,386);C=zeros(1,386);forn=1:1:386k1=(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;C(n)=k1;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);k3=(k/(36*nc^2*M))*Y(2*n)-ps0/(36*nc^2*M)-(k+ks)/(36*nc^2*M)*(y+det/2*yy+det^2/4*k1)-beta1/(6*nc*M)*(yy+det/2*k2);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);y=y+det*yy+det^2/6*(k1+k2+k3);A(n)=y;yy=yy+det/6*(k1+2*k2+2*k3+k4);B(n)=yy;endplot(x,A);holdon;plot(x,B*20)plot(x,C*20)31\n重庆大学本科学生毕业设计（论文）参考文献参考文献31