单自由度机械系统动力学——牛头刨床运动例题.doc

单自由度机械系统动力学作业
题目：
图1所示为一牛头刨床。各构件长度为：，，；尺寸，。导杆3重量，质心位于导杆中心，导杆绕的转动惯量。滑枕5的重量。其余构件重量均可不计。电动机型号为Y100L2-4，电动机轴至曲柄1的传动比，电动机转子及传动齿轮等折算到曲柄上的转动惯量。刨床的平均传动效率。空行程时作用在滑枕上的摩擦阻力，切削某工件时的切削力和摩擦阻力如图2所示。
1） 求空载启动后曲柄的稳态运动规律；
2） 求开始刨削工件的加载过程，直至稳态。

图1 牛头刨床 图2 牛头刨床加工某工件时的负载图

解：
（1）运动分析
可以用解析法列出各杆角速度、各杆质心速度的表达式。但为简便起见，现调用改自课本附录Ⅰ中的Matlab子程序来进行计算。图1中给出了构件和运动副的编号。先调用子程序crank分析点②的运动学参数，再调用子程序vosc进行滑块2—导杆3这一杆组的运动学分析，然后再调用子程序vguide进行小连杆4—滑枕5这一杆组的运动学分析。这一段的Matlab程序如下：
crank(1,2,L(1),TH(1),W(1));
vosc(2,3,4,L(3));
vguide(4,5,L(4));
其中：L(i)、TH(i)、W(i)分别表示第i个杆的长度、位置角、角速度。
（2） 等效转动惯量和等效力矩
取曲柄1为等效构件，等效转动惯量为
(a)
式中：g为重力加速度，为导杆3质心的速度，为滑枕的速度。
等效驱动力矩可由电动机机械特性导出，设、分别为电动机输出力矩和等效驱动力矩，两者有如下关系：
(b)
式中i为电动机轴和曲轴间的传动比。
电动机轴转速和曲柄转速间有如下关系：
(c)
将式(b)和式(c)代入电动机机械特性
(d)
可得
(e)
将传动比i的值和课本例题3.2.2中求出的系数a、b、c的值代入式(e)，得到等效驱动力矩
(f)
等效阻力矩中只计入滑枕上的摩擦阻力和切削阻力，以及导杆的重力：
(g)
式中为导杆3重心的y向速度。
等效力矩为
(h)
等效力矩和等效转动惯量均随机构位置而变化。需将曲柄运动周期分成k个等份（k可取为60），对每一机构位置计算等效力矩和等效转动惯量。
（3） 运动方程的求解
本题属于等效力矩同时为等效构件转角和角速度的函数，而等效力矩的表达式中与可以分离，即可以表达为两个函数的和，其中一个等效驱动力矩为角速度的函数，另一个等效阻力矩为转角的函数。这样采用能量形式的运动方程求解更为简便、快速。
已知等效驱动力矩的表达式为式（f），等效阻力矩的表达式为式（g），设，为已知量。由能量形式的运动方程，对从到的区间，可以写出
(i)
式中，、为与角相对应的位置的角速度和等效转动惯量。
用梯形公式求积分，式(i)可写为
(j)
用式（e）和代入，得到
(k)
这是一个以为未知数的一元二次方程，如果已知，便可很容易地求出。
同理，对第i个区间，即和之间的区间，可以有如下递推公式：
(l)
式中：

用此递推公式，当已知初始条件、时便可逐步求出各位置的角速度。
计算空载启动后的稳态响应不必取初值，为在计算中迅速收敛，可任意取一接近电动机额定角速度的初值，如取，则不到两周便求出稳态解，如图3所示。可以看出，在空载下速度波动很小。
开始刨削后的加载过程的初值可取空载稳态时的值。加载过程如图4所示，最后得到切削时的稳态响应如图5中曲线(2)所示，可以看出，负载的波动导致了较大的速度波动。将图5、图4与图2对比，可以很清楚地看出，工作循环中的两段有切削力的部分基本上与变化中两次降速的位置相对应。

图3 空载启动后曲柄的稳态运动规律

图4 开始刨削工件的加载过程
(2)
(1)

图5 空载与切削时的稳态响应

