陳玲，魏沛堂，余泳，劉懷舉，朱才朝，周浩

(1.重慶大學機械傳動國家重點實驗室，重慶，400044；2.貴州群建精密機械有限公司，貴州遵義，563003)

齒輪是機械傳動中的重要基礎件，廣泛應用于風電、航空、航天、汽車和船舶等領域裝備的運動及動力傳遞，被譽為工業(yè)的象征。隨著現(xiàn)代高端裝備對齒輪傳動高功率密度和高承載能力需求的日益提升，齒輪的接觸疲勞問題成為限制這些機械設備壽命和可靠性的重要原因。影響齒輪疲勞失效的機制復雜，出現(xiàn)早期微點蝕、點蝕、剝落、齒面斷裂等接觸疲勞失效形式，并且失效形式隨著材料、結構、工況等參數(shù)的不同而變化[1?3]。經(jīng)典疲勞研究認為，材料微觀結構特征從根本上決定齒輪等結構件的服役性能和疲勞壽命[4]。SANGID等[5?6]詳細綜述了疲勞裂紋形成過程中微觀結構的根本性作用，強調高周疲勞壽命主要受裂紋萌生過程影響，而裂紋萌生顯著受到微觀結構特征影響。因此，引入材料微觀結構特征研究齒輪疲勞失效機理具有重要意義。

接觸疲勞的數(shù)值模擬是預測結構/材料行為和優(yōu)化結構/材料力學性能的重要工具，將材料微觀結構信息納入數(shù)值模型是材料計算細觀力學的挑戰(zhàn)之一。為了解決這個問題，國內外學者提出了不同的方法來考慮微觀結構特征的影響。WARHADPANDE 等[7?8]利用Voronoi 剖分法模擬軸承鋼的微結構，開發(fā)了彈塑性Voronoi 有限元模型，預測的疲勞剝落形狀與實驗結果吻合較好。NOYEL 等[9]構建了考慮材料微觀結構的晶體，結合力有限元模型模擬軸承鋼滾動接觸疲勞，發(fā)現(xiàn)在三叉晶界處出現(xiàn)應力畸點，易形成微裂紋。WANG 等[10]通過修正Voronoi 多邊形法考慮實際晶粒尺寸、取向以及珠光體和鐵素體的比例，有效模擬了真實顯微組織結構的疲勞行為。PAUL[11]建立了基于顯微組織DP 590 鋼變形行為、塑性應變局部化和塑性失穩(wěn)的細觀力學模型。此外，多相有限元法在模擬材料微觀結構方面也得到了廣泛應用，該方法將不同相位特性分配給單元中的各個積分點，與單相有限元不同，多相有限元法的網(wǎng)格與材料的相結構無關，可以使用有限元網(wǎng)格來模擬復雜材料中微觀結構的行為[12]。

為了充分考慮多晶體的微觀結構，如晶粒和晶界的尺寸、形狀、空間取向、分布等特征對材料疲勞性能的影響，需要進行詳細的實驗研究，獲取材料的真實微觀幾何拓撲結構。HOLLISTER等[13]提出了基于數(shù)字圖像的建模技術，將骨的微觀形態(tài)影響納入生物工程有限元模擬中。丁友平等[14]將計算機圖像處理技術運用于晶粒度級別測定。KATANI 等[15]利用MATLAB 圖像處理軟件，實現(xiàn)了鐵素體和馬氏體的雙相識別。LEWIS等[16?17]利用連續(xù)切片技術和EBSD生成的微觀結構，重構了與材料真實微觀結構相似的三維拓撲幾何模型。LIAN 等[18]采用EBSD 表征方法與圖像處理技術區(qū)分鐵素體?馬氏體雙相，建立了考慮DP600鋼板真實微觀組織的二維有限元模型。在齒輪接觸疲勞壽命研究中，基于真實微觀結構的方法還鮮有報道。

