李炎，董慶兵，羅振濤，何東

（1.重慶大學(xué) a.機械傳動國家重點實驗室 b.機械工程學(xué)院，重慶 400044；2.重慶齒輪箱有限責(zé)任公司，重慶 402263）

接觸問題廣泛出現(xiàn)在航空航天、汽車等領(lǐng)域，對于粗糙表面接觸應(yīng)力的計算顯得尤為重要。表面粗糙形貌對零件的使用功能、壽命有重大的影響，因此對接觸部位上的應(yīng)力、接觸區(qū)域的大小進行分析具有十分重要的意義。

隨著計算機性能的提高，對于粗糙表面接觸研究，國內(nèi)外學(xué)者們提出了粗糙表面彈性接觸數(shù)值模型。Webster 等[1]提出了光滑無摩擦的實際粗糙表面模型，用于解決實際工程粗糙表面接觸問題，并發(fā)現(xiàn)表面平滑情況有很高的準確性。Venner 等[2]提出了采用多重網(wǎng)格方法來加速各種粗糙表面接觸問題的數(shù)值求解過程。Liu 等[3]首先將離散卷積快速傅里葉變換（Discrete Convolutional Fourier Transform, DC-FFT）運用到求解接觸問題中，提高了工作效率。接著，Liu等[4]進一步利用離散卷積快速傅里葉（DC-FFT）算法加速計算已知粗糙表面壓力和切向力下層狀半空間的應(yīng)力。Wang 等[5]進一步比較了直接求和、多重網(wǎng)格積分和離散卷積-快速傅里葉變換法（DC-FFT），結(jié)果顯示，DC-FFT 相比于其他兩種方法，在計算效率方面有很大的優(yōu)勢。Lin 等[6]提出了一種求解含摩擦各向異性材料二維等溫粗糙表面接觸問題的方法，采用共軛梯度法（Conjugate Gradient Method, CGM）計算接觸壓力，采用疊加法計算各向異性半空間在接觸壓力作用下的二維應(yīng)力。Scholz 等[7]提出了求解各向異性彈性等溫粗糙表面接觸問題的方法，通過快速傅里葉變換，得到各向異性彈性半空間在均勻壓力作用下的表面位移，利用卷積定理計算了各向異性半空間的法向接觸應(yīng)力和牽引接觸應(yīng)力引起的表面位移和應(yīng)力。徐迪初等[8]建立了綜合載荷下粗糙表面彈塑性接觸的確定性模型，考慮微凸峰在接觸時彈塑性變形情況，數(shù)值求解得到了實際接觸面積、壓力分布和微凸峰彈塑性變形。

近些年，有限元方法（Finite Element Method,FEM）已經(jīng)廣泛用于模擬光滑或粗糙表面的接觸問題，是分析粗糙接觸相互作用信息的有力工具。Kosior 等[9]結(jié)合區(qū)域分解法和邊界元法對庫侖摩擦的正常彈性接觸進行了建模。Kogut 等[10]研究了球形和半空間之間的真實接觸，得到了無量綱平均壓力和接觸面積的經(jīng)驗表達式，并將其應(yīng)用于粗糙接觸的研究。Liu 等人[11]建立了具有線性硬化規(guī)律的有限元彈塑性接觸模型來分析二維粗糙表面。Poulio 等[12]研究了具有真實表面形貌的無摩擦彈塑性接觸問題。An等[13]通過譜插值方法研究了粗糙表面之間的相互作用，討論了不同網(wǎng)格分辨率對面積-載荷關(guān)系的影響。有限元法能夠模擬復(fù)雜的材料行為和接觸形狀，然而，有限元模擬復(fù)雜實際粗糙表面的接觸問題會耗費大量的時間。

移動平均濾波是數(shù)字信號處理器中最常見的濾波器，能有效減少隨機噪聲，并保持尖銳的階躍響應(yīng)。Everitt 等[14]將移動平均濾波器引入表面粗糙度處理，對比SOSA 等[15]的實驗結(jié)果，驗證了方法的準確性。

