含硬涂層結構的表面微觀接觸新模型

李 玲, 李治強， 張 健, 蔡安江

(西安建筑科技大學 機電工程學院, 西安 710055)

制造業的不斷發展，對機械加工精度和加工效率提出了更高的要求，機床的動態性能在實際生產中就顯得非常重要。機床結合面的動態參數對機床整體動態特性起著決定性影響，它們決定機床的加工精度和穩定性。因此，開展結合面動態特性研究對揭示機床整體動態特性具有重要的理論意義和實際應用價值。

結合面的研究主要是粗糙表面間接觸行為的研究，而關于粗糙表面間接觸行為的研究以微觀表面微凸體的研究為主[1]。國內外大量學者對粗糙表面微凸體采用不同的研究方法，提出了不同的粗糙表面接觸模型。其中應用最為廣泛的表面接觸模型是Greenwood等[2]提出的統計接觸模型(Greenwood-Williamson, GW模型)。GW模型是一個根據Hertz理論提出的粗糙表面統計接觸模型，模型假設：① 微凸體的峰頂是球體；② 具有相同的曲率半徑；③ 微凸體高度服從高斯分布；④ 表面粗糙度各向同性；⑤ 微凸體只進行彈性接觸，且各個微凸體的變形相互獨立。大多學者基于GW模型，對微觀接觸模型進行了拓展，Kogut等[3]提出描述粗糙表面的彈性、彈塑性和全塑性形變的接觸模型；Chang等[4]根據微凸體塑性變形體積守恒原理，提出了塑性接觸模型；趙永武等[5-7]基于接觸力學理論提出一種新的粗糙表面彈塑性微觀接觸模型(簡稱ZMC模型)；溫淑花等[8-9]基于球體與平面的接觸理論和粗糙表面的接觸分形理論，提出接觸剛度和接觸阻尼模型；尤晉閩等[10-11]利用分形接觸理論，建立了結合面法向接觸分形模型；李玲等[12]根據在變形狀態轉變的臨界點處微凸體真實接觸面積與接觸載荷轉化均滿足連續條件，建立了結合面真實接觸面積、接觸載荷與接觸剛度模型；肖會芳等[13]采用GW模型描述粗糙表面形貌，建立了不同潤滑狀態下的滑動粗糙界面模型，研究了界面載荷分配、油膜厚度、摩擦因數和臨界速度隨界面形貌參數、法向載荷、潤滑油屬性參數的變化規律；楊紅平等[14]基于分形幾何理論和接觸力學理論，用分形理論表征粗糙表面微凸體參數，建立各變形階段微凸體的接觸剛度模型。

上述學者都通過單微凸體接觸來描述粗糙表面的變形，都假設粗糙表面的微觀形態是許多微凸體附著在一個基體上，且只考慮了接觸過程中微凸體的變形情況，并沒有考慮微凸體所附著的基體的變形。當表面微凸體的材料和基體材料相同時，GW接觸模型和它的擴展改進模型通常能有效地分析接觸行為。然而，對于一些有硬涂層的粗糙表面，其接觸特性主要由軟基體上的硬涂層的材料決定，較軟基體的變形對粗糙表面變形也有相當大的影響[15]。傳統的GW接觸模型它的擴展改進模型不能完全描述它的接觸行為。Shi等[16]測量了一種磁記錄材料的接觸剛度，測得剛度值大約是106N/m，而通過拓展GW模型計算得出的接觸剛度值為107～108N/m。測量得到的接觸剛度測量值比拓展GW模型得到的計算值低，這也證明了基體變形對粗糙表面接觸特性的影響。

為了準確地預測有硬涂層粗糙表面的接觸行為，接觸模型中應包含較軟的基體變形和較硬的微凸體變形。在各位學者研究的基礎上，通過分析較軟基體上有較硬涂層的粗糙表面接觸行為，提出了一個在GW接觸模型基礎上改進的，并同時考慮基體變形和微凸體變形的新彈性接觸模型。通過對新模型進行有限元模擬，并將新模型與有限元分析結果、Hertz接觸模型進行對比，發現新模型可以有效地描述有硬涂層粗糙表面的真實接觸情況。

1 接觸過程中的接觸變形和應力分布