齒輪疲勞失效體現(xiàn)為損傷逐漸累積的動態(tài)演變過程。隨著循環(huán)接觸次數(shù)增加，不斷累積的疲勞損傷造成材料力學屬性衰退，材料屬性演變導致應力?應變場的重新分布和損傷的非均勻化發(fā)展。LEMAITRE 等[19]提出了連續(xù)損傷力學(continuum damage mechanics，CDM)，通過循環(huán)載荷作用下的應力?應變響應定義損傷變量，實現(xiàn)材料本構和損傷累積的耦合，從而通過損傷描述材料的劣化狀態(tài)。RAJE 等[20?21]將損傷變量定義為剪應力幅值的相關函數(shù)，同時考慮材料性能退化研究軸承鋼滾動接觸疲勞，實現(xiàn)了接觸疲勞過程中材料力學屬性劣化演變的模擬；LOPES 等[22]采用Lemaitre連續(xù)損傷模型得到多軸非比例加載條件下的疲勞壽命；LI等[23]以軸承為研究對象，建立了耦合損傷累積和材料本構方程的疲勞分析數(shù)值模型，探究軸承循環(huán)接觸下的疲勞裂紋萌生規(guī)律；ZHAO等[24]考慮水靜應力和八面體剪應力對損傷累積的影響，開發(fā)了接觸疲勞連續(xù)損傷力學模型，并用其預測接觸疲勞裂紋萌生壽命；HE 等[25]綜合考慮齒輪循環(huán)受載過程中的彈性和塑性損傷的累積，基于CDM理論，預測風電齒輪損傷演化規(guī)律以及疲勞裂紋萌生壽命。相比基于經(jīng)驗型的應力/應變?壽命映射關系，基于連續(xù)損傷理論的高周疲勞損傷研究進一步揭示了疲勞失效的物理本質，然而這些研究難以描述材料的微觀結構特征對滾動接觸疲勞壽命的影響。

本文采用機器視覺圖像處理方法構建了齒輪節(jié)點處嚙合區(qū)域材料的真實微觀拓撲結構，結合損傷耦合彈塑性材料本構關系，建立了風力機齒輪接觸疲勞有限元數(shù)值模型，研究齒輪重復嚙合過程中的力學響應，分析不同微觀結構特征對損傷速率的影響。

1 齒輪材料微觀結構表征

圖1所示為風力發(fā)電機機齒輪箱模型。該齒輪副在工程實踐中經(jīng)常發(fā)生早期接觸疲勞失效。齒輪材料為18CrNiMo7-6，其基本參數(shù)如表1所示。在齒輪接觸疲勞模擬中，采用額定輸入轉矩代表齒輪運行的載荷工況條件。

圖1 風力發(fā)電機機齒輪箱模型Fig.1 Wind turbine gearbox model

表1 齒輪副參數(shù)Table 1 Parameters of gear pair

根據(jù)齒輪嚙合理論，齒輪在節(jié)點處嚙合區(qū)域的接觸疲勞失效風險最大。因此，選擇中間級輸出齒輪節(jié)點附近長×寬為0.1 mm×1.0 mm 的區(qū)域作為代表性體積單元(representative volume element，RVE)區(qū)域，生成微觀拓撲結構，如圖2(a)所示。RVE 可以體現(xiàn)材料微觀組織結構的幾何、分布等信息，并能代表材料微觀組織結構的基本特征[26]。采用樹脂熱鑲嵌方法鑲嵌金相試樣，在碳化硅耐水砂紙上加水濕磨，隨后在絨布上采用金剛石噴霧拋光劑拋光，拋光面光潔平整、無劃痕，用無水乙醇清洗、吹干后腐蝕，腐蝕液配比為100 mL熱水+苦味酸(4 g)+洗衣粉(十二烷基磺酸鈉)，腐蝕時間為15～20 min[27]。采用金相顯微鏡拍攝RVE 區(qū)域表征照片。利用圖像拼接方法，得到圖2中所示齒輪表面距心部深度方向為1 mm、沿齒面滾動方向為0.1 mm 的區(qū)域的顯微結構圖，圖2(b)給出了局部區(qū)域微觀結構的放大圖。

圖2 齒輪節(jié)點附近RVE區(qū)域微觀結構Fig.2 Microstructure near gear pitch point within RVE

基于齒輪材料微觀結構表征結果，通過圖3所示的灰度化處理、灰度變換增強、濾波處理、邊緣檢測、形態(tài)學操作、晶粒重構等步驟生成齒輪材料虛擬微觀結構，具體步驟如圖4所示。

圖3 圖像處理分析框圖Fig.3 Flow chart of image processing

圖4 圖像處理過程圖Fig.4 Image processing process

1)讀取齒輪材料微觀結構表征圖像，采用MATLAB 圖像處理工具箱(image processing toolbox，IPT)進行灰度化處理[28]，去除源圖片中冗余信息，減小計算量。

2)進一步對圖4(b)所示圖像進行灰度變化增強處理，調整對比度，采用中值濾波去除圖4(c)中的噪點，以精確識別出單個晶粒。

3)采用Sobel 算子檢測圖4(d)的晶粒邊界，得到的二值圖中白色像素即為檢測出的晶粒邊界。

