弯扭耦合振动分岔图

lazyeboy

某一偏心下，转子耦合振动的分岔图如附件中所示，大家帮我看看，为什么在最初的转速下振动幅值为0，这样的分岔图对不对，应该如何分析？w为转动角速度（rad/s）。所用的程序如下：

w=1:100:10100;
for  n=1:length(w);
      T=2*pi;
      ts=[0:T/100:100*T];                        
      zo=[z1;z2;z3;z4;z5;z6;z7;z8;z9;z10;
      z11;z12;z13;z14;z15;z16;z17;z18;z19;z20;
      z21;z22;z23;z24;z25;z26;z27;z28;z29;z30;];                 
      [tt,z]=ode45('wn3f1',ts,zo,[],w(n));
      
      figure(1)
      plot(w(n),z(5000:100:10000,1),'*');
      xlabel('\fontsize{18}\omega');
      ylabel('\fontsize{18}x');grid
      hold on

      figure(2)
      plot(w(n)*30/pi,z(5000:100:10000,3),'*');
      xlabel('\fontsize{18}\omega');
      ylabel('\fontsize{18}y');grid
      hold on

end

[ 本帖最后由 lazyeboy 于 2007-6-21 20:27 编辑 ]

lazyeboy

扭转振动分岔图

咕噜噜

m函数文件呢？

无水1324

化分岔图的时候不要用“*”，这个会产生很多毛刺。
应该是概周期进入混沌的过程
改为：plot(w(n)*30/pi,z(5000:100:10000,3),'k.');

lazyeboy

多谢指点！m函数文件是5个轮盘弯扭耦合振动方程，由于文件较大，就不上传了。还想问一句的是，我的这个程序应当没有问题吧

无水1324

刚才没有看清楚，发现你的周期应该跟你这里面的w有关系吧？
有关系就得改为：T=2*pi/w(n);

gghhjj

为什么在最初的转速下振动幅值为0，这样的分岔图对不对

理解有问题，这里的振动位移和你取点有关系
分叉图上的点仅表示庞加莱截面上的映射，其位移大小表示系统在庞加来截面上的位移，而非振动幅值

gghhjj

另外方程式有量纲的还是无量纲的

最好把wn3f1的代码贴出来或者给出方程

lazyeboy

方程是经过无量纲处理的，T=2*pi，就不应除以w(n)了吧，否则不同转速下每个周期的点数是不同的，就不能每隔100个点取值了。
wn3f1文件如下：
function zdot=wn3f1(tt,z,flag,w);
zdot=zeros(30,1);

省略参数的取值。。。。。

zdot(1)=z(2);                                                               % 轮盘1振动方程
zdot(2)=E1*((1+z(6))^2*cos(tt+z(5)+a1)+zdot(6)*sin(tt+z(5)+a1))...
        -K1*(z(1)-z(7))...
        -C1*(z(2)-z(8));
zdot(3)=z(4);
zdot(4)=E1*((1+z(6))^2*sin(tt+z(5)+a1)-zdot(6)*cos(tt+z(5)+a1))...
        -K1*(z(3)-z(9))...
        -C1*(z(4)-z(10))-D;
zdot(5)=z(6);                                                   
zdot(6)=(E1*(zdot(2)*sin(tt+z(5)+a1)-(zdot(4)+D)*cos(tt+z(5)+a1))...
         -Ct1*(z(6)-z(12))...
         -Kt1*(z(5)-z(11)))/(h1^2/2+E1^2);

         
zdot(7)=z(8);                                                               % 轮盘2振动方程
zdot(8)=E2*((1+z(12))^2*cos(tt+z(11)+a2)+zdot(12)*sin(tt+z(11)+a2))...
        -K2*(2*z(7)-z(1)-z(13))...
        -C2*(2*z(8)-z(2)-z(14))-Kz1*z(7);
zdot(9)=z(10);
zdot(10)=E2*((1+z(12))^2*sin(tt+z(11)+a2)-zdot(12)*cos(tt+z(11)+a2))...
        -K2*(2*z(9)-z(3)-z(15))...
        -C2*(2*z(10)-z(4)-z(16))-D-Kz1*z(9);
zdot(11)=z(12);                                                         
zdot(12)=(E2*((zdot(8))*sin(tt+z(11)+a2)-(zdot(10)+D)*cos(tt+z(11)+a2))...
         -Ct2*(2*z(12)-z(6)-z(18))...
         -Kt2*(2*z(11)-z(5)-z(17)))/(h2^2/2+E2^2);

         
zdot(13)=z(14);                                                             % 轮盘3振动方程
zdot(14)=E3*((1+z(18))^2*cos(tt+z(17)+a3)+zdot(18)*sin(tt+z(17)+a3))...
        -K3*(2*z(13)-z(7)-z(19))...
        -C3*(2*z(14)-z(8)-z(20));
