裂紋故障對輪齒時變嚙合剛度的影響分析

孟 宗 李佳松 潘作舟 龐修身 崔玲麗 樊鳳杰

1.燕山大學電氣工程學院，秦皇島，0660042.北京工業大學機械工程與應用電子技術學院，100124

0 引言

齒輪傳動系統在工業領域應用范圍極為廣泛[1],與人類生產生活息息相關。齒輪作為旋轉機械中改變轉速與傳遞扭矩的重要零部件之一，長期處于復雜、變化環境下，易產生磨損、裂紋、點蝕、剝落、斷齒等故障[2]。裂紋作為齒輪故障主要存在形式之一，產生的危害輕則影響設備正常工作，降低人類生產效益，重則威脅人類生產生活安全，甚至造成機毀人亡慘劇，因此，對齒輪傳動系統進行故障檢測與健康監測具有重要意義。齒輪振動來源于系統內部激勵與外部激勵，剛度激勵作為內部激勵的主要來源，當齒輪嚙合時，裂紋使齒輪時變嚙合剛度(time-varying mesh stiffness,TVMS)發生不規則變化，進而引起齒輪異常振動，產生較大噪聲，因此，準確計算TVMS是構建動力學模型的關鍵所在。針對裂紋故障情況下的齒輪TVMS的計算及動力學分析，國內外學者做了許多研究。萬志國等[3]考慮齒根圓與基圓大小不重合現象，對TVMS算法進行改進，并與ISO標準比對，結果表明誤差在5%以內。張振等[4]提出雙齒根裂紋模型，指出在兩輪齒裂紋深度之和相等的前提下，具有最大裂紋深度模型對TVMS劣化率有較大影響。YANG等[5]構建了考慮齒廓間隙與軸承間隙的四自由度動力學參數模型，研究了不同故障程度下齒輪系統的動力學響應，結果表明故障程度越高，系統龐加萊截面圖與相圖越混亂。李秀紅等[6]通過ABAQUS確定裂紋萌生位置，分析影響齒輪疲勞壽命的多種因素，最終選取合適齒頂修緣量與增大表面硬化方法來延長齒輪使用壽命。MA等[7]考慮了裂紋分別沿基體和輪齒拓展的兩種情況，并根據有限元法得出裂紋在相同深度情況下，對沿基體方向拓展的TVMS影響最大的結論。CHEN等[8]采用解析有限元法對復雜基體類型的TVMS進行求解，結果表明與輪緣厚度相比，腹板厚度對TVMS的影響更大。吳家騰等[9]通過實驗臺采集裂紋輪齒的齒根應變，將其與有限元法計算得到的應變代入反問題模型中，實現了對TVMS的反向求解。MENG等[10]采用兩種方法計算了在不同裂紋等級下振動信號的統計指標，結果表明峭度指標對早期故障特征較為敏感。雷敦財等[11]參照圓柱直齒輪有限元模型構建方法，得出面齒輪傳動載荷分布不均勻、TVMS呈周期性正弦曲線變化的結論。

作為TVMS的兩種基本計算方法，有限元法與解析法各有利弊。有限元法是將齒輪模型離散為若干微小單元，并通過施加相應載荷條件與邊界條件計算系統各部分發生的位移；解析法則是在MATLAB中構建TVMS參數模型并輸入相應齒輪參數進行求解計算。相比有限元法，解析法求解速度較快但精度較差。因此，通過研究TVMS的變化規律并對模型算法進行合理改進，可以保證解析法運算效率以及結果的準確性。

本文首先修正了雙齒嚙合區TVMS計算公式，減小了計算誤差；其次，分析“死區”的產生原因，通過有限元法確定輪齒易發生斷裂的位置；構建了分別沿裂紋深度與長度方向的拓展曲線，定量研究了不同裂紋深度及長度對TVMS的影響規律；最后討論了不同裂紋深度下的輪齒對承載能力的變化過程。

1 健康齒輪TVMS計算

齒輪在正常工況下工作時，輪齒接觸數量與嚙合位置發生周期性改變，并使TVMS呈現周期性(高-低-高)變化規律[12]。在單雙齒交替時刻，剛度發生突變產生嚙合沖擊并導致振動幅值改變。因此，準確評估TVMS有助于描述齒輪傳動系統動態響應，為動力學參數模型的確立提供必要條件。

