1 单自由度机械系统动力学作业题目：图 1 所示为一牛头刨床各构件长度为：1110Lmm，3540Lmm，4135Lmm；尺寸580Hmm，1380Hmm导杆 3 重量3200GN，质心3S位于导杆中心，导杆绕3S的转动惯量231.1Jkg m滑枕 5 的重量5700GN其余构件重量均可不计电动机型号为Y100L2-4 ，电动机轴至曲柄1 的传动比23.833i，电动机转子及传动齿轮等折算到曲柄上的转动惯量21133.3Jkg m刨床的平均传动效率0.85空行程时作用在滑枕上的摩擦阻力50fFN，切削某工件时的切削力和摩擦阻力如图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 为等效构件，等效转动惯量为精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 1 页，共 8 页 - - - - - - - - -2 2223335513111()()()SeJJJGvGvgg(a) 式中： g 为重力加速度，3Sv为导杆 3 质心的速度，5v为滑枕的速度等效驱动力矩可由电动机机械特性导出，设mM、deM分别为电动机输出力矩和等效驱动力矩，两者有如下关系：demMiM(b) 式中 i 为电动机轴和曲轴间的传动比电动机轴转速m和曲柄转速1间有如下关系：1mi(c) 将式 (b)和式 (c)代入电动机机械特性2mmmMabc(d) 可得23211deMaibici(e) 将传动比i 的值和课本例题3.2.2 中求出的系数a、b、c 的值代入式 (e)，得到等效驱动力矩211148466076.8580.26deM(f) 等效阻力矩reM中只计入滑枕上的摩擦阻力fF和切削阻力rF，以及导杆的重力3G：3351(| ()|)/ ()reS yfrMG vFFv(g) 式中3S yv为导杆 3 重心3S的 y 向速度。
等效力矩eM为edereMMM(h) 等效力矩和等效转动惯量均随机构位置而变化需将曲柄运动周期分成k 个等份（ k 可取为 60） ，对每一机构位置计算等效力矩eM和等效转动惯量eJ3）运动方程的求解本题属于等效力矩同时为等效构件转角和角速度的函数，而等效力矩eM的表达式中与可以分离，即可以表达为两个函数的和，其中一个等效驱动力矩deM为角速度的函数，另一个等效阻力矩reM为转角的函数这样采用能量形式的运动方程求解更为简便、精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 2 页，共 8 页 - - - - - - - - -3 快速已知等效驱动力矩deM的表达式为式（f） ，等效阻力矩reM的表达式为式（g） ，设( )rerMM，为已知量由能量形式的运动方程，对从1到2的区间，可以写出2112222111()d22eedereJJMM(i) 式中，i、eiJ为与角i相对应的位置的角速度和等效转动惯量用梯形公式求积分，式(i)可写为222211112211()()()()222eederedereJJMMMM(j) 用式（ e）和( )rerMM代入，得到222222121111()()()20erreJcbMMabcJ(k) 这是一个以2为未知数的一元二次方程，如果1已知，2便可很容易地求出。
同理，对第i 个区间，即i和1i之间的区间，可以有如下递推公式：2110iiiiiABC(l) 式中：1221()()()()2ieiiiririiieiicAJBMabcbCMJ用此递推公式，当已知初始条件1、1时便可逐步求出各位置的角速度计算空载启动后的稳态响应不必取初值10，为在计算中迅速收敛，可任意取一接近电动机额定角速度的初值，如取16.5/rads， 则不到两周便求出稳态解，如图 3 所示可以看出，在空载下速度波动很小开始刨削后的加载过程的初值可取空载稳态1360时的1值加载过程如图4 所示，最后得到切削时的稳态响应如图5 中曲线 (2)所示，可以看出，负载的波动导致了较大的速度波动将图5、图 4 与图 2 对比，可以很清楚地看出，工作循环中的两段有切削力的部分基本上与1变化中两次降速的位置相对应精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 3 页，共 8 页 - - - - - - - - -4 0501001502002503003504006.546.556.566.576.586.596.66.616.626.6311/(rad/s)图 3 空载启动后曲柄的稳态运动规律01002003004005006007008005.966.16.26.36.46.56.66.76.811/(rad/s)图 4 开始刨削工件的加载过程精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 4 页，共 8 页 - - - - - - - - -5 0501001502002503003504005.966.16.26.36.46.56.66.76.811/(rad/s)图 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 = 1e-4if n = 0 n=1;else w( 1)= w( 61 );% 更新原点比较值 n=n+1;endfor 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)+(-580.26 )* w( k) 2)+ Je ( k)* w( k ) 2); w( k+1)=(-B+sqrt( B2- 4* A* C)/(2* A);endendPhi =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);精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 6 页，共 8 页 - - - - - - - - -7 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);endPhi =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);endendPhi =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 精品学习资料 可选择p d f - - - - - - - - - - - - - - 第 7 页，共 8 页 - - - - - - - - -8 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( N。