连杆机构运动分析.doc

1. 构件上点的运动分析
构件上点的运动分析
A
j
l
rN
N
x
O
y
B
rA
q
w，a

函数文件（m文件）
格式：　function [ 输出参数 ] = 函数名（ 输入参数 ）
p_crank.m function [p_Nx,p_Ny]=p_crank(Ax,Ay,theta,phi,l1)
v_crank.m　function [v_Nx,v_Ny]=v_crank(l1,v_Ax,v_Ay,omiga,theta,phi)
a_crank.m　function [a_Nx,a_Ny]=a_crank(l1,a_Ax,a_Ay,alpha,omiga,theta,phi)

函数中的符号说明
函数中的符号
说明
p_Nx，p_Ny
v_Nx，v_Ny
a_Nx，a_Ny
l
Ax，Ay
v_Ax，v_Ay
a_Ax，a_Ay
theta
phi
omiga
alpha
N点位置在x、y轴方向的分量
N点速度在x、y轴方向的分量
N点加速度在x、y轴方向的分量
AN杆长度
A点位置在x、y轴方向的分量
A点速度在x、y轴方向的分量
A点加速度在x、y轴方向的分量
AB杆与x轴夹角
AN与AB杆夹角（构件结构参数）
AB杆角速度
AB杆角加速度

2. RRRⅡ级杆组运动分析

实线位置，m=1
虚线位置，m= -1
C
D
B
C ¢
C
j
d
rD
D
x
O
y
B
rB
q2
w2，a2
w3，a3
l3
3
2
q3
l2
a
RRRⅡ级杆组运动分析

函数文件（m文件）
格式：　function [ 输出参数 ] = 函数名（ 输入参数 ）
p_RRR.m function [cx,cy,theta2,theta3]=p_RRR(bx,by,dx,dy,l2,l3,m)
v_RRR.m　 function [vcx,vcy,omiga2,omiga3]=v_RRR(vbx,vby,vdx,vdy,cx,cy,bx,by,dx,dy)
a_RRR.m function [acx,acy,alpha2,alpha3]=a_RRR(abx,aby,adx,ady,cx,cy,bx,by,dx,dy,omiga2,omiga3)

函数中的符号说明
函数中的符号
说明
bx，by
cx，cy
dx，dy
vbx，vby
vcx，vcy
vdx，vdy
abx，aby
acx，acy
adx，ady
l2，l3
theta2，theta3
omiga2，omiga3
alpha2，alpha3
m
B点位置在x、y轴方向的分量
C点位置在x、y轴方向的分量
D点位置在x、y轴方向的分量
B点速度在x、y轴方向的分量
C点速度在x、y轴方向的分量
D点速度在x、y轴方向的分量
B点加速度在x、y轴方向的分量
C点加速度在x、y轴方向的分量
D点加速度在x、y轴方向的分量
BC杆和CD杆的长度
BC杆和CD杆与x轴的夹角
BC杆和CD杆的角速度
BC杆和CD杆的角加速度
杆组的装配模式系数

3. RRPⅡ级杆组运动分析
C
s²r
P
B
l2
s¢r
C ¢
d
C
s²r
P
B
l2
s¢r
C ¢
d
l2<d
实线位置，ÐBCP< 90°，m=1
虚线位置，ÐBC¢P> 90°，m=-1
l2>d
ÐBCP< 90°，ÐBC¢P> 90°，m=1
w3，a3
w2，a2
2
d
C
sr
rB
P
x
O
y
B
rP
q2
RRPⅡ级杆组运动分析

rC
q3
l2
3

