第二节 铰链四杆机构的力分析
在如图2-4所示的铰链四杆机构中，已知各构件的尺寸和质心的位置、各构件的质量和转动惯量、原动件1的方位角θ1和匀角速度ω1以及构件3的工作阻力矩Mr，求各运动副中的反力和原动件上的平衡力矩Mb。
图2-4 铰链四杆机构受力分析
一、数学模型的建立
1.惯性力和惯性力矩的计算
由第一章介绍的运动分析方法可求出铰链四杆机构各构件的位移、速度和加速度，并可进一步计算出各构件质心的加速度。
构件1质心S1的加速度
（2-4）
构件2质心S2的加速度
（2-5）
构件3质心S3的加速度
（2-6）
由构件质心的加速度和构件的角加速度可以确定其惯性力和惯性力矩
（2-7）
2.平衡方程的建立
由图2-4所示的铰链四杆机构的受力分析可知，该机构有4个运动副，每个运动副反力可分解为x、y方向的两个分力，另外还有一个待求的平衡力矩共9个未知量，需列出九个方程式求解。根据前述的方法，建立各个构件的力平衡方程。
构件1受惯性力、构件2和构件4对它的作用力以及平衡力矩。对其质心S1点取矩，根据∑
=0、∑Fx=0和∑Fy=0，写出如下平衡方程
（2-8）
同理，对构件2进行受力分析，并对其质心S2点取矩，写出如下平衡方程
（2-9）
同理，对构件3进行受力分析，并对其质心S3点取矩，写出如下平衡方程
（2-10）
根据以上九个方程式可以解出各运动副反力和平衡力矩等九个未知量，由于以上九个方程式都为线性方程，为便于MATLAB编程求解，将以上线性方程组合写成矩阵形式的平衡方程
（2-11）
式中，C为系数矩阵；FR为未知力列阵；D为已知力列阵。其中
二、计算实例
【例2-1】 如图2-4所示，已知铰链四杆机构各构件的尺寸为：l1=400mm，l2=1000mm，l3=700mm，l4=1200mm，各杆质心都在杆的中点处，各构件的质量为：m1=1.2kg，m2=3kg，m3=2.2kg，各构件的转动惯量为：J1=0.016kg·m2，J2=0.25kg·m2，J3=0.09kg·m2，构件3的工作阻力矩为Mr=100N·m，顺时针方向，其他构件外力及外力矩不计，构件1以匀角速度ω1=10rad/s逆时针方向转动，不计摩擦时，求各转动副中的反力及平衡力矩Mb。
三、程序设计
铰链四杆机构力分析程序crank_rocker_force文件。
********************************************************
%1.输入已知数据
clear；
l1=0.40；
l2=1；
l3=0.70；
l4=1.200；
las1=0.2
lbs2=0.5；
lds3=0.35
omega1=10；
hd=pi/180；
du=180/pi；
%l=67.2；
m1=1.2；
m2=3；
m3=2.2；
g=10；
Js1=0.016；
Js2=0.25；
Js3=0.09；
Mr=100
%2.铰链四杆机构运动分析
for n1=1：360
L=sqrt（14*14﹢11*11-2*11*14*cos（n1*hd））；
phi（n1）=asin（（11/L）*sin（n1*hd））；
beta（n1）=acos（（-12*12﹢13*13﹢L*L）/（2*13*L））；
if beta（n1）<0
beta（n1）=beta（n1）﹢pi；
end
theta3（n1）=pi-phi（n1）-beta（n1）； %theta3表示杆3转过角度
theta2（n1）=asin（（13*sin（theta3（n1））-11*sin（n1*hd））/12）； %theta2表示杆2转过角度
omega3（n1）=omega1*（11*sin（（n1*hd-theta2（n1））））/（13*sin（（theta3（n1）-theta2（n1））））；
%omega3表示杆3角速度
omega2（n1）=-omega1*（11*sin（（n1*hd-theta3（n1））））/（12*sin（（theta2（n1）-theta3（n1））））；
%omega2表示杆2角速度
alpha3（n1）=（omega1∧2*11*cos（（n1*hd-theta2（n1）））﹢omega2（n1）∧2*12-omega3（n1）∧2*13
*cos（（theta3（n1）-theta2（n1））））/（13*sin（（theta3（n1）-theta2（n1））））；
%alpha3表示杆3角加速度
alpha2（n1）=（-omega1∧2*11*cos（（n1*hd-theta3（n1）））﹢omega3（n1）∧2*13-omega2（n1）∧2*12
*cos（（theta2（n1）-thena3（n1））））/（12*sin（（theta2（n1）-thena（n1））））
%alpha2表示杆2角加速度
%3.铰链四杆机构力平衡计算
%计算质心速度
as1x（n1）=-las1*omega1∧2*cos（n1*hd）； %质心s1在x轴的加速度
as1y（n1）=-las1*omega1∧2*sin（n1*hd）； %质心s1在y轴的加速度
as2x（n1）=-l1*omega1∧2*cos（n1*hd）-lbs2*（alpha2（n1）*sin（theta2（n1））
﹢omega2（n1）∧2*cos（theta2（n1）））； %质心s2在x轴的加速度
as2y（n1）=-l1*omega1∧2*sin（n1*hd）﹢lbs2*（alpha2（n1）*cos（theta2（n1））
-omega2（n1）∧2*sin（theta2（n1）））； %质心s2在y轴的加速度
as3x（n1）=-lds3*（cos（theta3（n1））*omega3（n1）∧2﹢sin（theta3（n1））*alpha3（n1））；
%质心s3在x轴的加速度
as3y（n1）=-lds3*（sin（theta3（n1））*omega3（n1）∧2-cos（theta3（n1））*alpha3（n1））；
%计算构件惯性力和惯性力矩
F1x（n1）=-as1x（n1）*m1；
F1y（n1）=-as1y（n1）*m1；
F2x（n1）=-as2x（n1）*m2；
F2y（n1）=-as2y（n1）*m2；
F3x（n1）=-as3x（n1）*m3；
F3y（n1）=-as3y（n1）*m3；
M2（n1）=-alpha2（n1）*Js2； %作用于杆2的合力矩
M3（n1）=-alpha3（n1）*Js3-Mr； %作用于杆3的合力矩
%计算各个铰链点坐标， 计算各个质心点坐标
xa=0；
ya=0；
xb（n1）=l1*cos（n1*hd）；
yb（n1）=l1*sin（n1*hd）；
xc（n1）=l4﹢l3*cos（theta3（n1））；
yc（n1）=l3*sin（theta3（n1））；
xd=l4；
yd=0；
xs1（n1）=（xb（n1）﹢xa）/2；
ys1（n1）=（yb（n1）﹢ya）/2；
xs2（n1）=（xb（n1）﹢xc（n1））/2；
ys2（n1）=（yb（n1）﹢yc（n1））/2；
xs3（n1）=（xc（n1）﹢xd）/2；
ys3（n1）=（yc（n1）﹢yd）/2；
%未知力系数矩阵
A=zeros（9）；
A（1，1）=1；A（1，2）=-（ys1（n1）-ya）；A（1，3）=-（xa-xs1（n1））；
A（1，4）=-（ys1（n1）-yb（n1））；A（1，5）=-（xb（n1）-xs1（n1））；
A（2，2）=-1；A（2，4）=-1；
A（3，3）=-1；A（3，5）=-1；
A（4，4）=（ys2（n1）-yb（n1））；A（4，5）=（xb（n1）-xs2（n1））；
A（4，6）=-（ys2（n1）-yc（n1））；A（4，7）=-（xc（n1）-xs2（n1））；
A（5，4）=1；A（5，6）=-1；
A（6，5）=1；A（6，7）=-1；
A（7，6）=（ys3（n1）-yc（n1））；A（7，7）=（xc（n1）-xs3（n1））；
A（7，8）=-（ys3（n1）-yd）；A（7，9）=-（xd-xs3（n1））；
A（8，6）=1；A（8，8）=-1；
A（9，7）=1；A（9，9）=-1；
%已知力列
B=zeros（9，1）；
B（1）=0；
B（2）=-F1x（n1）；
B（3）=-F1y（n1）﹢m1*g；%B（3）=-F1y（n1）；
B（4）=-M2（n1）；
B（5）=-F2x（n1）；
B（6）=-F2y（n1）﹢m2*g； %B（6）=-F2y（n1）；
B（7）=-M3（n1）；
B（8）=-F3x（n1）；
B（9）=-F3y（n1）﹢m3*g；%B（9）=-F3y（n1）；
C=A\B；
Mb（n1）=C（1）；Fr14x（n1）=C（2）；Fr14y（n1）=C（3）；Fr12x（n1）=C（4）；Fr12y（n1）=C（5）；
Fr23x（n1）=C（6）；Fr23y（n1）=C（7）；Fr34x（n1）=C（8）；Fr34y（n1）=C（9）；
end；
%4.铰链四杆机构力分析图形输出
figure（4）；
n1=1：360；
subplot（2，2，1）； %绘运动副反力FR14曲线图
plot（n1， Fr14x，'b'）；
hold on
plot（n1，Fr14y，'k--'）；
legend（'F_R_1_4_x'，'F_R_1_4_y'）
title（'运动副反力F_R_1_4曲线图'）；
xlabel（'曲柄转角 \theta_1/\circ'）
ylabel（'F/N'）
grid on；
subplot（2，2，2）； %绘运动副反力FR23曲线图
plot（n1，Fr23x，'b'）；
hold on
plot（n1，Fr23y，'k--'）；
hold on
legend（'F_R_2_3_x'，'F_R_2_3_y'）
title（'运动副反力F_R_2_3曲线图'）；
xlabel（'曲柄转角 \theta_1/\circ'）
ylabel（'F/N'）
grid on；
subplot（2，2，3）； %绘运动副反力FR34曲线图
plot（n1，Fr34x，'b'）；
hold on
plot（n1，Fr34y，'k--'）；
hold on
legend（'F_R_3_4_x'，'F_R_3_4_y'）
title（'运动副反力F_R_3_4_x曲线图'）；
xlabel（'曲柄转角 \theta_1/\circ'）
ylabel（'F/N'）
grid on；
subplot（2，2，4）； %绘平衡力矩M曲线图
plot（n1，Mb）
title（'力矩Mb图'）
xlabel（'曲柄转角 \theta_1/\circ'）；
ylabel（'M/N.m'）
hold on；
grid on；
text（100，1.9*10∧6，'Mb'）
四、运算结果
图2-5为铰链四杆机构的力分析线图。
图2-5 铰链四杆机构力分析线图