1.1 輪齒齒廓方程的確定

根據材料力學相關理論，將圓柱直齒輪輪齒簡化為固定在齒根圓上截面不斷變化的懸臂梁。輪齒齒廓共由三部分組成，如圖1所示，分別為齒頂圓圓弧(AB段)、漸開線齒廓(BC段)和基圓與齒根圓之間的過渡曲線(CD段)[13]，輪齒齒廓沿x軸對稱分布，在BC段上輪齒相互接觸并發生嚙合過程，BC段漸開線曲線公式為

x=Rb[sin(α2+α)+cosα]

(1)

hx1=Rb[cos(α2+α)-sinα]

(2)

圖1 健康輪齒懸臂梁Fig.1 Healthy tooth cantilever beam

過渡曲線雖未因接觸發生形變，但對分析輪齒彎曲形變有重要意義。近年來，諸多學者將過渡曲線簡化為一條直線[14]，但實際加工采用展成法切削方式，漸開線與過渡曲線分別由標準齒條的直線部分與齒頂圓弧部分切割而成，過渡曲線表達式為

(3)

(4)

(5)

式中，r為分度圓半徑；a1為齒頂圓弧圓心距中線距離；rρ為齒頂圓弧半徑；b1為齒頂圓弧圓心距齒條齒槽中心線的距離；γ為變參數，變化范圍為α0≤γ≤π/2；α0為分度圓壓力角。

1.2 改進勢能法計算公式

TVMS作為表征振動特性的一個重要指標，具有重要研究價值。根據彈性力學相關理論，接觸齒對所具備的總勢能U包括輪齒具有的勢能Utooth與輪體具有的勢能Ubody，Utooth包括赫茲接觸勢能Uh、彎曲勢能Ub、剪切勢能Us、軸向壓縮勢能Ua。根據勢能與剛度的對應關系，可推導出赫茲接觸剛度Kh、彎曲剛度Kb、剪切剛度Ks、軸向壓縮剛度Ka及基體剛度Kf。隨著齒輪轉動，嚙合點沿著發生線移動，產生不同變形量，引起勢能與剛度改變，其中，總勢能及每部分剛度計算公式為

(6)

(7)

(8)

(9)

(10)

(11)

(12)

式中，F為嚙合力；E為彈性模量；L為齒輪寬度；μ為泊松比；d為嚙合點到齒根圓距離；h為嚙合點距輪齒中心線高度；G為剪切模量；hx1、hx2分別為距離齒根圓x處且分別位于漸開線與過渡曲線處的齒廓與x軸的距離；Ix、Ax分別為距離齒根圓x處的截面慣性矩與截面面積。

計算基體變形量時，假設輪齒是剛性的，而輪體為一個彈性圓環，當齒對參與嚙合時，認為輪齒不會因接觸而變形，但會導致齒輪體產生彈性形變，進而引起輪齒發生形變。在雙齒嚙合區，當同時有第i與第j對(i≠j)輪齒參與嚙合過程時，傳統方法僅考慮第i(j)對齒的嚙合力在輪體作用下對第i(j)對齒沿發生線的變形量，卻忽視第i(j)對齒的嚙合力對第j(i)對齒的變形量。因此，低估由輪體偏轉作用下引起的輪齒總變形量，從而高估雙齒嚙合區TVMS數值，這種不同齒對之間的變形量互相影響的現象稱為相鄰齒耦合效應。MA等[15]通過有限元法確定修正系數并對雙齒嚙合區的基體剛度進行修正。隨后XIE等[16]基于彈性圓環理論，假設在齒根圓處法向應力呈三次分布、切向應力呈拋物線分布，并對基體剛度進行改進。相鄰輪齒耦合效應示意圖見圖2a，輪齒i與輪齒j在耦合效應作用下的變形量計算公式分別為

(13)

(14)

