行星齒輪典型斷齒故障的動力學仿真*

楊文廣,蔣東翔

(清華大學電力系統國家重點實驗室 北京，100084)

行星齒輪典型斷齒故障的動力學仿真*

楊文廣,蔣東翔

(清華大學電力系統國家重點實驗室 北京，100084)

基于一種改進的行星齒輪箱集總參數模型，將斷齒故障等效到時變嚙合剛度中，建立了直齒行星齒輪的太陽輪斷齒故障、行星輪斷齒故障和內齒圈斷齒故障的動力學模型。以某單級行星齒輪為研究對象，考慮扭轉方向外部激勵，對其正常和各斷齒故障狀態下的動力學模型進行求解，對內齒圈垂直振動加速度信號進行分析。研究結果表明：只考慮扭轉方向的外部激勵時，正常狀態下內齒圈最高點處的垂直振動加速度信號頻譜中幅值最高的共振峰與模態能量分布相關；在各斷齒故障狀態下，內齒圈垂直振動加速度的有效值均大于正常值，時域波形中可看到不同形式的沖擊，包絡譜中可看到對應的故障特征頻率及其倍頻成分。

風力機；行星齒輪箱；斷齒故障；動力學模型；故障機理

引 言

齒輪箱是雙饋風力發電機組的關鍵部件之一。統計表明，風力機故障停機時間20%由齒輪箱故障引起，而且齒輪箱故障維修費用高，嚴重影響風力機的運行安全和經濟效益[1]。由于行星齒輪箱結構復雜，同時承受時變風載和電氣負載，其振動信號為非線性時變信號，因此研究齒輪箱在各種故障情況下的振動特性對于齒輪箱的振動故障診斷具有重要意義。

根據數據來源的不同，對齒輪箱振動故障特性的研究可大致分為動力學模型驅動和數據驅動。本研究內容屬于前者。文獻[2]將行星齒輪箱的動力學模型分為3類，即純扭轉模型、剛性多體動力學模型(也稱集總參數模型)和柔性多體動力學模型。龍泉[3]建立了一個一級行星傳動兩級平行傳動的齒輪箱的純扭轉模型，考慮剛度激勵和誤差激勵，研究了齒根裂紋和齒面點蝕故障下的齒輪箱動態響應。Chaari等[4]建立了齒面點蝕和齒根裂紋故障狀態下的行星齒輪集總參數模型，研究了齒輪在這兩種故障激勵下的振動響應。文獻[5]建立了嚙合剛度隨齒根裂紋的裂紋深度變化模型，采用集總參數模型研究不同裂紋深度下平行齒輪副的振動響應。陳裴[6]利用ADAMS軟件進行船用行星齒輪太陽輪斷齒故障的仿真，指出太陽輪斷齒故障時嚙合頻率周圍出現調制現象，調制頻率為其旋轉頻率乘以行星輪數。韓振南等[7]研究了平行傳動齒輪系統中存在齒輪剝離故障時的動力學特性。文獻[8]基于行星齒輪的振動響應特性和嚙合沖擊的傳遞路徑，綜合考慮各嚙合成分之間的相位差，建立了行星齒輪箱的振動信號仿真模型。與純扭轉模型和柔性多體動力學模型相比較，集總參數模型可以以相對較小的計算復雜度、獲得較準確的振動信號細節[9]，非常適合于齒輪箱故障機理研究。傳統的集總參數模型建立在以行星架轉速旋轉的旋轉坐標系中，存在殼體的旋轉自由度難以約束等問題。

筆者提出了一種改進的剛性多體動力學模型[10]。在此基礎上建立直齒行星齒輪的太陽輪斷齒故障、行星輪斷齒故障和內齒圈斷齒故障的動力學模型，考慮扭轉方向的外部激勵和時變剛度等內部激勵，采用隱式Newmark算法對模型求解。通過對仿真信號進行分析，從動力學角度揭示了行星齒輪箱在各種斷齒故障狀態下的動力學特性。

1 改進的行星齒輪箱集總參數模型

改進的集總參數模型直接選擇慣性坐標系作為除行星輪外的其他部件的運動參考坐標系，避免了殼體旋轉自由度難以約束的問題。對各行星輪，選擇其運動參考坐標系作為繞太陽輪公轉、同時以相同角速度自轉的動坐標系。與傳統模型相比，每個行星齒輪新增一個公轉方向的角位移自由度。在該模型中，各行星輪的公轉角位移與行星架的自轉角位移相同。

