高精密軋機(jī)支承輥輥系動(dòng)力特性有限元法分析

王煜 閆成琨 王婉晴 謝犁 王果 馬彥芝

（中國(guó)重型機(jī)械研究院股份公司 陜西 西安 710032）

1 前言

高精密軋機(jī)在軋制生產(chǎn)過(guò)程中，支承輥輥系處于復(fù)雜的應(yīng)力狀態(tài)。如果支承輥輥系的選材、設(shè)計(jì)、制作工藝等不合理，或軋制時(shí)卡鋼等造成局部發(fā)熱引起熱沖擊等，都易使支承輥疲勞失效［1］。軋輥的尺寸結(jié)構(gòu)、材質(zhì)、使用在相當(dāng)程度上決定了軋機(jī)的技術(shù)水平，此外隨著軋制技術(shù)的發(fā)展，支承輥的工作環(huán)境也越來(lái)越苛刻。輥系軸承是影響支承輥輥系單元?jiǎng)討B(tài)特性最主要的部件［2］。角接觸球軸承、圓錐滾子軸承及圓柱滾子軸承等均可用作支承輥輥系的支撐，角接觸球軸承因其高速、高精度、高剛度、長(zhǎng)壽命及能同時(shí)承受軸向與徑向載荷等特點(diǎn)在支承輥輥系裝配單元上的應(yīng)用越來(lái)越廣泛［3，4］，輥系軸承的動(dòng)靜剛度等動(dòng)力學(xué)性能直接影響支承輥輥系靜態(tài)和動(dòng)態(tài)回轉(zhuǎn)精度、抗振、溫升以及速度等性能［5］，高速時(shí)鋼球的離心力及陀螺力矩等慣性效應(yīng)導(dǎo)致內(nèi)外圈接觸角、鋼球接觸載荷以及接觸區(qū)潤(rùn)滑狀態(tài)發(fā)生改變，使球軸承運(yùn)動(dòng)狀態(tài)的剛度與靜止?fàn)顟B(tài)的剛度產(chǎn)生差異［6］，同時(shí)，對(duì)球軸承適當(dāng)?shù)念A(yù)載，消除內(nèi)外圈之間的游隙，使鋼球和內(nèi)外圈產(chǎn)生彈性接觸變形，可增加承載鋼球的數(shù)量及使各鋼球受力趨于均勻，有利于提高支承輥輥系支撐軸承的剛度、輥系單元的回轉(zhuǎn)精度及抗振等性能［7，8］。

鑒于支承輥輥系單元?jiǎng)討B(tài)特性對(duì)高精密軋機(jī)整體精度的重要性［9］，在支承輥輥系單元的設(shè)計(jì)過(guò)程中對(duì)支承輥輥系單元和支撐軸承的性能及影響因素進(jìn)行準(zhǔn)確分析是保證所設(shè)計(jì)軋機(jī)的生產(chǎn)效率、加工精度及可靠性等性能的有力手段。對(duì)支承輥輥系單元的動(dòng)態(tài)特性分析，主要有模態(tài)綜合法、傳遞矩陣法及有限元法等方法，包括固有頻率、振型的模態(tài)分析、不平衡響應(yīng)及周期性外載荷作用下的諧響應(yīng)分析等內(nèi)容［10-12］，其中，如何正確簡(jiǎn)化軸承彈性支撐是支承輥輥系單元?jiǎng)討B(tài)分析的關(guān)鍵的技術(shù)難點(diǎn)［13-15］。

