基于邊界元法求解三維彈性摩擦接觸問題

舒小敏， 李 堅

(中國航發(fā)湖南動力機械研究所，株洲 412002)

1 引 言

接觸問題廣泛存在于工程機械領域。在機械中各部件通常以接觸的形式傳遞運動或載荷，而接觸區(qū)常常應力集中，致使機械部件更早地出現疲勞斷裂等問題，降低部件使用壽命，增加安全隱患。為降低這些影響及危害，接觸問題得到大量研究。

早期主要通過理論解析法求解接觸問題[1-3]，但解析法局限于簡單幾何體、簡單加載及簡單邊界條件的接觸問題。實際工程中，接觸問題不僅模型復雜而且受力加載情況多變。工程中很少存在一個完全滿足解析法的理想模型。因此，尋求數值算法求解接觸問題，常用的接觸數值求解算法包括有限元法[4-8]和邊界元法[9-12]。

邊界元法求解三維摩擦接觸問題有30年的歷史。Kong等[13]采用數學規(guī)劃法求解，該方法將庫侖摩擦錐近似為多邊形的金字塔模型。Yamazaki等[14,15]采用罰函數法求解三維摩擦接觸問題。Garrido等[16]采用直接法并利用增量技術求解。此后，Leahy等[17-19]進一步采用直接法并利用二次單元求解了彈性及彈塑性摩擦接觸問題。國內也展開了相關算法研究，如數學規(guī)劃法[20-22]、罰函數法[23]和直接法[24，25]。

利用邊界元法求解三維摩擦接觸問題，其中一個關鍵點在于如何確定滑移方向。即當出現相對滑移時，滑移方向如何確定。不同于二維摩擦接觸問題，滑移方向在一條線上，只有兩個切向滑移方向。三維中切向滑移方向位于一個平面內，難以確定。

邊界元法求解摩擦接觸問題常采用兩種方式確定滑移方向。一種是利用切向面力得到滑移方向；另一種是利用切向相對位移得到滑移方向。在第一次迭代時，文獻[16,25]通常假設物體之間沒有相對滑動，根據切向面力得到滑移方向，而后針對滑移點采用切向相對位移得到滑移方向。文獻[16]采用該方法并不能保證收斂。同樣在文獻[17-19]中第一次迭代，假設物體之間沒有滑動，根據切向面力得到滑移方向，但卻沒有說明在以后的迭代中如何得到相應的滑移方向。此外，文獻[16-19]只考慮了由粘結接觸點變成滑移接觸點，卻沒有再對滑移接觸點進行判斷，也就是一旦接觸點進入滑移就不再對其判斷。針對這一點，文獻[25]采用了摩擦耗能原理來判斷是否仍是滑移點。

針對以上問題，本文采用直接法求解并利用滑移方向預測技術。該預測技術[26]在二維摩擦接觸問題也得到采用，不過其采用無摩擦預測技術。即首先將摩擦接觸問題假設為無摩擦接觸問題，計算得到相應的切向相對位移，該切向相對位移就是以后的滑移方向。不過，假設為無摩擦接觸可能導致預測的滑移方向失真。為此本文提出并采用無限摩擦預測技術，同時采用增量法及摩擦耗能原理識別非滑移點，避免一旦接觸點進入滑移狀態(tài)就不再對其進行判斷的問題。最后，利用數值算例驗證了本文算法的有效性和收斂性。

2 邊界積分方程

對于彈性問題，在不考慮體力的情況下，每個體的邊界積分方程可表示為[27]

((P,Q)∈Γ)

(1)

式中uj和tj分別為位移和面力分量，Ui j和Ti j分別為位移和面力核函數或基本解，ci j(P)為關于邊界布置系數矩陣，其中(i,j=1,2,3)。位移和面力基本解表達式為

(2)

(1-2ν)(njr,i-nir,j)

(3)

式中r為源點P到場點Q的距離，n為場點Q處的外法線，G和v分別為材料的剪切模量和泊松比。

3 三維庫侖摩擦定律

在摩擦接觸問題中，常采用庫侖摩擦定律來描述界面的摩擦現象，如圖1所示。