其中，δij(i≠j)為考慮相鄰齒耦合效應時，當第j個輪齒在嚙合變形時對第i個輪齒沿發生線造成的變形量；Fi、ui、αi定義方式如圖2b所示，i表示輪齒編號。Ai～Ii數值參考文獻[16]，參數Mi、Qi、Ti、Ui與Li、Pi、Ri、Si、Vi計算公式如下：

(15)

(16)

(a)相鄰齒耦合效應示意圖(b)輪齒受力放大圖圖2 相鄰齒耦合效應Fig.2 Coupling effect of adjacent teeth

在雙齒嚙合區，除了考慮由相鄰齒耦合效應引起的齒形變外，仍需考慮嚙合齒對i(j)受輪體彈性形變作用產生的變形量的影響，在兩者共同影響下引起基體變形，總的基體變形量公式為

δf1=F1δ11+F2δ12

(17)

δf2=F1δ21+F2δ22

(18)

(19)

(20)

(21)

則單、雙齒嚙合區TVMS的計算公式分別為

(22)

(23)

其中，δf1、δf2為在不考慮齒接觸變形時，輪齒沿嚙合線的總變形量；δij(i=j)為當第j個輪齒在嚙合變形時對第i個輪齒沿發生線造成的變形量；Lsr1、Lsr2為齒對負載分擔比。上標和下標p、g分別代表主、從動輪，i代表接觸齒對的數量，采用表1所示的齒輪參數，將基于改進勢能法的TVMS結果與其他方法進行對比，如圖3所示。

圖3 不同方法TVMS對比Fig.3 TVMS Comparison of different methods

由圖3可以看出，文獻[17]與其他兩種方法相比，在單雙齒嚙合區TVMS均相差過大，這是由于該方法視齒輪體為剛性，只考慮輪齒的形變作用，而忽略了輪體變形量的影響，使得TVMS計算結果過大。文獻[18]為傳統勢能法，與改進后勢能法相比結果相差不大，但在雙齒嚙合區未考慮由相鄰齒耦合效應引起的輪齒變形量，使得整體變形量偏小，而在單齒嚙合區由于齒輪變形量擬合方程系數差異性引起結果不同。因此，改進后勢能法與另外兩種方法相比，考慮范圍更為全面，能夠得到更準確的齒輪變形量，獲得更精確的TVMS數值。

2 裂紋齒輪TVMS計算

TVMS異常變化是齒輪故障及外界工況變化等諸多因素共同作用的結果。當旋轉機械長期在高強度的環境下工作時，在齒根附近會產生應力集中現象并生成微觀裂紋。隨著循環加載次數累積，形成宏觀裂紋并逐步傳播，從而導致TVMS隨之下降，使得振動信號故障特征愈發明顯，因此由裂紋故障引起的危害是不容忽視的問題。本節考慮裂紋在不同時期拓展情況，構建不同裂紋模型并分析不同裂紋深度與長度對TVMS的影響。

2.1 沿深度方向拓展裂紋模型的構建及TVMS計算

與健康輪齒相比，裂紋導致輪齒截面面積A(x)及截面慣性矩I(x)減小，其原因為在載荷作用下，這部分輪齒缺乏足夠機械強度，被稱為“死區”，表明該部分輪齒不但失去正常承載能力，而且會在彎曲應力持續作用下拓展范圍越來越大，直至產生斷齒故障并導致輪齒對的承載能力完全喪失。根據第一強度理論，最大拉應力為引起材料斷裂的最主要因素，以此為依據可確定裂紋萌生位置。文獻[19]指出裂紋萌生點位于漸開線齒廓上，結合文獻[20]的分析方法，在ANSYS Workbench靜力學模塊中，對單齒嚙合區一點施加500 N的載荷，并約束輪孔表面的6個自由度，主應力分析結果如圖4所示。結果表明齒輪最大拉應力點位于過渡曲線處(圖4中紅色區域)，而不是在漸開線齒廓上，因此，可確定裂紋萌生點多發生在過渡曲線附近[20]。

圖4 最大主應力分析結果Fig.4 Analysis results of maximum principal stress

