微動磨損工況下二維柱面對平面磨損形貌的數值計算與分析

唐正強，張慧杰，李 儼，張太華

（1.貴州大學機械工程學院，貴州 貴陽 550025；2.貴州大學現代制造技術教育部重點實驗室，貴州 貴陽 550025；3.貴州師范大學機械與電氣工程學院，貴州 貴陽 550001）

1 引言

微動是發生在兩接觸表面之間微、納米量級的運動[1]。這種微小幅度的振動容易被人們忽視，卻是機械零構件失效的主要原因[2]。因此，微動被認為是現代工業機械部件的癌癥[3]。

眾多學者對微動磨損做了大量的實驗研究，研究微動工況、材料屬性等因素對磨損的影響[4-8]。文獻[4]研究了304 不銹鋼在不同位移振幅下的微動磨損行為，結果表明304 不銹鋼的磨損系數和微動損傷均隨微動幅值的增大而增大。文獻[5]在鋯合金包殼的微動疲勞試驗中研究了滑移狀態下的磨損形貌、磨損體積和磨損深度等特性。文獻[6]研究了熱處理對316L 不銹鋼表面TiO2涂層磨損的影響，結果表明熱處理增強了涂層的抗磨損能力。文獻[7]通過掃面電鏡觀察了不同接觸載荷下鋼絲的磨損形貌，結果表明隨著接觸載荷的增加，鋼絲的穩定摩擦系數呈現降低趨勢。文獻[8]研究了偏心載荷作用下螺栓連接的松動行為，發現螺紋接觸表面的微動磨損將會導致螺栓受到的軸向力下降。微動磨損對構件的影響往往在高周循環載荷之后才能體現出來，有的構件甚至要在微動工況下服役幾年的時間。因此，微動磨損實驗無論是從時間還是經濟上都是很大的浪費，數值模擬則彌補了實驗的這些缺點。

近年來，基于Archard 磨損方程的微動磨損數值模擬研究受到眾多學者的關注[9-14]。文獻[9]是較早應用Archard 磨損方程來模擬二維柱面與平面磨損的學者，通過對比試驗研究驗證了其數值模擬的可靠性。文獻[10]推導了基于Archard 方程的局部磨損方程，并結合有限元法預測了材料表面磨損形貌。文獻[11]利用數值模擬方法對考慮累積塑性變形的平面-平面接觸問題進行了模擬研究，結果表明塑性變形量變大，最大接觸應力變小。文獻[12]建立了基于有限元法的二維平面-柱面微動磨損模型，并通過多層節點更新方法解決了磨損深度大于有限單元邊長的問題。文獻[13]基于二維柱面-平面研究了整體滑移與部分滑移對鋯合金磨損速率的影響，發現磨損速率隨應變載荷振幅增加而增加，磨損率在不同接觸狀態隨法向力變化的規律不同。文獻[14]利用有限元法研究了髖關節假體頭頸錐度連接處的微動磨損，發現錐角的類型和大小對材料損有很大影響。

從上述研究中可知，在微動磨損實驗研究進展的同時，越來越多的學者采用數值模擬的方法研究微動磨損。一方面可以節約實驗所需的時間和成本，另一方面可以通過接觸體磨損過程中的彈塑性變形、非對稱性應力應變等變化深入探索磨損機理，提出有效抑制磨損的方法。針對微動磨損工況下接觸表面磨損形貌問題，基于Archard 磨損方程推導了適應于任意微動振幅的局部磨損公式，建立了二維剛性柱面與彈性平面的微動磨損模型，研究了二維柱面半徑、微動載荷與柱面彈性模量對彈性平面磨損深度、磨損寬度和磨損面積的影響。

2 微動磨損模型與數值計算方法

2.1 二維柱面-平面有限元模型

為研究在接觸過程中不同參數對接觸副磨損輪廓的影響，根據三維圓柱與平板接觸模型，如圖1（a）所示。建立二維柱面與平面的有限元接觸模型，如圖1（b）所示。圓柱面半徑為R，長度為L，平板的厚度為t。模型通過文獻[15]進行網格劃分，不同的區域采用了不同的網格大小。區域Ⅰ和區域Ⅳ是柱面與平面的接觸區域，使用了邊長為0.5mm 的單元；遠離接觸的區域Ⅲ和區域Ⅵ，使用了邊長為2mm 的單元；區域Ⅱ與區域Ⅴ為網格的過渡區域，使用邊長為1mm 的單元。為了研究剛性柱面對彈性平面的微動磨損影響，剛性柱面按照文獻[16]中“理想剛體”的定義，設置二維柱面材料的彈性模量為平面材料的1000 倍。

圖1 微動磨損三維模型與二維有限元模型Fig.1 Three-Dimensional Model of Fretting Wear and FE Model

二維柱面與平面的接觸計算采用罰函數法[17]，平面通過約束底端全部節點x 和y 方向的自由度而固定，圓柱面頂端施加載荷P/L（P 為正壓力載荷，L 為圓柱的長度），同時圓柱面頂端節點上施加了平行于x 方向的往復式位移約束，如圖2 所示。

