安先龍，王佳，趙誠，高超

（中機(jī)科（北京）車輛檢測工程研究院有限公司，北京 102100）

0 引言

工程表面微觀上是凹凸不平的，其摩擦過程實質(zhì)為摩擦副粗糙表面上微凸體間的相互摩擦過程[1]。摩擦過程中，摩擦功產(chǎn)生大量的熱量，改變摩擦副摩擦磨損特性，同時由于壓力的作用，摩擦副間磨損的產(chǎn)生也會改變接觸表面的壓力、溫度場的分布狀態(tài)。因此，摩擦過程涉及復(fù)雜的應(yīng)力、溫度和磨損耦合現(xiàn)象[2，3]。

摩擦副粗糙面具有分形特征，Mandelbort[4]創(chuàng)立分形幾何學(xué)，并在Weierstrass函數(shù)的基礎(chǔ)上，提出了分形W-M函數(shù)，因其具有自相似性、自仿射性，因此，可被用來描述粗糙表面形態(tài)。后來有不少多學(xué)者對分形理論進(jìn)行了研究，黃健萌等[5]通過建立粗糙面滑動接觸模型，仿真分析粗糙實體的熱、應(yīng)力分布規(guī)律。丁雪興[6]通過分形理論建立了磨粒磨損的預(yù)測模型，并推導(dǎo)出磨損率和分形模型參數(shù)之間的關(guān)系。Wang等[7]通過分形理論表征分形粗糙表面，分析了粗糙表面微凸體處的摩擦溫升。Gong和Komvopoulos[8]通過建立球體彈塑性有限元模型，研究了球體幾何形狀對接觸壓力分布和表面應(yīng)力-應(yīng)變場的影響規(guī)律。

依據(jù)W-M分形函數(shù)，在Matlab、Catia軟件中聯(lián)合建立接觸式干滑動模型，然后基于應(yīng)力、溫度和磨損耦合方法，通過編寫Umeshmotion子程序代碼，在ABAQUS軟件中建立微凸體表面干滑動磨損仿真模型，探究滑動過程中粗糙面微凸體的接觸壓力、磨損和溫度場的分布影響規(guī)律。

1 計算模型建立

1.1 粗糙面模型建立

W-M分形函數(shù)模型[9]，其三維函數(shù)表達(dá)式：

式中，L為取樣長度；G為特征尺度系數(shù)；cosφm，n表示（0，2）間均勻分布的隨機(jī)相位；D為分形維數(shù)；M為重構(gòu)表面時疊加峰值的數(shù)量，一般選為10；γ為隨機(jī)輪廓的空間頻率；n為頻率系數(shù)，nmax為頻率的上限。分形維數(shù)和特征尺度系數(shù)是決定粗糙面形貌的主要參數(shù)[10]，為減少有限元仿真分析占用過多計算資源，僅對摩擦副上粗糙面采用分形理論建模，下粗糙面則采用理想平面實體。

分形接觸模型如圖1所示。摩擦塊尺寸：內(nèi)半徑（半徑）270 μm，外徑（半徑）420 μm，厚度80 μm，摩擦盤尺寸：內(nèi)徑（半徑）270 μm，外徑（半徑）430 μm，厚度為100 μm。

圖1 銷-盤模型

1.2 計算模型假設(shè)條件

此計算模型的假設(shè)條件如下：

（1）摩擦副材料（摩擦塊和摩擦盤）各向同性。

（2）摩擦過程符合庫倫定律，摩擦盤在轉(zhuǎn)動過程中，由于時間較短，默認(rèn)摩擦系數(shù)為定值。

（3）忽略摩擦過程中慣性力和振動造成的影響。

（4）忽略熱輻射。

1.3 計算模型熱邊界條件

熱邊界條件是溫度場周圍的表面熱交換條件，共有3種[11]。如下所示：

（1）第一類邊界條件表明邊界上的溫度值是已知的，故不存在第一類邊界條件。

（2）第二類邊界條件規(guī)定了邊界上熱流密度，其方程表示為：

式中，q為邊界熱流密度，μ摩擦系數(shù)，p為接觸壓力，?γ為相對滑動速度。

默認(rèn)接觸面產(chǎn)生熱量在摩擦塊和摩擦盤間自由分配，根據(jù)能量守恒和熱傳導(dǎo)的傅立葉定律，摩擦副熱傳導(dǎo)方程表示為：

