裂纹转子振动仿真

zhlong

最近看了一些裂纹转子故障诊断方面的论文，里面很多都利用到裂纹转子的 动力学模型，以产生裂纹转子振动仿真信号，然后再利用适当的信号处理方法如小波分析等从仿真信号中提取特征。动力学方程是一个两个二阶微分方程组成的方程组，但是其中有一个参数即裂纹模型f是随自变量变化的，即导致微分方程成为变参数微分方程。

附件的WORD文档是我主要根据文献《小波变换在横向裂纹转子升速过程状态监测中的应用》（作者：钱立军，等。中国电机工程学报）列出的裂纹转子动力学方程

下面是我编写的matlab程序，但是得不到正确的结果，请各位帮忙看看
=======================================================================
===================主程序================================================

1. %赋初值  这里的处置是随意赋予的，可能有些不符合常理
   
2. a=11;%加速度
   
3. m=1;%转子质量
   
4. k=80;%支撑刚度
   
5. c=7;%阻尼
   
6. bita=0.5;%偏心质量与裂纹法向夹角
   
7. k1=0.07;%裂纹导致的刚度变化量
   
8. e=0.00001;%偏心距
   
9. alfa=pi/6;%裂纹张角的一半
   
10. tspan=[0.0001,2*5*pi];%采样转子转动5周的数据
    
11. y0=[0,0,0,0];              %给定初值
    
12. options=odeset; options.RelTol=1e-6;
    
13. [xx,XX]=ode45('vib',tspan,y0,[],a,m,k,c,bita,k1,e,alfa);

复制代码

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝‘
＝＝＝＝＝＝＝振动微分方程＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝＝二阶微分方程降阶为一阶＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

1. function dX=vib(xita,X,flag,a,m,k,c,bita,k1,e,alfa)
   
2. %运动微分方程   将二阶微分方程降阶为一阶
   
3. %a——转子旋转角加速度；m——转子质量；k——转子支撑刚度；c——转子支撑阻尼；
   
4. %xita   偏心质量与旋转坐标原点O’连线与固定坐标系x轴之间的夹角===>用于表示转子在固定坐标中的角位置
   
5. %bita   偏心质量与旋转坐标原点O’连线与旋转坐标系e轴的夹角（e轴将裂纹垂直平分为两半），这个角度一旦设定之后在转子运行过程中是不变的
   
6. %xita+bita 表示裂纹中心与固定坐标系x轴之间的夹角
   
7. %k1 裂纹导致的刚度下降比例  （刚度变化量/无裂纹时的刚度）
   
8. %e  偏心质量与旋转坐标系O'
   
9. t=abs(sqrt(2*xita/a)); w=a*t; wn=sqrt(k/m); g=9.8; Xs=m*g/k; omiga=w/wn;ibuxi=c/(2*sqrt(m*k)); ee=e/Xs;
   
10. %t——时间；w——角速度；wn——固有频率；g——重力加速度；Xs——重力导致的转子变形；
    
11. 
    
12. dX=[X(2);
    
13.     -2*ibuxi/omiga*X(2)-X(1)/omiga^2+fcrack(xita,bita,alfa)*k1/(2*omiga^2)*((1+cos(2*(xita+bita)))*X(1)+sin(2*(xita+bita))...
    
14.     *X(3))+1/omiga^2+ee/w^2*(a*sin(xita)+w^2*cos(xita));
    
15.    X(4);
    
16.    -2*ibuxi/omiga*X(4)-X(3)/omiga^2+fcrack(xita,bita,alfa)*k1/(2*omiga^2)*((1-cos(2*(xita+bita)))*X(3)+sin(2*(xita+bita))...
    
17.    *X(1))+ee/w^2*(-a*cos(xita)+w^2*sin(xita))];

复制代码

=====================================================================
==================裂纹开闭函数==========================================
=================用的是综合模型==========================================

1. %裂纹开闭模型函数
   
2. function f=fcrack(xita,bita,alfa)
   
3. %采用综合模型描述裂纹的开闭
   
4. %xita   偏心质量与旋转坐标原点O’连线与固定坐标系x轴之间的夹角===>用于表示转子在固定坐标中的角位置
   
5. %bita   偏心质量与旋转坐标原点O’连线与旋转坐标系e轴的夹角（e轴将裂纹垂直平分为两半），这个角度一旦设定之后在转子运行过程中是不变的
   
6. %alfa   裂纹所占角度的一半
   
7. xita=xita-floor(xita/(2*pi))*2*pi;
   
8. f=((0<=xita & xita<=pi/2-alfa-bita)|(3*pi/2+alfa-bita<=xita & xita<=2*pi))+...
   
9.     (pi/2-alfa-bita<xita & xita<=pi/2+alfa-bita)*(1/2+1/2*cos((xita+bita-pi/2+alfa)./(2*(xita+bita))))+...
   
10.     (pi/2+alfa-bita<xita & xita<=3*pi/2-alfa-bita)*(0)+...
    
11.     (3*pi/2-alfa-bita<=xita & xita<3*pi/2+alfa-bita)*(1/2+1/2*cos((xita+bita-3*pi/2+alfa)./(2*(xita+bita))));

复制代码

beyondbear

你的裂纹开闭模型函数公式写错了吧和《小波变换在横向裂纹转子升速过程状态监测中的应用》对不上