如圖1所示，當兩個球體發生軸對稱接觸時，R1和R2為上下兩個半球體的曲率半徑；E1和E2為上下兩個半球體材料的彈性模量；υ1和υ2為上下兩個半球體材料的泊松比。假設沒有切向滑移，最大的剪切應力將出現在較軟材料球體的接觸表面上。在彈性變形范圍內，隨著載荷的增加，最大剪應力的位置會離接觸表面越來越遠。當接觸球體有較大的曲率半徑和較高的屈服強度或硬度時，最大剪應力值將更大，其位置也將距接觸表面更遠。如果材料表面涂有硬涂層，接觸時應力分布將會不均勻，簡單的Hertz接觸模型無法準確地描述這種接觸。

圖1 球體-球體接觸圖

為證明這個觀點，將硬涂層結構簡化為一個球面微凸體和一個平整基體結合，并與一個剛性平面進行接觸，具體形式如圖2所示，采用ABAQUS有限元軟件對新接觸模型進行分析模擬。圖2中微凸體的半徑都取R=15 μm，微凸體的高度都取ha=100 nm，不同的是，一個系統中微凸體和基體的材料是相同的，它們的彈性模量Ea=Eb=100 GPa；另一個系統中微凸體材料較硬于基體材料，即Ea=200 GPa，Eb=100 GPa。為了簡化對比，兩個系統中微凸體和基體材料的泊松比取υa=υb=0.3。

圖2 單微凸體/基體系統

假設微凸體和基體均是彈性體，將微凸體和基體模型都采用軸對稱八節點等參數單元方法進行有限元分析。為了加快分析速度，網格劃分時，接近接觸表面的區域使用較密的面積為1 nm×1 nm的網格，遠離測試區域的網格可以相對粗化一些。單微凸體/基體系統和它接觸的邊界條件如圖2所示。微凸體綁定在基體之上，微凸體和基體之間無切向滑移，微凸體和基體的x軸方向和基體的y軸方向是固定的。微凸體上方的剛性平面上下移動使微凸體和基體發生形變。為避免由邊界條件引起的誤差，將基體的體積設置得比微凸體體積大一些，這樣就形成了一個無限彈性半空間體。

圖3是有限元模擬時，當實際位移量δ=4.5 nm時，微凸體系統所產生的接觸變形和應力分布圖。其中，圖3(a)為微凸體和基體是同種材料的情況，從圖可以看出它的應力場是連續且均勻的散開在微凸體和基體之中。這種情況與傳統的微凸體包含在基體中的Hertz接觸模型一致，從而證明傳統GW模型的球面接觸使用Hertz接觸解法的正確性。圖3(b)為微凸體的材料較硬于基體的材料，從圖中可以看出，微凸體系統的應力場分布不均勻，且最大應力出現在微凸體和基體的交界處。由于微凸體和基體材料特性的差異，Hertz接觸模型無法解決這種復雜情況。因此，當接觸表面微凸體材料比基體材料有較高的強度和硬度時，應該考慮較軟基體對接觸分析的影響。

(a) Ea=Eb=100 GPa

(b) Ea=200 GPa， Eb=100 GPa

從圖3(b)可看出，有限元模擬時，由于假設微凸體和基體之間沒有發生相對滑移，且微凸體和基體的材料不同，所以在加載條件下，微凸體和基體之間的應力分布不連續。但是，本文研究對象由于所施加的負載或位移相對較小，接觸變形均是在彈性范圍內，所以將忽略微凸體和基體之間的相對滑移。

2 單微凸體接觸分析建模

第1節有限元分析結果表明，微凸體和基體的變形對粗糙表面微觀接觸特性的影響不可忽略。考慮到基體材料和基體變形對整體接觸的影響，可以分開計算微凸體和基體的變形量，分別建立微凸體和基體的剛度模型，然后將這兩個剛度模型耦合形成系統總的剛度模型，基體變形量與微凸體/基體系統總變形量之間關系，微凸體剛度與基體剛度之間關系如圖4所示。這種微凸體基體耦合模型得到的接觸結果并不是理論上的精確解，因為它不包括材料的切向接觸，接觸壓力主要取決于微凸體和基體的幾何和物理特性。當微凸體材料比基體材料硬時，這種把微凸體和基體剛度分開計算的方法是體現微凸體和基體變形的一種簡單方法。

圖4 接觸模型原理圖