式中，λ為熱導(dǎo)率，ρ為密度，c為比熱容。

（3）第三類邊界條件定義了流體溫度以及物體與流體之間的傳熱系數(shù)，流體與固體表面間的對流換熱[12]方程表示為：

式中，q為散熱的熱流密度；h為對流換熱系數(shù)；Ts為固體溫度；Tf為流體的溫度。

1.4 計算模型Archard磨損理論

使用Archard理論模擬表面磨損損傷[13]，其廣義增量磨損高度定義如下：

式中，△h為磨損高度增量，k為無量綱常數(shù)，是表征材料耐磨性的磨損系數(shù)，可由試驗獲得；p為接觸壓力；△s為相對滑動位移增量。

假設(shè)摩擦塊有限元模型接觸節(jié)點(diǎn)為I，積分步驟為j，則磨損增量表示為：

累計磨損表示為：

2 熱力耦合和磨損仿真分析

2.1 有限元模型建立

摩擦塊材料為機(jī)械化合金，摩擦盤材料為鑄鋼，摩擦系數(shù)設(shè)為0.3，材料性能參數(shù)見表1。摩擦塊和摩擦盤六面體單元屬性均為C3D8T，摩擦塊單元數(shù)為54000，摩擦盤單元數(shù)為19136。為提高仿真結(jié)果準(zhǔn)確性，摩擦塊接觸區(qū)域網(wǎng)格做加密處理，劃分后模型網(wǎng)格如圖2所示。

表1 材料屬性

圖2 有限元模型網(wǎng)格劃分

為明顯區(qū)分磨損所造成的影響，因此選取工況為：第一步為壓力分析步，即在40 μs內(nèi)將20 MPa法向壓力通過摩擦盤施加到摩擦塊上；第二步為滑動分析步，即在法向壓力的作用下摩擦盤以10 m/s的線速度單向勻速旋轉(zhuǎn)40 μs。在摩擦盤中心創(chuàng)建參考點(diǎn)RP-1，由RP-1來控制摩擦盤的轉(zhuǎn)動和載荷的施加。仿真過程中，摩擦盤僅保留Z方向的自由度，其他方向位移施加約束限制。

有限元仿真簡化流程如圖3所示，在ABAQUS軟件中完成材料賦予、接觸屬性的創(chuàng)建等。仿真過程中，設(shè)置摩擦塊網(wǎng)格為ALE自適應(yīng)網(wǎng)格區(qū)域，通過提取ABAQUS接觸壓力和滑動位移，并根據(jù)式（5）編寫的磨損子程序代碼，計算節(jié)點(diǎn)的磨損增量，并按照計算的偏移方向，不斷調(diào)整表面網(wǎng)格，仿真模擬磨損深度處于不斷變化的磨損過程，計算摩擦塊的累計磨損量?；诓粩喔伦兓谋砻嫘蚊?，在溫度-位移耦合分析步中，得到不斷更新的表面壓力、溫度場的分布狀態(tài)。

圖3 有限元仿真簡化流程

2.2 仿真分析結(jié)果與討論

圖4為粗糙實體摩擦塊接觸面不同時刻的接觸壓力分布，從圖中可以看出，滑動分析步初期，摩擦塊端面接觸壓力分布不均勻，微凸體尖端處產(chǎn)生了巨大的應(yīng)力集中，最大接觸壓力均達(dá)到1853 MPa，是法向壓力的92.65倍。

圖4 滑動過程接觸壓力分布

滑動初期，粗糙實體摩擦塊表面微凸體與摩擦盤接觸后相互摩擦，此時實際摩擦面積為微凸體接觸面積，隨著摩擦盤的轉(zhuǎn)動，摩擦?xí)r間的增長，在法向壓力持續(xù)作用下，摩擦表面微凸體逐漸磨損變形，摩擦表面微凸體不斷被磨平，摩擦面積不斷增大，從而使得接觸面整體的接觸壓力逐漸下降，并逐漸在微凸體的邊緣呈現(xiàn)出較高壓力值?；瑒雍笞罡呓佑|壓力下降到419.7 MPa，滑動后最高接觸壓力下降為288.9 MPa。

