matlab对牛头刨床的仿真(速度,加速度,力,运动动画,matlab可运行程序)

北华航天工业学院牛头刨床主体机构的分析《机械原理》课程设计班级_______________姓名_______________序号_______________2015年1月摘要：文中建立了用解析法作牛头刨床的运动分析，基于机构的位置方程式，将位置方程对时间求一次和二次倒数，得到机构的速度和加速的方程，然后就得到位移、速度、加速度，并通过复数矢量法和矩阵法，并且借助Matlab等软件对牛头刨床进行运动分析，并对牛头刨床的运动学规律进行了仿真，画出牛头刨床的运动动画,使设计的机构更加直白,为该机床的优化设计提供了理论支撑。

目录摘要 (1)目录 (2)1.引言 (3)2.机构结构分析 (3)3.机构运动分及程序设计 (3)4.计算机仿真及结果 (4)5.图解法验证 (17)6.guide的制作 (22)7.参考文献 (30)8.研究感想 (30)1.引言：对牛头刨床进行运动分析的任务是在已知机构尺寸及原动件运动规律的情况下，确定机构中其他构件上某些点的轨迹、位移、速度、加速度和构件的角位移、角速度及角加速度。

而这些内容，不论是设计新的机器，还是为了了解现有机械的运动性能，都是十分必要的，而且它还是研究机械动力性能的必要基础和前提。

对于牛头刨床来说，其刨刀在有效工作行程中应该接近等速运动，而回程速度却应高于前者速度，从而提高效率。

为了确定牛头刨床的设计是否满足要求，就必须对其进行运动分析。

本文建立了牛头刨床的力学模型，并推导了其运动微分方程，通过数值方法求得了牛头刨床的运动学规律并开发了一个用于牛头刨床的通用运动学分析程序。

2.机构结构分析：1.自由度计算：F=3n-2pl -phF=5*3-7*2=12.结构分析：5个可动构件,5个转动副2个移动副.3.机构运动分及程序设计3.1.如图先建立坐标系，已知l1=100mm, l3=500mm，l4=200mm,原动件1的方位角00～3600和等角速度w 1=2*pi rad/s 如图所示：解：3.1.1.位置分析：由于有四个位置量，为了求解需要建立两个封闭矢量方程。

