基于擬靜力學的軸承內圈滾道工作應力場研究

中圖分類號： TH133.33+1 DOI： 10.16578/j.issn.1004.2539.2025.05.016

0 引言

軸承滾道工作應力場是指軸承在受載運轉過程中，滾動體與滾道之間相互接觸所引起的、在軸承滾道表層所形成的工作應力場，由軸承的工況條件決定，包括軸向力、徑向力和傾覆力矩等。研究軸承滾道工作應力場，對提高滾動軸承的服役壽命與可靠性有重要意義。

為求解軸承滾道工作應力場，需要得到特定工況下滾動體與滾道之間的接觸角、接觸載荷和接觸區等參數。在求解這些參數的諸多模型中，擬靜力學模型的應用最為廣泛[1]。陳瑛琳等[2]12-16提出一種適用于球數較多的角接觸球軸承的擬靜力學簡化計算方法，得到了低速重載工況下滾動體與內圈滾道接觸角分布。王延忠等[3]I3-17在擬靜力學模型基礎上，引入對迭代初值的幾何約束，通過選取合適的初值，解決了因初值選取不當導致收斂困難的難題。王明凱等4建立5自由度擬靜力學模型，在此基礎上，計算得到高速角接觸球軸承在聯合載荷作用下滾動體與內外圈滾道的接觸角、接觸載荷等力學參數，并進一步探究了軸向力、徑向力等工況條件改變對軸承力學參數的影響。

軸承滾道工作應力場的分布狀態對軸承疲勞壽命有著十分重要的影響，最大切應力峰值及其所在深度是決定軸承疲勞壽命的關鍵因素。一般認為，疲勞裂紋萌生的特征深度位于最大切應力處5，且最大切應力峰值越大，萌生疲勞裂紋的概率就越大；峰值所在深度距滾道表面越近，疲勞裂紋擴展至滾道表面所需的時間就越短，導致軸承的疲勞壽命也就越短。因此，開展軸承滾道工作應力場研究，尤其是最大切應力分布的研究，是進行滾動軸承疲勞壽命分析的重要基礎。

本文以7008C軸承內圈為研究對象，通過分析滾動體與滾道之間的接觸與變形，建立并求解軸承擬靜力學方程組，得到滾動體與內圈滾道之間的接觸應力分布；在此基礎上，利用有限元軟件AnsysWorkbench建立7008C軸承內圈滾道工作應力場有限元模型，分析了工況條件對最大切應力分布的影響。

1滾動體與滾道接觸分析

高速精密角接觸球軸承7008C通常被用作機床主軸軸承，其幾何結構如圖1所示，具體的幾何參數如表1所示。

根據赫茲接觸理論，滾動體與內外圈滾道橢圓接觸區的長、短半軸長度 a 、b，以及彈性趨近量δ的

數值計算式[7]70-71分別為

式中， Q 為滾動體與內外圈滾道之間的法向接觸力； 分別為滾動體與滾道的彈性模量，對于全鋼軸承，均為 2.07×10⁵MPa^[7]73-74 分別為滾動體與滾道的泊松比，對于全鋼軸承，均為 0.3^[7]73-74 m_a，m_b 均為與橢圓偏心率 e 相關的系數； k 為接觸橢圓短半軸與長半軸的比。 （20 k 的表達式分別為

式（3）中的 K（e） 為第一類完全橢圓積分，表達式為

式中， φ 為角度變量。

式（4）＼～式（5）中的 L（e） 為第二類完全橢圓積分，表達式為

式（7）＼～式（8）中， e 的表達式為

式（1）＼～式（3）中， Σρ 為滾動體與內外圈滾道接觸的主曲率和，表達式為

Σρ=ρ_I1+ρ_I2+ρ_II1+ρ_II2

式中， 表示滾動體與內外圈滾道接觸時的主曲率。角接觸球軸承7008C的主曲率計算式[8]12-13與結果如表2所示。

表2滾動體與內外圈滾道接觸時主曲率的計算式與結果

表2中， γ 為引入的一個無量綱參數，表達式為

由于橢圓偏心率 e 未知，目前還無法直接計算m_a 、 m_b 和 K（e） 。為了計算它們，需定義滾動體與內外圈滾道接觸時的主曲率函數 F（ρ）^[8]12-13 ，表達式為

利用式（10）＼～式（12）和表2，可以直接計算出角接觸球軸承7008C的主曲率函數 F（ρ） 的值。根據赫茲接觸理論[8]1-13，主曲率函數 F（ρ） 還可以表示為

