完整的支承松动数值仿真程序
本帖最后由 牛小贱 于 2014-3-27 21:33 编辑%运动微分方程
function d=fun(t,y,w)
N=length(y);
w=2100;
m1=4;%两端滑动轴承处等效集中质量
m2=32.1; %转子圆盘等效集中质量
m3=50.0;%轴承支座处等效集中质量
g=9.81;
e=0.00005; %偏心距
k=2.5e7;%弹性轴刚度
delta2=0.6e-3;%初始间隙
c1=1050;%转子圆盘处阻尼系数
c2=2100;%转子在轴承处阻尼系数
k1=7.5e7;
k2=2.5e9;
cb1=350;
cb2=500;
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=oilx( ox1, oy1, odx1, ody1, w);
fy=oily( ox1, oy1, odx1, ody1, w);
fx1=oilx( ox3,oy3-oy4,odx3,ody3-ody4,w);
fy1=oily( ox3,oy3-oy4,odx3,ody3-ody4,w);
F=[ fx;
fy-m1*g;
m2*e*w^2*cos(t);
m2*e*w^2*sin(t)-m2*g;
fx1;
fy1-m1*g;
-fy1-m3*g ];
C=C/w;
K=K/w^2;
F=F/c/w^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
%x方向油膜力
function oilforce=oilx(x,y,dx,dy,wi)
R=0.025; L=0.012; miu=0.018;dfai=1; dert=0.00011;
e=sqrt(x*x+y*y);
delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
ppp1=(y+2.0*dx)/(x-2.0*dy);
sign1=sign(ppp1);
ppp2=y+2.0*dx;
sign2=sign(ppp2);
alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));
f1=sqrt(abs(abs((x-2.0*dy)*(x-2.0*dy))+abs((y+2.0*dx)*(y+2.0*dx))))/(1.0-abs(x*x)-abs(y*y));
fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);
oilforce=fx*delta;
%y方向油膜力
function oilforce=oily(x,y,dx,dy,wi)
R=0.025; L=0.012; miu=0.018;dfai=1; dert=0.00011;
e=sqrt(x*x+y*y);
delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
ppp1=(y+2.0*dx)/(x-2.0*dy);
sign1=sign(ppp1);
ppp2=y+2.0*dx;
sign2=sign(ppp2);
alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));
f1=sqrt(abs(abs((x-2.0*dy)*(x-2.0*dy))+abs((y+2.0*dx)*(y+2.0*dx))))/(1.0-abs(x*x)-abs(y*y));
fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);
oilforce=fy*delta;
%主分析程序
clear all
h=pi/256;
w=2100;
tf=300000*2*pi/w;
tspan=0:h:tf;
y0=;
options=odeset('RelTol',10^-6,'AbsTol',10^-6);
=ode45(@fun,tspan,y0);
figure(1)
subplot(2,2,1);
plot(t,y(:,1),'r',t,y(:,2),'b')
title('未松动端位移响应');xlabel('t');ylabel('x1-red,y1-blue');
subplot(2,2,2);
plot(t,y(:,3),'r',t,y(:,4),'b')
title('圆盘位移响应');xlabel('t');ylabel('x2-red,y2-blue');
subplot(2,2,3);
plot(t,y(:,5),'r',t,y(:,6),'b')
title('松动端轴心位移响应');xlabel('t');ylabel('x3-red,y4-blue');
subplot(2,2,4);
plot(t,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');
主要参考文献：李振平，罗跃纲等.转子系统支承松动的非线性动力学及故障特征.东北大学学报(自然科学版)，2002,23(11):1049-1051.
收藏了,谢谢楼主! 发现错误：
其中
fx1=oilx( ox3,oy3-oy4,odx3,odx3-ody4,w);
fy1=oily( ox3,oy3-oy4,odx3,odx3-ody4,w);
应该改为
fx1=oilx( ox3,oy3-oy4,odx3,ody3-ody4,w);
fy1=oily( ox3,oy3-oy4,odx3,ody3-ody4,w);
其他应该没什么问题了，仿真结果跟原论文差不多。
[ 本帖最后由 pwangeng 于 2007-5-25 12:21 编辑 ] 已经帮你在原帖修正
[ 本帖最后由 ChaChing 于 2010-3-1 10:13 编辑 ]
感触下先
好东西,学习的时候就是少了点程序,看理论还是有点远:lol运行出现错误
??? Undefined function or variable 'c'.Error in ==> fun at 92
F=F/c/w^2;
Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in ==> main_songdong at 11
=ode45(@fun,tspan,y0); 有没有相应的理论方法介绍？ 原帖由 huazaiwang 于 2007-10-21 13:38 发表
??? Undefined function or variable 'c'.
Error in ==> fun at 92
F=F/c/w^2;
Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
...
这的确是个问题，除非就是C，否则恐怕要楼主现身才能修正
回复 11楼 的帖子
程序里面的c就是dert=0.00011,也就是油膜力的平均厚度。个人感觉程序有点差异。C=C/w;
K=K/w^2;
F=F/c/w^2;
应该是对方程进行无量纲化，既然无量纲化了系统的周期也就变成了2pi呀！望和楼主讨论:@)
还有，前期的运动一般都不稳定，要去掉瞬态解。
刚接触非线性不久，说的不对的地方请楼主包涵:loveliness: {:{23}:}我太需要这个了！谢谢楼主！ 本帖最后由 mjtjiang 于 2010-9-28 19:11 编辑
pwangeng 发表于 2007-5-24 15:20 http://www.chinavib.com/static/image/common/back.gif
alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));
f1=sqrt(abs(abs((x-2.0*dy)*(x-2.0*dy))+abs((y+2.0*dx)*(y+2.0*dx))))/(1.0-abs(x*x)-abs(y*y));
其中用abs是什么意思呢？若不用abs会出现复数，其中的虚部有什么意义吗？ mjtjiang 发表于 2010-9-28 19:09 static/image/common/back.gif
其中用abs是什么意思呢？若不用abs会出现复数，其中的虚部有什么意义吗？
abs应该是绝对值的意思
我好像见过这个模型，应该是有这个绝对值得 回复 雪缘 的帖子
虽然还不是很懂，但先记下吧。
求这个非线性油膜力的时候，为什么前面要乘以-0.1啊，这是什么意义啊，书上没这么写啊。
fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);
还有
C=C/w;
K=K/w^2;
F=F/c/w^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
mjtjiang 发表于 2010-9-29 09:12 static/image/common/back.gif
回复 雪缘 的帖子
虽然还不是很懂，但先记下吧。
我查了一下，这是1991年法国人Capone提出来的短圆柱瓦轴承油膜力模型，建议你找这个模型对照一下就都清楚了
问题1：不是 -0.1 是 -1.0 ,这表示力的方向问题；
问题2：是无量纲化处理，两者的解理论上应该是完全等价的，但是实际上有可能由于计算过程中的截断误差而导致结果不一样，无量纲化目的有两个，一是避免引入过大的截断误差，二是为了定性分析；
问题3：这一部分是将二阶微分方程转换为一阶状态方程，得到ode45所能求解的微分方程标准格式。 本帖最后由 mjtjiang 于 2010-9-29 16:40 编辑
你说的那个原文 Capone G. Descrizion analitica del campo di forze fluidodinamiconei cuscinetti cilindrici cubrificati.L’EnergiaElettrica,1991,3(1):105-110. 我没找到，但我有这个表达式：
我不是学转子动力学的，其原理我都不懂，但老师要我做的课题的一部分是转子故障信号仿真，我就只是按照书上给的公式编程序，发现有好多问题
问题一：楼主的这个程序是在程序里将二阶微分方程转换为一阶状态方程，我是在纸上手动计算把他降阶，在程序里直接编入的就是无量纲性形式的一阶状态方程，我这样做可以吗？反正程序能运行，但输出总是与书上给出的不符，我算出的轴心轨迹是像正弦信号一样震荡的。我也不知道自己错哪儿了。
问题二：楼主的
tf=300000*2*pi/w;和
stepf=fs/N1;k=pi*2/882.5056;%882.5056转子固有频率sqrt(k/m2)
n=0:stepf*k:(fs/2-stepf)*k;这两句我不太懂，还请校长指教。