外圈故障引起圓柱滾子軸承IAS擾動動力學建模

摘要： 針對瞬時角速度（instantaneous angular speed，IAS）信號與軸承機械動力學之間的關系研究和工業機器人RV減速器轉臂軸承在低速工況下易發生失效的問題，提出了一種三自由度局部故障滾子軸承IAS擾動的動力學模型。模型基于Hertz線接觸理論，分析了滾子與滾道間相互作用產生的局部變形對IAS的影響機理，給出了耦合切向力與力矩的計算方法；將故障沖擊引入扭矩分析，計算故障區帶來的扭矩變化，增加角自由度，建立了法向力與切向力相耦合的三自由度軸承IAS擾動動力學模型。采用四階?龍庫塔數值積分法求解，并在近似條件下將仿真與實驗結果進行對比分析。結果表明，該模型能較好地解釋圓柱滾子軸承IAS擾動的來源，有效反映外圈故障對IAS帶來的影響，進一步完善了滾動軸承動力學理論。

關鍵詞： 轉子動力學； 圓柱滾子軸承； 瞬時角速度； Hertz線接觸

中圖分類號： O347.6；TH133.33""" 文獻標志碼： A""" 文章編號： 1004-4523（2025）02-0441-08

DOI：10.16385/j.cnki.issn.1004-4523.2025.02.023

收稿日期： 2023-03-17； 修訂日期： 2023-06-14

基金項目："國家自然科學基金資助項目（52165067）；云南省重大專項科技計劃資助項目（202002AC080001）

Dynamic modeling of instantaneous angular speed disturbances in cylindrical roller bearings caused by outer ring failure

HAN Qiao， GUO Yu， FAN Jiawei

（Faculty of Mechanical and Electrical Engineering，Kunming University of Science and Technology，Kunming 650500， China）

Abstract： In light of investigating the interplay between the instantaneous angular speed （IAS） signal and mechanical dynamics of bearings， alongside addressing the vulnerability of rotary arm bearings within industrial robot RV reducers to failure under low-speed conditions， this paper presents a novel three-degree-of-freedom dynamics model to explain the IAS perturbations resultant from localized roller bearing failures. The model， rooted in the Hertz line contact theory， dissects the influence mechanism of localized deformations stemming from roller-raceway interactions on the IAS. A comprehensive approach to computing the coupled tangential force and torque is outlined. The integration of failure-induced impacts into torque analysis computes torque variations introduced by the failure zone， thereby augmenting angular degrees of freedom. As a result， a three-degree-of-freedom bearing IAS disturbance dynamic model， coupling both normal and tangential forces， is established. Employing fourth-order Runge-Kutta numerical integration， the model is solved， and its results are meticulously compared and analyzed against experimental approximations under near-approximate conditions. The findings underscore the model’s ability to effectively expound upon the origins of IAS perturbations in cylindrical roller bearings while adeptly reflecting the ramifications of outer ring failures on IAS. This work contributes to the refinement of rolling bearing dynamics theory， advancing the comprehension of intricate mechanical systems.

Keywords： rotor dynamics；cylindrical roller bearings；instantaneous angular speed；Hertz line contact

RV減速器作為工業機器人的核心部件之一，結構緊湊，包含眾多旋轉零部件，也是工業機器人主要的故障源。軸承在RV減速器中起支承及傳遞力和扭矩的作用［1］，受力狀況復雜，RV減速器的失效大多由軸承失效引起。其中圓柱滾子軸承作為RV減速器的轉臂軸承，安裝在曲柄軸中間，外圈與擺線輪相連，用于支撐曲柄擺線輪，在低速重載的條件下，往往最先發生失效［2］。然而，工業機器人低速工況下圓柱滾子軸承傳遞能量低、負載高，傳統的低速振動響應分析通常無法檢測到故障發出的振動。瞬時角速度（instantaneous angular speed，IAS）與軸的旋轉嚴格同步，更少地依賴于故障區與傳感器之間的傳輸路徑，可以有效揭示滾子元件進入軸承圈座剝落區的動態行為，彌補低速故障軸承監測的不足。

