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

单自由度机械系统动力学牛头刨床运动例题单自由度机械系统动力学作业 题目： 图1所示为一牛头刨床。各构件长度为：L1=110mm，L3=540mm，L4=135mm；尺寸H=580mm，H1=380mm。导杆3重量G3=200N，质心S3位于导杆中心，导杆绕S3的转动惯量J3=1.1kg×m2。滑枕5的重量G5=700N。其余构件重量均可不计。电动机型号为Y100L2-4，电动机轴至曲柄1的传动比i=23.833，电动机转子及传动齿轮等折2算到曲柄上的转动惯量J1=133.3kg×m。刨床的平均传动效率h=0.85。空行程时作用在滑枕上的摩擦阻力Ff=50N，切削某工件时的切削力和摩擦阻力如图2所示。 1）求空载启动后曲柄的稳态运动规律； 2）求开始刨削工件的加载过程，直至稳态。 图1 牛头刨床 图2 牛头刨床加工某工件时的负载图 解： 运动分析 可以用解析法列出各杆角速度、各杆质心速度的表达式。但为简便起见，现调用改自课本附录中的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个杆的长度、位置角、角速度。 等效转动惯量和等效力矩 取曲柄1为等效构件，等效转动惯量为 1 Je=J1+J3(w32G3vS32G5v52)+ (a) w1gw1gw1式中：g为重力加速度，vS3为导杆3质心的速度，v5为滑枕的速度。 等效驱动力矩可由电动机机械特性导出，设Mm、Mde分别为电动机输出力矩和等效驱动力矩，两者有如下关系： Mde=iMm (b) 式中i为电动机轴和曲轴间的传动比。 电动机轴转速wm和曲柄转速w1间有如下关系： w1=将式(b)和式(c)代入电动机机械特性 wmi (c) 2Mm=a+bwm+cwm (d) 可得 Mde=ai+bi2w1+ci3w12 (e) 将传动比i的值和课本例题3.2.2中求出的系数a、b、c的值代入式(e)，得到等效驱动力矩 Mde=-14846+6076.8w1-580.26w12 (f) 等效阻力矩Mre中只计入滑枕上的摩擦阻力Ff和切削阻力Fr，以及导杆的重力G3： Mre=(-G3vS3y-|(Ff+Fr)v5|)/(wh1) (g) 式中vS3y为导杆3重心S3的y向速度。 等效力矩Me为 Me=Mde+Mre (h) 等效力矩和等效转动惯量均随机构位置而变化。需将曲柄运动周期分成k个等份，对每一机构位置计算等效力矩Me和等效转动惯量Je。 运动方程的求解 本题属于等效力矩同时为等效构件转角和角速度的函数，而等效力矩Me的表达式中w与j可以分离，即可以表达为两个函数的和，其中一个等效驱动力矩Mde为角速度w的函数，另一个等效阻力矩Mre为转角j的函数。这样采用能量形式的运动方程求解更为简便、2 快速。 已知等效驱动力矩Mde的表达式为式，等效阻力矩Mre的表达式为式，设Mre=Mr(j)，为已知量。由能量形式的运动方程，对从j1到j2的区间，可以写出 j2112Je2w2-Je1w12=ò(Mde+Mre)dj (i) j122式中，wi、Jei为与角ji相对应的位置的角速度和等效转动惯量。 用梯形公式求积分，式(i)可写为 11Dj2Je2w2-Je1w12=Mde(j1)+Mre(j1)+Mde(j2)+Mre(j2) (j) 222用式和Mre=Mr(j)代入，得到 2(Je2-Djc)w2-Djbw2-DjMr(j1)+Mr(j2)+2a+bw1+cw12-Je1w12=0 (k) 这是一个以w2为未知数的一元二次方程，如果w1已知，w2便可很容易地求出。 同理，对第i个区间，即ji和ji+1之间的区间，可以有如下递推公式： Aiwi2+1+Biwi+1+Ci=0 (l) 式中： Ai=(Je)i+1-DjcBi=-DjbCi=-DjMr(ji)+Mr(ji+1)+2a+bwi+cwi2+(Je)iwi2用此递推公式，当已知初始条件j=j1、w=w1时便可逐步求出各位置的角速度。 计算空载启动后的稳态响应不必取初值w1=0，为在计算中迅速收敛，可任意取一接近电动机额定角速度的初值，如取w1=6.5rad/s，则不到两周便求出稳态解，如图3所示。可以看出，在空载下速度波动很小。 开始刨削后的加载过程的初值可取空载稳态j1=360°时的w1值。加载过程如图4所示，最后得到切削时的稳态响应如图5中曲线(2)所示，可以看出，负载的波动导致了较大的速度波动。将图5、图4与图2对比，可以很清楚地看出，工作循环中的两段有切削力的部分基本上与w1变化中两次降速的位置相对应。 3 6.636.626.616.66.596.586.576.566.556.54050100150200250300350400f1图3 空载启动后曲柄的稳态运动规律 6.86.76.66.56.46.36.26.165.90100200300400500600700800f1图4 开始刨削工件的加载过程 4 6.86.7(1) 6.66.56.46.36.26.165.9050100150200250300350400(2) f1图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*W32+200/10*Vs32+700/10*VP(5,1)2; if(k>=33 && k<=43)|(k>=50 && k<=59) F=9500; 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(B2-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取值范围为160 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(B2-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(B2-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