摩擦系數(shù)對高副接觸磨損影響的有限元模擬

張優(yōu)云，陳 朝，朱永生，蔣翔俊

(西安交通大學(xué)現(xiàn)代設(shè)計與軸承轉(zhuǎn)子系統(tǒng)教育部重點實驗室，710049西安)

磨損占機械零件失效的60%～80%［1］，因此，研究磨損具有非常重要的意義.隨著計算機技術(shù)的發(fā)展，數(shù)值仿真技術(shù)已經(jīng)成為磨損研究的一種非常重要的手段.蔣翔俊等［2］利用數(shù)值仿真方法研究了不同交變載荷力作用下燃機端面齒的磨損規(guī)律.在數(shù)值仿真研究方法中，有限元仿真分析已經(jīng)在工程中得到了廣泛的應(yīng)用，并且在磨損分析方面具有很強的實用性［3］.馮偉［4］利用有限元方法進行滑動磨損試驗的摩擦學(xué)仿真研究.Sfantos等［5］、Andersson等［6］提出了應(yīng)用邊界元法及有限元方法的磨損預(yù)測方法.Rezaei等［7］利用有限元方法和Archard方程，建立了徑向滑動軸承磨損過程的模型.Cruzado等［8］利用商用有限元軟件ABAQUS，研究了鋼絲繩的微動磨損，并且進行了試驗驗證.宿月文［9］利用有限元方法進行了滑動磨損過程的仿真研究，對銷軸的磨損進行了預(yù)測.但是上述研究只能對長距離的滑動磨損進行仿真，而不能對小幅振動下的磨損進行有效的仿真.摩擦系數(shù)、接觸壓力和滑移量等特征參數(shù)是影響磨損的主要因素［10］.在磨損過程中，摩擦系數(shù)并不是一個恒定的值，而是隨著磨損循環(huán)次數(shù)的增加而逐漸增大，最終在一個較小的范圍內(nèi)波動［11］.在目前的磨損有限元仿真中，針對變化的摩擦系數(shù)的影響研究并不多見.

本文建立了一種基于有限元方法的磨損預(yù)測模型，利用有限元軟件ANSYS12.1對圓柱和平板相接觸的小幅振動磨損過程進行了仿真分析，將仿真結(jié)果與文獻［11］中的實驗數(shù)據(jù)進行了對比驗證，并利用該模型分析了摩擦系數(shù)的變化對磨損量、接觸壓力、滑移量、正應(yīng)力、剪應(yīng)力的影響.

1 磨損分析的有限元建模

1.1 接觸分析

圓柱和平板的接觸問題廣泛應(yīng)用于磨損的研究中，其接觸情況如圖1所示.根據(jù)赫茲接觸計算公式，可以得到接觸壓力的分布為

式中，接觸區(qū)域半寬度a為

最大接觸壓力p0為

其中

式中，E1和E2、ν1和ν2、R1和R2分別為兩個接觸體的彈性模量、泊松比和曲率半徑.

圖1 圓柱和平板的接觸

在有切向力T(t)的作用下，兩個接觸表面間會出現(xiàn)相對滑動.在整體滑移的情況下，切向力的最大值Tmax和摩擦力μp相等(μ為摩擦系數(shù)).如果Tmax＜μp，則會產(chǎn)生局部滑移.此時，接觸區(qū)域的中心會產(chǎn)生一個半寬度為c的粘著區(qū)，而在c≤｜x｜≤a的區(qū)域則會產(chǎn)生一個局部滑移區(qū).瞬時時刻的粘著區(qū)的寬度c*可以表示為

1.2 有限元模型

本文針對的材料與文獻［11］相同，為高強度合金鋼.由于圓柱和平板的接觸屬于線接觸，因此可以將接觸模型簡化成二維模型.圖2是磨損仿真的有限元模型.

圖2 圓柱和平板接觸的有限元模型

