跪求用Riccati传递矩阵法计算磁悬浮转子临界转速的Matlab的程序,越详 ...

发布网友发布时间：2024-06-01 16:54

共1个回答

热心网友时间：2024-07-20 16:43

%求解转子系统前三个临界转速和主振型的传递矩阵法
%等截面轴参数
l1=0.12; d=0.04; A=pi*d*d/4;
%轮盘参数
D=0.5; h=0.025; %盘轴材料参数（忽略轴的质量）
a=1; u=0.3; rou=7800; E=2.0e11; G=E/(2*(1+u)); I=pi*(d^4)/64; K1=2.0e7;
v1=6*E*I/(a*G*A*l1*l1);
mi=rou*pi*D^2/4*h;%轮盘的集质量
Jp=mi*D^2/8; Jd=Jp/2; Ji=Jp-Jd; %参数的数组形式
L=[l1 l1 l1 l1 l1 l1 l1 l1 l1 l1 l1 l1 l1 0 0]; M=[0 mi mi mi mi mi mi 0 0 0 0 0 mi mi 0]; K=[K1 0 0 0 0 0 0 K1 0 0 0 K1 0 0 0];
v=[v1 v1 v1 v1 v1 v1 v1 v1 v1 v1 v1 v1 v1 0 0]; J=[0 Ji Ji Ji Ji Ji Ji 0 0 0 0 0 Ji Ji 0]; k=0;
Tit=['第一阶频率的振型和弯矩图';'第二阶频率的振型和弯矩图';'第三阶频率的振型和弯矩图'];
for w=0:0.01:4000; for i=1:15;
T(:,:,i)=[1+(L(i)^3)*(1-v(i))*(M(i)*w^2-K(i))/(6*E*I) L(i)+L(i)^2*J(i)*w^2/(2*E*I)
L(i)^2/(2*E*I) L(i)^3*(1-v(i))/(6*E*I); (L(i)^2)*(M(i)*w^2-K(i))/(2*E*I) 1+L(i)*J(i)*w^2/(E*I) L(i)/(E*I) L(i)^2/(2*E*I);
L(i)*(M(i)*w^2-K(i)) J(i)*w^2 1 L(i); M(i)*w^2-K(i) 0 0 1]; end H=T(:,:,1);
for i2=2:15;
H=T(:,:,i2)*H; end

F=H(3,1)*H(4,2)-H(3,2)*H(4,1); if F*(-1)^k < 0 %求解临界转速 k=k+1; wi(k)=w; w=wi(k)
ni(k)=wi(k)*30/pi; end end
for i1=1:3; w=wi(i1);
for j=1:14;
T(:,:,j)=[1+(L(j)^3)*(1-v(j))*(M(j)*w^2-K(j))/(6*E*I)
L(j)+L(j)^2*J(j)*w^2/(2*E*I)
L(j)^2/(2*E*I) L(j)^3*(1-v(j))/(6*E*I);
(L(j)^2)*(M(j)*w^2-K(j))/(2*E*I) 1+L(j)*J(j)*w^2/(E*I) L(j)/(E*I) L(j)^2/(2*E*I); L(j)*(M(j)*w^2-K(j)) J(j)*w^2 1 L(j); M(j)*w^2-K(j) 0 0 1]; end
H=T(:,:,1); for j=2:15;
H=T(:,:,j)*H; end
b=-H(4,1)/H(4,2); X(:,1)=([1 b 0 0]'); for n=2:16;
X(:,n)=T(:,:,n-1)*X(:,n-1); %相邻两质点右边的传递关系 end for j1=1:15; y(j1)=X(1,j1); z(j1)=X(3,j1); x(j1)=(j1-1)*l1; end
y(16)=X(1,16);
x(16)=1.56; z(16)=X(3,16);
y=y/max(abs(y));%归一化 z=z/max(abs(z)); subplot(3,1,i1)
plot(x,y,'b-',x,z,'r:') title(Tit(i1,:))
xlabel('轴长'),ylabel('不平衡值') axis([0,1.56,-1.2,1.2]) grid on z;
end

热心网友时间：2024-07-31 06:52