多尺度粗糙接觸剛度及數(shù)值等效方法研究

馬奔奔 李源 趙偉 郭旭 劉暢

（1 東方電氣風(fēng)電股份有限公司，研發(fā)中心，德陽 618000；2工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連理工大學(xué)工程力學(xué)系，大連 116023）

0 引言

粗糙界面之間接觸剛度的量化對很多工程裝配體問題的研究都起到至關(guān)重要的作用[1-3]，如截面不對稱的大型汽發(fā)轉(zhuǎn)子槽楔附加剛度穩(wěn)定性分析、軸承座軸瓦之間接觸剛度對系統(tǒng)振動(dòng)的影響等。然而，界面之間接觸剛度的量化依然是較為困難的，主要是因?yàn)榻佑|界面上具有很多微觀的粗糙峰，且這些粗糙峰具有多尺度和分形的特征（圖1），即大粗糙峰上會(huì)分布著很多更加細(xì)小的粗糙峰。為解決這個(gè)難題，有學(xué)者利用仿真[4,5]和試驗(yàn)[6]方法對接觸剛度進(jìn)行研究。但是試驗(yàn)和仿真結(jié)果具有單一性，往往只適用于當(dāng)前粗糙樣件。目前主流的粗糙接觸理論主要有兩類：1）Greenwood和Williamson為代表的多峰接觸G-W模型[7]；2）Persson提出的基于放大率的壓力擴(kuò)散理論[8]。G-W模型具有清晰的物理表征，它假設(shè)所有的粗糙峰都可以簡化為高度符合高斯分布，且變形相互獨(dú)立、具有相等曲率半徑的球型粗糙峰，而整體接觸面上的接觸狀態(tài)可通過對每個(gè)粗糙峰的結(jié)果進(jìn)行求和得到。但是G-W模型由于受到單一尺度思想的限制（忽略了粗糙度的多尺度特征和峰間的變形耦合作用），所以被普遍認(rèn)為在小壓縮量的時(shí)候比較準(zhǔn)確[6,9-11]。相反，Persson提出了粗糙接觸理論則是基于放大率的壓力擴(kuò)散方程的一系列數(shù)學(xué)推導(dǎo)[8]，考慮了粗糙表面的多尺度特征。但是受到壓力擴(kuò)散方程中全接觸假設(shè)的影響，Persson模型往往在大壓縮量下被認(rèn)為比較精確[6,12-14]。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，粗糙接觸界面的仿真方法也得到了飛速的發(fā)展。Liu等人[5]采用（法向和切向）彈簧單元來模擬拉桿轉(zhuǎn)子軸承系統(tǒng)中兩個(gè)圓盤之間的界面接觸作用，彈簧單元?jiǎng)偠纫蕾嚠?dāng)前狀態(tài)下的接觸應(yīng)力。雖然Liu等人的彈簧單元等效方法是有創(chuàng)意且有效的，但不容忽視的是該方法對有限元模型網(wǎng)格劃分的質(zhì)量有著較強(qiáng)的要求，接觸界面網(wǎng)格必須嚴(yán)格匹配來保證彈簧單元法向方向不發(fā)生偏斜。嚴(yán)格的網(wǎng)格匹配條件還造成了大量的模型細(xì)節(jié)簡化，影響結(jié)果的準(zhǔn)確性。此外，大量的彈簧單元也意味著耗時(shí)的建模過程。

為了能夠在整個(gè)壓縮過程中預(yù)測粗糙界面之間的接觸狀態(tài)，近期郭旭等人提出了一個(gè)基于放大率的多峰接觸MBMA模型[15]。該模型融合了Persson的粗糙表面放大率量化思想和G-W的球型粗糙峰表征方式的優(yōu)點(diǎn)，考慮了粗糙表面所具有的多尺度和峰間變形耦合的特性，與較多試驗(yàn)[6,16-18]結(jié)果相吻合。MBMA模型對本文的研究提供了新的視角，本文將針對工程分析中常見的粗糙表面形貌特征，利用MBMA模型來對粗糙界面接觸剛度進(jìn)行量化研究，使相關(guān)工程設(shè)計(jì)人員通過簡單查圖表便可以得到較為精確的接觸剛度數(shù)值。為了能夠精確高效地模擬粗糙接觸問題，本文將粗糙接觸界面等效為一層正交各向異性材料薄層，通過調(diào)整薄層的材料屬性來精確控制接觸界面之間法向和切向接觸剛度，以此達(dá)到模擬真實(shí)粗糙界面之間接觸特性的目的。通過在薄層與部件之間施加“綁定”約束，可以進(jìn)一步地消除接觸非線性，極大提高計(jì)算效率。