4)設置面積閾值，將圖4(e)中的雜點顆粒刪除，采用形態(tài)學膨脹處理，填充不連續(xù)晶界間的空隙，避免晶粒幾何模型重構過程中出現(xiàn)錯誤[29]。再經(jīng)細化處理、剪枝運算，獲取單像素(單一寬度)圖像。

5)基于求得的晶粒邊界交點坐標將晶界曲線簡化為直線，并融合短邊[30]。最終建立了如圖4(g)所示的基于真實微觀結構的齒輪材料晶粒重構模型。

2 齒輪材料力學性能測試

在齒輪試樣橫截面上，通過顯微壓痕試驗測量RVE 區(qū)域內材料沿深度方向的維氏硬度，測點分布如圖5(a)所示。試驗采用加載力為9.8 N，保持時間為15 s。為避免試驗中測量誤差的影響，進行了4組顯微壓痕試驗。如圖5(b)所示為滲碳材料維氏硬度隨深度變化的測量結果。由于選取RVE區(qū)域深度范圍較小，在1 mm深度范圍內維氏硬度近似恒定，主要集中在650～690。取測量結果的平均值671.1。研究表明硬化鋼材料屈服強度σy與維氏硬度Hv之間存在線性關系[8]

圖5 齒輪材料力學性能測試Fig.5 Mechanical properties test of gear material

式中：σy為屈服強度，MPa；Hv為維氏硬度。由式(1)計算得屈服強度σy=1839.4 MPa。

3 損傷耦合彈塑性本構模型

3.1 損傷耦合彈塑性本構關系

由于環(huán)境引起的極端載荷條件，加上材料性能的惡化，風力發(fā)電機組齒輪將發(fā)生塑性變形。彈塑性材料本構方程能很好地反映循環(huán)載荷作用下齒輪材料的應力?應變響應。本文基于彈塑性本構方程，結合CDM理論，計算風電機組在循環(huán)接觸疲勞過程中的損傷演化。引入損傷變量D表征材料在循環(huán)載荷作用下力學性能的劣化。只考慮材料的各向同性時，損傷變量D為標量。當D達到臨界值DC=1 時，材料點被完全破壞，該材料點不能承受任何載荷，第1個微裂紋視為在該臨界位置萌生。將Von Mises 屈服準則與連續(xù)損傷理論相結合，考慮損傷變量的屈服函數(shù)f如下：

式中：J2為Von Mises等效應力；k為屈服面的初始半徑；S為偏應力張量；α為表示屈服面中心的背應力張量；符號:表示張量雙點積運算。

總應變ε等于塑性應變εp與彈性應變εe的總和，即

由彈性應變張量εe和疲勞損傷標量D共同決定應力張量σ，其表達式為

式中：C代表4階彈性張量矩陣。

齒輪材料采用與速率無關的塑性模型，塑性應變率計算公式如下：

式中：表示塑性系數(shù)；是塑性流動方向。

塑性系數(shù)計算如下

采用PRAGER[31]提出的線性運動硬化模型，并結合連續(xù)損傷力學，確定背應力速率如下：

式中：M為線性運動硬化模量，齒輪材料18CrNiMo7-6的初始值設為M=5%E。

根據(jù)ZHAO 等[24,32?33]研究結果，總損傷由2 部分組成：依賴于循環(huán)應力狀態(tài)的彈性損傷和由每個疲勞周期累積塑性應變控制的塑性損傷。因此，損傷演化速率為

式中：De和Dp分別代表彈性損傷和塑性損傷；N代表循環(huán)加載次數(shù)。

本研究采用的彈塑性損傷演化規(guī)律遵循WARHADPANDE 等[7]提出的形式。其中，彈性損傷演化方程為

式中：τR和m為材料常數(shù)；ΔτGB是整個循環(huán)加載過程中的晶界上的剪應力變化范圍。

塑性損傷演化方程為

式中：S和q為材料常數(shù)；σM,max為加載循環(huán)中的最大Von Mises應力；為一個加載循環(huán)中沿晶界分解的累積塑性應變率，可定義為2個相鄰加載循環(huán)之間的累積塑性應變增量。

3.2 疲勞損傷參數(shù)確定

為了求解齒輪材料的損傷演化方程，需要確定式(9)～(10)中滲碳齒輪材料的疲勞損傷參數(shù)τR，m，S和q。假設在每個循環(huán)加載期間，應力場和應變場不變。對式(9)和式(10)分別積分，得到ΔτGB與σM,max表達式如下：

