碰摩转子故障仿真程序

zhaowu

该程序在不加碰摩力时，得出结果很好。加入碰摩力，积分结果发散，得
不出仿真结果。

不平衡量在m2处，碰摩也发生在m2处。
  
仿真程序如下，matlab程序，newmark积分。
clear all
clc
%%----------------------------------
M=[3.5, 0,  0,   0,   0,   0                  
   0,   3.5,0,   0,   0,   0
   0,   0,  3.5, 0,   0,   0
   0,   0,  0,   3.5, 0,   0
   0,   0,  0,   0,   3.5, 0
   0,   0,  0,   0,   0,   3.5];
C=[1.05e3, 0,  0,   0,   0,   0
   0,   1.2e3,0,   0,   0,   0
   0,   0,  1.05e3, 0,   0,   0
   0,   0,  0,   1.05e3, 0,   0
   0,   0,  0,   0,   1.2e3, 0
   0,   0,  0,   0,   0,   1.05e3];
K=[3.1543e7, -3.0173e7,  1.2342e7,   0,   0,   0
   -3.0173e7,   4.3885e7,-3.0173e7,   0,   0,   0
   1.2342e7,   -3.0173e7,  3.1543e7, 0,   0,   0
   0,   0,  0,   3.1543e7, -3.0173e7,  1.2342e7
   0,   0,  0,   -3.0173e7,   4.3885e7,-3.0173e7
   0,   0,  0,   1.2342e7,   -3.0173e7,  3.1543e7];
%%------------------------------------            
u=zeros(6,1);
v=zeros(6,1);
x(:,1)=u;                                     %位移                       
x1(:,1)=v;                                    %速度
ee=0.13e-3;                                   %偏心距
h=0.1e-3;;                                    %动静间隙
kk=0.5e9;                                     %碰摩时接触刚度
miur=0.2;                                     %摩擦系数
delta0=0.0000001;                             %不对中
RR0=0.1;
delt=0.5;
alpha=0.25;
dt=2*pi/(1024*8);                            %积分步长
b0=1/(alpha*dt^2);
b1=delt/(alpha*dt);
b2=1/(alpha*dt);
b3=(1/(2*alpha))-1;
b4=(delt/alpha)-1;
b5=(dt/2)*((delt/alpha)-2);
b6=dt*(1-delt);
b7=delt*dt;      
w=0;
t_max=4;                 
ii=1;
t(1)=0;                 
f=86;
q=zeros(6,1);
w=0;
while t(ii)<t_max
     w=2*pi*f;                                %角速度
     fux(ii)=5.5*ee*(w*w*cos(w*t(ii)));       %不平衡力
     fuy(ii)=5.5*ee*(w*w*sin(w*t(ii)));       %不平衡力
     %----------------------碰摩---------------------------
     rrr(ii)=sqrt((x(2,ii))^2+(x(5,ii)-delta0)^2);
     ggg(ii)=(rrr(ii)-h);
     fai(ii)=((x(1,ii)*x1(5,ii)-(x(5,ii)-delta0)*x1(2,ii))/rrr(ii))+RR0*w;
       if fai(ii)>0
          thetar=1;
       elseif fai(ii)<0
          thetar=-1;
       else
          thetar=0;
       end
       if (ggg(ii))>0
        N(ii)=(kk*ggg(ii)/rrr(ii));
        frx(ii)=(-N(ii))*((x(2,ii))-miur*((x(5,ii))-delta0));
        fry(ii)=(-N(ii))*((miur*x(2,ii))+(x(5,ii)-delta0));
       else
        fry(ii)=0;
        frx(ii)=0;
       end
       Q2(:,ii)=zeros(6,1);
       Q2(2,ii)=Q2(2,ii)+frx(ii);                                            
       Q2(5,ii)=Q2(5,ii)+fry(ii);
       %------------------------------------------------------------
     Q(:,ii)=zeros(6,1);
     Q1(:,ii)=zeros(6,1);
     Q1(2,ii)=Q1(2,ii)+fux(ii);
     Q1(5,ii)=Q1(5,ii)+fuy(ii);
     Q(:,ii)=Q1(:,ii);
     Ccc=C;
     K1=K+b0*M+b1*Ccc ;
     if ii==1
        x2(:,ii)=inv(M)*(Q(:,ii)-K*x(:,ii)-Ccc*x1(:,ii));
        q(:,ii)=Q(:,ii)+M*(b0*x(:,ii)+b2*x1(:,ii)+b3*x2(:,ii))+Ccc*(b1*x(:,ii)+b4*x1(:,ii)+b5*x2(:,ii));
        x(:,ii+1)=K1\q(:,ii);
        x2(:,ii+1)=b0*(x(:,ii+1)-x(:,ii))-b2*x1(:,ii)-b3*x2(:,ii);
        x1(:,ii+1)=x1(:,ii)+b6*x2(:,ii)+b7*x2(:,ii+1);
      else
        q(:,ii+1)=Q(:,ii)+M*(b0*x(:,ii)+b2*x1(:,ii)+b3*x2(:,ii))+Ccc*(b1*x(:,ii)+b4*x1(:,ii)+b5*x2(:,ii));
        x(:,ii+1)=K1\q(:,ii+1);
        x2(:,ii+1)=b0*(x(:,ii+1)-x(:,ii))-b2*x1(:,ii)-b3*x2(:,ii);
        x1(:,ii+1)=x1(:,ii)+b6*x2(:,ii)+b7*x2(:,ii+1);            
      end
       ii=ii+1;
       t(ii)=t(ii-1)+dt;