Gentle等人對(duì)純軸向載荷作用下角接觸球軸承，采用擬動(dòng)力學(xué)方法對(duì)其穩(wěn)態(tài)工作特性進(jìn)行了理論分析，利用試驗(yàn)驗(yàn)證了潤(rùn)滑介質(zhì)產(chǎn)生的拖動(dòng)力模型［16］。Nelias和Sainsot等人考慮了保持架對(duì)鋼球運(yùn)動(dòng)規(guī)律的影響，分析了球軸承的運(yùn)動(dòng)規(guī)律及摩擦功率損耗，建立了聯(lián)合載荷作用下角接觸球軸承擬動(dòng)力學(xué)分析模型，通過(guò)分析得到的保持架速度與試驗(yàn)值吻合較好，從而驗(yàn)證了該模型的合理性［17］。Walters首先提出了鋼球與保持架分別為四、六個(gè)自由度運(yùn)動(dòng)規(guī)律的動(dòng)力學(xué)分析模型［18］。Gupta在此基礎(chǔ)上詳細(xì)分析了球和溝道的相互作用，推導(dǎo)出球軸承中各零件的運(yùn)動(dòng)分析式作為數(shù)學(xué)模型，開(kāi)發(fā)了ADOBE程序?qū)崟r(shí)分析球軸承動(dòng)力學(xué)性能［19］。王黎欽建立了航空發(fā)動(dòng)機(jī)高速球軸承擬動(dòng)力學(xué)模型，分析了工況參數(shù)和結(jié)構(gòu)參數(shù)對(duì)滾動(dòng)軸承動(dòng)態(tài)特性的影響規(guī)律，通過(guò)Gupta的動(dòng)力學(xué)分析程序算例驗(yàn)證了程序的準(zhǔn)確性［20］。蔣興奇等對(duì)角接觸球軸承進(jìn)行了擬動(dòng)力學(xué)分析，探討更加符合實(shí)際的高速球軸承分析研究方法［21］。李錦標(biāo)等最早采用動(dòng)力學(xué)方法對(duì)滾子軸承進(jìn)行了分析計(jì)算，綜合考慮了滾子軸承各組件之間的作用力，得到了比較滿意的分析結(jié)果［22］。

支承輥輥系單元是一種典型的加入復(fù)雜擾動(dòng)的軸承轉(zhuǎn)子系統(tǒng)［23-27］，對(duì)該復(fù)雜系統(tǒng)動(dòng)力學(xué)的主要分析方法可分為兩類：傳遞矩陣法和有限元法。傳遞矩陣法具有占用計(jì)算機(jī)資源少、計(jì)算速度快及數(shù)值穩(wěn)定性差等特點(diǎn)。

有限元法將無(wú)限自由度的連續(xù)系統(tǒng)簡(jiǎn)化為有限自由度的離散系統(tǒng)，將各種因素考慮在內(nèi)的有限元模型是比較精確的模型，計(jì)算結(jié)果精度高，但是占用更多的計(jì)算機(jī)資源，若根據(jù)具體問(wèn)題在不影響分析精度時(shí)，忽略不必要的非線性等因素，將非線性系統(tǒng)線性化將會(huì)明顯提高計(jì)算效率減小計(jì)算機(jī)資源占用。隨計(jì)算機(jī)的發(fā)展，有限元的建模方法使得列寫乃至求解復(fù)雜精密的軸承轉(zhuǎn)子結(jié)構(gòu)系統(tǒng)運(yùn)動(dòng)方程成為可能。

2 輥系裝配單元有限元法建模

輥系裝配單元是一種典型的加入復(fù)雜擾動(dòng)軸承轉(zhuǎn)子系統(tǒng)，通常軸承轉(zhuǎn)子系統(tǒng)可以沿軸線化分為在節(jié)點(diǎn)處聯(lián)結(jié)的離散圓盤、分布質(zhì)量的彈性軸段及軸承座等單元，單盤轉(zhuǎn)子系統(tǒng)如圖1所示。

圖1 單盤轉(zhuǎn)子系統(tǒng)

若A及B位置為軸承支撐，轉(zhuǎn)子軸的中心線通過(guò)圓盤中心并以等角速度Ω自轉(zhuǎn)，忽略轉(zhuǎn)子的扭轉(zhuǎn)變形，轉(zhuǎn)子的自轉(zhuǎn)角φ：

以A為坐標(biāo)原點(diǎn)，以軸承座中心線為As軸建立固定坐標(biāo)系A(chǔ)xys，圖中節(jié)點(diǎn)o'的坐標(biāo)為（x，y），產(chǎn)生微小變形節(jié)點(diǎn)處與轉(zhuǎn)子軸線垂直的橫截面的偏轉(zhuǎn)角為（θx，θy），在坐標(biāo)系A(chǔ)xys中，任一橫截面位移向量可描述如下：

