基于-matlab的六桿機構動力學分析及仿真

1、<p>　　六桿機構的動力學分析仿真</p><p><b>　　一 系統(tǒng)模型建立</b></p><p>　　為了對機構進行仿真分析，首先必須建立機構數(shù)學模型，即位置方程，然后利用MATLAB仿真分析工具箱Simulink對其進行仿真分析。圖3．24所示是由原動件(曲柄1)和RRR—RRP六桿機構。各構件的尺寸為r1=400mm，r2=1200mm，r3=

2、800mm,r4=1500mm，r5=1200mm；各構件的質心為rc1=200mm,rc2＝600mm,rc3=400mm，rc5=600mm；質量為m1=1．2kg，m2＝3kg，m3＝2．2kg；m5=3．6kg，m6=6kg; 轉動慣量為J1=0.016kg·m2，J2=0.25kg·m2；J3=0.09kg·m2，J5=0.45kg·m2；構件6的工作阻力F6=1000N，其他構件所受外

3、力和外力矩均為零，構件1以等角速度10 rad/s逆時針方向回轉，試求不計摩擦時，轉動副A的約束反力、驅動力矩、移動副F的約束反力。 </p><p><b>　　圖1-1</b></p><p>　　此機構模型可以分為曲柄的動力學、RRR II級桿組的動力學和RRP II級桿組的動力學，再分別對這三個模型進行相應參數(shù)的求解。</p><p>

4、　　圖1-2 AB構件受力模型</p><p>　　如上圖1-2對于曲柄AB由理論力學可以列出表達式：</p><p>　　由運動學知識可以推得：</p><p>　　將上述各式合并成矩陣形式有，</p><p><b>　　(1-21)</b></p><p>　　如圖1-3，對構件BC的約束

5、反力推導如下，</p><p>　　圖1-3 BC構件受力模型</p><p>　　如圖1-4，對構件BC的約束反力推導如下，</p><p>　　圖 1-4 CD構件受力模型</p><p>　　由運動學可以推導得，</p><p>　　將上述BC構件，CD構件各式合并成矩陣形式有，</p><

6、p><b>　　= (1-22)</b></p><p>　　如圖1-5 對構件5進行約束反力的推導如下，</p><p>　　圖1-5 CE桿件受力模型</p><p>　　如圖1-6 對滑塊進行受力分析如下，</p><p><b>　　滑塊受力模型</b></p><

7、;p><b>　　由運動學可推，</b></p><p><b>　?。?-23）</b></p><p><b>　　二 編程與仿真</b></p><p>　　利用MATLAB進行仿真分析，主要包括兩個步驟：首先是編制計算所需要的函數(shù)模塊，然后利用其仿真工具箱Simulink建立仿真系統(tǒng)框圖

8、，設定初始參數(shù)進行仿真分析。針對建立完成的數(shù)學模型，為了進行矩陣運算，根據(jù)以上式子編制M函數(shù)文件chengcrank.m ，chengrrr.m、chengcrankdy.m、chengrrrdy.m、chengrrp.m和chengrrpdy.m如下：</p><p>　　曲柄原動件M函數(shù)文件chengcrank.m：</p><p>　　function y=chengcrank(x)

9、</p><p>　　%%Function to compute the accleration of crank</p><p>　　%Input parameters</p><p>　　%x(1)=theta-1</p><p>　　%x(2)=dtheta-1</p><p>　　%x(3)=ddtheta-1

10、</p><p>　　%0utput parameters</p><p>　　%y(1)=Re[ddB]</p><p>　　%y(2)=Im[ddB]</p><p><b>　　r1=0.4;</b></p><p>　　ddB=[r1*x(3)*cos(x(1)+pi/2)+r1*x(2)

11、^2*cos(x(1)+pi);r1*x(3)*sin(x(1)+pi/2)+r1*x(2)^2*sin(x(1)+pi)];</p><p><b>　　y=ddB;</b></p><p>　　RRR II級桿組M函數(shù)文件chengrrr.m：</p><p>　　function y=chengrrr(x)</p><

