在如图2-6所示的曲柄滑块机构中，已知各构件的尺寸和质心的位置、各构件的质量和转动惯量、原动件1的方位角θ1和匀角速度ω1以及滑块3的水平工作阻力Fr，求各运动副中的反力和原动件上的平衡力矩Mb。

图2-6　曲柄滑块机构受力分析

由第一章介绍的运动分析方法可求出曲柄滑块机构各构件的位移、速度和加速度，并可进一步计算出各构件质心的加速度。

构件1质心S1的加速度

　　（2-12）

构件2质心S2的加速度

　　（2-13）

构件3质心S3的加速度

　　（2-14）

由构件质心的加速度和构件的角加速度可以确定其惯性力和惯性力矩

　　（2-15）

曲柄滑块机构有3个铰链，每个铰链受杆的作用分别有x、y方向的两个分力，另外还有一个移动副反力和一个待求的平衡力矩共8个未知量，需列出八个方程式求解。

如图2-6所示，对构件1进行受力分析，构件1受惯性力、构件2和构件4对它的作用力以及平衡力矩。对其质心S1点取矩，根据∑=0、∑Fx=0和∑Fy=0，写出如下平衡方程

　　（2-16）

同理，对构件2进行受力分析，并对其质心S2点取矩，写出如下平衡方程

　　（2-17）

同理，对构件3进行受力分析，根据∑Fx=0和∑Fy=0，写出如下平衡方程

　　（2-18）

根据以上八个方程式可以解出运动副反力和平衡力矩等八个未知量，由于以上八个方程式都为线性方程，为便于MATLAB编程求解，将以上线性方程组合写成矩阵形式的平衡方程

　　（2-19）

式中，C为系数矩阵；FR为未知力列阵；D为已知力列阵。其中

【例2-2】　在图2-6所示的曲柄滑块机构中，已知各构件的尺寸分别为l1=400mm，l2=1200mm， =200mm， =600mm，ω1=10rad/s，各构件的质量及转动惯量分别为：m1=1.2kg，m2=3.6kg，m3=6kg， =0.45kg·m2，滑块3上作用外力Fr=-1000N，求各运动副中的反力及原动件1的平衡力矩Mb。

曲柄滑块机构力分析程序slider_crank__force文件

********************************************************

%1.输入已知数据

clear；

l1=0.400；

l2=1.200；

las1=0.2；

lbs2=0.6；

m1=1.2；

m2=3.6；

m3=6；

g=10；

J2=0.45；

G1=m1*g；

G2=m2*g；

G3=m3*g；

Fr=-1000；

e=0；

hd=pi/180；

du=180/pi；

omega1=10；

alpha1=0；

%2.曲柄滑块机构力平衡计算

for n1=1：360

　theta1（n1）=（n1-1）*hd；

　% 调用函数 slider_crank 计算曲柄滑块机构位移，速度，加速度

［theta2（n1），s3（n1），omega2（n1），v3（n1），alpha2（n1），a3（n1）］=slider_crank（theta1（n1），omega1，alpha1，l1，l2，e）；

　% 计算各个质心点加速度

　as1x（n1）=-las1*cos（theta1（n1））*omega1∧2；

　as1y（n1）=-las1*sin（n1*hd）*omega1∧2；

　as2x（n1）=-l1*omega1∧2*cos（n1*hd）-lbs2*（omega2（n1）∧2*cos（theta2（n1））﹢alpha2（n1）*sin（theta2（n1）））；　%质心s2在x轴的加速度

　as2y（n1）=-l1*omega1∧2*sin（n1*hd）-lbs2*（omega2（n1）∧2*sin（theta2（n1））-alpha2（n1）*cos（theta2（n1）））；　%质心s2在y轴的加速度

　% 计算各构件惯性力和惯性力矩

　F1x（n1）=-as1x（n1）*m1；

　F1y（n1）=-as1y（n1）*m1；

　F2x（n1）=-as2x（n1）*m2；

　F2y（n1）=-as2y（n1）*m2；

　F3x（n1）=-a3（n1）*m3；

　F3y（n1）=0；

　FR43x（n1）=Fr；

　M2（n1）=-alpha2（n1）*J2；

　% 计算各个铰链点坐标， 计算各个质心点坐标

　xa=0；

　ya=0；

　xs1=las1*cos（n1*hd）；

　ys1=las1*sin（n1*hd）；

　xb=l1*cos（n1*hd）；

　yb=l1*sin（n1*hd）；

　xs2=xb﹢lbs2*cos（theta2（n1））；

　ys2=yb﹢lbs2*sin（theta2（n1））；

　xc=xb﹢l2*cos（theta2（n1））；

　yc=yb﹢l2*sin（theta2（n1））；

　　% 未知力系数矩阵

　　A=zeros（8）；

　　A（1，1）=1；A（1，2）=-（ys1-yb）；A（1，3）=-（xb-xs1）；A（1，4）=-（ys1-ya）；

　　A（1，5）=-（xa-xs1）；

　　A（2，2）=-1；A（2，4）=-1；

　　A（3，3）=-1；A（3，5）=-1；

　　A（4，2）=（ys2-yb）；A（4，3）=（xb-xs2）；A（4，6）=-（ys2-yc）；A（4，7）=-（xc-xs2）；

　　A（5，2）=1；A（5，6）=-1；

　　A（6，3）=1；A（6，7）=-1；

　　A（7，6）=1；

　　A（8，7）=1；A（8，8）=-1；

　　%已知力列阵

　　B=zeros（8，1）；

　　B（2）=-F1x（n1）；

　　B（3）=-F1y（n1）﹢G1；

　　B（4）=-M2（n1）；

　　B（5）=-F2x（n1）；

　　B（6）=-F2y（n1）﹢G2；

　　B（7）=-F3x（n1）﹢FR43x（n1）；

　　B（8）=-F3y（n1）；

　　C=A\B；

　　Mb（n1）=C（1）；Fr12x（n1）=C（2）；Fr12y（n1）=C（3）；Fr14x（n1）=C（4）；Fr14y（n1）=C（5）；

　　Fr23x（n1）=C（6）；Fr23y（n1）=C（7）；Fr34y（n1）=C（8）；

end

　　%3.曲柄滑块机构力分析图形输出

　figure（2）；

　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（n1），'b'）；

　hold on

　plot（n1，Fr23y（n1），'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，Fr34y，'b'）；

　hold on

　legend（'F_R_3_4_y'）

　title（'运动副反力F_R_3_4_y曲线图'）；

　xlabel（'曲柄转角 \theta_1/\circ'）

　ylabel（'F/N'）

　grid on；

　subplot（2，2，4）；　　%绘平衡力矩Mb曲线图

　plot（n1，Mb）

　title（'力矩Mb图'）

　xlabel（'曲柄转角 \theta_1/\circ'）；

　ylabel（'M/N.m'）

　hold on；

　grid on；

图2-7为曲柄滑块机构的力分析线图。

图2-7　曲柄滑块机构力分析线图