关于粘滑振动求助

byandby

大家好，我现在正在研究盘式制动器的粘滑振动，建立了动力学模型，用龙格库塔法仿真计算，但却出不了文献上的粘滑效果，请问大侠哪里出错了，我搞了20天了，找不到原因啊，请大家帮助分析，先谢了。
我的程序如下
m函数
function dz=nianhuamoban(t,z,c1m1,k1m1,c2m2,k2m2,Ffm1,Ffm2)
dz=[z(2);-c1m1*z(2)-k1m1*z(1);z(4);-c2m2*z(4)-k2m2*z(3)]+[0;Ffm1;0;-Ffm2];

主程序：
clear;
clc;
ep=0.005;
m1=1;m2=1;k1=1;k2=2;c1=0.11;c2=0.11;mus=0.6;
N=10;alpha=0.025;
v0=1;
Pn=mus*N;
Ff0=N*(mus-alpha*v0);
c1m1=c1/m1;k1m1=k1/m1;c2m2=c2/m2;k2m2=k2/m2;
a=1;                    
b=3;
x=[0 0 0 0;0 0.1 0 -0.1];     %质量刚度相等时，两个速度相等时或位移相等时收敛
Fs=k1*x(2,1)+c1*x(2,2)-k2*x(2,3)-c2*x(2,4);
vr=v0+x(2,4)-x(2,2);
Ff=(mus-alpha*vr)*N;
if abs(vr)<=ep
    Ff=min(abs(Fs),Pn)*sign(Fs)-Ff0;
end
if vr>ep
    Ff=(mus-alpha*vr)*N-Ff0;
end
if -vr>ep
    Ff=-(mus+alpha*vr)*N*sign(vr)-Ff0;
end
Fs1(1)=Fs;
h=1;
y(1,:)=x(2,:);
Ffm1=Ff/m1;
Ffm2=Ff/m2;
t0(1)=a*0.001;
while(1)
       h=h+1;
       tspan=(a:1:b)*0.01;
       [t,x]=ode45(@nianhuamoban,tspan,x(2,:),[],c1m1,k1m1,c2m2,k2m2,Ffm1,Ffm2);
       vr=v0+x(2,4)-x(2,2);
       Fs=k1*x(2,1)+c1*x(2,2)-k2*x(2,3)-c2*x(2,4);
       Fs1(h)=Fs;
       y(h,:)=x(2,:);
       if abs(vr)<ep
          Ff=min(abs(Fs),Pn)*sign(Fs)-Ff0;
       end
       if vr>=ep
          Ff=(mus-alpha*vr)*N-Ff0;
       end
       if -vr>=ep
          Ff=-(mus+alpha*vr)*N*sign(vr)-Ff0;
       end  
       t0(h)=t(2);   
       a=a+1;
       b=b+1;
       if a==1000
           break;
       end
end
t00=t0';
figure(1)
subplot(2,1,1);
plot(t00,y(:,1));
subplot(2,1,2);
plot(t00,y(:,2));
figure(2)
plot(y(:,1),y(:,2));
figure(3)
subplot(2,1,1);
plot(t00,y(:,3));
subplot(2,1,2);
plot(t00,y(:,4));
figure(4)
plot(y(:,3),y(:,4));

greatchina

你把动力学模型穿上来看看，但从程序看很难看出什么问题。

萨布江胡

请问这个问题最后解决了吗？