转子动力学传递矩阵法matlab计算程序

clear all;
clc;
%控制参数，用国际单位制
f=20000;                             %频率控制，单位为rpm
E=2.01e+011;                         %轴段材料的弹性模量单位为N/m^2
P=7850;                              %材料密度Kg/m^3
%安装轴承部位将轴同样分成两段，盘分成两段2个轴承，2个盘
l=[0.018;0.01;0.008;0.008;0.1;0.012;0.019;0.019;0.432;0.019;0.019;0.012;0.03;0.006;0.006;0.01;0.068;0.009];%分成18段，各轴段长度
d=[0.024;0.025;0.03;0.03;0.033;0.037;0.0419;0.0419;0.033;0.0419;0.0419;0.037;0.033;0.030;0.030;0.025;0.0235;0.018];%各轴段直径
num=[18;10;8;8;100;12;19;19;432;19;19;12;30;6;6;10;68;9];%各轴段分段数

m_shaft_total=0;
s=0;
m_shaft_L1=0;
for i=1:18
    m_shaft（i）=pi*d（i）^2/4*l（i）*P;        %各轴段质量
    m_shaft_total=m_shaft_total+m_shaft（i）;%轴总质量
    if s>3 && s<15
    m_shaft_L1=m_shaft_L1+m_shaft（s）;      %轴承中间的轴段质量
    end
    J（i）=pi*（d（i）^4）/64;%轴段截面惯性矩
    s=s+1;
end

Dl=0.26;                                                %盘L径
Bl=0.04;                                               %盘L厚
m_disk_L=pi*Dl^2/4*Bl*P-pi*d（7）^2/4*Bl*P;                   %盘L质量
Ip_L=1/2*m_disk_L*（Dl^2/4+d（7）^2/4）;                          %盘L极转动惯量
Id_L=1/12*m_disk_L*（3*（Dl^2/4+d（7）^2/4）+Bl^2）;                 %盘L直径转动惯量

Dr=0.28;                                                %盘R径
Br=0.04;                                               %盘R厚
m_disk_R=pi*Dr^2/4*Br*P-pi*d（14）^2/4*Br*P;                   %盘R质量
Ip_R=1/2*m_disk_R*（Dr^2/4+d（14）^2/4）;                          %盘R极转动惯量
Id_R=1/12*m_disk_R*（3*（Dr^2/4+d（14）^2/4）+Br^2）;                 %盘R直径转动惯量

M_total=m_shaft_total+m_disk_L+m_disk_R;                  %实验器总重量

bearing_k1=3e+006;
bearing_k2=1e+007;                  %支承刚度，根据实际情况给定bearing_k1为滚棒轴承支承刚度，bearing_k2为3#轴承支承刚度

P_matrix=[1000;
          0100;
          0010;
          0001];                     %初始矩阵

B_matrix1=[1000;
          0100;
          0010;
         -bearing_k1001]; 
     
B_matrix2=[1000;
          0100;
          0010;
         -bearing_k2001];

for i=1:1:f
    Omega=i/60*2*pi;   %角速度
    %第一大段，细分成18小段
    n=num（1）;
    N=1;
    L=l（N）/n;
     M_matrix=[1      000 ;
               0      100 ;
               0      010 ;
           m_shaft（N）/n*Omega^2 001]; 
    L_matrix= [1         L        （L^2）/（2*E*J（N））（L^3）/（6*E*J（N））;
               0         1           L/（E*J（N））      （L^2）/（2*E*J（N））;
               0         0           1                  L           ;
               0         0           0                  1          ];
    H_matrix=（M_matrix*L_matrix）^n*P_matrix;

    %第二大段，细分成10小段
    n=num（2）;
    N=2;
     L=l（N）/n;
     M_matrix=[1      000 ;
               0      100 ;
               0      010 ;
           m_shaft（N）/n*Omega^2 001]; 
    L_matrix= [1         L        （L^2）/（2*E*J（N））（L^3）/（6*E*J（N））;
               0         1           L/（E*J（N））      （L^2）/（2*E*J（N））;
               0         0           1                  L           ;