式中：Δεp是對應于σM,max的一個加載循環(huán)的循環(huán)應力?應變行為的塑性應變范圍，Δεp設定為0.1。

根據(jù)Basquin 規(guī)則[34]，S?N 曲線的一般方程采用冪律形式，擬合為

式中：A和B分別為疲勞強度系數(shù)和疲勞強度指數(shù)。A和B可以從扭轉S?N 曲線中獲得，根據(jù)LI等[35]提供的扭轉S ?N 曲線，計算結果為A=1 391.5 MPa；B=?0.097。

對于純扭轉情況，ΔτGB和σM,max的關系如下：

由上述方程，計算可得m=10.3；τR=3521.2 MPa；q=5.15；S=25.8 MPa。

模擬風力發(fā)電機齒輪的高周循環(huán)疲勞需要大量循環(huán)次數(shù)，在計算上不可能求解上述各循環(huán)載荷的損傷耦合本構方程。采用ZHAN 等[36]中的“Jump in cycle”方法模擬齒輪高周疲勞，可以降低計算時間。從初始的未損傷狀態(tài)到D=DC的疲勞破壞過程大致由一系列的有限個ΔN循環(huán)(代表加載塊)表示。假設每個加載塊的加載循環(huán)次數(shù)ΔN恒定。在每個循環(huán)加載塊內，RVE 區(qū)域內各單元的應力、應變和損傷率也假定不變。因此，損傷變量的增量ΔD可以表示為

圖6所示為齒輪接觸疲勞數(shù)值處理流程圖，其中i和j分別表示第i個加載塊和RVE區(qū)域內的單元編號數(shù)j。在循環(huán)載荷作用下，疲勞損傷逐漸累積，導致材料力學性能下降。本文考慮了彈性模量E、硬化模量M和初始屈服極限σy不斷退化，以及損傷變量D的更新。材料力學性能的更新如下

圖6 接觸疲勞損傷演變流程圖Fig.6 Calculation flow chart of contact fatigue damage evolution

4 齒輪接觸疲勞有限元模型構建

基于表1中的齒輪基本參數(shù)，采用商業(yè)有限元軟件ABAQUS 6-14，建立如圖7(a)所示中間并聯(lián)級齒輪的幾何模型。采用Python 編程將MATLAB中生成的基于真實齒輪材料的微觀幾何拓撲信息導入到ABAQUS，在中間級輸出齒輪節(jié)點RVE 區(qū)域處重新生成齒輪材料的微觀結構。RVE 區(qū)域內共有1 779個晶粒。連接每個晶粒的質心與其相應的頂點，將齒輪材料晶粒離散為更細的三角形網(wǎng)格，如圖7(b)所示。在RVE 的計算域內采用了非常精細的CPE3類型網(wǎng)格，單元平均尺寸3 μm，如圖7(c)所示。為節(jié)約計算成本，RVE以外區(qū)域采用過度網(wǎng)格，網(wǎng)格尺寸逐漸增大。對輸出齒輪施加210 kN·m 的扭矩(額定負載)和對輸入齒輪施加轉角，模擬齒輪的重復嚙合過程。

圖7 齒輪接觸疲勞有限元模型Fig.7 FE model of gear contact fatigue

利用用戶材料子程序(UMAT)將損傷耦合彈塑性本構方程和損傷演化規(guī)律應用到有限元程序ABAQUS 中，定義齒輪材料力學性能參數(shù)包括彈性模量、泊松比、硬化模量以及屈服極限，并與ABAQUS 進行數(shù)據(jù)交換。對于每個加載塊，可以通過由UMAT 編碼的相關方程計算預定義的狀態(tài)變量，如晶界剪切應力的范圍、最大Mises應力和等效塑性應變增量。重復該模擬程序，直到RVE區(qū)域內晶界單元達到臨界損傷(DC=1)。

5 齒輪接觸疲勞結果分析

圖8所示為接觸壓力分布隨循環(huán)加載次數(shù)的演化。由圖8可見：接觸壓力分布并不嚴格遵循理想的赫茲壓力分布，在接觸疲勞過程中，接觸壓力在靠近接觸中心線的位置不斷波動。在RVE 區(qū)域中，最大接觸壓力的分布形狀及位置幾乎保持不變，但由于疲勞劣化過程中材料的承載能力逐漸降低，導致最大接觸壓力隨著滾動循環(huán)次數(shù)增加而呈現(xiàn)逐漸減小的趨勢。

圖8 接觸壓力分布演化特征Fig.8 Evolution characteristics of pressure distribution