Matlab程序：
[main.m]
global P VP %各点位置与速度为全局变量
P=zeros(5,2);
VP=zeros(5,2);
P(3,2)=-0.38;
P(5,2)=0.2;
Je=zeros(1,61);
Mre=zeros(1,61);
Mre0=zeros(1,61);
DeltaPhi=pi/30;
%准备工作，先计算各个位置时的等效转动惯量Je，等效阻力矩Mre
%因为本题等效转动惯量与等效阻力矩均只与机构位置有关，与角速度无关，设曲柄角速度为1进行计算
for k=1:60
crank(1,2,0.11,2*pi-(k-1)*DeltaPhi,1);
W3=vosc(2,3,4,0.54);
vguide(4,5,0.135);
Vs3=sqrt(VP(4,1)^2+VP(4,2)^2)/2;
Je(k)=133.3+1.1*W3^2+200/10*Vs3^2+700/10*VP(5,1)^2;
if((k>=33 && k<=43)||(k>=50 && k<=59))
F=9500;
else
F=50;
end
Mre(k)=(-200*VP(4,2)/2-abs(F*VP(5,1)))/0.85; %刨削工件时的阻力矩
Mre0(k)=(-200*VP(4,2)/2-abs(50*VP(5,1)))/0.85; %空载阻力矩
end
Je(61)=Je(1); %第61点值与第1点值相同，只是为了方便后面的迭代计算
Mre(61)=Mre(1);
Mre0(61)=Mre0(1);

n=0; %记录迭代次数，其实没什么用
w=zeros(1,61);
w(1)=6.5; w(61)=1;
while abs(w(61)-w(1))/w(61)>=1e-4
if n==0
n=1;
else
w(1)=w(61); %更新原点比较值
n=n+1;
end
for k=1:60
A=Je(k+1)-DeltaPhi*(-580.26);
B=-DeltaPhi*6076.8;
C=-(DeltaPhi*(Mre0(k)+Mre0(k+1)+2*(-14846)+6076.8*w(k)+(-58 0.26)*w(k)^2)+Je(k)*w(k)^2);
w(k+1)=(-B+sqrt(B^2-4*A*C))/(2*A);
end
end
Phi=0:6:360;
plot(Phi,w); %绘制空载稳态
xlabel('\phi_1');
ylabel('\omega_1/(rad/s)');
figure(3); %用来绘制空载与工作状态稳态对比
plot(Phi,w); %绘制空载稳态
xlabel('\phi_1');
ylabel('\omega_1/(rad/s)');

w0=w(61);
w=zeros(1,121); %取加载后两个周期的数据
w(1)=w0;
for k=1:120
sk=mod(k-1,60)+1; %sk取值范围为1~60
A=Je(sk+1)-DeltaPhi*(-580.26);
B=-DeltaPhi*6076.8;
C=-(DeltaPhi*(Mre(sk)+Mre(sk+1)+2*(-14846)+6076.8*w(k)+(-580.26)*w(k)^2)+Je(sk)*w(k)^2);
w(k+1)=(-B+sqrt(B^2-4*A*C))/(2*A);
end
Phi=0:6:120*6;
figure;
plot(Phi,w); %绘制加载过程
xlabel('\phi_1');
ylabel('\omega_1/(rad/s)');

%计算加载后达到的稳态响应，其实之前的计算值已经满足要求精度了，所以while里的语句不会执行
w(1:61)=w(61:121);
while abs(w(61)-w(1))/w(61)>=1e-4
w(1)=w(61);
for k=1:60
A=Je(k+1)-DeltaPhi*(-580.26);
B=-DeltaPhi*6076.8;
C=-(DeltaPhi*(Mre(k)+Mre(k+1)+2*(-14846)+6076.8*w(k)+(-580.26)*w(k)^2)+Je(k)*w(k)^2);
w(k+1)=(-B+sqrt(B^2-4*A*C))/(2*A);
end
end
Phi=0:6:360;
figure(3);
hold on;
plot(Phi,w(1:61)); %绘制工作状态稳态,与空载对比

[crank.m]
function crank(N1,N2,R,TH,W)
global P VP
VP(N1,1)=0;
VP(N1,2)=0;
RX=R*cos(TH);
RY=R*sin(TH);
P(N2,1)=P(N1,1)+RX;
P(N2,2)=P(N1,2)+RY;
VP(N2,1)=-RY*W;
VP(N2,2)=RX*W;

[vosc.m]
function [W]=vosc(N1,N2,N3,R)
global P VP
TH=posc(N1,N2,N3,R);
R2=sqrt((P(N2,1)-P(N1,1))^2+(P(N2,2)-P(N1,2))^2);
W=((VP(N1,2)-VP(N2,2))*cos(TH)-(VP(N1,1)-VP(N2,1))*sin(TH))/R2;
VP(N3,1)=VP(N2,1)-W*R*sin(TH);
VP(N3,2)=VP(N2,2)+W*R*cos(TH);

[posc.m]
function [TH]=posc(N1,N2,N3,R)
global P
TH=atan2(P(N1,2)-P(N2,2),P(N1,1)-P(N2,1));
P(N3,1)=P(N2,1)+R*cos(TH);
P(N3,2)=P(N2,2)+R*sin(TH);

[vguide.m]
function vguide(N1,N2,R)
global P VP
TH=pi-asin((P(N2,2)-P(N1,2))/R);
W=-VP(N1,2)/(R*cos(TH));
VP(N2,1)=VP(N1,1)-R*W*sin(TH);

8