关 键 词：

机构 动态 仿真 说明书 仿单 matlab

资源描述：

内容简介：

计算机工程学(2002) 18: 109115所有权和版权2002年伦敦Springer-Verlag出版社对基于固体参数模型的汽车部件的结构最优化M. E. Botkin美国，沃伦，密歇根州，冈比亚研发中心摘要：参数塑造是用于修建汽车前部结构概念的某些模型的，而汽车前部结构概念则是利用碳纤维合成材料和对应的成型过程。超轻量级铝身前部结构被再设计成包括所有合成材料的前部结构。有两个概念被人们研究，这两个概念都把汽车前部结构描绘成一个固定的外型。闭合部分是由内外两块构成的。因此，参数塑造在构造和改善模型以用于最优化概念研究方面是一个很有用的工具。这些模型可以迅速被建立，拟订尺寸和定位尺寸尤其能用以做必须的调整，把各个主要部件安装在一起。然后参数模型必须被连接成一个几何体模型以获得一个表面滤网。这样 结构优化输入数据才可能迅速地被以参数模型为基础的有限元素模型创造从而开始交易学习。参数塑造与结构优化连接的完整过程被用于轻量级前部结构的设计研究中 。人们研究一些碳纤维化合物质以尽可能地大量减少某些比过去成本低的碳纤维材料。结构优化用于组合结构和固定的前部结构设计的对比。通过各种材料和合成铺垫对八个案例进行了研究。预计钢用量节省45%-64%。最合理的设计包括一个由相对低价的切割碳化纤维和机织碳化纤维合成物，在震动区域顶部使用了20毫米的轻质木材。这个设计最大厚度达到7毫米，可大量减少钢用量，大约达到62%。书信和选刊可寄至：美国沃伦密歇根州冈比亚研发中心，主要研究工程师Mark E. Botkin,邮箱编号: 480- 106-256, 30500 MI 48090-9055。电子邮件: mark.e.botkin_ .关键字： 自动化配合; 汽车; CAD 塑造; 综合; 有限元素塑造; 优化 1. 介绍结构被显示在图1 是几个先进材料组成的轻量级身体。设计这个身体被提供了在Prsa 1 并且是执行作为PNGV 2个的部份。用主要设计目标是重量减少, PNGV 身体近似地接近 70% 比克莱斯勒触毛钢质车身它是随意替换。主部位主要由碳纤维外壳的三明治形状的结构和蜂窝状铝结构构成。少量Kevlar 并且被使用了某一Nomex 核心。但是, 前面结构的负荷成份是铝。本文描述一个项目再设计所有碳纤维结构。这个概念被考虑了。概念被描述在本文里能被描绘对象一起形成闭合部分路轨的板料。概念被塑造，使用了参数塑造特点 UG 3 , 并且被显示在下个部分。滤网自动地被创造了使用UG 的假定容量。 图1. 碳纤维身体。110 M. E. Botkin 由于滤网被创造了根据了固体几何学模型, 分开的物产地区是自动地创造为各张表面面孔。这做优化数据创作逐步采用更加容易地比早先研究 4 。譬如 Patran5 可能自动地创造优化数据从适当地被创造的结构分析数据。本文同时描述了Nastran5 的用途解答200 执行设计研究为综合设计观念为前期结构。 2. 参数塑造参数塑造的过程从二维空间的参量横断面开始，过去经常创造固体模型通过挤压一条直线或清扫沿曲线。前线结构被显示在图1 被塑造了, 使用这参数过程, 作为保税一起的综合板料。模型组成, 主要作为三主要组分; 上部路轨, 更低的路轨, 和震动塔。其它次要盘区栓这些三组分一起。为例, 图2 展示外面片断的参数横剖面上部路轨。维度被显示是可变物那可能由设计师改变或者使用塑造的节目或自动地使用优化规程。在可变物被修改之后, 整个固体模型自动地被更新。图 3 个展示创造的过程三维几何从二维部分。作为部分被清扫沿曲线, 立体几何是自动地引起。当任何部下信息被修改, 如部件尺度, 曲线定位等, 立体几何被 图2. 上部路轨外面片断。自动地更新。相似的操作被使用为更低的路轨和震动塔获得完成将被使用的参数固体模型在这设计研究中, 被显示在图4 。图4 代表三个片断综合结构(鞋帮路轨、更低的路轨, 和震动tower/apron) 被加入与胶粘剂一起。为这设计研究的目的, 参数塑造只被使用了作为一个方便方法创造有限元素模型。最后, 然而, 一种更加全面的设计方法会准许参数模型是整体部分优化过程参数维度有时间被使用当设计可变物。这是可能的, 到大规模范围, 以版本UG16 和以后和被展示了在Botkin 6 。 3. 有限元素分析3.1. 滤网世代滤网被显示在图5 被创造了从UG Scenariousing fully-automatic 四方形捕捉世代能力使用名词性的词10 毫米元素大小。这是一个推进的前面技术特点那些发现了在最商业塑造节目。但是, 由于紧的综合化在参数modeler 和mesher 之间, 分明表面地区被显示在图4 被维护当物产地区在滤网和, 依照被显示在图6, 能自动地被指定作为设计可变物在优化模型。模型显示在图5 由11,710 个结组成, 12,225 元素, 58,877 程度自由, 和23 不同的元素物产。除11488 之外轰击元素, 737 个CBAR 元素被用于代表被使用加入的黏着性债券综合片断一起。 3.2. 分析模型图5 并且显示装载和边界条件使用在这项研究。最优化研究的目的是设计一个与图一中铝结构坚硬程度一致的合成材料的前部结构。图1中结构的扭转力坚硬程度是 1 10,246 N.m/deg 。对震动区域的分析发现它的坚硬度比垂直方向17,963 N/mm在侧向方向12,941 N/mm的总体部位高很多。111 车身组分的结构优化 图3. 逐渐变细了上部路轨。 图4. 最后的前面结构模型。 图5. 最后的滤网从UG/Scenario 。 图6. 物质物产地区。 图7. 结合物模型。那对应对055 毫米的限制的偏折为 1000 N 装载被显示在图5 和077 毫米为a 相似的侧面1000 N 装载。 112 M. E. Botkin 3.3. 联合建模依照被注意及早, 三单独综合片断将由对胶粘剂的用途加入。图7 展示塑造的方法被选择为这项研究。元素列沿边缘各个片断由对射线元素的用途加入为哪些物产是坚定的从胶粘剂的物质物产(被显示在表里 1) 射线物产是结合物的区域考虑对结A 和B 贡献(参见图 7) 和转动惯量。这些价值是50 mm2 和416.67 mm4, 分别为， 承担a 10 毫米(正方形的) 有名无实的元素大小。 表1 ，合成道具 4 4. 优化研究Nastran 解答200 被使用执行八材料的各种各样的组合的专题研究。物质物产被给在表1 。所有缺省Nastran 优化参量 7 被使用了除了近似优化被选择是凸面线性化的技术 8,9 。依照被注意在Botkin 4 , 凸面线性化, 是最保守略计方法可利用在Nastran, 被使用了由于困难接近失败限制。4.1. 宗旨和限制作用研究的目标函数是极小的。设计是, 依照前面提到, 正如当初提到的。限制被强分别加了给侧向和垂直的偏折0.077 和0.055 毫米。4.2. 设计可变物依照被早期提起的, 设计可变物能自动地被创造从物产地区显示在图6 。图8 显示盘区从被显示当设计研究的Patran 5 工具被使用。虽然有23 设计可变物, 不是所有被显示在一个唯一盘区。原始价值能被修改使用这个盘区和然后解答 200 数据输入(DESVAR & DVPREL1 卡片图象) 自动地被创造。它应该是 113 车身组分的结构优化 。 图8. 设计可变物盘区。指出厚度价值被显示在图 8 是10 毫米的原始价值为优化浇铸。30 毫米价值为PCOMP 词条包括一个10 毫米浮筒木筏核心和二10 毫米综合面对板料。虽然优化使用的方法在这中设计研究, Nastran 解答200, 是可胜任对待综合使用角度可变物- 象被展示了在Botkin 4 - 唯一固定的有机材料, 90个地方被使用了在努力减少制造费用。4.3. 设计观念选择虽然图4 指设计观念, 参数固体模型可能被重视作为a 几何概念与综合相对设计观念。这设计研究的目标是使用低成本碳纤维产品譬如砍纤维和较不昂贵的大拖曳产品。估计的物产为这些综合可能是看见在表1 与所有材料一起物产使用为这项研究。最低界面在所有设计可变物被设置了在被考虑的2 毫米是可能当前被铸造的最稀薄的零件使用液体造型过程。箱子被分离入三个小组由最高界面: 15, 10, 和 _ 10 毫米。这些厚度代表连续研究发生在努力对发现a 设计以一种适当地小最大厚度。所有案件和厚度各自的表现在总结在图9 和10。它会是最中意的为能力存在会确定最宜综合概念而不是必须比较数选择了概念。这样能力的存在不为创作者所知道。 图9. 案件许多总结。 图10. 一致设计可变物。4.4. 15 毫米最高界面案例 1第一案件将是那使用所有被砍碳材料没有三明治建筑。这会被认为最少昂贵的案件使用最便宜的材料和有核心。核心增加了材料加工成本和铸造难度。正如所有案件那里是23个设计可变物和位移限制将被安置在震动塔的上面维护目标僵硬提前先被显示。A 初步分析表示, 最初的设计违犯位移限制大约33% 。图9展示这个情况的主体历史.。虽然最初的僵硬目标被违犯了近似地 13.6 公斤, 优化器能发现设计那满足了僵硬目标在6.2 公斤, 或64% (6.2 kg/64%) 许多减少估计钢大量。考虑这是一相对地降低表现材料, 这些结果是卓越的。图10 显示设计易变的发行这个案件。所有设计可变物的最初的厚度是10.0 毫米。数可变物被终止 114 M. E. Botkin 在被考虑2.0 毫米的他们的最低界面是可能合理地是的最小的厚度铸造使用树脂调动造型(RTM) 或结构反应注坯模型(SRIM) 。几可变物, 然而, 任意地结束了在 15 毫米选上的最高界面。它应该是针对性的在那之外终止在他们的鞋帮的可变物区域是更加敏感的(在数学感觉) 对限制比大量。剩下的人可变物对大量是敏感的比限制和因此提供机会为许多减少。一种15 毫米厚度被考虑是相当不切实际的价值和导致a 材料在那些地区需要的结论或者更高执行或三明治建筑, 带领案例2 。案例 2这个案件增加一个10 毫米厚实的浮筒木筏核心1.0 g/cm3 (6.5 lb/cf) 在震动的上面的密度塔, 即, 可变物7 和10 至13 依照被显示在图6 。由于增加的僵硬, 最初设计在这种情况下是7.8% 僵硬比目标价值。图9 显示这个案件的许多历史。虽然最初的设计是更加僵硬的, 最初大量更高归结于核心的重量和增加的面孔板料。另外, 最初的设计违反了方案1 (6.62 kg/61.3%) 。依照被展示在Botkin 4 当外壳厚度是非常大, 三明治建筑不能是有效。虽然几个变数被迫使他们的上界,最其它变数是与以防方案1相比更小。案例 3这个案件修改表面皮肤在三明治盘区(震动塔上面) 对一更高执行材料。这将是里的被编织的材料在物产被显示在表1 。它被承担这材料由低成本将被做 50K 拖曳碳纤维。图9 显示结果这个案件。可得出，最终质量只稍微好了一点(6.65 kg/61.2%)，厚度也没有如愿降低。即使弹性模数在织物上大得多，剪切性能也还是太低了。4.5.10 毫米最高界面以下案件显示减少的作用最高界面在厚度对10 毫米。它被考虑那15 毫米是过份地厚实的是合理的建筑。如同有人可能期待优选的大量增加在所有案件。案例4 这个案件维护所有切好的碳但与一个核心。(9.35 kg/45.4%) 案例5这个案件使用被编织的材料在震动塔区域(9.09 kg/46.9%) 。 案例6。 这个案件使用两层机织材料。外层0-90，内层45，以改善机织材料的剪切性能。这种构造是要明显降低其他两个案例中的质量。 (7.95 kg/53.6%) 。 4.6. _ 10 毫米最高界面 案例7这个案件是案例6的扩展。密度是每立方厘米1.5克20毫米，材料厚度上限7毫米。MTC 身体并且使用了一个20 毫米核心但铝蜂窝。这个案件导致大量6.55 公斤或61.8% 许多减少结束钢。这是a 非常合理的设计以最大厚度在实用考虑之内的范围。 案例8这个案件增加三明治建筑(一10 毫米浮筒木筏核心) 设计易变的地区17 和18 显示在图6 。案例7 的结果被表明这些地区是在他们的最高界面。此外, 依照被显示在图1, 震动塔并且被使用三明治建筑在侧面墙。正象一样案例6 和7 被编织的材料被使用了与同样放置。这个案件的结果非常是鼓励因为6 毫米一个最高界面是获得。虽然优选的大量大于方案7 (6.95 kg/59.4%), 如果最大厚度是问题, 这个案件也许是更加中意的。 4.7. 力量限制 案例9在这种情况下, 二种负荷状态增加来包括力量限制对僵硬设计(案例8) 。第一负荷状态代表平均障碍装载从35 哩/时(56.3 km/hr) 撞击在汽车的躯体成分较低栏杆的上面(10,000 N)和115结构上最优化应用(44,000 N),在每一方面。第二装货情况代表一个常用3g 坑洼装载(6675 N) 垂直应用了震动塔。约束被把强加给压力把碳和失败索引缩短为编织材料是4。从拉伸和耐压强度切好的碳是不同的, 更多保守的极限150 Mpa 被使用了作为重音限制。为失败索引(Hoffman), a 价值的1.0 被使用了作为限制极限。这个案件的结果被显示作为案例9 在图9 和10 。由于严格的限制是非常严格作为这低僵硬材料, 严格限制没有充当在的一个非常重要角色设计(7.1 kg/58.6%) 与案例8 比较了。它是仍然, 但是, 重要确定如果严格只是设计现实的。 4.8. 案件总结案例7 看来是最佳的设计。这个案件使用较不昂贵的切好的碳在所有区域除了二震动塔的上面被有机物的或被缝的席子材料层数被使用契合以一个20 毫米浮筒木筏核心。这型增强可能并且被获得从更加便宜, 大拖曳碳纤维和被显示有相似的物产。许多减少是61.8% 钢结构和一68% 相比为 PNGV 前面结构, 依照被指出及早, 被制造了从非常昂贵航空航天象材料和过程。5. 总结和结论 参数模型和结构最优化被用于研究设计轻量级的主体前部结构。无需手动数据准备，一个紧密连接的参数结构设计过程被描述了。参数模型被用于建造汽车前部结构概念的几个模型，汽车前部结构概念利用碳纤维合成材料和相应的成型步骤。通常有两个概念被研究，本文中只提出了一个，这个概念把汽车前部结构看作一个紧密的外形集合体。闭合部分是由内外两块构成的。超轻量级PNGV前部结构被再设计成包括所有合成材料的前部结构。人们通过研究几个碳纤维合成材料来尽可能地大量降低比过去低成本的碳纤维。结构优化用于组合结构和固定的前部结构设计的对比。通过各种材料和合成铺垫对八个案例进行了研究。预计钢用量节省45%-64%。最合理的设计包括一个由相对低价的切割碳化纤维和机织碳化纤维合成物，在震动区域顶部使用了20毫米的轻质木材。这个设计最大厚度达到7毫米，可大量减少钢用量，大约达到62%。同时人们还发现增加的力度限制对设计没有太大影响。参考文献： 1. Prsa, J. (1999) Hybrid material automotive structuredevelopment: Phase 2, 99IBECA-27, InternationalBody Engineering Conference, Detroit, MI2. Partnership for a New Generation of Vehicle HomePage. /pngv/3. UNIGRAPHICS online documentation, Version 15./docs/unigraphics.html4. Botkin, M. E. (2000) Modeling and optimal design ofa carbon fiber reinforced composite automotive roof,Eng. with Comput. 16(1), 16235. MSC PATRAN online documentation, Version 70.5(1999) The MSC. Software Corporation6. Botkin, M. E. (2000) An Assessment of DesignOptimization in Unigraphics Version 16. GM R&DPublication R&D-91007. Moore, G. J., Design Sensitivity and Optimization.MSC/NASTRAN Users Guide8. Starnes, J. H., Haftka, R. T. (1979) Preliminary designof composite wings for buckling, strength, and displacementconstraints. J. Aircraft 16(8), 5645709. Fleury, C., Braibant, V. (1986) Structural optimization a new dual method using mixed variables. Int.J. Numer. Meth. Eng. 23(3), 4094282002) 18: 1091152002 E. &D to of an An to an as a of an to be a to in be to be as in to to in to on to of of to be in to of 564% of a of a 20 mm in of a mm of 2%. &D 48030500 d., 80901 is a of 1 as a of 2 it to is a of as as of a an be as 3, 1. . E. a 4. 5 500 to to or a 1 as is as As an 2 of be by or a is is of As is of is 2. of to be 4. a of as a a to be an to be as is to a 6 65 a 0 is an of in of 4 as in as 6, be as in 5 is 1,710 2,225 58,877 3In 1488737 to in of is to an to as 1. of 1 to 1 10,246 Nm/It of is of 7,963 N/mm in 3. 4. 5. G/6. 41 N/mm in to a 055 mm 5 077 mm 000 N . E. 7. to be by of of by of of of to to (of 00 . 4 (kg/(6000 701 78 143 178 143 502000 200 30 370 530 370 75000. 800. 630. 840. 34. 110. A 9395 4937. 842. 55. 96. 55. 96. 4070. 159 2240. 108. 00 to of of . 7to 8,9. As 4, is of of of as on be 6. 5 is is 3 on a be 00 is It 8. 0 mm 0 mm a 10 mm 0 mm in 00, is as 4 90, in an to 4 is to as a be of as as to a of is as of on mm is be by 15, 10, in an to a of 0,be a an to of a is to 9. of 10. 15 be of no be to be no to to As 3 a be at of 3%. of at to a .2 4%(6.2 4%) is a 0 of 0.0 of . E. .0 mm is be be at 5 It be at in a to to A 15 mm is be a to in be or of a a 10 mm .5 lb/to of 0 3 as 6. to in of to of In (4 0 of in of to a be a in . It is be As be is as is 10 of on 0 It is 5 mm is to be a As in (in of 90 1100645in to of of 0 mm (. is an of a 20 .5 g/.5 lb/an on a 20 mm of in a is of a 10 to 7 86. to be at as 1, in as of in an mm is (, if an be In to ). a 35 at 10,000 N) 44,000 N), a g 6675 N) at on 4 50 as .0 as of 0. a in 7.1 to . It is to if to be in of of or is a 20 mm be to a 68% as to on A no is to of an in an of an to an to of be in of 5% 4%of a of a 20 in of Thisde 毕业设计（论文）外文翻译 题目 对基于固体参数模型的汽车部件的结构最优化 专 业 名 称 机械设计 与 制造及其自动化 班 级 学 号 078105118 学 生 姓 名 施小康 指 导 教 师 朱保利 填 表 日 期 2011 年 3 月 20 日 南昌航空大学科技学院学士学位论文目 录1 绪论(2)1.1引言(2)1.2平面连杆机构及杆组概述(2)1.3 进行杆组系统仿真的意义(3)1.4 仿真软件的发展状况与应用(3)1.5 MATLAB概述(3)2 7R六杆级机构运动学仿真(5)2.1 曲柄原动件运动学分析(5)2.2 6R级杆组运动学分析(6)2.3 7R六杆级机构MATLAB仿真积分模块初值的确定(11)2.4 7R六杆级机构运动学仿真模型及结果(16)3 7R六杆级机构动力学仿真(23)3.1 曲柄原动件动力学数学模型的建立(23)3.2 6R级杆组动力学数学模型的建立(25)3.3 需要引用的函数(30)3.4 7R六杆级机构运动学仿真模型及结果(32)4 结论(39)参考文献(40)致谢(41)1 绪论1.1引言大学的四年生活，通过老师的讲解和我自己的学习，我收获了很多，我也深深的喜欢上了机械这个行业，对机械加工和制造方面尤为感兴趣，我觉得通过自己的努力和思考来改变工艺规程来提高生产效率，提高经济效益很有成就感。我所研究的课题就是给了这样的机会我可以通过我的努力来优化工艺规程，提高经济效益。此次毕业设计，是在我们学完了机械制造工艺学、工艺装备设计等课程，进行了生产实习之后，进行的一个重要的实践性环节。这要求我们把所学的工艺理论和实践知识，在实际的工艺、夹具设计中综合地加以运用，这有助与提高了我们分析和解决生产实际问题的能力，为以后从事相关的技术工作奠定的基础。1.2平面连杆机构及杆组概述平面连杆机构是将各构件用转动副或移动副联接而成的平面机构。最简单的平面连杆机构是由四个构件组成的，简称平面四杆机构。它的应用非常广泛，而且是组成多杆机构的基础。全部用回转副组成的平面四杆机构称为铰链四杆机构。仅能在某一角度摆动的连架杆，称为摇杆。对于铰链四杆机构来说，机架和连杆总是存在的，因此可按照连架杆是曲柄还是摇杆，将铰链四杆机构分为三种基本型式：曲柄摇杆机构、双曲柄机构和双摇杆机构。在实际机械中，平面连杆机构的型式是多种多样的，但其中绝大多数是在铰链四杆机构的基础上发展和演化而成。如曲柄滑块机构、导杆机构等。任何机构都是由原动件、机架和从动件构成的系统。由于机架的自由度为零，一般每个原动件的自由度为1，且根据运动链成为机构的条件可知，机构的自由度与原动件为应相等，所以，从动件系统的自由度数必为零。机构的从动件系统还可以进一步分解成若干个不可再分的自由度为零的构件组合，这种组合称为杆组。设n表示活动构件数，PL表示低副个数，根据n的取值不同，村级可分为级杆组和级杆组。其中级杆组分为5种：RRR级杆组、RRP级杆组、RPR级杆组、PRP级杆组以及RPP级杆组。任何机构都可以看作是由若干个基本村级依次联接于原动件和机架而构成的，这就是所谓机构的组成原理。通常，把由最高级别为级杆组的基本杆组构成的机构称为级机构；把最高级为级杆组的基本杆组构成的机构称为级机构。1.3 进行杆组系统仿真的意义系统仿真是建立在控制理论、相似理论、信息处理技术和计算机初等理论基础之上的，以计算机和其他专用物理效应设备为工具，利用系统模型对真实或假设的系统进行试验，并借助于专家的经验知识、统计数据和信息资料对实验结果进行分析研究，进而做出决策的一门综合的实验性学科。仿真技术是分析、研究各种系统，尤其是复杂系统的重要工具。随着机械行业的迅速发展，对研究、设计的机械设备越来越复杂，用于制造各种零件的材料价格越来越昂贵，不可能每一步都采取试制再修改的方法进行设计，采用仿真的方法可以在一定程度上克服这种不足的不足，降低研究成本，提高效率。而连杆机构作为常见的传动机构，对其进行运动学和动力学仿真，建立起基本杆组模块的仿真模型，无疑对日后的设计大有裨益。一般机构的运动分析，使用Quik BASIC语言或Fortran语言编写程序进行计算，其缺点“透明性”差，修改麻烦等而用MATLAB对机构进行运动仿真，利用MATLAB的simulink仿真模型的数据可视化的特点，就可以很容易观察到运动参数是如何变化的，极其简便同时，用MATLAB建立和修改仿真模型具有方便、快捷、很容易扩展等优点MATLAB仿真求解器提供很多解不同微分方程的方法，可以根据不同的微分方程类型选择相应的求解方法机构的动力学分析，由已知工作阻力，求出运动副的约束反力和驱动力（或力矩），为选择和设计轴承和零部件强度的计算及选择原动机提供理论依据。1.4 仿真软件的发展状况与应用早期的计算机仿真技术大致经历了几个阶段：20世纪40年代模拟计算机仿真；50年代初数字仿真；60年代早期仿真语言的出现等。80年代出现的面向对象仿真技术为系统仿真方法注入了活力。我国早在50年代就开始研究仿真技术了，当时主要用于国防领域，以模拟计算机的仿真为主。70年代初开始应用数字计算机进行仿真4。随着数字计算机的普及，近20年以来，国际、国内出现了许多专门用于计算机数字仿真的仿真语言与工具，如CSMP，ACSL， SIMNOM， MATLAB/Simulink， Matrix/System Build， CSMP-C等。1.5 MATLAB概述MATLAB是国际上仿真领域最权威、最实用的计算机工具。它是MathWork公司于1982年推出的一套高性能的数值计算和可视化数学软件，被誉为“巨人肩上的工具”。MATLAB是一种应用于计算技术的高性能语言。它将计算，可视化和编程结合在一个易于使用的环境中，此而将问题解决方案表示成我们所熟悉的数学符号，其典型的使用包括:.数学计算.运算法则的推导.模型仿真和还原.数据分析，采集及可视化.科技和工程制图.开发软件，包括图形用户界面的建立MATLAB是一个交互式系统，它的基本数据元素是矩阵，且不需要指定大小。通过它可以解决很多技术计算问题，尤其是带有矩阵和矢量公式推导的问题，有时还能写入非交互式语言如C和Fortran等。MATLAB的名字象征着矩阵库。它最初被开发出来是为了方便访问由LINPACK和EISPAK开发的矩阵软件，其代表着艺术级的矩阵计算软件。 MATLAB在拥有很多用户的同时经历了许多年的发展时期。在大学环境中，它作为介绍性的教育工具，以及在进阶课程中应用于数学，工程和科学。在工业上它是用于高生产力研究，开发，分析的工具之一。Simulink概述 Simulink是用于仿真建模及分析动态系统的一组程序包，它支持线形和非线性系统，能在连续时间，离散时间或两者的复合情况下建模。系统也能采用复合速率，也就是用不同的部分用不同的速率来采样和更新。Simulink提供一个图形化用户界面用于建模，用鼠标拖拉块状图表即可完成建模。在此界面下能像用铅笔在纸上一样画模型。相对于以前的仿真需要用语言和程序来表明不同的方程式而言有了极大的进步。Simulink拥有全面的库，如接收器，信号源，线形及非线形组块和连接器。同时也能自己定义和建立自己的块。模块有等级之分，因此可以由顶层往下的步骤也可以选择从底层往上建模。可以在高层上统观系统，然后双击模块来观看下一层的模型细节。这种途径可以深入了解模型的组织和模块之间的相互作用。2 7R六杆级机构运动学仿真2.1 曲柄原动件运动学分析2.1.1 曲柄原动件运动学数学模型的建立图1 曲柄的复数坐标系如图1所示，在复数坐标系中，曲柄AB复向量的模rj为常数、幅角j为变量，通过转动副A与机架连接，转动副A的复向量的模ri为常量、幅角i为常量，曲柄AB端点B的位移、速度和加速度的推导如下： （2.1）将方程2.1两边对时间t求两次导数得：（2.2）由式2.2写成矩阵形式有 （2.3）2.1.2 曲柄MATLAB运动学仿真模块M函数根据式(2.3)编写曲柄原动件MATLAB的M函数如下：function y=crank(x)%x(1)=rj 曲柄杆长%x(2)=thetaj曲柄与水平方向夹角%x(3)=dthetaj曲柄角速度%x(4)=ddthetaj曲柄角加速度%y(1)=ReddB转动副B加速度实轴分量%y(2)=ImddB转动副B加速度虚轴分量ddB=x(1)*x(4)*cos(x(2)+pi/2)+x(1)*x(3)2*cos(x(2)+pi); x(1)*x(4)*sin(x(2)+pi/2)+x(1)*x(3)2*sin(x(2)+pi); y=ddB;此函数模块用于计算转动副B的加速度的水平分量和垂直分量输入参数为曲柄的长度、角位移、角速度和角加速度；输出参数为曲柄端部（转动副B）的加速度的水平分量和垂直分量。2.2 6R级杆组运动学分析2.2.1 6R级杆组运动学数学模型的建立图2 6R级杆组的位置参数如图2所示，在复数坐标系中，由3个外转动副(B,C,D)和3个内转动副(E,F,G),4个构件(BE,CF,DG和EFG)构成1个6R级杆组，构件BE,CF,DG的长度分别为ri，rj，rk构件EFG的3个边为e，f，g，方向如图所示，规定所有复向量与实轴正方向逆时针夹角为，并用相应的下标来区别，用 B,C,D,E,F和G分别表示该转动副的复数坐标，则各个构件的运动参数推导如下： （2.4）（2.5）将式（2.4）代入式（2.5）并整理得： （2.6）将式(2.4)至式(2.6)三式合并成矩阵得： （2.7）将式(2.7)展开整理得：（2.8）点E，F，G的加速度分别为 （2.9）2.2.2 6R级杆组MATLAB运动学仿真模块M函数根据式(2.9)编写6R级杆组MATLAB的M函数如下：function y=R6ki(x)%x(1)=riBE杆长%x(2)=rjCF杆长%x(3)=rkDG杆长%x(4)=eFG杆长%x(5)=fGE杆长%x(6)=gEF杆长%x(7)=theta-iBE杆与水平方向夹角%x(8)=theta-jDF杆与水平方向夹角%x(9)=theta-kDG杆与水平方向夹角%x(10)=theta-eFG杆与水平方向夹角%x(11)=theta-fGE杆与水平方向夹角%x(12)=theta-gEF杆与水平方向夹角%x(13)=dtheta-iBE杆角速度%x(14)=dtheta-jEF杆角速度%x(15)=dtheta-kDG杆角速度%x(16)=dtheta-eGE杆角速度%x(17)=dtheta-fGE杆角速度%x(18)=dtheta-gEF杆角速度%x(19)=ReddB转动副B加速度实轴分量%x(20)=ImddB转动副B加速度虚轴分量%x(21)=ReddC转动副C加速度实轴分量%x(22)=ImddC转动副C加速度虚轴分量%x(23)=ReddD转动副D加速度实轴分量%x(24)=ImddD转动副D加速度虚轴分量%y(1)=ddtheta-iBE杆角加速度%y(2)=ddtheta-jCF杆角加速度%y(3)=ddtheta-kDG杆角加速度%y(4)=ddtheta-eFG杆角加速度%y(5)=ddtheta-fGE杆角加速度%y(6)=ddtheta-fEF杆角加速度%y(7)=ReddE转动副E加速度实轴分量%y(8)=ImddE转动副E加速度虚轴分量%y(9)=ReddF转动副F加速度实轴分量%y(10)=ImddF转动副F加速度虚轴分量%y(11)=ReddG转动副G加速度实轴分量%y(12)=ImddG转动副G加速度虚轴分量a=-x(1)*cos(x(7)+pi/2) x(2)*cos(x(8)+pi/2) 0 0 0 -x(6)*cos(x(12)+pi/2); -x(1)*sin(x(7)+pi/2) x(2)*sin(x(8)+pi/2) 0 0 0 -x(6)*sin(x(12)+pi/2); 0 -x(2)*cos(x(8)+pi/2) x(3)*cos(x(9)+pi/2) -x(4)*cos(x(10)+pi/2) 0 0; 0 -x(2)*sin(x(8)+pi/2) x(3)*sin(x(9)+pi/2) -x(4)*sin(x(10)+pi/2) 0 0; x(1)*cos(x(7)+pi/2) 0 -x(3)*cos(x(9)+pi/2) 0 -x(5)*cos(x(11)+pi/2) 0; x(1)*sin(x(7)+pi/2) 0 -x(3)*sin(x(9)+pi/2) 0 -x(5)*sin(x(11)+pi/2) 0;c=-x(1)*cos(x(7)+pi) x(2)*cos(x(8)+pi) 0 0 0 -x(6)*cos(x(12)+pi); -x(2)*sin(x(7)+pi) x(2)*sin(x(8)+pi) 0 0 0 -x(6)*sin(x(12)+pi); 0 -x(2)*cos(x(8)+pi) x(3)*cos(x(9)+pi) -x(4)*cos(x(10)+pi) 0 0; 0 -x(2)*sin(x(8)+pi) x(3)*sin(x(9)+pi) -x(4)*sin(x(10)+pi) 0 0; x(1)*cos(x(7)+pi) 0 -x(3)*cos(x(9)+pi) 0 -x(5)*cos(x(11)+pi) 0; x(1)*sin(x(7)+pi) 0 -x(3)*sin(x(9)+pi) 0 -x(5)*sin(x(11)+pi) 0;b1=c*x(13)2;x(14)2;x(15)2;x(16)2;x(17)2;x(18)2;b2=x(19)-x(21);x(20)-x(22);x(21)-x(23);x(22)-x(24);x(23)-x(19);x(24)-x(20);b=b1+b2;ddtheta=inv(a)*b;y(1)=ddtheta(1);y(2)=ddtheta(2);y(3)=ddtheta(3);y(4)=ddtheta(4);y(5)=ddtheta(5);y(6)=ddtheta(6);y(7)=x(19)+x(1)*ddtheta(1)*cos(x(7)+pi/2)+x(1)*x(13)2*cos(x(7)+pi);y(8)=x(20)+x(1)*ddtheta(1)*sin(x(7)+pi/2)+x(1)*x(13)2*sin(x(7)+pi);y(9)=x(21)+x(2)*ddtheta(2)*cos(x(8)+pi/2)+x(2)*x(14)2*cos(x(8)+pi);y(10)=x(22)+x(2)*ddtheta(2)*sin(x(8)+pi/2)+x(2)*x(14)2*sin(x(8)+pi);y(11)=x(23)+x(3)*ddtheta(3)*cos(x(9)+pi/2)+x(3)*x(15)2*cos(x(9)+pi);y(12)=x(24)+x(3)*ddtheta(3)*sin(x(9)+pi/2)+x(3)*x(15)2*sin(x(9)+pi);这个模块用于求级杆组中各杆的加速度的水平及垂直分量。输入参数为构件2、构件3的角位移和角速度，构件2、构件3和构件4的杆长，构件5的3个边长，构件2、构件3、构件4的角位移和角速度，构件5的3个边向量的角位移和3个转动副B,C,D的加速度；输出参数为构件2、构件3和构件4的角加速度，构件5的3个边向量的角加速度和转动副E,F,G的加速度。2.3 7R六杆级机构MATLAB仿真积分模块初值的确定2.3.1 运用牛顿辛普森法进行角位移分析 图3 7R六杆级机构图3所示是由原动件（曲柄1）和一个6R级杆组所组成的7R六杆级机构，复数向量坐标亦如图所示，各构件的尺寸为r1=120mm，r2=400mm，r3=300mm，r4=300mm，ReD=250mm，FG=450mm，GE=180mm，EF=350mm，ImD=350mm，ReC=700mm，ImC=350mm，构件1以等角速度10rad/s逆时针方向回转，试求构件2和构件3的位移、速度和加速度。由图3可列以下三个方程：AB+BEAD-DG EG=0，即r1+r2-AD-r4-f=0 (2.10)AB+BE+EF-AC-CF=0即r1+r2+g-AC-r3=0 (2.11)EF+FG+GE=0,即g+e+f=0 (2.12)由复向量坐标，可写出式(2.10)、式(2.11)及式(2.12)的角位移方程为： （2.13） (2.14) (2.15)将式(2.13)、式(2.14)、式(2.15)展开，整理得: (2.16)由式(2.16)求出雅可比矩阵为： (2.17)根据式(2.16)、式(2.17)，由牛顿辛普森求解方法得编制M函数如下：function y=r6posi(x)%x(1)=theta-1杆1与水平方向夹角%x(2)=theta-2 杆2与水平方向夹角（估计量）%x(3)=theta-3 杆3与水平方向夹角（估计量）%x(4)=theta-4 杆4与水平方向夹角（估计量）%x(5)=theta-e 杆e与水平方向夹角（估计量）%x(6)=theta-f 杆f与水平方向夹角（估计量）%x(7)=theta-g 杆g与水平方向夹角（估计量）%x(8)=theta-ACAC与水平方向夹角%x(9)=theta-ADAD与水平方向夹角%x(10)=r1杆1长度%x(11)=r2杆2长度%x(12)=r3杆3长度%x(13)=r4杆4长度%x(14)=e杆e长度%x(15)=f杆f长度%x(16)=g杆g长度%x(17)=AC杆AC长度%x(18)=AD杆AD长度%y(1)=theta-22杆与水平方向夹角%y(2)=theta-33杆与水平方向夹角%y(3)=theta-44杆与水平方向夹角%y(4)=theta-ee杆与水平方向夹角%y(5)=theta-ff杆与水平方向夹角%y(6)=theta-gg杆与水平方向夹角%theta2=x(2);theta3=x(3);theta4=x(4);theta5=x(5);theta6=x(6);theta7=x(7);%epsilon=1.0E-6;%f=x(10)*cos(x(1)+x(11)*cos(theta2)-x(18)*cos(x(9)-x(13)*cos(theta4)-x(15)*cos(theta6);x(10)*sin(x(1)+x(11)*sin(theta2)-x(18)*sin(x(9)-x(13)*sin(theta4)-x(15)*sin(theta6);x(10)*cos(x(1)+x(11)*cos(theta2)+x(16)*cos(theta7)-x(17)*cos(x(8)-x(12)*cos(theta3);x(10)*sin(x(1)+x(11)*sin(theta2)+x(16)*sin(theta7)-x(17)*sin(x(8)-x(12)*sin(theta3);x(14)*cos(theta5)+x(15)*cos(theta6)+x(16)*cos(theta7);x(14)*sin(theta5)+x(15)*sin(theta6)+x(16)*sin(theta7);%while norm(f) epsilon J= -x(11)*sin(theta2) 0 x(13)*sin(theta4) 0 x(15)*sin(theta6) 0; x(11)*cos(theta2) 0 -x(13)*cos(theta4) 0 -x(15)*cos(theta6) 0; -x(11)*sin(theta2) x(12)*sin(theta3) 0 0 0 -x(16)*sin(theta7); x(11)*cos(theta2) -x(12)*cos(theta3) 0 0 0 x(16)*cos(theta7); 0 0 0 -x(14)*sin(theta5) -x(15)*sin(theta6) -x(16)*sin(theta7); 0 0 0 x(14)*cos(theta5) x(15)*cos(theta6) x(16)*cos(theta7); dth=inv(J)*(-1.0*f); theta2=theta2+dth(1); theta3=theta3+dth(2); theta4=theta4+dth(3); theta5=theta5+dth(4); theta6=theta6+dth(5); theta7=theta7+dth(6);f=x(10)*cos(x(1)+x(11)*cos(theta2)-x(18)*cos(x(9)-x(13)*cos(theta4)-x(15)*cos(theta6);x(10)*sin(x(1)+x(11)*sin(theta2)-x(18)*sin(x(9)-x(13)*sin(theta4)-x(15)*sin(theta6);x(10)*cos(x(1)+x(11)*cos(theta2)+x(16)*cos(theta7)-x(17)*cos(x(8)-x(12)*cos(theta3);x(10)*sin(x(1)+x(11)*sin(theta2)+x(16)*sin(theta7)-x(17)*sin(x(8)-x(12)*sin(theta3);x(14)*cos(theta5)+x(15)*cos(theta6)+x(16)*cos(theta7);x(14)*sin(theta5)+x(15)*sin(theta6)+x(16)*sin(theta7);norm(f)end;y(1)=theta2;y(2)=theta3;y(3)=theta4;y(4)=theta5;y(5)=theta6;y(6)=theta7;此函数模块的输入量为杆1与水平方向的夹角、其余各杆与水平方向夹角的估计值以及杆的杆长及AC、AD两个收入参量的杆长和角度的估计值，输出参量为2，3，4，E,F,G杆与水平方向的夹角。7R六杆级机构在图1.3.1所示位置，估计构件2，3，4，e，f，g的角位移为2=5.9341rad，3=1.9199rad，4=1.9199rad，e=0，f=2.3562rad，g=0.3491。输入参数x=40*pi/180 340*pi/180 110*pi/180 110*pi/180 0 135*pi/180 20*pi/180 atan(350/700) atan(350/250) 120 400 300 300 450 180 350 sqrt(3502+7002) sqrt(3502+2502)，解得2=-0.3725rad，3=-1.2735rad，4=-1.2735rad，e=3.1416rad，f=-0.8040rad，g=0.3794rad。2.3.2 运用牛顿辛普森法进行角速度分析对式(2.10)-式(2.12)求导并展开成矩阵形式为： (2.18)根据式(2.18)编写M函数如下：function y=r6vel(x)%x(1)=theta-1杆1与水平方向夹角%x(2)=theta-2杆2与水平方向夹角%x(3)=theta-3杆3与水平方向夹角%x(4)=theta-4杆4与水平方向夹角%x(5)=theta-e杆e与水平方向夹角%x(6)=theta-f杆f与水平方向夹角%x(7)=theta-g杆g与水平方向夹角%x(8)=dtheta-1杆1角速度%x(9)=r1杆1长度%x(10)=r2杆2长度%x(11)=r3杆3长度%x(12)=r4杆4长度%x(13)=e杆e长度%x(14)=f杆f长度%x(15)=g杆g长度%y(1)=dtheta-2杆2角加速度%y(2)=dtheta-3杆3角加速度%y(3)=dtheta-4杆4角加速度%y(4)=dtheta-e杆e角加速度%y(5)=dtheta-f杆f角加速度%y(6)=dtheta-g杆g角加速度A= -x(10)*sin(x(2) 0 x(12)*sin(x(4) 0 x(14)*sin(x(6) 0; x(10)*cos(x(2) 0 x(12)*cos(x(4) 0 x(14)*sin(x(6) 0; -x(10)*sin(x(2) x(11)*sin(x(3) 0 0 0 -x(15)*sin(x(7); x(10)*cos(x(2) -x(11)*cos(x(3) 0 0 0 x(15)*cos(x(7); 0 0 0 -x(13)*sin(x(5) -x(14)*sin(x(6) -x(15)*sin(x(7); 0 0 0 x(13)*cos(x(5) x(15)*cos(x(6) x(15)*cos(x(7);B=x(10)*sin(x(1);-x(10)*cos(x(1);x(10)*sin(x(1);-x(10)*cos(x(1);0;0*x(8);y=inv(A)*B;此函数模块输入量为各杆与水平方向的角度以及杆1的角速度，输出参数为2，3，4，E,F,G杆的角加速度。图3所示机构，由位移分析计算出的各杆角度和曲柄1的角速度为10rad/s，则输入参数为：x= 40*pi/180 -0.3725 -1.2735 -1.2735 3.1416 -0.8040 0.3794 10 120 400 300 300 450 180 350，代入上面的M文件，求得2，3，4，e，f，g杆的角速度依次分别为 -3.49rad/s，-4.5298rad/s，-4.5298rad/s，0rad/s，0rad/s，0rad/s。2.4 7R六杆级机构运动学仿真模型及结果2.4.1 7R六杆级机构图3所示是由原动件（曲柄1）和1个6R级杆组所组成的7R六杆级机构，复数向量坐标如图1.4.1所示，各构件的尺寸为的r1=120mm，r2=400mm，r3=300mm，r4=300mm，ReD=250mm，FG=450mm，GE=180mm，EF=350mm，ImD=350mm，ReC=700mm，ImC=350mm，构件1以等角速度10rad/s逆时针方向回转，试求构件2和构件3的位移、速度、加速度。 2.4.2 7R六杆级机构MATLAB运动学仿真模型7R六杆级机构MATLAB运动学仿真模型如图2.4.2所示，在图2.4.2中各积分模块的初值是以曲柄1的幅角为0.72rad和角速度等于10rad/s逆时针方向回转时，相应各个构件的位移、速度的瞬时值，2个MATLAB函数模块分别为crank.m和r6ki.m，其中crank.m函数模块的输入参数为曲柄的长度、角位移、角速度和角加速度；输出参数为曲柄端部（转动副B）的加速度的水平分量和垂直分量;r6ki.m函数模块的输入参数为构件2、构件3的角位移和角速度，构件2、构件3和构件4的杆长，构件5的3个边长，构件2、构件3、构件4的角位移和角速度，构件5的3个边向量的角位移和3个转动副B,C,D的加速度；输出参数为构件2、构件3和构件4的角加速度，构件5的3个边向量的角加速度和转动副E,F,G的加速度。每个数据线上标注了相应变量，常量模块放置了各个构件的尺寸，长度分别为m，角度单位为rad。设置仿真时间为1s，仿真结果输出到工作空间变量simout中，输出格式为array,求解器选用ode45,步长选用变步长。2.4.3 7R六杆级机构MATLAB运动学仿真结果由于曲柄转速为10rad/s,因此每转动1周的时间是0.628s,用绘图命令plot(tout,simout(:,1)，plot(tout,simout(:,2)，plot(tout,simout(:,5)，plot(tout,simout(:,6)，plot(tout,simout(:,9)，plot(tout,simout(:,10)，plot(tout,simout(:,3)，plot(tout,simout(:,4)绘制出构件2和构件3的位移、速度、加速度，构件4的位移、向量e的位移，如图2.4.3所示。从该图中可以看出这些参数也都是周期变化的。从图3所给出的各个构件的尺寸可以看出，由构件3，4，5，6构成平行四边形，因此构件5就作平动，复数向量e,f,g的角位移为常量、角速度为零。从图6（g）中看出复数向量e仿真结果确实如此 。同时构件3，4的角位移、角速度、角加速度应该相同，比较图6（b）和（h）曲线也相同，这两组也说明了所推导的6R级杆组的运动学仿真公式及相应的M函数和仿真模型是正确的。 图6（b） 构件3的角位移 (纵坐标表示角位移的大小，单位为rad；横坐标表示时间，单位为s。) 图6（a） 构件2的角位移(纵坐标表示角位移的大小，单位为rad；横坐标表示时间，单位为s。) 图6（c） 构件3的角速度(纵坐标表示角速度的大小，单位为rad/s；横坐标表示时间，单位为s。) 图6（d） 构件2的角速度(纵坐标表示角位移的大小，单位为rad/s；横坐标表示时间，单位为s。) 图6（e） 构件2的角加速度(纵坐标表示角加速度的大小，单位为rad/s2；横坐标表示时间，单位为s。) 图6（f） 构件3的角加速度(纵坐标表示角加速度的大小，单位为rad/s2；横坐标表示时间，单位为s。) 图6（g） 向量e的角位移(纵坐标表示角位移的大小，单位为rad；横坐标表示时间，单位为s。) 图6（h） 构件4的角位移(纵坐标表示角位移的大小，单位为rad；横坐标表示时间，单位为s。)3 7R六杆级机构动力学仿真3.1 曲柄原动件动力学数学模型的建立3.1.1 曲柄原动件动力学数学分析 图7 曲柄的受力模型如图7所示，已知曲柄AB向量的模ri为常数，幅角i为变量，质心到转动副A的距离为rci,质量为mi，绕质心转动惯量为Ji，作用于质心上的外力为Fxi和Fyi、外力矩为Mi，曲柄与机架联接，转动副A的约束反力为Rxa和Rya，驱动力矩为Ml。 由理论力学可得： (3.1) (3.2) (3.3) 由运动学知识可推得： (3.4) (3.5)将式（3.4）、式（3.5）代入式（3.1）、式（3.2），并与式（3.3）合并得: (3.6)3.1.2 曲柄MATLAB动力学仿真模块M函数根据式(3.6)编写曲柄原动件MATLAB的M函数如下：function y=crankdy(x)%x(1)=theta-i曲柄与水平方向夹角%x(2)=dtheta-i曲柄角速度%x(3)=ddtheta-i曲柄角加速度%x(4)=RxB转动副B约束反力水平分量%x(5)=RyB转动副B约束反力垂直分量%y(1)=RxA转动副A约束反力水平分量%y(2)=RyA转动副A约束反力垂直分量%y(3)=M1转动副A的驱动力矩%g=9.8ri=0.4;rci=0.2;mi=1.2;Ji=0.016;Fxi=0;Fyi=0;Mi=0;ReddA=0;ImddA=0;y(1)=mi*ReddA+mi*rci*x(3)*cos(x(1)+pi/2)+mi*rci*x(2)2*cos(x(1)+pi)-Fxi+x(4);y(2)=mi*ImddA+mi*rci*x(3)*sin(x(1)+pi/2)+mi*rci*x(2)2*sin(x(1)+pi)-Fyi+x(5)+mi*g;y(3)=Ji*x(3)-y(1)*rci*sin(x(1)+y(2)*rci*cos(x(1)-x(4)*(ri-rci)*sin(x(1)+x(5)*(ri-rci)*cos(x(1)-Mi;此函数模块转用于求动副A的约束反力和曲柄上作用的驱动力矩。输入参数是曲柄原动件的角位移、角速度和角加速度 以及转动副B的反作用力；输出参数是转动副A的约束反力和曲柄上作用的驱动力矩。3.2 6R级杆组动力学数学模型的建立3.2.1 6R级杆组动力学数学分析由3个外转动副（B,C,D）和3个内转动副（E,F,G）组成的6R级杆组如图3.2.1(a)所示，其中具有2个转动副的3个构件：BE构件、FC构件和GD构件的受力情况如图3.2.1(b)、(c)、(d)所示。除不受驱动力矩外，其受力分析同曲柄的受力分析是一样的，直接给出公式为： (3.7) (3.8) 具有3个转动副EFG构件受力如图8(e)所示，可推导出： (3.9) (3.10) (3.11) 图8（d）GD杆件受力图 图8（c）FC杆件受力图 图8 (b)BE杆件受力图 图8（a）6R级杆组 图8（e）EFG杆件受力图 图8 各杆受力图将式(3.7)与式（3.11）整理并写成矩阵形式为： (3.12) 3.2.2 6R级杆组MATLAB动力学仿真模块M函数根据式(3.12)编写6R级杆组MATLAB的M函数如下：function y=R6dy(x)%x(1)=theta-iBE杆与水平方向夹角%x(2)=theta-jCF杆与水平方向夹角%x(3)=theta-kDG杆与水平方向夹角%x(4)=theta-eFG杆与水平方向夹角%x(5)=theta-fGE杆与水平方向夹角%x(6)=theta-gEF杆与水平方向夹角%x(7)=dtheta-iBE杆角速度%x(8)=dtheta-jCF杆角速度%x(9)=dtheta-kDG杆角速度%x(10)=dtheta-gFG杆角速度%x(11)=ddtheta-iBE杆角加速度%x(12)=ddtheta-jCF杆角加速度%x(13)=ddtheta-kDG杆角加速度%x(14)=ddtheta-gFG杆角加速度%x(15)=ReddB转动副B加速度实轴分量%x(16)=ImddB转动副B加速度虚轴分量%x(17)=ReddE转动副E加速度实轴分量%x(18)=ImddE转动副E加速度虚轴分量%y(1)=RxB转动副B约束反力水平分量%y(2)=RyB转动副B约束反力垂直分量%y(3)=RxC转动副C约束反力水平分量%y(4)=RyC转动副C约束反力垂直分量%y(5)=RxD转动副D约束反力水平分量%y(6)=RyD转动副D约束反力垂直分量%y(7)=RxE转动副E约束反力水平分量%y(8)=RyE转动副E约束反力垂直分量%y(9)=RxF转动副F约束反力水平分量%y(10)=RyF转动副F约束反力垂直分量%y(11)=RxG转动副G约束反力水平分量%y(12)=RyG转动副G约束反力垂直分量g=9.8;ri=0.4;rj=0.3;rk=0.3;rci=0.2;rcj=0.15;rck=0.15;rce=0.274;rcf=0.188;rcg=0.095;thgp=0.5612;thep=0.1931;thfp=0.5214;mi=8.5;mj=7.6;mk=7.6;ml=12;Ji=0.09;Jj=0.075;Jk=0.075;Jl=0.15;ReddC=0;ImddC=0;ReddD=0;ImddD=0;Fxi=0;Fyi=0;Fxj=0;Fyj=0;Fxk=0;Fyk=0;Fxl=0;Fyl=0;Mi=0;Mj=1000;Mk=1000;Ml=0;a=zeros(12);a(1,1)=1;a(1,7)=1;a(2,2)=1;a(2,8)=1;a(3,1)=rci*sin(x(1);a(3,2)=-rci*cos(x(1);a(3,7)=-(ri-rci)*sin(x(1);a(3,8)=(ri-rci)*cos(x(1);a(4,3)=1;a(4,9)=1;a(5,4)=1;a(5,10)=1;a(6,3)=rcj*sin(x(2);x(6,4)=-rcj*cos(x(2);a(6,9)=-(rj-rcj)*sin(x(2);a(6,10)=(rj-rcj)*cos(x(2);a(7,5)=1;a(7,11)=1;a(8,6)=1;a(8,12)=1;a(9,5)=rck*sin(x(3);a(9,6)=-rck*cos(x(3);a(9,11)=-(rk-rck)*sin(x(3);a(9,12)=(rk-rck)*cos(x(3);a(10,7)=-1;a(10,9)=-1;a(10,11)=-1;a(11,8)=-1;a(11,10)=-1;a(11,12)=-1;a(12,7)=-rcg*sin(x(6)+thgp);a(12,8)=rcg*cos(x(6)+thgp);a(12,9)=-rce*sin(x(4)+thep);a(12,10)=rce*cos(x(4)+thep);a(12,11)=-rcf*sin(x(5)+thfp);a(12,12)=rcf*cos(x(5)+thfp);b=zeros(12,1)b(1,1)=mi*x(15)+mi*rci*x(11)*cos(x(1)+pi/2)+mi*rci*x(7)2*cos(x(1)+pi)-Fxi;b(2,1)=mi*x(16)+mi*rci*x(11)*sin(x(1)+pi/2)+mi*rci*x(7)2*sin(x(1)+pi)-Fyi+mi*g;b(3,1)=Ji*x(11)-Mi;b(4,1)=mi*ReddC+mj*rcj*x(12)*cos(x(2)+pi/2)+mj*rcj*x(8)2*cos(x(2)+pi)-Fxj;b(5,1)=mi*ReddC+mj*rcj*x(12)*sin(x(2)+pi/2)+mj*rcj*x(8)2*sin(x(2)+pi)-Fyj+mj*g;b(6,1)=Jj*x(12)-Mj;b(7,1)=mi*ReddD+mk*rck*x(13)*cos(x(3)+pi/2)+mk*rck*x(9)2*cos(x(3)+pi)-Fxk;b(8,1)=mi*ReddD+mk*rck*x(13)*sin(x(3)+pi/2)+mk*rck*x(9)2*sin(x(3)+pi)-Fyk+mj*g;b(9,1)=Jk*x(13)-Mk;b(10,1)=ml*x(17)+ml*rcg*x(14)*cos(x(6)+thgp+pi/2)+ml*rcg*x(10)2*cos(x(6)+thgp+pi)-Fxl;b(11,1)=ml*x(18)+ml*rcg*x(14)*sin(x(6)+thgp+pi/2)+ml*rcg*x(10)2*sin(x(6)+thgp+pi)-Fyl;b(12,1)=Jl*x(14)-Ml;y=inv(a)*b;此函数模块用于求转动副B,C,D,E,F,G的约束反力。输入参数为构件2、构件3和构件4的角位移、角速度和角加速度、构件5的3个边向量的角位移及转动副B和E的加速度；输出参数是转动副B,C,D,E,F,G的约束反力。3.3 需要引用的函数 在机构的动力学仿真中，需要用到运动学中的两个函数，即crank.m及r6ki.m,同时，为了简化动力学simulink仿真模型，需要对这两个M文件进行细微改变，将各杆的杆长作为书籍量写入函数中，同时删去不必要的输出量，以达到简化模型的目的。经过改造后的函数分别命名为crank_1.m及r6ki_1.m。3.3.1 crank_1.m 此文件将曲柄的长度r作为已知量写入函数，得crank_1.m为：function y=crank_1(x)%x(1)=thetajCF杆与水平方向夹角%x(2)=dthetajCF杆角速度%x(3)=ddthetajCF杆角加速度%y(1)=ReddB转动副B加速度实轴分量%y(2)=ImddB转动副B加速度虚轴分量%r=0.12ddB=r*x(3)*cos(x(1)+pi/2)+r*x(2)2*cos(x(1)+pi); r*x(3)*sin(x(1)+pi/2)+r*x(2)2*sin(x(1)+pi);y=ddB;3.3.2 r6ki_1.m 此函数将2，3，4，e，f，g杆的杆长作为已知量写入函数，并删去输出量F和G的加速度，得r6ki_1.m为：function y=R6ki_1(x)%Input parameters%x(1)=theta-iBE杆与水平方向夹角%x(2)=theta-jCF杆与水平方向夹角%x(3)=theta-kDG杆与水平方向夹角%x(4)=theta-eFG杆与水平方向夹角%x(5)=theta-fGE杆与水平方向夹角%x(6)=theta-gEF杆与水平方向夹角%x(7)=dtheta-iBE杆角速度%x(8)=dtheta-jCF杆角速度%x(9)=dtheta-kDG杆角速度%x(10)=dtheta-eFG杆角速度%x(11)=dtheta-fGE杆角速度%x(12)=dtheta-gEF杆角速度%x(13)=ReddB转动副B加速度实轴分量%x(14)=ImddB转动副B加速度虚轴分量%y(1)=ddtheta-iBE杆角加速度%y(2)=ddtheta-jCF杆角加速度%y(3)=ddtheta-kDG杆角加速度%y(4)=ddtheta-eFG杆角加速度%y(5)=ddtheta-fGE杆角加速度%y(6)=ddtheta-gEF杆角加速度%y(7)=ReddE转动副E加速度实轴分量%y(8)=ImddE转动副E加速度虚轴分量%r2=0.4;r3=0.3;r4=0.3;re=0.15;rf=0.18;rg=0.35;a=-r2*cos(x(1)+pi/2) r3*cos(x(2)+pi/2) 0 0 0 -rg*cos(x(6)+pi/2); -r2*sin(x(1)+pi/2) r3*sin(x(2)+pi/2) 0 0 0 -rg*sin(x(6)+pi/2); 0 -r3*cos(x(2)+pi/2) r4*cos(x(3)+pi/2) -re*cos(x(4)+pi/2) 0 0; 0 -r3*sin(x(2)+pi/2) r4*sin(x(3)+pi/2) -re*sin(x(4)+pi/2) 0 0; r2*cos(x(1)+pi/2) 0 -r4*cos(x(3)+pi/2) 0 -rf*cos(x(5)+pi/2) 0; r2*sin(x(1)+pi/2) 0 -r4*sin(x(3)+pi/2) 0 -rf*sin(x(5)+pi/2) 0;c=-r2*cos(x(1)+pi) r3*cos(x(2)+pi) 0 0 0 -rg*cos(x(6)+pi); -r2*sin(x(1)+pi) r3*sin(x(2)+pi) 0 0 0 -rg*sin(x(6)+pi); 0 -r3*cos(x(2)+pi) r4*cos(x(3)+pi) -re*cos(x(4)+pi) 0 0; 0 -r3*sin(x(2)+pi) r4*sin(x(3)+pi) -re*sin(x(4)+pi) 0 0; r2*cos(x(1)+pi) 0 -r4*cos(x(3)+pi) 0 -rf*cos(x(5)+pi) 0; r2*sin(x(1)+pi) 0 -r4*sin(x(3)+pi) 0 -rf*sin(x(5)+pi) 0;b1=c*x(7)2;x(8)2;x(9)2;x(10)2;x(11)2;x(12)2;b2=x(13);x(14);0;0;-x(13);-x(14);b=b1+b2;ddtheta=inv(a)*b;y(1)=ddtheta(1);y(2)=ddtheta(2);y(3)