雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)實(shí)體有限元?jiǎng)恿Ψ治?/h1>
2014-05-16 07:01:40馬威猛王建軍
振動(dòng)與沖擊 2014年23期
馬威猛，王建軍
（北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京 100191）
雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)實(shí)體有限元?jiǎng)恿Ψ治?/p>
馬威猛，王建軍
（北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京 100191）
研究了同時(shí)考慮轉(zhuǎn)子和支承非對(duì)稱特征的轉(zhuǎn)子支承系統(tǒng)3D有限元建模和振動(dòng)特性分析方法。首先用ANSYS得到轉(zhuǎn)子在旋轉(zhuǎn)坐標(biāo)系下的整體質(zhì)量、剛度和科氏力矩陣；其次利用動(dòng)靜坐標(biāo)變換關(guān)系得到旋轉(zhuǎn)坐標(biāo)系下的時(shí)變支承剛度；最后利用節(jié)點(diǎn)編號(hào)與矩陣行號(hào)的對(duì)應(yīng)關(guān)系在整體矩陣中施加約束和時(shí)變支承剛度，得到系統(tǒng)的時(shí)變系數(shù)運(yùn)動(dòng)微分方程?；诟ヂ鍎P理論和赫爾無窮行列式方法對(duì)方程求解得到系統(tǒng)的頻率特性和穩(wěn)定性。利用文獻(xiàn)算例模型對(duì)該方法的有效性進(jìn)行了驗(yàn)證，并對(duì)實(shí)際工程轉(zhuǎn)子在考慮雙重非對(duì)稱特征時(shí)的振動(dòng)特性進(jìn)行了分析。
非對(duì)稱；轉(zhuǎn)子支承系統(tǒng)；3D有限元；時(shí)變系統(tǒng)；穩(wěn)定性
轉(zhuǎn)子支承系統(tǒng)是旋轉(zhuǎn)機(jī)械的核心部件，其振動(dòng)特性對(duì)旋轉(zhuǎn)機(jī)械的工作性能和可靠性具有重要的影響。轉(zhuǎn)子支承系統(tǒng)振動(dòng)特性研究對(duì)旋轉(zhuǎn)機(jī)械的結(jié)構(gòu)設(shè)計(jì)、維修保養(yǎng)、振動(dòng)控制、故障診斷等方面工作都具有重要的意義。
根據(jù)轉(zhuǎn)子與支承的軸對(duì)稱特征，轉(zhuǎn)子支承系統(tǒng)可以分為：軸對(duì)稱（axisymmetric）、支承非對(duì)稱（anisotropic）、轉(zhuǎn)子非對(duì)稱（asymmetric）和雙重非對(duì)稱（anisotropic asymmetric）轉(zhuǎn)子支承系統(tǒng)，其中雙重非對(duì)稱系統(tǒng)又稱為一般（general）轉(zhuǎn)子支承系統(tǒng)［1］。人們對(duì)軸對(duì)稱轉(zhuǎn)子支承系統(tǒng)開展了充分的研究工作，基于3D實(shí)體有限元的分析模型可以精確描述復(fù)雜截面轉(zhuǎn)子的振動(dòng)特性［2－4］。而在工程實(shí)際中，因加工、設(shè)計(jì)或故障破壞等原因，轉(zhuǎn)子或支承部件通常會(huì)包含一些非對(duì)稱特征，典型的如發(fā)電機(jī)轉(zhuǎn)子、裂紋轉(zhuǎn)子等。轉(zhuǎn)靜子上存在的非對(duì)稱特征使得轉(zhuǎn)子表現(xiàn)出獨(dú)特的振動(dòng)特性，需要進(jìn)行專門的研究，如參數(shù)振動(dòng)特性是雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的主要特征［5］。
20世紀(jì)60年代以來，國外學(xué)者針對(duì)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的振動(dòng)特性開展了廣泛的定性［6－8］、定量研究［9－12］。在建模分析方法方面，Nandi［12－14］在不同的框架下提出并開展了基于3D實(shí)體有限元的非對(duì)稱轉(zhuǎn)子系統(tǒng)動(dòng)力特性分析工作，Kim［15］仍在1D有限元的框架下考慮非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的部件耦合問題。沈松等［16－19］在非對(duì)稱轉(zhuǎn)子系統(tǒng)的振動(dòng)分析方面也做了出色的工作，但主要集中在非線性分析領(lǐng)域，對(duì)線性領(lǐng)域內(nèi)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的建模和振動(dòng)特性分析方法沒有加以關(guān)注。
3D實(shí)體單元建模不需對(duì)轉(zhuǎn)靜子結(jié)構(gòu)進(jìn)行過多簡(jiǎn)化，可以真實(shí)反映轉(zhuǎn)靜子的結(jié)構(gòu)特征，并且能夠準(zhǔn)確考慮轉(zhuǎn)靜子各部件之間的耦合效應(yīng)，是轉(zhuǎn)子動(dòng)力學(xué)建模方法的發(fā)展趨勢(shì)。目前，多數(shù)商用有限元分析軟件已支持基于3D實(shí)體單元建模的轉(zhuǎn)子動(dòng)力學(xué)分析，但其分析功能僅能考慮軸對(duì)稱或單一非對(duì)稱轉(zhuǎn)子支承系統(tǒng)，而不能進(jìn)行雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的振動(dòng)特性計(jì)算。
本文利用商用有限元軟件ANSYS的前處理功能，得到轉(zhuǎn)子3D有限元模型的質(zhì)量、剛度和科氏力矩陣，通過矩陣處理加入支承結(jié)構(gòu)在旋轉(zhuǎn)坐標(biāo)描述下的時(shí)變剛度，得到轉(zhuǎn)子支承系統(tǒng)的時(shí)變運(yùn)動(dòng)微分方程，然后基于時(shí)變系統(tǒng)分析的弗洛凱理論和無窮行列式方法將時(shí)變方程轉(zhuǎn)化為一定截?cái)嚯A次的時(shí)不變線性方程組，求解得到系統(tǒng)的頻率特性和穩(wěn)定性，實(shí)現(xiàn)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的3D有限元振動(dòng)特性分析功能。文中首先介紹相關(guān)的分析理論，包括坐標(biāo)系的選擇，運(yùn)動(dòng)微分方程的建立以及基于弗洛凱理論的時(shí)變參數(shù)系統(tǒng)的求解方法；然后以兩個(gè)算例對(duì)文中雙重非對(duì)稱轉(zhuǎn)子振動(dòng)特性分析方法的有效性和工程應(yīng)用價(jià)值進(jìn)行驗(yàn)證和說明。
1 坐標(biāo)系選擇
工程實(shí)際中，裂紋、鍵槽、非圓截面等非對(duì)稱轉(zhuǎn)子特征及各種形式的支承非對(duì)稱特征可以抽象為圖1所示的典型雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)［20］。圖1中的非圓截面軸代表前述多種因素引起的轉(zhuǎn)子剛度非對(duì)稱特征，非圓截面盤代表轉(zhuǎn)子的非對(duì)稱轉(zhuǎn)動(dòng)慣量特征，正交方向上不同的支承剛度和阻尼則表示支承系統(tǒng)的非對(duì)稱特征。
圖1 雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)示意圖Fig．1 Simple analysis model of anisotropic asymmetric rotor bearing system
對(duì)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)而言，采用旋轉(zhuǎn)或者固定坐標(biāo)系描述系統(tǒng)運(yùn)動(dòng)不像單一非對(duì)稱轉(zhuǎn)子支承系統(tǒng)那樣具有特殊的意義－選用合適的坐標(biāo)系可以把時(shí)變微分方程轉(zhuǎn)化為時(shí)不變微分方程的求解問題。雙重非對(duì)稱轉(zhuǎn)子系統(tǒng)在固定或轉(zhuǎn)動(dòng)坐標(biāo)系下描述均不能得到具有固定系數(shù)的運(yùn)動(dòng)微分方程。雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)在選取分析坐標(biāo)系時(shí)，考慮到：①目前ANSYS等商用有限元分析軟件僅支持非對(duì)稱轉(zhuǎn)子在旋轉(zhuǎn)坐標(biāo)系下的3D實(shí)體有限元分析；②若將靜子系統(tǒng)的支承剛度以參數(shù)形式給出，在旋轉(zhuǎn)坐標(biāo)系下表達(dá)的時(shí)變剛度系數(shù)僅涉及支承點(diǎn)附近的幾個(gè)相關(guān)自由度，表達(dá)形式簡(jiǎn)單。因此，一般選用隨轉(zhuǎn)子同速旋轉(zhuǎn)的旋轉(zhuǎn)坐標(biāo)系對(duì)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的運(yùn)動(dòng)進(jìn)行描述。
2 旋轉(zhuǎn)坐標(biāo)系下的運(yùn)動(dòng)微分方程
在旋轉(zhuǎn)坐標(biāo)系下，采用3D實(shí)體單元對(duì)非對(duì)稱轉(zhuǎn)子進(jìn)行建模能夠準(zhǔn)確反映非對(duì)稱幾何特征或材料特征對(duì)轉(zhuǎn)子截面彎曲剛度或轉(zhuǎn)動(dòng)慣量的影響。實(shí)體非對(duì)稱轉(zhuǎn)子在邊界節(jié)點(diǎn)附近由剛度單元支承，剛度單元在各個(gè)支承方向的剛度差異反映支承系統(tǒng)的非對(duì)稱特征。為方便描述，以下標(biāo)nb標(biāo)示轉(zhuǎn)子的支承節(jié)點(diǎn)。在固定坐標(biāo)系下，這些支點(diǎn)的剛度可以表示為［13］
轉(zhuǎn)子支點(diǎn)在旋轉(zhuǎn)坐標(biāo)系下的支承剛度可由式（1）經(jīng)坐標(biāo)變換得到
旋轉(zhuǎn)坐標(biāo)系下，無阻尼轉(zhuǎn)子支承系統(tǒng)在外力作用下的運(yùn)動(dòng)微分方程可以寫為
式中［Kr］是轉(zhuǎn)子支承系統(tǒng)的剛度矩陣，由轉(zhuǎn)子自身剛度和支承剛度組成，［Mr］是轉(zhuǎn)子系統(tǒng)在旋轉(zhuǎn)坐標(biāo)系下的質(zhì)量矩陣，［Cr］被稱為科氏力矩陣，是反對(duì)稱矩陣，｛U｝是轉(zhuǎn)子系統(tǒng)在旋轉(zhuǎn)坐標(biāo)系下的位移向量，｛f｝則是在旋轉(zhuǎn)坐標(biāo)系下表示的外部載荷向量。
在ANSYS前處理模塊中完成轉(zhuǎn)子結(jié)構(gòu)建模后，若不考慮支承的非對(duì)稱特征，利用ANSYS的轉(zhuǎn)子動(dòng)力學(xué)分析模塊可以直接得到方程（3）中各矩陣的元素?cái)?shù)值；而考慮支承非對(duì)稱時(shí)，剛度矩陣是時(shí)變的，需要進(jìn)一步編程處理才能得到相應(yīng)的矩陣。
考慮到支承剛度在旋轉(zhuǎn)坐標(biāo)系下的時(shí)變特征，方程（3）可以改寫為
式中［Kr］nb是不考慮支承剛度時(shí)，轉(zhuǎn)子系統(tǒng)在旋轉(zhuǎn)坐標(biāo)系下的剛度矩陣，可由ANSYS前處理模塊直接得到；［K0］b是支承剛度的非時(shí)變項(xiàng)，而［Kc］b和［Ks］b則是支承剛度時(shí)變項(xiàng)中的余弦和正弦分量，支承剛度各矩陣由式（2）經(jīng)坐標(biāo)變換得到。
考慮到三角函數(shù)與復(fù)數(shù)的變換關(guān)系，方程（4）可寫成復(fù)數(shù)形式，
3 時(shí)變參數(shù)系統(tǒng)的求解方法
時(shí)變參數(shù)系統(tǒng)的求解方法與一般時(shí)不變系統(tǒng)不同。弗洛凱理論和赫爾無窮行列式方法是分析時(shí)變參數(shù)系統(tǒng)常用的分析理論。
基于方程（5）的齊次方程，可對(duì)轉(zhuǎn)子支承系統(tǒng)的自由振動(dòng)特性進(jìn)行分析。根據(jù)弗洛凱理論，假設(shè)方程解的形式為
其中，Λ是系統(tǒng)的基礎(chǔ)特征頻率，Φ是與Λ相對(duì)應(yīng)的類模態(tài)矩陣。之所以稱Φ為類模態(tài)，是因?yàn)槠湫问诫m然與時(shí)不變系統(tǒng)的模態(tài)矩陣相似，但Φ具有時(shí)變的性質(zhì)，其時(shí)變周期與運(yùn)動(dòng)微分方程中時(shí)變系統(tǒng)的周期一致。
其中，lj為Φ的第j階階次展開模態(tài)。
將式（6）和式（7）代入方程（5）的齊次方程中，并整理可以得到
方程（8）中，令指數(shù)項(xiàng)對(duì)應(yīng)各階次系數(shù)等于零可以得到
式（9）在截?cái)嚯A次j＝∞時(shí)，與原時(shí)變微分方程是等價(jià)的，但在實(shí)際的計(jì)算中，需要進(jìn)行適當(dāng)階次的截?cái)?，得到原時(shí)變微分方程解的近似值。求解得到各展開階次的頻率及特征向量后，根據(jù)式（7）可得到原時(shí)變微分運(yùn)動(dòng)方程的近似特征值Λ和類模態(tài)Φ。
式（9）可寫成特征值求解問題的一般形式
上述特征值問題一般轉(zhuǎn)換到狀態(tài)空間進(jìn)行求解。式（10）在狀態(tài)空間中可以寫為
在狀態(tài)空間中求解式（11）的特征值問題，可得到相應(yīng)的特征頻率λ和特征向量Ψ，進(jìn)而得到轉(zhuǎn)子支承系統(tǒng)的近似特征頻率Λ和特征向量Φ。
根據(jù)式（7）的階次展開形式，可知每階類模態(tài)Φp對(duì)應(yīng)（2j＋1）個(gè)模態(tài)頻率，考慮到雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)時(shí)變剛度的變化圓頻率是2Ω，類模態(tài)Φp對(duì)應(yīng)的模態(tài)頻率中各階次頻率與基礎(chǔ)頻率的關(guān)系為
式中，ωj為第j階階次頻率，ω0為基礎(chǔ)頻率，等于j＝0時(shí)的階次頻率。
根據(jù)式（12）對(duì)計(jì)算得到的特征頻率按照基礎(chǔ)頻率分簇，并計(jì)算得到該基礎(chǔ)頻率簇對(duì)應(yīng)的類模態(tài)。類模態(tài)各階次頻率的實(shí)部作為判斷該階類模態(tài)穩(wěn)定性的依據(jù)。
4 算 例
這里用兩個(gè)算例對(duì)文中分析方法的有效性及工程實(shí)際應(yīng)用進(jìn)行說明。
算例一中的分析模型來源于一篇參考文獻(xiàn)［22］，文獻(xiàn)作者基于1D有限元方法對(duì)模型的振動(dòng)特性進(jìn)行了分析。通過將本文基于3D有限元方法得到的結(jié)果與參考文獻(xiàn)結(jié)果進(jìn)行對(duì)比，對(duì)文中分析方法的有效性進(jìn)行驗(yàn)證。
算例二的分析模型是某型航空發(fā)動(dòng)機(jī)的低壓轉(zhuǎn)子。該轉(zhuǎn)子由三級(jí)風(fēng)扇、兩級(jí)渦輪及連接軸組成，由三個(gè)支點(diǎn)進(jìn)行支撐，支點(diǎn)剛度包含機(jī)匣剛度和軸承剛度。分析模型中考慮了中間支點(diǎn)機(jī)匣支承剛度在水平和豎直方向上的差異以及轉(zhuǎn)子連接軸因裝配誤差引起的彎曲剛度差異?；谖闹械姆治龇椒ǎ瑢?duì)上述非對(duì)稱特征對(duì)轉(zhuǎn)子振動(dòng)特性的影響進(jìn)行分析。
4.1 算例一
根據(jù)參考文獻(xiàn)中提供的結(jié)構(gòu)數(shù)據(jù)，建立的轉(zhuǎn)子分析模型如圖2所示。
圖2 非對(duì)稱轉(zhuǎn)子的有限元模型Fig．2 Finite element model of asymmetric rotor
實(shí)體單元模型在支承結(jié)構(gòu)處存在多個(gè)節(jié)點(diǎn)，不便于施加支承剛度約束，通常的做法是將支承位置的整圈實(shí)體單元節(jié)點(diǎn)剛化到圓心節(jié)點(diǎn)上，對(duì)圓心節(jié)點(diǎn)施加相應(yīng)的約束［23］，文中沿用這種“剛化”處理方法。
利用商用有限元軟件ANSYS的矩陣生成及組集功能，可以得到轉(zhuǎn)子結(jié)構(gòu)模型的質(zhì)量、剛度和科氏力矩陣以及反映節(jié)點(diǎn)編號(hào)與矩陣行號(hào)對(duì)應(yīng)關(guān)系的mapping文件。利用這些矩陣及對(duì)應(yīng)關(guān)系，在MATLAB中對(duì)支承節(jié)點(diǎn)所對(duì)應(yīng)的矩陣位置添加時(shí)變支承剛度，形成整個(gè)系統(tǒng)的時(shí)變運(yùn)動(dòng)微分方程。需要注意的是，在ANSYS中生成矩陣時(shí)，需要施加除徑向支承剛度之外的其它約束條件，如軸向約束，轉(zhuǎn)速等。
4．1．1 剛性支承下的結(jié)構(gòu)驗(yàn)證
是基于1D框架給出的轉(zhuǎn)子模型描述數(shù)據(jù)，如盤的極轉(zhuǎn)動(dòng)慣量、直徑轉(zhuǎn)動(dòng)慣量等，而沒有轉(zhuǎn)子結(jié)構(gòu)的3D幾何信息。因此需對(duì)本文模型的結(jié)構(gòu)特征進(jìn)行驗(yàn)證。在轉(zhuǎn)子兩端支點(diǎn)處施加剛化約束，計(jì)算零轉(zhuǎn)速下該非對(duì)稱轉(zhuǎn)子的模態(tài)頻率，與文獻(xiàn)算例模型在相同邊界條件下的計(jì)算結(jié)果比較，可驗(yàn)證轉(zhuǎn)子結(jié)構(gòu)模型的正確性。本文模型與文獻(xiàn)模型的零轉(zhuǎn)速一階彎曲模態(tài)頻率列于表1中，兩者誤差較小，說明本文模型與文獻(xiàn)模型具有一致的結(jié)構(gòu)特征。
表1 本文與文獻(xiàn)模型在零轉(zhuǎn)速下的一階彎曲頻率Tab.1 First bending frequencies of the model for zero speed
4．1．2 彈性支承下的計(jì)算結(jié)果
在模型結(jié)構(gòu)一致的基礎(chǔ)上，考慮非對(duì)稱支承剛度對(duì)轉(zhuǎn)子振動(dòng)特性的影響，分析該轉(zhuǎn)子在旋轉(zhuǎn)坐標(biāo)系下的固有頻率隨轉(zhuǎn)速變化的情況。
圖3和圖4分別給出了解的傅里葉展開指數(shù)為p＝0和p＝1時(shí)，轉(zhuǎn)子在旋轉(zhuǎn)坐標(biāo)系下的頻率轉(zhuǎn)速變化曲線。
圖3中，截?cái)嚯A次p＝0，系統(tǒng)運(yùn)動(dòng)微分方程的等效線性方程組（9）僅包含A0項(xiàng)，不包含因支承剛度非對(duì)稱引起的時(shí)變系數(shù)Δ項(xiàng)，此時(shí)求解得到的結(jié)果是轉(zhuǎn)子在平均支承剛度下的基礎(chǔ)振動(dòng)頻率w0。a0表示的是轉(zhuǎn)子的一階彎曲基礎(chǔ)頻率的正進(jìn)動(dòng)曲線。在轉(zhuǎn)速范圍820～950 r／min內(nèi)，a0頻率等于零，相應(yīng)的特征根實(shí)部大于零，與參考文獻(xiàn)中給出的轉(zhuǎn)速不穩(wěn)定區(qū)范圍829～976 r／min一致。
圖4中，截?cái)嚯A次p＝1，計(jì)算得到的特征頻率除包含基礎(chǔ)振動(dòng)頻率w0外，還包含第一階階次振動(dòng)頻率w－1和w＋1。理論上，階次頻率與基礎(chǔ)頻率的相互關(guān)系由式（12）給出。圖4中，a0，a＋1和a－1分別是實(shí)際計(jì)算得到的基礎(chǔ)振動(dòng)頻率和＋1、－1階階次振動(dòng)頻率隨轉(zhuǎn)速的變化曲線。從圖4中可以看出，a0曲線始終位于a＋1和a－1中間，并與它們的頻率相差2Ω，與式（12）中給出的關(guān)系完全一致，并且，階次頻率變化曲線a＋1和a－1隨著轉(zhuǎn)速的變化發(fā)生耦合，產(chǎn)生新的轉(zhuǎn)速不穩(wěn)定區(qū)。上述結(jié)論與參考文獻(xiàn)一致，通過對(duì)簡(jiǎn)單轉(zhuǎn)子振動(dòng)特性的分析，文中分析方法的有效性得到驗(yàn)證。
圖3 p＝0時(shí)旋轉(zhuǎn)坐標(biāo)系下的頻率轉(zhuǎn)速圖Fig．3 Evolution of frequencies in the rotating coordinate system for p＝0
圖4 p＝1時(shí)旋轉(zhuǎn)坐標(biāo)系下的頻率轉(zhuǎn)速圖Fig．4 Evolution of frequencies in the rotating coordinate system for p＝1
4.2 算例二
圖5所示為某型航空發(fā)動(dòng)機(jī)低壓轉(zhuǎn)子的三維有限元模型。模型通過改變連接軸的部分單元材料參數(shù)模擬轉(zhuǎn)子在裝配過程中可能產(chǎn)生的在正交方向上的彎曲剛度差異，材料參數(shù)改變區(qū)域在圖5的局部放大圖中以紫色標(biāo)示。
圖5 某型航空發(fā)動(dòng)機(jī)轉(zhuǎn)子有限元模型Fig．5 Finite element rotor model of an aero-engine
該轉(zhuǎn)子采用1－1－1型支承方案，各支點(diǎn)位置如圖5中標(biāo)號(hào)1、2、3所示，其中2號(hào)支點(diǎn)在水平和豎直方向上具有不同的支承剛度，各支點(diǎn)的支承剛度列于表2中。
表2 某型發(fā)動(dòng)機(jī)低壓轉(zhuǎn)子各支點(diǎn)支承剛度Tab.2 Bearing stiffnesses of the low pressure rotor of the aero-engine
采用文中的雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)振動(dòng)特性分析方法，對(duì)該低壓轉(zhuǎn)子實(shí)體有限元模型在旋轉(zhuǎn)坐標(biāo)系下進(jìn)行模態(tài)特性分析，得到轉(zhuǎn)子模態(tài)頻率隨轉(zhuǎn)速的變化曲線。
圖6和圖7分別是考慮方程解的傅里葉展開階次p＝0和p＝1時(shí)的模態(tài)頻率隨轉(zhuǎn)速的變化曲線。
圖6中，截?cái)嚯A次p＝0，得到非對(duì)稱轉(zhuǎn)子在對(duì)稱支承剛度下的頻率轉(zhuǎn)速變化曲線。當(dāng)各階振型的正進(jìn)動(dòng)頻率（圖中a0，c0和e0）為0時(shí)，轉(zhuǎn)子進(jìn)入不穩(wěn)定區(qū)。需要說明的是，在轉(zhuǎn)速為0時(shí)，轉(zhuǎn)子僅在第四階振型上表現(xiàn)出明顯的頻率分離，前三階振型在該轉(zhuǎn)速下對(duì)轉(zhuǎn)子連接軸段的局部彎曲剛度變化不敏感。
圖6 p＝0時(shí)低壓轉(zhuǎn)子的頻率轉(zhuǎn)速變化曲線Fig．6 Evolution of nature frequencies of the low pressure rotor for p＝0
圖7 p＝1時(shí)低壓轉(zhuǎn)子的頻率轉(zhuǎn)速變化曲線Fig．7 Evolution of nature frequencies of the low pressure rotor for p＝1
圖7中，截?cái)嚯A次p＝1，除得到轉(zhuǎn)子的基礎(chǔ)頻率外，還得到轉(zhuǎn)子振動(dòng)類模態(tài)的一階階次頻率。從圖7中可以看出，考慮轉(zhuǎn)子的一階階次頻率后，轉(zhuǎn)子的頻率轉(zhuǎn)速變換曲線變得更加復(fù)雜。圖6中的每條基礎(chǔ)頻率w0，n曲線均擴(kuò)展為w－1，n，w0，n，w＋1，n三條曲線，并且不是簡(jiǎn)單的擴(kuò)展，從圖7中可以看到，階次頻率與基礎(chǔ)頻率以及階次頻率之間出現(xiàn)多處的耦合和頻率轉(zhuǎn)向現(xiàn)象。階次頻率的耦合增加了轉(zhuǎn)子的不穩(wěn)定區(qū)范圍。因此，與僅考慮單一非對(duì)稱的轉(zhuǎn)子不同，考慮轉(zhuǎn)子的雙重非對(duì)稱特征時(shí)，轉(zhuǎn)子的不穩(wěn)定區(qū)區(qū)間增多，轉(zhuǎn)子的安全工作轉(zhuǎn)速范圍變小。
另外需要說明的是，由于頻率轉(zhuǎn)向與階次頻率耦合，各階頻率與振型的對(duì)應(yīng)關(guān)系非常復(fù)雜，并且分析人員更加關(guān)注的是不穩(wěn)定區(qū)的范圍，進(jìn)行詳細(xì)的振型區(qū)分也沒有重要的意義，因此在圖7中用相同的符號(hào)標(biāo)示各階頻率變化曲線。
5 結(jié) 論
本文提出了基于3D實(shí)體有限單元建模的雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的振動(dòng)特性分析方法，詳細(xì)給出了系統(tǒng)時(shí)變運(yùn)動(dòng)微分方程的組集推導(dǎo)過程以及頻域特性的分析求解方法。文中以一個(gè)具有簡(jiǎn)單結(jié)構(gòu)的雙重非對(duì)稱轉(zhuǎn)子支承模型為例說明了方法的分析流程，對(duì)ANSYS模型分析矩陣的輸出和MATLAB中支承剛度及約束施加的方法進(jìn)行了必要的說明。對(duì)某型航空發(fā)動(dòng)機(jī)低壓轉(zhuǎn)子的振動(dòng)特性分析進(jìn)一步說明的文中所提方法的工程實(shí)用價(jià)值。
參考文獻(xiàn)
［1］Lee C W．Vibration analysis of rotors［M］．Dordrecht：Kluwer Academic Publishers，1993．
［2］Combescure D，Lazarus A．Refined finite element modelling forthevibrationanalysisoflargerotatingmachines：application to the gas turbine modular helium reactor power conversion unit［J］．Journal of Sound and Vibration，2008，318（4－5）：1262－1280．
［3］陳萌，洪杰，朱彬，等．基于實(shí)體單元的轉(zhuǎn)子動(dòng)力特性計(jì)算方法［J］．北京航空航天大學(xué)學(xué)報(bào)，2007，33（1）：10－13．
CHEN Meng，HONG Jie，ZHU Bin，et al．Rotor dynamic analysis based on solid element［J］．Journal of Beijing University of Aeronautics and Astronautics，2007，33（1）：10－13．
［4］Chatelet E，D'Ambrosio F，Jacquet-Richardet G．Toward global modelling approaches for dynamic analyses of rotating assemblies of turbomachines［J］．Journal of Sound and Vibration，2005，282（1－2）：163－178．
［5］Genta G．Dynamics of rotating systems［M］．springer，2005．
［6］Black H F．Parametrically excited lateral vibrations of an asymmetric slender shaft in asymmetrically flexible bearings ［J］．Journal of Mechanical Engineering Science，1969，11 （1）：57－67．
［7］Yamamoto T，Ota H，Kono K．On the unstable vibrations of a shaft with unsymmetrical stiffness carrying an unsymmetrical rotor［J］．Journal of Applied Mechanics，1968，35（2）：313 －321．
［8］Gladwell G M L，Stammers C W．On the stability of an unsymmetrical rigid rotor supported in unsymmetrical bearings ［J］．Journal of Sound and Vibration，1966，3：221－232．
［9］Genta G．Whirling of unsymmetrical rotors：a finite element approach based on complex coordinates［J］．Journal of Sound and Vibration，1988，124（1）：27－53．
［10］Lee C W，Suh J H，Hong S W．Modal analysis of rotor system using modulated co-ordinates for asymmetric rotor systemwithanisotropicstator［C］．8thInternational Conference on Vibrations in Rotating Machinery．University ofSwansea，UK：ProfessionalEngineeringPublishing （Institution of Mechanical Engineers），2004：575－584．
［11］Lee C W，Han D J，Hong S W．Modal analysis of periodically time-varying linear rotor systems using floquet theory［C］．Proceedings of the 7th IFToMM-Conference on Rotor Dynamics．Austria：International Conference on Rotor Dynamics，2006．
［12］Nandi A，Neogy S．An efficient scheme for stability analysis of finite element asymmetric rotor models in a rotating frame ［J］．Finite Elements in Analysis and Design，2005，41 （14）：1343－1364．
［13］Nandi A．Reduction of finite element equations for a rotor model on non-isotropic spring support in a rotating frame［J］．Finite Elements in Analysis and Design，2004，40（9）：935 －952．［14］Lazarus A，Prabel B，Combescure D．A 3d finite element model forthevibrationanalysisofasymmetricrotating machines［J］．Journal of Sound and Vibration，2010，329 （18）：3780－3797．
［15］Kim K T，Lee C W．Dynamic analysis of asymmetric bladedrotors supported by anisotropic stator［J］．Journal of Sound and Vibration，2012，331（24）：5224－5246．
［16］沈松，鄭兆昌，應(yīng)懷樵．非對(duì)稱轉(zhuǎn)子－軸承－基礎(chǔ)系統(tǒng)的非線性振動(dòng)［J］．振動(dòng)與沖擊，2004，23（4）：33－35．
SHENSong，ZHENGZhao-chang，YINGHuai-qiao．Nonlinear vibration of unsymmetrical rotor-bearing-foundation system［J］．Journal of Vibration and Shock，2004，23（4）：33－35．
［17］肖錫武，肖光華，楊叔子．不對(duì)稱轉(zhuǎn)子系統(tǒng)的參激強(qiáng)迫振動(dòng)［J］．振動(dòng)工程學(xué)報(bào)，2002（3）：71－74．
XIAO Xi-wu，XIAO Guang-hua，YANG Shu-zi．Parametrically excited and forced oscillations of an unsymmetrical rotor system［J］．Journal of Vibration Engineering，2002（3）：71－74．
［18］肖錫武，楊正茂，肖光華，等．不對(duì)稱轉(zhuǎn)子系統(tǒng)的非線性振動(dòng)［J］．華中科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2002（5）：81 －84．
XIAO Xi-wu，YANG Zheng-mao，XIAO Guang-hua，et al． Nonlinear oscillations of an unsymmetrical rotor system［J］．J．Huazhong Univ．of Sci．＆Tech．（Nature Science Edition），2002（5）：81－84．
［19］王培杰，滕弘飛，劉篤金，等．非對(duì)稱復(fù)雜轉(zhuǎn)子系統(tǒng)有限元?jiǎng)恿Ψ治觯跩］．固體力學(xué)學(xué)報(bào)，1991（4）：370－376．
WANG Pei-jie，TENG Hong-fei，LIU Du-jin，et al．The dynamic analysis for complex asymmetric rotor-bearing systems by finite element method［J］．Acta Mechanica Solida Sinica，1991（4）：370－376．
［20］Lee C W，Han D J，Suh J H，et al．Modal analysis of periodically time-varying linear rotor systems［J］．Journal of Sound and Vibration，2007，303（3－5）：553－574．
［21］Han D J．Modal analysis of periodically time-varying rotor systems using floquet theory and modulated coordinate transformation：applications to crack detection and modal balancing［D］．Daejeon：Korea Advanced Institute of Science and Technology，2006．
［22］Boru F E．Numerical and experimental response and stability investigations of anisotropic rotor-bearing systems［D］．Kassel：kassel university，2010．
［23］Rao J S，Sreenivas R．Dynamics of a three level rotor system using solid elements［C］．Proceedings of ASME Turbo Expo 2003．Atlanta，Georgia，USA：ASME，2003，601－606．
Dynamic analysis method of an anisotropic and asymmetric rotor-bearing system based on 3D FEM
MA Wei-meng，WANG Jian-jun
（School of Energy and Power Engineering，Beijing University of Aeronautics and Astronautics，Beijing 100191，China）
Dynamic analysis method of an anisotropic and asymmetric rotor-bearing system based on 3D finite element method was presented．The system mass，stiffness and Coriolis force matrices expressed with respect to a rotating coordinate system were generated with ANSYS．Time-varying bearing stiffness matrix with respect to the rotating coordinates was then obtained with a coordinate transformation．The differential equations of motion of the system were generated by applying constraints and time-varying bearing stiffness to the system matrices using the corresponding relationships between node number and row index of the matrices．Floquet theory and Hill infinite determinant method were adopted to solve the differential equations，frequency characteristics and stability of the system were obtained．The proposed method was validated based on a model in literature，and the proposed method was applied to analyze the dynamic behavior of an industrial rotor considering anisotropic and asymmetric characters．
asymmetric；rotor-bearing system；3D finite element；time-varying system；stability
V 231．96
A
10．13465／j．cnki．jvs．2014．23．002
2013－07－16 修改稿收到日期：2014－01－02
馬威猛男，博士生，1987年1月生
王建軍男，教授，博士生導(dǎo)師，1956年生
猜你喜歡
階次分析在驅(qū)動(dòng)橋異響中的應(yīng)用汽車實(shí)用技術(shù)(2019年24期)2019-12-27 03:52:46 基于Vold-Kalman濾波的階次分析系統(tǒng)設(shè)計(jì)與實(shí)現(xiàn)*組合機(jī)床與自動(dòng)化加工技術(shù)(2019年7期)2019-08-06 03:51:06 非對(duì)稱Orlicz差體數(shù)學(xué)年刊A輯(中文版)(2018年2期)2019-01-08 01:59:50 基于齒輪階次密度優(yōu)化的變速器降噪研究價(jià)值工程(2017年28期)2018-01-23 20:48:29 基于時(shí)變Copula的股票市場(chǎng)相關(guān)性分析智富時(shí)代(2017年4期)2017-04-27 17:08:47 煙氣輪機(jī)復(fù)合故障時(shí)變退化特征提取廣東石油化工學(xué)院學(xué)報(bào)(2016年6期)2016-05-17 05:17:43 點(diǎn)數(shù)不超過20的旗傳遞非對(duì)稱2-設(shè)計(jì)數(shù)學(xué)理論與應(yīng)用(2016年4期)2016-05-17 04:50:23 基于MEP法的在役橋梁時(shí)變可靠度研究中國鐵道科學(xué)(2015年4期)2015-06-21 06:46:08 非對(duì)稱負(fù)載下矩陣變換器改進(jìn)型PI重復(fù)控制電測(cè)與儀表(2015年4期)2015-04-12 00:43:04 自適應(yīng)兩級(jí)UKF算法及其在時(shí)變偏差估計(jì)中的應(yīng)用空間控制技術(shù)與應(yīng)用(2010年3期)2010-12-23 08:04:54
馬威猛，王建軍
（北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京 100191）
雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)實(shí)體有限元?jiǎng)恿Ψ治?/p>
馬威猛，王建軍
（北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京 100191）
研究了同時(shí)考慮轉(zhuǎn)子和支承非對(duì)稱特征的轉(zhuǎn)子支承系統(tǒng)3D有限元建模和振動(dòng)特性分析方法。首先用ANSYS得到轉(zhuǎn)子在旋轉(zhuǎn)坐標(biāo)系下的整體質(zhì)量、剛度和科氏力矩陣；其次利用動(dòng)靜坐標(biāo)變換關(guān)系得到旋轉(zhuǎn)坐標(biāo)系下的時(shí)變支承剛度；最后利用節(jié)點(diǎn)編號(hào)與矩陣行號(hào)的對(duì)應(yīng)關(guān)系在整體矩陣中施加約束和時(shí)變支承剛度，得到系統(tǒng)的時(shí)變系數(shù)運(yùn)動(dòng)微分方程?；诟ヂ鍎P理論和赫爾無窮行列式方法對(duì)方程求解得到系統(tǒng)的頻率特性和穩(wěn)定性。利用文獻(xiàn)算例模型對(duì)該方法的有效性進(jìn)行了驗(yàn)證，并對(duì)實(shí)際工程轉(zhuǎn)子在考慮雙重非對(duì)稱特征時(shí)的振動(dòng)特性進(jìn)行了分析。
非對(duì)稱；轉(zhuǎn)子支承系統(tǒng)；3D有限元；時(shí)變系統(tǒng)；穩(wěn)定性
轉(zhuǎn)子支承系統(tǒng)是旋轉(zhuǎn)機(jī)械的核心部件，其振動(dòng)特性對(duì)旋轉(zhuǎn)機(jī)械的工作性能和可靠性具有重要的影響。轉(zhuǎn)子支承系統(tǒng)振動(dòng)特性研究對(duì)旋轉(zhuǎn)機(jī)械的結(jié)構(gòu)設(shè)計(jì)、維修保養(yǎng)、振動(dòng)控制、故障診斷等方面工作都具有重要的意義。
根據(jù)轉(zhuǎn)子與支承的軸對(duì)稱特征，轉(zhuǎn)子支承系統(tǒng)可以分為：軸對(duì)稱（axisymmetric）、支承非對(duì)稱（anisotropic）、轉(zhuǎn)子非對(duì)稱（asymmetric）和雙重非對(duì)稱（anisotropic asymmetric）轉(zhuǎn)子支承系統(tǒng)，其中雙重非對(duì)稱系統(tǒng)又稱為一般（general）轉(zhuǎn)子支承系統(tǒng)［1］。人們對(duì)軸對(duì)稱轉(zhuǎn)子支承系統(tǒng)開展了充分的研究工作，基于3D實(shí)體有限元的分析模型可以精確描述復(fù)雜截面轉(zhuǎn)子的振動(dòng)特性［2－4］。而在工程實(shí)際中，因加工、設(shè)計(jì)或故障破壞等原因，轉(zhuǎn)子或支承部件通常會(huì)包含一些非對(duì)稱特征，典型的如發(fā)電機(jī)轉(zhuǎn)子、裂紋轉(zhuǎn)子等。轉(zhuǎn)靜子上存在的非對(duì)稱特征使得轉(zhuǎn)子表現(xiàn)出獨(dú)特的振動(dòng)特性，需要進(jìn)行專門的研究，如參數(shù)振動(dòng)特性是雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的主要特征［5］。
20世紀(jì)60年代以來，國外學(xué)者針對(duì)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的振動(dòng)特性開展了廣泛的定性［6－8］、定量研究［9－12］。在建模分析方法方面，Nandi［12－14］在不同的框架下提出并開展了基于3D實(shí)體有限元的非對(duì)稱轉(zhuǎn)子系統(tǒng)動(dòng)力特性分析工作，Kim［15］仍在1D有限元的框架下考慮非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的部件耦合問題。沈松等［16－19］在非對(duì)稱轉(zhuǎn)子系統(tǒng)的振動(dòng)分析方面也做了出色的工作，但主要集中在非線性分析領(lǐng)域，對(duì)線性領(lǐng)域內(nèi)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的建模和振動(dòng)特性分析方法沒有加以關(guān)注。
3D實(shí)體單元建模不需對(duì)轉(zhuǎn)靜子結(jié)構(gòu)進(jìn)行過多簡(jiǎn)化，可以真實(shí)反映轉(zhuǎn)靜子的結(jié)構(gòu)特征，并且能夠準(zhǔn)確考慮轉(zhuǎn)靜子各部件之間的耦合效應(yīng)，是轉(zhuǎn)子動(dòng)力學(xué)建模方法的發(fā)展趨勢(shì)。目前，多數(shù)商用有限元分析軟件已支持基于3D實(shí)體單元建模的轉(zhuǎn)子動(dòng)力學(xué)分析，但其分析功能僅能考慮軸對(duì)稱或單一非對(duì)稱轉(zhuǎn)子支承系統(tǒng)，而不能進(jìn)行雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的振動(dòng)特性計(jì)算。
本文利用商用有限元軟件ANSYS的前處理功能，得到轉(zhuǎn)子3D有限元模型的質(zhì)量、剛度和科氏力矩陣，通過矩陣處理加入支承結(jié)構(gòu)在旋轉(zhuǎn)坐標(biāo)描述下的時(shí)變剛度，得到轉(zhuǎn)子支承系統(tǒng)的時(shí)變運(yùn)動(dòng)微分方程，然后基于時(shí)變系統(tǒng)分析的弗洛凱理論和無窮行列式方法將時(shí)變方程轉(zhuǎn)化為一定截?cái)嚯A次的時(shí)不變線性方程組，求解得到系統(tǒng)的頻率特性和穩(wěn)定性，實(shí)現(xiàn)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的3D有限元振動(dòng)特性分析功能。文中首先介紹相關(guān)的分析理論，包括坐標(biāo)系的選擇，運(yùn)動(dòng)微分方程的建立以及基于弗洛凱理論的時(shí)變參數(shù)系統(tǒng)的求解方法；然后以兩個(gè)算例對(duì)文中雙重非對(duì)稱轉(zhuǎn)子振動(dòng)特性分析方法的有效性和工程應(yīng)用價(jià)值進(jìn)行驗(yàn)證和說明。
1 坐標(biāo)系選擇
工程實(shí)際中，裂紋、鍵槽、非圓截面等非對(duì)稱轉(zhuǎn)子特征及各種形式的支承非對(duì)稱特征可以抽象為圖1所示的典型雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)［20］。圖1中的非圓截面軸代表前述多種因素引起的轉(zhuǎn)子剛度非對(duì)稱特征，非圓截面盤代表轉(zhuǎn)子的非對(duì)稱轉(zhuǎn)動(dòng)慣量特征，正交方向上不同的支承剛度和阻尼則表示支承系統(tǒng)的非對(duì)稱特征。
圖1 雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)示意圖Fig．1 Simple analysis model of anisotropic asymmetric rotor bearing system
對(duì)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)而言，采用旋轉(zhuǎn)或者固定坐標(biāo)系描述系統(tǒng)運(yùn)動(dòng)不像單一非對(duì)稱轉(zhuǎn)子支承系統(tǒng)那樣具有特殊的意義－選用合適的坐標(biāo)系可以把時(shí)變微分方程轉(zhuǎn)化為時(shí)不變微分方程的求解問題。雙重非對(duì)稱轉(zhuǎn)子系統(tǒng)在固定或轉(zhuǎn)動(dòng)坐標(biāo)系下描述均不能得到具有固定系數(shù)的運(yùn)動(dòng)微分方程。雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)在選取分析坐標(biāo)系時(shí)，考慮到：①目前ANSYS等商用有限元分析軟件僅支持非對(duì)稱轉(zhuǎn)子在旋轉(zhuǎn)坐標(biāo)系下的3D實(shí)體有限元分析；②若將靜子系統(tǒng)的支承剛度以參數(shù)形式給出，在旋轉(zhuǎn)坐標(biāo)系下表達(dá)的時(shí)變剛度系數(shù)僅涉及支承點(diǎn)附近的幾個(gè)相關(guān)自由度，表達(dá)形式簡(jiǎn)單。因此，一般選用隨轉(zhuǎn)子同速旋轉(zhuǎn)的旋轉(zhuǎn)坐標(biāo)系對(duì)雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的運(yùn)動(dòng)進(jìn)行描述。
2 旋轉(zhuǎn)坐標(biāo)系下的運(yùn)動(dòng)微分方程
在旋轉(zhuǎn)坐標(biāo)系下，采用3D實(shí)體單元對(duì)非對(duì)稱轉(zhuǎn)子進(jìn)行建模能夠準(zhǔn)確反映非對(duì)稱幾何特征或材料特征對(duì)轉(zhuǎn)子截面彎曲剛度或轉(zhuǎn)動(dòng)慣量的影響。實(shí)體非對(duì)稱轉(zhuǎn)子在邊界節(jié)點(diǎn)附近由剛度單元支承，剛度單元在各個(gè)支承方向的剛度差異反映支承系統(tǒng)的非對(duì)稱特征。為方便描述，以下標(biāo)nb標(biāo)示轉(zhuǎn)子的支承節(jié)點(diǎn)。在固定坐標(biāo)系下，這些支點(diǎn)的剛度可以表示為［13］
轉(zhuǎn)子支點(diǎn)在旋轉(zhuǎn)坐標(biāo)系下的支承剛度可由式（1）經(jīng)坐標(biāo)變換得到
旋轉(zhuǎn)坐標(biāo)系下，無阻尼轉(zhuǎn)子支承系統(tǒng)在外力作用下的運(yùn)動(dòng)微分方程可以寫為
式中［Kr］是轉(zhuǎn)子支承系統(tǒng)的剛度矩陣，由轉(zhuǎn)子自身剛度和支承剛度組成，［Mr］是轉(zhuǎn)子系統(tǒng)在旋轉(zhuǎn)坐標(biāo)系下的質(zhì)量矩陣，［Cr］被稱為科氏力矩陣，是反對(duì)稱矩陣，｛U｝是轉(zhuǎn)子系統(tǒng)在旋轉(zhuǎn)坐標(biāo)系下的位移向量，｛f｝則是在旋轉(zhuǎn)坐標(biāo)系下表示的外部載荷向量。
在ANSYS前處理模塊中完成轉(zhuǎn)子結(jié)構(gòu)建模后，若不考慮支承的非對(duì)稱特征，利用ANSYS的轉(zhuǎn)子動(dòng)力學(xué)分析模塊可以直接得到方程（3）中各矩陣的元素?cái)?shù)值；而考慮支承非對(duì)稱時(shí)，剛度矩陣是時(shí)變的，需要進(jìn)一步編程處理才能得到相應(yīng)的矩陣。
考慮到支承剛度在旋轉(zhuǎn)坐標(biāo)系下的時(shí)變特征，方程（3）可以改寫為
式中［Kr］nb是不考慮支承剛度時(shí)，轉(zhuǎn)子系統(tǒng)在旋轉(zhuǎn)坐標(biāo)系下的剛度矩陣，可由ANSYS前處理模塊直接得到；［K0］b是支承剛度的非時(shí)變項(xiàng)，而［Kc］b和［Ks］b則是支承剛度時(shí)變項(xiàng)中的余弦和正弦分量，支承剛度各矩陣由式（2）經(jīng)坐標(biāo)變換得到。
考慮到三角函數(shù)與復(fù)數(shù)的變換關(guān)系，方程（4）可寫成復(fù)數(shù)形式，
3 時(shí)變參數(shù)系統(tǒng)的求解方法
時(shí)變參數(shù)系統(tǒng)的求解方法與一般時(shí)不變系統(tǒng)不同。弗洛凱理論和赫爾無窮行列式方法是分析時(shí)變參數(shù)系統(tǒng)常用的分析理論。
基于方程（5）的齊次方程，可對(duì)轉(zhuǎn)子支承系統(tǒng)的自由振動(dòng)特性進(jìn)行分析。根據(jù)弗洛凱理論，假設(shè)方程解的形式為
其中，Λ是系統(tǒng)的基礎(chǔ)特征頻率，Φ是與Λ相對(duì)應(yīng)的類模態(tài)矩陣。之所以稱Φ為類模態(tài)，是因?yàn)槠湫问诫m然與時(shí)不變系統(tǒng)的模態(tài)矩陣相似，但Φ具有時(shí)變的性質(zhì)，其時(shí)變周期與運(yùn)動(dòng)微分方程中時(shí)變系統(tǒng)的周期一致。
其中，lj為Φ的第j階階次展開模態(tài)。
將式（6）和式（7）代入方程（5）的齊次方程中，并整理可以得到
方程（8）中，令指數(shù)項(xiàng)對(duì)應(yīng)各階次系數(shù)等于零可以得到
式（9）在截?cái)嚯A次j＝∞時(shí)，與原時(shí)變微分方程是等價(jià)的，但在實(shí)際的計(jì)算中，需要進(jìn)行適當(dāng)階次的截?cái)?，得到原時(shí)變微分方程解的近似值。求解得到各展開階次的頻率及特征向量后，根據(jù)式（7）可得到原時(shí)變微分運(yùn)動(dòng)方程的近似特征值Λ和類模態(tài)Φ。
式（9）可寫成特征值求解問題的一般形式
上述特征值問題一般轉(zhuǎn)換到狀態(tài)空間進(jìn)行求解。式（10）在狀態(tài)空間中可以寫為
在狀態(tài)空間中求解式（11）的特征值問題，可得到相應(yīng)的特征頻率λ和特征向量Ψ，進(jìn)而得到轉(zhuǎn)子支承系統(tǒng)的近似特征頻率Λ和特征向量Φ。
根據(jù)式（7）的階次展開形式，可知每階類模態(tài)Φp對(duì)應(yīng)（2j＋1）個(gè)模態(tài)頻率，考慮到雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)時(shí)變剛度的變化圓頻率是2Ω，類模態(tài)Φp對(duì)應(yīng)的模態(tài)頻率中各階次頻率與基礎(chǔ)頻率的關(guān)系為
式中，ωj為第j階階次頻率，ω0為基礎(chǔ)頻率，等于j＝0時(shí)的階次頻率。
根據(jù)式（12）對(duì)計(jì)算得到的特征頻率按照基礎(chǔ)頻率分簇，并計(jì)算得到該基礎(chǔ)頻率簇對(duì)應(yīng)的類模態(tài)。類模態(tài)各階次頻率的實(shí)部作為判斷該階類模態(tài)穩(wěn)定性的依據(jù)。
4 算 例
這里用兩個(gè)算例對(duì)文中分析方法的有效性及工程實(shí)際應(yīng)用進(jìn)行說明。
算例一中的分析模型來源于一篇參考文獻(xiàn)［22］，文獻(xiàn)作者基于1D有限元方法對(duì)模型的振動(dòng)特性進(jìn)行了分析。通過將本文基于3D有限元方法得到的結(jié)果與參考文獻(xiàn)結(jié)果進(jìn)行對(duì)比，對(duì)文中分析方法的有效性進(jìn)行驗(yàn)證。
算例二的分析模型是某型航空發(fā)動(dòng)機(jī)的低壓轉(zhuǎn)子。該轉(zhuǎn)子由三級(jí)風(fēng)扇、兩級(jí)渦輪及連接軸組成，由三個(gè)支點(diǎn)進(jìn)行支撐，支點(diǎn)剛度包含機(jī)匣剛度和軸承剛度。分析模型中考慮了中間支點(diǎn)機(jī)匣支承剛度在水平和豎直方向上的差異以及轉(zhuǎn)子連接軸因裝配誤差引起的彎曲剛度差異?；谖闹械姆治龇椒ǎ瑢?duì)上述非對(duì)稱特征對(duì)轉(zhuǎn)子振動(dòng)特性的影響進(jìn)行分析。
4.1 算例一
根據(jù)參考文獻(xiàn)中提供的結(jié)構(gòu)數(shù)據(jù)，建立的轉(zhuǎn)子分析模型如圖2所示。
圖2 非對(duì)稱轉(zhuǎn)子的有限元模型Fig．2 Finite element model of asymmetric rotor
實(shí)體單元模型在支承結(jié)構(gòu)處存在多個(gè)節(jié)點(diǎn)，不便于施加支承剛度約束，通常的做法是將支承位置的整圈實(shí)體單元節(jié)點(diǎn)剛化到圓心節(jié)點(diǎn)上，對(duì)圓心節(jié)點(diǎn)施加相應(yīng)的約束［23］，文中沿用這種“剛化”處理方法。
利用商用有限元軟件ANSYS的矩陣生成及組集功能，可以得到轉(zhuǎn)子結(jié)構(gòu)模型的質(zhì)量、剛度和科氏力矩陣以及反映節(jié)點(diǎn)編號(hào)與矩陣行號(hào)對(duì)應(yīng)關(guān)系的mapping文件。利用這些矩陣及對(duì)應(yīng)關(guān)系，在MATLAB中對(duì)支承節(jié)點(diǎn)所對(duì)應(yīng)的矩陣位置添加時(shí)變支承剛度，形成整個(gè)系統(tǒng)的時(shí)變運(yùn)動(dòng)微分方程。需要注意的是，在ANSYS中生成矩陣時(shí)，需要施加除徑向支承剛度之外的其它約束條件，如軸向約束，轉(zhuǎn)速等。
4．1．1 剛性支承下的結(jié)構(gòu)驗(yàn)證
是基于1D框架給出的轉(zhuǎn)子模型描述數(shù)據(jù)，如盤的極轉(zhuǎn)動(dòng)慣量、直徑轉(zhuǎn)動(dòng)慣量等，而沒有轉(zhuǎn)子結(jié)構(gòu)的3D幾何信息。因此需對(duì)本文模型的結(jié)構(gòu)特征進(jìn)行驗(yàn)證。在轉(zhuǎn)子兩端支點(diǎn)處施加剛化約束，計(jì)算零轉(zhuǎn)速下該非對(duì)稱轉(zhuǎn)子的模態(tài)頻率，與文獻(xiàn)算例模型在相同邊界條件下的計(jì)算結(jié)果比較，可驗(yàn)證轉(zhuǎn)子結(jié)構(gòu)模型的正確性。本文模型與文獻(xiàn)模型的零轉(zhuǎn)速一階彎曲模態(tài)頻率列于表1中，兩者誤差較小，說明本文模型與文獻(xiàn)模型具有一致的結(jié)構(gòu)特征。
表1 本文與文獻(xiàn)模型在零轉(zhuǎn)速下的一階彎曲頻率Tab.1 First bending frequencies of the model for zero speed
4．1．2 彈性支承下的計(jì)算結(jié)果
在模型結(jié)構(gòu)一致的基礎(chǔ)上，考慮非對(duì)稱支承剛度對(duì)轉(zhuǎn)子振動(dòng)特性的影響，分析該轉(zhuǎn)子在旋轉(zhuǎn)坐標(biāo)系下的固有頻率隨轉(zhuǎn)速變化的情況。
圖3和圖4分別給出了解的傅里葉展開指數(shù)為p＝0和p＝1時(shí)，轉(zhuǎn)子在旋轉(zhuǎn)坐標(biāo)系下的頻率轉(zhuǎn)速變化曲線。
圖3中，截?cái)嚯A次p＝0，系統(tǒng)運(yùn)動(dòng)微分方程的等效線性方程組（9）僅包含A0項(xiàng)，不包含因支承剛度非對(duì)稱引起的時(shí)變系數(shù)Δ項(xiàng)，此時(shí)求解得到的結(jié)果是轉(zhuǎn)子在平均支承剛度下的基礎(chǔ)振動(dòng)頻率w0。a0表示的是轉(zhuǎn)子的一階彎曲基礎(chǔ)頻率的正進(jìn)動(dòng)曲線。在轉(zhuǎn)速范圍820～950 r／min內(nèi)，a0頻率等于零，相應(yīng)的特征根實(shí)部大于零，與參考文獻(xiàn)中給出的轉(zhuǎn)速不穩(wěn)定區(qū)范圍829～976 r／min一致。
圖4中，截?cái)嚯A次p＝1，計(jì)算得到的特征頻率除包含基礎(chǔ)振動(dòng)頻率w0外，還包含第一階階次振動(dòng)頻率w－1和w＋1。理論上，階次頻率與基礎(chǔ)頻率的相互關(guān)系由式（12）給出。圖4中，a0，a＋1和a－1分別是實(shí)際計(jì)算得到的基礎(chǔ)振動(dòng)頻率和＋1、－1階階次振動(dòng)頻率隨轉(zhuǎn)速的變化曲線。從圖4中可以看出，a0曲線始終位于a＋1和a－1中間，并與它們的頻率相差2Ω，與式（12）中給出的關(guān)系完全一致，并且，階次頻率變化曲線a＋1和a－1隨著轉(zhuǎn)速的變化發(fā)生耦合，產(chǎn)生新的轉(zhuǎn)速不穩(wěn)定區(qū)。上述結(jié)論與參考文獻(xiàn)一致，通過對(duì)簡(jiǎn)單轉(zhuǎn)子振動(dòng)特性的分析，文中分析方法的有效性得到驗(yàn)證。
圖3 p＝0時(shí)旋轉(zhuǎn)坐標(biāo)系下的頻率轉(zhuǎn)速圖Fig．3 Evolution of frequencies in the rotating coordinate system for p＝0
圖4 p＝1時(shí)旋轉(zhuǎn)坐標(biāo)系下的頻率轉(zhuǎn)速圖Fig．4 Evolution of frequencies in the rotating coordinate system for p＝1
4.2 算例二
圖5所示為某型航空發(fā)動(dòng)機(jī)低壓轉(zhuǎn)子的三維有限元模型。模型通過改變連接軸的部分單元材料參數(shù)模擬轉(zhuǎn)子在裝配過程中可能產(chǎn)生的在正交方向上的彎曲剛度差異，材料參數(shù)改變區(qū)域在圖5的局部放大圖中以紫色標(biāo)示。
圖5 某型航空發(fā)動(dòng)機(jī)轉(zhuǎn)子有限元模型Fig．5 Finite element rotor model of an aero-engine
該轉(zhuǎn)子采用1－1－1型支承方案，各支點(diǎn)位置如圖5中標(biāo)號(hào)1、2、3所示，其中2號(hào)支點(diǎn)在水平和豎直方向上具有不同的支承剛度，各支點(diǎn)的支承剛度列于表2中。
表2 某型發(fā)動(dòng)機(jī)低壓轉(zhuǎn)子各支點(diǎn)支承剛度Tab.2 Bearing stiffnesses of the low pressure rotor of the aero-engine
采用文中的雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)振動(dòng)特性分析方法，對(duì)該低壓轉(zhuǎn)子實(shí)體有限元模型在旋轉(zhuǎn)坐標(biāo)系下進(jìn)行模態(tài)特性分析，得到轉(zhuǎn)子模態(tài)頻率隨轉(zhuǎn)速的變化曲線。
圖6和圖7分別是考慮方程解的傅里葉展開階次p＝0和p＝1時(shí)的模態(tài)頻率隨轉(zhuǎn)速的變化曲線。
圖6中，截?cái)嚯A次p＝0，得到非對(duì)稱轉(zhuǎn)子在對(duì)稱支承剛度下的頻率轉(zhuǎn)速變化曲線。當(dāng)各階振型的正進(jìn)動(dòng)頻率（圖中a0，c0和e0）為0時(shí)，轉(zhuǎn)子進(jìn)入不穩(wěn)定區(qū)。需要說明的是，在轉(zhuǎn)速為0時(shí)，轉(zhuǎn)子僅在第四階振型上表現(xiàn)出明顯的頻率分離，前三階振型在該轉(zhuǎn)速下對(duì)轉(zhuǎn)子連接軸段的局部彎曲剛度變化不敏感。
圖6 p＝0時(shí)低壓轉(zhuǎn)子的頻率轉(zhuǎn)速變化曲線Fig．6 Evolution of nature frequencies of the low pressure rotor for p＝0
圖7 p＝1時(shí)低壓轉(zhuǎn)子的頻率轉(zhuǎn)速變化曲線Fig．7 Evolution of nature frequencies of the low pressure rotor for p＝1
圖7中，截?cái)嚯A次p＝1，除得到轉(zhuǎn)子的基礎(chǔ)頻率外，還得到轉(zhuǎn)子振動(dòng)類模態(tài)的一階階次頻率。從圖7中可以看出，考慮轉(zhuǎn)子的一階階次頻率后，轉(zhuǎn)子的頻率轉(zhuǎn)速變換曲線變得更加復(fù)雜。圖6中的每條基礎(chǔ)頻率w0，n曲線均擴(kuò)展為w－1，n，w0，n，w＋1，n三條曲線，并且不是簡(jiǎn)單的擴(kuò)展，從圖7中可以看到，階次頻率與基礎(chǔ)頻率以及階次頻率之間出現(xiàn)多處的耦合和頻率轉(zhuǎn)向現(xiàn)象。階次頻率的耦合增加了轉(zhuǎn)子的不穩(wěn)定區(qū)范圍。因此，與僅考慮單一非對(duì)稱的轉(zhuǎn)子不同，考慮轉(zhuǎn)子的雙重非對(duì)稱特征時(shí)，轉(zhuǎn)子的不穩(wěn)定區(qū)區(qū)間增多，轉(zhuǎn)子的安全工作轉(zhuǎn)速范圍變小。
另外需要說明的是，由于頻率轉(zhuǎn)向與階次頻率耦合，各階頻率與振型的對(duì)應(yīng)關(guān)系非常復(fù)雜，并且分析人員更加關(guān)注的是不穩(wěn)定區(qū)的范圍，進(jìn)行詳細(xì)的振型區(qū)分也沒有重要的意義，因此在圖7中用相同的符號(hào)標(biāo)示各階頻率變化曲線。
5 結(jié) 論
本文提出了基于3D實(shí)體有限單元建模的雙重非對(duì)稱轉(zhuǎn)子支承系統(tǒng)的振動(dòng)特性分析方法，詳細(xì)給出了系統(tǒng)時(shí)變運(yùn)動(dòng)微分方程的組集推導(dǎo)過程以及頻域特性的分析求解方法。文中以一個(gè)具有簡(jiǎn)單結(jié)構(gòu)的雙重非對(duì)稱轉(zhuǎn)子支承模型為例說明了方法的分析流程，對(duì)ANSYS模型分析矩陣的輸出和MATLAB中支承剛度及約束施加的方法進(jìn)行了必要的說明。對(duì)某型航空發(fā)動(dòng)機(jī)低壓轉(zhuǎn)子的振動(dòng)特性分析進(jìn)一步說明的文中所提方法的工程實(shí)用價(jià)值。
參考文獻(xiàn)
［1］Lee C W．Vibration analysis of rotors［M］．Dordrecht：Kluwer Academic Publishers，1993．
［2］Combescure D，Lazarus A．Refined finite element modelling forthevibrationanalysisoflargerotatingmachines：application to the gas turbine modular helium reactor power conversion unit［J］．Journal of Sound and Vibration，2008，318（4－5）：1262－1280．
［3］陳萌，洪杰，朱彬，等．基于實(shí)體單元的轉(zhuǎn)子動(dòng)力特性計(jì)算方法［J］．北京航空航天大學(xué)學(xué)報(bào)，2007，33（1）：10－13．
CHEN Meng，HONG Jie，ZHU Bin，et al．Rotor dynamic analysis based on solid element［J］．Journal of Beijing University of Aeronautics and Astronautics，2007，33（1）：10－13．
［4］Chatelet E，D'Ambrosio F，Jacquet-Richardet G．Toward global modelling approaches for dynamic analyses of rotating assemblies of turbomachines［J］．Journal of Sound and Vibration，2005，282（1－2）：163－178．
［5］Genta G．Dynamics of rotating systems［M］．springer，2005．
［6］Black H F．Parametrically excited lateral vibrations of an asymmetric slender shaft in asymmetrically flexible bearings ［J］．Journal of Mechanical Engineering Science，1969，11 （1）：57－67．
［7］Yamamoto T，Ota H，Kono K．On the unstable vibrations of a shaft with unsymmetrical stiffness carrying an unsymmetrical rotor［J］．Journal of Applied Mechanics，1968，35（2）：313 －321．
［8］Gladwell G M L，Stammers C W．On the stability of an unsymmetrical rigid rotor supported in unsymmetrical bearings ［J］．Journal of Sound and Vibration，1966，3：221－232．
［9］Genta G．Whirling of unsymmetrical rotors：a finite element approach based on complex coordinates［J］．Journal of Sound and Vibration，1988，124（1）：27－53．
［10］Lee C W，Suh J H，Hong S W．Modal analysis of rotor system using modulated co-ordinates for asymmetric rotor systemwithanisotropicstator［C］．8thInternational Conference on Vibrations in Rotating Machinery．University ofSwansea，UK：ProfessionalEngineeringPublishing （Institution of Mechanical Engineers），2004：575－584．
［11］Lee C W，Han D J，Hong S W．Modal analysis of periodically time-varying linear rotor systems using floquet theory［C］．Proceedings of the 7th IFToMM-Conference on Rotor Dynamics．Austria：International Conference on Rotor Dynamics，2006．
［12］Nandi A，Neogy S．An efficient scheme for stability analysis of finite element asymmetric rotor models in a rotating frame ［J］．Finite Elements in Analysis and Design，2005，41 （14）：1343－1364．
［13］Nandi A．Reduction of finite element equations for a rotor model on non-isotropic spring support in a rotating frame［J］．Finite Elements in Analysis and Design，2004，40（9）：935 －952．［14］Lazarus A，Prabel B，Combescure D．A 3d finite element model forthevibrationanalysisofasymmetricrotating machines［J］．Journal of Sound and Vibration，2010，329 （18）：3780－3797．
［15］Kim K T，Lee C W．Dynamic analysis of asymmetric bladedrotors supported by anisotropic stator［J］．Journal of Sound and Vibration，2012，331（24）：5224－5246．
［16］沈松，鄭兆昌，應(yīng)懷樵．非對(duì)稱轉(zhuǎn)子－軸承－基礎(chǔ)系統(tǒng)的非線性振動(dòng)［J］．振動(dòng)與沖擊，2004，23（4）：33－35．
SHENSong，ZHENGZhao-chang，YINGHuai-qiao．Nonlinear vibration of unsymmetrical rotor-bearing-foundation system［J］．Journal of Vibration and Shock，2004，23（4）：33－35．
［17］肖錫武，肖光華，楊叔子．不對(duì)稱轉(zhuǎn)子系統(tǒng)的參激強(qiáng)迫振動(dòng)［J］．振動(dòng)工程學(xué)報(bào)，2002（3）：71－74．
XIAO Xi-wu，XIAO Guang-hua，YANG Shu-zi．Parametrically excited and forced oscillations of an unsymmetrical rotor system［J］．Journal of Vibration Engineering，2002（3）：71－74．
［18］肖錫武，楊正茂，肖光華，等．不對(duì)稱轉(zhuǎn)子系統(tǒng)的非線性振動(dòng)［J］．華中科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2002（5）：81 －84．
XIAO Xi-wu，YANG Zheng-mao，XIAO Guang-hua，et al． Nonlinear oscillations of an unsymmetrical rotor system［J］．J．Huazhong Univ．of Sci．＆Tech．（Nature Science Edition），2002（5）：81－84．
［19］王培杰，滕弘飛，劉篤金，等．非對(duì)稱復(fù)雜轉(zhuǎn)子系統(tǒng)有限元?jiǎng)恿Ψ治觯跩］．固體力學(xué)學(xué)報(bào)，1991（4）：370－376．
WANG Pei-jie，TENG Hong-fei，LIU Du-jin，et al．The dynamic analysis for complex asymmetric rotor-bearing systems by finite element method［J］．Acta Mechanica Solida Sinica，1991（4）：370－376．
［20］Lee C W，Han D J，Suh J H，et al．Modal analysis of periodically time-varying linear rotor systems［J］．Journal of Sound and Vibration，2007，303（3－5）：553－574．
［21］Han D J．Modal analysis of periodically time-varying rotor systems using floquet theory and modulated coordinate transformation：applications to crack detection and modal balancing［D］．Daejeon：Korea Advanced Institute of Science and Technology，2006．
［22］Boru F E．Numerical and experimental response and stability investigations of anisotropic rotor-bearing systems［D］．Kassel：kassel university，2010．
［23］Rao J S，Sreenivas R．Dynamics of a three level rotor system using solid elements［C］．Proceedings of ASME Turbo Expo 2003．Atlanta，Georgia，USA：ASME，2003，601－606．
Dynamic analysis method of an anisotropic and asymmetric rotor-bearing system based on 3D FEM
MA Wei-meng，WANG Jian-jun
（School of Energy and Power Engineering，Beijing University of Aeronautics and Astronautics，Beijing 100191，China）
Dynamic analysis method of an anisotropic and asymmetric rotor-bearing system based on 3D finite element method was presented．The system mass，stiffness and Coriolis force matrices expressed with respect to a rotating coordinate system were generated with ANSYS．Time-varying bearing stiffness matrix with respect to the rotating coordinates was then obtained with a coordinate transformation．The differential equations of motion of the system were generated by applying constraints and time-varying bearing stiffness to the system matrices using the corresponding relationships between node number and row index of the matrices．Floquet theory and Hill infinite determinant method were adopted to solve the differential equations，frequency characteristics and stability of the system were obtained．The proposed method was validated based on a model in literature，and the proposed method was applied to analyze the dynamic behavior of an industrial rotor considering anisotropic and asymmetric characters．
asymmetric；rotor-bearing system；3D finite element；time-varying system；stability
V 231．96
A
10．13465／j．cnki．jvs．2014．23．002
2013－07－16 修改稿收到日期：2014－01－02
馬威猛男，博士生，1987年1月生
王建軍男，教授，博士生導(dǎo)師，1956年生
猜你喜歡
汽車實(shí)用技術(shù)(2019年24期)2019-12-27 03:52:46
組合機(jī)床與自動(dòng)化加工技術(shù)(2019年7期)2019-08-06 03:51:06
數(shù)學(xué)年刊A輯(中文版)(2018年2期)2019-01-08 01:59:50
價(jià)值工程(2017年28期)2018-01-23 20:48:29
智富時(shí)代(2017年4期)2017-04-27 17:08:47
廣東石油化工學(xué)院學(xué)報(bào)(2016年6期)2016-05-17 05:17:43
數(shù)學(xué)理論與應(yīng)用(2016年4期)2016-05-17 04:50:23
中國鐵道科學(xué)(2015年4期)2015-06-21 06:46:08
電測(cè)與儀表(2015年4期)2015-04-12 00:43:04
空間控制技術(shù)與應(yīng)用(2010年3期)2010-12-23 08:04:54