函数文件（m文件）
格式：　function [ 输出参数 ] = 函数名（ 输入参数 ）
p_RRP.m function [cx,cy,sr,theta2]=p_RRP(bx,by,px,py,theta3,l2,m)
v_RRP.m　　function [vcx,vcy,vr,omiga2]=v_RRP(bx,by,cx,cy,vbx,vby,vpx,vpy,theta2,theta3,l2,sr,omiga3)
a_RRP.m function [acx,acy,ar,alpha2]=a_RRP(bx,by,cx,cy,px,py,abx,aby,apx,apy,theta3,vr,omiga2,omiga3,alpha3)
函数中的符号说明
函数中的符号
说明
bx，by
cx，cy
px，py
vbx，vby
vcx，vcy
vpx，vpy
abx，aby
acx，acy
apx，apy
l2
theta2，theta3
omiga2，omiga3
alpha2，alpha3
sr，vr，ar
m
B点位置在x、y轴方向的分量
C点位置在x、y轴方向的分量
P点位置在x、y轴方向的分量
B点速度在x、y轴方向的分量
C点速度在x、y轴方向的分量
P点速度在x、y轴方向的分量
B点加速度在x、y轴方向的分量
C点加速度在x、y轴方向的分量
P点加速度在x、y轴方向的分量
BC杆的长度
BC杆和滑块导路与x轴的夹角
BC杆和滑块导路的角速度
BC杆和滑块导路杆的角加速度
滑块上C点相对于参考点P的位移、速度、加速度
杆组的装配模式系数

4. RPRⅡ级杆组运动分析
w3，a3
2
f
D
sr
d
N¢
x
O
y
B
q2
RRPⅡ级杆组运动分析

C
q3
l3
3
实线位置，m=1
虚线位置，m= -1

e
N

函数文件（m文件）
格式：　function [ 输出参数 ] = 函数名（ 输入参数 ）
p_RPR.m function [dx,dy,sr,theta3]=p_RPR(bx,by,cx,cy,e,l3,m)
v_RPR.m　　 function [vdx,vdy,omiga3,vr]=v_RPR(bx,by,cx,cy,dx,dy,vcx,vcy,vbx,vby,theta3)
a_RPR.m function [adx,ady,alpha3,ar]=a_RPR(bx,by,cx,cy,dx,dy,acx,acy,abx,aby,vr,omiga3,theta3)
函数中的符号说明
函数中的符号
说明
bx，by
cx，cy
dx，dy
vbx，vby
vcx，vcy
vdx，vdy
abx，aby
acx，acy
adx，ady
l3
e
theta3
omiga3
alpha3
sr，vr，ar
m
B点位置在x、y轴方向的分量
C点位置在x、y轴方向的分量
D点位置在x、y轴方向的分量
B点速度在x、y轴方向的分量
C点速度在x、y轴方向的分量
D点速度在x、y轴方向的分量
B点加速度在x、y轴方向的分量
C点加速度在x、y轴方向的分量
D点加速度在x、y轴方向的分量
ND杆的长度
内接移动副导路至B点的偏心距
ND杆与x轴的夹角
ND杆的角速度
ND杆的角加速度
滑块上C点相对于N点的位移、速度、加速度
杆组的装配模式系数

5. RRRⅡ级杆组力分析
RRRⅡ级杆组力分析

x
O
y
D
B
F2
2
3
S2
M2
C
S3
M3
F3
B
F2y
2
S2
M2
C
-R23x
-R23y
F2x
R12y
R12x
F3x
D
3
C
S3
M3
R23x
R23y
-R34x
-R34y
F3y

函数文件（m文件）
格式：　function [ 输出参数 ] = 函数名（ 输入参数 ）
F_RRR.m function [R12x,R12y,R23x,R23y,R34x,R34y]=F_RRR(bxy,cxy,dxy,s2,s3,m2,m3, Js2,Js3,M2,M3,F2,F3,as2,as3,alpha2,alpha3)
函数中的符号说明
函数中的符号
说明
bxy
cxy
dxy
s2，s3
m2，m3
Js2，Js3
M2，M3
F2，F3
as2，as3
alpha2，alpha3
R12x，R12y，R23x，R23y
R34x，R34y
B点位置在x、y轴方向的分量
C点位置在x、y轴方向的分量
D点位置在x、y轴方向的分量
构件2、构件3的质心位置在x、y轴方向的分量
构件2、构件3的质量
构件2、构件3绕其质心的转动惯量
作用在构件2、构件3质心上的主矩
作用在构件2、构件3质心上的主矢
构件2、构件3质心加速度在x、y轴方向的分量
构件2、构件3的角加速度
运动副B、C、D中的约束反力

6. RRPⅡ级杆组力分析
RRPⅡ级杆组力分析

x
O
y
C
F2
2
S2
M2
B
S3
M3
F3
3
-R23x
F2x
R12y
R12x
-R23y
F2y
C
2
S2
M2
B
q3
N
R23y
-R34
C
S3
M3
F3y
3
F3x
R23x

