吝水林，孫建亮，彭 艷

(燕山大學國家冷軋板帶裝備及工藝工程技術(shù)研究中心，河北，秦皇島 066004)

軋機滾動軸承是軋機裝備的關(guān)鍵零部件，其性能和質(zhì)量是保證軋機穩(wěn)定運行的關(guān)鍵。軋輥軸承通常在高速高溫的惡劣環(huán)境下工作，采用多列軸承來承受軋制負荷，輥身和輥頸由于軋制負荷過大發(fā)生彈性彎曲，徑向載荷、軸向載荷和力矩載荷的聯(lián)合作用導致各列滾動體和滾道單位接觸長度上載荷非均勻分布。由于服役工況下軋機部件的磨損，現(xiàn)代軋制工藝的改進，軋制力、彎輥力、軸向力的聯(lián)合作用，軸承的制造誤差、軸承座的安裝精度等綜合影響，導致軸承在總負載不變的情況下各列滾動體承載不均，此外，推力載荷進一步加劇了滾動體與滾道接觸應(yīng)力的非均勻性，導致軸承長期在偏載狀態(tài)下運行，大大降低了軸承的耐久性和疲勞壽命，嚴重影響軋機運行穩(wěn)定性和生產(chǎn)效率。因此，有必要研究高速板帶，軋機偏載狀態(tài)下軋機軸承力學性能，對保證軋機穩(wěn)定運行和提高軸承服役壽命有重要意義。

針對滾動軸承多種多樣的應(yīng)用領(lǐng)域及工況條件，研究者們從擬靜/動力學模型[1?4]、動力學模型[5?7]、故障力學模型[8?15]、接觸潤滑理論及力學性能方面開展研究，力求構(gòu)建軸承在復雜工況下運行的力學特性與服役狀態(tài)，更好地服務(wù)于軸承設(shè)計與應(yīng)用。

在擬靜力學、擬動力學、動力學模型方面，Hertz[16]提出了空間點接觸理論和平面無限長線接觸理論，建立了接觸問題的基本方程，該理論在軋機滾動軸承靜力學分析上得到了廣泛應(yīng)用。Jones[17]在球軸承力學模型中考慮了鋼球的離心力和陀螺力矩，并提出了套圈控制理論，通過Newton-Raphson 迭代法得到了滾動體的運動參數(shù)、載荷分布以及軸承剛度等重要參數(shù)。Rumbarger 和Filetti[18]將彈流潤滑理論引入擬動力學模型中，建立了軸承動力學分析模型。張彪等[19]針對森基米爾軋機軸承提出了考慮軸承套圈變形的載荷計算分析方法，對于多輥軋機支承輥軸承等載荷計算，以Hertz 理論為基礎(chǔ)的彈性變形理論難以滿足，必須考慮套圈彎曲變形的影響。史修江和王黎欽[20]以D1842926 航空發(fā)動機主軸滾子軸承典型工況為算例，基于擬動力學分析結(jié)果，獲得滾動體與套圈之間的接觸微區(qū)運動和受力狀態(tài)。陳芳華和汪久根[21]考慮接觸體表面的彈性變形效應(yīng)與潤滑劑的壓粘效應(yīng)，建立圓錐滾子軸承的彈性流體動力潤滑模型，并進行了完全數(shù)值解分析，得到了流體動壓力分布、油膜形狀與表面層內(nèi)Mises 應(yīng)力分布。

在求解軸承載荷分布結(jié)合有限元分析仿真方面，王思明等[22]在不考慮支撐結(jié)構(gòu)剛度的情況下根據(jù)受力和變形的關(guān)系推導出變槳軸承溝道上載荷分布的計算公式，并計算出變槳軸承在某一負游隙時具有較好的力學性能。劉光明等[23]采用ANSYS/LS-DYNA 有限元軟件，建立了包含支承輥軸承及軸承座組件、軋輥和軋件的四輥軋機三維實體同步耦合有限元分析模型，通過對比和分析咬入階段、穩(wěn)定軋制階段和拋鋼階段各列滾動體的受力情況，得到了軋制過程中軸承載荷的分布規(guī)律。高學海等[24]不考慮軸承支撐結(jié)構(gòu)和軸承套圈變形建立了雙列四點接觸轉(zhuǎn)盤軸承的載荷分布與接觸角分布模型，討論了軸承幾何參數(shù)對軸承承載能力的影響。聶重陽等[25]建立了線接觸載荷作用下的界面力學分析模型。趙興乾等[26]提出“徑向軸承載荷增量”這一概念實現(xiàn)對整體隔振系統(tǒng)中軸系靜態(tài)安全性定量分析，并計算了在不同螺旋槳推力載荷作用下及船舶縱傾工況下的軸承載荷增量。王興東等[27]利用有限元分析軟件建立低速重載軸承全模型，分析了軸承滾動體與溝道間的接觸應(yīng)力，軸承內(nèi)、外圈的應(yīng)力與應(yīng)變，為分析軸承溝道疲勞壽命的應(yīng)力計算奠定基礎(chǔ)。Singh 等[28]使用顯式動力學有限元軟件包LS-DYNA，數(shù)值求解了帶有外滾道缺陷的滾動軸承的動態(tài)非線性有限元模型。

在力學性能分析方面，曹宏瑞等[29]通過建立高速滾動軸承力學模型，基于外圈圓心位置不變假設(shè)，分析了滾動體受離心力、陀螺力矩、軸承內(nèi)、外圈受熱膨脹變形時，涉及軸承中滾珠變形狀態(tài)和方位角的軸承時變剛度。Cavallaro 等[30]建立了高速柔性套圈圓柱滾子軸承的分析模型，研究了套圈的變形、軸承內(nèi)部的受力分布、能量損耗以及滾子滑動速度與載荷分布的關(guān)系等。

上述研究從軸承的力學模型、力學性能、仿真分析等方面做了大量工作，但在軋機輥系運行過程中，未考慮配套軸承在惡劣的服役工況下以偏載狀態(tài)運行時的動態(tài)行為和應(yīng)力分布，本文針對機架輥系剛度偏差、來料楔形、非對稱軋制工藝等因素引起的高速板帶軋機輥系軸承的偏載力學行為，建立對應(yīng)力學模型，研究了滾動軸承在高速偏載下的力學性能。

1 Hertz 理想接觸模型

Hertz 接觸理論是研究兩物體因受壓相觸后產(chǎn)生的局部應(yīng)力和應(yīng)變分布規(guī)律的理論方法。對于滾子軸承，滾動體與內(nèi)、外圈的接觸可以看作是兩曲面物體之間的相互擠壓。Hertz 理論有如下假設(shè)：1)接觸區(qū)發(fā)生小變形；2)接觸面呈橢圓形；3)相接觸的物體被看作是彈性半空間，接觸面只有垂直分布壓力。

1.1 彈性接觸區(qū)和最大接觸應(yīng)力模型

圖1 為彈性半空間內(nèi)均布載荷分布，q 為作用在圓柱體半空間的均布載荷，R1和R2分別為兩個圓柱體半徑。由Hertz 線接觸理論可知，在q 的作用下，兩圓柱體的中心軸線相互趨近，在接觸界面會產(chǎn)生寬度為2a 的接觸區(qū)，如圖2 所示，P0為均布載荷下的最大接觸應(yīng)力，a 為接觸界面半寬度，實際運行中滾動體與內(nèi)、外滾道接觸時會分別產(chǎn)生上、下兩個接觸區(qū)，如圖3 所示，最大接觸應(yīng)力P1和P2、接觸半寬度a1與a2的計算如下。

圖 1 彈性半空間均布載荷分布Fig. 1 Uniform loads distribution in elastic half space

圖 2 接觸界面放大圖Fig. 2 Enlarged contact interface

圖 3 半橢圓接觸示意圖Fig. 3 Schematic of semi-elliptical contact

1.2 接觸載荷與彈性趨近量關(guān)系

彈性趨近量即為滾動體與內(nèi)、外圈接觸時的彈性變形量，圖3 中滾動體與內(nèi)外圈接觸時的彈性趨近量通過赫茲理論的基本方程無法求解，大多數(shù)情況下難以找到閉合形式的理論解，大多采用修正的Palmgren 經(jīng)驗公式，它考慮了滾子與滾道的直徑以及曲率的凹凸性，如式(3)所示：

2 偏載接觸效應(yīng)的形成與建模

2.1 偏載接觸效應(yīng)的形成

圖4 所示為軋制負荷大、來料楔形、非對稱軋制、軸向力、軋機各部件磨損、軸承制造誤差及非對稱因素等導致軸承內(nèi)、外圈產(chǎn)生角度為θ 的傾斜現(xiàn)象，滾動體在軸向方向應(yīng)力不同，各列滾動體也由于軸頸的傾斜承受不同的軋制負荷，導致軸承長期處于偏載狀態(tài)下運行。軸承承受徑向載荷的傳遞路徑為輥徑到內(nèi)圈，內(nèi)圈到滾動體，滾動體到外圈，外圈到軸承座，本節(jié)主要研究在偏載狀態(tài)下滾動體與滾道的接觸效應(yīng)。

圖 4 軸承套圈的傾斜Fig. 4 Tilt of the bearing ring

2.2 偏載下有限長線接觸模型及數(shù)值算法

圖5 為偏載狀態(tài)下采用條形單元法將載荷接觸區(qū)進行離散化處理。四列圓柱滾子軸承的滾子與滾道接觸長度有限，圓柱滾子與內(nèi)、外圈的滾道發(fā)生傾斜形成偏載效應(yīng)時，無限長圓柱體之間理想化接觸狀態(tài)不再適用，超出了赫茲線接觸理論的范圍，屬于非赫茲接觸問題。對于非赫茲接觸問題，假設(shè)滾動體與滾道接觸區(qū)在平行于軸承徑向平面被劃分為一定數(shù)量的條形單元，由于接觸變形很小，可以忽略條形單元之間的切應(yīng)力，僅僅考慮彈性變形。建立各個單元的力平衡方程和彈性變形方程，取滾動體與內(nèi)圈的接觸區(qū)進行分析。單元j 中，在x 軸上任一點上的接觸應(yīng)力為：

式中：aj為j 單元寬度；hj為單元長度的一半；Zi(yi)為第i 個單元的y 軸幾何中心的表面函數(shù)值；Q 為總體載荷；E′為等效彈性模量；Dij為影響系數(shù)。式(9)構(gòu)成 n+1線性方程組，當Q 和Zi已知時，可以解出Poj和 δ共n+1未知數(shù)。

圖 5 接觸區(qū)條形單元Fig. 5 Bar units of contact area

影響系數(shù)Dij表示j 單元上的力對i 單元位移的影響，完善了傳統(tǒng)切片法的各單元之間不相關(guān)的缺點，使每個單元之間具有了相關(guān)性。影響系數(shù)具體計算如下：

式中：yi、yj分別為i、j 單元的中心坐標；h 為每個單元長度的一半，h=l/2n。

式(11)的進一步計算必須采用數(shù)值積分的方法。此外，對于基本方程組，應(yīng)滿足應(yīng)力非負的約束條件，即 Poj≥0 (j=1,2,3···)。

采用一維方法，將滾子離散為n 個單元，每個單元直接采用柔度矩陣法進行關(guān)聯(lián)，將有線長接觸的基本方程化為n+1 個方程組，方程組需迭代求解，首先對彈性趨近量δ 和接觸半寬aj賦初值，計算影響系數(shù)Dij、最大接觸應(yīng)力Pi，若在迭代過程中，Poj為負值，說明該單元未發(fā)生接觸，進行排除處理，再利用力平衡方程，對所有Poj求和，判斷其是否滿足平衡條件，即是否收斂。若不收斂，則調(diào)整δ，重新迭代，直至其收斂，圖6為計算流程圖，其中am是中間變量，儲存上一次循環(huán) aj的值，ε1、ε2為誤差控制小量，調(diào)整δ 時，采用線性誤差控制方案，此外aj的調(diào)整方案采用均值調(diào)整方案。

圖 6 有限長數(shù)值解法計算流程Fig. 6 Calculation process of finite length numerical solution

2.3 滾動體與滾道接觸模型及數(shù)值算法

圖7 為單列軸承界面，可以清楚看到徑向力及滾子位置角，徑向載荷下，在位置角為φi的滾動體處，其與內(nèi)、外圈接觸的彈性變形量為：

由此，可確定Qri和αi，通過2.2 節(jié)圓柱滾子軸承的有限長數(shù)值解法，可求解位置角處滾動體的接觸應(yīng)力分布。圖8 為單個滾子在力和力矩綜合作用下平衡示意圖，對于單個滾動體，求解力與力矩。

圖 7 單列軸承截面Fig. 7 Cross section of single row bearing

圖 8 力與力矩平衡Fig. 8 Force and moment balance

式中：n 為沿滾子長度劃分的單元數(shù)；qij為單元上的集中載荷；hj為單元長度一半；aij為接觸半寬度；Pij為每單元內(nèi)的最大接觸應(yīng)力；Fr為徑向載荷；yi為單元中心坐標；yc為滾動體中心坐標；z 為滾子個數(shù)，力學模型可通過迭代求解。

將上述模型編制計算機程序，如圖9 所示，可計算總徑向力與力矩、單個滾子的接觸載荷、接觸應(yīng)力與彈性變形等，具體流程如下：首先給定軋輥和軸承的材料、尺寸參數(shù)及輥間載荷q，通過分析軋輥，確定軋輥作用與軸承的總徑向載荷Fr和輥頸傾斜角θ。進入外循環(huán)對四列圓柱滾子軸承上每列軸承的等效徑向載荷Frk賦初值，初值設(shè)定為總載荷的1/4。進入內(nèi)循環(huán)對每列軸承的徑向位移δrk賦初值，通過方程求解每個滾子的接觸載荷并確定約束條件。之后代入2.2 節(jié)單個滾子有限長接觸的數(shù)值計算程序中，求得單個滾子的接觸應(yīng)力Pij、接觸半寬aij、寬度hj與集中力qij，利用數(shù)值計算的結(jié)構(gòu)，計算單列滾子等效徑向力和力矩，通過誤差控制小量εr判斷是否收斂，不收斂采用線性誤差控制方案調(diào)整δrk，收斂穩(wěn)定，循環(huán)迭代結(jié)束。

圖 9 計算流程圖Fig. 9 Calculation flowchart

3 理想及偏載接觸模擬仿真

建立有限元分析模型時既要如實反映軸承實際結(jié)構(gòu)的重要力學特性，又要采用較少單元數(shù)量以及簡單的單元形態(tài)。本研究采用四列圓柱滾子軸承FC4666170，圖10(a)為建立的三維有限元模型，在整體模型基礎(chǔ)上提取內(nèi)圈和滾動體接觸模型，建立基于Hertz 理論的理想接觸模型，圖10(b)為整體網(wǎng)格劃分，圖10(c)為滾動體網(wǎng)格劃分。

1)劃分單元及確定接觸類型。理想接觸模型采用8 節(jié)點四面體單元，通過控制接觸區(qū)域的網(wǎng)格密度來保證計算精度，采用手動劃分網(wǎng)格的方法，單元數(shù)量25 988，節(jié)點數(shù)80 196。三維有限元模型整體采用8 節(jié)點四面體單元和8 節(jié)點六面體單元耦合形式，滾動體與滾道接觸區(qū)域網(wǎng)格細分，滾動體等規(guī)則體采用8 節(jié)點六面體單元，非規(guī)則體采用8 節(jié)點四面體單元。單元數(shù)量176 933，節(jié)點數(shù)506 146。采用柔體-柔體的點-面接觸類型定義內(nèi)外滾道與滾動體之間的接觸。

2)設(shè)置接觸對、邊界條件及外載荷。理想接觸模型設(shè)置滾動體與內(nèi)滾道的接觸對。有限元模型設(shè)置四個接觸對：滾動體與內(nèi)、外滾道、軸承座與外圈、內(nèi)圈與輥徑。在理想接觸模型四象限位置留硬點。在Mechanical 中設(shè)置2D Behavior 為平面應(yīng)力類型，內(nèi)圈與大地接觸線添加固定約束，滾動體及接觸區(qū)創(chuàng)建位移約束，X 方向為0，Y 方向自由。偏載接觸模型中軸承座兩個側(cè)表面施加X 方向約束，軸承座前表面施加Z 方向約束，使其軸向固定，支承輥中間處設(shè)置對稱約束。理想接觸模型滾動體上施加徑向載荷。偏載接觸模型外載荷通過軸承座施加于軸承。

圖 10 有限元分析模型Fig. 10 Finite element analysis model

4 算例與分析

4.1 計算參數(shù)

以生產(chǎn)線的850 mm 四輥冷軋機支承輥為例，采用配套的四列圓柱滾子軸承FC4666170，支承輥與軸承的模型參數(shù)如表1 所示，以實際工況為算例。

表 1 模型參數(shù)Table 1 Model parameters

4.2 結(jié)果計算與分析

4.2.1 理想接觸模型結(jié)果對比

對理想接觸模型模擬，取0°位置角滾動體與內(nèi)圈接觸模型，圖11 為采用理想接觸模型進行有限元計算的位移分布圖與應(yīng)力分布圖。從圖11 中可以看到接觸區(qū)的半橢圓分布形式，與Hertz 接觸理論假設(shè)相同，接觸應(yīng)力上、下對稱分布，最大接觸應(yīng)力出現(xiàn)在接觸區(qū)中線附近。本次分析只考慮表面及表面附近的應(yīng)力所產(chǎn)生的變形以及由接觸應(yīng)力所產(chǎn)生的變形，由于滾動元件的剛度特性，本文采用的鋼制軸承變形一般小于0.025 mm，由表2 可知，接觸半寬度、最大接觸應(yīng)力、彈性趨近量等關(guān)鍵參數(shù)的誤差控制在實際工程要求范圍內(nèi)，有限元解與解析解吻合程度較好。

圖 11 有限元結(jié)果分布云圖Fig. 11 Cloud diagram of finite element results

4.2.2 徑向力與彈性變形量

圖12、圖13 為每個滾動體的總徑向載荷和彈性變形量，總軋制力在操作側(cè)軸承和傳動側(cè)軸承兩側(cè)對稱分布，軸承第1 列～第4 列的順序為靠近輥頸側(cè)到遠離輥頸側(cè)。輥間傾斜角為0.01°，通過計算可知軸承各列的總徑向力Fr分別為246.15 kN、232.05 kN、217.95 kN、203.85 kN，總力矩Mm分別為46.51 N·m、46.08 N·m、45.63 N·m、45.16 N·m。整體結(jié)構(gòu)呈對稱分布，取右側(cè)滾動體分析，最上端滾動體為0°位置角，滾動體所在位置角分別為0°～180°。同一列滾動體，隨著其所在位置角的增大，徑向載荷和彈性變形量逐漸減小，90°變?yōu)?°，滾動體與滾道開始脫離接觸，不再有接觸應(yīng)力和接觸變形。相同位置角的不同列滾動體，承受徑向載荷和變形量依次減小。

表 2 數(shù)值解析解與有限元解對比Table 2 Comparison of numerical analytical results withfinite element results

圖 12 單個滾子的總徑向載荷Fig. 12 Total radial load of single roller

圖 13 單個滾子的彈性變形量Fig. 13 Amount of elastic deformation of single roller

4.2.3 滾動體接觸應(yīng)力

圖14 為單個滾動體沿長度方向的接觸應(yīng)力分布，對于每一個滾動體，單元接觸應(yīng)力沿滾子長度依次減小，接觸應(yīng)力分布呈“半船形”，最大接觸應(yīng)力達到1500 MPa 左右，且第1 列～第4 列最大接觸應(yīng)力依次減小，也間接反映了列間偏載的狀態(tài)，滾動體的邊緣接觸應(yīng)力比中心位置大，本模型考慮了滾動體端部發(fā)生應(yīng)力集中現(xiàn)象，彌補了Harris 的簡單切片法不能反映滾動體端部應(yīng)力集中的弊端。

4.2.4 偏載接觸模型結(jié)果對比

圖15(a)為有限元模擬的單列滾動體周向最大接觸應(yīng)力，圖15(b)為利用迭代進行數(shù)值求解的周向最大接觸應(yīng)力結(jié)果。滾動體端面與內(nèi)滾道產(chǎn)生接觸變形，最大接觸應(yīng)力隨著滾動體位置角急劇變化，為了使分析過程更具有一般性，三維模型中滾動體在周向方向裝配是隨機的，模擬中0°位置恰好沒有滾動體與套圈接觸，15°～30°位置角的滾動體，承載力最大，整體呈“M”形分布，超過90°后，接觸應(yīng)力急劇降為零，滾動體與套圈沒有發(fā)生有效的接觸，即沒有起到承載作用。圖15(b)的數(shù)值計算結(jié)果中設(shè)置同樣邊界條件，結(jié)果顯示在0°位置角附近沒有數(shù)值顯示，說明在15°～30°處應(yīng)力數(shù)值已達到頂峰，與有限元分析中在0°位置急劇降為零趨勢相同。有限元分析結(jié)果中第1 列和第2 列接觸應(yīng)力大于第 3 列和第4 列，整體趨勢與數(shù)值求解結(jié)果相同，整個偏載過程可以很清晰看到。由于有限元分析更接近實際模型和工況，計算結(jié)果更精確，且略小于數(shù)值求解結(jié)果。

圖 15 滾動體周向最大接觸應(yīng)力分布Fig. 15 Maximum contact stress distribution on the roller in the circumferential direction

4.2.5 不同載荷分布對偏載接觸效應(yīng)的影響

以上介紹了在均布載荷q=2×106N 時軸承整體及單個滾子的受力及變形情況，以下將討論不同載荷分布的影響。圖16 為每列軸承在不同分布載荷下承受的總徑向力，隨著工況載荷的增大，同一列的總徑向力在不斷增大，四列徑向力線性分布的斜率增大，隨著均布載荷的增加，軸承上各列載荷的偏載程度越來越嚴重，致使軋制軸承長期處于嚴重偏載的狀態(tài)運行，大大降低了使用壽命。

圖 16 不同載荷分布時單列軸承徑向力對比Fig. 16 Radial force comparison of single-row bearings under different load distributions

5 結(jié)論

本文針對機架輥系剛度偏差、來料楔形、非對稱軋制工藝等因素引起的高速板帶軋機輥系軸承的偏載力學行為，建立對應(yīng)力學模型，研究了滾動軸承在高速偏載下的力學性能。針對輥間不同載荷進行了模型數(shù)值求解，仿真模擬了理想及偏載接觸過程。研究為延長軸承使用壽命、保證軋機安全穩(wěn)定運行提供理論依據(jù)，具有重要工程意義。主要結(jié)論包括：

(1)模擬了Hertz 線接觸問題。對理想接觸模型模擬，取0°位置角滾動體與內(nèi)圈接觸模型，通過有限元計算的位移分布圖與應(yīng)力分布圖得到接觸區(qū)的半橢圓分布形式，接觸應(yīng)力上、下對稱分布，最大接觸應(yīng)力出現(xiàn)在接觸區(qū)中線附近，解析解與有限元解誤差在15%以內(nèi)。

(2)同一列滾動體，隨著其所在位置角的增大，徑向載荷和彈性變形量逐漸減小，90°變?yōu)?°，滾動體與滾道開始脫離接觸，不再有接觸應(yīng)力和接觸變形。相同位置角的不同列滾動體，承受徑向載荷和變形量依次減小。

(3)滾動體在周向隨機裝配，在15°～30°位置角的滾動體，承載力最大，整體呈“M”形分布，超過90°后，接觸應(yīng)力急劇降為零，滾動體與套圈沒有發(fā)生有效的接觸，沒有起到承載作用。

(4)對比了不同載荷分布對偏載接觸效應(yīng)的影響，隨著軋制載荷的增大，軸承同一列滾子的總徑向力不斷增大，并且四列徑向力線性分布的斜率增大，隨著軋制載荷的增加，軸承各列載荷的偏置程度越來越嚴重。