半圓和矩形的單元類型均為PLANE42，彈性模量為200 GPa，泊松比為0.3，摩擦系數(shù)為0.8.接觸區(qū)域的半寬度設(shè)置為0.3 mm，這一寬度是根據(jù)赫茲接觸計算所確定的.為了得到比較精確的計算結(jié)果并盡可能的減少節(jié)點數(shù)量，劃分網(wǎng)格時只將接觸區(qū)域的網(wǎng)格細化(網(wǎng)格寬度為5 μm).為了使網(wǎng)格劃分較細的區(qū)域和網(wǎng)格劃分較粗的區(qū)域能夠平穩(wěn)的過渡，該模型在兩個區(qū)域之間增加了一個過渡區(qū).節(jié)點總數(shù)為24 004個.本模型使用接觸單元CONTA172和目標(biāo)單元TARGE169來定義接觸對.但是ANSYS只能保存接觸單元的計算結(jié)果.為了能將相互接觸的兩個面的數(shù)據(jù)都提取出來，該模型定義了兩個接觸對，即以接觸區(qū)域中上方的圓弧為接觸面、下方的直線為目標(biāo)面的接觸對，和以直線為接觸面、圓弧為目標(biāo)面的接觸對.圓柱在平板上的運動方式與文獻［11］中的實驗相同，即以坐標(biāo)原點為中心，沿x方向作往復(fù)運動，幅值為25 μm.

1.3 磨損計算模型

Archard磨損模型是一種經(jīng)典且應(yīng)用非常廣泛的磨損計算模型.該模型可表述如下:

式中:V為總的體積磨損量，m3;P為法相載荷，單位為N;S為總的滑動距離，m;H為材料的布氏硬度，Pa;K為磨損系數(shù)，是一個無量綱的比例常數(shù).為了體現(xiàn)兩接觸表面的輪廓隨磨損循環(huán)進行的變化，研究磨損深度更具實際意義.假設(shè)在一個極小的接觸區(qū)域dA內(nèi)，有一段較小的滑動距離dS，此時應(yīng)用式(1)可以得到

方程兩邊同時除以dA可以得到

式中，dP/dA為局部的接觸壓力，它是一個接觸點水平位置的函數(shù)，記為p(x);dV/dA為局部磨損深度，記為dh;將K/H記為k1，則式(2)可寫為

如果將磨損過程離散為若干個子步來計算磨損量，則第i個節(jié)點在第j次磨損循環(huán)時的磨損增量為式中，k1與材料屬性、運行工況等因素有關(guān);Δsi，j是接觸區(qū)域的滑移量，通常小于實際的滑動距離;接觸壓力pi，j可通過數(shù)值方法，如有限元方法來得到.

1.4 磨損仿真過程

磨損仿真的過程是十分緩慢的，如果完全仿真實際磨損過程的每一步，將耗費大量的資源和時間，而計算精度卻不會有太大的提高.本文采用的磨損計算流程見圖3.為了提高磨損計算的效率，本文主要采用兩種方法:其一，如果將半圓往復(fù)運動1次視為1個磨損循環(huán)，并且假定n(本文的模型中n=20)次磨損循環(huán)內(nèi)的磨損量都是相等的，這樣n次磨損循環(huán)的磨損量就可以用1次磨損循環(huán)的磨損量的n倍得到;其二，1個磨損循環(huán)內(nèi)半圓滑動的距離是原點到右極限位置距離的4倍，假設(shè)1個磨損循環(huán)內(nèi)產(chǎn)生的磨損量是半圓動到右極限位置所產(chǎn)生的磨損量的4倍，這樣也可以在保證計算精度的情況下提高計算效率.在施加y方向向下的載荷時，如果直接將一個集中力施加在A點位置的節(jié)點上，則會導(dǎo)致受集中力作用的區(qū)域產(chǎn)生變形，且受力情況與實際不符.為了避免這一情況，本文用CERGID命令將A點處的節(jié)點和y坐標(biāo)值與其相等的節(jié)點耦合，形成一個剛性區(qū)域，這樣在半圓圓心的節(jié)點處施加集中力載荷時，就是通過整個剛性區(qū)域來傳遞載荷了.平板的底部和兩邊都添加了x和y方向上的位移約束，使之在仿真過程中固定不動.

圖3 磨損計算流程圖