zdot(15)=z(16);
zdot(16)=E3*((1+z(18))^2*sin(tt+z(17)+a3)-zdot(18)*cos(tt+z(17)+a3))...
        -K3*(2*z(15)-z(9)-z(21))...
        -C3*(2*z(16)-z(10)-z(22))-D;
zdot(17)=z(18);                                                            
zdot(18)=(E3*(zdot(14)*sin(tt+z(17)+a3)-(zdot(16)+D)*cos(tt+z(17)+a3))...
         -Ct3*(2*z(18)-z(12)-z(24))...
         -Kt3*(2*z(17)-z(11)-z(23)))/(h3^2/2+E3^2);  
     

zdot(19)=z(20);                                                             % 轮盘4振动方程
zdot(20)=E4*((1+z(24))^2*cos(tt+z(23)+a4)+zdot(24)*sin(tt+z(23)+a4))...
        -K4*(2*z(19)-z(13)-z(25))...
        -C4*(2*z(20)-z(14)-z(26))-Kz2*z(19);
zdot(21)=z(22);
zdot(22)=E4*((1+z(24))^2*sin(tt+z(23)+a4)-zdot(24)*cos(tt+z(23)+a4))...
        -K4*(2*z(21)-z(15)-z(27))...
        -C4*(2*z(22)-z(16)-z(28))-D-Kz2*z(21);
zdot(23)=z(24);                                                         
zdot(24)=(E4*((zdot(20))*sin(tt+z(23)+a4)-(zdot(22)+D)*cos(tt+z(23)+a4))...
         -Ct4*(2*z(24)-z(18)-z(30))...
         -Kt4*(2*z(23)-z(17)-z(29)))/(h4^2/2+E4^2);
                  
         
zdot(25)=z(26);                                                             % 轮盘5振动方程
zdot(26)=E5*((1+z(30))^2*cos(tt+z(29)+a5)+zdot(30)*sin(tt+z(29)+a5))...
         -K5*(z(25)-z(19))...
         -C5*(z(26)-z(20));
zdot(27)=z(28);
zdot(28)=E5*((1+z(30))^2*sin(tt+z(29)+a5)-zdot(30)*cos(tt+z(29)+a5))...
         -K5*(z(27)-z(21))...
         -C5*(z(28)-z(22))-D;      
zdot(29)=z(30);                                                      
zdot(30)=(E5*(zdot(26)*sin(tt+z(29)+a5)-(zdot(28)+D)*cos(tt+z(29)+a5))...
         -Ct5*(z(30)-z(24))...
         -Kt5*(z(29)-z(23)))/(h5^2/2+E5^2);

gghhjj

你这个程序能够运行？感觉很奇怪

你的函数要求输入的参数为flag,w两个

而

1. [tt,z]=ode45('wn3f1',ts,zo,[],w(n));

复制代码

中仅输入了w(n)一个量

gghhjj

另外，像这样长的程序，如果希望别人帮你看就应该给出相关参数

别人是没有这么多时间慢慢读程序的

通常都是运行以下程序，然后根据经验调整一些参数来发现问题

你这里也省略参数的取值那里也省略参数的取值，根本没办法运行

lazyeboy

首先谢谢指点，完整的程序如下：

function zdot=wn3f1(tt,z,flag,w);
zdot=zeros(30,1);


m=22.89240653;                         % 转子的总质量

m1=8.11334199;                         % 各轮盘或轴段简化成的轮盘的质量
m2=1.71254673;
m3=4.23564871;
m4=1.72335468;
m5=7.34341583;

g=9.8;

e1=0.001;                         % 轮盘的偏心距
e2=0.001;
e3=0.001;
e4=0.001;
e5=0.001;

R1=0.128;                              % 轮盘的半径
R2=0.075;
R3=0.08;
R4=0.075;
R5=0.139;

d=0.02;                                % 单位位移

E1=e1/d;                               % 轮盘的相对偏心率
E2=e2/d;
E3=e3/d;
E4=e4/d;
E5=e5/d;

h1=R1/d;
h2=R2/d;
h3=R3/d;
h4=R4/d;
h5=R5/d;


c1=800*pi;                              % 弯曲振动阻尼
c2=800*pi;
c3=800*pi;
c4=800*pi;
c5=800*pi;
C1=c1/(m1*w);                       % 无量纲化后弯曲阻尼
C2=c2/(m2*w);
C3=c3/(m3*w);
C4=c4/(m4*w);
C5=c5/(m5*w);


a1=30*180/pi;                           % 轮盘振动初相位
a2=30*180/pi;
a3=30*180/pi;
a4=30*180/pi;
a5=30*180/pi;


