有没有哪位高人做过Muszynska的密封力模型

fanhuasijin

我现在做的是密封力作用下的转子非线性研究，转子是单圆盘Jeffcott转子模型，密封力采用的是
Muszynska模型。用matlab仿真结果总是不收敛，有没有人遇到过这种问题，最后怎么解决的。请大家多指教。下面就是该模型的M文件

1. 
   
2. function dxdt=Muszynska(t,x)
   
3. global w;
   
4. Ke=1.0e+008;De=500;   %轴的刚度、转子外阻尼
   
5. u=2.0e-005;   %流体粘度  l=0.05;    %密封长度     c=0.001;       %密封间隙    R=0.5;    %转子半径   tao1=0.49; b=0.15;m=3000;    %转子质量   m0=-0.25;n0=0.066;v=50;   % 轴向速度  p=200;    % 流体密度 z=0.1;   %水力损失 n=2.5;    r=0.2;%转子的相对偏心距       %未标注的为具体的密封参数
   
6. Ra=p*v*c/u;    %轴向雷诺数
   
7. Rv=p*R*c*w/u;    %周向雷诺数
   
8. ac=Ra/Rv;
   
9. lam=n0*Ra^(m0)*(1+(1/(2*ac))^2)^((1+m0)/2);
   
10. o=lam*l/c;
    
11. T=1/v;
    
12. B1=1/(1+(2*ac)^2);
    
13. ap=(1+z+2*o)*p*v^2/2;    %密封前后压差
    
14. B=1+4*ac^2*B1*(1+m0);
    
15. E=(1+z)/(2*(1+z+2*o));
    
16. u0=2*o^2*E*(1-m0)/(1+z+2*o);
    
17. u1=2*o^2*(E/o+B*(1/6+E)/2)/(1+z+2*o);
    
18. u2=o*(1/6+E)/(1+z+2*o);
    
19. u3=pi*ap*R/lam;
    
20. mf=u2*u3*T^2;
    
21. M=m+mf;
    
22. K0=u0*u3;
    
23. D0=u1*u3*T;
    
24. tao=tao1*(1-r)^b;
    
25. K=K0*(1-r)^(-n);
    
26. D=D0*(1-r)^(-n);
    
27. K1=(Ke+K-tao^2*w^2*mf)/(M*w^2);
    
28. K2=D*tao/(M*w);
    
29. D1=(De+D)/(M*w);
    
30. D2=2*mf*tao/M;
    
31. G=9.8*m/(M*c*w^2);
    
32. dxdt=[x(2);-K1*x(1)-K2*x(3)-D1*x(2)-D2*x(4)+0.2*cos(t);x(4);K2*x(1)-K1*x(3)+D2*x(2)-D1*x(4)-G+0.2*sin(t)]
    

复制代码

做过的高手进来帮帮忙吧，急呀
以上所有参数值均来自一篇已发表的博士论文
下面是该模型的无量纲模型，但是后面缺了转子质量偏心力，在M文件里有

liliangbiao

你想要实现什么呢？

fanhuasijin

不平衡转子在密封力作用下的非线性运动，画它的分岔图。但是画出来的图和参考文献里的相去甚远，搞不明白什么原因

liliangbiao

那么这个非线性方程适用于频闪法吗？对初值敏感吗？另外，解方程是需要数值解还是迭代解，比如Rung-kutta法和常微分方程中可以根据系统的特征来假设解的格式，然后再逐步迭代的。我不知道你选用的是那种方法？

fanhuasijin

这个是用四阶Rung-kutta法解微分方程，我用最大值法画图，得到的结果是不收敛的，而文献里的结果即有分岔也有混沌现象而且无量纲位移在一定的范围里，所以我也搞不清楚这个是方程的问题还是程序的问题。不知道动力学方程除无量纲化以为还需要做什么处理吗。

lhy

把你的分岔图贴上来看看吧，这个本身没有什么问题，初值也很重要吧，你的初值赋给的是多少？

fanhuasijin

多谢各位了，我又找了个参数齐全的算例，并不是方程的问题也不是程序的问题，就是参数的值有些出入。谢谢了！