查看: 8014|回复: 57

[共享资源] 完整的支承松动数值仿真程序

  [复制链接]
发表于 2007-5-24 15:20 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
本帖最后由 牛小贱 于 2014-3-27 21:33 编辑
  1. %运动微分方程
  2. function d=fun(t,y,w)
  3. N=length(y);
  4. w=2100;
  5. m1=4;%两端滑动轴承处等效集中质量
  6. m2=32.1; %转子圆盘等效集中质量
  7. m3=50.0;%轴承支座处等效集中质量
  8. g=9.81;
  9. e=0.00005; %偏心距
  10. k=2.5e7;%弹性轴刚度
  11. delta2=0.6e-3;%初始间隙
  12. c1=1050;%转子圆盘处阻尼系数
  13. c2=2100;%转子在轴承处阻尼系数
  14. k1=7.5e7;
  15. k2=2.5e9;
  16. cb1=350;
  17. cb2=500;

  18. ox1=y(1);%未松动端竖直方向位移x1
  19. oy1=y(2);%未松动端竖直方向位移y1
  20. odx1=y(8);
  21. ody1=y(9);

  22. ox2=y(3);%圆盘位移x2
  23. oy2=y(4);%圆盘位移y2
  24. odx2=y(10);
  25. ody2=y(11);

  26. ox3=y(5);%松动端轴心位移x3
  27. oy3=y(6);%松动端轴心位移y3
  28. odx3=y(12);
  29. ody3=y(13);

  30. oy4=y(7);%质量m3在竖直方向位移y4
  31. ody4=y(14);

  32. if oy4<0
  33.     cb=cb2;
  34.     kb=k2;
  35. elseif  (oy4>=0)&(oy4<=delta2)
  36.         cb=0;
  37.         kb=0;
  38.     else
  39.          cb=cb1;
  40.          kb=k1;
  41. end

  42. M=[  m1   0     0    0    0     0    0;
  43.      0    m1    0    0    0     0    0;
  44.      0    0     m2   0    0     0    0;
  45.      0    0     0    m2   0     0    0;
  46.      0    0     0    0    m1    0    0;
  47.      0    0     0    0    0     m1   0;
  48.      0    0     0    0    0     0    m3;];

  49. C=[  c1   0     0    0    0     0    0;
  50.      0    c1    0    0    0     0    0;
  51.      0    0     c2   0    0     0    0;
  52.      0    0     0    c2   0     0    0;
  53.      0    0     0    0    c1    0    0;
  54.      0    0     0    0    0     c1   0;
  55.      0    0     0    0    0     0    cb;];
  56.   
  57. K=[  k   0     -k    0    0     0    0;
  58.      0    k    0    -k    0     0    0;
  59.      -k   0    2*k    0    -k    0    0;
  60.      0    -k   0    2*k    0    -k    0;
  61.      0    0    -k    0     k    0    0;
  62.      0    0    0    -k     0    k    0;
  63.      0    0    0     0     0    0    kb;];

  64. fx=oilx( ox1, oy1, odx1, ody1, w);
  65. fy=oily( ox1, oy1, odx1, ody1, w);
  66. fx1=oilx( ox3,oy3-oy4,odx3,ody3-ody4,w);
  67. fy1=oily( ox3,oy3-oy4,odx3,ody3-ody4,w);



  68. F=[   fx;
  69.       fy-m1*g;
  70.       m2*e*w^2*cos(t);
  71.       m2*e*w^2*sin(t)-m2*g;
  72.       fx1;
  73.       fy1-m1*g;
  74.       -fy1-m3*g ];
  75.    
  76. C=C/w;
  77. K=K/w^2;
  78. F=F/c/w^2;

  79. for i=1:1:N/2
  80.     y1(i,1)=y(i);
  81.     y2(i,1)=y(i+N/2);
  82. end


  83. yy2=inv(M)*(F-C*y2-K*y1);

  84. d=zeros(N,1);

  85. for i=1:1:N/2
  86.     d(i)=y2(i,1);
  87.     d(i+N/2)=yy2(i,1);
  88. end


  89. %x方向油膜力

  90. function oilforce=oilx(x,y,dx,dy,wi)

  91. R=0.025; L=0.012; miu=0.018;  dfai=1; dert=0.00011;

  92. e=sqrt(x*x+y*y);
  93. delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
  94. ppp1=(y+2.0*dx)/(x-2.0*dy);
  95. sign1=sign(ppp1);
  96.         
  97. ppp2=y+2.0*dx;
  98. sign2=sign(ppp2);

  99. alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
  100. alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
  101. fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
  102. fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
  103. fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));

  104. 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));

  105. fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
  106. fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);

  107. oilforce=fx*delta;



  108. %y方向油膜力

  109. function oilforce=oily(x,y,dx,dy,wi)

  110. R=0.025; L=0.012; miu=0.018;  dfai=1; dert=0.00011;

  111. e=sqrt(x*x+y*y);
  112. delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
  113. ppp1=(y+2.0*dx)/(x-2.0*dy);
  114. sign1=sign(ppp1);
  115.       
  116. ppp2=y+2.0*dx;
  117. sign2=sign(ppp2);

  118. alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
  119. alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
  120. fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
  121. fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
  122. fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));

  123. 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));

  124. fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
  125. fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);

  126. oilforce=fy*delta;



  127. %主分析程序

  128. clear all
  129. h=pi/256;
  130. w=2100;
  131. tf=300000*2*pi/w;
  132. tspan=0:h:tf;
  133. y0=[0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5];
  134. options=odeset('RelTol',10^-6,'AbsTol',10^-6);  
  135. [t,y]=ode45(@fun,tspan,y0);

  136. figure(1)
  137. subplot(2,2,1);
  138. plot(t,y(:,1),'r',t,y(:,2),'b')
  139. title('未松动端位移响应');xlabel('t');ylabel('x1-red,y1-blue');
  140. subplot(2,2,2);
  141. plot(t,y(:,3),'r',t,y(:,4),'b')
  142. title('圆盘位移响应');xlabel('t');ylabel('x2-red,y2-blue');
  143. subplot(2,2,3);
  144. plot(t,y(:,5),'r',t,y(:,6),'b')
  145. title('松动端轴心位移响应');xlabel('t');ylabel('x3-red,y4-blue');
  146. subplot(2,2,4);
  147. plot(t,y(:,7),'b')
  148. title('m3在竖直方向位移响应');xlabel('t');ylabel('y4');

  149. figure(2)
  150. subplot(2,2,1);
  151. plot(y(:,1),y(:,2))
  152. title('未松动端轴心轨迹');xlabel('x1');ylabel('y1');
  153. subplot(2,2,2);
  154. plot(y(:,5),y(:,6))
  155. title('松动端轴心轨迹');xlabel('x3');ylabel('y3');

  156. yy1=sqrt(y(:,1).^2+y(:,2).^2);%yy1=sqrt(x1^2+y1^2)
  157. yy2=sqrt(y(:,5).^2+y(:,6).^2);%yy2=sqrt(x3^2+y3^2)
  158. N1=256;fs=1024;stepf=fs/N1;k=pi*2/882.5056;%882.5056转子固有频率sqrt(k/m2)
  159. n=0:stepf*k:(fs/2-stepf)*k;
  160. YY1=abs(fft(yy1));
  161. YY2=abs(fft(yy2));
  162. figure(3)
  163. plot(n,abs(YY1(1:N1/2)));grid on;
  164. title('未松动情况下幅值谱');xlabel('频率比');ylabel('FFT幅值');
  165. figure(4)
  166. plot(n,abs(YY2(1:N1/2)));grid on;
  167. title('松动情况下幅值谱');xlabel('频率比');ylabel('FFT幅值');
  168.       
  169. fclose('all');