圖2 正壓力載荷與切向位移Fig.2 Applied Normal Load and Tangential Displacement

2.2 磨損的計算方法

為計算接觸平面的磨損輪廓，基于Archard 磨損方程，文獻[10]給出了式（1）來計算磨損深度：

式中：δ—設定微動振幅的幅值；n—每個微動循環子步數。微動過程中二維柱面和平面的接觸示意圖，如圖3 所示。

圖3 時間步j 到j+1 的柱面/平面的接觸示意圖Fig.3 Contact at the Cylinder-Flat Interface for Time Step j and j+1

每一個接觸節點在不同的時間步，其接觸壓強不同。從圖中可以看出，當時間步從j 增加到j+1 時，節點i 的接觸壓強從Pi，j變到Pi，j+1（其中i 為節點號，j 為時間步），則根據式（1），可得到節點i 從時間步j 到時間步j+1 的磨損深度為：

微動磨損數值模擬的具體計算流程，如圖4 所示。（1）在LSDyna 有限元軟件中建立二維柱面與平面的接觸模型；（2）將有限元模型保存為LS-Dyna 求解器可識別的K 文件，便于后續模型的更改；（3）調用LS-Dyna 的求解器求解K 文件，計算接觸壓強；（4）通過MATLAB 讀取結果文件，獲得接觸表面節點的坐標、接觸壓強等信息；（5）通過磨損深度計算式（6），計算在特定循環增量步內接觸節點的磨損深度，獲得磨損輪廓；（6）判斷微動循環次數是否達到總的循環次數，若未達到，則更新二維柱面與平面的幾何模型，獲得更新后的K 文件，若達到預設的循環次數，則輸出結果，計算結束。

圖4 微動磨損模擬流程圖Fig.4 Flowchart of the Fretting Wear Procedure

3 結果討論與分析

3.1 接觸壓強

為了驗證這里有限元模型的正確性，比較了Hertz 接觸理論與有限元模擬結果。Hertz 接觸對于二維的柱面-平面的接觸問題給出了接觸面正壓力分布的解析表達式[18]：

式中：p0—接觸區域的最大壓強；a—接觸區域長度；E*—接觸副材料的等效彈性模量；R—等效接觸體半徑；Ec、Ef—接觸副材料的彈性模量；vc、vf—接觸副材料的泊松比；Rc、Rf—接觸副的曲率半徑；“c”、“f”—柱面和平面。本模型中壓力載荷設為20MPa；Ec設置為200000GPa、Ef設置為200GPa；vc、vf均為0.3。將有限元計算得到的壓強與Hertz 接觸理論計算得到的結果進行對比，結果如圖5 所示。

圖5 接觸壓強分布的有限元結果與Hertz 理論比較Fig.5 Comparison between Contact Pressure from FEM Analysis and Analytical Solution of Hertz

從圖5 中可以看出，有限元結果和Hertz 接觸理論計算結果趨于一致，但有限元結果的接觸寬度值略大于Hertz 接觸理論，這是由兩接觸面接觸區域的網格單元大小引起的誤差。綜合考慮計算精度和磨損模型的計算時間，采取了該網格尺寸。

在不同工況下平面的壓強分布，如圖6 所示。從圖6（a）中可以看出，隨著柱面半徑的增大，最大接觸壓強降低，接觸長度增加；這是因為隨著柱面半徑的增大，圓柱面和平面的等效接觸半徑增大，導致最大接觸壓強降低；從圖6（b）中看出接觸壓強也隨著載荷的增大而增大；圖6（c）顯示了不同彈性模量下平面的接觸壓強變化，隨著平面彈性模量的增加，由式（10）可知等效彈性模量增大，再由式（9）可知接觸區域的最大壓強增大。

圖6 不同工況下平面的接觸壓強Fig.6 Effect of the Fretting Wear Conditions on Contact Pressure of the Flat

3.2 圓柱半徑對磨損的影響

磨損發生在接觸體表面，多數學者研究了接觸副中被動磨損件結構參數對磨損輪廓的影響，而主動磨損件結構參數對磨損輪廓的研究鮮有報道。通過改變二維柱面半徑研究二維柱面幾何尺寸對平面磨損輪廓的影響。微動工況下（壓力為20MPa，微動振幅為4μm，微動循環次數為1000 次）圓柱半徑對平面磨損形貌的影響，從圖中可發現，隨著圓柱半徑的增加，平面的磨損形貌呈現出變寬變淺的趨勢，如圖7（a）所示。

