声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3165|回复: 13

[分形与混沌] 求助,碰磨转子动力学模型分析,用龙格库塔方法求解没有结果??

  [复制链接]
发表于 2012-3-29 09:20 | 显示全部楼层 |阅读模式

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

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

x
我最近在研究碰磨转子的动力学模型,龙哥库塔法求解方程,并画分岔图和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文件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

微分方程

微分方程
回复
分享到:

使用道具 举报

发表于 2012-3-30 08:07 | 显示全部楼层
方程定义为:function  uu = equ_elastic_without_UMP(t,flag,u)  ,其中u是分岔参数。
 楼主| 发表于 2012-3-30 08:58 | 显示全部楼层
回复 2 # kezairenjian 的帖子

您好,你做过这方面的东西吗?,其他地方还有错误吗?
发表于 2012-3-30 15:16 | 显示全部楼层
碰摩转子没做过,不过matlab会一点。
发表于 2012-4-12 20:26 | 显示全部楼层
发表于 2012-4-19 21:33 | 显示全部楼层
你这个轴心轨迹图做出来了吗?
你做的是转子碰摩和轴端油膜力吧
发表于 2012-5-22 02:27 | 显示全部楼层
楼主你好,请问您的程序最后做好了吗,我也是在做碰摩转子,遇到很多问题,能否参考一下您的,万分感谢了!
发表于 2012-6-29 15:44 | 显示全部楼层
我想问下,你是参考哪篇论文中的方程编写的程序呀
发表于 2012-9-12 10:58 | 显示全部楼层
学习中,我也在研究中
发表于 2012-9-12 15:51 | 显示全部楼层
我还没找到方向转子-轴承要看哪些书呀?
发表于 2012-9-12 18:45 | 显示全部楼层
发表于 2012-11-23 22:57 | 显示全部楼层
我现在也在做这个,可以交流一下
发表于 2012-11-28 16:30 | 显示全部楼层
不是没有结果,而是有错误啊
发表于 2013-3-12 11:25 | 显示全部楼层
我也刚刚接触碰磨这块,新手上路啊~
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-23 20:10 , Processed in 0.075554 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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