本文以線接觸問題為基礎(chǔ)，建立了二維平面接觸的模型，并對初始粗糙表面進行平均光滑處理，能快速準確求解粗糙表面的接觸問題。同時也考慮了微接觸間的相互作用，能準確地計算局部接觸區(qū)域的應(yīng)力，進行粗糙表面的接觸應(yīng)力計算，并提出0 楊氏模量涂層的構(gòu)想，輸出平面應(yīng)力結(jié)果。文中將以光滑圓面和規(guī)則圓形微凸體、正弦形表面、粗糙表面及平均濾波光滑處理后的粗糙平面為例，將共軛梯度法應(yīng)用于模型中，計算接觸壓力，并通過修正離散卷積傅里葉變換，提高應(yīng)力計算效率，研究綜合載荷下接觸應(yīng)力的分布規(guī)律，研究不同摩擦系數(shù)下應(yīng)力的變化情況。

1 模型

圓柱和平面接觸，可以看成是當(dāng)量圓柱與平面的線接觸，接觸區(qū)域的長度遠大于接觸區(qū)域?qū)挾?，因此可以把接觸區(qū)域視為無限長，是典型的平面應(yīng)變問題。圓柱與實際粗糙表面接觸截面如圖1 所示。采用直角坐標系，x-y平面經(jīng)過粗糙表面尖峰所在最高點位置，圓柱面最低點與x-y平面相切，相切點為坐標系原點，z軸沿公法線指向彈性平面深度方向。通過計算二維接觸的力學(xué)行為，以x-z平面為研究對象。為準確描述表面壓力和亞表層應(yīng)力分布，采用1024×1024 網(wǎng)格離散化。y軸方向?qū)τ嬎憬Y(jié)果不產(chǎn)生任何影響，為節(jié)省計算時間，該方向網(wǎng)格數(shù)量取2。

圖1 圓面與粗糙表面接觸Fig.1 Cylindrical surface in contact with rough surface

假設(shè)圓柱是剛性體，圓面和彈性平面的當(dāng)量楊氏模量如式(1)所示。

式中：E為當(dāng)量楊氏模量；E1、E2分別為剛性圓柱和平面的楊氏模量；μ1、μ2分別為剛性圓柱和平面的泊松比。

根據(jù)赫茲彈性接觸理論，接觸半徑aH可用式(2)表示。

式中：R為曲率半徑；W為載荷。

在接觸區(qū)上的最大接觸壓力pH的計算如式(3)所示。

接觸面上其余各點的壓力呈半橢圓規(guī)律分布，壓力計算見式(4)。

圓柱與平面之間的間隙h由式(5)計算。

式中：h0為初始間隙；ν(x,y)為表面彈性變形；s(x,y)為表面粗糙度。

壓力分布的本構(gòu)方程如式(6)所示。在接觸間隙為0 的區(qū)域，壓力大于0；間隙大于0 的區(qū)域，壓力等于0[16]。采用改進的共軛梯度法進行迭代求解[17-18]。

應(yīng)力場用式(7)表示[4]。

計算接觸應(yīng)力的DC-FFT 方法中，將計算區(qū)域在x和y方向擴展至原區(qū)域的2 倍，壓力在擴展區(qū)域用0 來補充。針對本研究的圓柱與平面接觸，采用修正的DC-FFT 方法，在x方向使用DC-FFT 中的方法擴充補0，而在y方向復(fù)制擴充擴展區(qū)域[19-21]。針對文中二維平面應(yīng)變問題中的應(yīng)力分量遵循式(8)。

通過離散卷積傅里葉變換可求得各分量的值。通過計算圓面與平面接觸時的Mises 應(yīng)力，來表示應(yīng)力分布情況，如式(9)所示。

式中：σvm為Mises 應(yīng)力；τxy= 0；τyz= 0。

