牛头刨床运动分析实例.doc

例： 如图所示为一牛头刨床的机构运动简图。设已知各构件尺寸为：，，，原动件1的方位角和等角速度。试用矩阵法求该机构中各从动件的方位角、角速度和角加速度以及E点的位移、速度和家速度的运动线图。

解：先建立一直角坐标系，并标出各杆矢量及方位角。其中共有四个未知量、、及。为求解需建立两个封闭矢量方程，为此需利用两个封闭图形ABCA及CDEGC，由此可得，
(1-1)
写成投影方程为：
(1-2)
解上面方程组，即可求得、、及四个位置参数，其中。
将上列各式对时间取一次、二次导数，并写成矩阵形式，即可得以下速度和加速度方程式。速度方程式： (1-3)
机构从动件的位置参数矩阵：
机构从动件的的速度列阵：
机构原动件的位置参数矩阵：
：机构原动件的角速度
加速度方程式：
(1-4)
机构从动件的位置参数矩阵求导：
机构从动件的的加速度列阵：
机构原动件的位置参数矩阵求导： 主程序（matlab）：

%牛头刨床运动分析主程序
%x(1)——代表;
%x(2)—— 代表构件3的转角；
%x(3）——代表构件4的转角；
%x(4)——代表E点的线位移;
%x(5)——代表；
%x(6)——代表;
%x(7)——代表;
%x(8)——代表;
%x(9）——代表;
%x(10)——代表构件1的转角。
x=[ 0.302 65*pi/180 169*pi/180 0.1 0.125 0.6 0.15 0.275 0.575 0];%赋初值
dr=pi/180;%度转化为弧度
dth=10*dr;w1=1;%每10度计算一个点
for i=1:37
y=ntpc(x); %调用从动件位置方程求解函数ntpc(自编)
s3=y(1);theta3=y(2);theta4=y(3);se=y(4); %得到位置参数。
%将各位置参数用向量储存，便于后面绘图，角度用度表示
ss3(i)=y(1);th1(i)=x(10)/dr;th3(i)=y(2)/dr;th4(i)=y(3)/dr;sse(i)=y(4);

%进行速度分析
A=[cos(theta3) -s3*sin(theta3) 0 0;
sin(theta3) s3*cos(theta3) 0 0;
0 -x(6)*sin(theta3) -x(7)*sin(theta4) -1;
0 x(6)*cos(theta3) x(7)*cos(theta4) 0];%A机构从动件的位置参数矩阵
B=[-x(5)*sin(x(10));x(5)*cos(x(10));0;0];%B机构原动件的位置参数列阵
yy=w1*inv(A)*B;%公式1-3求解，yy表示机构从动件速度列阵，inv(A)是A的逆阵
vs3=yy(1);w3=yy(2);w4=yy(3);vse=yy(4);
%将各速度参数以向量的方式表示，以便后面绘图
dvs3(i)=yy(1);dw3(i)=yy(2);dw4(i)=yy(3);dvse(i)=yy(4);
%dA为从动件位置参数矩阵对时间一次求导

%进行角速度分析
dA=[-w3*sin(theta3) -vs3*sin(theta3)-s3*w3*cos(theta3) 0 0;
w3*cos(theta3) vs3*cos(theta3)-s3*w3*sin(theta3) 0 0;
0 -x(6)*w3*cos(theta3) -x(7)*w4*cos(theta4) 0;
0 -x(6)*w3*sin(theta3) -x(7)*w4*sin(theta4) 0];
%dB就是原动件位置参数列阵对时间一次求导
dB=[-x(5)*w1*cos(x(10));-x(5)*w1*sin(x(10));0;0];
KK=-dA*yy+w1*dB; %KK为公式1-4右端
ya=inv(A)*KK;%公式1-4求解，ya为从动件加速度列阵
%将各加速度以向量表示
as3(i)=ya(1);atheta3(i)=ya(2);atheta4(i)=ya(3);ase(i)=ya(4);
x(10)=x(10)+dth; %计算下一个点
x(1)=s3;
x(2)=theta3;
x(3)=theta4;
x(4)=se;
end

%绘制运动参数曲线
subplot(2,2,1); % 选择第1个子窗口
plot(th1,th3,th1,th4,th1,sse*1e3);%绘制位置线图
subplot(2,2,2);
plot(th1,dw3,th1,dw4,th1,dvse);%绘制速度线图
subplot(2,2,3);
plot(th1,atheta3,th1,atheta4,th1,ase); %绘制加速度线图
%这个函数是关于牛头刨床位置方程求解，可得：s3,theta3,theta4,sE四个运动变量。
子程序：
function y=ntpc(x)
%x(1)——代表;
%x(2)—— 代表构件3的转角；
%x(3）——代表构件4的转角；
%x(4)——代表E点的线位移;
%x(5)——代表；
%x(6)——代表;
%x(7)——代表;
%x(8)——代表;
%x(9）——代表;
%x(10)——代表构件1的转角。
%先赋初值；这些初值来自于主程序。
s3=x(1);
theta3=x(2);
theta4=x(3);
se=x(4);
epsilon=1e-6;%设置求解精度为10-6
%用矩阵的形式表示位置方程组（4x1的矩阵）
f=[s3*cos(theta3)-x(5)*cos(x(10));s3*sin(theta3)-x(5)*sin(x(10))-x(8);
x(6)*cos(theta3)+x(7)*cos(theta4)-se; x(6)*sin(theta3)+x(7)*sin(theta4)-x(9)];
%用牛顿-辛普森法求解
while norm(f)>epsilon
%J位置方程组的雅可比矩阵，即从动件位置参数矩阵
J=[cos(theta3) -s3*sin(theta3) 0 0;
sin(theta3) s3*cos(theta3) 0 0;
0 -x(6)*sin(theta3) -x(7)*sin(theta4) -1;
0 x(6)*cos(theta3) x(7)*cos(theta4) 0];
dth=inv(J)*(-1.0*f); %计算增量，进行迭代，inv(J) 为J的逆阵
s3=s3+dth(1);
theta3=theta3+dth(2);
theta4=theta4+dth(3);
se=se+dth(4);
f=[s3*cos(theta3)-x(5)*cos(x(10));s3*sin(theta3)-x(5)*sin(x(10))-x(8);
x(6)*cos(theta3)+x(7)*cos(theta4)-se; x(6)*sin(theta3)+x(7)*sin(theta4)-x(9)];
norm(f);%若未达精度，会继续迭代。
end
%输出4个参数
y(1)=s3;
y(2)=theta3;
y(3)=theta4;
y(4)=se;