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

答：

图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);

 

解题分析参见这里

参考书籍：

机械动力学（第二版），张策，高等教育出版社