求助 裂纹转子系统故障信号matlab仿真

mjtjiang

请大家帮我看看这段程序吧，原方程和降阶后的一阶方程我都贴出来了，程序也在下面，方程中的参数都是按书上给的输入的，可仿真结果跟书上差很多啊，不知道是我的方程编写错误还是求频谱的傅里叶变换程序部分出错了，希望大家指点啊！

无量纲参数：


裂纹转子系统运动微分方程的无量纲形式：



因为是要用ode45求解方程，所以将方程降阶：
令



为编程方便，又把系数简化了一下：



程序贴在楼下

mjtjiang

回复 1 # mjtjiang 的帖子

方程在楼上贴出了，以下是程序，就算是截取后半段稳定值作分析，运行结果也还是与书上的不相符，请大家帮忙指点

function out=liewenfun(t,z,flag)

omega=2000;           % 轴的转速
a=0.0135;                % 裂纹深度
r=0.015;                  % r转轴半径
L1=0.57;                 % L1转轴长度
R=0.025;                 % R轴承半径
L=0.012;                 % L轴承长度
b=0.05*10^(-3);    % b圆盘偏心率
m1=4.0;                 % m1轴承处的等效集中质
m2=32.1;               % m2转轴中央圆盘质量
c=0.11*10^(-3);    % c轴承间隙
mu=0.018;             % mu润滑油粘度
c1=1050.0;             % c1转子在轴承处的结构阻尼
c2=2100.0;             % c2转子在圆盘处的结构阻尼
k=2.5*10^7;           % k转轴刚度
g=9.81;
epsilon=3.3064e-007;
delta=3.0244e+006;     % epsilon与delta为仅与裂纹深度a有关的相关刚度参数
P=0.5*m2;             % P转子圆盘重量的一半

beta=0;        % beta裂纹方向与偏心之间的夹角
phi=0;          % phi转轴初始相位角
theta=omega*t+phi+beta; % theta裂纹截面转角
F=(1+cos(theta))/2;     % F裂纹开闭规律，描述的是裂纹转角对实际裂纹轴刚度的影响规律(裂纹开闭规律模型有好几种，我用的是比较简单的余弦波模型)

s=(mu*omega*R*L*((R/c)^2)*(L/(2*R))^2)/P; % s为Sommerfeld修正系数，c平均油膜厚度
xi1=c1/(m1*omega);
xi2=c2/(m2*omega);
eta1=k/(m1*omega^2);
eta2=2*k/(m2*omega^2);
G=g/(c*omega^2);
M1=(c*m1*omega^2)/(s*P);
b1=b/c;
tau=omega*t;
w=1-epsilon*delta*0.5*F;
u=epsilon*0.5*F.*cos(2*tau);
l=epsilon*0.5*F.*sin(2*tau); %时变系数

out=[z(5);...
z(6);...
z(7);...
z(8);...
(-eta1*w-eta1*u)*z(1)-eta1*l*z(2)+(eta1*w+eta1*u)*z(3)+eta1*l*z(4)-xi1*z(5)+fx(z(1),z(2),z(5),z(6))/M1;...
-eta1*l*z(1)+(-eta1*w+eta1*u)*z(2)+eta1*l*z(3)+(eta1*w-eta1*u)*z(4)-xi1*z(6)+fy(z(1),z(2),z(5),z(6))/M1-G;...
(eta2*w+eta2*u)*z(1)+eta2*l*z(2)+(-eta2*w-eta2*u)*z(3)-eta2*l*z(4)-xi2*z(7)+b1*cos(tau);...
eta2*l*z(1)+(eta2*w-eta2*u)*z(2)-eta2*l*z(3)+(-eta2*w+eta2*u)*z(4)-xi2*z(8)+b1*sin(tau)-G];

fx、fy是无量纲非线性油膜力

function fxout=fx(x,y,Dx,Dy)
if x-2*Dy~=0
alpha=atan((y+2*Dx)./(x-2*Dy))-0.5*pi*sign((y+2*Dx)./(x-2*Dy))-0.5*pi*sign(y+2*Dx);
else
alpha=atan2((y+2*Dx),(x-2*Dy))-0.5*pi*sign(y+2*Dx)-0.5*pi*sign(y+2*Dx);
end
G=(2.0/sqrt(abs(1.0-abs(x.^2)-abs(y.^2)))).*(pi/2+atan((y.*cos(alpha)-x.*sin(alpha))./sqrt(abs(1.0-abs(x.^2)-abs(y.^2)))));
V=(2.0+(y.*cos(alpha)-x.*sin(alpha)).*G)./(1.0-abs(x.^2)-abs(y.^2));
S=(x.*cos(alpha)+y.*sin(alpha))./(1.0-abs((x.*cos(alpha)+y.*sin(alpha)).^2));
fxout=-1.0*(sqrt(abs(abs((x-2*Dy).^2)+abs((y+2*Dx).^2)))./(1.0-abs(x.^2)-abs(y.^2))).*(3.0.*x.*V-sin(alpha).*G-2.0.*cos(alpha).*S);