12、p>　　%function to compute the acceleration for RRR bar group</p><p>　　%Input parameters</p><p>　　%x(1)=theta-2</p><p>　　%x(2)=theta-3</p><p>　　%x(3)=dtheta-2</p>

13、;<p>　　%x(4)=dtheta-3</p><p>　　%x(5)=Re[ddB]</p><p>　　%x(6)=Im[ddB]</p><p>　　%Output parameters</p><p>　　%y(1)=ddtheta-2</p><p>　　%y(2)=ddtheta-3<

14、;/p><p>　　%y(3)=Re[ddC]</p><p>　　%y(4)=Im[ddC]</p><p>　　r2=1.2; r3=0.8; ReddD=0; ImddD=0;</p><p>　　a=[r2*cos(x(1)+pi/2) -r3*cos(x(2)+pi/2); r2*sin(x(1)+pi/2) -r3*sin(x(2)+

15、pi/2)];</p><p>　　b=[-r2*cos(x(1)+pi) r3*cos(x(2)+pi); </p><p>　　-r2*sin(x(1)+pi) r3*sin(x(2)+pi)]*[x(3)^2;x(4)^2]+[ReddD-x(5);ImddD-x(6)];</p><p>　　ddth=inv(a)*b;</p><p

16、>　　y(1)=ddth(1);</p><p>　　y(2)=ddth(2);</p><p>　　y(3)=x(5)+r2*ddth(1)*cos(x(1)+pi/2)+r2*x(3)^2*cos(x(1)+pi);</p><p>　　y(4)=x(6)+r2*ddth(1)*sin(x(1)+pi/2)+r2*x(3)^2*sin(x(1)+pi);&

17、lt;/p><p>　　曲柄原動件動力學M函數(shù)文件 chengcrankdy.m：</p><p>　　function y=chengcrankdy(x)</p><p>　　%Function for Dyanmic analysis of crank</p><p>　　%%Input parameters</p><p

18、>　　%x(1)=theta-1</p><p>　　%x(2)=dtheta-1</p><p>　　%x(3)=ddtheta-1</p><p><b>　　%x(4)=RxB</b></p><p><b>　　%x(5)=RyB</b></p><p>　　%

19、%0utput parameters</p><p><b>　　%y(1)=RxA</b></p><p><b>　　%y(2)=RyA</b></p><p><b>　　%y(3)=M1 </b></p><p>　　g=9.8; %重力加速度</p>&

20、lt;p>　　r1=0.4; %曲柄長度</p><p>　　rc1=0.2;%質心離鉸鏈A的距離</p><p>　　m1=1.2;%曲柄質量</p><p>　　J1=0.016; %繞質心轉動慣量</p><p>　　Fx1=0; Fy1=0; MF=0;%作用于質心的外力和外力矩</p><p>　　Re

21、ddA=0; ImddA=0;%鉸鏈A的加速度</p><p>　　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);</p><p>　　y(2)=m1*ImddA+m1*rc1*x(3)*sin(x(1)+pi/2)+m1*rc1*x(2)^2*sin(x(1)+pi)-Fy1+x(

22、5)+m1*g;</p><p>　　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;</p><p>　　RRR II級桿組動力學M函數(shù)文件 chengrrrdy.m：</p><p>　　function y

23、=chengrrrdy(x)</p><p>　　%Function for Dyanmic analysis of RRR dayard group</p><p>　　%Input parameters</p><p>　　%x(1)=theta-2 </p><p>　　%x(2)=theta-3</p>&l

24、t;p>　　%x(3)=dtheta-2 </p><p>　　%x(4)=dtheta-3</p><p>　　%x(5)=ddtheta-2 </p><p>　　%x(6)=ddtheta-3</p><p>　　%x(7)=Re[ddB] </p><p>　　%x(8)=Im[ddB]&l

25、t;/p><p>　　%x(9)=Fx3 </p><p>　　%x(10)=Fy3 </p><p><b>　　%x(11)=M3</b></p><p>　　%0utput parameters</p><p><b>　　%y(1)=RxB<