至此，可利用迭代法數值求解 和 K（e） ，基本過程為： ① 利用式（12）計算主曲率函數 F（ρ） 的準確值； ② 給定 k 的初值，由式（6）可知， k 在0＼～1變化，本文中設置 k 的初值為0.01； ③ 利用式（13）獲得 F（ρ） 的計算值； ④ 計算主曲率函數 F（ρ） 的準確值與計算值之間的差值，若差值大于 10^-5 ，則 k 值增加 10^-5 ，再重新計算差值，如此迭代計算，直到差值小于 10^-5 ，最終得到 k 的終值； ⑤ 基于 k 的終值，分別利用式（4）、式（5）和式（7）計算 m_a 、 m_b 和 K（e） 。

根據赫茲接觸理論，滾動體與內外圈滾道橢圓接觸區內的壓力分布為橢球狀不均勻分布，壓力值隨坐標值不斷變化并在中心點處取得最大值 P₀ ，如圖2所示。

滾動體與內外圈滾道橢圓接觸區內的壓力分布計算式可表示為

橢圓接觸區中心點處的最大壓力為

為計算 、 δ 和 P₀ ，還需建立軸承的擬靜力學方程組，以求解滾動體與內外圈滾道之間的法向接觸力 Q 。這將在下一節詳細描述。

2 擬靜力學方程組

2.1擬靜力學方程組的建立

首先，對滾動體依據其方位角進行編號，如圖3所示。圖中， Δψ 為任意相鄰兩滾動體之間的方位角之差，且 Δψ=360^°/z ; ψ_j 為第 j 個滾動體的方位角。

在軸承受載之前，滾動體球心與內外圈滾道的曲率中心共線。在軸承受載后，在軸向力和徑向力的作用下，內圈滾道的曲率中心將發生軸向與徑向位移。在離心力和陀螺力矩的作用下，滾動體球心將發生偏移。由于外圈固定在支座上，故其曲率中心不發生偏移。軸承受載前后，滾動體球心與內外圈滾道曲率中心之間的位置關系940-41如圖4所示。

圖4中， O_e 為外圈滾道曲率中心； o 、 分別為滾動體球心初始位置和末位置； O_i 2 O_i^′ 分別為內圈滾道曲率中心初始位置和末位置； α_ij 人 α_ej 分別為軸承受載后滾動體與內外圈所成的接觸角； α₀ 為受載前滾動體與內外圈所成接觸角，即初始接觸角； X_1j X_2j 分別為外圈滾道曲率中心與滾動體球心末位置的軸向、徑向距離； A_1j 1 A_2j 分別為外圈滾道曲率中心與內圈滾道曲率中心末位置的軸向、徑向距離； l_1j， l_2j 分別為軸承受載前后內圈滾道曲率中心的軸向、徑向位移； l_ij?l_ej 分別為軸承受載之后內、外圈滾道曲率中心到滾動體球心距離； l_j 為在任意方位角滾動體j 處，軸承受載之前外圈滾道曲率中心與內圈滾道曲率中心之間連線距離。

根據圖4中的幾何關系，受載前后軸承的變形協調方程組[9]40-41為

其中， A_1j?A_2j 表達式分別為

式中， δ_a 人 δ_? 和 θ 分別為軸承內圈的軸向、徑向和角位移； R_i 為內圈滾道曲率中心軌跡圓半徑，其表達式為

l_j 的表達式為

式（16）中， l_ij?l_ej 的表達式分別為

式中， δ_?ij?δ_?ej 分別為滾動體 j 對內、外圈的彈性趨近量。

在軸承受載運行過程中，所有滾動體對內圈的作用力應與外部載荷平衡。由此可得，軸承內圈受力平衡方程組[9]40-41為

式中， F_a， F_r 和 M 分別為軸承所受的軸向力、徑向力和傾覆力矩； Q_ij 為滾動體 j 與內圈的法向接觸力。

任意方位角處滾動體 j 的受力分析如圖5所示。圖5中， F_cj 為滾動體受到的離心力； M_gj 為滾動體受到的陀螺力矩； ω_?mj，ω_?Rj 分別為滾動體的公轉、自轉角速度； ω_i 為軸承內圈運轉角速度； β_j 為第 j 個滾動體的姿態角； Q_ij 、 Q_ej 分別為滾動體與內、外圈的法向接觸力。

根據滾動體軸向與徑向受力平衡，滾動體受力平衡方程組[9]40-41為

