裂纹转子振动仿真

zhlong

最近看了一些裂纹转子故障诊断方面的论文，里面很多都利用到裂纹转子的 动力学模型，以产生裂纹转子振动仿真信号，然后再利用适当的信号处理方法如小波分析等从仿真信号中提取特征。

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

下面是我编写的matlab程序，但是得不到正确的结果，请各位帮忙看看
=======================================================================
===================主程序================================================
%赋初值  这里的处置是随意赋予的，可能有些不符合常理
a=11;%加速度
m=1;%转子质量
k=80;%支撑刚度
c=7;%阻尼
bita=0.5;%偏心质量与裂纹法向夹角
k1=0.07;%裂纹导致的刚度变化量
e=0.00001;%偏心距
alfa=pi/6;%裂纹张角的一半


tspan=[0.0001,2*5*pi];%采样转子转动5周的数据

y0=[0,0,0,0];              %给定初值

options=odeset; options.RelTol=1e-6;

[xx,XX]=ode45('vib',tspan,y0,[],a,m,k,c,bita,k1,e,alfa);

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝‘
＝＝＝＝＝＝＝振动微分方程＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝＝二阶微分方程降阶为一阶＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
function dX=vib(xita,X,flag,a,m,k,c,bita,k1,e,alfa)
%运动微分方程   将二阶微分方程降阶为一阶
%a——转子旋转角加速度；m——转子质量；k——转子支撑刚度；c——转子支撑阻尼；
%xita   偏心质量与旋转坐标原点O’连线与固定坐标系x轴之间的夹角===>用于表示转子在固定坐标中的角位置
%bita   偏心质量与旋转坐标原点O’连线与旋转坐标系e轴的夹角（e轴将裂纹垂直平分为两半），这个角度一旦设定之后在转子运行过程中是不变的
%xita+bita 表示裂纹中心与固定坐标系x轴之间的夹角
%k1 裂纹导致的刚度下降比例  （刚度变化量/无裂纹时的刚度）
%e  偏心质量与旋转坐标系O'
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;
%t——时间；w——角速度；wn——固有频率；g——重力加速度；Xs——重力导致的转子变形；

dX=[X(2);
    -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))...
    *X(3))+1/omiga^2+ee/w^2*(a*sin(xita)+w^2*cos(xita));
   X(4);
   -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))...
   *X(1))+ee/w^2*(-a*cos(xita)+w^2*sin(xita))];

=====================================================================
==================裂纹开闭函数==========================================
=================用的是综合模型==========================================
%裂纹开闭模型函数
function f=fcrack(xita,bita,alfa)
%采用综合模型描述裂纹的开闭
%xita   偏心质量与旋转坐标原点O’连线与固定坐标系x轴之间的夹角===>用于表示转子在固定坐标中的角位置
%bita   偏心质量与旋转坐标原点O’连线与旋转坐标系e轴的夹角（e轴将裂纹垂直平分为两半），这个角度一旦设定之后在转子运行过程中是不变的
%alfa   裂纹所占角度的一半

xita=xita-floor(xita/(2*pi))*2*pi;
if (0<=xita & xita<=pi/2-alfa-bita)|(3*pi/2+alfa-bita<=xita & xita<=2*pi)
    f=1;
elseif pi/2-alfa-bita<xita & xita<=pi/2+alfa-bita
    f=1/2+1/2*cos((xita+bita-pi/2+alfa)./(2*(xita+bita)));
elseif pi/2+alfa-bita<xita & xita<=3*pi/2-alfa-bita
    f=0;
elseif 3*pi/2-alfa-bita<=xita & xita<3*pi/2+alfa-bita
     f=1/2+1/2*cos((xita+bita-3*pi/2+alfa)./(2*(xita+bita)));   
end

[ 本帖最后由 zhlong 于 2006-12-15 16:09 编辑 ]

zhlong

微分方程

（附件已删除）

[ 本帖最后由 xinyuxf 于 2006-12-26 16:42 编辑 ]

zhlong

哪位做过裂纹转子仿真的帮忙看看好吗？
在你们也许看来很简单的问题，可是我已经花了好几个星期了也没做出来。

mriook

该去matlab版问吧

zhlong

谢谢mriook 的提醒，故障诊断版的高手啊
你是在工作还是学校的啊

zhlong

救命啊，都没人解答啊

mriook

我可不是高手，工作了
在这问是地方不大对，呵呵

zhlong

我现在学校读研，想做点正真有用的东西，但似乎不容易
根据你的工作经验，我们这些在校学生应该往哪些方面作点工作
才能产生实际效应啊。