clc
clear
close
%%%%%%%%%%%%%%%%%轴段的状态参数
m=[0.0571 10.7023 0.2499 0.1538 7.0869 0.0385 0.04699 7.202 3.692 0.04699];
l=[0.0762 0.1778 0.1524 0.0508 0.0508 0.0508 0.1524 0.0508];
J=[2.6367e-9 2.6367e-9 2.6367e-9 2.6367e-9 2.6367e-9 2.1935e-8 2.1935e-8 2.1935e-8];
E=2.068e11;
Ip=[0.0859 0.0678 0.0429 0.0271];Id=[0.04295 0.0339 0.02145 0.01355];
k11=2.6269e7;k12=1.7513e7;k21=1.7513e7;k22=2e7;
v(1)=0;v(2)=0;
w1=1047;w2=1571;
p=-4000:1:6000;
nuaa=1;
for p=-4000:1:6000
for i=1:2
g(i)=Ip(i)*w1*p-Id(i)*p*p;
end
for i=3:4
g(i)=Ip(i)*w2*p-Id(i)*p*p;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%内转子
for i=1:5
a11=l(i)^3/3/E/J(i);
a12=l(i)^2/2/E/J(i);
a22=l(i)/E/J(i);
Tz(:,:,i)=[1 l(i) a12 a12*l(i)-a11;
0 1 a22 a22*l(i)-a12;
0 0 1 l(i);
0 0 0 1 ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%外转子
for i=6:8
a11=l(i)^3/3/E/J(i);
a12=l(i)^2/2/E/J(i);
a22=l(i)/E/J(i);
Tz(:,:,i)=[1 l(i) a12 a12*l(i)-a11;
0 1 a22 a22*l(i)-a12;
0 0 1 l(i);
0 0 0 1 ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%质量点传递矩阵
Td(:,:,1)=[1 0 0 0;0 1 0 0;0 0 1 0; m(1)*p*p-k22 0 0 1];
Td(:,:,3)=[1 0 0 0;0 1 0 0;0 0 1 0; m(3)*p*p 0 0 1];
Td(:,:,4)=[1 0 0 0;0 1 0 0;0 0 1 0; m(4)*p*p 0 0 1];
Td(:,:,6)=[1 0 0 0;0 1 0 0;0 0 1 0; m(6)*p*p-k22 0 0 1];
Td(:,:,7)=[1 0 0 0;0 1 0 0;0 0 1 0; m(7)*p*p-k22 0 0 1];
Td(:,:,10)=[1 0 0 0;0 1 0 0;0 0 1 0;m(10)*p*p 0 0 1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%盘的传递矩阵
Tp(:,:,1)=[1 0 0 0;0 1 0 0;0 g(1) 1 0;m(2)*p*p 0 0 1];
Tp(:,:,2)=[1 0 0 0;0 1 0 0;0 g(2) 1 0;m(5)*p*p 0 0 1];
Tp(:,:,3)=[1 0 0 0;0 1 0 0;0 g(3) 1 0;m(8)*p*p 0 0 1];
Tp(:,:,4)=[1 0 0 0;0 1 0 0;0 g(4) 1 0;m(9)*p*p 0 0 1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T11=Td(:,:,4)*Tz(:,:,3)*Td(:,:,3)*Tz(:,:,2)*Tp(:,:,1)*Tz(:,:,1)*Td(:,:,1);
T12=Td(:,:,6)*Tz(:,:,5)*Tp(:,:,2)*Tz(:,:,4);
T21=Td(:,:,10)*Tz(:,:,8)*Tp(:,:,4)*Tz(:,:,7)*Tp(:,:,3)*Tz(:,:,6)*Td(:,:,7);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A11=T12(3,1)*T11(1,1)+T12(3,2)*T11(2,1)+T12(3,3)*T11(3,1)+T12(3,4)*T11(4,1);
A12=T12(3,1)*T11(1,2)+T12(3,2)*T11(2,2)+T12(3,3)*T11(3,2)+T12(3,4)*T11(4,2);
A21=T12(4,1)*T11(1,1)+T12(4,2)*T11(2,1)+T12(4,3)*T11(3,1)+T12(4,4)*T11(4,1);
A22=T12(4,1)*T11(1,2)+T12(4,2)*T11(2,2)+T12(4,3)*T11(3,2)+T12(4,4)*T11(4,2);
% A(:,:,1)=[0,0,0,0,0;
% 0,0,0,0,0;
% 0,0,0,0,0;
% 0,0,0,0,0;
% 0,0,0,0,0];
A=[A11 A12 0 0 T12(3,4);
A21 A22 0 0 T12(4,4);
0 0 T21(3,1) T21(3,2) 0;
0 0 T21(4,1) T21(4,2) -1;
k22*T11(1,1) k22*T11(1,2) -k22*T21(1,1) -k22*T21(1,2) 1];
h(nuaa)=det(A);
v(1)=v(2);
v(2)=h(nuaa);
if v(1)*v(2)<0
H=p;
nuaa=nuaa+1;
end
plot(p,h);
set(gca,'XGrid','on','YGrid','on');
axis([-5000,6000,-10e45,10e45]);
hold on
end |