弯扭耦合振动分岔图
某一偏心下，转子耦合振动的分岔图如附件中所示，大家帮我看看，为什么在最初的转速下振动幅值为0，这样的分岔图对不对，应该如何分析？w为转动角速度（rad/s）。所用的程序如下：w=1:100:10100;
forn=1:length(w);
T=2*pi;
ts=;
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;];
=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 编辑 ] 扭转振动分岔图 m函数文件呢？ 化分岔图的时候不要用“*”，这个会产生很多毛刺。
应该是概周期进入混沌的过程
改为：plot(w(n)*30/pi,z(5000:100:10000,3),'k.'); 多谢指点！m函数文件是5个轮盘弯扭耦合振动方程，由于文件较大，就不上传了。还想问一句的是，我的这个程序应当没有问题吧 刚才没有看清楚，发现你的周期应该跟你这里面的w有关系吧？
有关系就得改为：T=2*pi/w(n); 为什么在最初的转速下振动幅值为0，这样的分岔图对不对
理解有问题，这里的振动位移和你取点有关系
分叉图上的点仅表示庞加莱截面上的映射，其位移大小表示系统在庞加来截面上的位移，而非振动幅值 另外方程式有量纲的还是无量纲的
最好把wn3f1的代码贴出来或者给出方程 方程是经过无量纲处理的，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); 你这个程序能够运行？感觉很奇怪
你的函数要求输入的参数为flag,w两个
而=ode45('wn3f1',ts,zo,[],w(n)); 中仅输入了w(n)一个量 另外，像这样长的程序，如果希望别人帮你看就应该给出相关参数
别人是没有这么多时间慢慢读程序的
通常都是运行以下程序，然后根据经验调整一些参数来发现问题
你这里也省略参数的取值那里也省略参数的取值，根本没办法运行 首先谢谢指点，完整的程序如下：
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;
forn=1:length(w);
T=2*pi;
ts=; % 时间区间/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包含的三十个向量
=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
回复 #10 gghhjj 的帖子
这个是正确的，flag是个标示符号！[ 本帖最后由 无水1324 于 2007-6-23 09:29 编辑 ] 原帖由 lazyeboy 于 2007-6-23 09:18 发表 http://www.chinavib.com/forum/images/common/back.gif
首先谢谢指点，完整的程序如下：
function zdot=wn3f1(tt,z,flag,w);
zdot=zeros(30,1);
m=22.89240653; % 转子的总质量
m1=8.11334199; % 各轮盘 ...
你的程序应该是正确的。适当延长计算的时间，画图看看 好的，谢谢！虽然现在计算的时间比较长，只要程序不错就好，不过我会试试的。
页:
[1]
2