節(jié)點(diǎn)的位移向量定義為通過(guò)該節(jié)點(diǎn)且與轉(zhuǎn)子軸線垂直截面的位移向量，有限元中，任一瞬時(shí)轉(zhuǎn)子的位置用該單元節(jié)點(diǎn)的位移來(lái)表示，各節(jié)點(diǎn)的位移向量組成了軸承轉(zhuǎn)子系統(tǒng)的廣義坐標(biāo)。分析單元節(jié)點(diǎn)力與位移的關(guān)系，建立單元的運(yùn)動(dòng)方程，以節(jié)點(diǎn)位移為廣義坐標(biāo)綜合各單元的運(yùn)動(dòng)方程可得到系統(tǒng)的運(yùn)動(dòng)微分方程。

2.1 質(zhì)量圓盤的運(yùn)動(dòng)方程

對(duì)于剛性圓盤，若其質(zhì)量、直徑轉(zhuǎn)動(dòng)慣量及極轉(zhuǎn)動(dòng)慣量分別為m、Jd和Jp。圓盤所在單元的節(jié)點(diǎn)位移向量為：

求得剛性圓盤的動(dòng)能根據(jù)Lagrange方程可得剛性圓盤的運(yùn)動(dòng)微分方程：

式中： Md—圓盤的質(zhì)量矩陣；

Ω［J］—回轉(zhuǎn)矩陣；

Q1d、Q1d—廣義力。

圓盤的質(zhì)量矩陣及回轉(zhuǎn)矩陣如下：

若圓盤有微小的偏心距eξ及eη，剛性圓盤的不平衡力表示成廣義力為：

2.2 彈性軸段的運(yùn)動(dòng)方程

對(duì)于彈性軸段單元，如圖2所示，將該單元的兩端的節(jié)點(diǎn)位移作為其廣義坐標(biāo)：

圖2 彈性軸段單元

若軸段單元的單位長(zhǎng)度質(zhì)量、直徑轉(zhuǎn)動(dòng)慣量及極轉(zhuǎn)動(dòng)慣量分為為μ、jd及jp，通過(guò)求得該單元的動(dòng)能及Lagrange方程可得該彈性軸段的運(yùn)動(dòng)方程：

式中： Ms—質(zhì)量矩陣；

Ω［Js］—回轉(zhuǎn)矩陣；

Ks—?jiǎng)偠染仃嚕?/p>

Q1s、Q1s—廣義力向量。

剛度矩陣為：

式中： EI—彈性軸段抗彎截面系數(shù)；

l—彈性軸段單元長(zhǎng)度；

Ks—?jiǎng)偠染仃嚒?/p>

通常，在轉(zhuǎn)子動(dòng)力學(xué)中將無(wú)限多個(gè)自由度質(zhì)量連續(xù)分布的彈性轉(zhuǎn)子系統(tǒng)簡(jiǎn)化為沿軸線若干個(gè)集總質(zhì)量的多自由度系統(tǒng)，集總質(zhì)量的節(jié)點(diǎn)一般選擇在圓盤中心、軸頸中心及軸截面突變處等，如圖3所示。當(dāng)節(jié)點(diǎn)間的彈性軸段為等截面軸時(shí)，質(zhì)量、直徑轉(zhuǎn)動(dòng)慣量及極轉(zhuǎn)動(dòng)慣量集總?cè)缦拢?/p>

圖3 轉(zhuǎn)子質(zhì)量離散

式中：Mi、Jpi、Jdi—集總到節(jié)點(diǎn)i處質(zhì)量、極轉(zhuǎn)動(dòng)慣量及直徑轉(zhuǎn)動(dòng)慣量；

μ、Jpi、Jdi—單位長(zhǎng)度軸段的質(zhì)量、極轉(zhuǎn)動(dòng)慣量及直徑轉(zhuǎn)動(dòng)慣量；

l—軸段長(zhǎng)度。

若圓盤及彈性軸段集總到A和B兩點(diǎn)的質(zhì)量、直徑轉(zhuǎn)動(dòng)慣量及極轉(zhuǎn)動(dòng)慣量分別為mB、JdA、JdB、JpA及JpB，則質(zhì)量矩陣、回轉(zhuǎn)矩陣及剛度矩陣如下：

2.3 軸承座的運(yùn)動(dòng)方程