end
%%----------------------------------------
T=1/dt;
TT=T*t_max;
t=3000:TT;
subplot(2,2,(1:2))
plot(t/T,x(2:2,t)*1000)   
title('位移-1x')
grid on
subplot(2,2,3)
Nx=length(x(2:2,t));
Y=fft(x(2:2,t)*T,Nx);
pyy=Y.*conj(Y)/(1024*8);
f=T*(0:Nx/2-1)/Nx;
plot(f,pyy(1:(Nx/2)))  
title('幅值谱-节点1x')
grid on
subplot(224)
plot(x(2:2,t)*1000,x(5:5,t)*1000)
title('轴心轨迹')
grid on

zhaowu

得出结论了
上面的程序是对的,只是积分步长不够小,改小就ok

qijunshuai

你好，我想问一下，是不是用有限元方法列出的转子系统方程后，直接把几个系数矩阵带入就可以了？可是我建立一个双圆盘转子的运动微分方程，为什么结出的结果不正确？ 难道我求质量矩阵，刚度矩阵，阻尼矩阵出现了问题？
这位大侠能否留下联系方式，请教下问题

zhaowu

试着积分步长取小些
qq:120852933

songshuomjy

原帖由 zhaowu 于 2007-9-18 20:01 发表
试着积分步长取小些
qq:120852933

..................

yejet

原帖由 qijunshuai 于 2007-9-13 10:32 发表
你好，我想问一下，是不是用有限元方法列出的转子系统方程后，直接把几个系数矩阵带入就可以了？可是我建立一个双圆盘转子的运动微分方程，为什么结出的结果不正确？ 难道我求质量矩阵，刚度矩阵，阻尼矩阵出现 ...

直接把几个系数矩阵带入方程，形成一个二阶微分方程组，然后求解就可以了

chaijize

请问楼主,这个程序碰摩力加在那里?

zhaowu

Q2(2,ii)=Q2(2,ii)+frx(ii);                                            
       Q2(5,ii)=Q2(5,ii)+fry(ii);
       frx，fry表示x,y方向碰摩力

0810064

谢谢楼主正学习中好东东当然要收藏了

菜肴全

感谢楼主好东西哈