式中， λ_ij. λ_ej 分別為內外滾道控制系數，本文采用外滾道控制，取 λ_ij=0 ， λ_ej=2^[10]93-94 ; J 為轉動慣量； ?_m 為滾動體的質量； Q_ij 、 Q_ej 分別為滾動體對內、外圈的彈性趨近量，即

式中， K_ij 、 K_ej 分別為滾動體與內、外圈滾道的接觸載荷變形系數，具體計算方法見文獻[11]；@mi、@Ri以及 β_j 的計算方法見文獻 [10]^63-71

至此，聯立軸承變形協調方程組式（16）、軸承內圈受力平衡方程組式（21）和滾動體受力平衡方程組式（22），即可建立高速精密角接觸球軸承7008C的擬靜力學方程組。

2.2擬靜力學方程組的數值求解方法

由上述擬靜力學方程組可知，對于第 j 個滾動體而言，共有4個未知數，分別為 δ_ij δ_ej 、 X_1j 、 X_2j ；再加上軸承內圈的軸向位移 δ_a 、徑向位移 δ_r 和角位移 θ 對于整個軸承而言，一共有 4z+3 個未知數。為進一步化簡計算流程，由文獻[2]13-14可進一步得到以下化簡關系

通過上述關系式，將關于 δ_ij， δ_?ej ， X_1j， X_2j， δ_a ， δ_r 和 θ 的 4z+3 個未知量轉化為關于 α_ij?α_ej 的 2z+3 個未知量，使未知數的個數大大減少，可以使計算更加簡便且結果更容易收斂。

采用Newton-Raphson方法對建立的擬靜力學方程組進行求解，并引入對迭代初值的幾何約束[3]16-17，以提高數值計算的收斂速度。計算流程如圖6所示，圖6中，收斂精度 ε 取值為 10^-8 。

2.3 模型驗證

為驗證本文數值求解算法的準確性，在某一軸承參數下，將本文的計算結果與文獻 [12]^26-27 、[13]427-43進行對比，結果如表3所示。由表3可知，本文計算結果的最大誤差小于 3% ，這表明本文的數值求解算法具有較高的可靠性。具體的軸承參數為： z= 19； f_i=f_e=0.52 ； ； d_m=105mm ； α₀= 25^° . F_a=889.84N ；軸承工作轉速 n=15000r/min 。

在 F_a=500N 、 F_r=500N 、 M=5N?m 以及工作轉速 的工況條件下，滾動體對內圈接觸角 α_ij 以及與內圈接觸區內最大壓力 P₀ 的計算結果如圖7所示。由圖7可知，方位角為 0^° 處的1號滾動體受載最大。后續以其與內圈滾道的接觸為對象，分析軸承內圈滾道工作應力場，求解其與內圈滾道之間的橢圓接觸區、接觸載荷和接觸角，并作為工作應力場分析的必要條件。

3軸承內圈滾道工作應力場分析

3.1 有限元模型建立

采用AnsysWorkbench軟件建立7008C軸承內圈的幾何模型，并建立圖8所示方位角為 0^° 處滾動體與內圈滾道之間的橢圓接觸區，其中，長、短半軸a 、 b 和接觸角 α_i 分別為 0.55mm 、 0.087mm 和20.2^° 。在劃分網格時，對橢圓接觸區以 0.01mm 為單元尺寸進行網格細化，將其余區域的單元尺寸設置為 0.5mm ，網格增長率設置為1.1，并以六面體網格為主導進行劃分，如圖8所示。此時，共有180138個網格，653800個節點。網格獨立性驗證表明，進一步加密網格對仿真結果的影響不到 1.5% 。

軸承內圈材料為淬硬軸承鋼GCr15，工作應力場分析需要的力學性能參數：密度為 7830kg/m^3[14]2-8 彈性模量為 2.07×10⁵MPa^[7]73-74 ，泊松比為 0.3^[7]73-74 屈服強度為 1410.17MPa^[14]2-8 。在施加邊界條件時，考慮軸承實際工況下的安裝條件，限制軸承內圈兩端面的軸向位移，并限制內孔表面的周向與徑向位移，如圖9所示；然后將算得的壓力分布施加于橢圓接觸區。圖7中方位角為 0^° 處滾動體與內圈滾道之間算得的壓力分布如圖10所示。

3.2 工作應力場分布特征分析