目前，基于編碼器IAS信號的旋轉機械故障診斷研究逐漸引起國內外專家學者的關注，該技術在工業機器人RV減速器等應用場合較振動信號有著更大的優勢［3］，其不受振動信號面臨的時變傳遞路徑、界面間非線性能量耗散和多源多路徑振動信號疊加衰減等因素的影響，直接與旋轉機械動態特性相關，擁有更高的信噪比，便于故障特征的識別、提取和檢測。例如，MOUSTAFA等［4］介紹了一種基于IAS的技術，用于低速滾動軸承故障檢測和故障寬度估計。AIT SGHIR等［5］提出一個基于譜相關密度的指標，證明了IAS信號在低速狀態下診斷滾動軸承的適用性。為了方便在齒輪傳動和軸承診斷中了解缺陷的真正角度周期性，RéMOND等［6］提出了一些處理角域經典運動方程的建模方法，用于描述旋轉機器在非恒速運行條件下的行為，如風力渦輪機或直升機等。在考慮滾動元件動力學的研究中，SAWALHI等［7］建立了一種考慮滾動元件動力學特性的模型，該模型利用Hertz理論估計滾動單元所產生的非線性接觸力，并通過與其位移相關的時變剛度進行轉化，同時將局部缺陷添加到內圈（IR）、外圈（OR）和滾動元件中。BOURDON等［8］開發了一種在角域描述機械系統的運動方程，介紹了代表運動部件機械故障的角周期扭矩擾動的參數模型。其所針對的機械結構由兩個非線性方程組描述，可以用常規數值方法求解；GOMEZ等［9］對SAWALHI建立的模型進行了擴展，采用角度法建立了相應的有限元模型，其中設置了六個自由度，允許模擬非平穩條件及IAS現象。THIBAULT等［10］建立了高速軸向載荷下球軸承的動力學模型，對IAS的性能與局限性進行了分析。

雖然上述學者對滾動軸承開展了部分動力學建模研究，但大多數研究仍針對振動響應，針對IAS響應進行動力學建模的研究較少，對解釋軸承IAS擾動的原因以及故障區帶來的附加扭矩對IAS影響的研究不足，建立的動力學模型較為簡單，較少關注到IAS對于低速或超低速工況下軸承動態行為監測的有效性［4］。雖然GOMEZ等［9］通過有限元模型在滾動軸承動力學中引入了角度擾動，但是僅針對中高速工況，對故障區位移變化計算較為理想，且忽略了滾子進出故障區時帶來的附加扭矩，未能通過動力學模型對低速或超低速工況下故障軸承IAS擾動進行解釋。因此，對低速或超低速工況下滾子軸承IAS擾動現象進行深入研究并建立相應的局部故障滾子軸承IAS動力學模型具有重要意義。

滾子軸承的局部故障一般會導致異常的接觸力激勵和位移激勵，這些激勵會引起軸承瞬時扭矩的變化，從而引起IAS變化。據此，本文以低速工況下含外圈局部故障的NU204 ECP型軸承為研究對象，基于Hertz線接觸理論，增加內圈角自由度，分析無故障滾子軸承IAS擾動的原因，計算法向接觸力與切向力的耦合，考慮滾子進出故障區時帶來的附加扭矩影響，建立三自由度局部故障圓柱滾子軸承IAS擾動的動力學模型。通過仿真和實驗對比分析，驗證了所建模型的正確性，討論了故障對IAS信號的影響，進而為低速工況下更高精度的滾子軸承狀態監測提供理論依據。

1 無故障軸承受力分析

無故障滾子軸承在徑向載荷作用下，位于載荷區的滾子與滾道間相互作用，局部變形導致出現滾動阻力現象，使得切向力與接觸力相耦合，由此引起扭矩的變化，此即為無故障圓柱滾子軸承IAS擾動的起源。

1.1 滾子承受徑向力計算

通常情況下，位于徑向載荷作用線上的滾子所承受的載荷最大，其兩邊滾子所承受的載荷依次遞減，且滾子的載荷分布相對于載荷線對稱，如圖1所示。

由靜態平衡條件可知，作用在滾子上的徑向外載荷等于滾子所承受的載荷的垂直分量。滾子所受載荷的大小以及承受載荷的滾子數由軸承內部的幾何參數（游隙等）和外部載荷的大小共同決定。

載荷區內滾子在任意時刻所受載荷的大小可表示為［11］：

（1）

式中，表示載荷區圓心角度的一半；ε為載荷分布系數，它取決于承載區范圍大小；n為載荷?形變指數，對于滾子軸承n=10/9；Qmax表示滾子所受最大載荷數值，STRIBECK［12］提出其可表示為：