本文在第一章將首先明確接觸剛度的物理意義以及粗糙表面的形貌特征，并簡要介紹粗糙接觸力學(xué)MBMA模型。第二章將基于MBMA模型對整個(gè)壓縮過程中的界面接觸剛度的量化進(jìn)行研究。一種高效的受壓粗糙接觸界面等效方法及精確性驗(yàn)證將在第三章進(jìn)行介紹。本文的總結(jié)部分放在最后一章。

1 粗糙接觸剛度的理論求解

本章將首先明確法向和切向接觸剛度的物理意義和表征方式。隨后，通過引入粗糙接觸MBMA模型來量化接觸剛度的求解方法。最后，將針對粗糙表面的形貌特點(diǎn)及其表征方式，歸納影響接觸剛度的粗糙表面形貌參數(shù)。

1.1 接觸剛度的求解方法

對于如圖1所示的粗糙接觸問題。接觸界面法向剛度kn可以計(jì)算為[4,19,20]

圖1 粗糙表面接觸狀態(tài)示意圖 Fig.1 The schematic of rough contact problem

其中，F(xiàn)n為法向載荷；Δ為接觸界面間平均間距u的壓縮量，Δ與u的和為常數(shù)；A0為名義接觸面積；σ為平均接觸壓應(yīng)力。切向接觸剛度kt可以計(jì)算為[5]

其中，F(xiàn)t為切向載荷，δ為接觸界面間切向的相對位移，τ為平均接觸切應(yīng)力。實(shí)際上，切向接觸剛度kt對應(yīng)的是靜摩擦接觸狀態(tài)，沒有相對滑移發(fā)生?；诖植诒砻媸歉飨蛲缘淖苑律浞中伪砻孢@一假設(shè)，Persson等人[19]和Akarapu等人[20]都得到了下述結(jié)論

對于常見金屬材料，泊松比ν≈0.3，所以kt≈0.82kn。因此，隨后本文所提到的接觸剛度將都只是法向剛度kn，切向接觸剛度kt的數(shù)值讀者可以根據(jù)式(3)得到。

1.2 基于放大率的多峰接觸MBMA模型

根據(jù)1.1節(jié)可知，法向載荷Fn與接觸界面間的平均間距u直接關(guān)系到界面法向剛度kn的求解。為了求解Fn和u這兩個(gè)參數(shù)，我們這里采用的是郭旭等人提出的一個(gè)基于放大率的多峰接觸MBMA模型[15]。該模型在整個(gè)壓縮過程中，都可以精確地、定量地描述真實(shí)接觸狀態(tài)。

MBMA模型的建立基于真實(shí)的接觸是發(fā)生在一系列“接觸孤島”上的，在這些接觸孤島上分布著許多更細(xì)小的粗糙峰這一物理現(xiàn)象（如圖1所示）。所以在當(dāng)前壓縮量下，只有尺寸在某一特定范圍內(nèi)的粗糙峰對該接觸孤島的接觸狀態(tài)起到?jīng)Q定性作用。這是因?yàn)樾〕叽绲拇植诜鍖⒅苯颖粔浩?，而大尺寸的粗糙峰很難發(fā)生變形。根據(jù)Greenwood和Tripp[21]的理論便可以解出上述接觸孤島的接觸情況，然后整個(gè)的接觸狀態(tài)可以通過對每個(gè)接觸孤島進(jìn)行求和得出[15]。

整體法向載荷Fn可以寫為

式(4)和式(5)的具體求解方法和參數(shù)物理意義，建議參考論文[15]，我們這里不再詳細(xì)展開。最后，將式(1)、式(4)和式(5)聯(lián)立便可以解得接觸界面法向剛度kn。

1.3 粗糙表面的形貌特點(diǎn)及表征方式

式(4)中整體法向接觸載荷Fn的數(shù)值實(shí)際上與粗糙表面的形貌特征有著密切的聯(lián)系。假設(shè)粗糙表面的形貌特征是各向同性的，那么粗糙表面粗糙峰的密度Dsum，粗糙峰峰頂平均的球半徑R和粗糙峰頂點(diǎn)高度分布的標(biāo)準(zhǔn)差σ可以分別寫為