圖7 描述了不同柱面半徑下平面的磨損深度（圖7（b））、磨損寬度（圖7（c））以及磨損面積（圖7（d））隨微動循環次數的變化規律。磨損面積是指平面的磨損輪廓與原始輪廓圍成的面積。從圖中可以看出平面的磨損深度、磨損寬度和磨損面積隨著微動循環次數的增加而增加。其次，我們發現平面磨損深度隨著圓柱半徑的增大而減小，磨損寬度隨著圓柱半徑的增大而增大。由式（11）可知，假設平面的半徑Rf趨于無窮大，隨著圓柱半徑Rc的增大，等效接觸半徑R 增大，導致磨損區域接觸長度a 增大。同時，在載荷不變的情況下導致接觸壓強減小，因此，平面的磨損深度隨著Rc的增大而減小。其次，我們發現平面的磨損面積隨著微動循環次數的增加而呈線性增加，另外，圓柱半徑對平面的磨損面積沒有影響。這是因為隨著圓柱半徑的增加，平面磨損深度減小而磨損寬度增加共同作用的結果。

圖7 柱面半徑大小對平面磨損的影響Fig.7 Effect of the Cylindricalsurface Radius on Wear Performance of the Flat

3.3 載荷對磨損的影響

在1000 次微動循環次數下平面的磨損形貌，我們可以看出隨著載荷的增加，平面的磨損輪廓有顯著的區別，表明載荷對磨損形貌的變化有著顯著影響，如圖8（a）所示。

圖8（b）～圖8（d）描述了不同載荷作用下平面的磨損深度（圖8（b））、磨損寬度（圖8（c））和磨損面積（圖8（d））隨微動循環次數的變化規律。從圖中可以看出，隨著微動循環次數的增加，平面的磨損深度、磨損寬度和磨損面積均逐漸增加；同時隨著載荷的增加，平面的磨損深度、磨損寬度和磨損面積也大幅增加。

圖8 載荷對平面磨損的影響Fig.8 Effect of the Normal Load on Wear Performance of the Flat

3.4 彈性模量對磨損的影響

1000 次微動循環后的平面磨損輪廓，發現在不同的彈性模量下，平面的磨損輪廓改變較大，如圖9（a）所示。隨著彈性模量的增加，平面磨損輪廓從寬而淺變化到深而窄。

圖9（b）～圖9（d）顯示了不同彈性模量下平面的磨損深度（圖9（b））、磨損寬度（圖9（c））和磨損面積（圖9（d））隨著微動循環次數的變化規律。隨著微動循環次數的增加，平面的磨損深度、磨損寬度和磨損面積均逐漸增加；隨著彈性模量的增加，平面的磨損深度增加但磨損寬度減小。從圖9（d）中看出，雖然彈性模量對平面的磨損深度和磨損寬度呈現出相反的影響效果，但是隨彈性模量的增加平面的磨損面積增大。此外，我們發現不同彈性模量下磨損面積與微動循環次數均呈現出線性關系。

圖9 彈性模量對平面磨損的影響Fig.9 Effect of the Elastic Modulus on Wear Performance of the Flat

這里的模型基于有限元接觸壓強的計算和Archard 磨損公式的推導，通過迭代更新的方法計算磨損形貌。在初始接觸階段，有限元的接觸壓強計算結果應與Hertz 接觸理論保持一致。因此，我們試想直接通過Hertz 接觸理論和Archard 磨損方程來預測平面磨損面積隨著微動循環次數的變化規律，與模型計算結果進行對比，如圖10 所示。選取了圓柱半徑Rc=5mm，載荷為20MPa，微動振幅為4μm 的工況采用兩種計算方法進行磨損面積預測。從圖中可以看出，在微動循環次數較低時，根據Hertz 接觸理論和Archard 磨損方程預測的平面磨損面積與本論文中模型計算的結果差別不大，但隨著微動循環次數的增加，Hertz 理論預測結果小于這里模型計算的結果。這是因為在模型計算過程中，隨著微動循環次數的增加，彈性平面的磨損形貌實時更新，二維柱面與磨損平面的接觸不再符合經典的Hertz 接觸理論，加上微動振幅的作用，使得這里模型的接觸半徑大于Hertz 理論的接觸半徑，因而在磨損面積計算時這里模型計算的磨損面積值略大于Hertz 理論預測值。

圖10 磨損面積預測結果對比Fig.10 Wear Area Comparison between Hertz Theory and the Present Result

這里的模型忽略了材料的界面環境和磨粒（第三體）的參與，因此，隨著微動循環次數的增加，磨損面積呈線性增加。雖然該模型的假設與真實環境中的磨損現象有差別，但能通過模型快速預測材料表面的磨損趨勢，合理的選擇材料參數，設計較合理的微動工況，抑制或者減緩材料的磨損，增加零部件的使用壽命。

4 結論

基于改進的Archard 磨損方程，推導了適應于不同微動振幅下的局部磨損方程，建立了二維剛性柱面-平面的數值模型。分別從平面的磨損深度、磨損寬度和磨損面積三個特征量研究了二維柱面對平面磨損形貌的影響。研究結果表明隨著柱面半徑的增大，平面的磨損深度降低、磨損寬度增大、磨損面積保持不變；隨著載荷的增大，磨損深度、磨損寬度、磨損面積均呈現增大趨勢；而隨著平面彈性模量的增大，平面的磨損深度增大、磨損寬度降低，但磨損面積增大。該模型能快速預測材料表面的磨損趨勢，有助于設計較合理的微動工況，抑制或者減緩材料的磨損。