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

1、利用传递矩阵法和riccati传递矩阵法分析转子临界转速一、所需求解转子参数 将转子简化为如下所示：m2m1m3abcjj1j23m13.5kgm27kgm33kg三个盘的参数为：i p1i d10.016kgm2 0.012kgm2i p 2i d 20.05kgm20.025kgm2i p3i d 30.016kgm20.012kgm2另，阶梯轴的三段轴的截面惯性矩分别为：j1.7cm412j3.2cm43j0.9cm4三段轴的单位长度轴段的质量分别为：m12.45kg / m m23.063kg / mm31.587kg / m二、试算转轴的传递矩阵取试算转速 p1200rad/ s ;

2、就，各轴段的传递矩阵分别为：第 1 段l0.06 mj1.7108 m4m2.45 kg / m1.0006e+000 6.0007e-002 5.2943e-007 1.0588e-0083.7356e-002 1.0006e+000 1.7649e-005 5.2943e-007h 16.3506e+003 1.2701e+002 1.0006e+000 6.0007e-0022.1170e+005 6.3506e+003 3.7356e-0021.0006e+000第 2 段l0.15mj3.2108 m4m3.063kg / m1.0145e+000 1.5044e-001 1.759

3、5e-006 8.7927e-0083.8782e-001 1.0145e+000 2.3506e-005 1.7595e-006h 24.9669e+004 2.4821e+003 1.0145e+000 1.5044e-001第 3 段l6.6353e+005 4.9669e+004 3.8782e-0010.05m1.0145e+000j3.2108 m4m3.063kg / m1.0002e+000 5.0002e-002 1.9531e-007 3.2552e-0091.4358e-002 1.0002e+000 7.8128e-006 1.9531e-007h 35.5135e+0

4、03 9.1890e+001 1.0002e+000 5.0002e-002第 4 段2.2054e+005 5.5135e+003 1.4358e-0021.0002e+000l0.03mj3.2108 m4m3.063kg / m1.0000e+000 3.0000e-002 7.0313e-008 7.0313e-0103.1013e-003 1.0000e+000 4.6875e-006 7.0313e-008h 41.9848e+003 1.9848e+001 1.0000e+000 3.0000e-002第 5 段1.3232e+005 1.9848e+003 3.1013e-00

5、31.0000e+000l0.1mj0.9108 m4m1.587kg / m1.0053e+000 1.0011e-001 2.7788e-006 9.2607e-0082.1163e-001 1.0053e+000 5.5614e-005 2.7788e-006h 51.1430e+004 3.8094e+002 1.0053e+000 1.0011e-001第 6 段2.2877e+005 1.1430e+004 2.1163e-0011.0053e+000l0.06mj0.9108 m4m1.587kg / m1.0007e+000 6.0008e-002 1.0000e-006 2.

6、0000e-0084.5706e-002 1.0007e+000 3.3338e-005 1.0000e-006h 64.1137e+003 8.2272e+001 1.0007e+000 6.0008e-0021.3714e+005 4.1137e+003 4.5706e-0021.0007e+000此 6 段传递矩阵均采纳matlab编程求解， matlab的源文件为 h.m三、采纳传递矩阵法 进行各段轴的状态参数的传递初始参数列阵为 :x 0101m 01q01x 01011i pp i dmp2 x2dip01011令 x 011 ，就初始矩阵可化为：015760015.04e6以初始

7、矩阵乘第一轴段的传递矩阵，就可得第一段轴的终端状态参数：x k10.0630601+ 1.054k11.10201 +2.706m k1qk15890 016566.0 013.088e55.255e6由于考虑支座的支撑刚度系数变化从1*10 51*101 0 ，先取1*105 ，那么100001000010k001kk，此处 k105 ，就可得支座 a 后第 2 段的起始端参数阵为：x 020.0630601 + 1.054021.10201 +2.706m 025890.0 013.088*10 5q260.25.149* 1060201用第 2 段的传递矩阵乘此矩阵，可得第2 段终端参数