根據文獻[17-18]，假設粗糙表面的微凸體為球體，其半徑和高度可從粗糙表面的形貌測量得出。微凸體的剛度可以通過Hertz接觸理論獲得其表達式[19]

(1)

需要注意的是，微凸體的法向位移量δa和系統總的位移量δ并不是同一個量，系統總的位移量是微凸體位移量δa和基體位移量δb之和。當上部的微凸體受到載荷時，微凸體和基體的接觸區是一個半徑為rb的圓形區域，rb的值可通過微凸體的曲率半徑R和微凸體的高度ha求得，可表示為

(2)

基體圓形接觸區域內各點受到的壓力各不相同，接觸區域中心位置受到的壓力最大，最大壓力P0為

(3)

距離接觸區域中心的距離r不同，基體表面發生的法向位移不同，當基體受到的最大壓力P0在不斷變化時，基體表面的法向位移為

(4)

根據平衡條件，圓形區域實際接觸力F值和作用在微凸體上的法向接觸力相同。所以，基體與微凸體所接觸圓形觸區域中心點的位移為

(5)

通過式(3)和式(5)，可以得到基體的剛度計算公式為

(6)

大部分的微凸體的球形區域被認為是包含在基體中，所以式(1)中包括了基體接觸變形的影響。基體的影響分別出現在式(1)和式(6)中，因此當微凸體和基體應用相同材料的時候，提出的接觸模型計算出接觸變形量將高于實際接觸變形量。為了得到一個確切的接觸剛度解，需對式(1)進行修改，排除其中基體對微凸體剛度的影響。然而，當微凸體材料比基體的材料硬時，較軟基體的接觸變形在整個接觸過程中占主導地位，因此式(1)使用Hertz解的影響將會很小，最終導致的誤差可以忽略。為了簡化和方便應用，該微凸體的剛度模型使用Hertz解。

當微凸體和基體的剛度值是常數時，等效剛度k可以通過簡單的公式1/k=1/ka+1/kb確定。然而，如式(1)所示，由于ka不是一個定值，而是一個關于微凸體形變量δa的一個方程，所以在求等效剛度之前，應該事先確定獲得微凸體形變量δa的值。

根據力的相互作用原理，可得

k×δ=ka(δa)×δa=kb×δb

(7)

δ=δa+δb

(8)

將式(1)和式(6)代入式(7)式(8)，可以得到δa關于δ的公式為

(9)

通過化簡和轉化，可以將求解δa的函數為

δa=f(δ)

(10)

(11)

將式(11)中的δa替換為與δ、κ、ξ有關的函數g(δ)，式(11)為

(12)

利用不動點迭代法，可以得到式(12)的近似解，式(12)又為

(13)

理論上來說，起始點g0可以是任何值，但是為了避免符號運算造成的復雜性，假設初始點時微凸體和基體擁有相同的幾何形狀和不同材料特性，即ξ=1，從而得出g0為

(14)

因此迭代序列采用以下方式不斷迭代

(15)

式中：n=0，1，…，∞，迭代的次數越多得到的結果越精確。本文只考慮彈性變形的前提下，經過計算發現，只需迭代兩次得到的近似解就非常接近于計算得到的數值解，且誤差不超過5%。因此微凸體的形變量可以近似地表示為一個關于實際位移量δ的函數，其為

(16)

圖5描述了采用數值解法的式(10)和近似解法的式(16)求得的微凸體形變量在各種幾何和物理特性之下的對比。微凸體材料比基體材料軟時κ的值小于1，微凸體材料比基體材料硬時κ的值大于1，如圖5(a)所示，實際位移量不變的情況下，隨著微凸體材料不斷變硬，微凸體的形變量變得越來越小，基體的形變量變得越來越大，這也證明當微凸體材料硬于基體材料時，微凸體/基體系統中基體的變形量占主導地位。由于一個平滑的表面擁有較大的微凸體半徑R和較小的微凸體高度ha，通過式(11)可以得出，表面越光滑，表面粗糙參數ξ的值就越大。如圖5(b)所示，實際位移量不變的情況下，ξ值越大，微凸體的形變量變得越來越小，基體的形變量變得越來越大。從這些結果中，可以清楚地看到，表面涂層較硬時，基體的變形對接觸行為的影響很大，這和第1節中有限元模擬的分析結果一致。從圖5還可以清楚地看出，由式(16)得出的微凸體形變量近似解(虛線)與式(10)得出的微凸體形變量數值解(實線)是基本相同的，所以我們在建模時使用其簡單的近似解。

