单自由度机械系统动力学——牛头刨床运动例题
答:
图3 空载启动后曲柄的稳态运动规律
图4 开始刨削工件的加载过程
图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)+(-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; 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);
解题分析参见这里
参考书籍:
机械动力学(第二版),张策,高等教育出版社