對(duì)于軸承座單元，將軸承—轉(zhuǎn)子系統(tǒng)中軸承簡(jiǎn)化為彈簧支撐模型，即軸承支撐等效為阻止徑向位移的線性彈簧及限制轉(zhuǎn)動(dòng)的扭轉(zhuǎn)彈簧固定在其支承中心，如圖4所示。

圖4 軸承—轉(zhuǎn)子支承模型原理圖

若軸頸中心的編號(hào)為s（j），軸承中心及軸頸中心的坐標(biāo)為（xb、yb）和（xs（j）、ys（j）），軸承座的運(yùn)動(dòng)方程為：

假設(shè)基礎(chǔ)剛性較好，軸承可簡(jiǎn)化為各向同性且不計(jì)阻尼的等剛度彈性支撐，若其剛度系數(shù)為kb，其中kb與其承擔(dān)的載荷有關(guān)。則軸承作用在軸頸中心的廣義力為：

2.4 主軸單元的運(yùn)動(dòng)方程

將軸承轉(zhuǎn)子系統(tǒng)劃分為N個(gè)節(jié)點(diǎn)N-1個(gè)軸段組成的有限元模型，系統(tǒng)的位移向量為：

綜合式（6）-（7）和式（14）-（15），即綜合剛性圓盤與彈性軸段單元運(yùn)動(dòng)方程，將軸承的支撐廣義力式（23）并入轉(zhuǎn)子系統(tǒng)的剛度矩陣相應(yīng)元素中，可得軸承轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)方程：

其中，Q不包含支承軸承反力，k為除2s（j）-1行及2s（j）-1列處為kb外，均為零的2N×2N階矩陣，若假設(shè)：

轉(zhuǎn)子系統(tǒng)的整體質(zhì)量矩陣［M1］、回轉(zhuǎn)矩陣Ω［J1］及剛度矩陣［K1］均為2N×2N階矩陣，對(duì)于整體質(zhì)量矩陣［M1］，將式（15）-（17）代入式（18）得到集總到節(jié)點(diǎn)i處質(zhì)量矩陣，將其疊加到對(duì)角線上，表示圓盤及彈性軸段單元質(zhì)量矩陣對(duì)整體質(zhì)量矩陣的貢獻(xiàn)。回轉(zhuǎn)矩陣Ω［J1］的形成和

將式（28）-（32）代入式（26）-（27），則軸承轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)方程可表示為：整體質(zhì)量矩陣的形成方法類似。對(duì)于剛度矩陣［K1］和整體質(zhì)量矩陣的形成方法類似，由式（14）可得第N-1個(gè)節(jié)點(diǎn)與第N個(gè)節(jié)點(diǎn)間軸段的剛度矩陣，將對(duì)應(yīng)節(jié)點(diǎn)處的元素疊加到對(duì)角線上。

3 輥系單元的臨界轉(zhuǎn)速與不平衡響應(yīng)

綜合轉(zhuǎn)子系統(tǒng)圓盤、彈性軸段及軸承座等單元運(yùn)動(dòng)方程，建立軸承轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)微分方程，通過(guò)求解系統(tǒng)的運(yùn)動(dòng)方程的齊次解，可得轉(zhuǎn)子自轉(zhuǎn)角速度Ω時(shí)的渦動(dòng)頻率及轉(zhuǎn)子的臨界轉(zhuǎn)速及振型，采用Newmark-β方法求解方程的穩(wěn)態(tài)解可得轉(zhuǎn)子系統(tǒng)不平衡響應(yīng)。

系統(tǒng)的運(yùn)動(dòng)方程式（33）的齊次式為：

若：

由式（35）-（37）將二階微分方程式（34）化為8N個(gè)一階線性方程：

令：

將軸承轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速問(wèn)題轉(zhuǎn)化為矩陣［D］的標(biāo)準(zhǔn)特征值v問(wèn)題，即將式（39）-（40）代入式（38）可得：

矩陣［D］為8N×8N階矩陣，其特征值是與轉(zhuǎn)子的自轉(zhuǎn)角速度Ω有關(guān)的成對(duì)的共軛復(fù)數(shù)，它的虛部代表自轉(zhuǎn)角速度為Ω時(shí)的轉(zhuǎn)子正反渦動(dòng)頻率。

