帮忙看看这个转子轴承碰摩程序

秋月

function dy=youmofangcheng(t,x)
%参数
global w kc
m1=4;%轴承处集中质量
m2=32.1;%圆盘处集中质量
R=0.025;%轴承半径
L=0.012;%轴承长度
c=0.00011;%油膜厚度
mu=0.018;%油膜黏度
% w=900;
f=0.1;%摩擦系数
c1=1050;%轴承阻尼
c2=2100;%圆盘阻尼
k=2.5e7;%弹性轴刚度
% kc=0;
b=0.00005;%偏心距
dd=0.000016;%转子与定子间隙
% P=16*9.8;
g=9.8;
d=b/c;
G1=g/(c*w^2);
G2=g/(c*w^2);
e=sqrt(x(5)^2+x(7)^2);
%判断是否碰摩
if e<dd
    rud=0;
else
    rud=1;
end

%碰摩力
Px=-c*(e-dd)*kc/e*(x(5)-f*x(7));
Py=-c*(e-dd)*kc/e*(f*x(5)+x(7));


%非线性油膜力
s=mu*w*R*L*(R/c)^2*(L/(2*R))^2;
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign(x(3)+2*x(2));
G=2/sqrt(1-x(1)^2-x(3)^2)*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/sqrt(1-x(1)^2-x(3)^2)));
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cos(a)+x(3)*sin(a))^2);
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G)/(1-x(1)^2-x(3)^2);
fx=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(1)*V-sin(a)*G-2*cos(a)*S);
fy=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(3)*V+cos(a)*G-2*sin(a)*S);

%方程组
dy=[x(2);
    -c1/(w*m1)*x(2)-k/(w^2*m1)*(x(1)-x(5))+s/(c*w^2*m1)*fx;
    x(4);
    -c1/(w*m1)*x(4)-k/(w^2*m1)*(x(3)-x(7))+s/(c*w^2*m1)*fy-G1;
    x(6);
    -c2/(w*m2)*x(6)-2*k/(w^2*m2)*(x(5)-x(1))+rud*Px/(c*w^2*m2)+d*cos(t);
    x(8);
    -c2/(w*m2)*x(8)-2*k/(w^2*m2)*(x(7)-x(3))+rud*Py/(c*w^2*m2)+d*sin(t)-G2];

帮忙看看这个方程有什么问题?
怎么计算出结果啊??
参数有问题?还是??
谢谢了!

无水1324

G1=g/(c*w^2);
G2=g/(c*w^2);
怎么是一样的？

秋月

谢谢!请高手帮忙看看!
因为重力化为无量纲量把质量削掉了没有关系!

咕噜噜

说实话没看太明白你的程序，因为没有找到主程序

sssssxxxxx921

她那上边就没有主程序  找到就不对了   呵呵
她就是在问用什么方法求解了   我想数值方法还是用ode45好一点

咕噜噜

^_^，她说怎么计算出结果？呵呵，我还奇怪呢，没有主程序怎么出结果
用ode45完全可以得

octopussheng

global w kc
w=2;
kc=3000;
x0=[0 0 0 0 0 0 0 0];
t0=0;
tf=20;
[t,y]=ode45(@youmofangcheng,[t0 tf],x0)

上面的是求解程序，但是你程序里面的两个全局变量w kc我不知道取值多少，所以就试算了一下，发现解不出来！

秋月

kc的取值在10的6次方,w 取值 900

资料上说kc取0,w取值900能出现混沌

秋月

我也是用ode45做的
可是结果什么都不是
参数值有问题还是方程有问题呢?
帮忙看看
谢谢了!

秋月

clc
clear all
global w kc

w=900;
% kc=2.5e7;
kc=0;
x0=[0;0.1;0;0.1;0;0.1;0;0.1];
options=odeset;
options.reltol=1e-6;
T=2*pi/w;
[t,y]=ode45('youmofangcheng',[0:T/100:1000*T],x0);
plot(t,y(:,1))
figure
plot(y(:,5),y(:,7))
怎么不对啊
初值和步长有问题吗?:'(
高手帮帮忙!谢谢你们了!

octopussheng

和初值关系不大的，我算了一下，没有出来图形，呵呵，你还是修改主程序中的参数再试试吧！

秋月

有没有参考数据啊???
不知道哪个参数有问题!!:'(
谢谢了!

秋月

Warning: Failure at t=1.295358e-001.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.440892e-016) at time t.
> In ode45 at 355
  In youmojisuan at 12
>>
这个怎么会事啊?

octopussheng

这个就是参数不对导致积分最小误差不能满足，积分无法继续！论坛好些帖子都有讨论这个的，具体解决还需要修改参数，或者换个方法，如采用ode15等试试

无水1324

Warning: Divide by zero.
(Type "warning off MATLAB:divideByZero" to suppress this warning.)
> In C:\MATLAB6p5\work\youmofangcheng.m at line 33
  In C:\MATLAB6p5\toolbox\matlab\funfun\private\odearguments.m at line 104
  In C:\MATLAB6p5\toolbox\matlab\funfun\ode45.m at line 155

这是我运行的时候产生的Warning