有图形ABCA 可写出矢量方程:sl l 316=+可变为：e θθe e313i12πi6is l l =+ (1)θθθθ331163316s i n s i n c o s c o s s l l s l =+=sl ll l s3113221163θθ)θc o s(θc o s a r c c o s 11)s i n (=++=由封闭图形CDEGC 可写出矢量方程：3346's l ll +=+34πiii '2346θθe e eEll l s +=+ (2)得到方程：ll l l l '64433E4433θθs θθs i n s i n 0-c o s c o s =+=+由上式可得：'633443344s i n a r c s i n c o s c o s θθθθEl l ls l l -==+ 3.1.2.速度和加速度分析式(1)和(2)对时间取一次二次导数写成矩阵形式，即得一下速度方程：=加速度方程=+3.2对滑枕5的受力分析重力G=50*9.8;切削阻力f=10000;质量m=50;aE=2m s;13.16/4.计算机仿真及结果程序1(运动动画模拟):clear all;clc;clear;close all;global w1 l1 l3 l4 l6 l61 H k theta N s str=inputdlg({'请输入刨床速度(转/分):', ...'请输入曲柄长度l1(m):', ...'请输入刨床的最大行程H(m):', ...'请输入摆杆与滑枕的连接杆长度l4(m):', ...'请输入行程速比系数k:', ...'是否保存动画:是(Y)否(N)'}, ...'Please input!',1,{'60','0.1','0.5','0.2','2','N'}); w1=2*pi*str2num(str{1})/60;l1=str2num(str{2});H=str2num(str{3});l4=str2num(str{4});k=str2num(str{5});theta=(pi/180)*180*(k-1)/(k+1);l3=H/(2*sin(theta/2));l6=l1/sin(theta/2);l61=l4*sin(pi/6)+l3*cos(theta/2);syms Y N;s=char(str{6});g=50*9.8;f=10000;m=50;%w1=2*pi;l1=100;l3=500;l6=200;l61=100+500*cos(pi/6);l4=20 0;for m11=1:361theta1(m11)=pi*(m11-1)/180;%theta1theta3(m11)=atan((l6+l1*sin(theta1(m11)))/(l1*cos(theta1( m11))));%theta3if theta3(m11)>=0theta3(m11)=theta3(m11);else theta3(m11)=pi+theta3(m11);end;s3(m11)=(l1*cos(theta1(m11)))/cos(theta3(m11));o4(m11)=pi-asin((l61-l3*sin(theta3(m11)))/l4);%转角4 xe(m11)=l3*cos(theta3(m11))+l4*cos(o4(m11));%E点的横坐标if theta1(m11)==pi/2 theta3(m11)=pi/2;s3(m11)=l1+l6;endif theta1(m11)==3*pi/2 theta3(m11)=pi/2;s3(m11)=l6-l1;endA1=[cos(theta3(m11)),-s3(m11)*sin(theta3(m11)),0,0; ...sin(theta3(m11)),s3(m11)*cos(theta3(m11)),0,0; ...0,-l3*sin(theta3(m11)),-l4*sin(o4(m11)),-1; ...0,l3*cos(theta3(m11)),l4*cos(o4(m11)),0];B1=w1*[-l1*sin(theta1(m11));l1*cos(theta1(m11));0;0]; D1=A1\B1;E1(:,m11)=D1;ds(m11)=D1(1);w3(m11)=D1(2);w4(m11)=D1(3);ve(m11)=D1(4);A2=[cos(theta3(m11)),-s3(m11)*sin(theta3(m11)),0,0; ...sin(theta3(m11)),s3(m11)*cos(theta3(m11)),0,0; ... 0,-l3*sin(theta3(m11)),-l4*sin(o4(m11)),-1; ...0,l3*cos(theta3(m11)),l4*cos(o4(m11)),0];B2=-[-w3(m11)*sin(theta3(m11)),(-ds(m11)*sin(theta3(m11)) -s3(m11)*w3(m11)*cos(theta3(m11))),0,0; ...w3(m11)*cos(theta3(m11)),(ds(m11)*cos(theta3(m11))-s3(m11 )*w3(m11)*sin(theta3(m11))),0,0; ...0,-l3*w3(m11)*cos(theta3(m11)),-l4*w4(m11)*cos(o4(m11)),0 ; ...0,-l3*w3(m11)*sin(theta3(m11)),-l4*w4(m11)*sin(o4(m11)),0 ]*[ds(m11);w3(m11);w4(m11);ve(m11)]; ...C2=w1*[-l1*w1*cos(theta1(m11));-l1*w1*sin(theta1(m11));0; 0]; ...B=B2+C2;D2=A2\B;E2(:,m11)=D2;dds(m11)=D2(1);a3(m11)=D2(2) ;a4(m11)=D2(3);ae(m11)=D2(4);F(m11)=(f-m*ae(m11)*10^(-3))/cos(pi-o4(m11));Fn(m11)=g-F(m11)*sin(pi-o4(m11));end;theta1=theta1*180/pi;theta3=theta3*180/pi;theta4=o4*180/pi;figure(1);subplot(311);h11=plot(theta1,theta3);xlabel('\it\theta1');ylabel('\it\ theta3');grid on;title('位置曲线图');subplot(312);h12=plot(theta1,theta4);xlabel('\it\theta1');ylabel('\it\ theta4');grid on;subplot(313);h13=plot(theta1,xe);xlabel('\it\theta1');ylabel('xe');gri d on;figure(2);subplot(311)h21=plot(theta1,w3);xlabel('\it\theta1');ylabel('\it\omeg a3');grid on;title('速度曲线图');subplot(312)h22=plot(theta1,w4);xlabel('\it\theta1');ylabel('\it\omeg a4');grid on;subplot(313);h23=plot(theta1,ve);xlabel('\it\theta1');ylabel('Ve');gri d on;figure(3);subplot(311);h31=plot(theta1,a3);xlabel('\it\theta1');ylabel('\it\alph a3');grid on;title('加速度曲线图');subplot(312);h32=plot(theta1,a4);xlabel('\it\theta1');ylabel('\it\alph a4');grid on;subplot(313);h33=plot(theta1,ae);xlabel('\it\theta1');ylabel('aE');gri d on;%力学分析figure(4);plot(theta1,F);xlabel('\it\theta1');ylabel('F45');grid on; title('杆4对滑枕5的力');figure(5);plot(theta1,Fn);xlabel('\it\theta1');ylabel('Fn');grid on; title('滑枕5受到的支持力');%分析摆杆3的受力%动画hd=pi/180;du=180/pi;for n1=1:90;theta1(n1)=(n1-1)*hd;theta3(n1)=atan((l6+l1*sin(theta1(n1)))/(l1*cos(theta1(n1 ))));theta4(n1)=asin((l61-l3*sin(theta3(n1)))/l4);endfor n1=91:270;theta1(n1)=(n1-1)*hd;theta3(n1)=pi-atan((l6+l1*sin(theta1(n1)))/(l1*cos(pi-the ta1(n1))));theta4(n1)=asin((l61-l3*sin(theta3(n1)))/l4);endfor n1=272:361;theta1(n1)=(n1-1)*hd;theta3(n1)=atan((l6+l1*sin(theta1(n1)))/(l1*cos(theta1(n1 ))));theta4(n1)=asin((l61-l3*sin(theta3(n1)))/l4);endfigure(6);m11=moviein(16);j=0;hmax=-1.6107e+03;hmin=0;for n1=1:270;j=j+1;clf;x(1)=0;y(1)=0;%画A点x(2)=l1*cos(theta1(n1));y(2)=l1*sin(theta1(n1));%画B点 x(3)=0;y(3)=-l6;%画D点x(4)=l3*cos(theta3(n1));y(4)=l3*sin(theta3(n1))-l6;%画E点x(5)=l3*cos(theta3(n1))-l4*cos(theta4(n1));y(5)=l61-l6;%画S点x(6)=l1*cos(theta1(n1))+50*cos(theta3(n1))/1000+25*sin(th eta3(n1))/1000;%画滑块2的右上角点y(6)=l1*sin(theta1(n1))+50*sin(theta3(n1))/1000-25*cos(th eta3(n1))/1000;%画滑块2的右上角点x(7)=x(6)-100*cos(theta3(n1))/1000;%画滑块2的右下角y(7)=y(6)-100*sin(theta3(n1))/1000;%画滑块2的右下角点x(8)=x(7)-50*sin(theta3(n1))/1000;%画滑块2的左下角点y(8)=y(7)+50*cos(theta3(n1))/1000;%画滑块2的左下角点x(9)=x(6)-50*sin(theta3(n1))/1000;%画滑块2的左上角点y(9)=y(6)+50*cos(theta3(n1))/1000;%画滑块2的左上角点x(10)=x(5)+500/1000;y(10)=y(5);%画滑枕x(11)=x(5)-500/1000;y(11)=y(5);%画滑枕与刀头的连接点x(12)=x(5)-500/1000;y(12)=y(10)+20/1000;%画刀头右上角x(13)=x(5)-500/1000;y(13)=y(10)-80/1000;%画刀头右下角x(14)=x(5)-500/1000-20/1000;y(14)=y(10)+20/1000;%画刀头左上角x(15)=x(5)-500/1000-20/1000;y(15)=y(10)-80/1000;%画刀头左下角x(16)=x(5)-500/1000-20/1000-20/1000;y(16)=y(10)-110/1000; %画刀尖L(n1)=sqrt((x(5)-x(4))^2+(y(5)-y(4))^2);if hmin>=x(16)hmin=x(16);n11=n1;endk=1:2;plot(x(k),y(k));hold on;k=3:4;plot(x(k),y(k));hold on;k=4:5;plot(x(k),y(k));hold on;k=6:9;plot(x(k),y(k));hold on;line([x(1) x(3)],[y(1) y(3)]);line([x(6) x(9)],[y(6) y(9)]);line([x(10) x(11)],[y(10) y(11)]);line([x(12) x(13)],[y(12) y(13)]);line([x(12) x(14)],[y(12) y(14)]);line([x(14) x(15)],[y(14) y(15)]);line([x(15) x(16)],[y(15) y(16)]);line([x(13) x(16)],[y(13) y(16)]);grid on;plot(x(1),y(1),'o');plot(x(2),y(2),'o');plot(x(3),y(3),'o');plot(x(4),y(4),'o');plot(x(5),y(5),'o');axis([-1100/1000 800/1000 -500/1000 500/1000]);m11(j)=getframe;endfor n1=272:361;j=j+1;clf;x(1)=0;y(1)=0;%画A点x(2)=l1*cos(theta1(n1));y(2)=l1*sin(theta1(n1));%画B点 x(3)=0;y(3)=-l6;%画D点x(4)=l3*cos(theta3(n1));y(4)=l3*sin(theta3(n1))-l6;%画E点x(5)=l3*cos(theta3(n1))-l4*cos(theta4(n1));y(5)=l61-l6;%画S点x(6)=l1*cos(theta1(n1))+50*cos(theta3(n1))/1000+25*sin(th eta3(n1))/1000;%画滑块2的右上角点y(6)=l1*sin(theta1(n1))+50*sin(theta3(n1))/1000-25*cos(th eta3(n1))/1000;%画滑块2的右上角点x(7)=x(6)-100*cos(theta3(n1))/1000;%画滑块2的右下角y(7)=y(6)-100*sin(theta3(n1))/1000;%画滑块2的右下角点x(8)=x(7)-50*sin(theta3(n1))/1000;%画滑块2的左下角点y(8)=y(7)+50*cos(theta3(n1))/1000;%画滑块2的左下角点x(9)=x(6)-50*sin(theta3(n1))/1000;%画滑块2的左上角点y(9)=y(6)+50*cos(theta3(n1))/1000;%画滑块2的左上角点x(10)=x(5)+500/1000;y(10)=y(5);%画滑枕x(11)=x(5)-500/1000;y(11)=y(5);%画滑枕与刀头的连接点x(12)=x(5)-500/1000;y(12)=y(10)+20/1000;%画刀头右上角x(13)=x(5)-500/1000;y(13)=y(10)-80/1000;%画刀头右下角x(14)=x(5)-500/1000-20/1000;y(14)=y(10)+20/1000;%画刀头左上角x(15)=x(5)-500/1000-20/1000;y(15)=y(10)-80/1000;%画刀头左下角x(16)=x(5)-500/1000-20/1000-20/1000;y(16)=y(10)-110/1000; %画刀尖if hmax<=x(16)hmax=x(16);n12=n1;endL(n1)=sqrt((x(5)-x(4))^2+(y(5)-y(4))^2);k=1:2;plot(x(k),y(k));hold on;k=3:4;plot(x(k),y(k));hold on;k=4:5;plot(x(k),y(k));hold on;k=6:9;plot(x(k),y(k));hold on;line([x(6) x(9)],[y(6) y(9)]);grid on;plot(x(1),y(1),'o');plot(x(2),y(2),'o');plot(x(3),y(3),'o');plot(x(4),y(4),'o');plot(x(5),y(5),'o');line([x(1) x(3)],[y(1) y(3)]);line([x(6) x(9)],[y(6) y(9)]);line([x(10) x(11)],[y(10) y(11)]);line([x(12) x(13)],[y(12) y(13)]);line([x(12) x(14)],[y(12) y(14)]);line([x(14) x(15)],[y(14) y(15)]);line([x(15) x(16)],[y(15) y(16)]);line([x(13) x(16)],[y(13) y(16)]);axis([-1100/1000 800/1000 -500/1000 500/1000]);m11(j)=getframe;endh=hmax-hmin;movie(m11,2,100);disp '*******压力角是:'90-du*theta3(n12)if s==N;elsemovie2avi(m11,'ex_movie2avi','compression','none'); end运行界面:运行结果如下：05010015020025030035040050100150θ1θ3(°)位置曲线图050100150200250300350400140160180θ1θ4(°)050100150200250300350400-0.500.5θ1x e (m )050100150200250300350400-10010θ1ω3(r a d /s )速度曲线图050100150200250300350400-505θ1ω4(r a d /s )050100150200250300350400-505θ1V e (m /s )050100150200250300350400-1000100θ1α3(r a d 2/s )加速度曲线图050100150200250300350400-2000200θ1α4(r a d 2/s )050100150200250300350400-50050θ1a E (m 2/s )5.图解法验证如下图(已知曲柄在右侧极限位置时对应图上横坐标取330):050100150200250300350400-10010X: 330Y: -2.22e-15θ1ω3(r a d /s )速度曲线图050100150200250300350400-505X: 330Y: -3.205e-15θ1ω4(r a d /s )050100150200250300350400-505X: 330Y: 1.282e-15θ1V e (m /s )当横坐标都取300时，w3，w4，Ve 都为1510-在计算机中忽略为0，与图解法算的结果一致。