在添加好約束，并在A點施加了y方向向下的載荷后，進入后處理器，將接觸區(qū)節(jié)點的數(shù)據(jù)提取出來.提取的數(shù)據(jù)主要包括節(jié)點的x坐標(biāo)、y坐標(biāo)、壓力和滑移量.然后再在A點對圓柱施加x方向向右的位移載荷.為了提高計算結(jié)果的精度，這里不是直接將半圓移動到右極限位置，而是將這個過程分成若干個子步，每進行一個子步之后，提取節(jié)點的壓力和滑移量，根據(jù)Archard方程計算這一子步的磨損量.在所有的子步都進行完之后，將計算的磨損量疊加，就可以得到將圓柱從原點位置移動到右極限位置這一過程的磨損量.如前所述，將這一結(jié)果乘以4就可以得到一個磨損循環(huán)內(nèi)的磨損量.再乘以n，即可得到n次磨損循環(huán)的磨損量.然后用節(jié)點在經(jīng)歷磨損循環(huán)前的y坐標(biāo)值減去其對應(yīng)的磨損量，得到其磨損后的位置，再將節(jié)點移動到各自的新位置，就可以得到磨損后的表面輪廓.最后卸掉所有的載荷.使用循環(huán)命令反復(fù)進行上述求解操作，就可以得到指定循環(huán)次數(shù)后的磨損量和磨損后的表面輪廓.

2 有限元仿真結(jié)果

圖4是在185 N的載荷下，本文的有限元預(yù)測模型與赫茲接觸理論分別計算出的壓力分布.從圖中可以看出，這兩種計算方法所得到的壓力分布很接近，最大壓力值的誤差為4.67%.

圖4 有限元模型與赫茲接觸分別計算的壓力分布(載荷:185 N)

圖5是有限元預(yù)測模型在載荷為185 N，往復(fù)運動的幅值為25 μm的條件下，不同磨損循環(huán)次數(shù)后，圓柱和平板的表面輪廓.從圖中可以看出，隨著磨損次數(shù)的增加，兩接觸表面的磨損深度逐漸增加，接觸區(qū)域逐漸增大，并且兩個表面的接觸區(qū)的寬度非常接近.

圖5 不同磨損循環(huán)次數(shù)后的表面輪廓(載荷:185 N，幅值:25 μm)

圖6是在載荷為185 N，往復(fù)運動幅值為25 μm的條件下，文獻［11］的實驗數(shù)據(jù)與有限元磨損模型的計算結(jié)果對比圖.從圖中可以看出，在載荷為185 N的條件下，實驗結(jié)果和有限元模型計算的結(jié)果非常接近.

圖6 參考文獻實驗數(shù)據(jù)與本文的有限元磨損模型計算結(jié)果對比(載荷:185 N，磨損次數(shù):18 000次)

圖7是文獻［11］的實驗所測得的摩擦系數(shù).從圖中可以看出，在磨損過程中，摩擦系數(shù)并不是一個恒定的值，而是由初期一個較小的值逐漸增大并最終在一個小的范圍內(nèi)波動.因此，本文就不同的摩擦系數(shù)對磨損的影響進行了討論.

圖7 摩擦系數(shù)隨磨損次數(shù)的變化(往復(fù)運動幅值:25 μm，頻率:20 Hz)

圖8是有限元預(yù)測模型在載荷為500 N，往復(fù)運動幅值為25 μm的條件下，在不同的摩擦系數(shù)下磨損5 000次后圓柱的表面輪廓.從圖中可以看出，隨著摩擦系數(shù)的增大，磨損量呈現(xiàn)出減小的趨勢.根據(jù)本文的Archard模型，磨損量是接觸壓力和滑移量的函數(shù).從圖9中可以看出，摩擦系數(shù)對接觸壓力的影響并不是很明顯，但放大后仍可發(fā)現(xiàn)接觸壓力隨摩擦系數(shù)的增大而減小.而隨著摩擦系數(shù)增大，滑移量亦隨之減小(如圖10所示)，這都是導(dǎo)致磨損量隨摩擦系數(shù)增大而減小的原因.從圖10還可以看出，在500 N的載荷下的滑移均為整體滑移，并且滑移量均小于實際的滑動距離.越靠近中心區(qū)域，滑移量越小.