對(duì)于自由度比較多的軸承轉(zhuǎn)子系統(tǒng)，求解矩陣的特征值比較耗時(shí)。通常，沒(méi)有必要求解出所有的特征值而只需要求解出系統(tǒng)的前若干階。通過(guò)式（34），求解得轉(zhuǎn)子的各種自轉(zhuǎn)角速度Ω時(shí)的渦動(dòng)頻率ωr并將其按序排列，以自轉(zhuǎn)角速度Ω為橫坐標(biāo)，渦動(dòng)頻率ωr為縱坐標(biāo)，作出系列曲線，作直線Ω=±ω與曲線相交可得到軸承轉(zhuǎn)子系統(tǒng)的同步正向渦動(dòng)頻率及同步反向渦動(dòng)頻率。由于實(shí)際轉(zhuǎn)子系統(tǒng)不平衡質(zhì)量的激勵(lì)，轉(zhuǎn)子作同步正向渦動(dòng)，因此，轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速通常只考慮同步正向渦動(dòng)的臨界轉(zhuǎn)速。

對(duì)于求得的每一階臨界轉(zhuǎn)速，代入式（41），總可得8N個(gè)齊次線性方程組，令｛V0｝中的某元素為1，從方程組中任選（8N-1）個(gè)方程，求解其組成的方程組便可求得轉(zhuǎn)子系統(tǒng)相應(yīng)階的主振型即模態(tài)振型，主振型上可以直觀的看出不同位置相對(duì)振動(dòng)幅度分布規(guī)律。

若軸承轉(zhuǎn)子系統(tǒng)不計(jì)轉(zhuǎn)軸的偏心的影響，只考慮由圓盤偏心質(zhì)量引起不平衡響應(yīng)，將軸承轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)方程式（33）寫成如下形式：

式中： ¨Ut、˙Ut、Ut—t時(shí)刻的加速度、速度及位移；

Qt—t時(shí)刻的圓盤偏心質(zhì)量引起不平衡激勵(lì)。

在t+Δt瞬時(shí)，由式（42）可得：

假定在時(shí)間間隔［t，t+Δt］內(nèi)的加速度為：

將式（44）代入式（45），可得：

將 { Ut+Δt} 展開(kāi)為二階Taylor展開(kāi)式：

此時(shí)，又假定在時(shí)間間隔［t，t+Δt］內(nèi)的加速度為：

將式（48）代入式（47），可得：

若令：

通過(guò)選擇恰當(dāng)?shù)目刂茀?shù)β1及β2，在已知t時(shí)刻系統(tǒng)的位移、速度及加速度的情況下，理論上可以求得t+Δt系統(tǒng)的位移、速度及加速度：

4 數(shù)值分析步驟及實(shí)例驗(yàn)證

通過(guò)分析軸承的載荷，可得角接觸軸承的剛度系數(shù)，綜合含不平衡質(zhì)量激勵(lì)的剛性圓盤的運(yùn)動(dòng)微分方程式（6）-（7）、通過(guò)集總質(zhì)量模型建立彈性軸段的運(yùn)動(dòng)微分方程式（14）-（15）及考慮載荷對(duì)軸承剛度影響的軸承座的運(yùn)動(dòng)微分方程式（22），建立質(zhì)量陣［M］、剛度陣［K］及回轉(zhuǎn)矩陣［G］，得到轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)微分方程式（33）。

經(jīng)狀態(tài)變量替換，轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速及模態(tài)振型問(wèn)題，轉(zhuǎn)化為求解式（42）特征值及特征向量問(wèn)題，通過(guò)MATLAB程序中Eig函數(shù)可以得到轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速，將特征值代入方程求解其特征向量可得對(duì)應(yīng)階次的模態(tài)振型。通過(guò)Newmark法求解非齊次線性微分方程式（42）可得到在剛性圓盤的不平衡激勵(lì)作用下的不平衡響應(yīng)。

為驗(yàn)證臨界轉(zhuǎn)速計(jì)算的正確性，本文計(jì)算了《轉(zhuǎn)子動(dòng)力學(xué)》中的兩端簡(jiǎn)支的單盤轉(zhuǎn)子系統(tǒng)實(shí)例，如圖1所示，由于文獻(xiàn)中未考慮轉(zhuǎn)子質(zhì)量，故將其密度設(shè)為ρ=1kg／m3，支撐為忽略阻尼的剛性支撐，剛度為K=7×106N／mm，具體參數(shù)如表1所示。