筆者以直齒的行星輪和太陽輪嚙合的嚙合剛度矩陣為例，簡單介紹該模型的推導過程，詳細推導過程可參考文獻[10]。記太陽輪為第1個齒輪，行星輪為第2個齒輪，則這兩個齒輪間的嚙合變形協調方程為

(1)

其中：θrevp為行星輪的公轉角位移；φi為第i個齒輪對中心線到該齒輪運動參考坐標系x軸的夾角；xi,yi,θzi分別為第i個齒輪的水平振動、垂直振動和繞軸向(z軸)的扭轉角位移；ri為第i個齒輪的節圓半徑；rc為行星架的半徑；廣義自由度q=[x1,y1,θz1,x2,y2,θz2,θrevp]。

嚙合力在第1個齒輪的各運動方向上的投影Fx1，Fy1，Tz1的計算公式為

(2)

嚙合力在行星輪的各運動方向上的投影Fx1，Fy1，Tz1，Trevp的計算公式為

(3)

其中：Trevp為行星輪受到的嚙合力在其公轉自由度方向上的投影。

定義廣義載荷為

F=[Fx1，Fy1，Tz1，Fx2，Fy2，Tz2，Trevp]

(4)

利用式(5)，得到行星輪與太陽輪嚙合的嚙合剛度矩陣Ke，該矩陣的大小為7×7

(5)

其他的剛度矩陣和阻尼矩陣可采用類似方法推導得到。根據各個齒輪間的相互作用關系，進行單元裝配和約束處理，得到行星齒輪級動力學數學模型為

(6)

其中：M，C，K分別為總體質量、總體阻尼和總體剛度矩陣；F(t)為外部激勵。

當阻尼項C中各個阻尼取值比較困難時，采用瑞利阻尼近似

C=αM+βK

(7)

其中：α，β為瑞利阻尼系數，根據具體情況選擇。

為提高數值求解的穩定性，將式(6)進行變換

(8)

其中：Feq(t)為等效激勵；F(t)為外部激勵；K(t)q為內部激勵；Kman為時不變的綜合剛度矩陣。

2 斷齒故障時變嚙合剛度模型

齒輪副在一個時刻可能有多對齒參與嚙合，一對嚙合的直齒的嚙合剛度隨嚙合位置在不斷變化。對于行星齒輪系統，還需要考慮各對嚙合齒輪之間的相對相位關系。因此，采用一個簡化的行星齒輪時變嚙合剛度模型[10]。對于斷齒故障，當斷齒參與嚙合時，其提供的嚙合力為0，可通過在時變嚙合剛度模型中加入這一過程來模擬斷齒故障。

以使用的行星齒輪為仿真對象，假設行星架轉速恒為120 r/min，基于上述時變剛度模型仿真得到正常狀態、太陽輪單齒斷全齒故障、內齒圈單齒斷全齒故障以及1號行星輪單齒斷全齒故障的部分齒輪副的時變嚙合剛度，如圖1～4所示。將故障狀態下的時變嚙合剛度與正常狀態下的相比，可以看到斷齒會導致時變剛度中出現周期性的負沖擊。從圖2看到，當太陽輪發生斷齒故障時，太陽輪與不同行星輪的時變嚙合剛度中由于斷齒導致的時變剛度負沖擊之間存在120°相位差。從圖4發現，斷齒行星輪分別與太陽輪和行星輪嚙合時的時變嚙合剛度的負沖擊之間存在180°相位差。

圖1 正常行星齒輪中部分齒輪對的時變嚙合剛度Fig.1 Mesh stiffness of some gear pairs in the healthy planetary gear

圖2 太陽輪斷齒故障時部分齒輪對的時變嚙合剛度Fig.2 Mesh stiffness of some gear pairs in the planetary gear with sun tooth breakage fault

圖3 內齒圈斷齒故障時部分齒輪對的時變嚙合剛度Fig.3 Mesh stiffness of some gear pairs in the planetary gear with ring tooth breakage fault

圖4 行星輪斷齒故障時部分齒輪對的時變嚙合剛度Fig.4 Mesh stiffness of some gear pairs in the planetary gear with planet tooth breakage fault

3 行星齒輪參數及固有頻率分析