（2）

式中，Z為滾子個數；Fy為徑向負載大小；α為軸承接觸角，本研究中接觸角為0°；第k個滾子在任意時刻的角位置表示為：

（3）

式中，Φ0為第一個滾動體初始時刻的角位置；ωcage為保持架角速度，計算公式為：

（4）

式中，ω為轉軸角速度；Dc為滾子直徑；Dp為節圓直徑。

1.2 局部變形引起的耦合力計算

不同于球軸承滾子與滾道的點接觸問題，圓柱滾子軸承中滾子與滾道的接觸屬于線接觸，假設相接觸的兩個圓柱體光滑且長度均勻，軸線相平行，可表示為圖2（a）所示的兩彈性圓柱體的接觸。由于載荷的作用，兩圓柱體在接觸點將發生彈性變形，形成一個具有一定寬度的矩形接觸區，線接觸彈性理論研究的即為兩個彈性圓柱體受壓接觸后產生的局部應力和應變沿母線的分布規律。滾子?滾道Hertz線接觸模型如圖2（b）所示。

滾子與滾道的接觸面半寬為：

（5）

式中，Q為載荷；l為有效接觸長度；為兩接觸體材料的泊松比；r1、r2為接觸體曲率半徑，若滾子與滾道為外接觸（即兩凸面），則r1與r2之間用正號，反之（一凸面與一凹面接觸）為負號［13］；E為等效彈性模量，通過下式計算：

（6）

式中，E1、E2為兩接觸體材料的彈性模量。

外接觸的彈性變形量為：

（7）

內接觸的彈性變形量為：

（8）

滾子與內滾道的接觸剛度為：

（9）

滾子與外滾道的接觸剛度為：

（10）

在徑向載荷的加載下，彈性接觸總變形量可表示為：

（11）

圓柱滾子軸承總接觸剛度表示為：

（12）

進而基于Hertz線接觸理論的法向接觸力N估算為：

（13）

通過法向接觸力施加點的向前移動來模擬局部變形，產生滾動阻力現象進而出現切向力的耦合。如圖3所示，wc為質心角速度，V為質心線速度，前移距離bn用滾動阻力系數μ與滾子半徑rc的函數表示：

（14）

其中，μ直接關系到IAS擾動模型，其選擇應為滾子軸承恒定摩擦系數［9］，本研究中取為0.0015。

對滾子列出力與力矩平衡方程：

（15）

式中，mc為滾子質量；Ic為滾子等效慣量；θk為第k個滾子在xoy平面內的絕對旋轉角度；Tik，Tek分別為滾子與軸承內、外圈接觸產生的切向力；Nik，Nek分別為滾子與軸承內、外圈接觸產生的法向力，數值與N相等。其中，切向力為代數量，方向與滾動元件質心位移方向相同，所有分析忽略相互作用體的內部阻尼，內圈受力如圖4所示。

通過解上述方程，得到切向力關于法向接觸力的計算公式：

（16）

法向接觸力作用點的前移使得其作用線偏移軸心，進而對軸產生扭矩且方向與電機扭矩相反，切向力與軸相切，力臂即為內圈滾道半徑，扭矩方向與法向力矩相同。第k個滾子任意角位置對軸產生的扭矩變化用下式計算：

（17）

式中，ri為軸承內圈滾道半徑。

2 外圈故障引起IAS擾動建模

2.1 外圈故障引起的時變位移激勵模型

圓柱滾子軸承內圈旋轉，外圈固定在軸承座上，所以外圈剝落位置也固定。當外滾道局部剝落類型為故障長度大于滾子長度時，滾子落入故障區進而產生附加位移。隨著滾子經過外圈故障區，滾子與內外滾道的相對位置發生改變，因此需要構建相應的數學模型準確描述滾子與外圈剝落之間接觸位移的變化，來彌補滾子與正常滾道接觸中傳統接觸位移計算方法的不足。

從滾子進入故障區開始，相對位移量逐漸增大，在滾子運動觸碰到故障區中心位置時，相對內圈的位移量達到最大；而后，滾子在保持架的作用下沿滾道繼續公轉，從故障區正中心位置逐步退出故障區，滾子在內外圈之間的位移從最大值逐漸恢復到正常接觸狀態。基于上述滾子故障區運動分析，在滾子與無故障內圈和故障外圈接觸的整個過程中，建立以下時變位移激勵模型：

（18）