這里α被稱為是Nayak參數(shù)，并滿足

其中，mn稱為Nayak矩。通過假定粗糙表面是各向同性的自仿射分形表面[22,23]，Nayak矩可以根據(jù)式（8）和式（9）求得

和

在式(9)中，H是Hurst指數(shù)；hrms是表面粗糙度的均方根；q0是長程截?cái)嗖〝?shù)（與最大尺度的粗糙峰相關(guān)）；q1是短程截?cái)嗖〝?shù)（與試驗(yàn)可以觀測到的最小尺度的粗糙峰相關(guān)），其數(shù)值等于π除以當(dāng)前的取樣間距。實(shí)際測得的波數(shù)q的范圍依賴于當(dāng)前的取樣方案（如圖2所示）。根據(jù)上述分析可知，粗糙表面的形貌參數(shù)主要由H、hrms、q0和q1這四個(gè)變量決定。因此，上述四個(gè)變量將直接影響到界面接觸剛度kn的數(shù)值評估。

圖2 自仿射分形表面的功率譜密度函數(shù)[11, 24] Fig.2 Power spectral density function of a detected self-affine surface

2 數(shù)值結(jié)果

本章將基于MBMA模型對整個(gè)壓縮過程中的界面接觸剛度的量化進(jìn)行研究。根據(jù)1.3節(jié)分析可知，H、hrms、q0和q1這四個(gè)變量將直接影響接觸界面法向剛度kn的預(yù)測結(jié)果。因此，本章將為上述四個(gè)變量選取盡量大的取值范圍，涵蓋現(xiàn)實(shí)中常見的粗糙表面形貌特征，便于讀者從本文圖表直接量化其所需的接觸剛度數(shù)值。

2.1 q1對接觸剛度的影響

我們知道，粗糙表面上最小尺度的粗糙峰甚至可能為納米級[11,24]。但是，受到當(dāng)前測量儀器分辨率的限制，很難探測出粗糙面上所有尺度的粗糙峰。因此，本小節(jié)我們將分析能夠探測到的截?cái)嗖〝?shù)q1對接觸剛度的影響。這里粗糙表面的形貌參數(shù)定為：Hurst指數(shù)H=0.7；設(shè)定量q0hrms=0.1；假定能夠探測到的最大的波數(shù)q1分別為最小波數(shù)q0的300，600和900倍。在圖3中，我們給出了不同壓應(yīng)力σ下，粗糙接觸界面之間的接觸剛度?？梢钥闯瞿軌蛱綔y到的q1對應(yīng)的最小尺度粗糙峰幾乎不影響接觸剛度的數(shù)值。這是因?yàn)樾〕叨鹊拇植诜宄叽绾苄。欢鴫簯?yīng)力的數(shù)值卻又很高，所以越細(xì)小的粗糙峰越容易發(fā)生屈服而被壓平，對整體剛度起到很小的貢獻(xiàn)。

圖3 不同能夠探測到的最大的波數(shù)q1下，粗糙接觸界面之間的接觸剛度與壓應(yīng)力的關(guān)系曲線 Fig.3The relationship between contact stiffness and contact stress under different value of detected q1

2.2 H對接觸剛度的影響

常見粗糙表面的 Hurst指數(shù)一般滿足H≥0.7[23]。因此，我們這里粗糙表面的形貌參數(shù)設(shè)定為：q0hrms=0.1；能夠探測到的最大的波數(shù)q1為最小波數(shù)q0的300倍。Hurst指數(shù)H分別取0.7，0.75，0.8，0.85，0.9。不同壓應(yīng)力σ下，粗糙接觸界面之間的接觸剛度結(jié)果如圖4(a)和圖4(b)所示?？梢钥闯?，H取0.7時(shí)，接觸剛度的數(shù)值稍微高于H取0.9。實(shí)際上，分形維數(shù)Df與Hurst指數(shù)H有著下述關(guān)系，即Df=3-H。因此，H的數(shù)值越小，表明粗糙峰越“高聳”，因此接觸剛度也越大。