(a) 微凸體物理特性對形變量的影響

(b) 微凸體幾何特性對形變量的影響

Fig.5 Comparisons of single asperity deformation with applied displacement between the numerical solution and approximate solution

綜合式(1)、式(6)和式(16)，單微凸體系統的等效接觸剛度可以用一個與實際位移量有關的方程表示

(17)

接觸力可以通過式(18)確定

F(δ)=k(δ)×δ

(18)

3 模型驗證與分析

為了驗證第2節所提出接觸模型的正確性，利用有限元軟件模擬一個剛性平面同微凸體和基體的接觸。對于單微凸體模型，微凸體半徑為R，微凸體高度設為剛性面實際位移的三倍，即ha=3δ。其中基體上圓形接觸區域半徑rb可以通過式(3)得到。

表1 四種不同幾何和物理特性的微凸體

如表1所示，4組不同的微凸體高度ha和微凸體平均曲率半徑R代表現實中的不同粗糙表面。R值取5～20 μm，微凸體高度ha的范圍為3～750 nm，微凸體和基體彈性模量的比值κ取1～2，為了簡化計算，材料的泊松比取υa=υb=0.3，其余物理量rb，ξ可以通過公式計算得到。ξ的四個不同的取值代表了四種不同等級的接觸表面。

由于所有節點和單元的接觸的信息可以直接從有限元得到分析結果，所以微凸體和基體的接觸參數可以很容易地確定。將從有限元分析模擬的結果同提出的新的分析模型和Hertz接觸模型進行比較，以確定使用Hertz解決方案的局限性和所提出的模型的有效性。

當ξ=19.37×103m-1/2，κ=1，δ/rb=6.82×10-3時，微凸體/基體系統的應力分布圖如圖6(a)所示，由圖可知應力主要集中在基體內，形變量主要發生在微凸體上。圖6(b)既包含了Hertz模型解(實線)，新模型解(虛線)和有限元模擬的解(°點)，還包含了在不同的實際位移下，微凸體/基體系統應力應變。正如預期的那樣，在這種情況下，微凸體和基體的材料相同時(κ=1)，有限元分析結果與Hertz解一致，而新提出模型預測的接觸力與Hertz解不一致的原因是在第2節中討論過的位移的雙重計數的原因。

(a) δ/rb=6.82×10-3時的應力分布

(b) Hertz解、新模型解和有限元模擬解的比較

Fig.6 Stress distribution diagram and the comparison of the model results atξ=19.37×103m-1/2，κ=1

圖7(a)表明，當ξ=8.66×103m-1/2，κ=1.25，δ/rb=0.011時，微凸體/基體系統的應力分布圖，從圖中可以看出大部分的應力依舊集中在基體內，但是還是有一部分應力出現在微凸體內。從圖7(b)無量綱的位移和無量綱的接觸力的關系中可以看出，當δ/rb≤0.003 5時，代表有限元分析解，Hertz解，新模型解的三條曲線基本吻合；在δ/rb≥0.003 5時，有限元分析的解開始與Hertz解開始有差別，并慢慢趨近于新提出的分析模型。很顯然，在δ/rb的值較高時，新提出的模型優于Hertz模型。

(a) δ/rb=0.011時的應力分布

(b) Hertz解，新模型解，有限元模擬解的比較

Fig.7 Stress distribution diagram and the comparison of the model results atξ=8.66×103m-1/2，κ=1.25

(a) δ/rb=0.023時的應力分布

(b) Hertz解，新模型解，有限元模擬解的比較

Fig.8 Stress distribution diagram and the comparison of the model results atξ=3.36×103m-1/2，κ=1.5

圖8(a)表明，當ξ=3.66×103m-1/2，κ=1.5，δ/rb=0.023時，微凸體/基體系統的應力分布圖，從圖中可以看出，雖然最大的應力依舊集中在基體內，但是有相當大的應力出現在微凸體內。圖8(b)中的無量綱位移和無量綱接觸力關系圖也很好地顯示出了有限元分析解的趨勢，在δ/rb≤0.003時，有限元分析解，新模型解同Hertz解是基本相同的。但是δ/rb≥0.003之后，有限元分析解和新模型解開始同時偏離Hertz解，并且有限元分析解和新模型解的前進趨勢是幾乎相同。

