%Rotor critical speed calculation转子临界转速计算
%*****************************************************************************
%转子数据
E=2.06e11;%轴的弹性模量N/m2
l=[0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08];%轴总长为72cm,分成9段每段的长度相等都是0.08m
d=[0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025,0.025];%轴的直径m
m=[0.308,0.308,9.21,0.308,0.308,0.308,0.308,0.308,0];%各段的质量
Id=[0,0,0.03196,0,0,0,0,0,0];%直径转动惯量kg*m2
Io=[0,0,0.06392,0,0,0,0,0,0];%极转动惯量
J=[1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8,1.917e-8];%J为各段面积惯性矩，J=pi*d^4/64
% E=2.06e7 %轴的弹性模量N/cm2
% l=[8,8,8,8,8,8,8,8,8];%轴总长为72cm,分成9段每段的长度相等都是8cm
% d=[2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5];%轴的直径
% m=[0.308,0.308,9.21,0.308,0.308,0.308,0.308,0.308,0];%各段的质量
% Id=[0,0,319.6,0,0,0,0,0,0];%直径转动惯量
% Io=[0,0,639.2,0,0,0,0,0,0];%极转动惯量
% J=[1.917,1.917,1.917,1.917,1.917,1.917,1.917,1.917,1.917];%J为各段面积惯性矩，J=pi*d^4/64
%*****************************************************************************
%1-9段的传递矩阵L
clc
L=cell(1,9);
for i=1:9
L{1,i}=[1,l(i),(l(i))^2/(2*E*J(i)),(l(i))^3/(6*E*J(i));
0, 1, l(i)/(E*J(i)), (l(i))^2/(2*E*J(i));
0, 0, 1, l(i);
0, 0, 0, 1];
end
%1-9站的传递矩阵M
M=cell(1,9);
syms w;
for i=1:9
if i<=2||i>=4
M{1,i}=[1, 0, 0, 0;
0, 1, 0, 0;
0, 0, 1, 0;
m(i)*w^2, 0, 0, 1];
else
M{1,i}=[1, 0, 0, 0;
0, 1, 0, 0;
0, ((Io(i)/Id(i))-1)*Id(i)*w^2, 1, 0;
m(i)*w^2, 0, 0, 1];
end
end
%*****************************************************************************
%传递矩阵计算
%左右两端的边界条件,转子两边看成是是简支的
% BJ9=[0,Seta9,0,Q9];
% BJ0=[0,Seta0,0,Q0];
% BJ9=M{1,9}*L{1,9}*M{1,8}*L{1,8}*M{1,7}*L{1,7}*M{1,6}*L{1,6}*M{1,5}*L{1,5}
% *M{1,4}*L{1,4}*M{1,3}*L{1,3}*M{1,2}*L{1,2}*M{1,1}*L{1,1}*BJ0
H4x4=M{1,9}*L{1,9}*M{1,8}*L{1,8}*M{1,7}*L{1,7}*M{1,6}*L{1,6}*M{1,5}*L{1,5}*...
M{1,4}*L{1,4}*M{1,3}*L{1,3}*M{1,2}*L{1,2}*M{1,1}*L{1,1};
h12=H4x4(1,2);h14=H4x4(1,4);h32=H4x4(3,2);h34=H4x4(3,4);
%采用作图法，根据传递矩阵法的原理，画出Δ随Ω（轴的转速）
%变化的曲线，曲线与Δ=0时的水平轴交点频率就为临界转速。
Drta=h12*h34-h14*h32;
% w=200:1:1600;
% plot(w,subs(Drta,w))%作图观察第一阶临界转速的变化趋势和Drta=0轴的交点
% xlabel('转速/rad/s','FontSize',10.5);ylabel('Δ值','FontSize',10.5);
% grid on
%*****************************************************************************
% a=subs(Drta,w)
% Drta=vpa(a)
%采用试凑法，选定一系列的w值，代入Drta式子中进行计算，得到一系列的Drta值，当Drta值满足一定的
%的精度时就可以认为此时的W值为临界转数。
for w=200:1:400
a=subs(Drta,w);
if a>=-0.01&&a<=0.001%计算精度的选取很重要，精度不能太小，否则得不到准确结果
w1cor=w
break
end
end
%注意利用传递矩阵发计算转子的临界速度对MATLAB的计算精度有很高的要求，所以在计算中一定要
%采用双精度数，且各个计算物理量都要采用国际单位,才能得到正确的临界转数
%*****************************************************************************
评论2