26、/b></p><p><b>　　%Y(2)=RyB</b></p><p><b>　　%y(3)=RxC</b></p><p><b>　　%y(4)=RyC</b></p><p><b>　　%y(5)=RxD</b></p>

27、<p><b>　　%y(6)=RyD</b></p><p>　　g=9.8; %重力加速度</p><p>　　r2=1.2; r3=0.8; %兩桿的長度</p><p>　　rc2=0.6;rc3=0.4; %質心到鉸鏈B的距離 %質心到鉸鏈D的距離 </p><p>　　m2=3; m3=2.

28、2; %兩桿的質量</p><p>　　J2=0.25;J3=0.09;%兩桿的轉動慣量</p><p>　　ReddD=0;ImddD=0;</p><p>　　Fx2=0; Fy2=0; </p><p>　　M2=0; %2桿的外力和外力矩</p><p>　　a=zeros(6);</p><

29、;p><b>　　a(1,1)=1;</b></p><p>　　a(1,3)=1; </p><p>　　a(2,2)=1; </p><p><b>　　a(2,4)=1;</b></p><p>　　a(3,1)=rc2*sin(x(1)); </p><p>

30、;　　a(3,2)=-rc2*cos(x(1));</p><p>　　a(3,3)=-(r2-rc2)*sin(x(1));</p><p>　　a(3,4)=(r2-rc2)*cos(x(1));</p><p>　　a(4,3)=-1; </p><p><b>　　a(4,5)=1;</b></p>

31、<p>　　a(5,4)=-1; </p><p><b>　　a(5,6)=1;</b></p><p>　　a(6,3)=(r3-rc3)*sin(x(2));</p><p>　　a(6,4)=-(r3-rc3)*cos(x(2));</p><p>　　a(6,5)=rc3*sin(x(2));<