圖9所示為齒輪循環(huán)受載歷程不同階段時RVE區(qū)域內晶界剪應力ΔτGB分布云圖。由圖9可見：ΔτGB最大值主要集中在與齒輪接觸表面距離為0.3～0.5 mm 的范圍內。隨著齒輪循環(huán)加載，在臨界深度處的晶界剪應力大范圍由700 MPa 降至600 MPa。在同一深度范圍內ΔτGB也呈現(xiàn)出差異性，這主要是因為即使在同樣深度范圍內，晶界的角度也存在著較大差異。

圖9 不同接觸疲勞周期晶界剪應力ΔτGB分布情況Fig.9 Distributions of ΔτGB recorded at various number of loading cycles

為了進一步探討晶界角度θ與晶界剪應力ΔτGB的關系，本文選取關鍵區(qū)域0.3～0.5 mm 的深度范圍晶粒，提取該區(qū)域每個晶粒晶界法線方向與相應的ΔτGB數(shù)據(jù)，繪制了接觸疲勞進程中不同階段的晶界角度與晶界剪應力變化圖，如圖10所示。圖10中每個點代表1個選取關鍵區(qū)域的晶界。由圖10可以看出，在齒輪滾動接觸疲勞前期，晶界剪應力ΔτGB最大值出現(xiàn)在晶界角度θ為20°或110°附近。隨著齒輪接觸疲勞損傷不斷累積和材料力學性能退化，在θ為20°或110°附近，晶界剪應力ΔτGB由最初的最大值730 MPa 逐漸衰減至最小值450 MPa。

圖10 晶界角度θ與晶界剪應力ΔτGB隨循環(huán)加載次數(shù)的演化關系Fig.10 Interaction between θ and ΔτGB at various stages of fatigue process

圖11所示為不同接觸疲勞周期晶界損傷D分布情況。晶界損傷程度較高區(qū)域同樣主要集中在0.3～0.5 mm 深度范圍內，與晶界上剪應力變化范圍峰值的深度范圍相對應。隨著循環(huán)加載次數(shù)不斷增加，晶界的疲勞損傷逐漸累積，損傷開始出現(xiàn)局部化現(xiàn)象，并最終集中于個別晶界單元，在經(jīng)歷7×106次加載循環(huán)后，該單元的累積損傷達到臨界值，疲勞裂紋是在該位置處萌生。

圖11 不同接觸疲勞周期晶界損傷D分布情況Fig.11 Distributions of damage D recorded at various number of loading cycles

選取相同RVE區(qū)域，采用Voronoi剖分技術生成齒輪材料微觀結構，與機器視覺圖像處理技術識別的齒輪材料真實微觀結構進行對比研究。圖12所示為不同微觀結構特征下?lián)p傷速率隨加載循環(huán)次數(shù)演化情況。由圖12可見：對不同微觀結構，特征損傷速率的演變趨勢總體相同，隨著加載循環(huán)次數(shù)增加，損傷速率逐漸加快；由于微觀結構特征的差異性，在損傷前期，采用Voronoi 剖分法時損傷速率較圖像處理法略快，當接觸疲勞損傷不斷累積，采用圖像處理法的損傷速率顯著加快，而Voronoi剖分法的損傷速率增加趨勢較緩慢。

圖12 不同微觀結構特征下?lián)p傷速率隨加載循環(huán)次數(shù)演化情況Fig.12 Evolution of damage rate with increasing number of loading cycles under different microstructure characteristics

6 結論

1)最大接觸壓力隨著滾動循環(huán)次數(shù)增加而呈現(xiàn)逐漸減小的趨勢，并且在靠近接觸中心線的位置處，接觸壓力不斷波動。

2)晶界上剪應力峰值與晶界損傷程度較高區(qū)域均集中在0.3～0.5 mm范圍內。晶界的疲勞損傷隨著循環(huán)加載次數(shù)增加不斷累積，損傷逐漸出現(xiàn)局部化現(xiàn)象，最終集中于個別晶界單元。晶界角度θ為20°或110°時，晶界剪應力ΔτGB隨齒輪嚙合過程由最大值730 MPa逐漸衰減至最小值450 MPa。

3)對比機器視覺圖像處理技術與Voronoi 剖分技術構建的齒輪材料微觀結構，發(fā)現(xiàn)損傷速率的演變趨勢總體相同，但是兩者的損傷速率在損傷的不同階段表現(xiàn)出差異性，在損傷后期，前者損傷速率相比后者顯著加快。采用圖像處理法考慮了材料微觀結構內部的真實晶粒分布特征，更符合實際的損傷演化情況。