急！请各位高手帮忙修改一下滚动轴承松动故障仿真程序，谢谢！~

xufuze12

急，急，急！！！老师让我编一个滚动轴承支撑松动的程序~~（年前必须要弄好，现在好纠结啊）
这是我按照“含不平衡-碰摩-基础松动耦合故障转子-滚动轴承系统非线性动力响应分析” 陈果以及
“转子系统支撑松动的非线性动力学及故障特征” 李振平 2002两篇文章编的一个关于滚动轴承支撑松动故障的程序，自己也修改了好久，但是运行后响应图及轴心轨迹图有问题，出不来！
    请各位帮忙指点一下，看我程序哪个地方出错了，本人在此谢谢大家了！
%运动微分方程
function d=funa(t,y,w)
N=length(y);
m1=4;%两端滚动轴承处等效集中质量
m2=32.1; %转子圆盘等效集中质量
m3=50.0;%轴承支座处等效集中质量
g=9.8;
e=0.00005; %偏心距
k=2.5e7;%弹性轴刚度
delta2=0.6e-3;%初始松动间隙
c1=1050;%转子在轴承处阻尼系数
c2=2100;%转子圆盘处阻尼系数
r0=1.5e-6;%轴承间隙
Cb=13.34e9;%接触刚度;
k1=7.5e7;
k2=2.5e9;
cb1=350;
cb2=500;
w=2100;
ox1=y(1);%未松动端水平方向位移x1
oy1=y(2);%未松动端竖直方向位移y1
odx1=y(8);
ody1=y(9);
ox2=y(3);%圆盘位移x2
oy2=y(4);%圆盘位移y2
odx2=y(10);
ody2=y(11);
ox3=y(5);%松动端轴心位移x3
oy3=y(6);%松动端轴心位移y3
odx3=y(12);
ody3=y(13);
oy4=y(7);%质量m3在竖直方向位移y4
ody4=y(14);

%根据间隙大小判断
if oy4<0
    cb=cb2;
    kb=k2;
elseif  (oy4>=0)&(oy4<=delta2)
    cb=0;
    kb=0;
else
    cb=cb1;
    kb=k1;
end
M=[  m1   0     0    0    0     0    0;
    0    m1    0    0    0     0    0;
    0    0     m2   0    0     0    0;
    0    0     0    m2   0     0    0;
    0    0     0    0    m1    0    0;
    0    0     0    0    0     m1   0;
    0    0     0    0    0     0    m3;];
C=[  c1   0     0    0    0     0    0;
    0    c1    0    0    0     0    0;
    0    0     c2   0    0     0    0;
    0    0     0    c2   0     0    0;
    0    0     0    0    c1    0    0;
    0    0     0    0    0     c1   0;
    0    0     0    0    0     0    cb;];
K=[  k   0     -k    0    0     0    0;
    0    k    0    -k    0     0    0;
    -k   0    2*k    0    -k    0    0;
    0    -k   0    2*k    0    -k    0;
    0    0    -k    0     k    0    0;
    0    0    0    -k     0    k    0;
    0    0    0     0     0    0    kb;];

%未松动端轴承力fx
[Fx,Fy]=xufuze(t,ox1,oy1,w);
%松动端的
[Fx1,Fy1]=xufuze(t,ox3,oy3,w);
F=[   Fx;
    Fy-m1*g;
    m2*e*w^2*cos(w*t);
    m2*e*w^2*sin(w*t)-m2*g;
    Fx1;
    Fy1-m1*g;
    -Fy1-m3*g ];

%无量纲化
ox1=ox1/r0;
oy1=oy1/r0;
ox3=ox3/r0;
oy3=oy3/r0;
w0=((Cb*r0^1/2)/m1)^1/2;
t=w0*t;
oumiga=w/w0;
C=C/w0;%可能需要修改下，‘支承松动的转子-滚动轴承系统的分岔与混沌’
K=K/w0^2;
F=F/delta2/w0^2;

for i=1:1:N/2
    y1(i,1)=y(i);
    y2(i,1)=y(i+N/2);
end
yy2=inv(M)*(F-C*y2-K*y1);
d=zeros(N,1);
for i=1:1:N/2
    d(i)=y2(i,1);
    d(i+N/2)=yy2(i,1);