在裂紋拓展相關研究中，多數情況將裂紋簡化為對稱直線，但由于裂紋拓展到輪齒對稱線時，拓展速度加快、斜率更大，故將下半部分裂紋形狀修正為斜率更大的圓弧。如圖5所示，q1為未穿過x軸裂紋深度，此時裂紋尖端位于E點；隨著故障程度加劇，裂紋深度到達x軸(G點)并穿過x軸到達H點處，將穿過x軸裂紋深度記為q2，這兩部分深度之和記為總裂紋深度q。本文中，裂紋拓展角度取常數θ=80°，研究裂紋是否穿過x軸兩種情形，A(x)與I(x)計算公式如下。

圖5 裂紋輪齒懸臂梁Fig.5 Crack tooth cantilever beam

情形1：當裂紋深度在x軸上方，如圖5中E點所示，即裂紋總深度q≤q1max且q≥0時，有

(24)

(25)

hc=q1/sinθ

(26)

情形2：當裂紋深度在x軸下方，如圖5中H點所示，即裂紋總深度q=q1max+q2>q1max時，有

A(x)=

(27)

I(x)=

(28)

(29)

(30)

由式(24)～式(30)可知，在由齒根圓向齒頂圓嚙合過程中，情形1與情形2均會導致A(x)與I(x)減小。為更直觀地分析輪齒在不同裂紋深度下TVMS的變化規律，將裂紋深度取整，數值分別為0、1 mm、2 mm、3mm，且每種深度都有相應故障程度并對應表2，TVMS變化結果如圖6所示。

表2 裂紋深度參數

由圖6看出，裂紋使單雙齒嚙合區TVMS結果發生改變，當裂紋深度較小時，TVMS減小但變化不明顯，且隨著深度增加降幅逐漸變大。在第一個雙齒嚙合區比第二個雙齒嚙合區降幅要小，這是由于在初始嚙合階段健康齒對起主導作用。同時，負載分擔比作為衡量齒對承載能力變化的一個重要指標，能夠更加直觀地描述齒對在不同故障程度下承載能力的變化規律及輪齒損傷程度的變化過程。負載分擔比的變化如圖7所示。

圖6 不同裂紋深度TVMS變化Fig.6 TVMS changes with different crack depths

圖7 不同裂紋深度下負載分擔比Fig.7 Load sharing ratio for different crack depths

由圖7看出，在不同裂紋深度下，裂紋齒對的負載分擔比不同，并對第二個雙齒嚙合區結果影響較大。在單齒嚙合區負載分擔比為1，因為此時只有裂紋齒對參與嚙合；隨著裂紋深度的增加，雙齒嚙合區負載分擔比逐漸降低且變化幅度增大，而造成這一現象的原因為輪齒“死區”范圍拓展速度加快，使得裂紋齒對承載能力下降，直接印證了裂紋會縮短齒輪使用壽命并且降低齒輪系統傳動質量這一觀點。

2.2 沿長度方向拓展裂紋模型的構建及TVMS計算

裂紋在中早期拓展時，并不是等深度且完全貫穿的[5]，會在過渡曲線附近同時沿深度與長度方向拓展，直至完全貫穿整個齒寬。因此，本文針對這一現象，將中早期裂紋模型簡化為沿齒寬方向不等深度分布的裂紋模型，如圖8a所示，qs與qe分別為始端與尾端裂紋深度。沿齒寬方向將輪齒分割成無數小薄片，任取一小薄片如圖8b所示，當該小薄片厚度很小時，裂紋深度是確定的，但小薄片間裂紋深度不斷變化且都沿固定角度θ拓展。每一小薄片記為Δx，其厚度與剛度分別為ΔL、Kpiece，健康與故障輪齒所具有的剛度Kthealth與Ktcrack記為

(31)

(32)

(33)

N=L/ΔL

(34)

當只有裂紋齒對參與嚙合過程時，時變嚙合剛度K記為

(35)

式中，Ka(x)、Kb(x)、Ks(x)與ΔKa(x)、ΔKb(x)、ΔKs(x)分別為整個健康輪齒與裂紋輪齒任一小薄片的軸向壓縮剛度、彎曲剛度與剪切剛度；N為被分割后裂紋輪齒小薄片份數，求解時需同時保證結果精確性與運算速度，取N=60。