式中，Cd表示滾子進入故障區后的時變位移增量；表示故障區角度位置；mod（·）表示求余函數算子；Hm表示滾子撞擊故障區后緣那一刻的位移增量，計算式為：

（19）

設Θos表示故障區所在圓周上對應的圓心角弧度的一半，則：

（20）

式中，Ld為故障區沿軌道方向對應寬度；Do為軸承外圈滾道內徑。

2.2 瞬時沖擊力計算

當軸承外圈發生疲勞剝落時，滾子滾過剝落區接觸形式突變，徑向力將進行再分配，撞擊剝落區后緣的一瞬間產生沖擊力，從而影響扭矩的變化，誘發軸承IAS的明顯異常變化。

沖擊是動能向系統的傳遞，在相對較短的時間內發生。假設系統是保守的，則狀態1（進入故障區前一刻，如圖5所示）和狀態2（滾子撞擊故障區后緣時，如圖6所示）之間的機械能守恒為：

（21）

式中，g為重力加速度；λ表示滾子撞擊故障區后緣那一刻相對滾子質心的角度；V1、V2分別為狀態1與狀態2時滾子質心線速度；ω1、ω2分別為狀態1與狀態2時滾子質心角速度。剛進入故障區的圓柱滾子的慣量I1為：

（22）

由平行軸定理可以得到I2為：

（23）

質心的線速度V和角速度wc的關系為：

（24）

整理可得：

（25）

在進入剝落區到撞擊后緣時使用動量守恒定理可得：

（26）

式中，Vos表示滾子撞擊后緣時滾子的撞擊速度：

（27）

且可以知道撞擊后緣的時間為：

（28）

式中，Vout表示故障外圈前緣處的線速度，由于外圈固定不動，Vout=0。由式（26）～（28）可得滾子撞擊剝落區后緣的瞬時沖擊力為：

（29）

2.3 故障區帶來的時變附加力矩計算

由上述可知，當滾子進入故障區時起，滾子與無故障內圈和故障外圈的接觸變形逐漸變小，因此，滾子與內外圈之間的接觸力也逐漸變小，并在滾子與故障區后緣碰撞的前一刻減到最小值；滾子在保持架的帶動下撞擊故障區后緣的一瞬間產生瞬時沖擊力，接觸力迅速從最小值變到最大值；而后，滾子在保持架的作用下沿滾道繼續公轉，從故障區正中心位置逐步退出故障區的這一過程，接觸力也就從最大值逐步恢復到其正常的接觸力狀況。綜上所述，給出滾子經過故障區的時變接觸力模型：

（30）

因此，滾子經過故障區的瞬時沖擊力對內圈帶來的時變附加力矩為：

（31）

其中，沖擊力矩MT方向與電機施加的外扭矩方向相反。

2.4 外圈故障引起IAS擾動的動力學模型

為分析滾子軸承健康及故障引起的IAS擾動機理，需建立相應的動力學模型。將滾子與內、外圈間的接觸簡化為彈簧?質量系統［14］，如圖7所示，且滿足以下假設：（1）內圈、外圈與滾子之間沒有相對滑動；（2）假設軸承元件的接觸滿足Hertz接觸條件；（3）僅計算滾子與滾道間的彈性接觸變形，忽略軸承套圈的整體變形，將滾子視為非線性彈性接觸元件；（4）忽略保持架動力學、缺陷形狀和彈性流體動力潤滑（EHL）的非線性影響。

基于剛性套圈假設［15］，經以上簡化，可以得到滾子軸承三自由度的IAS擾動動力學方程為：

（32）

式中，m為軸承內圈和轉軸的等效總質量；I為軸承內圈和轉軸的等效總慣量；Cs為系統的阻尼；CT為內圈與空氣間的等效黏性阻尼；x、y、θ分別為內圈在水平和豎直方向的振動位移以及角位移；、分別為內圈在水平和豎直方向的振動速度；為內圈角速度；和分別為內圈在水平和豎直方向的振動加速度；為內圈角加速度；Qx和Qy分別為內圈在水平和豎直方向上所受外加徑向力；Mz為電機帶來的外加扭矩；δod表示每個滾子的時變接觸變形，計算式為：

（33）

式中，c為軸承初始徑向游隙；γ為用于判斷滾子是否進入故障區的參數：

（34）

η為用于判斷滾子是否在載荷區與內外滾道發生接觸變形的參數：

