求助，碰磨转子动力学模型分析，用龙格库塔方法求解没有结果？？

林天

我最近在研究碰磨转子的动力学模型，龙哥库塔法求解方程，并画分岔图和poincare截面图，但是用四阶龙格库塔法求解没有结果，不知道为什么，求大侠们帮我看看，万分感谢！！
下面是我的程序：
clc
clear
omega=900;            %想求系统随omega（转速）的变化分岔图和截面图，这里想用一个定值试试
period=2*pi;
tspan = (1:period/100:500*period);
y0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001];
[t,u] = ode45(@equ_elastic_without_UMP,tspan,y0);                          
%%%%%%运算没有结果，不知道为什么，错在哪？？？？，下面是微分方程
%%%%%微分方程（文献中的方程）%%%%%%%%%%%%%%%%%%%%%%
function  uu = equ_elastic_without_UMP(t,u)   %有问题
  m1 = 4;               %转子质量 kg
  m2 = 32.1;               %轴颈质量 kg
  R = 50*10^-3;         %轴承半径 m
  L = 20*10^-3;        %轴承长 m         
  miu=0.11;             %平均油膜厚度
  sigam=0.2;             %转静间隙
  mu=0.018;              %润滑油粘度
  c1 = 1050;             %转子阻尼 N.s/m
  c2 = 2100;             %轴承阻尼 N.s/m
  f=0.1;                 %摩擦系数  
  kc = 6.2*10^6;          %静子刚度 N/m
  k=2.5*0^7;              %弹性轴刚度
  e = 0.05;        %转子质量偏心 m                    
                                          
%delta_0 = 4.5*10^-3;   %均匀气隙大小 m                     
%mu = 18*10^-3;         %绝对润滑油粘度 无单位
%omega = 13;            %大轴旋转角速度
%global cz            %轴承间隙
global omega
  %%% 非线性油膜力表达式 fx  fy  %%%%
  M=m2/2;
  sigma = mu * omega * R * L/M * (R/miu)^2 * (L/2/R)^2;        %Sommerfeld修正系数%%%%%%%%改动(*2/m2)
  
  A1 = u(7) + 2*u(6);
  A2 = u(5) - 2*u(8);
  
  alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1);
  
E = sqrt(u(5).^2 +u(7).^2);                     %轴颈轴心轨迹
E_minus = sqrt(1 -E.^2);
B1 = u(7)*cos(alpha) - u(5)*sin(alpha);        
B2 = u(5)*cos(alpha) + u(7)*sin(alpha);
  
  G = 2./E_minus.^2 * ( pi/2 + atan(B1./ E_minus) );     %此处错误，应该是E_minus.^2;
  V = (2 + B1 * G) ./  E_minus.^2;
  S = B2 ./ ( 1 - B2.^2 );
  
  F_oil_same = -sqrt( A2.^2 + A1.^2 ) ./ E_minus.^2;
  
  f_x = F_oil_same * ( 3* u(5) *V - sin(alpha) * G - 2 *cos(alpha)*S );
  f_y = F_oil_same * ( 3* u(7) *V + cos(alpha) * G - 2 *sin(alpha)*S );     %油膜力的无量纲表达式
  
  fx = sigma * f_x;
  fy = sigma * f_y;              %非线性油膜力
  %%%摩擦力计算%%%%%%%%%%%%%%%%%%%%%%%%%自己加的，注意会出问题
  r=sqrt(u(5).^2+u(7).^2);
  PX=-kc* miu.*(1-sigam/r).*(u(5)-f.*u(7));
  PY=-kc* miu.*(1-sigam/r).*(f.*u(5)+u(7));
  g=9.8;
  G=g/(miu*omega^2);
  b=e/miu;
  tao=omega*t;
  uu = zeros(8,1);               %不明白什么意思
  uu(1) = u(2);
  uu(2) = -c1/(m1*omega) * u(2) - k/(m1*omega^2) * (u(1)-u(5)) + sigma*M*fx/(m1*omega^2*miu);
  uu(3) = u(4);
  uu(4) = -c1/(m1*omega) * u(4) - k/(m1*omega^2) * (u(3)-u(7)) + sigma*M*fy/(m1*omega^2*miu)-G;
  uu(5) = u(6);
  uu(6) = -c2/(m2*omega) * u(6) - 2*k/(m2*omega^2)* (u(5)-u(1)) + PX/(m2*omega^2*miu) +b*cos(tao);
  uu(7) = u(8);
  uu(8) = -c2/(m2*omega) * u(8) - 2*k /(m2*omega^2) * (u(7)-u(3)) + PY/(m2*omega^2*miu) +b*sin(tao)-G;
%%%%%%%%%%%%%%%%%%%%%%%equ_elastic_without_UMP.m文件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

kezairenjian

方程定义为：function  uu = equ_elastic_without_UMP(t,flag,u)  ，其中u是分岔参数。

林天

回复 2 # kezairenjian 的帖子

您好，你做过这方面的东西吗？，其他地方还有错误吗？

kezairenjian

碰摩转子没做过，不过matlab会一点。

xiaoshihanlan

伤痕累累

你这个轴心轨迹图做出来了吗？
你做的是转子碰摩和轴端油膜力吧

boom

楼主你好，请问您的程序最后做好了吗，我也是在做碰摩转子，遇到很多问题，能否参考一下您的，万分感谢了！

paopaotaiqiang

我想问下，你是参考哪篇论文中的方程编写的程序呀

nishoulong

学习中，我也在研究中

yuanrenlian2011

我还没找到方向转子-轴承要看哪些书呀？

funi

yuanrenlian2011 发表于 2012-9-12 15:51

                               
登录/注册后可看大图

我还没找到方向转子-轴承要看哪些书呀？

转子动力学

牛行天下

我现在也在做这个，可以交流一下

gghhjj

不是没有结果，而是有错误啊

andongny

我也刚刚接触碰磨这块，新手上路啊~