本文假設(shè)圓柱和平面區(qū)域之間的未接觸間隙空白區(qū)域存在一種楊氏模量為0 的特殊涂層材料，平面的其他部分和材料之間不存在相互作用。將涂層材料區(qū)域的Mises 應(yīng)力賦值為0，得到求解區(qū)域內(nèi)實際平面形狀的應(yīng)力分布情況。本文將引入進行信號處理的數(shù)字處理器中的數(shù)字濾波方法[22]，采用最常用的移動平均濾波處理實際粗糙表面的隨機高頻波峰波谷，并模擬表面的磨損過程[14]。

2 結(jié)果及分析

在本文中，平面的楊氏模量E2=100 GPa，泊松比μ2=0.3。施加在剛性圓柱上的法向載荷W=172.615 N，赫茲半徑aH=0.2 mm。

2.1 圓面與光滑平面接觸

圓面與光滑平面接觸，不考慮摩擦力，是典型的赫茲接觸。x方向求解區(qū)域為－1.2aH≤x≤1.2aH，z軸方向為2.0aH。圓面與光滑平面接觸的Mises 應(yīng)力分布和赫茲壓力如圖2 所示。在圖2a 中，應(yīng)力分布與文獻[23]中一致。通過對比CGM 迭代計算和公式(4)計算的壓力曲線（見圖2b），表明模型計算的準確性，可用來計算圓面與規(guī)則微凸體粗糙表面和實際粗糙表面接觸時的應(yīng)力分布情況。

圖2 圓面與光滑平面接觸的應(yīng)力分布及模型和解析解的壓力對比Fig.2 a) stress distribution of cylindrical surface in contact with smooth plane and b) comparison of pressure between solutions by the present model and analytical solutions

2.2 圓面與有單個圓形凸體平面接觸

采用光滑圓面與圓形微凸體的接觸，模擬圓面與粗糙表面規(guī)則尖峰輪廓接觸時的應(yīng)力分布情況。圖3為圓面與不同半徑圓形凸體平面的Mises 應(yīng)力分布，x方向的求解區(qū)域為–0.5aH≤x≤0.5aH，z軸方向為0.5aH。

圖3 不同半徑圓形凸體Mises 應(yīng)力云圖Fig.3 The Mises stress contour for a circular asperity with different radii: a) 0.1 times Hertz radius; b) 0.2 times Hertz radius

在圖3 中，圓面與有圓形凸體的平面接觸時，主要與平面上的圓形凸體部分接觸。由于圓形凸體接觸面積小，彈性變形遠遠大于在光滑表面上接觸產(chǎn)生的彈性變形，Mises 應(yīng)力更大，且應(yīng)力主要集中在接觸區(qū)內(nèi)，導(dǎo)致平面表層出現(xiàn)應(yīng)力集中現(xiàn)象。圓形凸體半徑增大，Mises 應(yīng)力減小。

2.3 圓面與有圓形凸體平面接觸模型驗證

為了驗證模型的準確性，用有限元軟件建立剛性圓面與有圓形凸起平面的二維模型，平面區(qū)域網(wǎng)格劃分如圖4 所示。為看清模型網(wǎng)格劃分，圖4 中對實際計算的平面網(wǎng)格劃分進行了稀疏處理。對比有限元仿真軟件和數(shù)值模型沿z軸方向Mises 應(yīng)力計算結(jié)果，如圖5 所示。有限元仿真軟件計算時，網(wǎng)格細小，網(wǎng)格疏密程度由接觸區(qū)域到非接觸區(qū)域逐漸稀疏，并且覆蓋整個平面區(qū)域，能考慮到整個平面區(qū)域的彈性變形及單元間的相互作用，還考慮了微凸體中的網(wǎng)格過度扭曲變形和尖角處應(yīng)力集中問題。因此，數(shù)值模型和有限元模型結(jié)果對比有偏差，但計算結(jié)果和仿真結(jié)果基本符合。

圖4 有限元模型局部區(qū)域網(wǎng)格細化Fig.4 Finite element model with locally refined meshes

圖5 有限元仿真和模型計算結(jié)果沿z 軸方向應(yīng)力分布Fig.5 Comparison of stresses along the z axis between the present model and finite element method

2.4 圓面與有多個圓形凸體平面接觸

