利用傳遞矩陣法和Riccati傳遞矩陣法分析轉(zhuǎn)子臨界轉(zhuǎn)速

利用傳遞矩陣法和Riccati傳遞矩陣法分析轉(zhuǎn)子臨界轉(zhuǎn)速所需求解轉(zhuǎn)子參數(shù)將轉(zhuǎn)子簡化為如下所示：三個盤的參數(shù)為：QUOTEm1=m2=另，階梯軸的三段軸的截面慣性矩分別為：QUOTEJ1=1.7cm4J2三段軸的單位長度軸段的質(zhì)量分別為：QUOTEm1=2.47kgmm試算轉(zhuǎn)軸的傳遞矩陣取試算轉(zhuǎn)速Q(mào)UOTEp=ω=1000rads;則，各軸段的傳遞矩陣分別為：第1段QUOTEl=0.06mJ=1.7×10-8m4m=2.47第2段QUOTEl=0.15mJ=3.2×10-8m4第3段第4段第5段第6段此6段傳遞矩陣均采用MATLAB編程求解，MATLAB的源文件為HYPERLINKcalculate.m在MATLAB中使用線性插值法尋找最佳p值使得?逼近于0。其程序為rotor.m經(jīng)計算，考慮支撐剛度變化之間時，取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為因此隨著剛度的增加，一階臨界轉(zhuǎn)速的值越來越大，而當(dāng)不考慮支座的剛度變化，假設(shè)為完全剛性的話，一階臨界轉(zhuǎn)速值為，因此當(dāng)取時一階臨界轉(zhuǎn)速值已相當(dāng)接近完全剛性的情況。采用Riccati傳遞矩陣法進(jìn)行各段軸的狀態(tài)參數(shù)的傳遞根據(jù)Riccati傳遞矩陣法的原理，只需在傳遞矩陣法的基礎(chǔ)之上求得各截面的Riccati傳遞矩陣。將轉(zhuǎn)子截面的狀態(tài)參數(shù)分組：,QUOTEfi={MQ}因為左端和右端均為自由端，故QUOTEf0=0，QUOTEe0≠0；QUOTEfN=0，QUOTEeN≠0;所以，我們可得到左端截面的Riccati傳遞矩陣QUOTES0=0；根據(jù)第i+1截面、QUOTEf和e之間的Riccati變換公式：可得（同樣，試算轉(zhuǎn)速選為QUOTEp=1000rads）：左盤右邊截面的Riccati傳遞矩陣：QUOTES1=0250007e60第1軸段末的Riccati傳遞矩陣：QUOTES2=4.063e5-3.23e38.07e6-4.063e5剛性支承在此處的處理因為涉及到剛度，取QUOTE∞的情況，所以在獲取其支承左邊的Riccati傳遞矩陣后，需將QUOTEf2=S2e2轉(zhuǎn)換為QUOTEe2=S2-1f2(即第二類Riccati變換)，然后再代入到普通傳遞矩陣的式子：;可得：。最終可得：,即QUOTE{e}3=00-3.0962e-40f3，即S3-1=00則，第2段軸段末的Riccati傳遞矩陣：可得：QUOTE[S]4=-8.309e51.270e5-5.134e6通過中間圓盤后，可得圓盤右邊的Riccati傳遞矩陣：;則，第3軸段右邊的Riccati傳遞矩陣為：;到達(dá)第二個剛性支承處，同樣采用第二類Riccati變換，并帶入Riccati傳遞矩陣公式：即可得：QUOTES7-1=00繼續(xù)進(jìn)行傳遞，第4軸段末的Riccati傳遞矩陣為：可得：QUOTE[S]8=-2.337e76.604e5-9.146e8由于，第4軸段末是球頭聯(lián)軸器，故，在此也要進(jìn)行另外的推導(dǎo)，由于球頭聯(lián)軸器的力矩剛性系數(shù)趨近于0，則在此利用彈性鉸鏈的傳遞矩陣：；；把QUOTEfi=SieiQUOTESi+1=Si(I+K則，可得：QUOTES9=00-8.765e7再根據(jù)第5軸段的傳遞矩陣，可得第5軸段末的Riccati傳遞矩陣：第5軸段末為一剛性支承，則同樣采用第二種Riccati變換，可得剛性支承右端的Riccati傳遞矩陣：同樣，根據(jù)第6軸段的傳遞矩陣，可得第6軸段末的Riccati傳遞矩陣：最后，通過最后一個圓盤的傳遞，可以得到轉(zhuǎn)子最右端的Riccati傳遞矩陣：即，最右邊截面上應(yīng)該滿足：QUOTEf13=S13e13；又，由于QUOTEf13={00}；所以，只有當(dāng)QUOTEdetS13=0時，才能取到QUOTEe13的值，所以只有臨界轉(zhuǎn)速值p才能使QUOTEdetS13盡量接近于0。此處QUOTEdetS13=1.2e12；此處的算法MATLAB程序為calculate2.m利用線性插值法來求取最終的臨界轉(zhuǎn)速值p，只需將第一次求解的rotor.m的程序中的算法程序calculate.m改為calculate2.m即可?？紤]支撐剛度變化之間時，取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為取時，一階臨界轉(zhuǎn)速值為因此隨著剛度的增加，一階臨界轉(zhuǎn)速的值越來越大，而當(dāng)不考慮支座的剛度變化，假設(shè)為完全剛性的話，一階臨界轉(zhuǎn)速值為，因此當(dāng)取時一階臨界轉(zhuǎn)速值已相當(dāng)接近完全剛性的情況。附錄：H.mfunctionoutput=H(in1,in2,in3,in4)%in1=lin2=Jin3=min4=psymsSTUVkaEJmpl;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=[ST/aU/(a^2*E*J)V/(a^3*E*J);a*VST/(a*E*J)U/(a^2*E*J);a^2*E*J*Ua*E*J*VST/a;a^3*E*J*Ta^2*E*J*Ua*VS];Calculate.mfunctionoutput=calculate(in)symspabp=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=[1a-(1-0.016/0.012)*0.012*p^2*a3.5*p^2]';o1=H1*o1;KK=[1000;0100;0010;-10^8001];o2=KK*o1;o2=H2*o2;K1=[1000;0100;00.025*p^210;7*p^2001];o3=K1*o2;o3=H3*o3;o4=KK*o3;o4=H4*o4;d=solve(o4(3));o5=subs(o4,'a',d);o5(2)=o5(2)+b;o5=H5*o5;o6=KK*o5;o6=H6*o6;K2=[1000;0100;0-(1-0.016/0.012)*0.012*p^210;3*p^2001];oo=K2*o6;RC=solve(oo(4));ee=subs(oo(3),'b',RC);output=double(ee);ee=vpa(ee);Calculate2.mfunctionoutput=calculate2(in)symspk;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);S=[00;00];K1=[1000;0100;0-(1-0.016/0.012)*0.012*p^210;3.5*p^2001];T22=[K1(1,1:2);K1(2,1:2)];T21=[K1(1,3:4);K1(2,3:4)];T12=[K1(3,1:2);K1(4,1:2)];T11=[K1(3,3:4);K1(4,3:4)];S1=(T11*S+T12)*inv(T21*S+T22);T22=[H1(1,1:2);H1(2,1:2)];T21=[H1(1,3:4);H1(2,3:4)];T12=[H1(3,1:2);H1(4,1:2)];T11=[H1(3,3:4);H1(4,3:4)];S2=(T11*S1+T12)*inv(T21*S1+T22);iS2=inv(S2);KK=[00;-10^50];I=eye(2);iS3=iS2*inv(I+KK*iS2);T11=[H2(1,1:2);H2(2,1:2)];T12=[H2(1,3:4);H2(2,3:4)];T21=[H2(3,1:2);H2(4,1:2)];T22=[H2(3,3:4);H2(4,3:4)];S4=(T21*iS3+T22)*inv(T11*iS3+T12);K2=[1000;0100;00.025*p^210;7*p^2001];T22=[K2(1,1:2);K2(2,1:2)];T21=[K2(1,3:4);K2(2,3:4)];T12=[K2(3,1:2);K2(4,1:2)];T11=[K2(3,3:4);K2(4,3:4)];S5=(T11*S4+T12)*inv(T21*S4+T22);T22=[H3(1,1:2);H3(2,1:2)];T21=[H3(1,3:4);H3(2,3:4)];T12=[H3(3,1:2);H3(4,1:2)];T11=[H3(3,3:4);H3(4,3:4)];S6=(T11*S5+T12)*inv(T21*S5+T22);iS6=inv(S6);iS7=iS6*inv(I+KK*iS6);T11=[H4(1,1:2);H4(2,1:2)];T12=[H4(1,3:4);H4(2,3:4)];T21=[H4(3,1:2);H4(4,1:2)];T22=[H4(3,3:4);H4(4,3:4)];S8=(T21*iS7+T22)*inv(T11*iS7+T12);KKK=[00;1/k0];S9=S8*inv(I+KKK*S8);S9=limit(S9,k,0);T22=[H5(1,1:2);H5(2,1:2)];T21=[H5(1,3:4);H5(2,3:4)];T12=[H5(3,1:2);H5(4,1:2)];T11=[H5(3,3:4);H5(4,3:4)];S10=(T11*S9+T12)*inv(T21*S9+T22);iS10=inv(S10);iS11=iS10*inv(I+KK*iS10);T11=[H6(1,1:2);H6(2,1:2)];T12=[H6(1,3:4);H6(2,3:4)];T21=[H6(3,1:2);H6(4,1:2)];T22=[H6(3,3:4);H6(4,3:4)];S12=(T21*iS11+T22)*inv(T11*iS11+T12);K3=[1000;0100;0-(1-0.016/0.012)*0.012*p^210;3*p^2001];T22=[K3(1,1:2);K3(2,1:2)];T21=[K3(1,3:4);K3(2,3:4)];T12=[K3(3,1:2);K3(4,1:2)];T11=[K3(3,3:4);K3(4,3:4)];S13=(T11*S12+T12)*inv(T21*S12+T22);ee=det(S13);output=double(ee);rotor.msymsab;a=1200;b=2300;a1=a+0.5*(b-a);f1=calculate(a);f2=calculate(b);fa=calculate(a1);fori=1:10000ifabs(fa)<1e-2disp('thecriticalspeedis');disp(a1);break;elseiff1*fa<0f2=fa;b=a1;a1=a