复制代码
主要参考文献:李振平,罗跃纲等.转子系统支承松动的非线性动力学及故障特征[J].东北大学学报(自然科学版),2002,23(11):1049-1051.

点评

赞成: 4.0
赞成: 4
不错的帖子!!  发表于 2014-3-27 21:32

评分

3

查看全部评分

回复
分享到:

使用道具 举报

发表于 2007-5-24 20:21 | 显示全部楼层
收藏了,谢谢楼主!
 楼主| 发表于 2007-5-25 12:19 | 显示全部楼层
发现错误:
其中
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 编辑 ]

点评

赞成: 4.0
赞成: 4
  发表于 2014-3-27 21:34

评分

1

查看全部评分

发表于 2007-5-25 13:30 | 显示全部楼层
已经帮你在原帖修正

[ 本帖最后由 ChaChing 于 2010-3-1 10:13 编辑 ]
发表于 2007-5-27 22:17 | 显示全部楼层

感触下先

好东西,学习的时候就是少了点程序,看理论还是有点远:lol
发表于 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.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> main_songdong at 11
[t,y]=ode45(@fun,tspan,y0);

评分

1

查看全部评分

发表于 2007-10-21 19:10 | 显示全部楼层
有没有相应的理论方法介绍?
发表于 2007-10-21 23:25 | 显示全部楼层
原帖由 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,否则恐怕要楼主现身才能修正

