高精密軋機支承輥輥系動力特性有限元法分析

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

（中國重型機械研究院股份公司 陜西 西安 710032）

1 前言

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

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

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

支承輥輥系單元是一種典型的加入復雜擾動的軸承轉子系統［23-27］，對該復雜系統動力學的主要分析方法可分為兩類：傳遞矩陣法和有限元法。傳遞矩陣法具有占用計算機資源少、計算速度快及數值穩定性差等特點。

有限元法將無限自由度的連續系統簡化為有限自由度的離散系統，將各種因素考慮在內的有限元模型是比較精確的模型，計算結果精度高，但是占用更多的計算機資源，若根據具體問題在不影響分析精度時，忽略不必要的非線性等因素，將非線性系統線性化將會明顯提高計算效率減小計算機資源占用。隨計算機的發展，有限元的建模方法使得列寫乃至求解復雜精密的軸承轉子結構系統運動方程成為可能。

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

輥系裝配單元是一種典型的加入復雜擾動軸承轉子系統，通常軸承轉子系統可以沿軸線化分為在節點處聯結的離散圓盤、分布質量的彈性軸段及軸承座等單元，單盤轉子系統如圖1所示。

圖1 單盤轉子系統

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

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

節點的位移向量定義為通過該節點且與轉子軸線垂直截面的位移向量，有限元中，任一瞬時轉子的位置用該單元節點的位移來表示，各節點的位移向量組成了軸承轉子系統的廣義坐標。分析單元節點力與位移的關系，建立單元的運動方程，以節點位移為廣義坐標綜合各單元的運動方程可得到系統的運動微分方程。

2.1 質量圓盤的運動方程

對于剛性圓盤，若其質量、直徑轉動慣量及極轉動慣量分別為m、Jd和Jp。圓盤所在單元的節點位移向量為：

求得剛性圓盤的動能根據Lagrange方程可得剛性圓盤的運動微分方程：

式中： Md—圓盤的質量矩陣；

Ω［J］—回轉矩陣；

Q1d、Q1d—廣義力。

圓盤的質量矩陣及回轉矩陣如下：

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

2.2 彈性軸段的運動方程

對于彈性軸段單元，如圖2所示，將該單元的兩端的節點位移作為其廣義坐標：

圖2 彈性軸段單元

若軸段單元的單位長度質量、直徑轉動慣量及極轉動慣量分為為μ、jd及jp，通過求得該單元的動能及Lagrange方程可得該彈性軸段的運動方程：

式中： Ms—質量矩陣；

Ω［Js］—回轉矩陣；

Ks—剛度矩陣；

Q1s、Q1s—廣義力向量。

剛度矩陣為：

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

l—彈性軸段單元長度；

Ks—剛度矩陣。

通常，在轉子動力學中將無限多個自由度質量連續分布的彈性轉子系統簡化為沿軸線若干個集總質量的多自由度系統，集總質量的節點一般選擇在圓盤中心、軸頸中心及軸截面突變處等，如圖3所示。當節點間的彈性軸段為等截面軸時，質量、直徑轉動慣量及極轉動慣量集總如下：

圖3 轉子質量離散

式中：Mi、Jpi、Jdi—集總到節點i處質量、極轉動慣量及直徑轉動慣量；

μ、Jpi、Jdi—單位長度軸段的質量、極轉動慣量及直徑轉動慣量；

l—軸段長度。

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

2.3 軸承座的運動方程

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

圖4 軸承—轉子支承模型原理圖

若軸頸中心的編號為s（j），軸承中心及軸頸中心的坐標為（xb、yb）和（xs（j）、ys（j）），軸承座的運動方程為：

假設基礎剛性較好，軸承可簡化為各向同性且不計阻尼的等剛度彈性支撐，若其剛度系數為kb，其中kb與其承擔的載荷有關。則軸承作用在軸頸中心的廣義力為：

2.4 主軸單元的運動方程

將軸承轉子系統劃分為N個節點N-1個軸段組成的有限元模型，系統的位移向量為：

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

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

轉子系統的整體質量矩陣［M1］、回轉矩陣Ω［J1］及剛度矩陣［K1］均為2N×2N階矩陣，對于整體質量矩陣［M1］，將式（15）-（17）代入式（18）得到集總到節點i處質量矩陣，將其疊加到對角線上，表示圓盤及彈性軸段單元質量矩陣對整體質量矩陣的貢獻?；剞D矩陣Ω［J1］的形成和