圖1 三維摩擦庫侖定律Fig.1 3D coulomb friction

其中庫侖摩擦定律對接觸區(qū)切向面力進行了約束，該約束可表示為

(4)

式中μ為摩擦系數，tτ為切向面力tτ1和tτ2合力，tn為法向面力，如圖2所示。只要切向面力tτ足夠小，小于摩擦極限，則出現粘結現象，如圖1的點α。其相對切向位移為零，

(5)

圖2 接觸方向定義Fig.2 Contact problem definition

當摩擦力大于摩擦極限，則會出現滑移現象。此時，切向面力tτ的大小等于μ|tn|，其方向與切向相對運動方向相反，如圖1的點β。切向摩擦力與切向相對運動方向相反，可表示為

(6)

此外，摩擦存在時，需滿足摩擦耗能條件為

(7)

4 接觸約束方程

由邊界積分方程(1)可知，在三維中每個節(jié)點可以配置3個邊界積分方程(i,j=1,2,3)。對于非接觸的節(jié)點，由于邊界條件的存在，只存在3個未知量，3個邊界積分方程，正好可解。然而對于可能接觸區(qū)的節(jié)點，由于面力和位移均為未知量，共有6個變量，卻僅有3個邊界積分方程，系統(tǒng)方程不可解。需在每個接觸節(jié)點上補充3個接觸約束方程，保證系統(tǒng)方程可解。由于摩擦的存在，計算結果和加載歷史相關，本文采用增量法及點對面(node -to -surface)[16,28,29]的方式施加接觸約束。如圖2所示，假設物體A上的點a和物體B上的點b接觸，點b位于單元Q1Q2Q3Q4中(圖3)。則可建立如下的接觸約束補充方程。

圖3 點對面施加接觸約束Fig.3 Contact constraints by node -to -surface method

對于接觸體A的粘結接觸節(jié)點a，接觸約束條件為法向間隙為零，切向力相對位移為零。

(8a)

(8b)

(8c)

對于接觸體A的滑移接觸節(jié)點a，接觸約束條件為法向間隙為零，切向力滿足庫侖摩擦定律。

(9a)

(9b)

(9c)

對于接觸體B的接觸節(jié)點b，接觸約束條件為法向和切向面力相等，

(10a)

(10b)

(10c)

對于非接觸節(jié)點，不管是接觸體A的節(jié)點a，還是接觸體B的節(jié)點b,約束條件皆為所有方向上面力為零。

(11a)

(11b)

(11c)

(12a)

(12b)

(12c)

從式(8～12)可看出，任意接觸狀態(tài)的節(jié)點都可以補充三個約束方程，從而保證最后的系統(tǒng)方程可解。這種根據接觸狀態(tài)補充約束方程，進而直接求解的算法[16-19,24-26]，本文稱為直接法。

5 滑移方向預測技術

采用滑移方向預測技術，可得到一個預測的滑移方向。在以后的迭代過程中若發(fā)生滑移，滑移方向一直取預測技術中得到的滑移方向。本文介紹兩種預測技術。

5.1 無摩擦預測技術

該方法將摩擦接觸問題假設為無摩擦接觸問題，然后將所有載荷進行加載并計算求解。根據計算結果，可求出切向方向的相對位移，該相對位移就是滑移方向。該方法[26]在二維摩擦接觸問題中應用，數值算例取得良好的結果。其中滑移方向或滑移角可表示為

(13)

不過，該方法存在一個問題，即在接觸面為平面時可能會失效。因為無摩擦接觸時，只有法向方向施加位移約束式(9a)，而切向約束式(9b，9c)為摩擦力為零約束，結果導致接觸體位移約束不足，計算結果失真。如圖4上面的接觸體，只有法向位移(z方向)受到位移約束，切向沒有位移約束(x和y方向)。即使最后方程可解，得到的位移也將失真，導致滑移方向預測失真。因此，該方法的應用受到一定的限制，必須要保證位移約束條件足夠充分(限制剛體位移)才能使用。

圖4 接觸面為平面的接觸問題Fig.4 Contact problem with a flat contact surface