在 F_a=500N. F=500N. ， M=5N?m 及 的工況條件下，受載最大滾動體與內圈滾道接觸區附近的工作應力場分布如圖11所示。由圖11（a）＼～圖11（c）可知，周向、軸向和徑向應力均為壓應力，且最大值分布在滾道表面，隨深度增加而逐漸減少；徑向應力最大，軸向應力次之，周向應力最小；接觸區以外應力數值較小，且為拉應力。由圖11（d）可知，最大切應力隨深度呈現先增大后減小的趨勢，在距滾道表面一定深度處取得最大值。

3.3工況條件對最大切應力分布的影響

軸承內圈滾道工作應力場最大切應力分布隨軸向力、徑向力和工作轉速的變化如圖12所示。由圖12（a）和圖12（b）可知，最大切應力峰值隨軸向力及徑向力的增加均呈現增大的趨勢，峰值所在深度均隨軸向力及徑向力的增加緩慢向滾道表面移動。在增加相同作用力的前提下，徑向力對最大切應力峰值的影響大于軸向力。

由圖12（c）可知，軸承內圈滾道工作應力場中，最大切應力峰值隨工作轉速的增加呈現減小的趨勢，并且峰值所在深度逐漸遠離滾道表面。這是因為隨著工作轉速的增大，滾動體的離心力也增大，滾動體與內圈滾道之間的法向接觸力減小。

4結論

1）通過分析滾動體與滾道之間的接觸與變形，建立并求解軸承擬靜力學方程組，得到了滾動體與內圈滾道之間的接觸應力分布。研究表明，方位角為 0^° 處的1號滾動體受載最大。

2）利用有限元軟件AnsysWorkbench仿真計算了滾動體與內圈滾道接觸區工作應力場分布。結果表明，滾動體與滾道接觸區附近的周向、軸向和徑向工作應力均為壓應力，且最大值分布在滾道表面，隨深度增加而逐漸減小；最大切應力隨深度增加呈現先增大后減小的趨勢。

3）隨著軸向力和徑向力的增加，最大切應力峰值呈現增大的趨勢，峰值所在深度逐漸靠近滾道表面。隨著工作轉速的增大，滾動體的離心力增大，滾動體與內圈滾道之間的法向接觸力減小，從而導致最大切應力峰值減小，峰值所在深度逐漸遠離滾道表面。

參考文獻

