利用传递矩阵法和Riccati传递矩阵法分析转子临界转速

利用传递矩阵法和Riccati传递矩阵法分析转子临界转速一、所需求解转子参数将转子简化为如下所示：ABC1m2m3m1J2J3J三个盘的参数为：1232221232221233.5730.0160.050.0160.0120.0250.012PPPdddmkgmkgmkgIkgmIkgmIkgmIkgmIkgmIkgm另，阶梯轴的三段轴的截面惯性矩分别为：4142431.73.20.9JcmJcmJcm三段轴的单位长度轴段的质量分别为：1232.45/3.063/1.587/mkgmmkgmmkgm二、试算转轴的传递矩阵取试算转速1200/prads;则，各轴段的传递矩阵分别为：第1段840.061.7102.45/lmJmmkgm11.0006e+0006.0007e-0025.2943e-0071.0588e-0083.7356e-0021.0006e+0001.7649e-0055.2943e-0076.3506e+0031.2701e+0021.0006e+0006.0007e-0022.1170e+0056.3506e+0033.7356e-002H1.0006e+000第2段840.153.2103.063/lmJmmkgm21.0145e+0001.5044e-0011.7595e-0068.7927e-0083.8782e-0011.0145e+0002.3506e-0051.7595e-0064.9669e+0042.4821e+0031.0145e+0001.5044e-0016.6353e+0054.9669e+0043.8782e-001H1.0145e+000第3段840.053.2103.063/lmJmmkgm31.0002e+0005.0002e-0021.9531e-0073.2552e-0091.4358e-0021.0002e+0007.8128e-0061.9531e-0075.5135e+0039.1890e+0011.0002e+0005.0002e-0022.2054e+0055.5135e+0031.4358e-002H1.0002e+000第4段840.033.2103.063/lmJmmkgm41.0000e+0003.0000e-0027.0313e-0087.0313e-0103.1013e-0031.0000e+0004.6875e-0067.0313e-0081.9848e+0031.9848e+0011.0000e+0003.0000e-0021.3232e+0051.9848e+0033.1013e-003H1.0000e+000第5段840.10.9101.587/lmJmmkgm51.0053e+0001.0011e-0012.7788e-0069.2607e-0082.1163e-0011.0053e+0005.5614e-0052.7788e-0061.1430e+0043.8094e+0021.0053e+0001.0011e-0012.2877e+0051.1430e+0042.1163e-001H1.0053e+000第6段840.060.9101.587/lmJmmkgm61.0007e+0006.0008e-0021.0000e-0062.0000e-0084.5706e-0021.0007e+0003.3338e-0051.0000e-0064.1137e+0038.2272e+0011.0007e+0006.0008e-0021.3714e+0054.1137e+0034.5706e-002H1.0007e+000此6段传递矩阵均采用MATLAB编程求解，MATLAB的源文件为H.m三、采用传递矩阵法进行各段轴的状态参数的传递初始参数列阵为:0101010120101012011PddXXIIpMpIQmpx令011X，则初始矩阵可化为：0101157605.046e以初始矩阵乘第一轴段的传递矩阵，则可得第一段轴的终端状态参数：1011011011010.06306+1.0541.102+58903.088526566.05..7062556kkkkeeXMQ由于考虑支座的支撑刚度系数变化从5101*101*10，先取51*10，那么100001000010001KKk，此处510k，则可得支座A后第2段的起始端参数阵为：0201020102015602010.06306+1.0541.102+5890.03.088*10260.25.149*2.01076XMQ用第2段的传递矩阵乘此矩阵，可得第2段终端参数：201201201201660.2402+2.4721.282+19.11900.01.147*1099133.06.170477*1kkkkXMQ用中间圆盘的传递矩阵乘第2段终端参数阵，即可得第3段起始端参数：03010306103016307010.2402+2.4721.282+1958022.01.848*102.52*103.11*10.47XMQ用第3段传递矩阵乘其始端参数矩阵：'013'013'013'63560170.3238+3.9092.231+401.855*103.419*102.581*103.178*10.02kkkkXMQ用上式乘以支座刚度矩阵，得其终端参数：01305667130130130.3238+3.9092.231+40.01.855*103.419*102.549*103.1392*10kkkkXMQ则，根据可得：,则可得支座B后第4段的起始端参数阵为：560104010401040104670.3238+3.9092.231+40.01.855*103.419*102.549*103.139*102XMQ同上，用此段轴的传递函数乘其起始端的状态参数，可得：40140140140156670.4056+5.3723.281+582.626*104.369*102.597*103..172*20kkkkXMQ则，根据40kM可得：01-16.64则，可得第5段的起始参数矩阵：050750505-1.3753.69601.12*10XMQ其中，5为铰链处的转角。用第5段的传递矩阵乘此参数矩阵，即得第5段的终端参数：55555655750.1001*-2.051.005*-2380.9*1.135*1011433.07.*1.11753*0kkkkXMQ用上式乘以支座刚度矩阵可得第6段的初始状态参数阵：065065065076650.1001*-2.051.005*-2380.9*1.135*101420.0*1.17.733*10XMQ则，用第6段的传递矩阵乘此状态参数即可得其终端的参数阵：65656565670.1609-5.0751.025-7960.91.827*1019300.01.178*610.99kkkkXMQ根据最右边盘得传递矩阵，可得转子终端的状态参数：67665660.1609*b-5.0751.025*b-76863.0*2.27*107.144*10*3.371*106.99kkkkbbXMQ则根据终端的自由状态，则应该60kM;60kQ通过令60kQ解出cR，并将其带入到6kM的表达式中，可得：-1.9463e+006；此处使用的MATLAB源程序为calculate.m在MATLAB中使用线性插值法寻找最佳p值使得逼近于0。其程序为rotor.m经计算，考虑支撑刚度变化5101*101*10之间时，取51*10时，一阶临界转速值为1.0344e+002/rads取61*10时，一阶临界转速值为3.2390e+002/rads取71*10时，一阶临界转速值为9.3859e+002/rads取81*10时，一阶临界转速值为1.8631e+003/rads取91*10时，一阶临界转速值为2.1843e+003/rads取101*10时，一阶临界转速值为2.2241e+003/rads因此随着刚度的增加，一阶临界转速的值越来越大，而当不考虑支座的刚度变化，假设为完全刚性的话，一阶临界转速值为2228.7/rads，因此当取101*10时一阶临界转速值已相当接近完全刚性的情况。四、采用Riccati传递矩阵法进行各段轴的状态参数的传递根据Riccati传递矩阵法的原理，只需在传递矩阵法的基础之上求得各截面的Riccati传递矩阵。将转子截面的状态参数分组：iiMfQ,iiXe因为左端和右端均为自由端，故00f，00e；0nf，0ne;所以，我们可得到左端截面的Riccati传递矩阵00S；根据第i+1截面f、e之间的Riccati变换公式：1111221221iiiSTSTTST可得（同样，试算转速选为1200/prads）：左盘右边截面的Riccati传递矩阵：615.04*1005760.00S；第1轴段末的Riccati传递矩阵：52653.273*105.826*10-13403.0.100273*S；刚性支承5101*101*10在此处的处理因为涉及到刚度，取51*10的情况，所以在获取其支承左边的Riccati传递矩阵2S后，需将222fSe转换为1222eSf(即第二类Riccati变换)，然后再代入到普通传递矩阵的式子：32100001000100001fkfee32232ffKeee;000Kk可得：1113322eSIKSf。最终可得：(5)(7)(35)1320.0001871.072*104.3865*101.072*10eSf,即(5)(7)1(5)30.00018761.072*104.385*101.072*10S;此处的13S即为刚性支承右端截面的第二类Riccati传递矩阵。则，第2段轴段末的Riccati传递矩阵：11121221112433SHSHHSH可得：545658.213*101.632*103.97*108.213*10S通过中间圆盘后