vr是变化的,是两物体之间的相对速度,x是初始值,x=[0 0 0 0;0 0.01 0 0.01];我自己编了一个如下:
function dz=nianhuamoban(t,z,c1m1,k1m1,Ffm1,c2m2,k2m2,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=1;c1=0.01;c2=0.01;mus=0.6;
N=15;alpha=0.012;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=5;
x=[0 0 0 0;0 0.01 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);
if abs(vr)<=ep
Ff=min(abs(Fs),Pn)*sign(Fs)-Ff0;
end
if vr>ep
Ff=(mus-alpha*vr)*N*sign(vr)-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.1;
while(1)
h=h+1;
tspan=(a:1:b)*0.1;
[t,x]=ode45(@nianhuamoban,tspan,x(2,:),[],c1m1,k1m1,Ffm1,c2m2,k2m2,Ffm2);
Fs=k1*x(2,1)+c1*x(2,2)-k2*x(2,3)-c2*x(2,4);
Fs1(h)=Fs;
vr=v0+x(2,4)-x(2,2);
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*sign(vr)-Ff0;
end
if -vr>=ep
Ff=-(mus+alpha*vr)*N*sign(vr)-Ff0;
end
% g(h)=y(h,2)-y(h-1,2);
% if abs(g(h))<=0.001
% y(h,2)=y(h-1,2);
% end
% G(h)=y(h,4)-y(h-1,4);
% if abs(G(h))<=0.001
% y(h,4)=y(h-1,4);
% end
Ffm1=Ff/m1;
Ffm2=Ff/m2;
t0(h)=t(2);
a=a+1;
b=b+1;
if a==10000
break;
end
end
t00=t0';
% y=y/w1*u0;
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));
但是出来结果跟文献里的不相符啊,麻烦无水前辈指教
[ 本帖最后由 byandby 于 2009-12-12 16:47 编辑 ] |