[1］吳偉開.基于擬靜力學的滾動軸承動態性能分析技術與系統研 究[D].無錫：江南大學，2023：3-4. WU Weikai. Research on dynamic performance analysis technology and system of rolling bearing based on pseudo-statics[D]. Wuxi：Jiangnan University，2023：3-4.

[2]陳瑛琳，史文華，陳江義，等.RV減速器主軸承擬靜力學模型簡 化分析方法[J].軸承，2021（4）：12-17. CHENYinglin，SHIWenhua，CHENJiangyi，etal.Simplified analytical method for quasi-static model of main bearing forRV reducer[J].Bearing，2021（4）：12-17.

[3]王延忠，鄂世元，賈彥蓉，等.聯合載荷下改進的角接觸球軸承 擬靜力學分析模型[J].軸承，2022（11）：13-17. WANG Yanzhong，E Shiyuan，JIA Yanrong，etal.Improvedquasistatic analysis model of angular contact ball bearing under combined load[J].Bearing，2022（11）：13-17.

[4]王明凱，閆柯，馬帥軍，等.聯合載荷作用下精密角接觸球軸承 動態特性演變規律[J].航空動力學報，2022，37（6）：1134-1149. WANG Mingkai，YAN Ke，MA Shuaijun，et al. Evolutionary regularityof dynamic characteristics for precision angular contact ball bearing under combined loads[J]. Journal of Aerospace Power， 2022，37（6）：1134-1149.

[5]SANTUS C，BEGHINI M，BARTILOTTAI，et al. Surface and subsurface rolling contact fatigue characteristic depths and proposal of stressindexes[J].International Journal of Fatigue，2O12，45： 71-81.

[6]SADEGHI F，JALALAHMADI B，SLACK T S，et al. A review of rolling contact fatigue[J].Journal of Tribology，20o9，131（4）： 041403.

[7］周建男，陳長征，周劬勛.軋鋼機械滾動軸承[M].北京：冶金工 業出版社，2001：70-85. ZHOU Jiannan，CHEN Changzheng，ZHOU Quxun.Rolling bearings for steel rolling machinery[M].Beijing：Metallurgical Industry Press，2001：70-85.

[8］梁曉.高速角接觸陶瓷球軸承滾道接觸應力場分析[D].濟南： 山東大學，2015：12-30. LIANG Xiao.Contact stress field analysis of raceway of highspeed angular ceramic ballbearings[D]. Jinan： Shandong University，2015：12-30.

[9]萬超，張亢，田澤宇，等.高轉速下角接觸球軸承接觸應力分布 特征與疲勞剝落關系研究[J].機械傳動，2021，45（8）：38-44. WAN Chao，ZHANG Kang，TIAN Zeyu，et al. Study on the relationship between contact stress distribution characteristic and fatiguespallingof angular contact ball bearingat high speed[J]. Journal ofMechanical Transmission，2021，45（8）：38-44.

[10]HARRISTA，KOTZALASMN.Advanced conceptsof bearing

technology：rollingbearinganalysis[M].5th ed.Boca Raton：CRC Press，2006：63-94.

[11]汪智慧.航空滾動軸承力學特性分析研究[D].南京：南京航空 航天大學，2010：13-16. WANG Zhihui.Analytical research on mechanical propertiesof aerial rolling bearings[D].Nanjing：NanjingUniversity ofAeronautics and Astronautics，2010：13-16.

[12]徐立暉.高速電主軸角接觸球軸承力學特性研究[D].杭州：浙 江大學，2015：26-27. XULihui.Research on mechanical properties of angular contact ball bearings in high-speed motorized spindle[D].Hangzhou：ZhejiangUniversity，2015：26-27.

[13]朱益利，金超武，許磊，等.滾珠軸承力學模型的數值求解方法 研究[J].中國機械工程，2013，24（4）：427-431. ZHUYili，JINChaowu，XULei，etal.Numericalsolutionmethods of ball bearing mechanics model[J].China Mechanical Engineering，2013，24（4）：427-431.

[14]GUOYB，LIU CR.Mechanical properties of hardened AISI 52100 steel in hard machiningprocesses[J].Journal ofManufacturingScienceandEngineering，2002，124（1）：1-9.

Research on working stress field ofbearing inner ring raceway based on quasi-statics

GUOHaoyu'WANGDexiang12 GUO Feng12LI Xinming1，2 （1.SchoolofMechanicalandAutomotiveEngineering，Qingdao UniversityofTechnology，Qingdao26652o，China） （2.Shandong Zhilian Community Bearing Technology Co.，Ltd.，Liaocheng 252664，China）

Abstract：[Objective]Theresearchoftheworking stressfieldofthebearingracewayisanimportantbasis forthe fatigue lifeanalysisofrollingbearings.[Methods]Thecontacttressdistributionbetweentherollngelementandtheinnerringraceway of7008Cbearingwasobtainedbyanalyzing thecontactanddeformationbetween therollngelementandtheraceway，andby establishingandsolvingthebearingquasi-staticsequations.Onthisbasis，thefiniteelementmodeloftheworkingstresfieldof theracewaywas establishedusing theAnsys Workbench.Theinfluenceofworkingconditionsonthe maximumshearstresdistributionWas analyzed.[Results] The results show that the No.1 rolling element with anazimuth of 0^° has the maximum load. Withtheincreaseofaxialforceandradialfore，thepeakvalueofthemaximumshearstresshowsanincreasing trend，andthe depthofthepeakvaluegraduallyapproachestheracewaysurface.Withtheincreaseoftheworkingspeed，thepeakvaluedecreases，and the depth gradually moves away from the raceway surface.

Keywords：Bearing;Quasi-statics;Bearingraceway;Working stress field

（上接第92頁）

Imbalance response analysis of permanent magnet eddy current couplings based on the transfer matrix method

WANG JianXU Chengxi （School ofAutomation，Nanjing InstituteofTechnology，Nanjing211167，China）

Abstract：[Objective] BasedontheRiccati transfer matrix method，thedynamic lumped model fortheanalysisof the imbalancevibrationcharacteristicsofpermanent magneteddycurrentcouplings was proposed.[Methods]Bycalculatingand analyzingtherotorofthecouplingwithmultiplepresetimbalances，thevibrationcharacteristicswereexploredunderdifferent imbalancedistributions.Thevibrationcharacteristicsweresimulatedbyusingthre-dimensionalfiniteelementmethods，andthe analysisresultswerecompared with theresultsobtained fromthe transfer matrix method.[Results]Theresults indicate that withinthetestspeedrange，comparingtheimbalanceefectsbetweenthecasewheretheimbalanceaffcts therotoronthesame sideoftheobservationandthatontheoppositesideoftheobservation，therotoronthesamesidegenerateslargerimbalance vibration.

KeyWords：Permanent magnet eddy current coupling; Analysis of imbalance;Transfer matrix method