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

fanhuasijin

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

1. function dxdt=Muszynska(t,x)
   
2. global w;
   
3. Ke=1.0e+008;De=500;   %轴的刚度、转子外阻尼
   
4. 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;%转子的相对偏心距       %未标注的为具体的密封参数
   
5. Ra=p*v*c/u;    %轴向雷诺数
   
6. Rv=p*R*c*w/u;    %周向雷诺数
   
7. ac=Ra/Rv;
   
8. lam=n0*Ra^(m0)*(1+(1/(2*ac))^2)^((1+m0)/2);
   
9. o=lam*l/c;
   
10. T=1/v;
    
11. B1=1/(1+(2*ac)^2);
    
12. ap=(1+z+2*o)*p*v^2/2;    %密封前后压差
    
13. B=1+4*ac^2*B1*(1+m0);
    
14. E=(1+z)/(2*(1+z+2*o));
    
15. u0=2*o^2*E*(1-m0)/(1+z+2*o);
    
16. u1=2*o^2*(E/o+B*(1/6+E)/2)/(1+z+2*o);
    
17. u2=o*(1/6+E)/(1+z+2*o);
    
18. u3=pi*ap*R/lam;
    
19. mf=u2*u3*T^2;
    
20. M=m+mf;
    
21. K0=u0*u3;
    
22. D0=u1*u3*T;
    
23. tao=tao1*(1-r)^b;
    
24. K=K0*(1-r)^(-n);
    
25. D=D0*(1-r)^(-n);
    
26. K1=(Ke+K-tao^2*w^2*mf)/(M*w^2);
    
27. K2=D*tao/(M*w);
    
28. D1=(De+D)/(M*w);
    
29. D2=2*mf*tao/M;
    
30. G=9.8*m/(M*c*w^2);
    
31. 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文件里有