函数文件（m文件）
格式：　function [ 输出参数 ] = 函数名（ 输入参数 ）
F_RRP.m function [R12x,R12y,R23x,R23y,R34x,R34y,lcn]=F_RRP(bxy,cxy,s2,s3,m2,m3, Js2,Js3,M2,M3,F2,F3,theta3,as2,as3,alpha2,alph3)
函数中的符号说明
函数中的符号
说明
bxy
cxy
s2，s3
m2，m3
Js2，Js3
M2，M3
F2，F3
theta3
as2，as3
alpha2，alpha3
R12x，R12y，R23x，R23y
R34x，R34y
lcn
B点位置在x、y轴方向的分量
C点位置在x、y轴方向的分量
构件2、构件3的质心位置在x、y轴方向的分量
构件2、构件3的质量
构件2、构件3绕其质心的转动惯量
作用在构件2、构件3质心上的主矩
作用在构件2、构件3质心上的主矢
滑块导路与x轴的夹角
构件2、构件3质心加速度在x、y轴方向的分量
构件2、构件3的角加速度
运动副B、C中的约束反力
移动副的约束反力
力R34的作用点位置

7. RPRⅡ级杆组力分析
RPRⅡ级杆组力分析

x
O
y
D
C
F3
S2
M3
B
S3
M2
F2
2
-R34
3
N
D
F3x
M3
B
S3
-R34x
3
-R34y
R23y
R23x
-R35x
-R35y
F3y
q3
R12x
N
S2
M2
F2x
2
-R23
F2y
R12y
N
N
C

函数文件（m文件）
格式：　function [ 输出参数 ] = 函数名（ 输入参数 ）
F_RPR.m function [R12x,R12y,R23x,R23y,R35x,R35y,lcn]=F_RRP(bxy,cxy,dxy,s2,s3, m2,m3,Js2,Js3,M2,M3,F2,F3,R34,theta3,as2,as3,alpha3)
函数中的符号说明
函数中的符号
说明
bxy
cxy
dxy
s2，s3
m2，m3
Js2，Js3
M2，M3
F2，F3
R34
theta3
as2，as3
alpha3
R12x，R12y，R35x，R35y
R23x，R23y
lcn
B点位置在x、y轴方向的分量
C点位置在x、y轴方向的分量
D点位置在x、y轴方向的分量
构件2、构件3的质心位置在x、y轴方向的分量
构件2、构件3的质量
构件2、构件3绕其质心的转动惯量
作用在构件2、构件3质心上的主矩
作用在构件2、构件3质心上的主矢
杆组外构件作用在构件3上的约束反力
ND杆与x轴的夹角
构件2、构件3质心加速度在x、y轴方向的分量
构件3的角加速度
运动副C、B中的约束反力
滑块的约束反力
力R23的作用点位置

8. 作用有平衡力的构件力分析
N¢
A
F1y
1
S1
M1
F1x
R01y
R01x
作用有平衡力的构件力分析

x
O
y
-R12x
B
Mb
-R12y

函数文件（m文件）
格式：　function [ 输出参数 ] = 函数名（ 输入参数 ）
F_Bar.m function [R01x,R01y,Mb]=F_Bar(axy,bxy,s1,m1,Js1,M1,F1,R12,as1,alpha1)
函数中的符号说明
函数中的符号
说明
axy
bxy
s1
m1
Js1
M1
F1
as1
alpha1
R12x，R12y
R01x，R01y
Mb
A点位置在x、y轴方向的分量
B点位置在x、y轴方向的分量
构件1的质心位置在x、y轴方向的分量
构件1的质量
构件1绕其质心的转动惯量
作用在构件1质心上的主矩
作用在构件1质心上的主矢
构件1质心加速度在x、y轴方向的分量
构件1的角加速度
作用在运动副B中的约束反力
运动副A的约束反力
平衡力矩