ct1=1.2*pi;                             % 扭转振动阻尼
ct2=1.2*pi;
ct3=1.2*pi;
ct4=1.2*pi;
ct5=1.2*pi;
Ct1=ct1/(m1*w*d^2);                 % 无量纲化后扭转阻尼
Ct2=ct2/(m2*w*d^2);
Ct3=ct3/(m3*w*d^2);
Ct4=ct4/(m4*w*d^2);
Ct5=ct5/(m5*w*d^2);

D=g/(w^2*d);

k1=3.7e7;                               % 弯曲振动刚度
k2=3.7e7;
k3=3.7e7;
k4=3.7e7;
k5=3.7e7;
K1=k1/(m1*w^2);                     % 无量纲化后弯曲振动刚度
K2=k2/(m2*w^2);
K3=k3/(m3*w^2);
K4=k4/(m4*w^2);
K5=k5/(m5*w^2);


kz1=1e9;                               % 支撑刚度
kz2=1e9;
Kz1=kz1/(m2*w^2);                   % 无量纲化后的支撑刚度
Kz2=kz2/(m4*w^2);


kt1=8.45e5;                             % 扭转振动刚度
kt2=8.45e5;
kt3=8.45e5;                           
kt4=8.45e5;
kt5=8.45e5;                           
Kt1=kt1/(m1*w^2*d^2);               % 无量纲化后扭转振动刚度
Kt2=kt2/(m2*w^2*d^2);
Kt3=kt3/(m3*w^2*d^2);
Kt4=kt4/(m4*w^2*d^2);
Kt5=kt5/(m5*w^2*d^2);


zdot(1)=z(2);                                                               % 轮盘1振动方程
zdot(2)=E1*((1+z(6))^2*cos(tt+z(5)+a1)+zdot(6)*sin(tt+z(5)+a1))...
        -K1*(z(1)-z(7))...
        -C1*(z(2)-z(8));
zdot(3)=z(4);
zdot(4)=E1*((1+z(6))^2*sin(tt+z(5)+a1)-zdot(6)*cos(tt+z(5)+a1))...
        -K1*(z(3)-z(9))...
        -C1*(z(4)-z(10))-D;
zdot(5)=z(6);                                                               % 只考虑连接轴段的扭转刚度和扭转阻尼(1+z(6))
zdot(6)=(E1*(zdot(2)*sin(tt+z(5)+a1)-(zdot(4)+D)*cos(tt+z(5)+a1))...
         -Ct1*(z(6)-z(12))...
         -Kt1*(z(5)-z(11)))/(h1^2/2+E1^2);

         
zdot(7)=z(8);                                                               % 轮盘2振动方程
zdot(8)=E2*((1+z(12))^2*cos(tt+z(11)+a2)+zdot(12)*sin(tt+z(11)+a2))...
        -K2*(2*z(7)-z(1)-z(13))...
        -C2*(2*z(8)-z(2)-z(14))-Kz1*z(7);
zdot(9)=z(10);
zdot(10)=E2*((1+z(12))^2*sin(tt+z(11)+a2)-zdot(12)*cos(tt+z(11)+a2))...
        -K2*(2*z(9)-z(3)-z(15))...
        -C2*(2*z(10)-z(4)-z(16))-D-Kz1*z(9);
zdot(11)=z(12);                                                             % 只考虑连接轴段的扭转刚度和扭转阻尼(1+z(12))
zdot(12)=(E2*((zdot(8))*sin(tt+z(11)+a2)-(zdot(10)+D)*cos(tt+z(11)+a2))...
         -Ct2*(2*z(12)-z(6)-z(18))...
         -Kt2*(2*z(11)-z(5)-z(17)))/(h2^2/2+E2^2);

         
zdot(13)=z(14);                                                             % 轮盘3振动方程
zdot(14)=E3*((1+z(18))^2*cos(tt+z(17)+a3)+zdot(18)*sin(tt+z(17)+a3))...
        -K3*(2*z(13)-z(7)-z(19))...
        -C3*(2*z(14)-z(8)-z(20));
zdot(15)=z(16);
zdot(16)=E3*((1+z(18))^2*sin(tt+z(17)+a3)-zdot(18)*cos(tt+z(17)+a3))...
        -K3*(2*z(15)-z(9)-z(21))...
        -C3*(2*z(16)-z(10)-z(22))-D;
zdot(17)=z(18);                                                             % 只考虑连接轴段的扭转刚度和扭转阻尼(1+z(18))
zdot(18)=(E3*(zdot(14)*sin(tt+z(17)+a3)-(zdot(16)+D)*cos(tt+z(17)+a3))...
         -Ct3*(2*z(18)-z(12)-z(24))...
         -Kt3*(2*z(17)-z(11)-z(23)))/(h3^2/2+E3^2);  
     