圖9(a)表明，當ξ=1.24×103m-1/2，κ=2，δ/rb=0.055 3時，微凸體/基體系統的應力分布圖，從圖中可以看出，當實際位移量最大時，最大應力同時出現在微凸體和基體中，但是從圖9(b)中三張代表變形過程的圖可以看出，最大應力首先出現在微凸體中，隨著位移量的不斷增大，最大應力慢慢向下傳遞到了基體中。這說明將微凸體和基體的變形分開研究對于粗糙表面接觸研究非常重要。從圖9(b)無量綱位移和無量綱接觸力關系可知，粗糙度R=15 μm與R=5 μm、10 μm和15 μm 3種情況類似，在δ/rb值很小的時候，Hertz解和新模型解都很好的與有限元分析結果相吻合，隨著δ/rb值的變大之后，有限元分析解和新模型解開始同時偏離Hertz解。

(a) δ/rb=0.055 3時的應力分布

(b) Hertz解，新模型解，有限元模擬解的比較

Fig.9 Stress distribution diagram and the comparison of the model results atξ=1.24×103m-1/2，κ=2

通過4種情況分析對比，可以發現當實際位移δ的值很小時，Hertz解和新提出的模型解都和有限元分析解吻合，所以可以推測當基體的彈性模量或者微凸體和基體的復合彈性模量很小時，Hertz模型解法可以代替新提出的模型。當實際位移δ的值不斷增大時，Hertz解法就無法準確的表示真實地接觸情況。新提出的模型可以有效地描述表層比基體硬的接觸面的真實接觸情況。

4 結 論

(1) 當粗糙表面有硬涂層或硬薄膜時，基體材料和微凸體材料的不同將會影響粗糙表面的接觸特性。研究表明，基體的接觸變形在微凸體接觸研究中不可忽略。

(2) 基于GW模型和Hertz接觸理論，分別計算微凸體和基體的變形量，從而建立微凸體和基體的剛度模型，通過耦合的剛度模型提出了考慮基體變形的微凸體彈性接觸模型。

(3) 利用Hertz接觸理論求得了微凸體變形量與微凸體/基體系統總變形量之間的函數關系，并應用不動點迭代法驗證得到誤差不超過5%的簡單函數關系。實際位移量不變的情況下，隨著微凸體材料不斷變硬，微凸體的形變量越來越小；隨著微凸體不斷變小，微凸體的形變量越來越小。

(4) 利用Hertz模型和有限元分析結果驗證了新模型的正確性。在變形量很小的時候，Hertz模型解和新模型解都很好的與有限元分析結果相吻合，隨著變形量的變大，有限元分析解和新模型解開始同時偏離Hertz模型解，且有限元分析解和新模型解的前進趨勢相同。證明新提出的接觸模型適用于有硬涂層的接觸表面。

(5) 對有限元模型仿真結果進行分析，表明當基體材料和微凸體材料不同時，微凸體/基體系統的應力分布會不均勻，微凸體表面的接觸力比材料相同時的接觸力小，最大應力比材料相同時的最大應力大。

[1] BARBER J R，CIAVARELLA M. Contact mechanics[J]. International Journal of Solids and Structures，2000，37：29-43.

[2] GREENWOOD J A，WILLIAMSON J B P. Contact of nominally flat surfaces[J]. Mathematical and Physical Sciences，1966，295(1442)：300-319.

[3] KOGUT L，ETSION I. Elastic-plastic contact analysis of a sphere and a rigid flat[J]. Journal of Applied Mechanics, 2002，69(5)：657-662.

[4] CHANG W R，ETSION I，BOGY D B. An elastic-plastic model for the contact of rough surfaces[J]. ASME Journal of Tribology，1987，109(2)：257-263.

[5] 趙永武，呂彥明，蔣建忠. 新的粗糙表面彈塑接觸力學模型[J]. 機械工程學報，2007，43(3):95-101.

ZHAO Yongwu，Lü Yanming，JIANG Jianzhong. New elastic-plastic model for the contact of rough surfaces[J]. Chinese Journal of Mechanical Engineering，2007，43(3)：95-101.