選擇某單級行星齒輪進行研究[11]，該行星齒輪系統有3個行星輪，等間隔布置，詳細參數如表1所示。其中：下標s，p，c，r分別為太陽輪、行星輪、行星架和內齒圈；Ksn為太陽輪與行星輪的綜合嚙合剛度；Kpw為行星輪軸承的扭轉剛度，忽略軸承的扭轉剛度。通過計算可知，該級行星輪的傳動比約為4.56，輸入側等效轉動慣量。根據經驗選擇瑞利阻尼系數α=0．012 56，β=2．12×10-7。

表1 行星齒輪的參數Tab.1 Parameters of the planetary gear

將正常及斷齒故障狀態下的時變剛度模型分別帶入到改進的行星齒輪集總參數模型中，得到該行星齒輪在正常和各斷齒故障狀態下的動力學模型。

對該行星齒輪進行模態分析，得到的固有頻率如表2所示。該行星齒輪系統共6個旋轉固有頻率和6個平動固有頻率。由于模型中各軸承的水平和垂直方向的剛度相同，因此12個平動固有頻率兩兩相同。

表2 固有頻率計算結果Tab.2 Results of the modal frequency analysis

4 正常狀態仿真與信號分析

只考慮扭轉方向上的外部激勵，取行星架上驅動扭矩F1(t)=10 Nm,太陽輪上負載扭矩F2為

(9)

其中：Tr為齒輪箱的傳動比；v(t)為太陽輪的當前轉速；C=120 r/min；t0為太陽輪轉速首次達到轉速C時的時間。

定性分析可知，t0時刻以前傳動鏈空載運行，t0時刻負載接入，之后齒輪箱上的總扭矩達到動態平衡，轉速將在轉速C上下波動。可以估算t0≈CJ/F1=0．156s。選擇定步長隱式Newmark算法[12]對模型進行求解，仿真步長取1×10-6s，時長為5 s，采樣頻率為25 kHz。

圖5 正常狀態下各部件轉速Fig.5 The rotational speed of some components in the healthy planetary gear

進行正常狀態下行星齒輪的瞬態動力學分析，得到部分組件轉速的前2 s的時域波形如圖5所示。仿真得到t0≈0．158 s，與理論估計值之間誤差約為1.3%。取各轉速最后1 s的數據進行平均，作為各部件的平衡轉速，如表3所示。其中，各個理論值是以行星架轉速120 r/min為基礎，通過傳動比關系計算得到。可以看到，各仿真值與理論值之間的誤差在1%以內。

表3 轉速的計算值與理論值對比Tab.3 Comparisons between the calculation rotational speeds and theoretical rotational speeds

實際振動傳遞器采集到的振動信號應對應仿真的垂直振動加速度和旋轉振動加速度信號的合成，為簡便起見，筆者取內齒圈最高點的垂直振動加速度1～5s的信號進行信號分析，該點處旋轉運動在垂直方向上投影為0。其時域波形、頻譜和包絡譜分別如圖6所示。由于仿真模型只考慮了扭轉方向的外部激勵以及內齒圈柔性等因素，因此仿真得到的垂直振動加速度很小。可以看到，最大的峰值在1103 Hz處，其次為743，1894和2277 Hz處，高頻段也有兩個共振區。與表2中的固有頻率分析結果比較發現，這些頻率值與橫振固有頻率相吻合。

圖7為包絡譜幅值的局部放大圖,可以看到微弱的192 Hz嚙合頻率成分，低頻區中的成分比較雜亂。這說明該頻譜中的主要能量集中在固有頻率附近，驗證了行星齒輪振動信號中的共振調制現象。

圖6 正常狀態下內齒圈垂直振動的仿真結果Fig.6 The simulated signal of the ring vertical vibration acceleration in normal state

圖7 正常狀態下內齒圈垂直加速度頻譜和包絡譜局部放大圖Fig.7 Partial enlarged view of spectrum and envelop amplitude spectrum of the ring vertical vibration acceleration in normal state

對頻譜中的共振調制現象進行進一步分析。圖8為系統各個平動模態的能量隨組件的分布。模態能量的定義可參考文獻[13]，圖(a)～(f)分別對應第1個到最后一個平動固有頻率。每一副子圖均有4列，第1～4列分別對應該固有頻率的模態能量在行星架、內齒圈、太陽輪和所有行星輪上的分布。可以發現，1103 Hz(對應著子圖(b))的模態能量在內齒圈上分布的能量最大，其次為743 Hz。這說明內齒圈振動信號上的各固有頻率附近的能量值與該固有頻率對應的模態分布在該部件上的能量大小具有一定的關系。