圖4 不同Hurst指數(shù)H下，粗糙接觸界面之間的接觸剛度與壓應(yīng)力的關(guān)系曲線 Fig.4 The relationship between contact stiffness and contact stress under different value of H

但總的來看，對于Hurst指數(shù)H≥0.7的情況，接觸剛度數(shù)值變化較小，可以近似認(rèn)為Hurst指數(shù)H基本上不影響粗糙表面之間的接觸剛度。

2.3 q0hrms對接觸剛度的影響

經(jīng)過前面兩個(gè)小節(jié)的分析，我們知道q1和H基本上不影響或者很小影響粗糙接觸界面之間的接觸剛度，所以我們這里將對表面粗糙度的均方根hrms和長程的截?cái)嗖〝?shù)q0的影響進(jìn)行分析。此外，經(jīng)過MBMA理論的研究發(fā)現(xiàn)，粗糙界面接觸狀態(tài)往往與q0hrms的乘積有關(guān)[15]。因此，我們這里將給定q0hrms的乘積為0.01，0.05，0.1，0.5和1，基本涵蓋了常見粗糙表面的參數(shù)跨度。

不同q0hrms乘積下，粗糙接觸界面之間的接觸剛度與壓應(yīng)力的關(guān)系曲線如圖5所示?？梢钥闯鰍0hrms的乘積會(huì)顯著影響接觸剛度的數(shù)值。鑒于針對大多數(shù)經(jīng)過剖光處理的粗糙表面，hrms一般是微米級的。所以，從圖5可以看出當(dāng)壓應(yīng)力1MPa≤σ≤10MPa時(shí)，接觸剛度k的量級約在104～105MPa/mm之間；當(dāng)壓應(yīng)力10MPa≤σ時(shí)，接觸剛度k的量級約在105～106MPa/mm之間。劉意等人[4,5]的有限元仿真結(jié)果與上述接觸剛度k量級區(qū)間的結(jié)論相吻合。

2.4 材料對接觸剛度的影響

實(shí)際上，不僅僅MBMA理論，其他粗糙接觸力學(xué)理論，如G-W模型、Persson理論等均研究發(fā)現(xiàn)，粗糙界面之間的接觸狀態(tài)與無量綱化后的σ/E*密切相關(guān)。其中，1/E*=(1-ν12)/E1+(1-ν22)/E2，為接觸等效楊氏模量[15]。E1、E2和ν1、ν2分別為兩個(gè)接觸體的彈性模量和泊松比。所以對于其他的本文沒有考慮的接觸材料，建議根據(jù)圖5的解析結(jié)果對接觸壓應(yīng)力σ進(jìn)行無量綱化后替換即可。

經(jīng)過本章分析可以得到以下結(jié)論：能夠探測到的最小尺度粗糙縫q1、Hurst指數(shù)H基本上不影響或者很小影響接觸剛度的數(shù)值；然而，粗糙度的均方根為hrms、長程截?cái)嗖〝?shù)q0會(huì)顯著影響接觸剛度的數(shù)值。接觸材料屬性不同，接觸剛度也會(huì)顯著不同，但讀者可以根據(jù)圖5的解析結(jié)果對接觸壓應(yīng)力σ進(jìn)行無量綱化后得到所需接觸材料的接觸剛度?？偟膩碚f，圖5基本涵蓋了工程分析中常見粗糙表面的接觸剛度數(shù)值，設(shè)計(jì)人員可以根據(jù)圖5采用插值法獲得其想要的接觸剛度。

圖5 不同q0hrms乘積下，粗糙接觸界面之間的接觸剛度與壓應(yīng)力的關(guān)系曲線 Fig.5 The relationship between contact stiffness and contact stress under different value of q0hrms

3 受壓粗糙接觸界面等效方法

本章將受壓粗糙接觸界面等效為一層正交各向異性材料薄層，通過調(diào)整薄層的材料屬性來精確控制接觸界面之間法向和切向接觸剛度，以此達(dá)到模擬真實(shí)粗糙界面之間接觸特性的目的。通過在薄層與部件之間施加“綁定”約束，可進(jìn)一步地消除接觸非線性，極大提高計(jì)算效率。

3.1 等效正交各向異性薄層法

對于如圖6所示的受壓粗糙接觸問題，接觸界面可以等效為一層正交各向異性材料薄層。根據(jù)單晶正交各向異性材料的特點(diǎn)，該薄層的本構(gòu)方程寫為

