机械原理大作业--牛头刨床仿真分析(matlab)全自动版本
机械原理大作业——牛头刨床
机械原理大作业——牛头刨床大作业，一，平面连杆机构的运动分析题号: 6班级 : 姓名 : 学号 : 同组者 :成绩 :完成时间 :目录题目、原始数据及要求 ..................................................................... .......................1 一平面连杆机构运动分析方程 ..................................................................... . (1)1.1速度计算公式 ..................................................................... .. (2)1.2加速度计算公式 ..................................................................... ..............2 二程序 ..................................................................... (3)2.1计算程序框图 ..................................................................... (3)2.2计算源程序 ..................................................................... .........................4 三 3.1 (一组数据 Lab =200mm)计算结果 (9)3.2运动线图 ..................................................................... . (10)3.3 体会 ..................................................................... .................................... 12 四 4.1(第二组数据 Lab =150mm)计算结果 . (12)4.2 运动线图 ..................................................................... .. (13)4.3 体会 ..................................................................... .................................... 15 五 5.1(第三组数据 Lab =220mm)计算结果 . (16)5.2 运动线图 ..................................................................... (17)5.3 体会 ..................................................................... ...................................... 21 六参考资料 ..................................................................... (21)题目、原始数据及要求:图所示为一牛头刨床(?级机构)。
牛头刨床MATLAB运动仿真
牛头刨床机构简图二、运动方程的求解(1)对如右图所示的矢量三角形矢量方程为：.321R R R =+(*)将上述矢量方程分别向水平、竖直方向投影，得：⎩⎨⎧=+=.sin sin ;cos cos 332213322θθθθr r r r r （1-1）上式分别对时间t 求导，得：⎪⎩⎪⎨⎧+=-=-∙∙.cos sin cos ;sin cos sin 3333322233333222θωθθωθωθθωr r r r r r （1-2）上式写成关于3ω、∙3r 的矩阵方程，为：.cos sin sin cos cos sin 22222233333333⎪⎪⎭⎫⎝⎛-=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛-∙θωθωωθθθθr r r r r (1-3)将(2)式分别对时间t 求导，得：⎪⎩⎪⎨⎧-++=----=--∙∙∙∙∙∙∙∙.sin cos cos 2sin sin cos ;cos sin sin 2cos cos sin 33233333333322222223323333333332222222θωθαθωθθωθαθωθαθωθθωθαr r r r r r r r r r r r (1-4)上式写成关于3α、∙∙3r的矩阵方程，为：一、牛头刨床机构简介：牛头刨床是一种常用的平面切削加工机床。
电动机经皮带传动、齿轮传动最后带动曲柄1（见右图）转动。
刨床工作时，是由导杆机构1-2-3-4-5带动刨头和刨刀做往复直线运动刨头5右行时，刨头切断，称为工作行程，此时要求速度较低且均匀；刨头左行时，不进行切削，称为空行程，此时速度较高，以节省时间提高生产率。
为此刨床采用有急回作用的导杆机构。
这里给出各构件的长度，以供下面的仿真使用。
.88.210,92.702,72315102.133===DE CD AB l l l.sin cos 2sin cos cos sin 2cos sin sin cos cos sin 332333322222223323333222222233333333⎪⎪⎭⎫ ⎝⎛+--++--=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛-∙∙∙∙∙∙θωθωθωθαθωθωθωθααθθθθr r r r r r r r r r r (1-5)(2)对如右图所示的矢量三角形矢量方程为：.7654R R R R +=+(**)将上述矢量方程分别向水平、竖直方向投影，得：⎩⎨⎧=+=+.sin sin ;cos cos 6554475544r r r r r r θθθθ（2-1）上式分别对时间t 求导，得：⎪⎩⎪⎨⎧=+=--∙.0sin cos ;sin sin 5554447555444θωθωθωθωr r r r r （2-2）将(2)式分别对时间t 求导，得：⎪⎩⎪⎨⎧=-+-=----∙∙.0sin cos sin cos ;cos sin cos sin 55255554424444755255554424444θωθαθωθαθωθαθωθαr r r r r r r r r (2-3)上式写成关于5α、∙∙7r 的矩阵方程，为：.sin sin cos cos sin 0cos 1sin 552544244444424444755555⎪⎪⎭⎫ ⎝⎛--+=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--∙∙θωθωθαθωθααθθr r r r r r r r (2-4)三、构建SIMULINK 方框图a)构建速度方框图并仿真1.首先对方程(1-3)的常变量分析：已知常量有：2ω,2r ；变量有：2θ,3θ,3ω,3r ,3∙r 。
机械原理课程设计——牛头刨床.
一：课程设计题目、内容及其目的课题：牛头刨床内容1．对机构进行运动分析已知：曲柄每分钟转数错误!未找到引用源。
,各构件尺寸及质心位置。
作机构1~2个位置的速度多边形和加速度多边形，作滑块的运动线图，以上内容与后面动态静力分析一起画在1号图纸上。
2．对机构进行动态静力分析已知：各构件的重量G（曲柄1、滑块2、和连杆5的重量都可以忽略不计），导杆3的转动惯量错误!未找到引用源。
及切削力错误!未找到引用源。
变化规律如下图。
确定构件一个位置的各运动副反力及应加于曲柄上的平衡力矩。
3、用UG进行模拟运动仿真校核机构运动分析和动态静力分析的结果4、电动机功率的确定与型号的选择5、齿轮减速机构设计目的：1：学会机械运动见图设计的步骤和方法；2：巩固所学的理论知识，掌握机构分析与综合的基本方法；3：培养学生使用技术资料，计算作图及分析与综合能力；4：培养学生进行机械创新的能力。
二：牛头刨床简介和机构的要求1：牛头刨床简介牛头刨床是一种用于平面切削加工的机床，如图1。
电动机经皮带和齿轮传动，经过减速机构减速从而带动曲柄1。
刨床工作时，由导杆3 经过连杆4 带动刨刀5 作往复运动。
刨头左行时，刨刀进行切削，称工作行程，此时要求速度较低并且均匀，以减少电动机容量和提高切削质量，刨头右行时，刨刀不切削，称空行程，此时要求速度较高，以提高生产率。
为此刨床采用有急回作用的导杆机构。
刨刀每切削完一次，利用空回行程的时间，通过棘轮带动螺旋机构（图中未画），使工作台连同工件作一次进给运动，以便刨刀继续切削。
刨头在工作行程中，受到很大的切削阻力（在切削的前后各有一段约0.05H 的空刀距离），而空回行程中只有摩擦阻力。
因此刨头在整个运动循环中，受力变化是很大的，这就影响了主轴的匀速运转，故需安装飞轮来减小主轴的速度波动，以提高切削质量和减小电动机容量。
2：机构的要求牛头刨床的主传动的从动机构是刨头，在设计主传动机构时，要满足所设计的机构要能使牛头刨床正常的运转，同时设计的主传动机构的行程要有急回运动的特性，刨削速度尽可能为匀速运动，以及很好的动力特性。
机械原理课程设计——牛头刨床
机械能变化曲线：
飞轮设计：
Ｖ
A4
＝
A2 A4 A2
速度图解法：
V1A+V12=V 2A VF+VFB=V 2B V2B=βV 2A Β为常数比
加速度图解分析： a4An+a4Ar+a24Ar+ak24A =a2A 大小 方向
a4b+aF4Br=aF a4A=βV 4B
进给凸轮机构设计
主体机构设计
牛头刨床主体机构
主体结构设计
设计要求
（1）刨刀工作行程要求速度比较平稳，空回行程时 刨刀快速退回，机构行程速比系数在1.4左右。 （2）刨刀行程H=300mm或H=150mm。曲柄转速、 切削力、许用传动角等见表1，每人选取其中一组数据。 （3）切削力P大小及变化规律如图1所示，在切削行 程的两端留出一点空程。具体数据如下：
主体机构
电机转速n(r/mi n)
切削力P(N)
75
许用传动角[γ]
H=150mm
4500N
45°
刨刀行程：H=150 速比系数：K=1.4
主体机构（方案一）
方案一： 摆动导杆机构与摇杆滑块机构组合机构
机构简图：
计算机构的自由度 F=3×5－2×7=1
主体机构（方案一）
机构尺寸的计算：
在满足压力角条件确定基圆半径，摆杆中心间的中心距。
• 推程许用压力角为[α]= 38°； • 回程许用压力角为[α’]= 65°； • 试凑法：对照摆杆长度为L，赋值基圆半径， 中心距a=90,r0=50;经试验符合要求
滚子半径rf：rf<ρ mi n －3（mm）及rf<0.8ρ mi n（mm） 方法1用图解法确定凸轮理论廓线上某点A的曲率半径R： 以A点位圆心,任选较小的半径r 作圆交于廓线上，在圆A 两边分别以理论廓线上的B、C为圆心，以同样的半径r 画圆，三个小圆分别交于E、F、H、M四个点处。过E、 F H、M O点 O点近似为凸轮廓线上A OA。并且曲率中心肯定在曲线过A 点的法线上。可以通 过法线与直线EF或HM的交点求曲率中心。
机械原理课程设计牛头刨床讲解
机械原理课程设计——牛头刨床(速度分析与受力分析附于最后)说明书姓名：分析点：4,10点组号：第3组2011 年 7 月 15日工作原理 (3)一.设计任务 (4)二.设计数据 (4)三.设计要求 (4)1、运动方案设计 (4)2、确定执行机构的运动尺寸 (4)3、进行导杆机构的运动分析 (5)4、对导杆机构进行动态静力分析 (5)四.设计方案选定 (5)五. 机构的运动分析 (6)1. 4点速度分析,加速度分析 (7)2. 10点速度,加速度分析 (9)六.机构动态静力分析 (11)七.数据总汇并绘图 (13)九.参考文献 (16)工作原理牛头刨床是一种用于平面切削加工的机床，如图a）所示。
电动机经过皮带和齿轮传动，带动曲柄2和固结在其上的凸轮8。
刨床工作时，由导杆机构2－3－4－5－6带动刨头6和刨刀7作往复运动。
刨头左行时，刨刀不切削，称为空回行程，此时要求速度较高，以提高生产率。
为此刨床采用有急回运动的导杆机构。
刨刀每切削完一次，利用空回行程的时间，凸轮8通过四杆机构1－9－10－11与棘轮带动螺旋机构（图中未画），使工作台连同工件作一次进给运动，以便刨刀继续切削。
刨头在工作过程中，受到很大的切削阻力（在切削的前后各有一段0.05H的空刀距离，见图b），而空回行程中则没有切削阻力。
因此刨头在整个运动循环中，受力变化是很大的，这就影响了主轴的匀速转动，故需安装飞轮来减小主轴的速度波动，以提高切削质量和减少电动机容量。
(a) (b)图d一.设计任务1、运动方案设计。
2、确定执行机构的运动尺寸。
3、进行导杆机构的运动分析。
4、对导杆机构进行动态静力分析。
5、汇总数据画出刨头的位移、速度、加速度线图以及平衡力矩的变化曲线。
二.设计数据本组选择第六组数据表1表2三.设计要求1、运动方案设计根据牛头刨床的工作原理，拟定1～2个其他形式的执行机构（连杆机构），给出机构简图并简单介绍其传动特点。
2、确定执行机构的运动尺寸根据表一对应组的数据，用图解法设计连杆机构的尺寸，并将设计结果和步骤写在设计说明书中。
第2-1节 牛头刨床运动仿真
第2讲牛头刨床运动仿真实例一、启动pro/e并设置工作目录1.点击【开始】→【所有程序】→【pro engineer】，启动pro/e软件,如图1所示。
图1 启动pro/e2.设置工作目录：选择【文件】→【设置工作目录】，选择桌面《运动仿真牛头刨床》文件夹为工作目录，点击该图框右下方的【确定】键，如图2、图3所示。
图2 设置工作目录图3 选择工作目录二、设置运动仿真文件1.点击【文件】→【新建】命令，选择类型为【组件】，子类型选择【设计】，将名称改为【niutoubaochuang】，将【使用缺省模板】前的对勾去掉，,如图4所示。
单击【确定】后弹出【新文件选项】对话框，【模板】选mmns asm design，单击【确定】完成任务的新建，如图5所示。
图4 新建文件图5 选择模板三、零件的装配1.预览整个机构效果图单击工作窗口右边工具栏中的【装配】命令，在【打开】命令中单击【niutoubaoch.asm】项。
单击【预览】可在框图中查看整个凸轮机构组装好后的效果图，如图6所示，后面的组装皆可按此标准来进行。
图6 机构效果图2. 安装基座（1）单击工作窗口右边的【装配】命令，在【打开】命令中选择ground.prt文件，即机座（也可以现【预览】确认一下），点击【打开】如图7所示。
图7 选择基座文件（2）导入机座后单击图8中所示【自动】命令右边小三角形，将机座设置为【缺省】模式，再点击面板右端的对勾，完成机座的放置，如图8所示。
图8 设置基座3. 安装圆盘（1）点击【装配】选中并打开打开第二个文件yuanpan.prt，即圆盘机构然后预览并单击【打开】,如图9所示，点击工具栏旁的【用户定义】下拉菜单，选择【销钉】连接类型,如图10所示。
图9 选择圆盘文件图10 设置约束类型（2）先进行轴对齐，即选择机座机架上矮圆柱的中心线A-2，再选择圆盘上圆孔的中心线A-1，完成轴对齐，如图11所示。
图11 轴对齐（3）再面匹配，光标移到圆盘上选中表面，然后选择矮圆柱的表面，点击工具栏上的对勾完成凸轮的装配，如图12所示。
机械原理大作业牛头刨床 matlab编程
l4*cosd(x(2))+l3*cosd(x(3))-x(4);
l4*sind(x(2))+l3*sind(x(3))-h];
主程序：
theta1=0:360;
s3=0:360;
theta3=0:360;
-omega3(i)*cosd(theta3(i)) l4*omega4(i)*sind(theta4(i)) -v3(i)*cosd(theta3(i))+s3(i)*omega3(i)*sind(theta3(i)) 0;
0 l4*omega4(i)*cosd(theta4(i)) l3*omega3(i)*cosd(theta3(i)) 0;
xlabel('{\theta}_1(^{ o })');
ylabel('s_c (mm)');
figure(5)
plot(theta1,vc,'k');
gridon;
xlabel('{\theta}_1(^{ o })');
ylabel('v_c (m/s)');
figure(6)
plot(theta1,ac,'k');
gridon;
xlabel('{\theta}_1(^{ o })');
ylabel('{\alpha}_3 {\alpha}_4 (r a d / s^2)');
legend('{\alpha}_3','{\alpha}_4');
(完整word版)牛头刨床运动仿真matlab程序
附录牛头刨床主运动机构MATLAB程序由主程序six_bar_main 和子函数six_bar两部分组成。
1。
主程序six_bar_main文件％1.输入已知数据clear；l1=0.125；l3=0。
600;l4=0。
150;l6=0.275;l61=0.575;omega1=1；alpha1=0；hd=pi/180；du=180/pi；％2。
调用子函数six_bar计算牛头刨床机构位移，角速度，角加速度for n1=1：459；theta1(n1)=—2*pi+5。
8199+(n1—1)＊hd;ll=［l1，l3，l4,l6，l61]；[theta，omega,alpha］=six_bar（theta1(n1),omega1,alpha1，ll);s3（n1)=theta（1）;theta3（n1)=theta（2)；theta4(n1）=theta(3）;sE（n1）=theta（4）；v2(n1)=omega(1);omega3（n1）=omega(2)；omega4（n1)=omega(3)；vE(n1）=omega（4)；a2(n1)=alpha（1）;alpha3（n1)=alpha（2);alpha4(n1）=alpha（3);aE(n1）=alpha(4）;end%3.位移、角速度、角加速度、和牛头刨床图形输出figure(3);n1=1：459；t=(n1-1）＊2*pi/360；subplot（2，2，1)；％绘角位移及位移线图plot(t,theta3＊du，’r-。
’）；grid on；hold on；axis auto；［haxes，hline1,hline2］=plotyy（t，theta4＊du，t,sE)；grid on;hold on;xlabel（’时间/s')axes（haxes(1));ylabel（'角位移/\circ’）axes(haxes(2)）；ylabel(’位移/m’）hold on;grid on;text(1。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
牛头刨床机构运动分析一、设计小组人员构成二、设计任务1、机构结构分析2、机构运动分析建立数学模型，解析法进行运动分析；程序编写；上机调试程序；位移、速度和加速度运动曲线图与分析；三、设计参数四、设计分工五、 设计内容如右图，建立直角坐标系，并标出各杆矢量及方位角。
利用两个封闭图形ABCA 及CDEGC ，由此可得:l 6⃑⃑⃑ +l 3⃑⃑⃑ =s 3⃑⃑⃑ , l 3⃑⃑⃑ +l 4⃑⃑⃑ =l 6′⃑⃑⃑ +S E ⃑⃑⃑⃑投影方程式为:s 3 cos θ3=l 1 cos θ1 s 3sin θ3=l 6+l 1sin θ1 l 3cos θ3+l 4cos θ4−s E =0l 3sin θ3+l 4sin θ4=l 6′有以上各式即可求得s 3 、θ3 、θ4 四个运动那个变量，二滑块2的方位角θ2=θ3。
然后，分别将上式对时间去一次、二次导数，并写成矩阵形式，即得以下速度和加速度方程式：[cos θ3−s 3sin θ3sin θ3s 3cos θ30 00 00 −l 3sin θ30l 3cos θ3−l 4sin θ4−1l 4cos θ40] [s 3ω3ω4v E ]=ω1[−l 1sin θ1l 1cos θ100][cos θ3−s 3sin θ3sin θ3s 3cos θ30 00 00 −l 3sin θ30l 3cos θ3−l 4sin θ4−1l 4cos θ40] [s 3α3α4αE ] =- [−ω3sin θ3 −s 3sin θ3−s 3ω3sin θ3ω3cos θ3 s 3cos θ3−s 3ω3sin θ30 00 00 −l 3ω3cos θ30 −l 3ω3sin θ3−l 4ω4cos θ4 0−l 4ω4sin θ4 0] [s 3ω3ω4v E ]+ ω1[−l 1ω1cos θ1−l 1ω1sin θ100]而ω2=ω3 、α2=α3 。
根据以上各式，将已知参数带入，编写程序，利用计算机辅助计算，求得的数值列于下表。
并根据数值做出机构运动的位置线图、速度线图、家速度线图。
六、结论A、计算数据θ1θ3θ4s Eω3ω4v Eα3α4αE0 65.5560 168.9382 0.1011 0.1712 0.2888 -0.1018 0.2477 0.2927 -0.138710.0000 67.4668 172.0273 0.0814 0.2093 0.3239 -0.1227 0.1908 0.1172 -0.102620.0000 69.7125 175.3266 0.0585 0.2386 0.3320 -0.1383 0.1472 -0.0185 -0.078030.0000 72.2163 178.5980 0.0333 0.2612 0.3192 -0.1504 0.1132 -0.1241 -0.061440.0000 74.9188 181.6560 0.0062 0.2785 0.2900 -0.1601 0.0861 -0.2076 -0.050350.0000 77.7722 184.3541 -0.0225 0.2915 0.2477 -0.1681 0.0640 -0.2742 -0.042360.0000 80.7379 186.5757 -0.0524 0.3010 0.1951 -0.1749 0.0454 -0.3264 -0.035570.0000 83.7830 188.2297 -0.0835 0.3075 0.1346 -0.1805 0.0290 -0.3646 -0.028680.0000 86.8791 189.2495 -0.1154 0.3113 0.0687 -0.1848 0.0142 -0.3882 -0.020390.0000 90.0000 189.5941 -0.1479 0.3125 0.0000 -0.1875 0 -0.3962 -0.0099100.0000 93.1209 189.2495 -0.1807 0.3113 -0.0687 -0.1881 -0.0142 -0.3882 0.0030 110.0000 96.2170 188.2297 -0.2134 0.3075 -0.1346 -0.1863 -0.0290 -0.3646 0.0183 120.0000 99.2621 186.5757 -0.2456 0.3010 -0.1951 -0.1816 -0.0454 -0.3264 0.0357 130.0000 102.2278 184.3541 -0.2766 0.2915 -0.2477 -0.1738 -0.0640 -0.2742 0.0544 140.0000 105.0812 181.6560 -0.3061 0.2785 -0.2900 -0.1626 -0.0861 -0.2076 0.0737 150.0000 107.7837 178.5980 -0.3332 0.2612 -0.3192 -0.1481 -0.1132 -0.1241 0.0929 160.0000 110.2875 175.3266 -0.3575 0.2386 -0.3320 -0.1302 -0.1472 -0.0185 0.1114 170.0000 112.5332 172.0273 -0.3785 0.2093 -0.3239 -0.1092 -0.1908 0.1172 0.1289 180.0000 114.4440 168.9382 -0.3955 0.1712 -0.2888 -0.0852 -0.2477 0.2927 0.1464 190.0000 115.9198 166.3668 -0.4080 0.1217 -0.2190 -0.0580 -0.3229 0.5144 0.1669 200.0000 116.8285 164.6992 -0.4155 0.0571 -0.1069 -0.0263 -0.4224 0.7754 0.1980 210.0000 116.9955 164.3854 -0.4168 -0.0275 0.0518 0.0126 -0.5522 1.0375 0.2539 220.0000 116.1941 165.8703 -0.4103 -0.1375 0.2504 0.0649 -0.7135 1.2133 0.3538 230.0000 114.1449 169.4393 -0.3929 -0.2775 0.4619 0.1392 -0.8905 1.1564 0.5061 240.0000 110.5470 174.9633 -0.3600 -0.4461 0.6286 0.2423 -1.0255 0.6624 0.6684 250.0000 105.1832 181.5495 -0.3071 -0.6259 0.6559 0.3651 -0.9935 -0.4598 0.6995 260.0000 98.1324 187.2640 -0.2337 -0.7742 0.4416 0.4682 -0.6461 -1.9991 0.4257 270.0000 90.0000 189.5941 -0.1479 -0.8333 -0.0000 0.5000 -0.0000 -2.8172 -0.0704 280.0000 81.8676 187.2640 -0.0639 -0.7742 -0.4416 0.4515 0.6461 -1.9991 -0.4435 290.0000 74.8168 181.5495 0.0072 -0.6259 -0.6559 0.3598 0.9935 -0.4598 -0.5742 300.0000 69.4530 174.9633 0.0612 -0.4461 -0.6286 0.2589 1.0255 0.6624 -0.5677 310.0000 65.8551 169.4393 0.0980 -0.2775 -0.4619 0.1646 0.8905 1.1564 -0.5068 320.0000 63.8059 165.8703 0.1194 -0.1375 -0.2504 0.0832 0.7135 1.2133 -0.4244 330.0000 63.0045 164.3854 0.1279 -0.0275 -0.0518 0.0168 0.5522 1.0375 -0.3369 340.0000 63.1715 164.6992 0.1261 0.0571 0.1069 -0.0348 0.4224 0.7754 -0.2561 350.0000 64.0802 166.3668 0.1165 0.1217 0.2190 -0.0734 0.3229 0.5144 -0.1893 360.0000 65.5560 168.9382 0.1011 0.1712 0.2888 -0.1018 0.2477 0.2927 -0.1387B、运动线图C、Matlab程序源代码>> w1=1;l1=0.125;l3=0.6;l6=0.275;l61=0.575;l4=0.15;for m=1:3601o1(m)=pi*(m-1)/1800;o31(m)=atan((l6+l1*sin(o1(m)))/(l1*cos(o1(m))));if o31(m)>=0o3(m)=o31(m);else o3(m)=pi+o31(m);end;s3(m)=(l1*cos(o1(m)))/cos(o3(m));o4(m)=pi-asin((l61-l3*sin(o3(m)))/l4);se(m)=l3*cos(o3(m))+l4*cos(o4(m));if o1(m)==pi/2o3(m)=pi/2; s3(m)=l1+l6;endif o1(m)==3*pi/2o3(m)=pi/2; s3(m)=l6-l1;endA1=[cos(o3(m)),-s3(m)*sin(o3(m)),0,0;sin(o3(m)),s3(m)*cos(o3(m)),0,0;0,-l3*sin(o3(m)),-l 4*sin(o4(m)),-1;0,l3*cos(o3(m)),l4*cos(o4(m)),0];B1=w1*[-l1*sin(o1(m));l1*cos(o1(m));0;0];D1=A1\B1;E1(:,m)=D1;ds(m)=D1(1);w3(m)=D1(2);w4( m)=D1(3);ve(m)=D1(4);A2=[cos(o3(m)),-s3(m)*sin(o3(m)),0,0;sin(o3(m)),s3(m)*cos(o3(m)),0,0;0,-l3*sin(o3(m)),-l 4*sin(o4(m)),-1;0,l3*cos(o3(m)),l4*cos(o4(m)),0];B2=-[-w3(m)*sin(o3(m)),(-ds(m)*sin(o3(m))-s3(m)*w3(m)*cos(o3(m))),0,0;w3(m)*cos(o3(m)),( ds(m)*cos(o3(m))-s3(m)*w3(m)*sin(o3(m))),0,0;0,-l3*w3(m)*cos(o3(m)),-l4*w4(m)*cos(o4(m)) ,0;0,-l3*w3(m)*sin(o3(m)),-l4*w4(m)*sin(o4(m)),0]*[ds(m);w3(m);w4(m);ve(m)];C2=w1*[-l1*w1*cos(o1(m));-l1*w1*sin(o1(m));0;0];B=B2+C2;D2=A2\B;E2(:,m)=D2;dds(m)=D2(1); a3(m)=D2(2);a4(m)=D2(3);ae(m)=D2(4);end;o11=o1*180/pi;y=[o3*180/pi;o4*180/pi];w=[w3;w4];a=[a3;a4];figure;subplot(221);h1=plotyy(o11,y,o11, se); axis equal;title('位置线图');xlabel('\it\theta1');ylabel('\it\theta3,\theta4,Se');text(315,215,'SE');text(180,180,'θ4');text(70,105,'θ3');subplot(222);h2=plotyy(o11,w,o11,ve);title('速度线图');text(300,0.62,'VE');text(165,0.35,'ω3');text(90,0.1,'ω4');xlabel('\it\theta1');ylabel('\it\omega3,\omega4,Ve');subplot(212);h3=plotyy(o11,a,o11,ae);title('加速度线图');text(185,-0.5,'α3');text(100,-0.3,'α4');text(100,0.7,'αE');。