关于转子支承松动的仿真计算，为什么会有不同的结果（附程序）

肥振

main.txt (1.43 KB, 下载次数: 27)

2011-8-15 09:38 上传

点击文件名下载附件
原主程序

fun.txt (2.2 KB, 下载次数: 26)

2011-8-15 09:38 上传

点击文件名下载附件
原子程序

oilx.txt (835 Bytes, 下载次数: 24)

2011-8-15 09:38 上传

点击文件名下载附件
不变子程序

oily.txt (836 Bytes, 下载次数: 23)

2011-8-15 09:38 上传

点击文件名下载附件
不变子程序

main1.txt (1.46 KB, 下载次数: 31)

2011-8-15 09:38 上传

点击文件名下载附件
现主程序

fun11.txt (1.65 KB, 下载次数: 29)

2011-8-15 09:38 上传

点击文件名下载附件
现主程序

大家好，这帖子主要是：完整的支承松动数值仿真程序
http://forum.vibunion.com/forum-viewthread-tid-44094-fromuid-171095.html 一文的再讨论，小弟近来对原帖的程序做了小改动，但是出来的结果竟然会完全的不同，不知各位对此有什么看法？

附件是包含的matlab程序文件，分别是(main+fun+oilx+oily)和(main1+fun11+oilx+oily), 其中主要改动是在子函数fun11，各位可以试试运行。
在此对原帖的作者和原程序的编写者致以敬意，感谢为我们提供了思考提高的平台。
谢谢大家！期盼大家的回复！

肥振

主要改动是下面：
原子函数：
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;%初始间隙
c=delta2;
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; %%c=C?

for i=1:1:N/2 %这一部分是将二阶微分方程转换为一阶状态方程，得到ode45所能求解的微分方程标准格式
    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
t

改动后fun11子函数：
function d=fun11(t,y)

d=zeros(14,1);
%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;%初始间隙
c=delta2;
c1=1050;%转子圆盘处阻尼系数
c2=2100;%转子在轴承处阻尼系数
k1=7.5e7;
k2=2.5e9;
cb1=350;
cb2=500;

ox1=y(1);%未松动端竖直方向位移x1
odx1=y(2);
oy1=y(3);%未松动端竖直方向位移y1
ody1=y(4);

ox2=y(5);%圆盘位移x2
odx2=y(6);
oy2=y(7);%圆盘位移y2
ody2=y(8);

ox3=y(9);%松动端轴心位移x3
odx3=y(10);
oy3=y(11);%松动端轴心位移y3
ody3=y(12);

oy4=y(13);%质量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

c1=c1/w;
c2=c2/w;
cb=cb/w;
k=k/w^2;
kb=kb/w^2;
m1=m1/c/w^2;
m2=m2/c/w^2;
m3=m3/c/w^2;

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

fx=fx/c/w^2;
fy=fy/c/w^2;
fx1=fx1/c/w^2;
fy1=fy1/c/w^2;


d(1)=odx1;     %主要是这里改动
d(2)=-(c1/m1)*odx1-(k/m1)*(ox1-ox2)+(fx/m1);
d(3)=ody1;
d(4)=-(c1/m1)*ody1-(k/m1)*(oy1-oy2)+(fy/m1)-g;
d(5)=odx2;
d(6)=-(c2/m2)*odx2-(k/m2)*(ox2-ox1)-(k/m2)*(ox2-ox3)+e*w^2*cos(w*t);
d(7)=ody2;
d(8)=-(c2/m2)*ody2-(k/m2)*(oy2-oy1)-(k/m2)*(oy2-oy3)+e*w^2*sin(w*t)-g;
d(9)=odx3;
d(10)=-(c1/m1)*odx3-(k/m1)*(ox3-ox2)+(fx1/m1);
d(11)=ody3;
d(12)=-(c1/m1)*ody3-(k/m1)*(oy3-oy2)+(fy1/m1)-g;
d(13)=ody4;
d(14)=-(cb/m3)*ody4-(kb/m3)*oy4-(fy1/m3)-g;

t % 显示时间
end

肥振

不能让这帖子沉啊。。。各位帮帮忙

肥振

一个暑假，又开学了。。。这个帖子还是没人看。。。

肥振

一个寒假过去了，还是没人看。。。

feizhongxiu

谢谢分享

wangfengt

谢谢楼主

伤痕累累

这个我做过了，李振平发表在东北大学学报上的，程序跟原文的出入很大。没搞明白。难道是参数原因？一模一样的参数啊

xufuze12

请问楼主上面那程序的问题怎么样了啊？？解决了没啊？

Rotor2014

楼主解决问题了吗

1713573225

下载下来看看怎么回事