圖8 各平動模態的能量隨組件的分布Fig.8 Energy distribution over components of each modal

5 斷齒故障仿真與信號分析

采用相同的外部激勵，對各斷齒狀態下的行星齒輪進行仿真，仿真時長為5 s，得到的各個轉速時域波形圖與圖5類似。行星架轉速為120 r/min時，該行星齒輪太陽輪斷齒故障特征頻率fs約為7.11Hz，內齒圈斷齒故障特征頻率fr約為2 Hz，該頻率同時也是行星架的旋轉頻率fc，行星輪斷齒故障特征頻率fp約為5.48 Hz。仍采用1～5 s的內齒圈垂直振動加速度信號進行信號分析，各狀態下加速度信號的時域波形如圖9所示。對比發現：各斷齒故障產生的沖擊導致該加速度的峰值明顯增大；不同的齒輪斷齒故障沖擊造成的波形不同。

圖9 各個狀態下該加速度信號的時域波形Fig.9 Comparisons among the waveforms of in different states

以太陽輪斷齒故障為例進行信號分析。該狀態下內齒圈垂直振動加速度的時域波形、頻譜和包絡譜如圖10所示。通過與正常狀態下該信號的對比發現，頻譜上1 103 Hz頻率處的振動幅值增長最多，導致高頻部分的成分被淹沒。圖11為故障狀態下該信號的頻譜和包絡譜的局部放大圖。與正常狀態下的信號對比發現，故障狀態下信號的頻譜中存在明顯的共振調制現象。從包絡譜的局部放大圖可以看到，除行星架的轉頻外，被調制的頻率在低頻段主要為7.11Hz以及其倍頻，而且其3倍頻較為突出，與行星輪的個數一致。

內齒圈斷齒故障和行星輪斷齒故障的內齒圈垂直振動加速度信號中也存在類似的共振調制現象。圖12為不同狀態下包絡譜局部放大圖的對比。可以發現，不同故障狀態下被調制的頻率成分不同，太陽輪和行星輪斷齒故障狀態下，被調制的主要頻率是對應的故障特征頻率，在其倍頻成分中，3倍頻成分較為突出，與行星輪的個數一致。行星輪斷齒故障狀態下，故障特征頻率主要為其故障特征頻率的偶數倍頻，其中2倍頻成分最突出，這是由于斷齒的行星輪在一個旋轉周期中先后與太陽輪和內齒圈嚙合的原因。

圖10 太陽輪斷齒故障狀態下內齒圈垂直振動加速度仿真結果Fig.10 The simulated signal of ring vertical vibration acceleration in the sun tooth break fault state

圖11 太陽輪斷齒故障狀態下內齒圈垂直振動頻譜和包絡譜的局部放大圖Fig.11 Partial enlarged view of spectrum and envelop spectrum of ring vertical vibration acceleration in the sun tooth break fault state

圖12 內圈圈垂直振動加速度的包絡譜局部放大圖Fig.12 Partial enlarge view of the envelop spectra of the ring vertical vibration acceleration

共振解調分析(demodulated resonance analysis，簡稱DRA)是齒輪箱故障信號分析常用的方法，筆者選擇高頻區5 000～10 000 Hz的頻帶進行共振解調分析，得到各狀態下包絡譜信號的局部放大圖如圖13所示。與包絡譜對比可以發現，共振解調得到的頻譜更加干凈，故障特征頻率更加突出；對行星輪斷齒故障，共振解調分析得到的故障特征頻率的偶數倍頻成分的幅值相對較大。這些現象與文獻[14]的定性分析和實驗結論一致。

圖13 內圈圈垂直振動加速度的共振解調分析局部放大圖Fig.13 Partial enlarge view of the DRA of the ring vertical vibration acceleration

6 結 論

1) 只考慮扭轉方向的外部激勵時，齒輪箱最高點處垂直振動加速度信號中共振頻率處的能量與其對應的固有頻率的模態能量分布存在一定的關系。

2) 正常狀態下，該振動加速度的調制頻率中存在嚙合頻率，但較為微弱，低頻區中的信號無明顯規律。

3) 當出現斷齒故障時，該振動加速度的有效值增大，頻譜中出現了明顯的共振調制現象。包絡譜分析顯示，被調制的頻率成分為對應的故障特征頻率及其倍頻。

