用Matlab对四连杆运动模拟引言四连杆机构因其结构灵活、能够传递动力并有效地实现预定动作，在很多领域得到了广泛应用。
进行连杆机构运动分析，传统方法主要是图解法或分析法，无论设计精度还是设计效率都相对低下，无法满足现代机械高速高精度的要求。
随着计算机技术的飞速发展，特别是以MATLAB为代表的数值计算软件的出现，为进行机构分析提供了有力的工具。
1、四连杆介绍1.1、四连杆介绍与分类所有运动副均为转动副的四杆机构称为铰链四杆机构，它是平面四杆机构的基本形式，其他四杆机构都可以看成是在它的基础上演化而来的。
选定其中一个构件作为机架之后，直接与机架链接的构件称为连架杆，不直接与机架连接的构件称为连杆，能够做整周回转的构件被称作曲柄，只能在某一角度范围内往复摆动的构件称为摇杆。
如果以转动副连接的两个构件可以做整周相对转动，则称之为整转副，反之称之为摆转副。
铰链四杆机构中，按照连架杆是否可以做整圆周转动，可以将其分为三种基本形式，即曲柄摇杆机构，双曲柄机构和双摇杆机构。
曲柄摇杆机构，两连架杆中一个为曲柄一个为摇杆的铰链四杆机构。
双曲柄机构，具有两个曲柄的铰链四杆机构称为双曲柄机构。
其特点是当主动曲柄连续等速转动时，从动曲柄一般做不等速转动。
在双曲柄机构中，如果两对边构件长度相等且平行，则成为平行四边形机构。
这种机构的传动特点是主动曲柄和从动曲柄均以相同的角速度转动，而连杆做平动。
双摇杆机构。
双摇杆机构是两连架杆均为摇杆的铰链四杆机构。
1.2、格拉霍夫定理杆长之和条件：平面四杆机构的最短杆和最长杆的长度之和小于或者等于其余两杆长度之和。
在铰链四杆机构中，如果某个转动副能够成为整转副，则它所连接的两个构件中，必有一个为最短杆，并且四个构件的长度关系满足杆长之和条件。
在有整装副存在的铰链四杆机构中，最短杆两端的转动副均为整转副。
此时，如果取最短杆为机架，则得到双曲柄机构；若取最短杆的任何一个相连构件为机架，则得到曲柄摇杆机构；如果取最短杆对面构件为机架，则得到双摇杆机构。
如果四杆机构不满足杆长之和条件，则不论选取哪个构件为机架，所得到机构均为双摇杆机构。
上述系列结论称为格拉霍夫定理。
运用条件分析：1.Lmax+Lmin＞其余两杆之和,------此铰链四杆机构为双摇杆机构;2.Lmax+Lmin≤其余两杆之和,要具体分析:(1)Lmin为机架------为双曲柄机构(2)Lmin为连架杆------为曲柄摇杆机构(Lmin为曲柄)(3)Lmin为连杆------为双摇杆机构1.3、理论分析对图1的四连杆机构列出方程组。
其中L1、L2、L3、L4为四根杆，L4为机架。
L1×cosθ1+L2×cosθ2=L3×cosθ1+L4L1×sinθ1+L2×sinθ2=L3sinθ3由此两个方程消去θ2，便可得到一个关于θ1、θ3的函数F(θ1, θ3)=（L3×cosθ3+L4−L1×cosθ1）2+（L3×sinθ3−L1×sinθ1）2−L22然后用隐函数求出θ3=θ3(θ1)另解出角速度：θ1=w1×tw2=d(θ3)/d(t)2、Matlab的实现(一)数据分析给出了3组数据：a)L1=8;L2=25;L3=20;L4=20;b)L1=7;L2=9;L3=11;L4=4;c)L1=7;L2=6.2;L3=10;L4=6;用格拉霍夫定理进行分析得到：a为曲柄摇杆机构;b为双曲柄机构；c双摇杆机构；(二)Matlab求解并进行动画模拟a)L1=8;L2=25;L3=20;L4=20; （曲柄摇杆机构）θ3=θ3(θ1)图像实时动画模拟（部分截屏）W2的图像说明：其中w2(1)为先符号求导再赋值计算，w2(2)为直接数值求导。
（下同）b)L1=7;L2=9;L3=11;L4=4; （双曲柄机构）θ3=θ3(θ1)图像实时动画模拟（部分截屏）W2的图像c)L1=7;L2=6.2;L3=10;L4=6; （双摇杆机构）θ3=θ3(θ1)图像实时动画模拟（部分截屏）W2的图像3、结果分析对θ3的角速度w2的求值得两种方法中，分析图像知：符号求导再代入数值更容易出现大的偏差（w2的图像中的红线）。
直接用数值求导得到的解更平稳。
4、程序1、函数文件dd.mfunction F=dd(y)global L1 L2 L3 L4 px=p;F=((L3*cos(y)+L4-L1*cos(x))^2)+((L3*sin(y)-L1*sin(x))^2)-(L2^2);2、计算θ3的文件y_figure.mglobal L1 L2 L3 L4 p %传递参数% L1=8;L2=25;L3=20;L4=20; %曲柄摇杆机构% L1=7;L2=9;L3=11;L4=4; %双曲柄机构L1=7;L2=6.2;L3=10;L4=6; %双摇杆机构% L1=9;L2=5;L3=10;L4=7; %双摇杆机构Ltotal=L1+L2+L3+L4;Lmax=max([L1 L2 L3 L4]);[Lmin,I]=min([L1 L2 L3 L4]);m=100; b=0;while (Lmax+Lmin)>=(Ltotal/2)str='双摇杆机构';b=1;breakendwhile (Lmax+Lmin)<(Ltotal/2)if I==4str='双曲柄机构';else if I==1||I==3str='曲柄摇杆机构';else if I==2str='双摇杆机构';b=1;endendendbreakendswitch bcase 0m=100;x=linspace(0,2*pi,m); %确定隐函数自变量的范围y0=0.7341; %第一个方程的初值case 1a=0.58;m=80;x=linspace(a,(2*pi-a),m);y0=1.2465;endy=[];f=[];for k=1:mp=x(k);[y1,fval,exitflag,output] = fzero('dd',y0);y0=y1;y=[y,y1];f=[f,fval];endfigureplot(x(1:m),y(1:m),'r.-'), %绘制隐函数图形grid on3、计算w2的文件w_figure.mglobal L1 L2 L3 L4 w1w1=100;syms xx yy%xx为输入杆转角，yy为输出角转角F=L1*cos(xx)+L2*(1-((L3*sin(yy)-L1*sin(xx))/(L2))^2)^(1/2)-L3*cos(yy)-L4; w=(-1)*(diff(F,xx)/diff(F,yy))*w1;w21=subs(w,{xx,yy},{x,y});%用隐函数求导t=x/w1;w22=diff(y)./diff(t);%直接数值求导plot(w21,'r')hold onplot(w22,'k')4、实时动画文件siliangan.mpic=figure('name',str);set(pic,'color','white');axis equalgrid onaxis([-10,30,-10,20]);axis onglobal L1 L3 L4xa=0;ya=0; %A点坐标xd=L4;yd=0; %D点坐标xb=L1*cos(x(1));yb=L1*sin(x(1));%B点坐标xc=L3*cos(y(1))+L4;yc=L3*sin(y(1));%C点坐标stick_1=line([xa;xb],[ya;yb],'color','red','linewidth',5,'linestyle','-');stick_2=line([xb;xc],[yb;yc],'color','g','linewidth',5,'linestyle','-');stick_3=line([xc;xd],[yc;yd],'color','b','linewidth',5,'linestyle','-');stick_0=line([xd;xa],[yd;ya],'color','y','linewidth',5,'linestyle','-');dot_a=line(xa,ya,'color','r','linestyle','.','markersize',30);dot_b=line(xb,yb,'color','black','linestyle','.','markersize',30);dot_c=line(xc,yc,'color','r','color',[0.1 0.7 0.3],'linestyle','.','markersize',30); dot_d=line(xd,yd,'color','r','linestyle','.','markersize',30);dt=2*pi/m;k=1;while 1if k>m;k=1;endxb=L1*cos(x(k))+xa;yb=L1*sin(x(k))+ya;xc=L3*cos(y(k))+L4;yc=L3*sin(y(k));set(stick_1,'xdata',[xa;xb],'ydata',[ya;yb]);set(stick_2,'xdata',[xc;xb],'ydata',[yc;yb]);set(stick_3,'xdata',[xd;xc],'ydata',[yd;yc]);set(dot_b,'xdata',xb,'ydata',yb);set(dot_c,'xdata',xc,'ydata',yc);set(gcf,'doublebuffer','on');drawnow;pausek=k+1;end。