单自由度机械系统动力学--牛头刨床运动例题
单自由度机械系统动力学作业
题目：
图1所示为一牛头刨床。各构件长度为：L 1=110mm ，L 3=540mm ，L 4=135mm ；尺寸H =580mm ，H 1=380mm 。导杆3重量G 3=200N ，质心S 3位于导杆中心，导杆绕
S 3的转动惯量J 3=1.1kg ⋅m 2。滑枕5的重量G 5=700N 。其余构件重量均可不计。电动
机型号为Y100L2-4，电动机轴至曲柄1的传动比i =23.833，电动机转子及传动齿轮等折算到曲柄上的转动惯量J 1=133.3kg ⋅m 2。刨床的平均传动效率η=0.85。空行程时作用在滑枕上的摩擦阻力F f =50N ，切削某工件时的切削力和摩擦阻力如图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为等效构件，等效转动惯量为
1
J e =J 1+J 3(
ω32G 3v S 32G 5v 52
) +() +() (a) ω1g ω1g ω1
式中：g 为重力加速度，v S 3为导杆3质心的速度，v 5为滑枕的速度。
等效驱动力矩可由电动机机械特性导出，设M m 、M de 分别为电动机输出力矩和等效驱动力矩，两者有如下关系：
M de =iM m (b)
式中i 为电动机轴和曲轴间的传动比。
电动机轴转速ωm 和曲柄转速ω1间有如下关系：
ω1=
将式(b)和式(c)代入电动机机械特性
ωm
i
(c)
2
(d) M m =a +b ωm +c ωm
可得
M de =ai +bi 2ω1+ci 3ω12 (e)
将传动比i 的值和课本例题3.2.2中求出的系数a 、b 、c 的值代入式(e)，得到等效驱动
力矩
M de =-14846+6076.8ω1-580.26ω12 (f)
等效阻力矩M re 中只计入滑枕上的摩擦阻力F f 和切削阻力F r ，以及导杆的重力G 3：
M re =(-G 3v S 3y -|(F f +F r ) v 5|)/(ωη1) (g)
式中v S 3y 为导杆3重心S 3的y 向速度。 等效力矩M e 为
M e =M de +M re (h)
等效力矩和等效转动惯量均随机构位置而变化。需将曲柄运动周期分成k 个等份（k 可取为60），对每一机构位置计算等效力矩M e 和等效转动惯量J e 。
（3）运动方程的求解
本题属于等效力矩同时为等效构件转角和角速度的函数，而等效力矩M e 的表达式中ω与ϕ可以分离，即可以表达为两个函数的和，其中一个等效驱动力矩M de 为角速度ω的函数，另一个等效阻力矩M re 为转角ϕ的函数。这样采用能量形式的运动方程求解更为简便、
2
快速。
已知等效驱动力矩M de 的表达式为式（f ），等效阻力矩M re 的表达式为式（g ），设
M re =M r (ϕ) ，为已知量。由能量形式的运动方程，对从ϕ1到ϕ2的区间，可以写出
ϕ21122
J e 2ω2-J e 1ω1=⎰(M de +M re )d ϕ (i)
ϕ1
22
式中，ωi 、J ei 为与角ϕi 相对应的位置的角速度和等效转动惯量。 用梯形公式求积分，式(i)可写为
11∆ϕ2
J e 2ω2-J e 1ω12=[M de (ϕ1) +M re (ϕ1) +M de (ϕ2) +M re (ϕ2)] (j) 222
用式（e ）和M re =M r (ϕ) 代入，得到
2
(J e 2-∆ϕc ) ω2-∆ϕb ω2-∆ϕ[M r (ϕ1) +M r (ϕ2) +2a +b ω1+c ω12]-J e 1ω12=0 (k)
这是一个以ω2为未知数的一元二次方程，如果ω1已知，ω2便可很容易地求出。 同理，对第i 个区间，即ϕi 和ϕi +1之间的区间，可以有如下递推公式：
A i ωi 2+1+B i ωi +1+C i =0 (l)
式中：
A i =[(J e ) i +1-∆ϕc ]B i =-∆ϕb
C i =-{∆ϕ[M r (ϕi ) +M r (ϕi +1) +2a +b ωi +c ωi 2]+(J e ) i ωi 2}
用此递推公式，当已知初始条件ϕ=ϕ1、ω=ω1时便可逐步求出各位置的角速度。 计算空载启动后的稳态响应不必取初值ω1=0，为在计算中迅速收敛，可任意取一接近电动机额定角速度的初值，如取ω1=6.5rad /s ，则不到两周便求出稳态解，如图3所示。可以看出，在空载下速度波动很小。
开始刨削后的加载过程的初值可取空载稳态ϕ1=360︒时的ω1值。加载过程如图4所示，最后得到切削时的稳态响应如图5中曲线(2)所示，可以看出，负载的波动导致了较大
的速度波动。将图5、图4与图2对比，可以很清楚地看出，工作循环中的两段有切削力的部分基本上与ω1变化中两次降速的位置相对应。
3
)
s d /r a (/1ω0
[***********]50400
φ1
图3 空载启动后曲柄的稳态运动规律
6.86.76.66.56.4
)
s /r a d (6.3/1ω6.26.165.90
100
200
300
400
500
600
700
800
φ1
图4 开始刨削工件的加载过程
4
6.86.7
6.66.56.4
ω1/(r a d /s )
6.36.26.165.90
50
100
150
200
250
300
350
400
φ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=50 && k
5
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);
6
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]
7
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