在磨損的過程中，往往同時也伴隨著疲勞.在疲勞的產(chǎn)生過程中，正應(yīng)力和剪應(yīng)力的變化是非常重要的現(xiàn)象［12］.圖11是有限元預(yù)測模型在載荷為500 N，往復(fù)運動幅值為25 μm的條件下，在不同的摩擦系數(shù)下磨損5 000次后的x方向正應(yīng)力和剪應(yīng)力分布.從圖中可以看出，不同摩擦系數(shù)下的x方向正應(yīng)力和剪應(yīng)力均在磨損邊界處出現(xiàn)峰值，并且隨著摩擦系數(shù)的增大，這個峰值也逐漸增大.因此，在其他條件相同的情況下，較大的摩擦系數(shù)會比較容易導(dǎo)致疲勞的發(fā)生.

圖8 不同的摩擦系數(shù)下磨損5 000次后圓柱的表面輪廓(載荷:500 N，幅值:25 μm)

圖9 不同摩擦系數(shù)的壓力分布(載荷:500 N)

圖10 不同摩擦系數(shù)下移動2.5 μm后的滑移量(載荷:500 N)

圖11 不同摩擦系數(shù)下磨損5 000次后應(yīng)力分布(載荷:500 N，幅值:25 μm)

3 結(jié)論

1)利用有限元的方法，將連續(xù)的磨損過程離散化，選取適當(dāng)?shù)哪p放大因子來縮短計算時間，采用有限元分析軟件來進行磨損過程的仿真，可以為實際工程應(yīng)用提供一種摩擦磨損分析工具.

2)由仿真計算結(jié)果發(fā)現(xiàn)，較大的摩擦系數(shù)會降低接觸面滑移量和接觸壓力，從而降低磨損量.較大的摩擦系數(shù)同時也會導(dǎo)致接觸區(qū)域正應(yīng)力和剪應(yīng)力上升，進而加速接觸疲勞的產(chǎn)生.

［1］溫詩鑄，黃平.摩擦學(xué)原理［M］.第2版.北京:清華大學(xué)出版社，2002:301-302.

［2］蔣翔俊，張優(yōu)云，袁淑霞.燃機端面齒微動磨損分析及數(shù)值仿真［J］.哈爾濱工業(yè)大學(xué)學(xué)報，2011，43(3):75-79.

［3］汪選國.銷盤滑動磨損試驗的仿真方法研究［D］.武漢:武漢理工大學(xué)能源與動力工程學(xué)院，2009:1-7.

［4］馮偉.滑動磨損試驗的有限元法數(shù)字仿真研究［D］.武漢:武漢理工大學(xué)能源與動力工程學(xué)院，2005:25-33.

［5］SFANTOS G K，ALIABADI M H.Wear simulation using an incremental sliding boundary element method［J］.Wear，2006，260(9/10):1119-1128.

［6］PODRA P，ANDERSSON S.Finite element analysis wear simulation of a conical spinning contact considering surface topography［J］.Wear，1999，224(1):13-21.

［7］REZAEI A，VAN PAEPEGEM W，DE BAETS P，et al.Adaptive finite element simulation of wear evolution in radial sliding bearings［J］.Wear，2012，196(1/2):660-671.

［8］CRUZADO A，URCHEGUI MA，GOMEZ X.Finite element modeling and experimental validation of fretting wear scars in thin steel wires［J］.Wear，2012，289:26-38.

［9］宿月文，陳渭，朱愛斌，等.滑動磨損過程有限元分析及銷磨損預(yù)測［J］.中國機械工程，2009，20(13):1573-1581.

［10］MOHD TOBI A L，DING J，BANDAK G，et al.A study on the interaction between fretting wear and cyclic plasticity for Ti-6Al-4V［J］.Wear，2009，267(1/2/3/4):270-282.

［11］MCCOLL IR，DING J，LEEN SB.Finite element simulation and experimental validation of fretting wear［J］.Wear，2004，256(11/12):1114-1127.

［12］DING J，LEEN S B，MCCOLL I R.The effect of slip regime on fretting wear-induced stress evolution［J］.International Journal of Fatigue，2004，26(5):521-531.