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

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

利用传递矩阵法和Riccati传递矩阵法分析转子临界转速 一、 所需求解转子参数 将转子简化为如下所示： 三个盘的参数为： 另，阶梯轴的三段轴的截面惯性矩分别为： 三段轴的单位长度轴段的质量分别为： 二、 试算转轴的传递矩阵 取试算转速 ; 则，各轴段的传递矩阵分别为： 第1段 第2段 第3段 第4段 第5段 第6段 此6段传递矩阵均采用MATLAB编程求解，MATLAB的源文件为H.m 三、 采用传递矩阵法进行各段轴的状态参数的传递 初始参数列阵为: 令，则初始矩阵可化为： 以初始矩阵乘第一轴段的传递矩阵，则可得第一段轴的终端状态参数： 由于考虑支座的支撑刚度系数变化从，先取，那么，此处，则可得支座A后第2段的起始端参数阵为： 用第2段的传递矩阵乘此矩阵，可得第2段终端参数： 用中间圆盘的传递矩阵乘第2段终端参数阵，即可得第3段起始端参数： 用第3段传递矩阵乘其始端参数矩阵： 用上式乘以支座刚度矩阵，得其终端参数： 则，根据可得： ,则可得支座B后第4段的起始端参数阵为： 同上，用此段轴的传递函数乘其起始端的状态参数，可得： 则，根据可得： 则，可得第5段的起始参数矩阵： 其中，为铰链处的转角。 用第5段的传递矩阵乘此参数矩阵，即得第5段的终端参数： 用上式乘以支座刚度矩阵可得第6段的初始状态参数阵： 则，用第6段的传递矩阵乘此状态参数即可得其终端的参数阵： 根据最右边盘得传递矩阵，可得转子终端的状态参数： 则根据终端的自由状态，则应该; 通过令解出，并将其带入到的表达式中，可得： ； 此处使用的MATLAB源程序为calculate.m 在MATLAB中使用线性插值法寻找最佳p值使得∆逼近于0。其程序为rotor.m 经计算，考虑支撑刚度变化之间时，取时，一阶临界转速值为 取时，一阶临界转速值为 取时，一阶临界转速值为 取时，一阶临界转速值为 取时，一阶临界转速值为 取时，一阶临界转速值为 因此随着刚度的增加，一阶临界转速的值越来越大，而当不考虑支座的刚度变化，假设为完全刚性的话，一阶临界转速值为，因此当取时一阶临界转速值已相当接近完全刚性的情况。 四、 采用Riccati传递矩阵法进行各段轴的状态参数的传递 根据Riccati传递矩阵法的原理，只需在传递矩阵法的基础之上求得各截面的Riccati传递矩阵。 将转子截面的状态参数分组：, 因为左端和右端均为自由端，故，；，; 所以，我们可得到左端截面的Riccati传递矩阵； 根据第i+1截面、之间的Riccati变换公式： 可得(同样，试算转速选为)： 左盘右边截面的Riccati传递矩阵：； 第1轴段末的Riccati传递矩阵：； 刚性支承在此处的处理因为涉及到刚度，取的情况，所以在获取其支承左边的Riccati传递矩阵后，需将转换为(即第二类Riccati变换)，然后再代入到普通传递矩阵的式子： ; 可得：。 最终可得：,即;此处的即为刚性支承右端截面的第二类Riccati传递矩阵。 则，第2段轴段末的Riccati传递矩阵： 可得： 通过中间圆盘后，可得圆盘右边的Riccati传递矩阵： ; 则，第3轴段右边的Riccati传递矩阵为： ; 到达第二个刚性支承处，同样采用第二类Riccati变换，并带入Riccati传递矩阵公式： 即可得： 继续进行传递，第4轴段末的Riccati传递矩阵为： 可得： 由于，第4轴段末是球头联轴器，故，在此也要进行另外的推导，由于球头联轴器的力矩刚性系数趋近于0，则在此利用弹性铰链的传递矩阵： ；； 把带入上式中，可得： 则，可得： 再根据第5轴段的传递矩阵，可得第5轴段末的Riccati传递矩阵： 第5轴段末为一刚性支承，则同样采用第二种Riccati变换，可得刚性支承右端的Riccati传递矩阵： 同样，根据第6轴段的传递矩阵，可得第6轴段末的Riccati传递矩阵： 最后，通过最后一个圆盘的传递，可以得到转子最右端的Riccati传递矩阵： 即，最右边截面上应该满足：；又，由于；所以，只有当时，才能取到的值，所以只有临界转速值p才能使尽量接近于0。 此处； 此处的算法MATLAB程序为calculate2.m 利用线性插值法来求取最终的临界转速值p，只需将第一次求解的rotor.m的程序中的算法程序calculate.m改为calculate2.m即可。 考虑支撑刚度变化之间时，取时，一阶临界转速值为 取时，一阶临界转速值为 取时，一阶临界转速值为 取时，一阶临界转速值为 取时，一阶临界转速值为 取时，一阶临界转速值为 因此随着刚度的增加，一阶临界转速的值越来越大，而当不考虑支座的刚度变化，假设为完全刚性的话，一阶临界转速值为，因此当取时一阶临界转速值已相当接近完全刚性的情况。 附录： H.m function output = H(in1,in2,in3,in4) % in1=l in2=J in3=m in4=p syms S T U V k a E J m p l; E=2e11; l=in1; J=in2; m=in3; p=in4; k=l*(m*(p^2)/(E*J))^0.25; S=0.5*(0.5*(exp(k)+exp(-k))+cos(k)); T=0.5*(0.5*(exp(k)-exp(-k))+sin(k)); U=0.5*(0.5*(exp(k)+exp(-k))-cos(k)); V=0.5*(0.5*(exp(k)-exp(-k))-sin(k)); a=k/l; output=[S T/a U/(a^2*E*J) V/(a^3*E*J);a*V S T/(a*E*J) U/(a^2*E*J);a^2*E*J*U a*E*J*V S T/a;a^3*E*J*T a^2*E*J*U a*V S]; Calculate.m function output=calculate(in) syms p a b p=in; H1=H(0.06,1.7e-8,2.45,p); H2=H(0.15,3.2e-8,3.063,p); H3=H(0.05,3.2e-8,3.063,p); H4=H(0.03,3.2e-8,3.063,p); H5=H(0.1,0.9e-8,1.587,p); H6=H(0.06,0.9e-8,1.587,p); o1=[1 a -(1-0.016/0.012)*0.012*p^2*a 3.5*p^2] ; o1=H1*o1; KK=[1 0 0 0;0 1 0 0;0 0 1 0;-10^8 0 0 1]; o2=KK*o1; o2=H2*o2; K1=[1 0 0 0;0 1 0 0;0