关于支承松动数值分析的Matab程序的一些问题

pwangeng

运动微分方程

function d=fun(t,y,w)

N=length(y);
w=700;



m1=4;         m2=32.1;   

e=0.00005;
c=0.00011;
k=2.5e7;
delta2=0.4;
c1=1050;
c2=2100;
k1=7.5e7;
k2=2.5e9;
cb1=350;
cb2=500;

ox1=y(1);
oy1=y(2);
odx1=y(5);
ody1=y(6);

ox2=y(3);
oy2=y(4);
odx2=y(7);
ody2=y(8);

if oy1<0
    cb=cb2;
    kb=k2;
elseif  (oy1>=0)&(oy1<=delta2)
        cb=0;
        kb=0;
    else
         cb=cb1;
         kb=k1;
end

M=[  m1    0     0    0;
      0     m1    0    0;
      0     0     m2   0;
      0     0     0    m2; ];

C=[  c1    0     -c1   0;
      0     c1+cb   0    -c1;
      -c2    0     c2   0;
      0    -c2    0    c2; ];
  
K=[   k       0       -k      0;
      0       k+kb      0      -k;
      -k    0       k     0;
      0     -k      0     k; ];
  


fx=oilx( ox1, oy1, odx1, ody1, w);
fy=oily( ox1, oy1, odx1, ody1, w);



F=[   fx;
      fy-m1*9.81;
      m2*e*w^2*cos(t);
      m2*e*w^2*sin(t)-m2*9.81; ];
   
C=C/w;
K=K/w^2;
F=F/c/w^2;

for i=1:1:N/2
    y1(i,1)=y(i);
    y2(i,1)=y(i+N/2);
end


yy2=inv(M)*(F-C*y2-K*y1);

for i=1:1:N/2
    d(i)=y2(i,1);
    d(i+N/2)=yy2(i,1);
end

x方向上油膜力函数

unction oilforce=oilx(x,y,dx,dy,wi)

R=0.025; L=0.012; miu=0.018;  dfai=1; dert=0.00011;

e=sqrt(x*x+y*y);
delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
ppp1=(y+2.0*dx)/(x-2.0*dy);
sign1=sign(ppp1);
   
   
   
ppp2=y+2.0*dx;
sign2=sign(ppp2);



alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));

f1=sqrt(abs(abs((x-2.0*dy)*(x-2.0*dy))+abs((y+2.0*dx)*(y+2.0*dx))))/(1.0-abs(x*x)-abs(y*y));

fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);

oilforce=fx*delta;

y方向上油膜力函数



function oilforce=oily(x,y,dx,dy,wi)


  R=0.025; L=0.012; miu=0.018;  dfai=1; dert=0.00011;

e=sqrt(x*x+y*y);
delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
ppp1=(y+2.0*dx)/(x-2.0*dy);
sign1=sign(ppp1);
   
   
   
ppp2=y+2.0*dx;
sign2=sign(ppp2);



alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));

f1=sqrt(abs(abs((x-2.0*dy)*(x-2.0*dy))+abs((y+2.0*dx)*(y+2.0*dx))))/(1.0-abs(x*x)-abs(y*y));

fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);

oilforce=fy*delta;

主程序

%clc


t0=0;
tf=1000
y0=[0.05  0.05  0.05  0.05  0.05  0.05  0.05  0.05;]';  

[t,y]=ode45(@fun,[t0,tf],y0);
plot(t,y)

fclose('all');

%crack3

我试图用ode45去求解运动微分方程，但运行结果为：

tf =

        1000

??? Error using ==> funfun\private\odearguments
FUN must return a column vector.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...

Error in ==> rub_time at 8
[t,y]=ode45(@fun,[t0,tf],y0);

谁能够帮我修改下主程序 求解微分方程 我对matlab确实不是很了解 但毕业设计时间紧迫 所以求助下l

zhlong

你是做故障诊断中的支撑松动模拟吗，把你的背景知识介绍一下吧。
这么长的程序这么看很花时间啊，如果又不了解你的课题背景。
还有你确定是主程序问题吗？

