基于MATLAB的六桿機構(gòu)動力學(xué)分析和仿真

...wd......wd......wd...六桿機構(gòu)的動力學(xué)分析仿真一系統(tǒng)模型建設(shè)為了對機構(gòu)進展仿真分析，首先必須建設(shè)機構(gòu)數(shù)學(xué)模型，即位置方程，然后利用MATLAB仿真分析工具箱Simulink對其進展仿真分析。圖3．24所示是由原動件(曲柄1)和RRR—RRP六桿機構(gòu)。各構(gòu)件的尺寸為r1=400mm，r2=1200mm，r3=800mm,r4=1500mm，r5=1200mm；各構(gòu)件的質(zhì)心為rc1=200mm,rc2＝600mm,rc3=400mm，rc5=600mm；質(zhì)量為m1=1．2kg，m2＝3kg，m3＝2．2kg；m5=3．6kg，m6=6kg;轉(zhuǎn)動慣量為J1=0.016kg·m2，J2=0.25kg·m2；J3=0.09kg·m2，J5=0.45kg·m2；構(gòu)件6的工作阻力F6=1000N，其他構(gòu)件所受外力和外力矩均為零，構(gòu)件1以等角速度10rad/s逆時針方向回轉(zhuǎn)，試求不計摩擦?xí)r，轉(zhuǎn)動副A的約束反力、驅(qū)動力矩、移動副F的約束反力。圖1-1此機構(gòu)模型可以分為曲柄的動力學(xué)、RRRII級桿組的動力學(xué)和RRPII級桿組的動力學(xué)，再分別對這三個模型進展相應(yīng)參數(shù)的求解。圖1-2AB構(gòu)件受力模型如上圖1-2對于曲柄AB由理論力學(xué)可以列出表達式：由運動學(xué)知識可以推得：將上述各式合并成矩陣形式有，(1-21)如圖1-3，對構(gòu)件BC的約束反力推導(dǎo)如下，圖1-3BC構(gòu)件受力模型如圖1-4，對構(gòu)件BC的約束反力推導(dǎo)如下，圖1-4CD構(gòu)件受力模型由運動學(xué)可以推導(dǎo)得，將上述BC構(gòu)件，CD構(gòu)件各式合并成矩陣形式有，=(1-22)如圖1-5對構(gòu)件5進展約束反力的推導(dǎo)如下，圖1-5CE桿件受力模型如圖1-6對滑塊進展受力分析如下，滑塊受力模型由運動學(xué)可推，〔1-23〕二編程與仿真利用MATLAB進展仿真分析，主要包括兩個步驟：首先是編制計算所需要的函數(shù)模塊，然后利用其仿真工具箱Simulink建設(shè)仿真系統(tǒng)框圖，設(shè)定初始參數(shù)進展仿真分析。針對建設(shè)完成的數(shù)學(xué)模型，為了進展矩陣運算，根據(jù)以上式子編制M函數(shù)文件chengcrank.m，chengrrr.m、chengcrankdy.m、chengrrrdy.m、chengrrp.m和chengrrpdy.m如下：曲柄原動件M函數(shù)文件chengcrank.m：functiony=chengcrank(x)%%Functiontocomputetheacclerationofcrank%Inputparameters%x(1)=theta-1%x(2)=dtheta-1%x(3)=ddtheta-1%0utputparameters%y(1)=Re[ddB]%y(2)=Im[ddB]r1=0.4;ddB=[r1*x(3)*cos(x(1)+pi/2)+r1*x(2)^2*cos(x(1)+pi);r1*x(3)*sin(x(1)+pi/2)+r1*x(2)^2*sin(x(1)+pi)];y=ddB;RRRII級桿組M函數(shù)文件chengrrr.m：functiony=chengrrr(x)%functiontocomputetheaccelerationforRRRbargroup%Inputparameters%x(1)=theta-2%x(2)=theta-3%x(3)=dtheta-2%x(4)=dtheta-3%x(5)=Re[ddB]%x(6)=Im[ddB]%Outputparameters%y(1)=ddtheta-2%y(2)=ddtheta-3%y(3)=Re[ddC]%y(4)=Im[ddC]r2=1.2;r3=0.8;ReddD=0;ImddD=0;a=[r2*cos(x(1)+pi/2)-r3*cos(x(2)+pi/2);r2*sin(x(1)+pi/2)-r3*sin(x(2)+pi/2)];b=[-r2*cos(x(1)+pi)r3*cos(x(2)+pi);-r2*sin(x(1)+pi)r3*sin(x(2)+pi)]*[x(3)^2;x(4)^2]+[ReddD-x(5);ImddD-x(6)];ddth=inv(a)*b;y(1)=ddth(1);y(2)=ddth(2);y(3)=x(5)+r2*ddth(1)*cos(x(1)+pi/2)+r2*x(3)^2*cos(x(1)+pi);y(4)=x(6)+r2*ddth(1)*sin(x(1)+pi/2)+r2*x(3)^2*sin(x(1)+pi);曲柄原動件動力學(xué)M函數(shù)文件chengcrankdy.m：functiony=chengcrankdy(x)%FunctionforDyanmicanalysisofcrank%%Inputparameters%x(1)=theta-1%x(2)=dtheta-1%x(3)=ddtheta-1%x(4)=RxB%x(5)=RyB%%0utputparameters%y(1)=RxA%y(2)=RyA%y(3)=M1g=9.8;%重力加速度r1=0.4;%曲柄長度rc1=0.2;%質(zhì)心離鉸鏈A的距離m1=1.2;%曲柄質(zhì)量J1=0.016;%繞質(zhì)心轉(zhuǎn)動慣量Fx1=0;Fy1=0;MF=0;%作用于質(zhì)心的外力和外力矩ReddA=0;ImddA=0;%鉸鏈A的加速度y(1)=m1*ReddA+m1*rc1*x(3)*cos(x(1)+pi/2)+m1*rc1*x(2)^2*cos(x(1)+pi)-Fx1+x(4);y(2)=m1*ImddA+m1*rc1*x(3)*sin(x(1)+pi/2)+m1*rc1*x(2)^2*sin(x(1)+pi)-Fy1+x(5)+m1*g;y(3)=J1*x(3)-y(1)*rc1*sin(x(1))+y(2)*rc1*cos(x(1))-x(4)*(r1-rc1)*sin(x(1))+x(5)*(r1-rc1)*cos(x(1))-MF;RRRII級桿組動力學(xué)M函數(shù)文件chengrrrdy.m：functiony=chengrrrdy(x)%FunctionforDyanmicanalysisofRRRdayardgroup%Inputparameters%x(1)=theta-2%x(2)=theta-3%x(3)=dtheta-2%x(4)=dtheta-3%x(5)=ddtheta-2%x(6)=ddtheta-3%x(7)=Re[ddB]%x(8)=Im[ddB]%x(9)=Fx3%x(10)=Fy3%x(11)=M3%0utputparameters%y(1)=RxB%Y(2)=RyB%y(3)=RxC%y(4)=RyC%y(5)=RxD%y(6)=RyDg=9.8;%重力加速度r2=1.2;r3=0.8;%兩桿的長度rc2=0.6;rc3=0.4;%質(zhì)心到鉸鏈B的距離%質(zhì)心到鉸鏈D的距離m2=3;m3=2.2;%兩桿的質(zhì)量J2=0.25;J3=0.09;%兩桿的轉(zhuǎn)動慣量ReddD=0;ImddD=0;Fx2=0;Fy2=0;M2=0;%2桿的外力和外力矩a=zeros(6);a(1,1)=1;a(1,3)=1;a(2,2)=1;a(2,4)=1;a(3,1)=rc2*sin(x(1));a(3,2)=-rc2*cos(x(1));a(3,3)=-(r2-rc2)*sin(x(1));a(3,4)=(r2-rc2)*cos(x(1));a(4,3)=-1;a(4,5)=1;a(5,4)=-1;a(5,6)=1;a(6,3)=(r3-rc3)*sin(x(2));a(6,4)=-(r3-rc3)*cos(x(2));a(6,5)=rc3*sin(x(2));a(6,6)=-rc3*cos(x(2));b=zeros(6,1);b(1,1)=m2*rc2*x(5)*cos(x(1)+pi/2)+m2*x(7)+m2*rc2*x(3)^2*cos(x(1)+pi)-Fx2;b(2,1)=m2*rc2*x(5)*sin(x(1)+pi/2)+m2*x(8)+m2*rc2*x(3)^2*sin(x(1)+pi)-Fy2+m2*g;b(3,1)=J2*x(5)-M2;b(4,1)=m3*rc3*x(6)*cos(x(2)+pi/2)+m3*ReddD+m3*rc3*x(4)^2*cos(x(2)+pi)-x(9);b(5,1)=m3*rc3*x(6)*sin(x(2)+pi/2)+m3*ImddD+m3*rc3*x(4)^2*sin(x(2)+pi)-x(10)+m3*g;b(6,1)=J3*x(6)-x(11);y=inv(a)*b;RRPII級桿組M函數(shù)文件：functiony=chengrrp(x)%functiontocomputetheaccelerationforRRPbargroup%Inputparameters%x(1)=theta-5%x(2)=dtheta-5%x(3)=Re[ddC]%x(4)=Im[ddC]%x(5)=ds%Outputparameters%y(1)=ddtheta-5%y(2)=ddsr5=1.2;th6=0;ReddD=0;ImddD=0;a=[r5*cos(x(1)+pi/2)-cos(th6);r5*sin(x(1)+pi/2)-sin(th6)];b=[-r5*cos(x(1)+pi)0;-r5*sin(x(1)+pi)0]*[x(2)^2;x(5)]+[ReddD-x(3);ImddD-x(4)];y=inv(a)*b;RRPII級桿組動力學(xué)M函數(shù)文件：functiony=chengrrpdy(x)%FunctionforDyanm5canalysisofRRPdayardgroup%Inputparameters%x(1)=theta-5%x(2)=dtheta-5%x(3)=ddtheta-5%x(4)=dds-6%x(5)=Re[ddC]%x(6)=Im[ddC]%0utputparameters%y(1)=RxC%Y(2)=RyC%y(3)=RxE%y(4)=Ry%y(5)=RF%移動副的約束反力g=9.8;%重力加速度r5=1.2;%桿的長度rc5=0.6;%質(zhì)心到鉸鏈B的距離m5=3.6;m6=6;%桿、塊的質(zhì)量J5=0.45;Fx5=0;Fy5=0;Fx6=1000;Fy6=0;M5=0;th6=0;a=zeros(5);a(1,1)=1;a(1,3)=1;a(2,2)=1;a(2,4)=1;a(3,1)=rc5*sin(x(1));a(3,2)=-rc5*cos(x(1));a(3,3)=-(r5-rc5)*sin(x(1));a(3,4)=(r5-rc5)*cos(x(1));a(4,3)=-1;a(4,5)=-sin(th6);a(5,4)=-1;a(5,5)=cos(th6);b=zeros(5,1);b(1,1)=m5*x(5)+m5*rc5*x(3)*cos(x(1)+pi/2)+m5*rc5*x(2)^2*cos(x(1)+pi)-Fx5;b(2,1)=m5*x(6)+m5*rc5*x(3)*sin(x(1)+pi/2)+m5*rc5*x(2)^2*sin(x(1)+pi)-Fy5+m5*g;b(3,1)=J5*x(3)-M5;b(4,1)=m6*x(4)*cos(th6)-Fx6;b(5,1)=m6*x(4)*sin(th6)-Fx6+m6*g;y=inv(a)*b;三系統(tǒng)仿真框圖進入MATLAB，在命令欄中鍵入Simulink進入仿真界面，根據(jù)信息傳遞的邏輯關(guān)系，建設(shè)仿真系統(tǒng)框圖如圖3-1.然后設(shè)定各環(huán)節(jié)的初始參數(shù)，即可以對機構(gòu)進展運動學(xué)仿真分析，再利用MATLAB的plot命令根據(jù)需要繪制曲線。圖3-1四仿真的實現(xiàn)再設(shè)計完成仿真框圖之后，為了進展仿真還必須設(shè)定初始參數(shù)值。連桿機構(gòu)桿長已經(jīng)在simulink框圖中給定，如果設(shè)定初始夾角為62，=10rad/s，曲柄1作勻速轉(zhuǎn)動〔即〕，接下來要確定桿2，3的角位移和角速度，桿5的角位移和角速度，滑塊的速度?？梢岳眯疗丈椒ā苍贛ATLAB命令框中輸入M函數(shù)為rrrposi〕求得=0.3612rad/s，=1.8101rad/s，再