zdot(19)=z(20);                                                             % 轮盘4振动方程
zdot(20)=E4*((1+z(24))^2*cos(tt+z(23)+a4)+zdot(24)*sin(tt+z(23)+a4))...
        -K4*(2*z(19)-z(13)-z(25))...
        -C4*(2*z(20)-z(14)-z(26))-Kz2*z(19);
zdot(21)=z(22);
zdot(22)=E4*((1+z(24))^2*sin(tt+z(23)+a4)-zdot(24)*cos(tt+z(23)+a4))...
        -K4*(2*z(21)-z(15)-z(27))...
        -C4*(2*z(22)-z(16)-z(28))-D-Kz2*z(21);
zdot(23)=z(24);                                                             % 只考虑连接轴段的扭转刚度和扭转阻尼(1+z(24))
zdot(24)=(E4*((zdot(20))*sin(tt+z(23)+a4)-(zdot(22)+D)*cos(tt+z(23)+a4))...
         -Ct4*(2*z(24)-z(18)-z(30))...
         -Kt4*(2*z(23)-z(17)-z(29)))/(h4^2/2+E4^2);
                  
         
zdot(25)=z(26);                                                             % 轮盘5振动方程
zdot(26)=E5*((1+z(30))^2*cos(tt+z(29)+a5)+zdot(30)*sin(tt+z(29)+a5))...
         -K5*(z(25)-z(19))...
         -C5*(z(26)-z(20));
zdot(27)=z(28);
zdot(28)=E5*((1+z(30))^2*sin(tt+z(29)+a5)-zdot(30)*cos(tt+z(29)+a5))...
         -K5*(z(27)-z(21))...
         -C5*(z(28)-z(22))-D;      
zdot(29)=z(30);                                                             % 只考虑连接轴段的扭转刚度和扭转阻尼(1+z(30))
zdot(30)=(E5*(zdot(26)*sin(tt+z(29)+a5)-(zdot(28)+D)*cos(tt+z(29)+a5))...
         -Ct5*(z(30)-z(24))...
         -Kt5*(z(29)-z(23)))/(h5^2/2+E5^2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=1:100:10100;           %转子的转速是以rad/s为单位的转速


%位移、速度、转角、转角速度的初值
z1=1.4573958e-005;   z2=1.3450717e-005;   z3=-1.5054583e-5;     z4=3.8508189e-005;      z5=2.1710014e-005;   z6=-3.6692404e-5;    z7=-1.4573958e-005;   z8=1.3450717e-005;      
z9=1.5054583e-5;     z10=4.8508189e-005;  z11=2.1716014e-005;   z12=-4.6692404e-5;     
z13=3.1713014e-5;    z14=-2.6692404e-5;   z15=2.4573958e-005;   z16=3.3450717e-005;     
z17=2.113958e-005;   z18=2.3330717e-005;  z19=-3.2454583e-5;    z20=4.4508189e-005;     
z21=3.2310014e-005;  z22=-3.3192404e-5;   z23=-1.2873958e-005;  z24=2.6650717e-005;     
z25=2.1054583e-5;    z26=3.2508189e-005;  z27=1.3716014e-005;   z28=-3.21692404e-5;  z29=3.56513014e-5;   z30=-3.71192404e-5;                                                



for  n=1:length(w);
      T=2*pi;
      ts=[0:T/100:100*T];                          % 时间区间/w(n)
      zo=[z1;z2;z3;z4;z5;z6;z7;z8;z9;z10;
      z11;z12;z13;z14;z15;z16;z17;z18;z19;z20;
      z21;z22;z23;z24;z25;z26;z27;z28;z29;z30;];       % xo包含的三十个向量           
      [tt,z]=ode45('wn3f1',ts,zo,[],w(n));
      
      figure(1)
      plot(w(n),z(5000:100:10000,1),'*');
      xlabel('\fontsize{18}\omega');
      ylabel('\fontsize{18}x');grid
      hold on
      
      figure(2)
      plot(w(n),z(5000:100:10000,3),'*');
      xlabel('\fontsize{18}\omega');
      ylabel('\fontsize{18}y');grid
      hold on
end

无水1324

这个是正确的，flag是个标示符号！

[ 本帖最后由 无水1324 于 2007-6-23 09:29 编辑 ]

无水1324

原帖由 lazyeboy 于 2007-6-23 09:18 发表

                               
登录/注册后可看大图

首先谢谢指点，完整的程序如下：

function zdot=wn3f1(tt,z,flag,w);
zdot=zeros(30,1);


m=22.89240653;                         % 转子的总质量

m1=8.11334199;                         % 各轮盘 ...

你的程序应该是正确的。适当延长计算的时间，画图看看

lazyeboy

好的，谢谢！虽然现在计算的时间比较长，只要程序不错就好，不过我会试试的。