5.2 無限摩擦預測技術

該方法將摩擦接觸問題假設為無限摩擦接觸問題，然后將所有載荷進行加載并計算求解。此時法向和切向都是位移約束條件,即式(8a～8c)，不會出現位移約束不足的情況。根據求解結果，可得到相應的切向合力方向，該方向的負方向就是相對滑移方向。其中滑移方向或滑移角方向可表示為

(14)

本文采用該方法進行滑移方向預測。

6 摩擦接觸問題算法流程

6.1 摩擦接觸迭代收斂條件

(1) 對于可能接觸區(qū)的非接觸點，需滿足非貫穿條件，即法向間隙gn≥0。否則非接觸點變?yōu)榻佑|點。

(2) 對于可能接觸區(qū)的接觸點，需滿足壓應力條件，即法向壓力要求非正，即tn≤0。否則接觸點變?yōu)榉墙佑|點。

只有上面四條同時滿足，不出現違反的情況，摩擦接觸求解才算收斂。對于第(4)條的收斂判斷準則，文獻[16-19]均沒有采用。即認為一旦判斷為滑移接觸點，不管以后如何變化，一直為滑移接觸點。事實上這種假設值得懷疑，在二維摩擦接觸中，發(fā)現有部分滑移狀態(tài)[26]，指的是在增量加載過程中，一開始為滑移點，隨后又變?yōu)檎辰Y點。因此，本文將耗能摩擦條件也加入到收斂性判斷的準則中。

6.2 算法流程

三維摩擦接觸算法流程如圖5所示。在該算法中，根據壓應力條件、非貫穿條件及接觸點摩擦狀態(tài)(粘結或滑移)是否變化，來判斷當前增量加載是否收斂。由于增量步加載，當前步長不一定保證收斂，需判斷是否需要細分。圖5的迭代次數是否超過n，就是為了解決該問題。當迭代次數超過一定次數時，進行增量步細分。當所有增量載荷加載完，整個計算完成。本文在每個增量步加載前，采用無限摩擦預測技術得到滑移方向。

圖5 三維摩擦接觸算法流程Fig.5 Flowchart of 3D friction contact algorithm

7 數值算例

為驗證本文數值算法的有效性和收斂性，與有限元法進行比較。在有限元法中采用 surfure -to -surface接觸及拉格朗日乘子法求解[8]。

7.1 立方塊與基座的摩擦接觸

如圖6所示，彈性立方塊和彈性基座接觸。立方塊沖頭的尺寸為10×10×10，基座的尺寸為16×16×10。兩者采用相同的材料參數為彈性模量E=200 GPa，泊松比v=0.3?；牡酌婀潭?，沖頭的上表面受到一個均布壓力p=200 MPa。在邊界元法中采用線性四邊形單元進行離散，有限元法中采用線性六面體單元，并采用不同的摩擦系數μ=0.05，0.1和0.2。

圖6 網格離散模型Fig.6 Discrete model

從圖7～圖9的接觸壓力可以看出，兩種計算方法得到壓力的最大值都隨著摩擦系數的增大而增大，最大相對誤差為4.79%，如圖9所示。此外圖中壓力趨勢基本一致，表明接觸壓力計算結果的可靠性。

圖7 摩擦系數μ=0.05時的接觸壓力Fig.7 Contact pressure when friction coefficient μ=0.05

圖8 摩擦系數μ=0.1時的接觸壓力Fig.8 Contact pressure when friction coefficient μ=0.1

圖9 摩擦系數μ=0.2時的接觸壓力Fig.9 Contact pressure when friction coefficient μ=0.2

隨著摩擦系數從0.05增加到0.2，從圖10～圖12可以看出，粘結的區(qū)域逐漸變大(有限元法的紅色部分，邊界元法的紅點)。

圖10 摩擦系數μ=0.05時的滑移和粘結區(qū)域Fig.10 Slip and stick area when friction coefficient μ=0.05

圖11 摩擦系數μ=0.1時的滑移和粘結區(qū)域Fig.11 Slip and stick area when friction coefficient μ=0.1