8、：x k 20.240201+ 2.472k 2m k 21.282 0111900.0 01+ 19.471.147*10 6q99133.06.177*1 06k 201用中间圆盘的传递矩阵乘第2 段终端参数阵，即可得第3 段起始端参数：x 03030.24021.28201 + 2.47201 + 19.47m 03q58022.0 012.52*10 61.848*10 63.11*1070301用第 3 段传递矩阵乘其始端参数矩阵：x'k 30.3238 01 + 3.909k 3'2.231 01 + 40.02'm k 31.855*10 53.419*1

9、0 6q01k 301'2.581*1063.178*10 7用上式乘以支座刚度矩阵，得其终端参数：x k 3k 3m k 3qk 30.32382.2311.855*10 52.549*10 601 + 3.90901 + 40.02013.419*10 6013.139*10 7就，依据可得：, 就可得支座 b 后第 4 段的起始端参数阵为：x 0404m 04q040.32382.2311.855*10 52.549*10 601 + 3.90901 + 40.02013.419*10 6013.139*10 7同上，用此段轴的传递函数乘其起始端的状态参数，可得：x k 4k 4

10、m k 4qk 40.40563.2812.626*10 52.597*10 601 + 5.37201 + 58.27014.369*10 6013.2* 107就，依据m k 40 可得：01-16.64就，可得第 5 段的起始参数矩阵：x 0505m 05q05-1.3753.69601.12*10 7其中，5 为铰链处的转角；用第 5 段的传递矩阵乘此参数矩阵，即得第5 段的终端参数：x k 5k 50.1001*1.005*5 - 2.055 - 27.7m k 5qk 5380.9*511433.0*51.135*10 61.153* 107用上式乘以支座刚度矩阵可得第6 段的初始

11、状态参数阵：x 06060.1001*1.005*5 - 2.055 - 27.7m 06380.9*51.135*10 6q061420.0*1.133*10 75就，用第 6 段的传递矩阵乘此状态参数即可得其终端的参数阵：x k 6k 60.1609 51.025 5- 5.075- 76.99m k 6960.9 51.827*10 6qk 619300.0 51.178* 107依据最右边盘得传递矩阵，可得转子终端的状态参数：x k 6k 60.1609*b - 5.0751.025*b - 76.99m k 6qk 66863.0* b7.144*10 5 * b2.27*10 63

12、.371*10 7就依据终端的自由状态，就应当m k 60 ;qk 60通过令qk 60 解出rc ，并将其带入到m k 6 的表达式中，可得：-1.9463e+006；此处使用的 matlab源程序为 calculate.m在 matlab中使用线性插值法查找正确p 值使得靠近于 0；其程序为 rotor.m经运算，考虑支撑刚度变化1*10 51*10 10 之间时，取1*105 时，一阶临界转速值为 1.0344e+002rad / s取1*106 时，一阶临界转速值为3.2390e+002rad / s取1*107 时，一阶临界转速值为9.3859e+002rad / s取1*108 时

13、，一阶临界转速值为1.8631e+003rad / s9取1*10时，一阶临界转速值为2.1843e+003rad / s取1*10 10 时，一阶临界转速值为2.2241e+003rad / s因此随着刚度的增加， 一阶临界转速的值越来越大，而当不考虑支座的刚度变化，假设为完全刚性的话， 一阶临界转速值为 2228.7rad/ s ，因此当取1*10 10 时一阶临界转速值已相当接近完全刚性的情形；四、采纳 riccati传递矩阵法 进行各段轴的状态参数的传递依据 riccati传递矩阵法的原理， 只需在传递矩阵法的基础之上求得各截面的 riccati传递矩阵；将转子截面的状态参数分组：fi