[ 本帖最后由 eight 于 2007-5-19 21:28 编辑 ]

pwangeng

我做的是故障诊断中的支承松动数值仿真，课题为：转子系统支承松动故障特性研究
力学模型是简化的Jeffcott模型，运动学方程是根据选定的模型写的，油膜力基本是按公式写了函数应该没问题的。我的程序是从一个碰摩的程序改过来的原来的程序自己写了个RK法，我觉得用ode45有了运动方程用ode45比较简单，所以想直接用ode45来解微分方程，但就是老不对。

[ 本帖最后由 eight 于 2007-5-19 21:28 编辑 ]

zhlong

ode45是变步长的，得出来的结果会不会不符合故障诊断中的信号一般是等时间间隔采样的要求

pwangeng

现在的问题不在于结果是否有错，而是根本得不出结果，程序运行错误，但我找不到错在哪里。

无水1324

你定义的 函数fun有问题，只能定义为：
function d=fun(t,y,flag,w)

调用的时候改为：
[t,y]=ode45(@fun,[t0,tf],y0,[],w);%w为数值量

pwangeng

运动方程是
m1x1''+c1(x1'-x2')+k(x1-x2)=fx
m1y1''+c1(y1'-y2')+k(y1-y2)+cby1'+kby1=fy-m1g
m2x2''+c2(x2'-x1')+k(x2-x1)=m2ew^2coswt
m2y2''+c2(y2'-y1')+k(y2-y1)=m2ew^2sinwt-m2g

我令
ox1=y(1);
oy1=y(2);
odx1=y(5);
ody1=y(6);

ox2=y(3);
oy2=y(4);
odx2=y(7);
ody2=y(8);

其中 ox1 oy1 odx1 ody1分别代表 x1 y1 x1' y1' 同理ox2...

我的思路是先把二阶微分方程降为一阶
y1'=y2
y2'=inv(M)*(F-Cy2-Ky1)

然后在采用积分方法求解 我的问题就卡在写不出求解程序了

再问下无水为什么函数fun只能你说的那样定义 其中flag代表了什么 调用时的[]又为什么

zhlong

??? Error using ==> funfun\private\odearguments
FUN must return a column vector.

根据这个提示，在fun函数最后加个语句d=d';

试了一下，ok

你检查下结果看看

pwangeng

这样改后程序是能运行，但新的问题出现了
第一：加上d=d';对其含义不明
第二：我所要求的是方程组的解 也就是x1,y1,x1',y1',x2,y2,x2',y2'的数值解 而这样运行出来的结果是：
tf =

        1000

pwangeng

运行了半天才把图也给画出来了  我如果想把x1,y1,x1',y1',x2,y2,x2',y2'的图形单独一个一个画改怎么改呢

pwangeng

刚理解错误了把d=d';想成了d的微分，若是d的转置就知道是什么了，谢谢大家的帮忙啊

eight

原帖由 pwangeng 于 2007-5-19 22:59 发表

                               
登录/注册后可看大图

刚理解错误了把d=d';想成了d的微分，若是d的转置就知道是什么了，谢谢大家的帮忙啊

建议多阅读基础书

pwangeng

其实谁都知道基础的很重要，花一定的时间和精力进去有很多事情都会变得简单。但有时候问题往往就在于根本没那么多充足的时间让我们从最最基础的开始。就比如说我们做毕业设计要用到Matlab，但我们不可能先去把Matlab学精通了再开始做，那样不现实。我们所能做到的只是针对自己所要用到的那一块进行一定程度的学习，然而这样的学习往往会在某些细节上过不去，其实这也是我们的无奈。

eight

原帖由 pwangeng 于 2007-5-20 03:32 发表

                               
登录/注册后可看大图

其实谁都知道基础的很重要，花一定的时间和精力进去有很多事情都会变得简单。但有时候问题往往就在于根本没那么多充足的时间让我们从最最基础的开始。就比如说我们做毕业设计要用到Matlab，但我们不可能先去把M ...

迟点我会发文讨论这个问题，暂时不便争论

0810064

学习—精通—实用是一个慢长的过程，也是一个折磨人的过程