[6] JENG Y R，PENG S R. Elastic-plastic contact behavior considering asperity interactions for surfaces with various height distributions[J]. ASME Journal of Tribology，2006，128(2)：665-697.

[7] ALMQVIST A，SAHLIN F，LARSSON R，et al. On the dry elastic-plastic contact of nominally flat surfaces[J]. ASME Journal of Tribology，2007，40(10)：574-579.

[8] 溫淑花，張學良，武美先，等. 結合面法向接觸剛度分形模型建立與仿真[J]. 農業機械學報, 2009，40(11)： 197-202.

WEN Shuhua，ZHANG Xueliang，WU Meixian，et al. Fractal model and simulation of normal contact stiffness of joint interfaces and its simulation[J]. Transactions of the Chinese Society for Agricultural Machinery，2009，40(11)：197-202.

[9] 張學良，溫淑花，徐格寧，等. 結合部切向接觸剛度分形模型研究[J]. 應用力學學報，2003，20(1)：70-72.

ZHANG Xueliang，WEN Shuhua，XU Gening，et al. Fractal model of the tangential contact stiffness of machined surfaces in contact[J]. Chinese Journal of Applied Mechanics，2003，20(1)：70-72.

[10] 尤晉閩，陳天寧. 結合面法向動態參數的分形模型[J].西安交通大學學報，2009，43(9)：91-94.

YOU Jinmin，CHEN Tianning. Fractal model for normal dynamic parameters of joint surfaces[J]. Journal of Xi’an Jiao tong University，2009，43(9)：91-94.

[11] 尤晉閩，陳天寧. 基于分形接觸理論的結合面法向接觸參數預估[J]. 上海交通大學學報，2011，45(9)： 1275-1280.

YOU Jinmin，CHEN Tianning. Estimation for normal parameters of joint surfaces based on fractal theory[J]. Journal of Shanghai Jiaotong University，2011，45(9)：1275-1280.

[12] 李玲，蔡安江，蔡力鋼，等. 螺栓結合面微觀接觸模型[J]．機械工程學報，2016，52(7): 205-212．

LI Ling，CAI Anjiang，CAI Ligang，et al．Micro-contact model of bolted-joint interface[J]. Journal of Mechanical Engineering，2016，52(7): 205-212．

[13] 肖會芳，楊荃，邵毅敏，等. 潤滑狀態下線接觸滑動粗糙界面的動摩擦特性研究[J]. 振動與沖擊，2016，35(1)：188-194.

XIAO Huifang，YANG Quan，SHAO Yimin，et al. Dynamic friction characteristics of sliding rough interfaces in line contact under lubrication[J]. Journal of Vibration and Shock，2016，35 (1): 188-194.

[14] 楊紅平，傅衛平，王雯，等. 基于分形幾何與接觸力學理論的結合面法向接觸剛度計算模型[J]. 機械工程學報，2013，49(1)：102-107.

YANG Hongping，FU Weiping，WANG Wen，et al. Calculation model of the normal contact stiffness of joints based on the fractal geometry and contact theory[J].Journal of Mechanical Engineering， 2013， 49(1)：102-107.

[15] IIDA K，ONO K. Design consideration of contact/near-contact sliders based on a rough surface contact model[J]. ASME Journal of Tribology， 2003，125(3)：1518-1525.

[16] SHI X，ANDREAS A. Investigation of contact stiffness and contact damping for magnetic storage head-disk interfaces[J]. ASME Journal of Tribology，2008，130(2)：843-854.

[17] ZHAO Y W，DAVID D M，CHANG L. An asperity micro contact model incorporating the transition from elastic deformation to fully plastic flow[J]. ASME Journal of Tribology，2000，122(2)：86-93.

[18] 田紅亮，鐘先友，秦紅玲，等. 依據各向異性分形幾何理論的固定結合部法向接觸力學模型[J]. 機械工程學報，2013，49(21)：108-122.

TIAN Hongliang，ZHONG Xianyou，QIN Hongling，et al. Normal contact mechanics model of fixed joint interface adopting anisotropic fractal geometrical theory[J]. Journal of Mechanical Engineering，2013，49(21)：108-122.

[19] JOHNSON K L. Contact mechanics[M]. Cambridge: Cambridge University Press，1985.