9. 平面连杆机构运动分析算例
曲柄摇杆机构
l1
x
A
y
B
D
C
l2
l4
l3
例1图示曲柄摇杆机构，已知l1=150mm，l2=220mm，l3=250mm，l4=300mm，曲柄以n1=100r/min逆时针匀速转动，分析该机构的运动。
主程序
% 曲柄摇杆机构运动分析
clc,clear
l1=150;% 曲柄长度
l2=220;% 连杆长度
l3=250;% 摇杆长度
l4=300;% 机架长度
n=100;% 曲柄转速
m=1;% RRR II级杆组装配模式系数
omiga1=2*pi*n/60;alpha1=0;% 曲柄角速度、角加速度
ax=0;ay=0;vax=0;vay=0;aax=0;aay=0;% A点位置、速度、加速度
dx=l4;dy=0;vdx=0;vdy=0;adx=0;ady=0;% D点位置、速度、加速度
phi=0;% 曲柄的结构参数
theta1=0:10:360;% 曲柄转角(每隔10°计算一次)
theta1=theta1*pi/180;
% 调用crank函数，计算B点运动参数
[bx,by]=p_crank(ax,ay,theta1,phi,l1);
[vbx,vby]=v_crank(l1,vax,vay,omiga1,theta1,phi);
[abx,aby]=a_crank(l1,aax,aay,alpha1,omiga1,theta1,phi);
% 调用RRR函数，计算BC杆和CD杆以及C点运动参数
[cx,cy,theta2,theta3]=p_RRR(bx,by,dx,dy,l2,l3,m);
[vcx,vcy,omiga2,omiga3]=v_RRR(vbx,vby,vdx,vdy,cx,cy,bx,by,dx,dy);
[acx,acy,alpha2,alpha3]=a_RRR(abx,aby,adx,ady,cx,cy,bx,by,dx,dy,omiga2,omiga3);
% 绘制运动线图
theta1=theta1*180/pi;
figure(1)
subplot(3,1,1);
plot(theta1,cx,'-',theta1,cy,':r'),grid on
xlabel('曲柄转角( ° )');ylabel('位移(mm/s)');
legend('C点x方向位移','C点y方向位移');
subplot(3,1,2);
plot(theta1,vcx,'-',theta1,vcy,':r'),grid on
xlabel('曲柄转角( ° )');ylabel('速度(mm/s)');
legend('C点x方向速度','C点y方向速度');
subplot(3,1,3);
plot(theta1,acx,'-',theta1,acy,':r'),grid on
xlabel('曲柄转角( ° )');ylabel('加速度(mm/s^2)');
legend('C点x方向加速度','C点y方向加速度');
figure(2)
subplot(2,1,1);
plot(theta1,theta3*180/pi),grid on
xlabel('曲柄转角( ° )');ylabel('摇杆角位移( ° )');
subplot(2,1,2);
plot(theta1,omiga3),grid on
xlabel('曲柄转角( ° )');ylabel('摇杆角速度(rad/s)');
subplot(2,1,3);
plot(theta1,alpha3),grid on
xlabel('曲柄转角( ° )');ylabel('摇杆角加速度(/rad/s^2)');

