% 曲柄摇杆机构运动分析
% (1)-计算连杆的输出角th3和摇杆的输出角th4
% 设定各杆的长度(单位:毫米)
rs(1)=65.0; % 设定机架1长度
rs(2)=21.0; % 设定曲柄2长度
rs(3)=132.0; % 设定连杆3长度
rs(4)=110; % 设定摇杆4长度
dr=pi/180.0;% 角度与弧度的转换系数
% 设定初始推测的输入
% 机构的初始位置
th(1)=0.0; % 设定曲柄2初始位置角是0度(与机架1共线)
th(2)=45*dr; % 连杆3的初始位置角是 45度
th(3)=135*dr; % 摇杆4的初始位置角是135度
% 摇杆4的初始位置角可以用三角形的正弦定理确定
th(3)=pi-asin(sin(th(2))*rs(3)/rs(4));
dth=5*dr; % 循环增量
% 曲柄输入角从0度变化到360度,步长为5度,计算th34
for i=1:72
[th3,th4]=ntrps(th,rs); % 调用牛顿—辛普森方程求解机构位置解非线性方程函数文件
% Store results in a matrix-th34,in degrees
% 在矩阵th34中储存结果,以度为单位;(i,:)表示第i行所有列的元素;(:,i)表示第i列所有行的元素
th34(i,:)=[th(1)/dr th3/dr th4/dr]; % 矩阵[曲柄转角 连杆转角 摇杆转角] th(1)=th(1)+dth; % 曲柄转角递增
th(2)=th3; % 连杆转角中间计算值
th(3)=th4; % 摇杆转角中间计算值
end
% 绘制输出角th(2)与th(3)—输入角th(1)的关系曲线
subplot(2,2,1) % 选择第1个子窗口
plot(th34(:,1),th34(:,2),th34(:,1),th34(:,3))
axis([0 360 0 170])
grid % 网格线
ylabel('从动件角位移/deg')
title('角位移线图')
text(110,110,'摇杆4角位移')
text(50,35,'连杆3角位移')
% (2)-计算连杆的角速度om3和摇杆的角速度om4
% Setting initial conditions
% 设置初始条件
om2=125; % 曲柄角速度(等速输入)
T=2*pi/om2; % 机构周期-曲柄旋转1周的时间(秒) % 曲柄输入角从0度变化到360度,步长为5度,计算om34
for i=1:72
ct(2)=i*dth;
A=[-rs(3)*sin(th34(i,2)*dr) rs(4)*sin(th34(i,3)*dr);
rs(3)*cos(th34(i,2)*dr) -rs(4)*cos(th34(i,3)*dr)];
B=[om2*rs(2)*sin(ct(2));-om2*rs(2)*cos(ct(2))];
om=inv(A)*B; % 输出角速度矩阵
om3=om(1);
om4=om(2);
om34(i,:)=[i om3 om4]; % 矩阵[序号 连杆角速度 摇杆角速度] t(i)=i*T/72;
end
% 绘制连杆的角速度om3和摇杆的角速度om4—时间Times的关系曲线
subplot(2,2,2) % 选择第2个子窗口
plot(t,om34(:,2),t,om34(:,3))
axis([0 0.05 -190 210])
grid % 网格线
title('角速度线图')
ylabel('从动件角速度/rad/s')
text(0.001,170,'摇杆4角速度')
text(0.013,130,'连杆3角速度')
% (3)-计算连杆的角加速度a3和摇杆的角加速度a4
a2=0; % 曲柄角速度是等速,角加速度a2=dom2/dt=0
% 曲柄输入角从0度变化到360度,步长为5度,计算a34
for i=1:72
c(2)=i*dth;
C=[-rs(3)*sin(th34(i,2)*dr) rs(4)*sin(th34(i,3)*dr);
rs(3)*cos(th34(i,2)*dr) -rs(4)*cos(th34(i,3)*dr)];
D(1)=
a2*rs(2)*sin(c(2))+om2^2*rs(2)*cos(c(2))+om34(i,2)^2*rs(3)*cos(th34(i,2)*dr)-om34(i,3)^2*rs(4)*cos(th34(i,3)*dr);
D(2)=-a2*rs(2)*cos(c(2))+om2^2*rs(2)*sin(c(2))+om34(i,2)^2*rs(3)*sin(th34(i,2)*dr)-om34(i,3)^2*rs(4)*sin(th34(i,3)*dr);
a=inv(C)*D'; % 输出角加速度矩阵
a3=a(1);
a4=a(2);
a34(i,:)=[i a3 a4]; % 矩阵[序号 连杆角加速度 摇杆加角速度] t(i)=i*T/72;
end
% 绘制连杆的角加速度a3和摇杆的角加速度a4—时间Times的关系曲线
subplot(2,2,3) % 选择第3个子窗口
plot(t,a34(:,2),t,a34(:,3))
axis([0 0.05 -6*1e4 8*1e4])
grid % 网格线
title('角加速度线图')
xlabel('时间/s')
ylabel('从动件加速度/rad/s^{2}')
text(0.003,6.2*1e4,'摇杆4角加速度')
text(0.010,3.3*1e4,'连杆3角加速度')
%
% 输出1：四杆机构运动周期(0:5:360),时间,角位移,角速度,角加速度数据
disp ' 曲柄转角 连杆转角-摇杆转角-连杆角速度-摇杆角速度-连杆加速度-摇杆加速度' ydcs=[th34(:,1),th34(:,2),th34(:,3),om34(:,2),om34(:,3),a34(:,2),a34(:,3)];
disp (ydcs)
% 输出参数的数量级必须一致
%
% (4)-运动误差分析
% 闭环矢量方程:r2+r3-r4-r1=0
% 误差矢量E=r2+r3-r4-r1的模是表示仿真有效程度的标量(ex和ey是误差分量)
ex=rs(2)*cos(th34(:,1)*dth)+rs(3)*cos(th34(:,2)*dth)-rs(4)*cos(th34(:,3)*dth)-rs(1);
ey=rs(2)*sin(th34(:,2)*dth)+rs(3)*sin(th34(:,2)*dth)-rs(4)*sin(th34(:,3)*dth);
ee=norm([ex ey]); % 计算误差矢量矩阵的范数(模)
%
% 输出2：四杆机构运动周期(0:5:360),时间,X向误差分量,Y向误差分量
disp ' 曲柄转角 时间(秒) X向误差 Y向误差'
wc=[th34(:,1),t(:),ex(:,1),ey(:,1)];
disp (wc)
fprintf (1,' 误差矢量矩阵的模 ee = %3.4f \n',ee)
%
% 绘制均方根相容性误差曲线
subplot(2,2,4) % 选择第4个子窗口
plot(t,ex(:,1),t,ey(:,1))
axis([0 0.05 -800 600])
grid % 网格线
title('均方根误差曲线')
xlabel('时间/s')
ylabel('均方根误差')
text(0.012,350,'X向误差分量')
text(0.003,-600,'Y向误差分量')