2.平面机械的动态静力分析.doc

第2章 平面连杆机构的运动分析和动态静力分析为了适应现代技术的发展，机械原理课程设计方法应由传统的图解法转向以解析法为主，以计算机为主要设计手段。目前，机械原理教材中对图解法没有详细的论述。为了进行使学生能够熟练地应用解析法，有必要编制一些常用机构分析与设计的应用子程序。本指导教程主要有平面连杆机构运动分析、凸轮机构设计标准函数使用方法和实例，供课程设计使用。所用编程语言为目前流行的数学计算软件matlab，因此，要求学生必须具有matlab基本运算和编程基础。考虑到目前有关matlab的参考书较多，本书不再作详细的讲解。本教程中在平面连杆机构运动分析和动态分析，凸轮机构分析和设计中都应用了复数矢量法，为了便于学生理解、应用这一方法解决机械分析和设计问题，此处对复数矢量法的一些基本知识和运算方法作一简单介绍。2.1复数矢量的基本知识解析法对平面机构运动分析和设计的方法很多，目前国内教材中经常应用的有矩阵法、杆组法和矢量方程法等。根据相关文献的介绍，采用复数矢量方程是解决机构运动学问题的最简单有效的方法，尤其是平面问题。为了便于理解和应用，这里对复数矢量的一些知识以及复数矢量在matlab中的运算方法作一简单介绍。2.1.1 平面矢量的复数表达式 设矢量a在坐标轴x、y上的投影为a、b（图2-1），模为r，幅角为j，以弧度表示。则矢量r可表达成复数形式或 上式也可以表示a点在坐标轴上的位置，在速度和加速度分析中可用来表示某一点的速度、加速度矢量。图2-1对于一个矢量，若a=3,b=4,在matalab中其表达式为r=3+4*i或r=5*exp(0.34*i) 2.1.2 矢量的运算l 矢量的加（减）：两个矢量的加（减）法运算是它们的实部和虚部分别相加（减）。若，则在matlab中的表达式为：r1=a1+b1*i; r2=a2+b2i; r=r1+r2; l 矢量的乘积和商：两矢量和的乘积是将第一个矢量伸长倍，再转动角度而得到的。若，；则，； 注:作为分母项的矢量不能为零。复数的微分：matlab中的相关计算函数：若r1=3*exp(4*i)， r2=4*exp(0.2i)；两矢量的乘除：rm=r1.*r2；rd=r1./r2；注：两矢量的点乘、点除为两矢量对应项的乘、除。如果已知矢量a，则可以通过下列运算求出复数的相关量，复数的模： abs(a)；复数的幅角： angle(a);复数a的实部：real(a): 复数的虚部： imag(a)；复数a的共轭复数：conj(rd)；2.2刚体和级杆组的运动分析由机构组成原理可知，任何平面连杆机构都可以看作是由若干个基本杆组依次连接于原动件和机架上而构成。因此对常见的基本杆组进行运动分析并建立相应的子函数库，在进行机构运动分析时，就可以根据机构的组成情况，将机构分解为原动件、机构和若干个基本杆组，编制一个调用该机构杆组子函数的主程序，即可实现对整个机构的运动分析。考虑到工程实际中大多数机构是级机构，本节主要介绍同一构件上两点间的运动分析和最常见的rrr型、rrp型和rpr型级杆组的运动分析的子函数，并通过一个典型的实例说明用该方法进行机构运动分析的方法和步骤。2.2.1 同一构件上的点的运动分析已知构件上一点的位置矢量p1,速度v1、加速度a1, p1p2两点间距离r 和位置角q(以x轴为起始线逆时针方向为正)，角速度w、角加速度。求构件上另一点n2点的位置、速度和加速度。图2-2由图可得位置的复数矢量方程为 (2-1)式(2-1)对时间t求导,即得n2点速度的复数矢量方程 (2-2) 同样,式(2-2)对时间t求导可得p2点的加速度复数矢量方程 (2-3)根据以上矢量方程编制的求p2点位置、速度、加速度的函数为：function a2,v2,p2=rigid(a1,v1,p1,w,e1,r,theta)p2=p1+r*exp(i*theta);v2=v1+i*r*w.*exp(i*theta); a2=a1-(w.2-i*e1).*r.*exp(i*theta);其中：p,v,a是以复数形式给出的相应点的位置、速度和加速度矢量。qtheta，角度；ww角速度；e1e，角加速度当p1为固定铰链点时，求解时只需使v1=0，a1=0即可。这样就可用来求主动件曲柄上某一点的运动参数。图2-32.2.2 rrr基本杆组运动分析rrr级基本杆组由两个构件和三个转动副构成，如图2-3所示。已知：两构件杆长分别为rl和r2、，两外接运动副p1、p2的运动矢量参数(p1,v1,a1和p2、v2、a2)。求内接转动副p3点位置、速度、加速度及两构件的角速度、角加速度。1. 位置分析间的矢量为 (2-4)求出矢量d的模d和幅角d=abs(d), =angle(d) (2-5)根据它们之间的几何关系，若或,则rrr杆组不能装配。在程序运行中应进行判断。与之间的夹角为 (2-6)而 (2-7)式（2-7）前的正负号决定于杆组起始位置的安装形式。如图2-3中实线方式装配(构成顺时针方向)时，取正号；若按图中虚线方式装配(构成逆时针方向)时，则取负号.节点的位置为 (2-8) (2-9)2 速度分析式(2-8)对时间求导得的速度矢量方程为 (2-10)上两式相减得: (2-11)令： ， 则式(2-11)可写成： (2-12)求解得 (2-13) (2-14)则： (2-15)3加速度分析式(2-10)对时间t求导后,可得点的加速度复数矢量方程: (2-16)上两式相减得 (2-17)令： ， 则上式可表示为 (2-18)由此解得 (2-19) (2-20)则p3点的加速度 (2-21)根据以上运动分析的过程编制成matlab函数krrr如下 function p3,the1,the2,v3,w1,w2,a3,e1,e2=krrr(p1,p2,v1,v2,a1,a2,r1,r2,m)% detremine the acceleration of rrr component, known p1,p2,r1,r2;p12=p2-p1;d=abs(p12);phi=angle(p12);test1=(r1+r2)-d;test2=d-abs(r1-r2);if(test1=0|test2=0) cosin=(r12+d.2-r22)./(2*r1*d); alph=atan2(sqrt(1-cosin.2),cosin); the1=phi+m*alph; p3=p1+r1*exp(i*the1); the2=angle(p3-p2);else disp(rrr can not be assembled!);end % velocity v21m=abs(v2-v1); vthe=angle(v2-v1);w2=v21m.*cos(the1-vthe)./(r2*sin(the2-the1);w1=-v21m.*cos(the2-vthe)./(r1*sin(the1-the2);v3=v1+i*w1.*(p3-p1);% acceleration analysisa21m=abs(a2-a1); theta=angle(a2-a1);e1=-(a21m.*cos(the2-theta)+r1*w1.2.*cos(the2-the1)-r2*w2.2)./(r1*sin(the1-the2);e2=(a21m.*cos(the1-theta)-r2*w2.2.*cos(the1-the2)+r1*w1.2)./(r2*sin(the2-the1);a3=a1+(i*e1-w1.2).*(p3-p1);变量说明：m装配模式，当p1、p2、p3点的顺序为逆时针方向时，即在图中实线所示位置时，m=1;反之则m=-1。the1,the2构件杆长矢量、的幅角。2.2.3 rrp级基本杆组运动分析rrp级基本杆组由两个构件和两个转动副一个转动副构成，如图2-4所示。已知构件杆长为rl，两外接运动副p1、p2的运动矢量参数(p1,v1,a1和p2、v2、a2)及导杆的转动参数（，,）。求内接转动副p3点位置p3、速度v3、加速度a3和构件1的角速度、角加速度以及滑块相对于导杆的相对运动参数(、)。图2-4 图2-5 1位置分析由图2-4可得矢量d (2-22)则d的模与幅角为 , (2-23)因此，由图2-4有 等式两端同乘，得 （2-24）取上式中复数的实部和虚部，并整理得（2-24a） (2-24b) (2-24c)式中根号前的正负号决定于杆组的初始装配形式。如按图中p3实线位置()装配时,取正号；若为虚线位置()装配,则取负号.上述转动导杆级组的两种装配形式，实际上只是时,以为圆心, 为半径作圆与导杆相交得的两种可能位置。而若时，则转动导杆级组如图（2-5)所示，这时点仅有一个可能位置，而相应公式中根号前应取正号。若，则转动导杆级组不能装配。节点的位置为 (2-25) （2-26）2) 速度分析由位置分析知以上二式分别对时间求导数，得p3点的速度矢量为 （227） （228）以上两式相减，并令: ，则式(2-28)可写成： (2-29)求解得 (2-30) (2-31) 将式（2-30）和式（2-31）代入式（2-27）或（2-28）即可求得v33)加速度分析式(2-27)、（2-28）对时间t求导，得点得加速度矢量方程为 （2-32a） (2-32b)将二式相减，并令： ， 展开整理后求解得（2-33a） (2-33b)以上运动分析的过程编制成matlab函数krrp， 源程序：function p3,the1,r2,v3,w1,vr2,a3,e1,ar2=krrp(p1,p2,v1,v2,vbet,a1,a2,abet,r1,bet,m)% 求解rrp杆组的运动分析程序，已知：p1,p2,r1,bet，v1,v2,a1,a2,vbet,abet,m;p21=p2-p1;d=abs(p21);phi=angle(p21);test1=r12-d.2.*(sin(phi-bet).2;if(test1=0) sq=sqrt(test1); mode=m; if(r1=d) mode=1; end r2=mode*sq-d.*cos(phi-bet); p3=p2+r2.*exp(i*bet); the1=angle(p3-p1);else disp(rrp can not be assembled!);end % 求解速度volicityv21m=abs(v2-v1);thev=angle(v2-v1);w1=(v21m.*sin(thev-bet)+r2.*vbet)./(r1.*cos(the1-bet);vr2=-(v21m.*cos(thev-the1)-r2.*vbet.*sin(bet-the1)./(cos(bet-the1);v3=v1+i*w1.*(p3-p1);%求解加速度a21m=abs(a2-a1);thea=angle(a2-a1);e1=(a21m.*sin(thea-bet)+r1*w1.2.*sin(the1-bet)+r2.*abet+2*vbet.*vr2)./(r1*cos(the1-bet);ar2=-a21m.*cos(thea-bet)+r2.*vbet.2-r1*e1.*sin(the1-bet)-r1*w1.2.*cos(the1-bet);a3=a1+(i*e1-w1.2).*r1.*exp(i*(the1);变量说明：m装配模式。当p1p3与p2p3间的角度为锐角时，m=1；否则m=1。bet、vbet、abet导杆的位置角、转动角速度和角加速度。其余变量的含义同前。 ,2.2.4 rpr级基本杆组运动分析rpr级基本杆组两外接副为转动副，内接副为移动副，如图26所示。为了便于对与该基本杆组相连的其他基本杆组进行运动分析，在该基本杆组导杆上增设了p3点。已知两外接运动副p1、p2的运动参数(p1,v1,a1和p2、v2、a2)和构件尺寸c和r3；求: 导杆的位置角、角速度、角加速度及滑块在导杆上的相对位置r2、相对速度vr2，相对加速度ar2和p3点的位置p3，速度v3和加速度a3.。1)位置分析由图2-6可得 (2-34) 则d的模与幅角为 （2-35） (2-36)图2-6若,则摆动滑块级杆组不能装配。 由图2-6可得 (2-37) (2-38) (2-39)式(2-39)中前得正负号决定于杆组的初始装配位置，如图（2-6）所示。若按实线位置装配，取正号；若按虚线位置装配，则取负号。 则点的位置矢量为 (2-40) 2)速度分析由于的位置矢量可表示为（2-41）上式对时间t 求一次导数，得的速度矢量方程为 (2-42 )令 ， 则式(2-42 )可写成： （2-43）上式两边同乘，并展开成分量式，经整理后得 （2-44） (2-45)式(2-40)对时间求一次导数，得点的速度矢量为 (2-46)将求出的各量代入上式，求得p3点的速度。3)加速度分析式(2-42)对时间t求一次导数,得点的加速度矢量方程 (2-47)令： ， 代入式（247），并展开成分量形，整理后得 (2-48a) (2-48b)式(2-46)对时间求t求一次导数,得点得加速度矢量为（2-49）以上运动分析过程编写出matlab函数：function p3,the,r2,v3,w,vr2,a3,e,ar2=krpr(p1,p2,v1,v2,a1,a2,r3,c,m)% detremine the acceleration of rrp（摆动滑块杆组） component, known p1,p2,r1,bet;% determine the position of rpr component, known p1,p2,r3,er;p21=p2-p1;d=abs(p21);r2=sqrt(d.2-c.2);alf=atan(c./r2);phi=angle(p21);dtest=max(d);if(dtest=c) the=phi+m*alf; p3=p1+(r3-i*c).*exp(i*the);else disp(krpr linkage can not be assembled!);end v21m=abs(v2-v1); thev=angle(v2-v1);dm=abs(p2-p1); phi=angle(p2-p1);w=v21m.*sin(thev-the)./r2;vr2=v21m.*cos(thev-the)-c*w;v3=v1+w.*(c+i*r3).*exp(i*the);% acceleration a21m=abs(a2-a1); thea=angle(a2-a1);dm=abs(p2-p1); phi=angle(p2-p1);e=(a21m.*sin(thea-the)-2*vr2.*w-c*w.2)./r2;ar2=a21m.*cos(thea-the)+w.2.*r2-c.*e;a3=a1+e.*(c+i*r3).*exp(i*the)+w.2.*(i*c-r3).*exp(i*the);变量说明：m装配模式。当p3在p1p2方向的右侧时，m=1；否则m=1。r2、vr2、ar2导杆的相对位置、速度和加速度。其余变量的含义同前。 ,2.2.5 连杆机构运动分析实例 1机构运动分析步骤机械运动分析一般按下述步骤进行：（1） 画出机构运动简图，并把各构件和各有关节点编号或标注字母，以方便编程。（2） 把机构分解为主动件、机架的基本杆组。（3） 根据已知条件，先对主动件进行分析，求出其与其它构件连接点的运动参数。图2-7（4） 从与主动件连接的构件开始，根据其它已知点的运动参数，依次对各杆组进行分析，求出有关的运动参数；若杆组还有其它待求点（质心、与其它构件的连接点），应用求解刚体运动分析函数求解。（5） 根据分析过程，编制运动分析主程序，调试程序并输出相应的运动参数和及其变化曲线。2实例分析例2.1图2-7所示机构为牛头刨床的机构运动简图。给定尺寸=120，l4315，h616；（长度单位为mm）。主动曲柄以顺时针方向转动，其角速度=-3.14159rads(常数)。试求滑块(刨刀)的位移s、速度v5和加速度a5。解：图示机构可分解为图2-8所示的原动件1，rpr转动导杆组2、3，和一个移动导杆组4、5。依次调用刚体运动分析函数rigid，rpr杆组运动分析函数krpr，就可求出b点的运动参数，它的位置、速度和加速度在水平方向的分量就是构件5相应的运动参数。图2-8clear r1=0.12;h=0.616; p1=0.315i; v1=0;a1=0; w1=-3.14159;e1=0; p3=0;v3=0;a3=0;t=0:360; phi1=t*1*pi/180; th1=pi/2-phi1;a2,v2,p2=rigid(a1,v1,p1,w1,e1,r1,th1);m=1;er=0;p4,th3,r2,v4,w3,vr2,a4,ep3,ar2=krpr(p3,p2,v3,v2,a3,a2,h,er,m);p5=real(p4);v5=real(v4); a5=real(a4); phid=t*1; %将弧度转化为度；figure(1); plot(phid,p5);title(刨刀位移曲线)xlabel(phi(度);ylabel(s(m)figure(2); plot(phid,v5)title(刨刀速度曲线); xlabel(phi(度)ylabel(v(m/s)figure(3); plot(phid,a5)title(刨刀加速度曲线)xlabel(phi(度); ylabel(a(m/s/s)根据以上程序作出的刨刀位移、速度、加速度图见下图图29刨刀运动线图2.3连杆机构动态静力分析用复数矢量法对机构进行动态分析也是以杆组为基础的，以下是对平面连杆机构常用的三种级杆组进行动态静力分析过程。2.3.1 rrr级杆组rrr级杆组的受力分析见图（2-10,a为此杆组的受力图，j1、j2分别为组成杆组的两构件标号；p1、p2、p3为三个运动参考点，位置矢量分别为；s1、s2分别为两构件的质心，其位置矢量为。代表构件j的角加速度，质心的加速度矢量用表示，及分别代表构件j的质量及转动惯量。已知，质心s1、s2处分别受重力 、；外力 、；两构件上所受到的外力偶矩分别为、。求解节点p1，p2上受外接运动副反力 、及杆组内部节点p3受运动副反力。运动副反力的符号规定为：某一构件与下一个构件相连接的运动副处，其反力符号为“+” ；而与前一个构件连接处，运动副反力符号为“”。力矩一律以逆时针为正。(图2-10)由构件的质心加速度和质量，分别求出两构件质心所受惯性力为： （2-50a） (2-50b)由构件的转动惯量和角加速度，分别求出两构件所受惯性力矩： （2-51a） (2-51b)作用在构件质心处的惯性力、已知外力和重力的合力为：（2-52a） (2-52b)设：p1s1两点间的矢量为；p2s2两点间的矢量为；p1p3间的矢量为；p2p3间的矢量为；根据力系平衡方程，构件上的力对任一点的合力矩为零，可得： (2-53a) (2-53b)上式是关于、的线性方程组，解之，得 (2-54a) (2-54b)式中： ， 根据力的平衡方程，作用在构件1上的合力为零，得平衡方程 (2-55)则p1运动副上的反力为：（2-56a）同理，可求得p2点的运动副反力为： (2-56b)rrr杆组的动态静力分析子函数为frrr，源程序：function f1,f2,f3=frrr(p1,p2,p3,ps1,ps2,m1,m2,js1,js2,f01,f02,m01,m02,as1,as2,e1,e2)% 求解rrr二级杆组受力分析的子函数 % 已知外接节点p1,p2,p3的位置.构件1,2的质心位置ps1,ps2% 质量m1,m2(kg)，对质心的转动惯量js1,js2. 构件1,2在质心处所受的外力f01,f02,m01,m02.% 质心加速度as1,as2、角加速度e1,e2 fg1=-as1.*m1; fg2=-as2.*m2; mg1=-js1*e1; mg2=-js2*e2; fs1=fg1+f01-m1*9.8i; fs2=fg2+f02-m2*9.8i;ms1=m01+mg1; ms2=m02+mg2; p1s=ps1-p1; ma1=real(p1s).*imag(fs1)-imag(p1s).*real(fs1)+ms1; p2s=ps2-p2; mb1=real(p2s).*imag(fs2)-imag(p2s).*real(fs2)+ms2; p13=p3-p1; p23=p3-p2; p13x=real(p13); p13y=imag(p13); p23x=real(p23); p23y=imag(p23); detc=p13y.*p23x-p23y.*p13x; f3x=(ma1.*p23x+mb1.*p13x)./detc; f3y=(ma1.*p23y+mb1.*p13y)./detc;f3=f3x+f3y*i; f1=f3+fs1; f2=-fs2+f3; 2.3.2 rrp级杆组图2-11（a）为一转动导杆级杆组。分别取两构件j1、j2作为分离体画出受力图，如图2-11b所示。j1、j2分别为组成杆组的两构件标号；p1、p2、p3为三个运动参考点，其位置矢量为、，s1、s2为两构件质心，其位置矢量为、。代表构件j的角加速度，质心的加速度矢量用表示，及分别代表构件j的质量(kg)及转动惯量。已知两构件质心s1、s2分别受重力、；构件受外部载荷 、，外力偶矩、的作用。所要求的节点p1，p3处受运动副反力分别为-、，在外接移动副上c2处构件j2受到的反力为，其方向为p2p3绕c2点逆时针转动90，力的作用点c2距结点p3沿导杆方向的距离为。 (图2-11)惯性力、惯性力矩、作用于质心处的合力以及两构件上的合力矩的计算同式（2-50a）(2-52b)。设：p1s1矢量为；p3s2矢量为；p1p3矢量为；p2p3矢量为；上述矢量的分量在下标上加注x、y。在构件j2上，根据力矩平衡条件，构件j2上载荷p3点的力矩： (2-57)式中：，在整个杆组上，根据力矩平衡条件,得 (2-58)式中， ， ，由式(2-57)(2-58)，得 (2-59)解方程(2-59),得 (2-60)代入式(2-60),得 (2-61)则 (2-62)根据力的平衡方程，由件j2上力的平衡方程得：(2-63)对构件j1上列平衡方程：得， (2-64)根据以上过程编制出转动导杆级组动态静力分析子程序frrp。function f1,f3,fc2,rc2=frrp(p1,p2,p3,ps1,ps2,r2,m1,m2,js1,js2,f01,f02,m01,m02,as1,as2,e1,e2)% 求解frpr类级杆组(转动导杆组rrp）受力分析的子函数 % 已知外接节点p1,p2,p3的位置p1,p2，p3,构件1,2的质心位置ps1,ps2，构件位置r2,% 质量m1,m2(kg)，对质心的转动惯量js1,js2. 构件1,2在质心处所受的外力f01,f02,m01,m02.% 质心加速度as1,as2、角加速度e1,e2 fg1=-as1.*m1; fg2=-as2.*m2; % 惯性力 mg1=-js1*e1; mg2=-js2*e2; % 惯性力矩 fs1=fg1+f01-m1*9.8i; fs2=fg2+f02-m2*9.8i; % 作用在质心上的合力 ms1=m01+mg1; ms2=m02+mg2; % 作用在质心上的合力偶矩 p1s=ps1-p1; p3s=ps1-p3; p13=p3-p1; m1s1=real(p1s).*imag(fs1)-imag(p1s).*real(fs1); % 构件j1上fs1对p1之力矩 m3s2=real(p3s).*imag(fs2)-imag(p3s).*real(fs2); % 构件fs2对结点p3之力矩 m1s2=real(p1s2).*imag(fs2)-imag(p1s2).*real(fs2); % fs2对p1点之矩 bet=angle(p3-p2); % 矢量p3-p2的幅角 det=imag(p13).*sin(bet)+real(p13).*cos(bet); fc2m=-(ms1+m1s1+m1s2-m3s2)./det; % 构件2受到的外来作用力的大小 rc2=-(ms2+m3s2)./fc3m; % 构件2受到的外来作用力的作用点距c2点距离 fc2=fc2m.*exp(i*(pi/2+bet); % 构件2受到的外接副作用力的矢量 f3=fs2+fc2; f1=f3+fs1; 2.3.3 rpr级组的动态静力分析图2-12（a）为摆动滑块级杆组的受力图.j1、j2分别为组成杆组的两构件标号；p1、p2、p3为三个运动参考点，其位置矢量分别为、。s1、s2为两构件质心，其位置矢量分别为、。已知,构件的角加速度，质心的加速度矢量，构件j的质量(kg)及转动惯量、作用在节点p3点的作用力f3。构件在质心s1、s2分别受重力 、，作用在构件上的外力 、外力偶矩、。求节点p1，p2上受运动副反力 、。构件j1和构件j2间的作用力的大小和作用点c2的位置。两构件内移动副间的相互作用力的方向为移动导路p2p3方向绕c2点逆时针方向转动90，c2点沿移动方向距p1点的距离为。为了使图形清晰，图中未画出构件的重力、惯性力、角加速度和加速度。(图2-12)图2-12（b）、（c）为组成此杆组两构件的受力图。惯性力、惯性力矩的计算同式（2-50a）(2-51b)。作用在构件质心处的惯性力、已知外力和重力和节点p3处的作用力的合力：（2-65a） (2-65b)作用于质心处的合力以及两构件上的合力矩设：p1s1两点间的矢量为；p2s2间的矢量为；p1p3间的矢量为；p2p3间的矢量为；在以下的计算中矢量的分量在下标上加注、。 根据力系平衡条件，构件上的力对构件上一点的力矩和为零，在构件1上的力矩平衡方程，(2-66a)构件2上对节点p2之合力矩 (2-66b)式(2-66a)-(2-66b)，得： (2-67)求解得 （2-68）式中： （2-69a） (2-69b)式（2-68）代入式(2-66a)，得 (2-70)则： (2-71)由构件j1的力平衡方程： 得p1点的运动副上的反力为： (2-72a)同理，可得p2点的运动副反力为： (2-72a)据此可编制出摆动滑块杆组的动态静力分析子程序frpr,源程序如下：function f1,f2,fc2,rc=frpr(p1,p2,p3,ps1,ps2,c,r2,m1,m2,js1,js2,f3,f01,f02,m01,m02,as1,as2,e1,e2)% 求解rpr二级杆组(摆动滑块组）受力分析的子函数% 已知外接节点p1,p2，p3的位置p1,p2，p3,构件1,2的质心位置ps1,ps2，构件尺寸c,位置尺寸r2,% 两构件的质量m1,m2(kg)，对质心的转动惯量js1,js2. % 构件1,2在质心处所受的外力f01,f02,m01,m02.节点p3处的作用力f3% 质心加速度as1,as2、角加速度e1,e2 fg1=-as1.*m1; fg2=-as2.*m2; %惯性力 mg1=-js1*e1; mg2=-js2*e2; %惯性力矩 fs1=fg1+f01-m1*9.8i+f3; fs2=fg2+f02-m2*9.8i; %质心处的合力 ms1=m01+mg1; ms2=m02+mg2; % 合力矩 p1s=ps1-p1; p13=p3-p1; m1=real(p1s).*imag(fs1)-imag(p1s).*real(fs1)+ms1+real(p13).*imag(f3)-imag(p13).*real(f3); % 构件1上的力对点p1之合力矩 p2s=ps2-p2; m2=real(p2s).*imag(fs2)-imag(p2s).*real(fs2)+ms2; %构件2上对点p2之合力矩 fc2m=-(m1+m2)./r2; % 构件1,2间的相互作用力的大小 rc=-m1/fc2m; % 构件1,2作用力的作用点距p1点距离 theta=angle(p3-p2); % 导路移动方向的幅角 fc2=fc2m.*exp(i*(pi/2+theta); % 构件1,2间的相互作用力的矢量 f1=fc2+fs1; f2=-fs2+fc2;2.3.4 作用有待求平衡力的构件的力分析图2-13以作用有待求平衡力偶矩的刚体为例，图2-13为一刚体的受力图，p1、p2为二个运动参考点，其位置矢量为、；s1为构件质心，其位置矢量为。代表构件的角加速度，质心的加速度矢量用表示，及分别代表构件的质量(kg)及转动惯量。已知，质心s1处分别受重力、外力和外力偶矩分别为，参考点p2处的作用力，求解运动副p1的反力及构件上的平衡力偶。由构件的质心加速度和质量，求出构件质心所受惯性力为： （2-73）由构件的转动惯量和角加速度，求出构件所受惯性力矩： （2-74）作用在构件质心处的惯性力、已知外力和重力的合力：（2-75）设：p1s1两点间的矢量为；p1p2间的矢量为；根据力系平衡方程，构件上对任一点的合力矩为零，可得： (2-76)上式求解，得（2-77）根据力的平衡方程，作用在构件1上的合力为零，得平衡方程 (2-78)则p1运动副上的反力为：（2-79）根据以上过程编制出求待求平衡力的构件的力分析子程序fbody。function f1,m1=fbody(p1,p2,ps1,m1,js1,f2,f01,m01,as1,e1)% 求解待求平衡力的构件的力分析的子函数 % 已知节点p1,p2的位置p1,p2.构件的质心位置ps1% 质量m1(kg)，对质心的转动惯量js1. 构件在质心处所受的外力f01,m01.% 质心加速度as1、角加速度e1 fg1=-as1.*m1; mg1=-js1*e1; fs1=fg1+f01-m1*9.8i; ms1=m01+mg1; p1s=ps1-p1;p12=p2-p1; m1=-(real(p12).*imag(f2)-imag(p12).*real(f2)+real(p1s).*imag(fs1)-imag(p1s).*real(fs1)+ms1); f1=f2+fs1;2.3.5 平面级机构的动态静力分析实例1分析步骤平面连杆机构动态静图片分析的步骤如下：（1）