|
|
楼主 |
发表于 2007-9-18 10:41
|
显示全部楼层
- function dq=Rotors_System_Sub_Func(t,q)
- global BN Nb1 Nb2 w1 w2 Ro1 Ri1 Ro2 Ri2
- %-------------------------------------------------------------------------
- r01=0.00002;
- m1=7.86;
- m2=11.93;
- W1=m1*9.8;
- W2=m2*9.8;
- Kb1=4.6346e10;
- Kb2=3.4953e9;
- %-------------------------------------------------------------------------
- F=1485;
- F1 = 1485/2;
- delta_i=1.2443*w2^2*1e-11;
- delta_e=2.2183*w1^2*1e-11;
- Dm=(117+133.3)/2;
- wc=(w2*(1-8/Dm)+(1+8/Dm)*w1)/2;
- delta_o=4.912*wc^2*1e-14;
- delta_F=8.10*F1^0.925/7.8^0.85*1e-8;
- r02=0.00005-(delta_i-delta_o-delta_e-delta_F);
- %-------------------------------------------------------------------------
- Fx11=0;
- Fy11=0;
- Fx22=0;
- Fy22=0;
- for i=1:Nb1
-
- sita(i)=2*pi/Nb1*(i-1)+BN/Nb1*w1*2*pi*t;
- Deformation(i,1)=1e-6*q(1,1)*cos(sita(i))+1e-6*q(2,1)*sin(sita(i))-r01;
- if Deformation(i)<=0
- Deformation(i)=0;
- end
-
- fx11=Kb1*(Deformation(i))^1.5*cos(sita(i));
- fy11=Kb1*(Deformation(i))^1.5*sin(sita(i));
-
- Fx11=Fx11+fx11;
- Fy11=Fy11+fy11;
- end
- Fx12=Fx11;
- Fy12=Fy11;
- %-------------------------------------------------------------------------
- Fx21=0;
- Fy21=0;
- for j=1:Nb2
- sita(j)=2*pi/Nb2*(j-1)+(Ro2/(Ro2+Ri2)*w1+Ri2/(Ro2+Ri2)*w2)*2*pi*t;
- Deformation1(j,1)=1e-6*q(3)*cos(sita(j))+1e-6*q(4)*sin(sita(j))-r02;
- if Deformation1(j)<=0
- Deformation1(j)=0;
- end
- fx21=Kb2*(Deformation1(j))^1.1*cos(sita(j));
- fy21=Kb2*(Deformation1(j))^1.1*sin(sita(j));
-
- Fx21=Fx21+fx21;
- Fy21=Fy21+fy21;
- end
- %-------------------------------------------------------------------------
- for k=1:Nb1
- sita(k)=2*pi/Nb1*(k-1)+BN/Nb1*w2*2*pi*t;
- Deformation2(k,1)=1e-6*q(3,1)*cos(sita(k))+1e-6*q(4,1)*sin(sita(k))-r01;
- if Deformation2(k)<=0
- Deformation2(k)=0;
- end
-
- fx22=Kb1*(Deformation2(k))^1.5*cos(sita(k));
- fy22=Kb1*(Deformation2(k))^1.5*sin(sita(k));
-
- Fx22=Fx22+fx22;
- Fy22=Fy22+fy22;
- end
- %-------------------------------------------------------------------------
- P=1470;
- Cx11=1500;Cy11=1500;
- Cx12=2500;Cy12=2500;
- Cx21=2000;Cy21=2000;
- Cx22=7000;Cy22=7000;
- %-------------------------------------------------------------------------
- dq(5:8,1)=[-1e6/m1*((Cx11+Cx12)*q(5,1)*1e-6+Fx11+Fx12-(W1+Fx21+Cx21*q(7,1)*1e-6));...
- -1e6/m1*((Cy11+Cy12)*q(6,1)*1e-6+Fy11+Fy12-(Fy21+Cy21*q(8,1)*1e-6));...
- -1e6/m2*((Cx21+Cx22)*q(7,1)*1e-6+Fx21+Fx22-(P+W2));...
- -1e6/m2*((Cy21+Cy22)*q(8,1)*1e-6+Fy21+Fy22)];
- dq(1:4,1)=q(5:8,1);
主程序：
- function Rotors_System_Func
- clear
- clc
- global BN Nb1 Nb2 w1 w2 Ro1 Ri1 Ro2 Ri2
- %-------------------------------------------------------------------------
- n_one_T=100;
- n_T=400;
-
- Ro1=4.885e-2;
- Ri1=6.865e-2;
- Nb1=10;
-
- Ro2=66.515e-3;
- Ri2=58.5e-3;
- Nb2=34;
- w2_min=100;
- w2_max=10400;
- w2_step=100;
- w1_rpm=11000;
- q_initial(1:8,1)=1e-11;
- BN=Ri1/(Ri1+Ro1)*Nb1;
- %-------------------------------------------------------------------------
- tic
- w2_rpm=10000;
- w2=w2_rpm/60;
- w1=w1_rpm/60;
- r=Ro2/(Ro2+Ri2);
-
- w_vc_1=w1*BN
-
- w_cage_3=w1*r+w2*(1-r);
- w_vc_3=w_cage_3*Nb2
-
-
- w_cage=w2*BN;
- w_vc=w_cage
-
- T_vc=2*pi/w_vc;
- dt=T_vc/n_one_T;
- time=n_T*T_vc;
- n=round(time/dt);
- t_span(1:n)=linspace(0,time,n);
- [t,q]=ode45('Rotors_System_Sub_Func',t_span,q_initial);
- %-------------------------以下为绘制系统频谱图程序--------------------------
- xc =q(end-8192:4:end,1);
- yc =q(end-8192:4:end,2);
- xcc=q(end-8192:4:end,3);
- ycc=q(end-8192:4:end,4);
- X1 = fft(xc, 2048);
- X2 = fft(xcc,2048);
- X1(1) = [];
- X2(1) = [];
- Pxx = X1.* conj(X1) / 2048;
- Pyy = X2.* conj(X2) / 2048;
- f = 0.25*1/dt*(0:1024)/2048;
- subplot(221);plot(t_span(end-500:end),q(end-500:end,1))
- xlabel('t (s)')
- ylabel('x1 (m)')
- subplot(222);plot(f(1:1024),Pxx(1:1024))
- xlabel('频率(Hz)')
- ylabel('功率')
- subplot(223);plot(t_span(end-500:end),q(end-500:end,3))
- xlabel('t (s)')
- ylabel('x2 (m)')
- subplot(224);plot(f(1:1024),Pyy(1:1024))
- xlabel('频率(Hz)')
- ylabel('功率')
- %-------------------------------------------------------------------------
- toc
|
|
显身卡
发表于 2007-9-18 14:47
发表于 2007-9-18 16:50