(a)非均勻分布裂紋懸臂梁模型 (b)切割后裂紋輪齒小薄片圖8 非均勻分布裂紋模型Fig.8 Non-uniform distribution crack model

當裂紋沿齒寬方向投影長度小于齒寬時，稱為非貫穿型裂紋，反之稱為貫穿型裂紋，需同時討論這兩種情形下沿齒寬方向的裂紋深度分布情況，如圖8a與圖9所示。在裂紋深度(A-A)截面上，建立以初始裂紋萌生點O為坐標原點，OC與OB方向分別為X軸與Y軸的平面直角坐標系，將裂紋深度分布曲線按照拋物線處理，沿長度方向拓展的解析方程計算公式如下。

(1)貫穿型裂紋模型：

(36)

圖9 非均勻分布裂紋軌跡Fig.9 Non-uniform distribution crack trajectories

(2)非完全貫穿型裂紋模型：

(37)

在本文中，研究在中度故障程度前提下，分析貫穿型與非貫穿型裂紋對TVMS的變化規律，裂紋參數與其對應的結果分別見表3與圖10。

表3 非均勻分布模型裂紋參數

圖10 不同裂紋長度TVMS變化Fig.10 TVMS changes with different crack lengths

由圖10可知，在初始裂紋深度相同前提下，由不同裂紋尾端深度與裂紋投影長度引起的TVMS不同，且TVMS會隨著兩者的增加而減小。對比圖10結果與圖6裂紋深度結果可知，由裂紋長度與尾端深度引起的TVMS變化并不明顯；并且對比裂紋深度和裂紋長度對TVMS影響可知，裂紋深度對TVMS影響更大，說明當故障發生在中早期時，由故障產生的剛度激勵變化并不明顯，導致由輪齒故障引起的動態行為不易被監測。

3 有限元結果驗證

3.1 齒輪模型建立及網格劃分

以表1參數為基準，在SolidWorks中構建不同故障程度的圓柱直齒輪裝配模型，相應故障程度的輪齒部位放大圖見圖11，模型保存為x_t文件并通過SolidWorks與Workbench接口將齒輪模型導入Workbench瞬態動力學模塊中，模型被分割成輪齒與輪體兩部分，對三對嚙合輪齒的網格進行加密，如圖12所示。主動輪和從動輪只具備自身對地的旋轉自由度[21]，旋轉副采用MPC184單元，接觸區采用CONTA174單元和TARGE170單元，齒輪副為四面體與六面體混合網格，即輪體與輪齒網格分別為六面體與四面體。由于輪齒的故障程度不同，采用四面體網格能夠適應裂紋復雜形貌，使其求解時具備更好的適應性，此外假定材料為線彈性材料，其余網格設置不變。

(a)q=0 (b)q=1 mm

(c)q=2 mm (d)q=3 mm圖11 嚙合輪齒放大圖Fig.11 Enlargement of meshing gear teeth

圖12 有限元齒輪網格Fig.12 Finite element gear mesh

3.2 有限元法TVMS的求解

為更準確地分析齒輪副動態行為，并對2.1節結果進行驗證，在接觸設置中，選擇互相接觸的上下表面分別為接觸面與目標面，建立圖13所示的無摩擦力兩自由度齒輪副動力學參數模型[22]，并將齒輪副簡化為帶有質量的質心與形心重合的兩個圓盤。求解前，需將齒輪模型旋轉并使第二對輪齒恰好在初始嚙合位置上。求解結束后，提取對應時刻主從動輪輪孔相對旋轉角，分別記為θp(t)與θg(t)，由動力學方程可得

(38)

(39)

x(t)=rpbθp(t)-rgbθg(t)-e

(40)

式中，Jp、Jg分別為主、從動輪慣性矩；Tp、Tg分別為驅動力矩與負載力矩；rpb、rgb分別為主、從動輪基圓半徑；C(t)、Kt分別為嚙合阻尼與嚙合剛度；x(t)為輪齒沿嚙合線的變形量；2e為齒輪總間隙。

當齒輪副旋轉速度較慢時，由齒輪副轉動引起的慣性效應及嚙合阻尼很小[21]，可剔除式(38)、式(39)與轉角相關的導數項，簡化后的嚙合剛度Kt為

(41)