圖12 摩擦系數μ=0.2時的滑移和粘結區(qū)域Fig.12 Slip and stick area when friction coefficient μ=0.2

同時兩種方法得到粘結區(qū)的面積基本相同，表明粘結區(qū)和滑移區(qū)的大小判斷有效。

此外圖10～圖12中，箭頭方向表示沖頭接觸面中接觸點的相對滑移方向?？梢钥闯?，不管摩擦系數如何變化，其具有高度的對稱性。這是由于 圖6 的計算模型高度對稱，必然要求滑移方向具有高度的對稱性。圖10～圖12邊界元法滑移方向的對稱性，表明了本文無限摩擦滑移方向預測技術的有效性。

7.2 圓弧形沖頭與基座的摩擦接觸

彈性圓弧形沖頭和彈性基座模型如圖13所示。沖頭的幾何尺寸為16×16×5，基座的幾何尺寸為16×16×8。兩接觸體采用相同的材料參數，彈性模量E=200 GPa，泊松比v=0.3?；牡酌婀潭ǎ瑳_頭的上表面受到一個均布壓力p=200 MPa。在邊界元法中采用線性四邊形單元進行離散，有限元法中采用線性六面體單元離散(圖14)，并采用不同的摩擦系數μ=0.05，0.1和0.2。

圖13 圓弧形沖頭和基座幾何模型Fig.13 Geometric model of arc-shaped punch and base

圖14 網格離散模型Fig.14 Discrete model

從圖15～圖17的接觸壓力可以看出，兩種計算方法得到壓力最大值的相對誤差不超過0.52%。此外，圖中壓力趨勢基本一致，表明接觸壓力計算結果的可靠性。

隨著摩擦系數從0.05增加到0.2，從圖18～圖20可以看出，粘結的區(qū)域逐漸變大(有限元法的紅色部分，邊界元法的紅點)。同時兩種方法得到粘結區(qū)的面積基本相同，表明粘結區(qū)和滑移區(qū)的大小判斷有效。

圖15 摩擦系數μ=0.05時的接觸壓力Fig.15 Contact pressure when friction coefficient μ=0.05

圖16 摩擦系數μ=0.1時的接觸壓力Fig.16 Contact pressure when friction coefficient μ=0.1

此外圖18～圖20中，箭頭方向表示沖頭接觸面中接觸點的相對滑移方向?？梢钥闯?，不管摩擦系數如何變化，其都具有對稱性。這是由于 圖13 的計算模型對稱，必然要求計算滑移方向具有對稱性。圖18～圖20中邊界元法滑移方向的對稱性，表明本文無限摩擦滑移方向預測技術的有效性。

圖17 摩擦系數μ=0.2時的接觸壓力Fig.17 Contact pressure when friction coefficient μ=0.2

圖18 摩擦系數μ=0.05時的滑移和粘結區(qū)域Fig.18 Slip and stick area when friction coefficient μ=0.05

圖19 摩擦系數μ=0.1時的滑移和粘結區(qū)域Fig.19 Slip and stick area when friction coefficient μ=0.1

圖20 摩擦系數μ=0.2時的滑移和粘結區(qū)域Fig.20 Slip and stick area when friction coefficient μ=0.2

8 結 論

邊界元法基于直接法求解三維彈性摩擦接觸問題，其收斂性一直難以保證。針對這一問題，本文采用增量法求解，并利用無限摩擦滑移方向預測技術。該預測技術先將摩擦接觸問題假設為無限摩擦接觸問題，物體之間不發(fā)生滑動，計算得到相應的切向面力，該切向面力的反方向即為滑移方向。在以后的迭代過程中若發(fā)生滑移，滑移方向一直取預測技術中滑移方向。同時采用摩擦耗能條件對進入滑移狀態(tài)的接觸點進行判斷，避免接觸點一旦進入滑移狀態(tài)就不再對其進行判斷的問題。不同摩擦系數下的數值計算結果(接觸壓力、粘結區(qū)滑移區(qū)的大小及滑移方向)都表明本文算法的有效性和收斂性及滑移方向預測技術的有效性。