function fyout=fy(x,y,Dx,Dy)
if x-2*Dy~=0
alpha=atan((y+2*Dx)./(x-2*Dy))-0.5*pi*sign((y+2*Dx)./(x-2*Dy))-0.5*pi*sign(y+2*Dx);
else
alpha=atan2((y+2*Dx),(x-2*Dy))-0.5*pi*sign(y+2*Dx)-0.5*pi*sign(y+2*Dx);
end
G=(2.0/sqrt(abs(1.0-abs(x.^2)-abs(y.^2)))).*(pi/2+atan((y.*cos(alpha)-x.*sin(alpha))./sqrt(abs(1.0-abs(x.^2)-abs(y.^2)))));
V=(2.0+(y.*cos(alpha)-x.*sin(alpha)).*G)./(1.0-abs(x.^2)-abs(y.^2));
S=(x.*cos(alpha)+y.*sin(alpha))./(1.0-abs((x.*cos(alpha)+y.*sin(alpha)).^2));
fyout=-1.0*(sqrt(abs(abs((x-2*Dy).^2)+abs((y+2*Dx).^2)))./(1.0-abs(x.^2)-abs(y.^2))).*(3.0.*y.*V+cos(alpha).*G-2.0.*sin(alpha).*S);

% 主程序
clear all

h=pi/256;
w=2000;
tf=100000*2*pi/w;
tspan=0:h:tf;
z0=[0 0 0 0 0 0 0 0]';
[t,z]=ode45('liewenfun',tspan,z0);

figure(1)
subplot(2,2,1);
plot(t,z(:,1))
title('裂纹转子左端轴心径向位移x');
xlabel('t');ylabel('x1');
subplot(2,2,2);
plot(t,z(:,2))
title('裂纹转子左端轴心径向位移y');
xlabel('t');ylabel('y1');
subplot(2,2,3);
plot(z(:,1),z(:,2))
title('裂纹转子左端轴心轨迹');xlabel('x1');ylabel('y1');
subplot(2,2,4);
plot(z(:,1),z(:,5))
title('裂纹转子左端轴心相轨迹');xlabel('x1');ylabel('Dx1');

N1=length(z(:,1));fs=1024;stepf=fs/N1;
n=0:stepf:(fs/2-stepf);
Z1=abs(fft(z(:,1)));
figure(2)
plot(n(1:end),abs(Z1(1:N1/2)));grid on;
axis([-1 10 0 10000]);
title('转子左端轴心位移幅值谱');xlabel('频率');ylabel('FFT幅值');

happy

估计别人帮你查程序的可能性比较小
把书上的结果和你计算的结果都贴出来看看，是否能从结果上出什么端倪

mjtjiang

回复 3 # happy 的帖子

谢谢happy教授提醒

我查的一篇论文里给出的ω=1540 rad/s、A=0.5时的响应形式，和ω=1800 rad/s、A=0.9时的响应形式分别如下图
从上到下分别为时序波形，相轨迹和功率谱图：





而我做出的图形是这样的：
ω=1540 rad/s、A=0.5时的响应形式：
时序波形和相轨迹：


频谱图截图：


ω=1800 rad/s、A=0.9时的响应形式：
时序波形和相轨迹：


频谱图截图：


真的是弄不清楚是怎么回事了，在论坛里把所有关于转子的帖子都看过了，信号处理也看了不少，可就是弄不明白了，只有请大家帮助了！

happy

一个很明显的问题就是你的计算结果中瞬态量没有消除

个人感觉你对这类系统的仿真模拟还没有入门
建议一开始不要做这样复杂的模型

mjtjiang

回复 5 # happy 的帖子

请问happy教授，怎样才能入门啊？需要学什么呢？我现在是被老师赶鸭子上架做这个仿真，之前一点基础都没有的，请happy教授指点！

happy

mjtjiang 发表于 2010-12-4 09:56

                               
登录/注册后可看大图

回复 5 # happy 的帖子

请问happy教授，怎样才能入门啊？需要学什么呢？我现在是被老师赶鸭子上架做这个仿 ...

多看多练

woshiaq

我想问一下图中的A=0.9是什么参数呢   在你的公式中也没看到这个参数

mjtjiang

回复 8 # woshiaq 的帖子

A是无量纲裂纹深度，是实际裂纹深度与转轴半径的比值

0810064

我运行了楼主的程序 ，总是提示有错误，甚至连楼主的曲线都没得出，正学习中，研究两天了也没修正好程序

ChaChing

回复 10 # 0810064 的帖子

:@)求助完整格式：出错代码和出错提示

0810064

谢谢楼上的朋友。还是在学习中

0810064

总算得出和楼主一样的仿真图了。尽管不理想，总算是对两天的工作有一个交待，昨天推了半天的大矩阵。谢谢Chaching前辈的鼓励。刚入学，正处在学习的初级阶段。

tqc700322

学习学习，

伤痕累累

请问你最后做出来了吗？