舒小敏， 李 堅(jiān)

(中國(guó)航發(fā)湖南動(dòng)力機(jī)械研究所，株洲 412002)

1 引 言

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

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

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

利用邊界元法求解三維摩擦接觸問(wèn)題，其中一個(gè)關(guān)鍵點(diǎn)在于如何確定滑移方向。即當(dāng)出現(xiàn)相對(duì)滑移時(shí)，滑移方向如何確定。不同于二維摩擦接觸問(wèn)題，滑移方向在一條線上，只有兩個(gè)切向滑移方向。三維中切向滑移方向位于一個(gè)平面內(nèi)，難以確定。

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

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

2 邊界積分方程

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

((P,Q)∈Γ)

(1)

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

(2)

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

(3)

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

3 三維庫(kù)侖摩擦定律

在摩擦接觸問(wèn)題中，常采用庫(kù)侖摩擦定律來(lái)描述界面的摩擦現(xiàn)象，如圖1所示。

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

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

(4)

式中μ為摩擦系數(shù)，tτ為切向面力tτ1和tτ2合力，tn為法向面力，如圖2所示。只要切向面力tτ足夠小，小于摩擦極限，則出現(xiàn)粘結(jié)現(xiàn)象，如圖1的點(diǎn)α。其相對(duì)切向位移為零，

(5)

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

當(dāng)摩擦力大于摩擦極限，則會(huì)出現(xiàn)滑移現(xiàn)象。此時(shí)，切向面力tτ的大小等于μ|tn|，其方向與切向相對(duì)運(yùn)動(dòng)方向相反，如圖1的點(diǎn)β。切向摩擦力與切向相對(duì)運(yùn)動(dòng)方向相反，可表示為

(6)

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

(7)

4 接觸約束方程

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

圖3 點(diǎn)對(duì)面施加接觸約束Fig.3 Contact constraints by node -to -surface method

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

(8a)

(8b)

(8c)

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

(9a)

(9b)

(9c)

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

(10a)

(10b)

(10c)

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

(11a)

(11b)

(11c)

(12a)

(12b)

(12c)

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

5 滑移方向預(yù)測(cè)技術(shù)

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

5.1 無(wú)摩擦預(yù)測(cè)技術(shù)

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

(13)

不過(guò)，該方法存在一個(gè)問(wèn)題，即在接觸面為平面時(shí)可能會(huì)失效。因?yàn)闊o(wú)摩擦接觸時(shí)，只有法向方向施加位移約束式(9a)，而切向約束式(9b，9c)為摩擦力為零約束，結(jié)果導(dǎo)致接觸體位移約束不足，計(jì)算結(jié)果失真。如圖4上面的接觸體，只有法向位移(z方向)受到位移約束，切向沒(méi)有位移約束(x和y方向)。即使最后方程可解，得到的位移也將失真，導(dǎo)致滑移方向預(yù)測(cè)失真。因此，該方法的應(yīng)用受到一定的限制，必須要保證位移約束條件足夠充分(限制剛體位移)才能使用。

圖4 接觸面為平面的接觸問(wèn)題Fig.4 Contact problem with a flat contact surface

5.2 無(wú)限摩擦預(yù)測(cè)技術(shù)

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

(14)

本文采用該方法進(jìn)行滑移方向預(yù)測(cè)。

6 摩擦接觸問(wèn)題算法流程

6.1 摩擦接觸迭代收斂條件

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

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

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

6.2 算法流程

三維摩擦接觸算法流程如圖5所示。在該算法中，根據(jù)壓應(yīng)力條件、非貫穿條件及接觸點(diǎn)摩擦狀態(tài)(粘結(jié)或滑移)是否變化，來(lái)判斷當(dāng)前增量加載是否收斂。由于增量步加載，當(dāng)前步長(zhǎng)不一定保證收斂，需判斷是否需要細(xì)分。圖5的迭代次數(shù)是否超過(guò)n，就是為了解決該問(wèn)題。當(dāng)?shù)螖?shù)超過(guò)一定次數(shù)時(shí)，進(jìn)行增量步細(xì)分。當(dāng)所有增量載荷加載完，整個(gè)計(jì)算完成。本文在每個(gè)增量步加載前，采用無(wú)限摩擦預(yù)測(cè)技術(shù)得到滑移方向。

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

7 數(shù)值算例

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

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

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

圖6 網(wǎng)格離散模型Fig.6 Discrete model

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

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

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

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

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

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

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

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

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

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

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

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

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

圖14 網(wǎng)格離散模型Fig.14 Discrete model

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

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

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

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

此外圖18～圖20中，箭頭方向表示沖頭接觸面中接觸點(diǎn)的相對(duì)滑移方向?？梢钥闯觯还苣Σ料禂?shù)如何變化，其都具有對(duì)稱(chēng)性。這是由于 圖13 的計(jì)算模型對(duì)稱(chēng)，必然要求計(jì)算滑移方向具有對(duì)稱(chēng)性。圖18～圖20中邊界元法滑移方向的對(duì)稱(chēng)性，表明本文無(wú)限摩擦滑移方向預(yù)測(cè)技術(shù)的有效性。

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

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

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

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

8 結(jié) 論

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