圖13 兩自由度動力學參數模型Fig.13 2-DOF dynamic parameter model

在前處理階段，需設置如齒輪轉速與扭矩等相關邊界條件，本文以第二對輪齒開始嚙合到退出嚙合的整過程為研究對象，同時結合分析前后兩組齒對結果的影響，如圖14紅色虛線框所示，規定第二對齒開始進入嚙合(i=1)時角位移為0，單雙齒區角位移范圍計算公式如下：

(42)

(43)

θt=2π(Cr-1)/Z1

(44)

(45)

式中，θs、θd分別為主動輪位于單齒嚙合區與雙齒嚙合區輪齒角位移范圍；θt為雙齒嚙合區角位移；Z1為主動輪齒數；Cr為重合度；rai、rbi、ri(i=1, 2)分別為主從動輪齒頂圓、基圓、分度圓半徑；α為分度圓壓力角。

經式(42)～式(45)計算，確定主動輪旋轉角度為0.51 rad，轉速為0.51 rad/s，從動輪扭矩為200 N·m，求解時間為1 s，共分為50子步，即針對四種不同的裂紋深度情況，在0～1 s內每隔0.02 s分別采集主、從動輪的輪孔相對旋轉角。在后處理階段，將得到的輪孔相對旋轉角導入MATLAB中并根據式(41)計算嚙合剛度Kt，結果如圖15所示。

圖14 輪齒角位移變化Fig.14 Change of tooth angle displacement

圖15 有限元分析結果Fig.15 Finite element analysis results

由圖15可知，TVMS隨裂紋深度增加而降低，且裂紋對單雙齒嚙合區均有影響，并對第二個雙齒嚙合區結果影響較大，這與改進的勢能法所得結論一致；當裂紋深度達3 mm時，TVMS降幅變大，由于裂紋深度超過輪齒中心線，輪齒上半部分完全失去承載能力。值得注意的是，此時有限元結果存在一個“轉接區”，因為在嚙合時輪齒處于單雙齒嚙合交替的過渡階段，使得齒輪總變形量位于單齒變形量與雙齒變形量之間，能夠更準確地描述齒輪傳動的動態過程。

為研究有限元法與改進勢能法結果的差異性，將兩者的平均TVMS進行對比，對比結果見表4。當裂紋深度分別為0、1 mm、2 mm、3 mm時，與有限元法所得結果的相對誤差分別為5.9%、5.0%、2.6%、3.5%，結果在可接受誤差范圍內，為齒輪動力學模型的確立提供了一定參考。

表4 勢能法與有限元法對比結果

4 結論

(1)本文將輪齒簡化為固定在齒根圓變截面的懸臂梁，研究了由齒頂圓弧、漸開線齒廓及過渡曲線構成的完整齒廓曲線。針對傳統勢能法的局限性，分析了雙齒嚙合區相鄰齒耦合效應對TVMS的影響，將所得結果與另外兩種傳統算法進行了對比，結果表明傳統方法低估了輪齒總變形量，改進勢能法能夠實現雙齒區TVMS的降低。

(2)明確了“死區”出現的原因及危害性，探索了最大拉應力點(材料失效部位)出現的位置。研究了裂紋對截面面積及慣性矩的影響并推導了相應的計算公式，提出了一種沿深度方向拓展的裂紋路徑，分析了不同深度下的TVMS及負載分擔比的變化規律。研究了中早期故障程度下裂紋同時沿深度和長度方向拓展的情況，構建了一種非均勻分布裂紋模型并提出相應的TVMS計算方法，結果表明TVMS對裂紋深度較為敏感，而對裂紋長度的敏感性較差。

(3)構建了兩自由度齒輪傳動系統動力學參數模型，在ANSYS Workbench瞬態動力學模塊提取主從動輪輪孔相對旋轉角，并求得TVMS大小。結果表明有限元法與勢能法結果較為吻合，誤差均在6%的范圍內。

本文僅考慮了裂紋對TVMS的影響規律，下一步研究由裂紋故障引起的動力學響應并分析相應的故障特征，進而討論多種評價指標(時域、頻域)的變化趨勢，從而實現裂紋故障早期故障診斷的目的。