MATLAB机械分析.docx
《MATLAB机械分析.docx》由会员分享，可在线阅读，更多相关《MATLAB机械分析.docx（25页珍藏版）》请在冰豆网上搜索。
MATLAB机械分析
作业一凸轮轮廓曲线设计
凸轮机构是一种高副机构，广泛应用于各种机械的运动控制装置中，其功能类似于连杆机构中的函数发生，只不过凸轮可以实现的输入和输出函数关系更加复杂多样，但凸轮的制造费用相对比较高、承载能力比较低。
凸轮机构的运动设计主要包括从动件运动规律的确定和凸轮轮廓线的设计等。
通常是先确定从动件的运动规律，然后根据从动件的运动规律确定凸轮的轮廓曲线。
1.1凸轮从动件运动规律确定
常用的凸轮从动件运动规律主要有多项式运动规律、等速运动规律、等加速等减速运动规律、余弦加速度运动规律、正弦加速度运动规律等。
本文中凸轮从动件的运动规律是推程为五次多项式运动规律，回程是正弦加速度运动规律。
其中推程角为80°，远休止角为40°，回程角为150°，近休止角为90°。
1.2凸轮轮廓曲线设计
已知从动件运动规律的前提下，利用从动件的反转运动可以得到从动滚子中心的运动轨迹，该轨迹称为凸轮的理论轮廓曲线。
如图1所示，假设B0为滚子初始位置的中心，则r=OB0为凸轮的基圆半径。
当凸轮逆时针转过角度δ时，从动件逆时针摆动角度φ，此时从动件滚子中心到达点B1。
利用反转运动规律，将从动件以凸轮回转中心O为中心顺时针转动角度δ，到达图中A′B′1位置。
以凸轮回转中心O为原点建立坐标系XOY，利用几何法可解出图中点B′1的坐标，点B′1即为凸轮理论轮廓线上的一点。
令凸轮转动角度δ分别取到0°到360°中的每一个值，则可以得到凸轮理论轮廓曲线在坐标系XOY中的坐标描述。
已知条件：
基圆半径ro，摆杆转动中心A与凸轮转动中心O的距离d，摆杆AB的长度l，滚子半径rr，摆杆转动中心A和凸轮转动中心O的连线OA与水平X轴的夹角α。
图1摆动滚子从动件盘形凸轮轮廓曲线设计
摆杆初始摆角
所以点
的坐标为
当凸轮顺时针转动时，从动件逆时针反转运动，此时
1.3实例分析
利用MATLAB软件编写凸轮轮廓线绘制程序，要求用户输入参数如下：
请输入凸轮基圆半径（mm）：
80
输入从动件滚子半径（mm）：
30
请输入从动摆杆长度（mm）：
85
请输入凸轮转动中心与摆杆转动中心距离（mm）：
100
请输入摆杆转动中心和凸轮转动中心的连线与水平X轴的夹角（度）：
30
请选择凸轮的转动方向（顺时针为0，逆时针为1）：
1
根据输入的条件和凸轮轮廓曲线计算公式，MATLAB软件输出的凸轮理论轮廓曲线如图2所示。
图2凸轮理论轮廓曲线
利用SolidWorks三维建模软件建立凸轮的三维模型，如图3所示。
图3凸轮三维模型
1.4凸轮系
在凸轮设计时，经常需要设计一系列凸轮，这些凸轮之间有相位差。
假设第二个凸轮与第一个凸轮之间的相位角为60°。
利用MATLAB编制程序，将凸轮轮廓线上的每一个点坐标取出，乘以转换矩阵T，就得到凸轮二的轮廓线的各店的坐标。
最后绘制凸轮二的理论轮廓线如图4。
图4凸轮二的理论轮廓曲线
1.5附录1MATLAB程序
（1）从动件运动规律MATLAB程序
functionphi=yundongguilv（delta）%%定义从动件的运动规律函数
phimax=pi/4;%%从动件的最大摆角为45度
ifdelta<=80*pi/180%%推程
phi=phimax*（10*（delta/（80*pi/180））^3-15*（delta/（80*pi/180））^4+6*（delta/（80*pi/130
80））^5）;%%推程角为80度，推程为五次多项式运动规律
elseifdelta<=120*pi/180%%远休止
phi=phimax;
elseifdelta<=270*pi/180%%回程
phi=phimax*（1-（delta-120*pi/180）/（150*pi/180）+sin（2*pi*（delta-120*pi/180）/（150*pi/180））/（2*pi））;%%回程角为150度，回城为正弦加速度运动规律
else%%近休止
phi=0;
end
end
（2）凸轮理论轮廓线绘制MATLAB程序
clear;
clc;
ro=input（'请输入凸轮基圆半径（mm）：
'）;%%要求用户输入凸轮基圆半径
rr=input（'输入从动件滚子半径（mm）：
'）;%%要求用户输入从动件滚子半径
l=input（'请输入从动摆杆长度（mm）：
'）;%%要求用户输入从动摆杆长度
d=input（'请输入凸轮转动中心与摆杆转动中心距离（mm）：
'）;%%要求用户输入凸轮转动中心与摆杆转动中心距离
alpha=input（'请输入摆杆转动中心和凸轮转动中心的连线与水平X轴的夹角（度）：
'）;%%要求用户输入摆杆转动中心和凸轮转动中心的连线与水平X轴的夹角
alpha=alpha*pi/180;
beta=acos（（l^2+d^2-（ro+rr）^2）/（2*l*d））;%%摆杆初始摆角
direction=input（'请选择凸轮的转动方向（顺时针为0，逆时针为1）：
'）;%%要求用户选择凸轮的转动方向
switchdirection
case0
fori=0:
360
delta=i*pi/180;
phi=yundongguilv（delta）;
h=（l^2+d^2-2*l*d*cos（phi+beta））^0.5;
eta=acos（（h^2+d^2-l^2）/（2*h*d））;
theta=eta+alpha+delta;
lunkuo（i+1,1）=h*cos（theta）;
lunkuo（i+1,2）=h*sin（theta）;
lunkuo（i+1,3）=0;
end
case1
fori=0:
360
delta=i*pi/180;
phi=yundongguilv（delta）;
h=（l^2+d^2-2*l*d*cos（phi+beta））^0.5;
eta=acos（（h^2+d^2-l^2）/（2*h*d））;
theta=eta+alpha-delta;
lunkuo（i+1,1）=h*cos（theta）;
lunkuo（i+1,2）=h*sin（theta）;
lunkuo（i+1,3）=0;
end
end
save'F:
\lunkuo.txt'lunkuo-ascii;
plot（lunkuo（:
1）,lunkuo（:
2））;%%绘制凸轮一的轮廓曲线
fori=1:
361
p1（1,1）=lunkuo（i,1）;
p1（2,1）=lunkuo（i,2）;
p1（3,1）=lunkuo（i,3）;
gamma=60*pi/180;
T=[cos（gamma）-sin（gamma）0;sin（gamma）cos（gamma）0;001];
p2=T*p1;
lunkuo2（i,1）=p2（1,1）;
lunkuo2（i,2）=p2（2,1）;
lunkuo2（i,3）=p2（3,1）;
end
save'F:
\lunkuo2.txt'lunkuo2-ascii;
figure
（2）;
plot（lunkuo2（:
1）,lunkuo2（:
2）,'r'）%%绘制凸轮二的轮廓曲线
1.6附录2凸轮二维CAD图纸
见附页。
作业二ABBIRB660机械手运动学分析
2.1IRB660机械手简介
IRB660是一款专用堆垛机器人，其速度、到达距离和有效载荷在同类产品的市场上一枝独秀。
超高速4轴运行机构、3.15米到达距离加上250kg的有效载荷，使IRB660成为袋、盒、板条箱和瓶子等包装材料的理想堆垛工具。
IRB660机械手的外形和尺寸如图1所示。
图1IRB660机械手示意图
ABBIRB660机械手主要由固定底座、回转台、大臂、小臂、工作头、随动油缸和连杆等部分组成。
ABBIRB660机械手是一种四自由度的平行四边形机构的工业机器人，四个电机分别驱动回转台回转、大小臂的俯仰以及工作头的回转。
四个转动部件的工作范围见表1。
表1机械手各部件工作范围
2.2机械手工作空间分析
图2IRB660机械手机构简图与坐标系建立示意图
不考虑回转台的回转和工作台的转动，可以将机械手简化成机构简图，如图2所示。
杆件1和杆件2可以分别绕着关节A和关节B转动，杆件1的转动范围为
，杆件2的转动范围为
，如图2所示，其中逆时针转动为正，顺时针转动为负。
根据图2建立的各杆件坐标系，可以写出IRB660机械手的连杆参数，如表2所示。
表2IRB660机械手的连杆参数
连杆i
连杆长度ai-1
连杆转角αi-1
连杆偏距di
关节角θi
变量范围
1
0
0°
0
β1+90°
θ1
2
1280
0°
0
-90°-β1+β2
θ2
3
1350
0°
0
0
0
关节角θ1的转动范围是5°~132°，考虑到限位块的限位作用，关节角θ2的转动范围为-160°~-20°。
为了避免杆件之间的相互干涉关系，关节角θ1和关节角θ2之间存在相互约束关系，如图3所示。
图3关节角θ1和关节角θ2之间的角度约束关系
根据机械手的各连杆参数，建立连杆之间的D-H变换矩阵：
将表2中的数据代入可以得到各连杆坐标系之间的变换矩阵为
假设坐标系3中一点P，其坐标描述为
。
该点在坐标系0中描述为
，则有
取连杆3的关节点C为研究对象，点C在坐标系3中的坐标为
，利用转换矩阵
，可以找出点C在坐标系0（即固定坐标系）中的坐标
。
则机械手工作头中心D在坐标系0中的坐标为
。
利用MATLAB绘制机械手工作头中心D的工作范围，如图4所示。
图4机械手工作范围平面示意图
将MATLAB绘制的机械手工作范围（图4）与图1中的机械手工作范围对比，发现两者基本吻合，说明工作空间的仿真分析正确。
从图4和图1中可以看出，IRB660码垛机械手的举升重物为大臂回转中心水平方向0.5m到2.8m，竖直方向-1.5m到1.5m的范围内，举升范围大，能力强。
2.3机械手静力学分析
不考虑机械手各杆件的质量，当机械手工作头满载时，大臂铰点A和小臂铰点B的受力如图5所示。
图5机械手大臂和小臂静力学分析
根据图5，可以得到：
通过上面几个式子可以发现，MA和MB都是大臂转角β1和小臂转角β2的函数。
利用MATLAB函数编制程序绘制铰点A和铰点B的力矩图，如图5所示。
通过MATLAB可以得到大臂铰点A受到的最大力矩为7.2128kNm，小臂铰点B受到的最大力矩为4.0250kNm。
图6大臂和小臂铰点静力学分析
2.4附录MATLAB程序
（1）关节角1和关节角2转动角度关系计算：
forbeta1=-85:
42
forbeta2=-120:
20
if（-beta1+beta2-90<=-20&&-beta1+beta2-90>=-160）
theta1=90+beta1;
theta2=-beta1+beta2-90;
plot（theta1,theta2,'*'）;
holdon
end
end
end
title（'关节角θ1和θ2之间的约束关系'）;
xlabel（'关节角θ1'）;
ylabel（'关节角θ2'）;
（2）机械手工作范围计算：
clc;
clear;
fori=1:
50000
beta1=-85+127*rand;
beta2=-120+140*rand;
if（-beta1+beta2-90<=-20&&-beta1+beta2-90>=-160）
theta1=（90+beta1）*pi/180;
theta2=（-beta1+beta2-90）*pi/180;
T1=[cos（theta1）-sin（theta1）00;sin（theta1）cos（theta1）00;0010;0001];
T2=[cos（theta2）-sin（theta2）01280;sin（theta2）cos（theta2）00;0010;0001];
T3=[1001350;0100;0010;0001];
C3=[0;0;0;1];
C0=T1*T2*T3*C3;
P=[C0（1,1）+260,C0（2,1）-247];
plot（P（1,1）,P（1,2）,'*'）;
holdon;
end
end
（3）机械手静力学分析
clc;
clear;
a=linspace（-85*pi/180,42*pi/180,500）;
b=linspace（-120*pi/180,20*pi/180,500）;
[aa,bb]=meshgrid（a,b）;
Ma=2500*（-1.28*sin（aa）+1.35*cos（bb）+0.26）;
Mb=2500*（1.35*cos（bb）+0.26）;
beta1=linspace（-85,42,500）;
beta2=linspace（-120,20,500）;
subplot（1,2,1）;
mesh（beta1,beta2,Ma）;
title（'大臂铰点A力矩分析'）;
xlabel（'大臂转角β1'）;
ylabel（'小臂转角β2'）;
subplot（1,2,2）;
mesh（beta1,beta2,Mb）;
title（'小臂铰点B力矩分析'）;
xlabel（'大臂转角β1'）;
ylabel（'小臂转角β2'）;
作业三结构动态仿真
某三层钢筋混凝土结构（如图1），结构的各层特性参数分别为
，
，
。
地震波采用美国El-Centro南北方向地震波，延续时间为20s，分别计算三层结构在地震波作用下的位移、速度和加速度。
图1结构力学简化图
3.1建立系统微分方程
对第一层结构，有
（1）
对第二层结构，有
（2）
对第三层结构，有
（3）
将
（1）、
（2）和（3）式用矩阵方式表示为
（4）
其中：
，
，
，
，
，
整理式（4），得
（5）
由于
，因此
（I为三阶单位矩阵），
，所以可以将式（5）整理为：
（6）
3.2Simulink仿真分析
图2Simulink系统仿真模型
如图2，利用MATLAB自带的Simulink工具箱模拟El-Centro地震波对建筑结构的作用。
将三层结构在地震作用下的位移、速度和加速度混频输出到矩阵xa中。
3.3Simulink仿真结果
根据仿真结果，绘制出三层结构的位移图、速度图和加速度图，如图3、4、5所示。
表1是每一层结构的位移、速度、加速度的最大值。
图3各层建筑结构位移图
图4各层建筑结构速度图
图5各层建筑结构加速度图
表1每一层结构的位移、速度、加速度的最大值
结构
位移（m）
速度（m/s）
加速度（m/s2）
第一层
0.04618
0.43777
5.80162
第二层
0.08311
0.73147
9.72524
第三层
0.10593
0.95163
12.70139
通过表1可以看出，随着结构高度的增大，在地震中结构产生的位移、速度和加速度都在逐渐增大。
因此可以知道，在地震中，楼层越高，震感会越强烈，这与实际体验是一致的。
3.4附录MATLAB程序
（1）输入仿真参数：
clc
m=2923.39;
c=6.37374e3;
c1=c;c2=c;c3=c;
k=1391.06e3;
k1=k;k2=k;k3=k;
C=[c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3];
K=[k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3];
（2）结果分析程序
figure
（1）;
subplot（3,1,1）;
plot（tout,xa（:
1）,'b-','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第一层结构位移x1/m'）;
subplot（3,1,2）;
plot（tout,xa（:
2）,'r:
','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第二层结构位移x2/m'）;
subplot（3,1,3）;
plot（tout,xa（:
3）,'k-.','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第二层结构位移x3/m'）;
figure
（2）;
subplot（3,1,1）;
plot（tout,xa（:
4）,'b-','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第一层结构速度v1（m/s）'）;
subplot（3,1,2）;
plot（tout,xa（:
5）,'r:
','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第二层结构速度v2（m/s）'）;
subplot（3,1,3）;
plot（tout,xa（:
6）,'k-.','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第二层结构速度v3（m/s）'）;
figure（3）;
subplot（3,1,1）;
plot（tout,xa（:
7）,'b-','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第一层结构加速度a1（m/s^2）'）;
subplot（3,1,2）;
plot（tout,xa（:
8）,'r:
','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第二层结构加速度a2（m/s^2）'）;
subplot（3,1,3）;
plot（tout,xa（:
9）,'k-.','linewidth',2）;grid;
xlabel（'时间t/s'）;
ylabel（'第二层结构加速度a3（m/s^2）'）;
x1_max=max（abs（xa（:
1）））
x2_max=max（abs（xa（:
2）））
x3_max=max（abs（xa（:
3）））
v1_max=max（abs（xa（:
4）））
v2_max=max（abs（xa（:
5）））
v3_max=max（abs（xa（:
6）））
a1_max=max（abs（xa（:
7）））
a2_max=max（abs（xa（:
8）））
a3_max=max（abs（xa（:
9）））
升级会员 