秋月 发表于 2007-9-28 10:46

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

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 发表于 2007-9-28 11:25

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

秋月 发表于 2007-9-28 14:07

回复 #2 无水1324 的帖子

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

咕噜噜 发表于 2007-9-28 15:30

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

sssssxxxxx921 发表于 2007-9-28 15:41

回复 #4 咕噜噜 的帖子

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

咕噜噜 发表于 2007-9-28 15:59

回复 #5 sssssxxxxx921 的帖子

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

octopussheng 发表于 2007-9-28 16:12

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

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

秋月 发表于 2007-9-28 17:16

回复 #7 octopussheng 的帖子

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

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

秋月 发表于 2007-9-28 17:19

回复 #5 sssssxxxxx921 的帖子

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

秋月 发表于 2007-9-28 17:55

这是主程序

clc
clear all
global w kc

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

octopussheng 发表于 2007-9-28 18:12

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

秋月 发表于 2007-9-28 18:49

回复 #11 octopussheng 的帖子

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

秋月 发表于 2007-9-28 18:52

回复 #11 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 发表于 2007-9-28 19:01

回复 #13 秋月 的帖子

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

无水1324 发表于 2007-9-28 19:08

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
页: [1] 2 3 4
查看完整版本: 帮忙看看这个转子轴承碰摩程序