考慮表面形貌的接觸應力分析

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

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

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

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

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

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

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

1 模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中：σvm為Mises 應力；τxy= 0；τyz= 0。

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

2 結果及分析

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

2.1 圓面與光滑平面接觸

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

圖2 圓面與光滑平面接觸的應力分布及模型和解析解的壓力對比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 圓面與有單個圓形凸體平面接觸

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

圖3 不同半徑圓形凸體Mises 應力云圖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 中，圓面與有圓形凸體的平面接觸時，主要與平面上的圓形凸體部分接觸。由于圓形凸體接觸面積小，彈性變形遠遠大于在光滑表面上接觸產生的彈性變形，Mises 應力更大，且應力主要集中在接觸區內，導致平面表層出現應力集中現象。圓形凸體半徑增大，Mises 應力減小。

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

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

圖4 有限元模型局部區域網格細化Fig.4 Finite element model with locally refined meshes

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

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

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

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

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

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

圖7 不同摩擦系數時0.1 倍赫茲半徑圓形凸體的應力分布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

摩擦系數為0.15 時（見圖7a），切向摩擦會導致單個圓形凸體的Mises 應力略有增加，應力云圖形狀基本不變，但在應力集中區域，云圖會產生偏移。摩擦系數為0.3 時（見圖7b），對Mises 應力產生的影響更大，平面的Mises 應力增大，圓形凸體頂部應力集中區域的云圖偏移增大，且形狀不再完整，有部分應力轉移到圓形凸體與圓面接觸的表面區域。

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

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

2.6 圓面與正弦表面接觸

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

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

圖10 不同摩擦系數時正弦形表面應力分布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 所示。將原有齒輪漸開線輪廓去除后，得到表面的真實粗糙度峰值用于計算。用一維移動平均低通濾波平滑表面輪廓，以去除測量中的高頻噪聲。首先在水平方向上應用濾波，然后在豎直方向上根據當前節點的值再次濾波。為了解釋快速磨合過程[15,24]，應用另一個移動平均低通濾波。這是一個2D 過濾方法，使用水平和豎直方向上當前節點和鄰近節點的值。該濾波方法僅適用于z=0 以上的部分，以限制濾波的最高極限[14]。

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

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

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

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

圖13 Abbot 曲線Fig.13 Abbot curve

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

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

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

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

圖16 圓面與不同摩擦系數實際表面接觸時的應力分布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 結論

1）圓面與粗糙平面接觸時，在微凸體區域產生應力集中，摩擦力會導致應力發生偏移。本研究搭建的應力快速計算模型可為剪切應力導致的多軸疲勞失效提供量化的應力結果。

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

3）本文構建了考慮表面微觀形貌的應力分布情況，由于形貌特征產生的鏤空區域近似為不影響應力分布的涂層結構，從而準確描述了具體形貌表面下的應力分布。

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