32、;/p><p>　　a(6,6)=-rc3*cos(x(2));</p><p>　　b=zeros(6,1);</p><p>　　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;</p><p>　　b(2,1)=m2*rc2*x(5)*sin(

33、x(1)+pi/2)+m2*x(8)+m2*rc2*x(3)^2*sin(x(1)+pi)-Fy2+m2*g;</p><p>　　b(3,1)=J2*x(5)-M2;</p><p>　　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);</p><p>　　b(5

34、,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;</p><p>　　b(6,1)=J3*x(6)-x(11); </p><p>　　y=inv(a)*b;</p><p>　　RRP II級桿組M函數(shù)文件：</p><p>　　fu

35、nction y=chengrrp(x)</p><p>　　%function to compute the acceleration for RRP bar group</p><p>　　%Input parameters</p><p>　　%x(1)=theta-5</p><p>　　%x(2)=dtheta-5</p>

36、;<p>　　%x(3)=Re[ddC]</p><p>　　%x(4)=Im[ddC]</p><p><b>　　%x(5)=ds</b></p><p>　　%Output parameters</p><p>　　%y(1)=ddtheta-5</p><p><b&g

37、t;　　%y(2)=dds</b></p><p>　　r5=1.2; th6=0; ReddD=0; ImddD=0;</p><p>　　a=[r5*cos(x(1)+pi/2) -cos(th6); r5*sin(x(1)+pi/2) -sin(th6)];</p><p>　　b=[-r5*cos(x(1)+pi) 0; -r5*sin(x(

38、1)+pi) 0]*[x(2)^2; x(5)]+[ReddD-x(3); ImddD-x(4)];</p><p>　　y=inv(a)*b;</p><p>　　RRP II級桿組動力學M函數(shù)文件：</p><p>　　function y=chengrrpdy(x)</p><p>　　%Function for Dyanm5c ana

39、lysis of RRP dayard group</p><p>　　%Input parameters</p><p>　　%x(1)=theta-5 </p><p>　　%x(2)=dtheta-5 </p><p>　　%x(3)=ddtheta-5</p><p>　　%x(4)=dds-6 <

40、;/p><p>　　%x(5)=Re[ddC] </p><p>　　%x(6)=Im[ddC]</p><p>　　%0utput parameters</p><p>　　%y(1)=RxC </p><p><b>　　%Y(2)=RyC</b></p><p>

41、　　%y(3)=RxE </p><p><b>　　%y(4)=RyE</b></p><p>　　%y(5)=RF %移動副的約束反力</p><p>　　g=9.8; %重力加速度</p><p>　　r5=1.2; %桿的長度</p><p>　　rc5=0.6; %質心到鉸鏈

42、B的距離 </p><p>　　m5=3.6; m6=6; %桿、塊的質量</p><p><b>　　J5=0.45;</b></p><p>　　Fx5=0; Fy5=0;</p><p>　　Fx6=1000; Fy6=0;</p><p><b>　　M5=0;</b&g

43、t;</p><p><b>　　th6=0;</b></p><p>　　a=zeros(5); </p><p>　　a(1,1)=1; </p><p>　　a(1,3)=1; </p><p><b>　　a(2,2)=1;</b></p><

44、;p><b>　　a(2,4)=1;</b></p><p>　　a(3,1)=rc5*sin(x(1));</p><p>　　a(3,2)=-rc5*cos(x(1));</p><p>　　a(3,3)=-(r5-rc5)*sin(x(1)); </p><p>　　a(3,4)=(r5-rc5)*cos(x

45、(1)); </p><p>　　a(4,3)=-1;</p><p>　　a(4,5)=-sin(th6); </p><p>　　a(5,4)=-1;</p><p>　　a(5,5)=cos(th6);</p><p>　　b=zeros(5,1);</p><p>　　b(1,1)

46、=m5*x(5)+m5*rc5*x(3)*cos(x(1)+pi/2)+m5*rc5*x(2)^2*cos(x(1)+pi)-Fx5;</p><p>　　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;</p><p>　　b(3,1)=J5*x(3)-M5; </p>

47、;<p>　　b(4,1)=m6*x(4)*cos(th6)-Fx6; </p><p>　　b(5,1)=m6*x(4)*sin(th6)-Fx6+m6*g;</p><p>　　y=inv(a)*b;</p><p><b>　　三 系統(tǒng)仿真框圖</b></p><p>　　進入MATLAB，在命令欄

48、中鍵入Simulink進入仿真界面，根據(jù)信息傳遞的邏輯關系，建立仿真系統(tǒng)框圖如圖3-1. 然后設定各環(huán)節(jié)的初始參數(shù)，即可以對機構進行運動學仿真分析，再利用MATLAB的plot命令根據(jù)需要繪制曲線。</p><p><b>　　圖3-1</b></p><p><b>　　四 仿真的實現(xiàn)</b></p><p>　　再設計

49、完成仿真框圖之后，為了進行仿真還必須設定初始參數(shù)值。連桿機構桿長已經(jīng)在simulink框圖中給定，如果設定初始夾角為62，=10 rad/s，曲柄1作勻速轉動（即），接下來要確定桿2，3的角位移和角速度，桿5的角位移和角速度，滑塊的速度。</p><p>　　可以利用辛普森方法（在MATLAB命令框中輸入M函數(shù)為rrrposi）求得=0.3612rad/s，=1.8101rad/s，再利用MATLAB(在命令框輸

50、入rrrvel)求出W2=-2.2345, W3=3.3250,再利用桿3的角位移和角速度、桿5的角位移求得（在MATLAB命令框中輸入M函數(shù)為compvel）W5=0.6962，ds=-3.1323。對仿真框圖中各積分器設定參數(shù)變量x并在matlab命令框輸入變量 x=[62*pi/180 10 0.3612 1.8101 -2.2345 3.3250 -41*pi/180 0.6962 -3.1323];其中初始數(shù)值分別對應：the

51、ta-1、omega-1、theta-2、omega-2、omega-3、theta-5,omega-5 ds，以及仿真時間為1s，后進行仿真，利用MATLAB中的plot繪圖命令把角速度曲線分別繪制出來。在MATLAB命令中鍵入：plot（tout,simout(:,1)），plot（tout,simout(:,2)），plot（tout,simout(:,3)) ，plot（tout,simout</p><p&