假設(shè)規(guī)則的圓形凸體是平面上的一個表面粗糙度輪廓，以連續(xù)多個凸體形狀近似模擬實際表面，可以更好地反映出圓面與粗糙表面規(guī)則輪廓接觸時的應(yīng)力分布情況。沿z軸方向的應(yīng)力變化如圖6 所示。在相同的圓形凸體半徑下，圓面分別與有3、5 個圓形凸體的平面接觸。對比圓面和單個圓形凸體接觸，由于接觸區(qū)域增大，應(yīng)力大幅減小，曲線偏向x坐標軸正方向。在圓形凸體數(shù)相同時，圓形凸體半徑由0.1 倍赫茲半徑變?yōu)?.2 倍赫茲半徑，Mises 應(yīng)力會出現(xiàn)較大程度的減小。

圖6 圓面與粗糙表面規(guī)則輪廓接觸時沿z 軸應(yīng)力分布Fig.6 Stress distribution along the z axis

2.5 圓面與考慮切向力的圓形凸體平面接觸

在實際接觸情況中，粗糙表面都存在摩擦力。圓面與有圓形微凸體的平面接觸時，考慮切向摩擦力的情況，更能反映真實情況。圓面與不同半徑圓形凸體平面的Mises 應(yīng)力分布情況如圖7 所示。圓形凸體半徑為0.1 倍赫茲半徑，x方向求解區(qū)域為–0.5aH≤x≤0.5aH。

圖7 不同摩擦系數(shù)時0.1 倍赫茲半徑圓形凸體的應(yīng)力分布Fig.7 Stress distribution for a circular asperity with a radius of Hertz under different friction coefficients: a) coefficient of friction is 0.15; b) coefficient of friction is 0.3

摩擦系數(shù)為0.15 時（見圖7a），切向摩擦?xí)?dǎo)致單個圓形凸體的Mises 應(yīng)力略有增加，應(yīng)力云圖形狀基本不變，但在應(yīng)力集中區(qū)域，云圖會產(chǎn)生偏移。摩擦系數(shù)為0.3 時（見圖7b），對Mises 應(yīng)力產(chǎn)生的影響更大，平面的Mises 應(yīng)力增大，圓形凸體頂部應(yīng)力集中區(qū)域的云圖偏移增大，且形狀不再完整，有部分應(yīng)力轉(zhuǎn)移到圓形凸體與圓面接觸的表面區(qū)域。

在考慮切向力的情況下，沿z軸方向的應(yīng)力變化如圖8 所示。摩擦系數(shù)為0.15 時，與單個圓形凸體接觸時不同，在圓面與圓形凸體頂點處接觸時有最大的應(yīng)力點，應(yīng)力沿z軸迅速降低到較小值，再增長到峰值，曲線在峰值頂點處與摩擦系數(shù)為0.3 時的曲線重合。當(dāng)摩擦系數(shù)為0.3 時，這個現(xiàn)象非常明顯，對應(yīng)于云圖（圖7）中應(yīng)力形狀發(fā)生變化，在初始接觸點的應(yīng)力最大，部分較大應(yīng)力出現(xiàn)在圓形凸體的下表層近處。切向力會影響圓形微凸體內(nèi)Mises 應(yīng)力的分布，摩擦系數(shù)的變化主要影響應(yīng)力云圖的偏移，以及初始接觸頂點和主要應(yīng)力集中區(qū)域峰值點之間的應(yīng)力大小。

圖8 考慮切向摩擦力時沿z 軸的應(yīng)力分布Fig.8 Stress distribution along the z-axis

2.6 圓面與正弦表面接觸

圓面與具有正弦形粗糙表面的平面接觸，假設(shè)連續(xù)規(guī)則的正弦形狀是平面上的表面粗糙度輪廓，可以反映出圓面與連續(xù)周期性規(guī)則粗糙表面輪廓接觸時的應(yīng)力分布情況。選取正弦函數(shù)0.01cos (2πx/0.2)，圓面與正弦形狀的粗糙平面接觸時的應(yīng)力分布情況如圖9 所示。當(dāng)考慮切向摩擦力時，在圓面與正弦形接觸部分，出現(xiàn)與單個圓形凸體類似的現(xiàn)象，如圖10 所示。