14、m,exiqii由于左端和右端均为自由端，故f00 ，e 00 ；fn0 ， e n0 ;所以，我们可得到左端截面的riccati传递矩阵s 00 ；依据第 i+1 截面f、 e之间的 riccati变换公式：s i 1t11st12 i1t21st22 i可得（同样，试算转速选为p1200rad/ s ）：05760.0左盘右边截面的riccati传递矩阵：s 15.04*10 60；3.273*10 5-13400.0第 1 轴段末的 riccati传递矩阵：s 265；5.826*103.273* 10刚性支承 1*10 51*10 10 在此处的处理由于涉及到刚度，取1*105 的情形

15、，所以在猎取其支承左边的riccati传递矩阵s 2 后，需将f2s 2e转换为21e 2s 2f2 即其次类 riccati变换 ，然后再代入到一般传递矩阵的式子：1001000000fk0fe310e201f3f2ke 200; ke 3可得：e 3e 2k0111s 2iks 2f3 ；11.072*10 54 .385*107最终可得：e 3s20.00018761.072*10 5f3,即11.072*10 54.385*10 7s此处的1s即为刚性支承右端截面的其次30.0001876类 riccati传递矩阵；1.072*10 5;3就，第 2 段轴段末的 riccati传递矩阵

16、：1s 4h 21s 3h 221h 11 s 31h128.213*10 51.632*10 5可得： s 43.97*10 68.213*105通过中间圆盘后，可得圆盘右边的riccati传递矩阵：8.213*105s56.11*1061.992*10 5;8.213*10 5就，第 3轴段右边的 riccati传递矩阵为：-47877.090100.0s 67.642*1 0647877.0 ;到达其次个刚性支承处， 同样采纳其次类 riccati变换，并带入riccati传递矩阵公式：1111s i 1s iiks i1即可得：s7.019*10 81.322*10 771.106*1

17、0 57.019*10 8连续进行传递，第 4轴段末的 riccati传递矩阵为：1s 8h 21s 7h 221h11 s 71h 121.615*10 560100.0可得：s87.709*10 61.615*1 05由于，第 4轴段末是球头联轴器，故，在此也要进行另外的推导，由于球头联轴器的力矩刚性系数趋近于0，就在此利用弹性铰链的传递矩阵：1000e0110ek1fi 1f0010i0001e i 1e ikfi ； k001； kfi 1fik0k1f把is i带入上式中，可得：eis i 11s iiks i00就，可得：s 98.143* 1060再依据第 5轴段的传递矩阵，可得

18、第5轴段末的 riccati传递矩阵：1.589*10 61.586*10 5s 101.584*10 71.589*10 6第5轴段末为一刚性支承， 就同样采纳其次种 riccati变换，可得刚性支承右端的riccati传递矩阵：s10.00041854.177*1 0 5110.004199-0.0004185同样，依据第 6 轴段的传递矩阵，可得第6 轴段末的 riccati传递矩阵：2.501*10540211.0s121.473*1 062.501*10 5最终，通过最终一个圆盘的传递， 可以得到转子最右端的riccati传递矩阵：s 132.501*1052.847*10 6459

19、99.02.501* 105即，最右边截面上应当满意：f13s 13；又，由于e130f130；所以，只有当dets 130 时，才能取到e 13 的值，所以只有临界转速值p 才能使det s 13 尽量接近于 0；此处；dets 13-1.9345e+011此处的算法 matlab程序为 calculate2.m利用线性插值法来求取最终的临界转速值p，只需将第一次求解的rotor.m的程序中的算法程序calculate.m改为 calculate2.m即可；考虑支撑刚度变化1*10 51*10 10 之间时，取1*10 5 时，一阶临界转速值为1.0344e+002rad / s取1*106

20、 时，一阶临界转速值为3.2390e+002rad / s取1*107 时，一阶临界转速值为9.3859e+002rad / s8取1*10时，一阶临界转速值为1.8631e+003rad / s取1*109 时，一阶临界转速值为2.1843e+003rad / s取1*10 10 时，一阶临界转速值为2.2241e+003rad / s因此随着刚度的增加， 一阶临界转速的值越来越大，而当不考虑支座的刚度变化，假设为完全刚性的话， 一阶临界转速值为 2228.7rad/ s ，因此当取1*10 10 时一阶临界转速值已相当接近完全刚性的情形；附录：h.mfunctionoutput = hin