点评

从后面给出的模型来看,就是C 这是对激振力的无量钢化  发表于 2010-10-18 05:47
发表于 2008-5-12 10:34 | 显示全部楼层

回复 11楼 的帖子

程序里面的c就是dert=0.00011,也就是油膜力的平均厚度。个人感觉程序有点差异。
C=C/w;
K=K/w^2;
F=F/c/w^2;
应该是对方程进行无量纲化,既然无量纲化了系统的周期也就变成了2pi呀!望和楼主讨论:@)
还有,前期的运动一般都不稳定,要去掉瞬态解。
刚接触非线性不久,说的不对的地方请楼主包涵:loveliness:

点评

赞成: 4.0
赞成: 4
  发表于 2014-3-27 21:35

评分

2

查看全部评分

发表于 2010-9-27 21:08 | 显示全部楼层
我太需要这个了!谢谢楼主!
发表于 2010-9-28 19:09 | 显示全部楼层
本帖最后由 mjtjiang 于 2010-9-28 19:11 编辑
pwangeng 发表于 2007-5-24 15:20
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会出现复数,其中的虚部有什么意义吗?
发表于 2010-9-28 19:36 | 显示全部楼层
mjtjiang 发表于 2010-9-28 19:09
其中用abs是什么意思呢?若不用abs会出现复数,其中的虚部有什么意义吗?

abs应该是绝对值的意思
我好像见过这个模型,应该是有这个绝对值得
发表于 2010-9-29 09:12 | 显示全部楼层
回复 雪缘 的帖子

虽然还不是很懂,但先记下吧。

求这个非线性油膜力的时候,为什么前面要乘以-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
发表于 2010-9-29 09:31 | 显示全部楼层
mjtjiang 发表于 2010-9-29 09:12
回复 雪缘 的帖子

虽然还不是很懂,但先记下吧。

我查了一下,这是1991年法国人Capone提出来的短圆柱瓦轴承油膜力模型,建议你找这个模型对照一下就都清楚了

问题1:不是 -0.1 是 -1.0 ,这表示力的方向问题;
问题2:是无量纲化处理,两者的解理论上应该是完全等价的,但是实际上有可能由于计算过程中的截断误差而导致结果不一样,无量纲化目的有两个,一是避免引入过大的截断误差,二是为了定性分析;
问题3:这一部分是将二阶微分方程转换为一阶状态方程,得到ode45所能求解的微分方程标准格式。

点评

赞成: 4.0
赞成: 4
  发表于 2014-3-27 21:36

评分

1

查看全部评分

发表于 2010-9-29 10:59 | 显示全部楼层
本帖最后由 mjtjiang 于 2010-9-29 16:40 编辑

你说的那个原文 Capone G. Descrizion analitica del campo di forze fluidodinamiconei cuscinetti cilindrici cubrificati[J].L’Energia  Elettrica,1991,3(1):105-110. 我没找到,但我有这个表达式:
clip_image002.jpg
clip_image002.jpg
我不是学转子动力学的,其原理我都不懂,但老师要我做的课题的一部分是转子故障信号仿真,我就只是按照书上给的公式编程序,发现有好多问题
问题一:楼主的这个程序是在程序里将二阶微分方程转换为一阶状态方程,我是在纸上手动计算把他降阶,在程序里直接编入的就是无量纲性形式的一阶状态方程,我这样做可以吗?反正程序能运行,但输出总是与书上给出的不符,我算出的轴心轨迹是像正弦信号一样震荡的。我也不知道自己错哪儿了。
问题二:楼主的
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;
这两句我不太懂,还请校长指教。


您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表