牛头刨床的连杆机构运动分析0 前言机构运动分析的任务是对于结构型式及尺寸参数已定的具体机构，按主动件的位置、速度和加速度来确定从动件或从动件上指定点的位置、速度和加速度。
许多机械的运动学特性和运动参数直接关系到机械工艺动作的质量，运动参数又是机械动力学分析的依据，所以机构的运动分析是机械设计过程中必不可少的重要环节。
以计算机为手段的解析方法，由于解算速度快，精确度高，程序有一定的通用性，已成为机构运动分析的主要方法。
连杆机构作为在机械制造特别是在加工机械制造中主要用作传动的机构型式，同其他型式机构特别是凸轮机构相比具有很多优点。
连杆机构采用低副连接，结构简单，易于加工、安装并能保证精度要求。
连杆机构可以将主动件的运动通过连杆传递到与执行机构或辅助机构直接或间接相连的从动件，实现间歇运动，满足给定的运动要求，完成机器的工艺操作。
牛头刨床是一种利用工作台的横向运动和纵向往复运动来去除材料的一种切削加工机床。
工作台的纵向往复运动是机床的主运动，实现工件的切削。
工作台的横向运动即是进给运动，实现对切削精度的控制。
本文中只分析纵向运动的运动特性。
牛头刨床有很多机构组成，其中实现刨头切削运动的六连杆机构是一个关键机构。
刨床工作时，通过六杆机构驱动刨刀作往复移动。
刨刀右行时，当刨刀处于工作行程时;要求刨刀的速度较低且平稳，以减小原动机的容量和提高切削质量。
当刨刀处于返回行程时，刨刀不工作，称为空行程，此时要求刨刀的速度较高以提高生产率。
由此可见，牛头刨床的纵向运动特性对机床的性能有决定性的影响。
1 牛头刨床的六连杆机构牛头刨床有很多机构组成，其中实现刨头切削运动的六杆机构是一个关键机构。
图1所示的为一牛头刨床的六连杆机构。
杆1为原动件，刨刀装在C点上。
假设已知各构件的尺寸如表1所示，原动件1以等角速度ω1=1rad/s沿着逆时针方向回转，要求分析各从动件的角位移、角速度和角加速度以及刨刀C点的位移、速度和加速度的变化情况。
图1 牛头刨床的六连杆机构简图表1 六连杆机构的尺寸参数（单位：mm ）l 1 l 3 l 4 h h 1 h 2 1809601609004601102 六连杆机构的运动分析方程杆件1为主动件，六杆机构的运动随杆件1的位置变化而发生周期性变化。
在一个变化周期中，可以把杆件1的角位置分成36等分（θ1取362 πn ,其中n 取整数0～35,对应时间s 362πn t =），分别研究θ1在不同取值下杆件机构的位置参数和运动参数的变化。
3s 的长度与刨刀的运动行程成正比，因此可以用5S 表征刨刀的行程，用5S 关于时间的一阶导数5s来表征刨刀的运动速度，用5S 关于时间的二阶导数5s 来表征刨刀的加速度。
1）位置方程由图可知2θ=3θ，故未知量有3θ、4θ（x 轴与4l所成的的角度）、3S （直线BD 的长度）、5S （直线GC 的长度）。
利用两个封闭图形ABDEA 和EDCGE ，建立两个封闭矢量方程，由此可得：⎭⎬⎫+=++=+' s l l s l l l l56431643 (1) 把(1)式分别向x 轴、y 轴投影得：⎪⎪⎭⎪⎪⎬⎫=+=++=++=+ h l l s l l l h s l l h s l 33445334411133441123344sin sin cos cos sin sin sin cos cos cos θθθθθθθθθθ(2) 在（2）式中包含3s 、5s 、3θ、4θ四个未知数，消去其中三个可得到只含4θ一个未知数的方程：[][]{}[][]0sin sin sin 2sin cos cos sin sin 244111234424224244112244111=-+--+-++-+ θθθθθθθθl l h l hl h l l l h l l h (3)当1θ取不同值时，用牛顿迭代法解(3)式，可以求出每个4θ的值，再根据方程组(2)可以求出其他杆件的位置参数3s 、5s 、3θ的值：⎪⎪⎭⎪⎪⎬⎫-+=+=-= 3441113334453443sin sin sin cos cos )sin arcsin(θθθθθθθl l h s l l s l l h (4)2）速度方程对(2)式对时间求一次导数并把结果写成矩阵的形式得：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----00cos sin 0cos cos 01sin sin 00cos cos sin 0sin sin cos 1111143443344334433344333θθωωωθθθθθθθθθθl l v v l l l l l s l s C e B (5)其中C v 为刨刀的水平速度，v eB 为滑块2相对于杆3的速度。
由于每个1θ对应的3s 、3θ、4θ已求出，方程组式（5）的系数矩阵均为常数，采用按列选主元的高斯消去法可求解（式Ⅳ）可解得角速度ω3、ω4、eB v 、C v 。
3）加速度方程把式Ⅳ对时间求导得矩阵式：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--+⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡----------=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----00sin cos 0sin sin 00cos cos 00sin sin cos cos 0cos cos sin sin 0cos cos 01sin sin 00cos cos sin 0sin sin cos 11111114344433344433344433333344433333343443344334433344333θωθωωωωθωθωθωθωθωθωθθωθωθωθθωααl l v v l l l l l s v l s v a a θl θl θl θl θl θs θθl θs θC eBe B eBC e B（6）同样采用按列选主元的高斯消去法可求解（6）可得角加速度3α、4α、eB a 、C a 。
3 程序流程图4 Matlab程序编写%主程序开始l1=180;l3=960;l4=160;h=900;h1=460;h2=110;theta1=linspace(0,35*pi/18,36);%定义常量和已知参数，l1代表杆1的长度，l3代表杆3的长度，l4代表杆4的长度,h表示EG 的长度，h1表示AE的竖直距离,h2表示AE的水平距离，theta1表示角θ1的不同值。
theta3=zeros(1,36);theta4=zeros(1,36);s3=zeros(1,36);s5=zeros(1,36);t est=zeros(1,36);vBe=zeros(1,36);vc=zeros(1,36);omega1=ones(1,36);omeg a3=zeros(1,36);omega4=zeros(1,36);aBe=zeros(1,36);ac=zeros(1,36);alph a1=zeros(1,36);alpha3=zeros(1,36);alpha4=zeros(1,36);A=zeros(4,4);dA=zeros(4,1);%定义最终的结果数据，当θ1取不同值时,theta3表示θ3的值，theta4表示θ4的值，s3表示BD的长度，s5表示GC的长度，vBe表示B点在杆3上运动的速度，vc表示杆5的运动速度，即牛头刨刀的速度，omega3表示杆3的转动角速度，omega4表示杆4的转动角速度，aBe表示B点在杆3上运动的角加速度，ac表示杆5的加速度，即牛头刨刀的加速度，alpha3表示杆3的角加速度，alpha4表示杆4的角加速度，矩阵A,dA表示线性方程组的系数矩阵。
i=0; %i为循环变量，在循环结构中使用。
syms THETA1THETA4 %定义符号变量，为以下计算做准备。
fun1=((h1+l1*sin(THETA1)-l4*sin(THETA4))^2+(h2+l1*cos(THETA1)-l4*cos(THETA4))^2)*(l4^2*sin(THETA4)^2+h^2-2*h*l4*sin(THETA4))-l3^2*(h1+ l1*sin(THETA1)-l4*sin(THETA4))^2;%定义迭代法中要求解的关于THETA4的方程。
x0=0; %定义在牛顿迭代法中的变量THEA4的初值。
for i=1:36 %用循环结构求当theta1取不同值时，theta3值。
fun2=subs(fun1,THETA1,theta1(i));%把不同的THETA1的值代入要求解的方程。
[theta4(i),EA,it]=NEWTON(fun2,'THETA4',x0,0.0001,1000);%用牛顿迭代法求得THEATA4，并赋值到theta4的数组中。
x0=theta4(i); %把这次计算的解作为下一次计算的初值。
endfor i=1:36%用循环结构求当theta1的值取不同值时，theta3、s3、s5的取值。
因为theta3的值可能的取值范围为[0,π]，对theta3求解时应分以下两种情况讨论。
if sign(h2+l1*cos(theta1(i))-l4*cos(theta4(i)))>0 %theta3<π/2 theta3(i)=asin((h-l4*sin(theta4(i)))/l3);else theta3(i)=pi-asin((h-l4*sin(theta4(i)))/l3); %theta3>π/2 endtest(i)=h1+l1*sin(theta1(i))-l4*sin(theta4(i));s5(i)=l4*cos(theta4(i))+l3*cos(theta3(i));s3(i)=(h1+l1*sin(theta1(i))-l4*sin(theta4(i)))/sin(theta3(i)); endfor i=1:36%用循环结构求当theta1的值取不同值时vBe、omega3、omega4、vc的值。
A(1,1)=cos(theta3(i));A(1,2)=-s3(i)*sin(theta3(i));A(1,3)=-l4*sin(theta4(i));A(2,1)=sin(theta3(i));A(2,2)=s3(i)*cos(theta3(i));A(2,3)=l4*cos(theta4(i));A(3,2)=-l3*sin(theta3(i));A(3,3)=-l4*sin(theta4(i));A(4,2)=l3*cos(theta3(i));A(4,3)=l4*cos(theta4(i));dA(1,1)=-omega1(1,1)*l1*sin(theta1(i));dA(2,1)=omega1(1,2)*l1*cos(theta1(i));x=gauss(A,dA);%用按列选主元的高斯消去法求解。