圖9 未考慮切向力時正弦形表面應(yīng)力分布Fig.9 Tangential distribution is not considered a sinusoidal surface stress distribution

圖10 不同摩擦系數(shù)時正弦形表面應(yīng)力分布Fig.10 Sinusoidal surface stress distribution under different friction coefficients: a) coefficient of friction is 0.15; b) coefficient of friction is 0.30

2.7 圓面與實際表面接觸

實驗測得的齒面三維形貌如圖11a 所示，去除宏觀形狀誤差得到的實測表面三維形貌如圖11b 所示。將原有齒輪漸開線輪廓去除后，得到表面的真實粗糙度峰值用于計算。用一維移動平均低通濾波平滑表面輪廓，以去除測量中的高頻噪聲。首先在水平方向上應(yīng)用濾波，然后在豎直方向上根據(jù)當(dāng)前節(jié)點的值再次濾波。為了解釋快速磨合過程[15,24]，應(yīng)用另一個移動平均低通濾波。這是一個2D 過濾方法，使用水平和豎直方向上當(dāng)前節(jié)點和鄰近節(jié)點的值。該濾波方法僅適用于z=0 以上的部分，以限制濾波的最高極限[14]。

圖11 齒輪表面三維形貌Fig.11 3D topography of gear surface: a) initial gear surface; b) gear surface after removing macro shapes

平均濾波處理前后的齒輪表面二維形貌如圖12所示。表面的大量微小尖峰已被過濾，表面微峰降低了峰值，表面得到光滑處理。首先，在齒面接觸線方向進行一維平均，再在垂直接觸線方向進行一維平均，均基于當(dāng)前節(jié)點兩邊的值（本文選取9 個點），去除測量中出現(xiàn)的尖峰。最后，對大于0 的峰值點及環(huán)繞它的點（本文選取9 個點）進行二維平均，去除快速磨合過程中尖峰消失的部分。

圖12 齒輪表面二維形貌Fig.12 2D surface topography: a) two-dimension tooth topography before average; b) two-dimension tooth topography after average

在表面粗糙度測量中，支承長度率是唯一用于評定表面形狀特性的參數(shù)。將對應(yīng)于各個水平截距的tp值（支承長度率）繪制成曲線，得到支承長度率曲線，即Abbot 曲線[25-26]。Abbot 曲線可以很好地說明原始粗糙表面和磨合表面之間經(jīng)過平均濾波處理后的變化情況[14]。經(jīng)過平均濾波處理后，表面最高峰峰值30%的Abbot 曲線如圖13 所示。Deepak 等[27]研究了磨合和齒輪滾動接觸疲勞測試時表面形貌的變化，結(jié)果表明，磨合是一個快速的過程，其中變化最大的是凹凸不平的塑性變化使表面形貌發(fā)生的變化。因此，對于穩(wěn)定磨合期的表面形貌，在進行數(shù)值計算時，需采用合適的方法模擬表面的磨合過程。Everitt 等人[14]在論文中，采用移動平均濾波方法對噴丸和磨削粗糙表面進行處理，取得了較好的濾波效果，峰值降低到原始高度的70%左右，與SOSA 等[15]的齒輪表面磨合實驗結(jié)果一致。他們還將粗糙度數(shù)值用于研究表面熱彈性流體動力潤滑和疲勞模擬，并對齒輪點蝕出現(xiàn)位置進行了解釋。本文中的表面粗糙度來源于實際磨削齒輪表面，從AGMA 典型齒輪磨合測試中的實驗數(shù)據(jù)可知[28]，對于磨削加工而成的齒面，磨合后的齒輪表面粗糙度降低20%～40%。從圖13 中可以看出，平均濾波使峰值降低到原始高度的70%左右，特別是對于圖中坐標原點最高微凸峰的降低最為明顯。經(jīng)過處理后的表面形貌可以代表齒輪表面磨合后的形貌。

圖13 Abbot 曲線Fig.13 Abbot curve

