基于MATLAB的六杆机构动力学分析与仿真.docx
《基于MATLAB的六杆机构动力学分析与仿真.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的六杆机构动力学分析与仿真.docx（16页珍藏版）》请在冰豆网上搜索。
基于MATLAB的六杆机构动力学分析与仿真
六杆机构的动力学分析仿真
一 系统模型建立
为了对机构进行仿真分析，首先必须建立机构数学模型,即位置方程，然后利用MATLAB仿真分析工具箱Sｉmulink对其进行仿真分析。
图３.24所示是由原动件（曲柄1）和RRR—RＲP六杆机构。
各构件的尺寸为r1=４００mm，r2＝1200mｍ，r3=800mm,r4=1500mm,r５=1200mm;各构件的质心为rｃ1=2０0mm，rｃ2＝60０mｍ,rc3=400mm，rｃ５=600mm;质量为ｍ1＝1.2kｇ,ｍ2=3ｋg,m3=2．2ｋg;m５=3．6kg，m6＝6kg;转动惯量为Ｊ1=０.０１6kg·ｍ2，J２=0.25kg·ｍ2；Ｊ3=0.0９kg·ｍ２，Ｊ５=0.45kg·m2;构件6的工作阻力F６＝1000N,其他构件所受外力和外力矩均为零,构件1以等角速度１0rad/s逆时针方向回转,试求不计摩擦时，转动副Ａ的约束反力、驱动力矩、移动副Ｆ的约束反力。
图1－１
此机构模型可以分为曲柄的动力学、RRRII级杆组的动力学和RRＰ ＩＩ级杆组的动力学，再分别对这三个模型进行相应参数的求解。
图1-2ＡB构件受力模型
如上图1-2对于曲柄AB由理论力学可以列出表达式：
由运动学知识可以推得：
将上述各式合并成矩阵形式有，
（1－２1）
如图1-3,对构件BＣ的约束反力推导如下,
图1-３ BＣ构件受力模型
如图１－４,对构件BC的约束反力推导如下，
图1-4CD构件受力模型
由运动学可以推导得,
将上述ＢC构件,CＤ构件各式合并成矩阵形式有，
=
（1－22）
如图１－５对构件５进行约束反力的推导如下,
图1-5 ＣE杆件受力模型
如图1-６对滑块进行受力分析如下，
滑块受力模型
由运动学可推，
（1－2３）
二编程与仿真
利用MATLAB进行仿真分析,主要包括两个步骤：
首先是编制计算所需要的函数模块，然后利用其仿真工具箱Sｉmulｉnk建立仿真系统框图，设定初始参数进行仿真分析。
针对建立完成的数学模型,为了进行矩阵运算,根据以上式子编制M函数文件ｃhengcrank.m,chengrｒr.m、chengｃrａｎｋｄｙ．ｍ、chｅnｇｒrrdy．ｍ、chenｇrｒp.m和chengrｒpdy．ｍ如下:
曲柄原动件Ｍ函数文件chengcｒank.m：
fuｎctｉｏnｙ=cheｎgcrank（ｘ）
%％Fuｎctiontｏcｏmputetheaｃcleraｔionofcrank
％Inputpaｒaｍeterｓ
％x（１）=ｔheta-1
%ｘ
（2）=dtheta-1
%x（3）=dｄtheta-1
%0utput parａmeteｒs
％y
（1）=Re［ｄdB］
%y
（2）=Ｉm［ddＢ]
ｒ1=0．4；
ddB=[r1*x（3）＊cｏｓ（ｘ（１）+pｉ/2）+r1*x
（2）^2*coｓ（ｘ（１）＋pi）;r1*x（3）*ｓiｎ（x
（1）+ｐi/2）+r1*x
（2）^2*sｉｎ（x
（1）+pi）］;
ｙ＝dｄＢ;
RRRIＩ级杆组M函数文件chｅｎgｒｒr.m:
ｆuncｔｉony=cheｎｇrrr（x）
％ｆuｎｃtion to computethｅ aｃｃｅleratｉonforＲＲRbａｒ gｒoup
％Iｎpuｔpａrameters
%x
（1）＝ｔheｔa-2
％x
（2）=ｔhetａ-３
%x（3）=dthetａ-2
%ｘ（4）=dｔhｅta-3
%x（５）=Re［ddB]
%ｘ（6）＝Ｉm［ｄdB］
％Oｕtput pａｒameteｒs
%y（１）=ddtｈｅtａ-2
％ｙ
（2）=ddtheｔa－3
%y（３）=Re[ｄdC］
％y（4）＝Iｍ[ｄdC］
r2=1．2；r3=0.８;ReddＤ=0; ＩmddD＝0;
a＝［r２*ｃos（x
（1）+ｐｉ／2）-r３＊cｏs（x（２）＋pi/2）;r2*sin（x
（1）+pi/2）-r3*sin（x
（2）+pi/2）];
b=[-r２*coｓ（x
（1）＋pi）r３*cｏs（x
（2）+pｉ）；
－ｒ２*ｓｉｎ（x
（1）+ｐi）ｒ３*sin（x
（2）+ｐｉ）]*[x（3）^2;x（4）^2]+[RｅｄdD-ｘ（5）;IｍdｄD-ｘ（６）]；
ddtｈ=inｖ（a）＊b;
ｙ
（1）＝ｄdｔh（１）;
y
（2）＝ddth
（2）;
y（３）=ｘ（5）+ｒ2*ddtｈ
（1）*ｃｏs（x
（1）+pi/２）+r2*ｘ（3）^2*ｃos（x（１）+ｐi）;
y（４）=x（6）+ｒ２＊ddｔh（１）*ｓin（x
（1）＋ｐi／２）＋r2＊ｘ（3）＾2*sin（x
（1）+ｐi）；
曲柄原动件动力学M函数文件ｃｈenｇcrankdｙ.ｍ:
ｆunｃtiony＝chengcｒankdy（x）
%FuｎctionforDｙanmｉcanａlyｓisofcrank
%％Iｎput parａmｅterｓ
％ｘ
（1）＝tｈetａ-1
%x
（2）＝ｄｔｈｅta－1
%ｘ（３）=ｄｄｔheta-1
%x（4）=RxB
%ｘ（5）=RyB
%%0uｔputparameters
%y
（1）=RxA
％y（２）＝ＲyA
％y（３）＝M1
g=９.8; ％重力加速度
ｒ1=０.4;%曲柄长度
rc１=0．２;%质心离铰链A的距离
m1＝1.2;%曲柄质量
Ｊ1=０.0１６;%绕质心转动惯量
Fｘ１=０; Fy1=0; ＭF＝0；%作用于质心的外力和外力矩
ReｄｄA=０； ImｄｄA=0;%铰链A的加速度
y
（1）=m１*ReddA+m1＊ｒc１*ｘ（3）＊coｓ（x
（1）＋ｐi/２）+m1*rc１*x
（2）^2*cｏs（x
（1）+pｉ）－Fｘ1+ｘ（4）;
y（２）=m1*IｍdｄA＋m１*rc１＊x（3）*siｎ（x
（1）+ｐi/2）+m1＊rc１*x（２）^２*ｓｉｎ（x
（1）＋pi）-Fy1+x（5）+m1*g;
y（3）=J1*x（3）－y
（1）＊rc１*sin（x
（1））+y
（2）*rｃ1*coｓ（x
（1））－x（4）＊（r1－rc1）＊sin（x（１））+x（５）*（r1-rc1）*coｓ（x
（1））－MF；
ＲRRII级杆组动力学Ｍ函数文件 ｃhｅnｇrｒrdｙ.m:
fｕnctionｙ=chengrｒｒdy（x）
％FuncｔｉoｎforDyanmiｃaｎａlysisofRRRdayａrdgrｏup
%Iｎpuｔｐａrameteｒs
％x
（1）＝thetａ-2
%ｘ
（2）=ｔheta-3
%x（3）=dtheta-2
%x（4）=dtｈeｔａ－3
%x（5）＝ddtｈeｔａ-2
%x（6）＝ｄdｔｈeｔa-3
%ｘ（7）＝Re[ｄdＢ]
％x（8）=Im[ｄdB］
％x（９）＝Fx3
％x（10）=Ｆy３
%x（11）=Ｍ3
%０utputｐarameters
%y
（1）=RxB
％Y
（2）=RyＢ
%y（3）=RxＣ
%y（4）=RｙC
%ｙ（5）=ＲxD
%y（6）=RyD
g=9．8;%重力加速度
r2＝１.2; r3＝0.8;%两杆的长度
rｃ2=0．6;rc3=0.4;％质心到铰链Ｂ的距离%质心到铰链Ｄ的距离
ｍ２=3;ｍ3=2.2;%两杆的质量
J2＝０.25;J３＝0.0９;％两杆的转动惯量
ReddＤ＝０;ＩmddD=0;
Fｘ2＝0;Fy2=０;
M2=0;%2杆的外力和外力矩
ａ=zeｒｏs（６）;
a（1,1）=1；
a（１,３）=1;
a（2,２）=1；
a（2,4）=1;
a（３,1）＝rc2*sin（x
（1））;
ａ（3,2）=-rc2*cos（ｘ
（1））;
ａ（3,3）=-（ｒ2-ｒc2）*siｎ（ｘ
（1））；
a（3,４）=（r2－rc２）*cｏs（x（１））;
a（4,3）=-1;
a（4,5）＝1;
a（５,4）=-１;
ａ（５,６）=１;
a（６,３）＝（ｒ３-rc3）*sｉn（x
（2））；
ａ（６,4）=-（r3-rc3）*cｏs（x
（2））;
a（６,5）=rc3*siｎ（x
（2））;
ａ（6,６）=-rc3*ｃos（x
（2））；
b＝zeros（６，１）;
ｂ（1，1）=m2*rc2*ｘ（5）*ｃｏs（x
（1）＋pi/2）+m2＊x（７）+ｍ2＊ｒc2＊x（3）^2*coｓ（ｘ（１）+pi）-Fx２;
b（２,1）＝m2*rc2*x（5）*sｉn（x
（1）+pｉ／2）+ｍ2*x（8）+ｍ２*rc2*ｘ（3）^2*sin（x
（1）+pi）－Fy2+ｍ2*ｇ;
b（3,1）=J2*x（５）－M2;
ｂ（4,1）=ｍ３*ｒc3*ｘ（6）*coｓ（x
（2）+ｐi／２）+m3*RｅｄdＤ+ｍ3*rc3＊x（4）＾2*coｓ（x（２）＋pｉ）-x（９）;
b（5,1）=m3*ｒc3*ｘ（６）*sin（x
（2）+pi/２）+m3*IｍｄdD+m3*rｃ3*x（4）^２*sin（x（２）＋pi）-ｘ（1０）+ｍ3＊g;
b（６,１）=J3*x（6）-x（11）;
y=inv（ａ）*ｂ;
RRＰ II级杆组Ｍ函数文件:
fuｎｃtiony=ｃhengrrp（x）
%functiontｏcomputeｔheａｃceleratioｎfｏrRRP bar grouｐ
%Inpｕtparaｍeters
％x
（1）=thｅta-5
%x
（2）＝dｔhｅta-5
％x（3）＝Re［dｄC]
%x（４）=Im[dｄC]
％x（5）=ｄs
%Oｕtputｐarametｅｒs
％ｙ
（1）=ｄdtheta-5
%y
（2）=dｄs
r5=1.2;tｈ6=0;ReddD=0;ImｄｄD=0;
a=[r5*cｏs（ｘ
（1）＋pｉ／２）－cｏs（th6）; r5*sin（ｘ
（1）+pi/2） －sin（th6）];
ｂ=[－r5＊cｏs（ｘ
（1）+pi）0；-r5*siｎ（ｘ
（1）+pｉ）0]＊[ｘ（２）＾2;x（5）]+[ＲeｄdD－ｘ（３）; ImｄdD-x（4）］;
y=inv（a）＊b;
ＲRＰＩI级杆组动力学M函数文件:
fｕnctiony=chengｒrpdy（x）
％FuｎctｉｏnfｏrDｙanｍ5cａnalysｉsofRRPdayaｒdgroup
%Inputpaｒameteｒs
%x
（1）=tｈeｔa-５
％x
（2）=dthetａ-5
%ｘ（3）=ddthｅta-５
%ｘ（4）=dds－6
％ｘ（5）=Ｒｅ[dｄC］
%x（6）＝Ｉｍ[dｄC]
%0uｔputparaｍｅterｓ
％ｙ（１）=RxC
％Y
（2）=RyC
%y（３）＝ＲxE
%y（4）＝RyE
%y（5）=ＲF %移动副的约束反力
g=9.8; ％重力加速度
r５=１．2;%杆的长度
ｒc5=0．6;%质心到铰链B的距离
m5=3.6; m6＝6；%杆、块的质量
J5=0.45；
Fx5=0;Fy5=0；
Fx6＝1000； Fy６=0;
M5=0;
ｔh6＝0;
a=zｅros（5）;
a（１,1）=1;
ａ（１，3）=1；
a（2,２）=1;
a（２,4）=1;
a（3，1）=ｒc5*sin（ｘ
（1））;
a（3，２）＝-rｃ5*cos（x
（1））;
a（3,３）＝-（r5-ｒc5）*siｎ（x
（1））；
a（3,4）＝（ｒ５－ｒc5）*cｏｓ（x
（1））；
a（4,３）=-1;
ａ（4,５）=-sin（th6）;
ａ（5,４）=-1;
a（5，５）＝cｏs（th６）;
b=ｚeｒos（5，1）;
b（1,1）=m5＊x（5）+m5*rc5*ｘ（3）＊cos（ｘ（１）+pi／2）+ｍ5*rc５＊x
（2）^2*cｏs（x
（1）+pｉ）-Fx5;
b（2,１）=ｍ5*x（6）+m５*rc５*x（3）*siｎ（ｘ
（1）+pｉ/2）+m5*ｒc5*ｘ
（2）^2＊sin（x
（1）+pi）-Fy5＋m5*ｇ;
b（３,1）=Ｊ5*x（３）-M５;
b（４,1）=m6＊x（4）*cｏs（tｈ6）－Fx６;
b（５,1）=m６*x（４）＊sｉn（th６）－Fx6+ｍ6*ｇ；
y＝inv（a）*ｂ;
三 系统仿真框图
进入MATＬAB,在命令栏中键入Simulｉｎｋ进入仿真界面，根据信息传递的逻辑关系,建立仿真系统框图如图3-1.然后设定各环节的初始参数，即可以对机构进行运动学仿真分析,再利用ＭAＴLＡB的ｐlot命令根据需要绘制曲线。
图3-1
四仿真的实现
再设计完成仿真框图之后,为了进行仿真还必须设定初始参数值。
连杆机构杆长已经在ｓimｕlink框图中给定,如果设定
初始夹角为62
=１0rad／ｓ，曲柄1作匀速转动（即
），接下来要确定杆2，3的角位移和角速度，杆5的角位移和角速度,滑块的速度。
可以利用辛普森方法（在MATＬAB命令框中输入M函数为rrｒpｏsi）求得
=0.３61２rａd/ｓ，
=１．8１01rａd/s，再利用MＡTＬＡB（在命令框输入rrrｖeｌ）求出W2=-2.2３45, Ｗ3=3.３250,再利用杆3的角位移和角速度、杆5的角位移求得（在MATLAB命令框中输入Ｍ函数为compvｅｌ）W5=0.６９62，ds=－3.132３。
对仿真框图中各积分器设定参数变量ｘ并在matlab命令框输入变量x＝［62*pｉ／180１0０．3612１.810１－2．23453.3250-41＊ｐi/１8０0．６96２－３.１３23];其中初始数值分别对应:
tｈｅtａ-1、oｍega－1、thetａ－2、oｍegａ－2、omｅga－３、tｈetａ-5,ｏmega-5ｄs，以及仿真时间为1s，后进行仿真,利用ＭＡTＬＡB中的ploｔ绘图命令把角速度曲线分别绘制出来。
在ＭＡTＬAB命令中键入：
pｌoｔ（tout,sｉｍouｔ（:
1））,plｏｔ（ｔouｔ,siｍｏut（:
2））,pｌot（ｔｏｕt,sｉmoｕt（：
，３））,ｐlot（tｏｕt，siｍout（：
，4））,plｏt（ｔouｔ,siｍout2（:
，5））,即可得到点A的水平方向、垂直方向的约束反力、驱动力矩M1及其所作功W1的变化曲线,如图所示。
转动副A水平方向力 转动副A垂直方向力
曲柄上作用的力矩M１ 曲柄力矩所做的功W1
移动副F的约束反力
升级会员 