將式（28）-（32）代入式（26）-（27），則軸承轉子系統的運動方程可表示為：整體質量矩陣的形成方法類似。對于剛度矩陣［K1］和整體質量矩陣的形成方法類似，由式（14）可得第N-1個節點與第N個節點間軸段的剛度矩陣，將對應節點處的元素疊加到對角線上。

3 輥系單元的臨界轉速與不平衡響應

綜合轉子系統圓盤、彈性軸段及軸承座等單元運動方程，建立軸承轉子系統的運動微分方程，通過求解系統的運動方程的齊次解，可得轉子自轉角速度Ω時的渦動頻率及轉子的臨界轉速及振型，采用Newmark-β方法求解方程的穩態解可得轉子系統不平衡響應。

系統的運動方程式（33）的齊次式為：

若：

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

令：

將軸承轉子系統的臨界轉速問題轉化為矩陣［D］的標準特征值v問題，即將式（39）-（40）代入式（38）可得：

矩陣［D］為8N×8N階矩陣，其特征值是與轉子的自轉角速度Ω有關的成對的共軛復數，它的虛部代表自轉角速度為Ω時的轉子正反渦動頻率。

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

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

若軸承轉子系統不計轉軸的偏心的影響，只考慮由圓盤偏心質量引起不平衡響應，將軸承轉子系統的運動方程式（33）寫成如下形式：

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

Qt—t時刻的圓盤偏心質量引起不平衡激勵。

在t+Δt瞬時，由式（42）可得：

假定在時間間隔［t，t+Δt］內的加速度為：

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

將 { Ut+Δt} 展開為二階Taylor展開式：

此時，又假定在時間間隔［t，t+Δt］內的加速度為：

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

若令：

通過選擇恰當的控制參數β1及β2，在已知t時刻系統的位移、速度及加速度的情況下，理論上可以求得t+Δt系統的位移、速度及加速度：

4 數值分析步驟及實例驗證

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

經狀態變量替換，轉子系統的臨界轉速及模態振型問題，轉化為求解式（42）特征值及特征向量問題，通過MATLAB程序中Eig函數可以得到轉子系統的臨界轉速，將特征值代入方程求解其特征向量可得對應階次的模態振型。通過Newmark法求解非齊次線性微分方程式（42）可得到在剛性圓盤的不平衡激勵作用下的不平衡響應。

為驗證臨界轉速計算的正確性，本文計算了《轉子動力學》中的兩端簡支的單盤轉子系統實例，如圖1所示，由于文獻中未考慮轉子質量，故將其密度設為ρ=1kg／m3，支撐為忽略阻尼的剛性支撐，剛度為K=7×106N／mm，具體參數如表1所示。

表1 單盤轉子系統參數

基于MATLAB編寫轉子系統的有限元法分析程序，將轉子系統等分為76個節點，圓盤位于第26個節點上，給定轉子自轉角速度Ω，通過求解式（41）的特征值，得到對應于該轉速的正進動頻率ωF1、ωF2及反進動頻率ωB1、ωB2，如表2所示。

表2 單盤轉子系統進動頻率 （rad／s）

表2 （續） 單盤轉子系統進動頻率（rad／s）

以轉子自轉角速度Ω作橫坐標，對應的正、反進動角速度為縱坐標，即將表2作圖，如圖5所示。

圖5 單盤轉子系統正反進動頻率

作直線Ω=±ω與曲線ωF1、ωB1及ωB2相交，得到按絕對值大小排序的轉子系統的臨界轉速，其中第二階為同步正向渦動時的臨界轉速，如表3所示。從表可以得出，本文臨界轉速的分析計算與文獻中的計算結果誤差很小，計算結果很精確，表明程序的分析計算結果可信。

單盤轉子系統對應于表3前三階的臨界轉速的模態振型通過求解式（41）的確定，如圖6所示。圖中前三階的臨界轉速下的模態振型在形狀和趨勢上與理論上的振型基本一致。

圖6 單盤轉子系統主振型

表3 單盤轉子系統臨界轉速 （rad／s）

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

圖7 單盤轉子系統不平衡響應

5 結論

（1）詳細介紹了主軸單元有限元建模分析方法，首先建立軸承轉子系統的有限元模型，求解系統的質量矩陣、回轉矩陣與剛度矩陣，得到系統的運動微分方程。

（2）使用MATLAB編寫轉子系統的有限元法分析程序，利用MATLAB程序中的Eig函數求解系統臨界轉速、Newmark-β方法求解系統運動微分方程不平衡響應。

（3）以兩端簡支的單盤轉子系統為實例驗證了MATLAB數值分析的正確性。