本文選取圖12b 中一條水平數(shù)值線上經(jīng)過平均處理后的粗糙度數(shù)據(jù)進行接觸計算，處理前后該粗糙度的曲線如圖14 所示。經(jīng)過處理后的粗糙度形貌為穩(wěn)定磨合階段的形貌，以該階段的粗糙度數(shù)值進行計算能準確地反映實際情況中的接觸情況。圓面和有實際表面粗糙度的平面接觸，和具有規(guī)則微凸體形狀的平面相比，更能反映實際接觸問題中的應(yīng)力分布情況。實驗測量區(qū)域的寬度為1.1 mm，文中x方向求解區(qū)域為–2.75aH≤x≤2.75aH，z軸方向求解區(qū)域為0.01aH。

圖14 一條水平線的粗糙表面輪廓Fig.14 Roughness profile of a horizontal line

在圓面與實際粗糙表面接觸的Mises 應(yīng)力云圖（圖15）中，圓面主要和平面的不同尖峰接觸（表面輪廓高度在–3～2 μm 內(nèi)，圖中是對平面部分區(qū)域進行輸出，–1.0aH≤x≤1.0aH，z軸方向0.01aH）。云圖中較大Mises 應(yīng)力出現(xiàn)在圓面與粗糙接觸表面輪廓的尖峰處。圓面和原始粗糙表面，由于存在很小的尖峰，導(dǎo)致在局部會形成巨大的應(yīng)力集中。從圖15b 中可以看出，經(jīng)過光滑處理的粗糙表面，由于表面微凸峰尖銳程度減小，Mises 應(yīng)力大幅減少，應(yīng)力集中程度減輕。圓面與平均處理后的粗糙表面接觸時，考慮切向摩擦力下的應(yīng)力云圖見圖16。當(dāng)摩擦系數(shù)為0.15 時，對比沒有考慮摩擦力的情況（圖15b），應(yīng)力云圖略有偏移，形狀有變化，應(yīng)力略有增加；當(dāng)摩擦系數(shù)為0.3 時，應(yīng)力云圖形狀變化較大，偏移程度加大，應(yīng)力集中區(qū)域擴大，應(yīng)力增加。摩擦力會對粗糙表面的應(yīng)力分布產(chǎn)生影響，摩擦系數(shù)越大，應(yīng)力云圖形狀偏移增大，應(yīng)力也增大。

圖15 圓面與粗糙表面接觸應(yīng)力分布Fig.15 Stress distribution of circular surface and rough surface: a) initial actual surface; b) average rough surface

圖16 圓面與不同摩擦系數(shù)實際表面接觸時的應(yīng)力分布Fig.16 Stress distribution of cylindrical surface in contact with the actual surface at different friction coefficient: a) friction coefficient of 0.15; b) friction coefficient of 0.3

3 結(jié)論

1）圓面與粗糙平面接觸時，在微凸體區(qū)域產(chǎn)生應(yīng)力集中，摩擦力會導(dǎo)致應(yīng)力發(fā)生偏移。本研究搭建的應(yīng)力快速計算模型可為剪切應(yīng)力導(dǎo)致的多軸疲勞失效提供量化的應(yīng)力結(jié)果。

2）磨合作用能夠有效降低粗糙界面尖端的應(yīng)力集中。采用移動平均濾波對實際粗糙表面進行光滑處理，在保留表面特征的基礎(chǔ)上，能準確地描述表面經(jīng)過磨合后的粗糙情況。

3）本文構(gòu)建了考慮表面微觀形貌的應(yīng)力分布情況，由于形貌特征產(chǎn)生的鏤空區(qū)域近似為不影響應(yīng)力分布的涂層結(jié)構(gòu)，從而準確描述了具體形貌表面下的應(yīng)力分布。

4）本模型能夠快速精確計算出接觸區(qū)域的近場應(yīng)力以及應(yīng)力集中分布，可以為應(yīng)力集中產(chǎn)生的裂紋萌生、擴展等破壞提供理論基礎(chǔ)，同時為齒輪等傳動部件的表面完整性制造和優(yōu)化提供依據(jù)。