（35）

3 動力學仿真與實驗分析

3.1 動力學仿真

仿真研究以圓柱滾子軸承NU204 ECP為研究對象，其幾何尺寸參數如表1所列。設m=1.3 kg，系統等效阻尼Cs=400 N·s/m，扭轉阻尼系數與軸轉速有關［9］，且CT=1.52 N·m·s/rad，軸承受到沿x軸和y軸方向的外部徑向作用力分別為Fx=0 N，Fy=400 N，電機扭矩Mz=9.54 N·m；假設故障區沿外圈滾道方向寬度Ld=1 mm，故障長度大于滾子長度且位于外圈滾道正下方，軸承內圈在水平和豎直方向的初始位移分別為x=10-9 m和y=10-9 m，初始速度，初始角位移θ=10-9 rad，初始角速度；時間步長Δt=10-6 s，仿真時長t=10"s，內圈轉速為60 r/min，運用四階定步長龍格?庫塔數值積分法［16］求解式（32）。

由于所列動力學方程為時域積分求解，得到的序列需進行角域重采樣，進而得到IAS信號。無故障軸承以及外圈故障軸承的IAS信號如圖8所示。從圖8（a）中可看出無故障軸承呈現出的IAS擾動，擾動周期與角域故障周期基本相同，約為0.237（即1/0.237≈4.22）。圖8（b）顯示出滾子經過外圈故障區時IAS信號的明顯沖擊現象，并且觀察到IAS信號存在明顯周期性，擾動及沖擊特征的周期約為0.237（即1/0.237≈4.22）。角域上的外圈故障頻率理論值通過下式計算：

（36）

計算可得理論值為4.22 ev/r，即仿真結果與其基本吻合。對IAS信號進行階比譜分析得到如圖9所示的IAS階比譜，可發現在4.22×處出現明顯峰值，與外圈故障頻率理論值基本吻合，且存在明顯倍頻成分（2倍頻8.44×和3倍頻12.66×等）。

圖10為動力學模型仿真出的故障區寬度尺寸估計，可以看出從滾子進入故障區到退出故障區軸轉角為0.12 rad，故障區寬度估計的計算式為：

（37）

式中，θshaft為滾子通過故障區時軸轉過的角距離；fcage、fshaft分別為保持架和軸的轉速。

估算可得，仿真故障區寬度約為0.96 mm，與假設寬度1 mm基本吻合，模型可以有效識別沿滾道方向的缺陷長度。

3.2 實驗驗證分析

為驗證模型的正確性，以圓柱滾子軸承NU204 ECP外圈局部剝落為例，利用電火花方法在滾子軸承外圈滾道上加工出一個寬Ld=1 mm的損傷區域來模擬剝落故障（如圖11所示）。

使用皮帶輪給故障軸承提供徑向載荷，在旋轉機械故障模擬基礎試驗臺（如圖12所示）進行實驗。

圖13為轉速在60 r/min條件下采集到的滾子軸承外圈故障所對應的IAS擾動響應，顯示故障區沖擊周期約為0.237（即1/0.237≈4.22），對IAS帶來的影響與仿真結果基本一致。對比分析驗證了本文所建立的三自由度外圈故障引起IAS擾動動力學模型的正確性。

4 結" 論

本文基于Hertz線接觸理論，通過分析滾子在徑向載荷作用下，不同角位置與內、外圈滾道間的彈性變形量，將缺陷區帶來的附加力矩以及附加位移考慮進扭矩擾動中，并使切向力與法向接觸力相耦合，結合圓柱滾子軸承結構參數以及故障區尺寸，建立包含兩個位移自由度和一個旋轉自由度的軸承三自由度IAS動力學模型。在相同實驗條件下進行仿真與實測數據的對比分析研究，研究結果較好地解釋了外圈故障狀態下滾子軸承IAS擾動的來源，正確估計出故障區寬度尺寸，有效地反映了外圈故障對IAS帶來的影響。若進一步考慮內圈或者滾動體局部缺陷產生的附加力矩對系統的扭矩擾動，可以揭示軸承內圈、滾動體故障對IAS帶來的影響。本文所建動力學模型為低速工況下有效運用IAS進行故障診斷提供了理論支持。

參考文獻：