表1 單盤轉(zhuǎn)子系統(tǒng)參數(shù)

基于MATLAB編寫轉(zhuǎn)子系統(tǒng)的有限元法分析程序，將轉(zhuǎn)子系統(tǒng)等分為76個(gè)節(jié)點(diǎn)，圓盤位于第26個(gè)節(jié)點(diǎn)上，給定轉(zhuǎn)子自轉(zhuǎn)角速度Ω，通過(guò)求解式（41）的特征值，得到對(duì)應(yīng)于該轉(zhuǎn)速的正進(jìn)動(dòng)頻率ωF1、ωF2及反進(jìn)動(dòng)頻率ωB1、ωB2，如表2所示。

表2 單盤轉(zhuǎn)子系統(tǒng)進(jìn)動(dòng)頻率 （rad／s）

表2 （續(xù)） 單盤轉(zhuǎn)子系統(tǒng)進(jìn)動(dòng)頻率（rad／s）

以轉(zhuǎn)子自轉(zhuǎn)角速度Ω作橫坐標(biāo)，對(duì)應(yīng)的正、反進(jìn)動(dòng)角速度為縱坐標(biāo)，即將表2作圖，如圖5所示。

圖5 單盤轉(zhuǎn)子系統(tǒng)正反進(jìn)動(dòng)頻率

作直線Ω=±ω與曲線ωF1、ωB1及ωB2相交，得到按絕對(duì)值大小排序的轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速，其中第二階為同步正向渦動(dòng)時(shí)的臨界轉(zhuǎn)速，如表3所示。從表可以得出，本文臨界轉(zhuǎn)速的分析計(jì)算與文獻(xiàn)中的計(jì)算結(jié)果誤差很小，計(jì)算結(jié)果很精確，表明程序的分析計(jì)算結(jié)果可信。

單盤轉(zhuǎn)子系統(tǒng)對(duì)應(yīng)于表3前三階的臨界轉(zhuǎn)速的模態(tài)振型通過(guò)求解式（41）的確定，如圖6所示。圖中前三階的臨界轉(zhuǎn)速下的模態(tài)振型在形狀和趨勢(shì)上與理論上的振型基本一致。

圖6 單盤轉(zhuǎn)子系統(tǒng)主振型

表3 單盤轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速 （rad／s）

單盤轉(zhuǎn)子系統(tǒng)剛性圓盤引起的不平衡響應(yīng)，可由軸承轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)微分方程式（42）通過(guò)Newmark解法確定。通常，實(shí)際轉(zhuǎn)子在剛性圓盤不平衡質(zhì)量的激勵(lì)下將作同步正向渦動(dòng)，因此，更加重視同步正向渦動(dòng)的臨界轉(zhuǎn)速，不考慮同步反進(jìn)動(dòng)的臨界轉(zhuǎn)速，如圖7所示，轉(zhuǎn)子在臨界轉(zhuǎn)速為ω=244.35rad／s附近幅值出現(xiàn)了峰值。

圖7 單盤轉(zhuǎn)子系統(tǒng)不平衡響應(yīng)

5 結(jié)論

（1）詳細(xì)介紹了主軸單元有限元建模分析方法，首先建立軸承轉(zhuǎn)子系統(tǒng)的有限元模型，求解系統(tǒng)的質(zhì)量矩陣、回轉(zhuǎn)矩陣與剛度矩陣，得到系統(tǒng)的運(yùn)動(dòng)微分方程。

（2）使用MATLAB編寫轉(zhuǎn)子系統(tǒng)的有限元法分析程序，利用MATLAB程序中的Eig函數(shù)求解系統(tǒng)臨界轉(zhuǎn)速、Newmark-β方法求解系統(tǒng)運(yùn)動(dòng)微分方程不平衡響應(yīng)。

（3）以兩端簡(jiǎn)支的單盤轉(zhuǎn)子系統(tǒng)為實(shí)例驗(yàn)證了MATLAB數(shù)值分析的正確性。