圖6 受壓粗糙接觸界面等效為一層正交各向異性材料薄層 Fig.6 The compressive rough contact interface is equivalent to a thin layer of orthotropic material

其中，Cijkl為柔度矩陣，εij和σkl分別為應(yīng)變矢量和應(yīng)力矢量。對于單晶材料，三個(gè)主軸方向上材料的彈性特性相等，故其柔度矩陣可以寫為

對于圖1所示的工況，

聯(lián)立式(1)、式(2)和式(12)，可得等效的正交各向異性材料薄層的法向接觸剛度kn和切向接觸剛度kt與材料參數(shù)E和G的關(guān)系為

將式（13）形式改寫，可得等效的正交各向異性材料薄層的材料參數(shù)E和G為

通過在薄層與部件之間施加“綁定”約束，可以進(jìn)一步消除接觸非線性，極大提高計(jì)算效率。

3.2 等效薄層法數(shù)值驗(yàn)證

為了驗(yàn)證粗糙接觸界面的等效正交各向異性薄層法向和切向剛度是否精確可靠，我們基于ABAQUS有限元分析軟件建立了一個(gè)厚度僅為0.005mm，長度和寬度遠(yuǎn)大于其厚度值的薄板。在板的厚度方向，僅劃分了一層線性單元網(wǎng)格（C3D8單元）。假設(shè)已知要模擬粗糙接觸界面法向和切向剛度為kn/A0=kt/A0=5×105 MPa/mm，經(jīng)過式(14)，可以算出正交各向異性材料薄層此時(shí)對應(yīng)的材料參數(shù)E=G=2500MPa。

當(dāng)將薄板底面固定，上端面分別施加壓應(yīng)力和切向應(yīng)力為σ=τ=250MPa，根據(jù)式(12)可以解得此時(shí)薄板上下端面的法向壓縮量和切向相對移動(dòng)量分別為Δ=δ=5×10-4 mm。

有限元仿真的結(jié)果如圖7所示，可以看出仿真結(jié)果與理論預(yù)期相吻合，證明了該等效方法的精確性和可靠性。

圖7 薄板下端面固定：a) 上端面施加250MPa壓應(yīng)力；b) 上端面施加250MPa切應(yīng)力 Fig.7 Bottom surface of the thin plate is fixed

4 結(jié)論

本文基于多尺度粗糙接觸力學(xué)MBMA模型解析地對粗糙界面的接觸剛度進(jìn)行了量化研究。量化結(jié)果涵蓋了較寬的粗糙界面形貌參數(shù)范圍，能夠滿足工程設(shè)計(jì)人員獲取接觸壓力從0MPa到1000MPa的接觸剛度預(yù)測需求。通過接觸剛度對粗糙表面形貌參數(shù)的敏感性分析，發(fā)現(xiàn)不改變接觸材料的前提下，能夠探測到的最小尺度粗糙縫q1、Hurst指數(shù)H基本上不影響或很小影響接觸剛度的數(shù)值，可以不予考慮；粗糙度的均方根為hrms、長程截?cái)嗖〝?shù)q0顯著影響接觸剛度的數(shù)值。進(jìn)一步還發(fā)現(xiàn)，接觸材料屬性不同，接觸剛度也會(huì)顯著不同，可以根據(jù)圖5的解析結(jié)果對接觸壓應(yīng)力進(jìn)行無量綱化后得到所需的接觸材料的接觸剛度。

本文圖5所示的接觸剛度數(shù)值基本涵蓋了工程分析中常見粗糙表面的形貌特征，對工程師解決相關(guān)接觸問題研究將提供一定的幫助。隨后本文提出了一種粗糙接觸界面的等效實(shí)現(xiàn)方法。該方法的核心是將粗糙接觸界面等效為一層正交各向異性材料薄層，通過調(diào)整薄層的材料屬性來精確控制接觸界面之間法向和切向接觸剛度。隨后通過在薄層與部件之間施加“綁定”約束，可進(jìn)一步地消除接觸非線性，將非線性接觸問題簡化為線性問題，極大提高計(jì)算效率。最后，本文采用ABAQUS仿真軟件對上述等效方法進(jìn)行了驗(yàn)證，并得到只需一層線性積分單元便可以精確模擬該正交各向異性材料薄層的結(jié)論。