end
function [Fx,Fy]=xufuze(t,ox1,oy1,w)
r0=1.5e-6;%轴承间隙
Cb=13.34e9;%接触刚度;
R=63.9e-3;
r=40.1e-3;
Nb=8;
x=ox1;
y=oy1;
Fx=0;
Fy=0;
for j=1:1:Nb
%     H=0;
    BN=3.08;
    %角度计算
    w_cage=w*r/(R+r);
    aj=w_cage*t+2*pi/Nb*(j-1);
    %法向接触变形量
    deltaj=x*cos(aj)+y*sin(aj)-r0;
    if(deltaj>0)
        H=1;
    else
        H=0;
    end
    %接触压力计算
    fjx=Cb*((deltaj)^(3/2))*H*cos(aj);
    fjy=Cb*((deltaj)^(3/2))*H*sin(aj);

    Fx=Fx+fjx;
    Fy=Fy+fjy;
end

%主分析程序
clear all
clc
h=pi/256;
w=2100;
tf=300000*2*pi/w;
tspan=0:h:tf;
y0=[0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5];
options=odeset('RelTol',10^-6,'AbsTol',10^-6);  
[t,y]=ode45(@funa,tspan,y0);
figure(1)
subplot(2,2,1);
plot(t,y(:,1),'r',t,y(:,2),'b')
%plot(y(:,1),'r')
title('未松动端位移响应');xlabel('t');ylabel('x1-red,y1-blue');
subplot(2,2,2);
plot(t,y(:,3),'r',t,y(:,4),'b')
% plot(y(:,3),'r')
title('圆盘位移响应');xlabel('t');ylabel('x2-red,y2-blue');
subplot(2,2,3);
plot(t,y(:,5),'r',t,y(:,6),'b')
%plot(y(:,5),'r')
title('松动端轴心位移响应');xlabel('t');ylabel('x3-red,y4-blue');
subplot(2,2,4);
plot(y(:,7),'b')
title('m3在竖直方向位移响应');xlabel('t');ylabel('y4');
figure(2)
subplot(2,2,1);
plot(y(:,1),y(:,2))
title('未松动端轴心轨迹');xlabel('x1');ylabel('y1');
subplot(2,2,2);
plot(y(:,5),y(:,6))
title('松动端轴心轨迹');xlabel('x3');ylabel('y3');
yy1=sqrt(y(:,1).^2+y(:,2).^2);%yy1=sqrt(x1^2+y1^2)
yy2=sqrt(y(:,5).^2+y(:,6).^2);%yy2=sqrt(x3^2+y3^2)
N1=256;fs=1024;stepf=fs/N1;k=pi*2/882.5056;%882.5056转子固有频率sqrt(k/m2)
n=0:stepf*k:(fs/2-stepf)*k;
YY1=abs(fft(yy1));
YY2=abs(fft(yy2));
figure(3)
plot(n,abs(YY1(1:N1/2)));grid on;
title('未松动情况下幅值谱');xlabel('频率比');ylabel('FFT幅值');
figure(4)
plot(n,abs(YY2(1:N1/2)));grid on;
title('松动情况下幅值谱');xlabel('频率比');ylabel('FFT幅值');

fclose('all');

iwangtofly

个人对程序不太懂，但支持一下，还请高手指导~~
我不知道楼主这个程序是不是用matlab编写的，但我用matlab运行了一下，出现了下面的错误：

function d=funa(t,y,w)貌似定义的有问题~~

xufuze12

谢谢啦~！主程序在下面~~
%主分析程序
clear all
clc
h=pi/256;

xufuze12

祝各位新年快乐~！年前没弄好的程序，现在继续~不知哪位高手能不能帮忙修改下~！谢谢啦~！

东门吸雨

主程序显示错误是，Index exceeds matrix dimensions
YY1=abs(fft(yy1)); 是一个数值，只有一项
abs(YY1(1:N1/2)) 要调用 YY1的第1项到n/2项，所以报错了
我不懂机械，所以只能从编程的角度分析下。

xufuze12

东门吸雨 发表于 2013-2-19 07:24

                               
登录/注册后可看大图

主程序显示错误是，Index exceeds matrix dimensions
YY1=abs(fft(yy1)); 是一个数值，只有一项
abs(YY1( ...

非常感谢您的回复，这个地方是个问题，主要是我现在求的这个力太大了好像没收敛，如果力小一点就好了~我试着改了下 “Cb=13.34e9;%接触刚度;”把这个改小一点就不会有楼主所说的那种情况了，但是这个程序主要不是那里的问题，是我如何让这个力收敛，也就是在一定的范围内周期的变化~而不是一直增大

璃珞疏疏

请问楼主程序改好了吗？可以共享一下吗？毕业急需，谢谢