曲柄滑块机构
e
x
A
y
B
O
C
q2
sr
l2
q1
例2　图示曲柄滑块机构，已知l1=150mm，l2=150mm，e=25mm，曲柄以n1=955r/min逆时针匀速转动，分析该机构的运动。
主程序
% 曲柄滑块机构运动分析
clc,clear
l1=50;% 曲柄长度
l2=150;% 连杆长度
e=25;% 偏距
n=955;% 曲柄转速
m=1;% RRP II级杆组装配模式系数
omiga1=2*pi*n/60;alpha1=0;% 曲柄角速度、角加速度
ax=0;ay=e;vax=0;vay=0;aax=0;aay=0;% A点位置、速度、及速度
px=0;py=0;%　滑块导路上一定点（选为O点）的坐标
vpx=0;vpy=0;apx=0;apy=0;% 滑块导路上一定点的速度、加速度
theta3=0;omiga3=0;alpha3=0;% 滑块导路角位置、角速度、角加速度
phi=0;% 曲柄的结构参数
theta1=0:30:360;% 曲柄转角(每隔10°计算一次)
theta1=theta1*pi/180;
% 调用crank函数，计算B点运动参数
[bx,by]=p_crank(ax,ay,theta1,phi,l1);
[vbx,vby]=v_crank(l1,vax,vay,omiga1,theta1,phi);
[abx,aby]=a_crank(l1,aax,aay,alpha1,omiga1,theta1,phi);
% 调用RRP函数，计算BC杆和滑块的运动参数
[cx,cy,sr,theta2]=p_RRP(bx,by,px,py,theta3,l2,m);
[vcx,vcy,vr,omiga2]=v_RRP(bx,by,cx,cy,vbx,vby,vpx,vpy,theta2,theta3,l2,sr,omiga3);
[acx,acy,ar,alpha2]=a_RRP(bx,by,cx,cy,px,py,abx,aby,apx,apy,theta3,vr,omiga2,omiga3,alpha3);
% 绘制运动线图
crankx=l1.*cos(theta1);cranky=e+l1.*sin(theta1);
theta1=theta1*180/pi;
figure(1)
subplot(3,1,1);
plot(theta1,cx,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('滑块位移(mm/s)');
subplot(3,1,2);
plot(theta1,vcx,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('滑块速度(mm/s)');
subplot(3,1,3);
plot(theta1,acx,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('滑块加速度(mm/s^2)');

3
D
摆动导杆机构
h
x
A
B
C
1
q1
sr
l1
q3
2
例3 图示摆动导杆机构，已知l1=280mm，h=380mm，lCD=840mm，曲柄1以w1=18.012rad/逆时针匀速转动，分析该机构的运动。
主程序
% RPR II级杆组运动分析程序检验
clc,clear
l1=280;% 曲柄长度
h=380;% 机架长度
lcd=840;% 导杆长度
e=0;% 内接移动副导路至C点的偏心距
m=1;% RPR II级杆组装配模式系数
phi=0;% 曲柄的结构参数
omiga1=18.021;alpha1=0;% 曲柄角速度、角加速度
ax=0;ay=h;vax=0;vay=0;aax=0;aay=0;% A点位置、速度、加速度
cx=0;cy=0;vcx=0;vcy=0;acx=0;acy=0;% C点位置、速度、加速度
theta1=(0:5:360)*pi/180;% 曲柄转角(每隔5°计算一次)
% 调用crank函数，计算B点运动参数
[bx,by]=p_crank(ax,ay,theta1,phi,l1);
[vbx,vby]=v_crank(l1,vax,vay,omiga1,theta1,phi);
[abx,aby]=a_crank(l1,aax,aay,alpha1,omiga1,theta1,phi);
% 调用RPR函数，计算BC杆运动参数
[dx,dy,sr3,theta3]=p_RPR(cx,cy,bx,by,e,lcd,m);
[vdx,vdy,omiga3,vr3]=v_RPR(cx,cy,bx,by,dx,dy,vbx,vby,vcx,vcy,theta3);
[adx,ady,alpha3]=a_RPR(cx,cy,cx,by,dx,dy,abx,aby,acx,acy,vr3,omiga3,theta3);
% 绘制运动线图
theta1=theta1*180/pi;
subplot(3,1,1);
plot(theta1,theta3*180/pi,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('摇杆角位移( ° )');
subplot(3,1,2);
plot(theta1,omiga3,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('摇杆角速度(rad/s)');
subplot(3,1,3);
plot(theta1,alpha3,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('摇杆角加速度(rad/s^2)');

D
E
F
y
5
B
六杆机构
x
A
C
6
w1
H2
H1
4
3
2
1
例4 图示六杆机构，已知lAB=80mm，lBC=260mm，lDE=400mm，lCD=300mm，lEF=460mm，H1=90mm，H2=170mm，曲柄1逆时针匀速转动，w1=40rad/s，分析滑块5的运动。
主程序
% 平面六杆机构运动分析程序检验
clc,clear
l1=80;l2=260;l3=300;l31=400;l4=460;% 各构件长度
h1=90;h2=170;
e=0;% 内接移动副导路至C点的偏心距
m1=1;m2=1;% RRR、RRP II级杆组装配模式系数
phi1=0;phi2=0;% 曲柄和构件DE的结构参数
omiga1=40;alpha1=0;% 曲柄角速度、角加速度
ax=0;ay=h1;vax=0;vay=0;aax=0;aay=0;% A点位置、速度、加速度
dx=h2;dy=0;vdx=0;vdy=0;adx=0;ady=0;% D点位置、速度、加速度
px=0;py=h1;vpx=0,vpy=0;apx=0;apy=0;% 滑块导路上一定点（取为A点）运动参数
theta5=0;omiga5=0;alpha5=0;% 滑块导路角运动参数
theta1=(0:10:360)*pi/180;% 曲柄转角(每隔10°计算一次)
% 调用crank函数，计算B点运动参数
[bx,by]=p_crank(ax,ay,theta1,phi1,l1);
[vbx,vby]=v_crank(l1,vax,vay,omiga1,theta1,phi1);
[abx,aby]=a_crank(l1,aax,aay,alpha1,omiga1,theta1,phi1);
% 调用RRR函数，计算C点运动参数
[cx,cy,theta2,theta3]=p_RRR(bx,by,dx,dy,l2,l3,m1);
[vcx,vcy,omiga2,omiga3]=v_RRR(vbx,vby,vdx,vdy,cx,cy,bx,by,dx,dy);
[acx,acy,alpha2,alpha3]=a_RRR(abx,aby,adx,ady,cx,cy,bx,by,dx,dy,omiga2,omiga3);
% 调用crank函数，计算E点运动参数
[ex,ey]=p_crank(dx,dy,theta3,phi2,l31);
[vex,vey]=v_crank(l31,vdx,vdy,omiga3,theta3,phi2);
[aex,aey]=a_crank(l31,adx,ady,alpha3,omiga3,theta3,phi2);
% 调用RRP函数，计算F点运动参数
[fx,fy,sr,theta4]=p_RRP(ex,ey,px,py,theta5,l4,m2);
[vfx,vfy,vr,omiga4]=v_RRP(ex,ey,fx,fy,vex,vey,vpx,vpy,theta4,theta5,l4,sr,omiga5);
[afx,afy,ar,alpha4]=a_RRP(ex,ey,fx,fy,px,py,aex,aey,apx,apy,theta4,vr,omiga4,omiga5,alpha5);
% 绘制运动线图
theta1=theta1*180/pi;
subplot(3,1,1);
plot(theta1,sr,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('滑块5位移(mm)');
subplot(3,1,2);
plot(theta1,vr,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('滑块5速度(mm/s)');
subplot(3,1,3);
plot(theta1,ar,'r'),grid on
xlabel('曲柄转角( ° )');ylabel('滑块5加速度(mm/s^2)');

10. 平面连杆机构力分析算例
y
5
6
E
H1
4
1
3
D
平面六杆机构
x
A
B
C
n1
2
H
图示平面六杆机构，已知l1=0.1m，l3=0.7m，l4=0.45m，H=0.65 m，H1=0.5 m，曲柄1逆时针匀速转动，n1=100 r/min，构件质量m3=20 kg，m4=15 kg，m5=62 kg，构件3和构件4的质心分别位于各自构件长度的一半处，且构件3、4绕其质心的转动惯量分别为Js3=0.11 kg×m2，Js4=0.18 kg×m2，，构件1和构件2的质量忽略不计，该机构在工作行程时滑块受生产阻力Fr =110 N。分析滑块5的运动规律及在工作中应加在曲柄1上的平衡力矩。
主程序
% 六杆机构力分析
clc,clear
lab=0.1; % AB杆长度
lcd=0.7; % CD杆长度
lac=0.5; % 机架长度
lde=0.45; %DE杆长度
phi=0;% 曲柄的结构参数
ax=0;ay=lac;vax=0;vay=0;aax=0;aay=0;% 转动副A位置、速度、加速度
H=0.65; %　滑块5至铰链C垂直距离
cx=0;cy=0;vcx=0;vcy=0;acx=0;acy=0;% C点位置、速度、加速度
px=0;py=H;vpx=0;vpy=0;apx=0;apy=0;% 滑块导路上一定点P的位置、速度、加速度
s1=[0,0];as1=[0,0];% 曲柄质心位置和加速度
theta5=0;omiga5=0;alpha5=0;% 滑块导路的角位置、角速度、角加速度
m1=1;% RPR II级杆组装配模式系数
m2=1;% RRP II级杆组装配模式系数
e=0;% 内接移动副导路至C点的偏心距
n1=100; % 曲柄转速
omiga1=pi*n1/30;alpha1=0;% 曲柄角速度、角加速度
% 曲柄AB的运动参数
theta1=0:10:360; 曲柄转角(每隔10°计算一次)
theta1=theta1*pi/180;
% 调用crank函数，计算B点运动参数
[bx,by]=p_crank(ax,ay,theta1,phi,lab);
[vbx,vby]=v_crank(lab,vax,vay,omiga1,theta1,phi);
[abx,aby]=a_crank(lab,aax,aay,alpha1,omiga1,theta1,phi);
% 调用RPR函数，计算CD杆和D点运动参数
[dx,dy,sr3,theta3]=p_RPR(cx,cy,bx,by,e,lcd,m1);
[vdx,vdy,omiga3,vr3]=v_RPR(cx,cy,bx,by,dx,dy,vbx,vby,vcx,vcy,theta3);
[adx,ady,alpha3,ar3]=a_RPR(cx,cy,bx,by,dx,dy,abx,aby,acx,acy,vr3,omiga3,theta3);
% 调用RRP函数，计算E点运动参数
[ex,ey,sr,theta4]=p_RRP(dx,dy,px,py,theta5,lde,m2);
[vex,vey,vr,omiga4]=v_RRP(dx,dy,ex,ey,vdx,vdy,vpx,vpy,theta4,theta5,lde,sr,omiga5);
[aex,aey,ar,alpha4]=a_RRP(dx,dy,ex,ey,px,py,adx,ady,apx,apy,theta5,vr,omiga4,omiga5,alpha5);
% 调用crank函数，计算CD杆和DE杆质心点的运动参数
[s3x,s3y]=p_crank(dx,dy,theta3,phi,lcd/2);
[vs3x,vs3y]=v_crank(lcd/2,vcx,vcy,omiga3,theta3,phi);
[as3x,as3y]=a_crank(lcd/2,acx,acy,alpha3,omiga3,theta3,phi);
[s4x,s4y]=p_crank(dx,dy,theta4,phi,lcd/2);
[vs4x,vs4y]=v_crank(lcd/2,vdx,vdy,omiga4,theta4,phi);
[as4x,as4y]=a_crank(lcd/2,adx,ady,alpha4,omiga4,theta4,phi);
% 力参数
m1=0;m2=0;m3=20;m4=15;m5=62;% 各构件质量
Js1=0;Js2=0;Js3=0.11;Js4=0.18;Js5=0;% 各构件绕其质心的转动惯量
M1=0;M2=0;M3=0;M4=0;M5=0;% 作用于构件1、2、3、4、5的主矩
F1=[0,0];F2=[0,0];F3=[0,0];F4=[0,0];% 作用于构件1、2、3、4的主矢
fr=110;fr5=(omiga3>0)*fr;% 生产阻力
for i=1:length(theta1)
F5=[fr5(i),0];% 作用于构件5的主矢
% 力分析
% 调用RRP II级杆组力分析函数，计算运动副D、E中的约束反力
[R34x(i),R34y(i),R45x(i),R45y(i),R56x(i),R56y(i),lcn]=F_RRP([dx(i),dy(i)],[ex(i),ey(i)],...
[s4x(i),s4y(i)],[ex(i),ey(i)],m4,m5,Js4,Js5,M4,M5,F4,F5,theta5,...
[as4x(i),as4y(i)],[aex(i),aey(i)],alpha4(i),alpha5);
% 调用RPR II级杆组力分析函数，计算运动副B、C中的约束反力
[R12x(i),R12y(i),R23x(i),R23y(i),R36x(i),R36y(i),lcn(i)]=F_RPR([cx,cy],[bx(i),by(i)],...
[dx(i),dy(i)],[bx(i),by(i)],[s3x(i),s3y(i)],m2,m3,Js2,Js3,M2,M3,F2,F3,...
[R34x(i),R34y(i)],theta3(i),[abx(i),aby(i)],[as3x(i),as3y(i)],alpha3(i));
% 调用单杆力分析函数，计算运动副A中的约束反力及平衡力矩
[R16x(i),R16y(i),Mb(i)]=F_Bar([ax,ay],[bx(i),by(i)],s1,m1,Js1,M1,F1,...
[R12x(i