4) 對太陽輪和內齒圈斷齒故障，除1倍頻外，其3倍頻成分也較為突出，與行星輪的個數一致；對行星輪斷齒故障，故障特征頻率的2倍頻成分最為突出，這些結論與相關文獻中的分析一致。

5) 建立的行星齒輪典型斷齒故障的動力學模型可推廣用于復雜外部激勵作用下的故障機理研究以及多級行星齒輪箱和風力機傳動鏈的故障機理研究。

[1] 陳雪峰，李繼猛，程航，等.風力發電機狀態監測和故障診斷技術的研究與進展[J].機械工程學報，2011,47(9):45-52.

Chen Xuefeng,Li Jimeng,Cheng Hang,et al.Research and application of condition monitoring and fault diagnosis technology in wind turbines[J].Journal of Mechanical Engineering,2011,47(9):45-52.(in Chinese)

[2] Peeters J.Simulation of dynamic drive train loads in a wind turbine[D].Belgium:Katholieke Universiteit Leuven,2006.

[3] 龍泉.風電機組齒輪傳動系統動態特性及故障診斷方法研究[D].保定：華北電力大學,2012.

[4] Chaari F,Fakhfakh T,Haddar M.Dynamic analysis of a planetary gear failure caused by tooth pitting and cracking[J].Journal of Failure Analysis and Prevention，2006,6(2):73-78.

[5] Wu Siyan,Zuo Mingjian,Parey A.Simulation of spur gear dynamics and estimation of fault growth[J].Journal of Sound and Vibration,2008,317(3-5):608-624.

[6] 陳裴.船用行星齒輪斷齒故障機理研究[D].上海：上海交通大學,2014.

[7] 韓振南，孫文婷，高建新.含輪齒剝落的齒輪系統動力學故障模擬[J].振動、測試與診斷，2012,32(1):101-104.

Han Zhennan,Sun Wenting,Gao Jianxin.Dynamics fault simulation of gear transmission system including spalling[J].Journal of Vibration,Measurement &Diagnosis,2012,32(1):101-104.(in Chinese)

[8] Lei Yaguo,Tang Wei,Kong Detong,et al.Vibration signal simulation and fault diagnosis of planetary gearboxes based on transmission mechanism analysis[J].Journal of Mechanical Engineering,2014,50(17):61-68.

[9] Ambarisha V K,Parker R G.Nonlinear dynamics of planetary gears using analytical and finite element models[J].Journal of Sound and Vibration,2007,302(3):577-595.

[10] Yang Wenguang,Jiang Dongxiang.An improved rigid multibody model for the dynamic analysis of the planetary gearbox in a wind turbine[J].Shock and Vibration,2016,2016(8):1-18.

[11] Lin Jian,Parker R G.Analytical characterization of the unique properties of planetary gear free vibration[J].Journal of Vibration &Acoustics，1999,121:316-321.

[12] Gasmi A,Sprague M A,Jonkman J M,et al.Numerical stability and accuracy of temporally coupled multi-physics modules in wind-turbine CAE tools[R].[S.l.]:Office of Scientific &Technical Information Technical Reports,2013.

[13] Lin Jian,Parker R G.Sensitivity of planetary gear natural frequencies and vibration modes to model parameters[J].Journal of Sound and Vibration,1999,228(1):109-128.

[14] 馮志鵬，趙鐳鐳，褚福磊.行星齒輪箱齒輪局部故障振動頻譜特征[J].中國電機工程學報,2013,33(5):119-127.

Feng Zhipeng,Zhao Leilei,Chu Fulei.Vibration spectral characteristics of localized gear fault of planetary gearboxes [J].Proceedings of the CSEE,2013,33(5):119-127.(in Chinese)

10.16450/j.cnki.issn.1004-6801.2017.04.019

* 國家自然科學基金資助項目(51174273)

2015-12-16；

2016-02-23

TK83;TH17

楊文廣，男，1988年9月生，博士生。主要研究方向為風力機故障診斷。曾發表《An improved rigid multibody model for the dynamic analysis of the planetary gearbox in a wind turbine》(《Shock and Vibration》2016,No.8)等論文。 E-mail:ywg16@mail.tsinghua.edu.cn

蔣東翔，男，1963年2月生，教授、博士生導師。主要研究方向為動力系統設備故障診斷技術研究與應用。 E-mail:jiangdx@mail.tsinghua.edu.cn