21、1,in2,in3,in4% in1=l in2=j in3=m in4=psymss tu v kae jm pl ; e=2e11;l=in1; j=in2; m=in3; p=in4;k=l*m*p2/e*j0.25; s=0.5*0.5*expk+exp-k+cosk;t=0.5*0.5*expk-exp-k+sink;u=0.5*0.5*expk+exp-k-cosk;v=0.5*0.5*expk-exp-k-sink; a=k/l;output=s t/a u/a2*e*j v/a3*e*j;a*v s t/a*e*j u/a2*e*j;a2*e*j*u a*e*j*v s t/a

22、;a3*e*j*t a2*e*j*u a*v s;calculate.mfunctionoutput=calculatein symspabp=in;h1=h0.06,1.7e-8,2.45,p;h2=h0.15,3.2e-8,3.063,p;h3=h0.05,3.2e-8,3.063,p;h4=h0.03,3.2e-8,3.063,p;h5=h0.1,0.9e-8,1.587,p;h6=h0.06,0.9e-8,1.587,p;o1=1 a -1-0.016/0.012*0.012*p2*a 3.5*p2'o1=h1*o1;kk=1 0 0 0;0 1 0 0;0 0 1 0;-10

23、8 0 0 1;o2=kk*o1; o2=h2*o2;k1=1 0 0 0;0 1 0 0;0 0.025*p2 1 0;7*p2 0 0 1;o3=k1*o2; o3=h3*o3; o4=kk*o3; o4=h4*o4; d=solveo43;o5=subso4,'a',d;o52=o52+b;o5=h5*o5; o6=kk*o5; o6=h6*o6;k2=1 0 0 0;0 1 0 0;0 -1-0.016/0.012*0.012*p2 1 0;3*p2 0 0 1;oo=k2*o6;rc=solveoo4;ee=subsoo3,'b',rc;output=

24、doubleee; ee=vpaee;calculate2.mfunctionoutput=calculate2in symspk ;p=in;h1=h0.06,1.7e-8,2.45,p;h2=h0.15,3.2e-8,3.063,p;h3=h0.05,3.2e-8,3.063,p;h4=h0.03,3.2e-8,3.063,p;h5=h0.1,0.9e-8,1.587,p;h6=h0.06,0.9e-8,1.587,p; s=0 0;0 0;k1=10 0 0;01 0 0;0-1-0.016/0.012*0.012*p21 0;3.5*p200 1; t22=k11,1:2;k12,1:

25、2;t21=k11,3:4;k12,3:4;t12=k13,1:2;k14,1:2;t11=k13,3:4;k14,3:4; s1=t11*s+t12*invt21*s+t22; t22=h11,1:2;h12,1:2;t21=h11,3:4;h12,3:4;t12=h13,1:2;h14,1:2;t11=h13,3:4;h14,3:4; s2=t11*s1+t12*invt21*s1+t22;is2=invs2;kk=0 0;-105 0;i=eye2; is3=is2*invi+kk*is2; t11=h21,1:2;h22,1:2;t12=h21,3:4;h22,3:4;t21=h23,

26、1:2;h24,1:2;t22=h23,3:4;h24,3:4; s4=t21*is3+t22*invt11*is3+t12;k2=1 0 0 0;0 1 0 0;0 0.025*p2 1 0;7*p2 0 0 1;t22=k21,1:2;k22,1:2;t21=k21,3:4;k22,3:4;t12=k23,1:2;k24,1:2;t11=k23,3:4;k24,3:4; s5=t11*s4+t12*invt21*s4+t22; t22=h31,1:2;h32,1:2;t21=h31,3:4;h32,3:4;t12=h33,1:2;h34,1:2;t11=h33,3:4;h34,3:4; s6=t11*s5+t12*invt21*s5+t22;is6=invs6; is7=is6*invi+kk*is6; t11=h41,1:2;h42,1:2;t12=h41,3:4;h42,3: