偏載狀態下軋機軸承接觸力學模型研究

吝水林，孫建亮，彭 艷

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

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

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

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

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

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

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

1 Hertz 理想接觸模型

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

1.1 彈性接觸區和最大接觸應力模型

圖1 為彈性半空間內均布載荷分布，q 為作用在圓柱體半空間的均布載荷，R1和R2分別為兩個圓柱體半徑。由Hertz 線接觸理論可知，在q 的作用下，兩圓柱體的中心軸線相互趨近，在接觸界面會產生寬度為2a 的接觸區，如圖2 所示，P0為均布載荷下的最大接觸應力，a 為接觸界面半寬度，實際運行中滾動體與內、外滾道接觸時會分別產生上、下兩個接觸區，如圖3 所示，最大接觸應力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 接觸載荷與彈性趨近量關系

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

2 偏載接觸效應的形成與建模

2.1 偏載接觸效應的形成

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

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

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

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

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

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

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

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

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

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

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

2.3 滾動體與滾道接觸模型及數值算法

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

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

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

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

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

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

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

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

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

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

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

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

4 算例與分析

4.1 計算參數

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

表 1 模型參數Table 1 Model parameters

4.2 結果計算與分析

4.2.1 理想接觸模型結果對比

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

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

4.2.2 徑向力與彈性變形量

圖12、圖13 為每個滾動體的總徑向載荷和彈性變形量，總軋制力在操作側軸承和傳動側軸承兩側對稱分布，軸承第1 列～第4 列的順序為靠近輥頸側到遠離輥頸側。輥間傾斜角為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。整體結構呈對稱分布，取右側滾動體分析，最上端滾動體為0°位置角，滾動體所在位置角分別為0°～180°。同一列滾動體，隨著其所在位置角的增大，徑向載荷和彈性變形量逐漸減小，90°變為0°，滾動體與滾道開始脫離接觸，不再有接觸應力和接觸變形。相同位置角的不同列滾動體，承受徑向載荷和變形量依次減小。

表 2 數值解析解與有限元解對比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 滾動體接觸應力

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

4.2.4 偏載接觸模型結果對比

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

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

4.2.5 不同載荷分布對偏載接觸效應的影響

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

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

5 結論

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

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

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

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

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