求助，碰磨转子动力学模型分析，用龙格库塔方法求解没有结果？？
我最近在研究碰磨转子的动力学模型，龙哥库塔法求解方程，并画分岔图和poincare截面图，但是用四阶龙格库塔法求解没有结果，不知道为什么，求大侠们帮我看看，万分感谢！！下面是我的程序：
clc
clear
omega=900; %想求系统随omega（转速）的变化分岔图和截面图，这里想用一个定值试试
period=2*pi;
tspan = (1:period/100:500*period);
y0=;
= ode45(@equ_elastic_without_UMP,tspan,y0);
%%%%%%运算没有结果，不知道为什么，错在哪？？？？，下面是微分方程
%%%%%微分方程（文献中的方程）%%%%%%%%%%%%%%%%%%%%%%
functionuu = 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
%%% 非线性油膜力表达式 fxfy%%%%
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文件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
方程定义为：functionuu = equ_elastic_without_UMP(t,flag,u)，其中u是分岔参数。 回复 2 # kezairenjian 的帖子
您好，你做过这方面的东西吗？，其他地方还有错误吗？ 碰摩转子没做过，不过matlab会一点。 {:{20}:} 你这个轴心轨迹图做出来了吗？
你做的是转子碰摩和轴端油膜力吧 楼主你好，请问您的程序最后做好了吗，我也是在做碰摩转子，遇到很多问题，能否参考一下您的，万分感谢了！ 我想问下，你是参考哪篇论文中的方程编写的程序呀 学习中，我也在研究中 我还没找到方向{:{19}:}转子-轴承要看哪些书呀？ yuanrenlian2011 发表于 2012-9-12 15:51 static/image/common/back.gif
我还没找到方向转子-轴承要看哪些书呀？
转子动力学 我现在也在做这个，可以交流一下{:{39}:} 不是没有结果，而是有错误啊 我也刚刚接触碰磨这块，新手上路啊~
页:
[1]