粗糙面滑動摩擦熱-應力-磨損耦合仿真分析

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

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

0 引言

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

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

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

1 計算模型建立

1.1 粗糙面模型建立

W-M分形函數模型[9]，其三維函數表達式：

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

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

圖1 銷-盤模型

1.2 計算模型假設條件

此計算模型的假設條件如下：

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

（2）摩擦過程符合庫倫定律，摩擦盤在轉動過程中，由于時間較短，默認摩擦系數為定值。

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

（4）忽略熱輻射。

1.3 計算模型熱邊界條件

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

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

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

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

默認接觸面產生熱量在摩擦塊和摩擦盤間自由分配，根據能量守恒和熱傳導的傅立葉定律，摩擦副熱傳導方程表示為：

式中，λ為熱導率，ρ為密度，c為比熱容。

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

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

1.4 計算模型Archard磨損理論

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

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

假設摩擦塊有限元模型接觸節點為I，積分步驟為j，則磨損增量表示為：

累計磨損表示為：

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

2.1 有限元模型建立

摩擦塊材料為機械化合金，摩擦盤材料為鑄鋼，摩擦系數設為0.3，材料性能參數見表1。摩擦塊和摩擦盤六面體單元屬性均為C3D8T，摩擦塊單元數為54000，摩擦盤單元數為19136。為提高仿真結果準確性，摩擦塊接觸區域網格做加密處理，劃分后模型網格如圖2所示。

表1 材料屬性

圖2 有限元模型網格劃分

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

有限元仿真簡化流程如圖3所示，在ABAQUS軟件中完成材料賦予、接觸屬性的創建等。仿真過程中，設置摩擦塊網格為ALE自適應網格區域，通過提取ABAQUS接觸壓力和滑動位移，并根據式（5）編寫的磨損子程序代碼，計算節點的磨損增量，并按照計算的偏移方向，不斷調整表面網格，仿真模擬磨損深度處于不斷變化的磨損過程，計算摩擦塊的累計磨損量。基于不斷更新變化的表面形貌，在溫度-位移耦合分析步中，得到不斷更新的表面壓力、溫度場的分布狀態。

圖3 有限元仿真簡化流程

2.2 仿真分析結果與討論

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

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

滑動初期，粗糙實體摩擦塊表面微凸體與摩擦盤接觸后相互摩擦，此時實際摩擦面積為微凸體接觸面積，隨著摩擦盤的轉動，摩擦時間的增長，在法向壓力持續作用下，摩擦表面微凸體逐漸磨損變形，摩擦表面微凸體不斷被磨平，摩擦面積不斷增大，從而使得接觸面整體的接觸壓力逐漸下降，并逐漸在微凸體的邊緣呈現出較高壓力值。滑動后最高接觸壓力下降到419.7 MPa，滑動后最高接觸壓力下降為288.9 MPa。

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

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

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

圖6 節點磨損高度隨時間變化曲線

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

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

圖8 節點溫度隨時間變化曲線

3 結論

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

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

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

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