［1］"""" 姚燦江， 魏領會， 王海龍. RV減速器滾動軸承動態接觸應力的有限元分析［J］. 北方工業大學學報， 2016， 28（3）： 60-65.

YAO Canjiang， WEI Linghui， WANG Hailong. Finite element analysis of dynamic contact stress of rolling bearing of RV reducer［J］. Journal of North China University of Technology， 2016， 28（3）： 60-65.

［2］"""" 聶傲男， 李迎春， 沈文亮， 等. RV減速器曲柄支承軸承和轉臂軸承受力的變化規律研究［J］. 軸承， 2022（5）： 9-15.

NIE Aonan， LI Yingchun， SHEN Wenliang， et al.Research on variation laws of forces on crank support bearings and turning arm bearings of RV reducer［J］. Bearing， 2022（5）： 9-15.

［3］"""" ZENG Q， FENG G J， SHAO Y M， et al. Planetary gear fault diagnosis based on an instantaneous angular speed measurement system with a dual detector setup［J］. IEEE Access， 2020， 8： 66228-66242.

［4］"""" MOUSTAFA W， COUSINARD O， BOLAERS F， et al. Low speed bearings fault detection and size estimation using instantaneous angular speed［J］. Journal of Vibration and Control， 2016， 22（15）： 3413-3425.

［5］"""" AIT SGHIR K， COUSINARD O， EL BADAOUI M， et al. Diagnostic of rolling element bearing in low speed regime using the instantaneous angular speed and cyclo-stationary tools［C］∥International Conference on Noise and Vibration Engineering?International Conference on Uncertainty in Structural Dynamics（ISMA-USD 2014）. Leuven， Belgium： 2014.

［6］"""" RéMOND D， ANTONI J， RANDALL R B. Editorial for the special issue on instantaneous angular speed （IAS） processing and angular applications［J］. Mechanical Systems and Signal Processing， 2014， 44（1-2）：1-4.

［7］"""" SAWALHI N， RANDALL R B. Simulating gear and bearing interactions in the presence of faults： Part Ⅰ. The combined gear bearing dynamic model and the simulation of localised bearing faults［J］. Mechanical Systems and Signal Processing， 2008， 22（8）： 1924-1951.

［8］"""" BOURDON A， ANDRé H， RéMOND D. Introducing angularly periodic disturbances in dynamic models of rotating systems under non-stationary conditions［J］. Mechanical Systems and Signal Processing， 2014， 44（1-2）： 60?71.

［9］"""" GOMEZ J L， BOURDON A， ANDRé H， et al. Modelling deep groove ball bearing localized defects inducing instantaneous angular speed variations［J］. Tribology International， 2016， 98： 270?281.

［10］""" THIBAULT N， BOURDON A， RéMOND D， et al. Dynamic model of a deep grooves ball bearing dedicated to the study of instantaneous angular speed of rotating assemblies［J］. Tribology International， 2022， 174： 107753.

［11］""" 羅茂林， 郭瑜， 伍星. 考慮沖擊力的球軸承外圈剝落缺陷雙沖擊現象動力學建模［J］. 振動與沖擊， 2019， 38（14）： 48-54.

LUO Maolin， GUO Yu， WU Xing. Dynamic modeling of the dual-impulse behavior produced by a spall on the outer race of a ball bearing considering impact forces［J］. Journal of Vibration and Shock， 2019， 38（14）： 48-54.

［12］""" STRIBECK R. Ball bearings for various loads［J］. Transactions of the ASME， 1907 （29）： 420-463.

［13］""" PEREIRA C M， RAMALHO A L， AMBRóSIO J A. A critical overview of internal and external cylinder contact force models［J］. Nonlinear Dynamics， 2011， 63（4）： 681-697.

［14］""" PATIL M S， MATHEW J， RAJENDRAKUMAR P K， et al. A theoretical model to predict the effect of localized defect on vibrations associated with ball bearing［J］. International Journal of Mechanical Sciences， 2010， 52（9）： 1193-1201.

［15］""" JONES A B. Ball motion and sliding friction in ball bearings［J］. Journal of Basic Engineering， 1959， 81（1）： 1-12.

［16］""" DORMAND J R， PRINCE P J. A family of embedded Runge-Kutta formulae［J］. Journal of Computational and Applied Mathematics， 1980， 6（1）： 19-26.

第一作者： 韓" 譙（1998―），女，碩士研究生。E-mail： hertaq817@163.com

通信作者： 郭" 瑜（1971―），男，博士，教授。E-mail： kmgary@163.com