微動(dòng)磨損工況下二維柱面對(duì)平面磨損形貌的數(shù)值計(jì)算與分析

唐正強(qiáng)，張慧杰，李 儼，張?zhí)A

（1.貴州大學(xué)機(jī)械工程學(xué)院，貴州 貴陽(yáng) 550025；2.貴州大學(xué)現(xiàn)代制造技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，貴州 貴陽(yáng) 550025；3.貴州師范大學(xué)機(jī)械與電氣工程學(xué)院，貴州 貴陽(yáng) 550001）

1 引言

微動(dòng)是發(fā)生在兩接觸表面之間微、納米量級(jí)的運(yùn)動(dòng)[1]。這種微小幅度的振動(dòng)容易被人們忽視，卻是機(jī)械零構(gòu)件失效的主要原因[2]。因此，微動(dòng)被認(rèn)為是現(xiàn)代工業(yè)機(jī)械部件的癌癥[3]。

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

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

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

2 微動(dòng)磨損模型與數(shù)值計(jì)算方法

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

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

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

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

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

2.2 磨損的計(jì)算方法

為計(jì)算接觸平面的磨損輪廓，基于Archard 磨損方程，文獻(xiàn)[10]給出了式（1）來(lái)計(jì)算磨損深度：

式中：δ—設(shè)定微動(dòng)振幅的幅值；n—每個(gè)微動(dòng)循環(huán)子步數(shù)。微動(dòng)過(guò)程中二維柱面和平面的接觸示意圖，如圖3 所示。

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

每一個(gè)接觸節(jié)點(diǎn)在不同的時(shí)間步，其接觸壓強(qiáng)不同。從圖中可以看出，當(dāng)時(shí)間步從j 增加到j(luò)+1 時(shí)，節(jié)點(diǎn)i 的接觸壓強(qiáng)從Pi，j變到Pi，j+1（其中i 為節(jié)點(diǎn)號(hào)，j 為時(shí)間步），則根據(jù)式（1），可得到節(jié)點(diǎn)i 從時(shí)間步j(luò) 到時(shí)間步j(luò)+1 的磨損深度為：

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

圖4 微動(dòng)磨損模擬流程圖Fig.4 Flowchart of the Fretting Wear Procedure

3 結(jié)果討論與分析

3.1 接觸壓強(qiáng)

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

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

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

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

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

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

3.2 圓柱半徑對(duì)磨損的影響

磨損發(fā)生在接觸體表面，多數(shù)學(xué)者研究了接觸副中被動(dòng)磨損件結(jié)構(gòu)參數(shù)對(duì)磨損輪廓的影響，而主動(dòng)磨損件結(jié)構(gòu)參數(shù)對(duì)磨損輪廓的研究鮮有報(bào)道。通過(guò)改變二維柱面半徑研究二維柱面幾何尺寸對(duì)平面磨損輪廓的影響。微動(dòng)工況下（壓力為20MPa，微動(dòng)振幅為4μm，微動(dòng)循環(huán)次數(shù)為1000 次）圓柱半徑對(duì)平面磨損形貌的影響，從圖中可發(fā)現(xiàn)，隨著圓柱半徑的增加，平面的磨損形貌呈現(xiàn)出變寬變淺的趨勢(shì)，如圖7（a）所示。

圖7 描述了不同柱面半徑下平面的磨損深度（圖7（b））、磨損寬度（圖7（c））以及磨損面積（圖7（d））隨微動(dòng)循環(huán)次數(shù)的變化規(guī)律。磨損面積是指平面的磨損輪廓與原始輪廓圍成的面積。從圖中可以看出平面的磨損深度、磨損寬度和磨損面積隨著微動(dòng)循環(huán)次數(shù)的增加而增加。其次，我們發(fā)現(xiàn)平面磨損深度隨著圓柱半徑的增大而減小，磨損寬度隨著圓柱半徑的增大而增大。由式（11）可知，假設(shè)平面的半徑Rf趨于無(wú)窮大，隨著圓柱半徑Rc的增大，等效接觸半徑R 增大，導(dǎo)致磨損區(qū)域接觸長(zhǎng)度a 增大。同時(shí)，在載荷不變的情況下導(dǎo)致接觸壓強(qiáng)減小，因此，平面的磨損深度隨著Rc的增大而減小。其次，我們發(fā)現(xiàn)平面的磨損面積隨著微動(dòng)循環(huán)次數(shù)的增加而呈線(xiàn)性增加，另外，圓柱半徑對(duì)平面的磨損面積沒(méi)有影響。這是因?yàn)殡S著圓柱半徑的增加，平面磨損深度減小而磨損寬度增加共同作用的結(jié)果。

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

3.3 載荷對(duì)磨損的影響

在1000 次微動(dòng)循環(huán)次數(shù)下平面的磨損形貌，我們可以看出隨著載荷的增加，平面的磨損輪廓有顯著的區(qū)別，表明載荷對(duì)磨損形貌的變化有著顯著影響，如圖8（a）所示。

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

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

3.4 彈性模量對(duì)磨損的影響

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

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

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

這里的模型基于有限元接觸壓強(qiáng)的計(jì)算和Archard 磨損公式的推導(dǎo)，通過(guò)迭代更新的方法計(jì)算磨損形貌。在初始接觸階段，有限元的接觸壓強(qiáng)計(jì)算結(jié)果應(yīng)與Hertz 接觸理論保持一致。因此，我們?cè)囅胫苯油ㄟ^(guò)Hertz 接觸理論和Archard 磨損方程來(lái)預(yù)測(cè)平面磨損面積隨著微動(dòng)循環(huán)次數(shù)的變化規(guī)律，與模型計(jì)算結(jié)果進(jìn)行對(duì)比，如圖10 所示。選取了圓柱半徑Rc=5mm，載荷為20MPa，微動(dòng)振幅為4μm 的工況采用兩種計(jì)算方法進(jìn)行磨損面積預(yù)測(cè)。從圖中可以看出，在微動(dòng)循環(huán)次數(shù)較低時(shí)，根據(jù)Hertz 接觸理論和Archard 磨損方程預(yù)測(cè)的平面磨損面積與本論文中模型計(jì)算的結(jié)果差別不大，但隨著微動(dòng)循環(huán)次數(shù)的增加，Hertz 理論預(yù)測(cè)結(jié)果小于這里模型計(jì)算的結(jié)果。這是因?yàn)樵谀Ｐ陀?jì)算過(guò)程中，隨著微動(dòng)循環(huán)次數(shù)的增加，彈性平面的磨損形貌實(shí)時(shí)更新，二維柱面與磨損平面的接觸不再符合經(jīng)典的Hertz 接觸理論，加上微動(dòng)振幅的作用，使得這里模型的接觸半徑大于Hertz 理論的接觸半徑，因而在磨損面積計(jì)算時(shí)這里模型計(jì)算的磨損面積值略大于Hertz 理論預(yù)測(cè)值。

圖10 磨損面積預(yù)測(cè)結(jié)果對(duì)比Fig.10 Wear Area Comparison between Hertz Theory and the Present Result

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

4 結(jié)論

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