圖5為滑動摩擦塊的節(jié)點(diǎn)磨損分布云圖。從圖中可以看出，接觸壓力分布云圖與磨損分布云圖具有較高的一致性，此現(xiàn)象符合Archard磨損定律，即磨損與接觸壓力成正相關(guān)?；瑒雍?，摩擦塊形成了大量塊狀磨損區(qū)，微凸體中心處具有較大的磨損高度。最高微凸體節(jié)點(diǎn)，在滑動后磨損高度為6.51。選取磨損最為嚴(yán)重的A、B、C、D節(jié)分析點(diǎn)節(jié)點(diǎn)磨損、溫度隨時間的變化規(guī)律，所在位置如圖5所示。

圖5 滑動后摩擦塊磨損深度分布

圖6為A、B、C、D四個節(jié)點(diǎn)磨損高度隨滑動時間變化的曲線。四個節(jié)點(diǎn)的磨損高度變化規(guī)律相似性較高，可分為兩個階段。第一個階段為壓力分析步階段（0～40μs），微凸體節(jié)點(diǎn)經(jīng)過快速偏移到緩慢偏移的過程。這是因為接觸初期，接觸面積為微凸體接觸面積，壓力的增加致使微凸體節(jié)點(diǎn)快速沿z方向偏移，但隨著接觸面積的增大，節(jié)點(diǎn)應(yīng)力降低，節(jié)點(diǎn)偏移速率減緩。除此之外，接觸初期，A節(jié)點(diǎn)為最高微凸體，接觸后迅速向z方向偏移，但B、C、D節(jié)點(diǎn)在未和摩擦盤接觸時，仍發(fā)生輕微節(jié)點(diǎn)偏移，考慮為A節(jié)點(diǎn)受到壓力的作用，致使摩擦塊材料發(fā)生彈塑性變形，從而使得B、C、D節(jié)點(diǎn)沿壓力方向發(fā)生輕微偏移所致。第二個階段滑動分析步（40 μs～80 μs），此時表面微凸體節(jié)點(diǎn)的偏移速率由公式（5）archard磨損定律計算所得，滑動初期，磨損迅速增加，隨著接觸面積的進(jìn)一步增大，接觸壓力進(jìn)一步下降，磨損速率變得緩慢。

圖6 節(jié)點(diǎn)磨損高度隨時間變化曲線

滑動后模型的溫度場分布如圖7所示。溫度場分布以塊狀熱區(qū)為主，滑動后最高溫度為76.34℃，A、B、C、D 4個節(jié)點(diǎn)的溫度隨時間變化曲線如圖8所示，滑動初期，微凸體節(jié)點(diǎn)摩擦產(chǎn)生的熱量來不及傳遞，溫度聚集在微凸體內(nèi)部[14]，節(jié)點(diǎn)溫度迅速增加。以溫升較大的節(jié)點(diǎn)B為例，除接觸壓力較大影響外，其位置靠近摩擦塊外徑，相對滑動速度較大也是造成溫升較大的原因。節(jié)點(diǎn)B在3.18內(nèi)迅速溫升至171.847℃，隨著滑動的進(jìn)行，節(jié)點(diǎn)B微凸體被磨平，其他微凸體陸續(xù)進(jìn)入摩擦過程，接觸面積不斷增大，接觸壓力不斷下降，磨損消耗了部分熱量導(dǎo)致節(jié)點(diǎn)產(chǎn)生的熱量減少，同時加上節(jié)點(diǎn)熱傳導(dǎo)的散熱作用，溫度大幅下降后趨于穩(wěn)定。

圖7 滑動后摩擦塊溫度場分布

圖8 節(jié)點(diǎn)溫度隨時間變化曲線

3 結(jié)論

通過在ABAQUS軟件構(gòu)建粗糙面微凸體的干滑動磨損仿真模型，基于熱-應(yīng)力-磨損耦合分析方法，研究分析了滑動過程中磨損行為對接觸端面應(yīng)力、溫度場分布的影響規(guī)律。得到以下結(jié)論：

（1）摩擦過程中的磨損行為對摩擦面接觸壓力及溫度場分布影響較大。

（2）隨著滑動的進(jìn)行，摩擦塊表面微凸體在接觸后產(chǎn)生磨損，接觸表面逐漸被磨平，接觸面積不斷增大，接觸壓力逐漸降低，磨損速率變得緩慢；

（3）隨著滑動的進(jìn)行，磨損表面逐漸平滑，溫度場分布逐漸變均勻，主要以塊狀熱區(qū)分布為主，微凸體溫度出現(xiàn)閃溫后迅速下降，之后趨于平穩(wěn)。