董玉文,任青文,蘇琴(.重慶交通大學(xué)省部共建水利水運(yùn)工程教育部重點(diǎn)實(shí)驗(yàn)室,重慶400074;.河海大學(xué)土木工程學(xué)院,南京0098)

接觸摩擦問題的擴(kuò)展有限元數(shù)值模擬方法

董玉文1,任青文2,蘇琴1
(1.重慶交通大學(xué)省部共建水利水運(yùn)工程教育部重點(diǎn)實(shí)驗(yàn)室,重慶400074;2.河海大學(xué)土木工程學(xué)院,南京210098)

研究了采用擴(kuò)展有限元法分析接觸摩擦問題的方法。該方法基于單位分解思想,在位移插值形式中引入不連續(xù)函數(shù)和裂尖漸近位移場函數(shù)分別反映接觸面的不連續(xù)性和接觸面端點(diǎn)的奇異性;根據(jù)虛功原理建立有限元支配方程,接觸面不連續(xù)函數(shù)積分采用分區(qū)高斯積分的方法。擴(kuò)展有限元法在分析接觸問題時(shí)不需要引入接觸面單元,接觸面可以位于單元內(nèi)任意位置。因此該方法在分析接觸摩擦問題以及壓剪裂紋的開裂擴(kuò)展時(shí)有顯著的優(yōu)勢。

擴(kuò)展有限元;接觸摩擦;數(shù)值模擬

接觸問題廣泛存在于機(jī)械、土木、水利等工程領(lǐng)域。如混凝土壩縱縫和橫縫面的接觸,壩體與壩基的接觸,面板堆石壩中鋼筋混凝土面板與墊層之間的接觸,巖體節(jié)理面或裂縫面的工作狀態(tài)等。接觸非線性的特點(diǎn)和難點(diǎn)在于接觸邊界和接觸力均事先未知,而初始間隙和摩擦效應(yīng)使問題變得更加困難起來。所以,摩擦接觸行為的分析被認(rèn)為是固體力學(xué)領(lǐng)域中最具有挑戰(zhàn)性的問題之一。目前求解接觸問題的方法以有限元法為主,主要算法有拉氏乘子法[1]、罰函數(shù)法[2,3]、間隙元法[4,5]、數(shù)學(xué)規(guī)劃法[6]等。這些方法都需要引入層狀單元、節(jié)理單元等界面單元來模擬斷層、節(jié)理、裂隙等不連續(xù)面的局部接觸特性,并且在模擬界面接觸演化問題時(shí),隨著界面的演化,需要不斷更新界面單元,重新劃分網(wǎng)格,工作相當(dāng)煩瑣。

擴(kuò)展有限元法(XFEM)是新興的分析不連續(xù)問題的數(shù)值方法,在斷裂力學(xué)中得到了迅速發(fā)展。由于XFEM不要求界面與單元邊界一致,因而可以用來分析接觸問題。與常規(guī)有限元法相比,用XFEM分析接觸問題,有兩個(gè)優(yōu)點(diǎn):①不需要在接觸面上布設(shè)界面單元;②避免了隨著界面的發(fā)展演化重新劃分網(wǎng)格的麻煩。目前,用XFEM分析接觸摩擦問題的研究并不多,只有John Dolbow[7]和A.R. Khoei[8]等分別采用不同的迭代算法研究接觸摩擦問題。本文研究用XFEM結(jié)合增廣Lagrange乘子算法分析接觸問題的方法,為進(jìn)一步采用XFEM分析壓剪斷裂問題奠定了基礎(chǔ)。

1 接觸問題的虛功方程與離散

設(shè)有2個(gè)接觸體Ⅰ和Ⅱ,其可能的接觸邊界為Γ,兩接觸體在該界面上的點(diǎn)對可能存在的初始間隙量設(shè)為g0n(法向的相對距離),界面摩擦系數(shù)設(shè)為μ。該接觸問題的虛功方程為

其中:u為位移;δu為虛位移;b為體力;ˉt為面力; n為接觸面上的法向矢量;tN,tT分別為接觸面的法向與切向接觸力。

方程(1)的形式在具體的離散過程中可以簡化為更簡單的形式。假定2個(gè)接觸體中,Ⅰ為從接觸體(slave body),Ⅱ?yàn)橹鹘佑|體(master body),則方程(1)可以寫為

由于在式(2)中,等式右邊的項(xiàng)在未發(fā)生接觸的邊界上積分為零,而在發(fā)生接觸的邊界上接觸力大小相等,方向相反,因此有

其中,Γc表示兩接觸體發(fā)生接觸的邊界部分。

對求解問題進(jìn)行有限元離散,則上式(3)右邊積分可以用從接觸體上所有接觸節(jié)點(diǎn)的和來表示:

其中,nc表示發(fā)生接觸的節(jié)點(diǎn)數(shù)目;gNi表示節(jié)點(diǎn)i的法向間隙函數(shù);gTi表示節(jié)點(diǎn)i的切向間隙函數(shù)(點(diǎn)xi在時(shí)間步內(nèi)的切向相對運(yùn)動距離)。

2 接觸問題的XFEM模擬

2.1 擴(kuò)展有限元法位移模式的構(gòu)造

常規(guī)有限元法所表示的位移近似解為

其中:Ni為節(jié)點(diǎn)i的插值形函數(shù);ui為節(jié)點(diǎn)i的自由度向量。對于域內(nèi)任意一點(diǎn)(x,y),插值形函數(shù)應(yīng)該滿足單位分解:(x,y)≡1。

Duarte和Oden[9]發(fā)現(xiàn)基于插值形函數(shù)的單位分解性質(zhì),可在常規(guī)有限元表達(dá)式中增加能夠描述局部不連續(xù)特性的附加函數(shù),間接反映不連續(xù)面的存在。基于單位分解思想,Belytschko等[10,11]提出了適合于描述含不連續(xù)面的近似位移插值函數(shù),對于彈性問題,典型的表達(dá)式為

其中:N為網(wǎng)格中所有離散節(jié)點(diǎn)的集合;NΓ是被不連續(xù)面穿過、但不包含其端點(diǎn)的單元節(jié)點(diǎn)集合(圖1中小方框表示的節(jié)點(diǎn));N∧表示含不連續(xù)面端點(diǎn)的單元節(jié)點(diǎn)集合(圖1中小圓圈表示的節(jié)點(diǎn));H(x)為Heaviside函數(shù)[11],是反映不連續(xù)面位移不連續(xù)的附加函數(shù);ai為相應(yīng)的附加自由度,在不連續(xù)面上、下側(cè)H(x)分別取+1和-1,以反映位移不連續(xù)情況,即

其中,式(7)中?(x)表示不連續(xù)面的直線方程。

Φα為不連續(xù)面端點(diǎn)(裂尖)漸進(jìn)位移場附加函數(shù),反映端點(diǎn)附近的應(yīng)力奇異性,bαi為相應(yīng)的裂尖附加自由度;Φα通過線彈性斷裂力學(xué)中裂尖漸進(jìn)位移場提取得到。建立以裂尖為坐標(biāo)原點(diǎn)的(r,θ)極坐標(biāo)系(圖2),則Φα由下列4個(gè)基函數(shù)組成:

圖1 XFEM附加函數(shù)節(jié)點(diǎn)選取Fig.1 Nodes selection with XFEM additional function

圖2 裂尖局部坐標(biāo)系Fig.2 Local coordinate system at the crack tip

2.2 有限元支配方程的形成

對離散位移表達(dá)式(6)求導(dǎo)數(shù),可以得到應(yīng)變?yōu)?/p>

將式(9)和式(6)代入式(1),則可以得到以下有限元支配方程:

式中,d為域內(nèi)單元節(jié)點(diǎn)位移列向量集合,包括節(jié)點(diǎn)常規(guī)自由度和附加自由度,對于常規(guī)單元,= {};對于被不連續(xù)面貫穿的單元節(jié)點(diǎn)i,={,};對于含端點(diǎn)的單元節(jié)點(diǎn)j,=,。f為域內(nèi)單元節(jié)點(diǎn)力列向量集合,由單元節(jié)點(diǎn)力組裝而成,f的形式為KC為式(1)中等式左邊第一項(xiàng)積分后得到的整體剛度矩陣,由單元?jiǎng)偠染亟M裝而成,的形式為

式(15)中矩陣元素的上標(biāo)表示對應(yīng)的附加自由度。式(15)和式(14)中出現(xiàn)的子矩陣和外力矢量分別計(jì)算如下:

KI為式(1)中右邊的含不連續(xù)面Γ項(xiàng)積分后得到的附加剛度矩陣,是由含接觸面的不連續(xù)單元的組裝而成。得到各單元的和以后,按單元節(jié)點(diǎn)編號對號入座,進(jìn)行疊加就可以得到含接觸面紋單元的單元?jiǎng)偠染仃嘖e,然后組裝單元?jiǎng)偠染仃嚲涂梢缘玫秸w剛度矩陣,

2.3 接觸面的積分方法

要計(jì)算接觸面上的積分,必須對接觸面進(jìn)行離散。由于接觸面與有限元網(wǎng)格是相互獨(dú)立的,所以將接觸面分解成許多一維單元,具體的分解方法如圖所示,每個(gè)一維單元由不連續(xù)面與單元的兩個(gè)交點(diǎn)構(gòu)成(如圖3)。為了對式(1)或(4)進(jìn)行數(shù)值積分,在每個(gè)一維單元的兩側(cè)布置2個(gè)Gauss積分點(diǎn)(如圖4所示),組成一對接觸點(diǎn)對,接觸面兩側(cè)任意Gauss點(diǎn)的δu,tN,tT通過離散插值得到。

對于某一接觸點(diǎn)對,由式(6)可以得到兩接觸點(diǎn)之間的相對位移為

圖3 不連續(xù)面數(shù)值積分時(shí)一維子單元生成方法Fig.3 Generation of 1-D subelement on the discontinuous interface for integral

圖4 不連續(xù)面兩側(cè)Gauss積分點(diǎn)布置法Fig.4 Gauss quadrature points on each side of the discontinuous interface

其中Gauss點(diǎn)的位移采用母單元(四邊形單元)節(jié)點(diǎn)位移插值得到,一維單元上Gauss點(diǎn)在母單元中的局部坐標(biāo)的計(jì)算方法為:由一維單元上gauss點(diǎn)的局部坐標(biāo)ξ1D計(jì)算出其整體坐標(biāo)x,然后由整體坐標(biāo)x反推出該gauss點(diǎn)在四邊形母單元中的局部坐標(biāo)(ξ,η),由式(21)就可以得到接觸點(diǎn)對的相對位移。

2.4 計(jì)算接觸力的Lagrange乘子法

將式(21)代入式(4),有

式(22)中的接觸力tN,tT屬于未知量,需要采用增量迭代法計(jì)算。本文采用用增廣Lagrange乘子算法[12],即將tN,tT用Lagrange乘子代替并進(jìn)行迭代。在每一級加載中,都進(jìn)行Lagrange乘子迭代計(jì)算,得到真實(shí)的接觸狀態(tài)與接觸力,并作為下一級加載的初始接觸狀態(tài)。

3 算例分析

圖5所示為一相互接觸的彈性體,尺寸為長L=1.0 m,寬b=1.0 m,彈性模量為E=2×1011Pa,泊松比μ=0.3,摩擦系數(shù)為0.3。彈性體下部固定,上部邊界施加法向分布壓力σ=100 kPa,切向施加分布剪力τ。XFEM模型如圖6所示,采用四節(jié)點(diǎn)等參元,共劃分單元20×25=500個(gè),節(jié)點(diǎn)546個(gè),其中接觸面貫穿單元內(nèi)部。為了分析界面發(fā)生滑移的過程,先施加法向壓力σ,再分級施加切向剪力τ,直到兩彈性體發(fā)生完全滑移結(jié)束加載。

圖5 接觸體模型Fig.5 Model of a contact body

圖6接觸體擴(kuò)展有限元網(wǎng)格Fig.6 XFEM mesh of contact bodies

圖7 是τ=25.2 kPa時(shí)的σx,σy應(yīng)力分布云圖,可以看出接觸面上的切向方向的正應(yīng)力不連續(xù),而法向正應(yīng)力連續(xù),符合接觸分析的接觸面壓應(yīng)力條件;左下端受拉,右下端受壓,與實(shí)際受力情況吻合。

圖7正應(yīng)力分布圖(Pa)(τ=25.2 kPa)Fig.7 Normal and tangential stress contours

圖8 是不同分布剪力載荷下的剪應(yīng)力分布圖,從圖中可以得出以下結(jié)論:①彈性體左下部剪應(yīng)力在剪力載荷較小時(shí)均為負(fù),為壓剪狀態(tài),右下部均為正,為壓剪狀態(tài),與實(shí)際情況吻合;②在剪力荷載較小時(shí),接觸面處于粘結(jié)狀態(tài),接觸面上最大剪應(yīng)力位于接觸面中部,隨著剪力荷載的增大,接觸面逐漸發(fā)生滑移,首先發(fā)生滑移的部位在接觸面剪應(yīng)力較大的部位,即中部,發(fā)生滑移以后,接觸面上最大剪應(yīng)力的位置向后移動。

圖9給出了不同剪力載荷下的接觸面摩擦力分布與接觸狀態(tài)?？梢钥闯黾袅^小時(shí)接觸面處于完全粘結(jié)狀態(tài),最大摩擦力位于接觸面中部;當(dāng)τ逐漸增大時(shí),接觸面會發(fā)生滑移,隨著剪力荷載的增大,發(fā)生滑移的部位逐漸向后推移,最大摩擦力也逐漸向后推移;最終兩接觸體發(fā)生完全滑移。

圖8 不同分布剪力荷載作用下的剪應(yīng)力分布云圖(Pa)Fig.8 The shear stress contours corresponding to different shear loadings

圖9 不同分布剪力荷載作用下的接觸面摩擦力分布于接觸狀態(tài)(Pa)Fig.9 The contact frictional stresses and contact statuses on the interface corresponding to different shear loadings

4 結(jié)語

本文將擴(kuò)展有限元法與增廣Lagrange乘子法結(jié)合,研究了用擴(kuò)展有限元法分析摩擦接觸問題的方法,通過典型算例分析,驗(yàn)證了方法程序的正確性。該方法與常規(guī)有限元法的不同之處是,不需要采用接觸面單元,通過引入附加函數(shù)和附加自由度,使得不連續(xù)面可以位于單元內(nèi)任意位置,避免了常規(guī)有限元分析接觸問題時(shí)對網(wǎng)格劃分的限制。通過Lagrange乘子的迭代,使得乘子逐漸接近于真實(shí)的接觸力和真實(shí)的接觸狀態(tài),達(dá)到接觸分析的目的。

水利工程中巖體、混凝土的破壞形式很多屬于壓剪性開裂破壞,如巖質(zhì)邊坡、拱壩壩肩、重力壩壩踵、壩基深層滑動、板塊內(nèi)部地震和水庫誘發(fā)地震等。由于壓剪性裂紋擴(kuò)展涉及到接觸非線性和斷裂力學(xué)結(jié)合的問題,非常復(fù)雜,實(shí)際上壓剪性裂紋的開裂擴(kuò)展可以分為接觸摩擦和開裂擴(kuò)展兩個(gè)階段,一旦裂紋面的摩擦力過大直至超過材料的抗剪強(qiáng)度,即進(jìn)入開裂擴(kuò)展階段,本文方法為進(jìn)一步采用擴(kuò)展有限元法分析壓剪性裂紋的開裂擴(kuò)展奠定了基礎(chǔ)。

[1]王水林,鄧建輝,葛修潤.改進(jìn)的拉氏乘子法在接觸摩擦問題中的應(yīng)用[J].巖土工程學(xué)報(bào),1998,20(5):64 -67.

[2]KIKUHCI N.A Smoothing Technique for Reduced Integration Penalty Methods in Contact Problems[J].International Journal for Numerical Methods in Engineering, 1982,(18):343-350.

[3]PERIC D,OWEN DRJ.Computational Model for 3D Contact Problems with Friction Based on the Penalty Method[J].International Journal for Numerical Methods in Engineering,1992,(35):1289-1309.

[4]雷曉燕,SWOBODE G,杜慶華.接觸摩擦單元的理論及其應(yīng)用[J].巖土工程學(xué)報(bào),1994,16(3):23-32.

[5]邵偉,金鋒,王光倫.用于接觸面模擬的非線性薄層單元[J].清華大學(xué)學(xué)報(bào),1999,39(2):34-38.

[6]HAUG D,SAXCE G.Frictionless Contact of Elastic Bodies by Finite Element Method and Mathematical Programming Technique[J].Computers Structures,1980, (11):55-67.

[7]DOLBOWJ,MO¨ESN,BELYTSCHKO T.An Extended Finite Element Method for Crack Growth with Frictional Contact[J].Computer Methods in Applied Mechanics and Engineering,2001,(190):6825-6846.

[8]KHOEI A R,NIKBAKHT M.An Enriched Finite Element Algorithm for Numerical Computation of Contact Friction Problems[J].International Journal of Mechanical Sciences,2007,(49):183-199.

[9]DUARTE C A,ODEN J T.An H-P Adaptive Method Using Clouds[J].Computer Methods in Applied Mechanics and Engineering,1996,(139):237-262.

[10]BELYTSCHKO T,BLACK T.Elastic Crack Growth in Finite Elements with Minimal Remeshing[J].International Journal for Numerical Method in Engineering, 1999,(45):601-620.

[11]MO¨ES N.DOLBOW J,BELYTSCHKO T.A Finite Element Method For Crack Growth Without Remeshing [J].International Journal for Numerical Method in Engineering,1999,(46):131-150.

[12]SIMOJC,LAURSEN T A.An Augmented Lagrangian Treatment of Contact Problems Involving Friction[J]. Computers and Structures,1992,42(1):97-116.

(編輯:趙衛(wèi)兵)

Extended Finite Element Method of Frictional Contact Problem

DONG Yu-wen1,REN Qing-wen2,SU Qin1
(1.Key Laboratory of Hydro-Science and Engineering,Chongqing Jiaotong University,Chongqing 400074, China;2.College of Civil Engineering,Hohai University,Nanjing 210098,China)

A new method for modeling frictional contact problem by extended finite element method(XFEM)is researched.On the basis of the idea of partition of unity,a discontinuous function and the two-dimensional asymptotic crack-tip displacement fields are added to the displacement interpolation formula to account for discontinuity of the contact surface and stress singularity near the contact surface tip respectively.Governing equations are deduced by virtual work principle and the integral of the discontinuous functions are calculated by Gauss integral through sub-elements cut by discontinuous surface.Thus it is unnecessary for using contact element in contact analysis by XFEM compared with FEM,and the element can be meshed without considering the location of the contact surface.So XFEM is promising in modeling frictional contact problem and the crack growth of compression-shear mixed mode crack.

extended finite element method;frictional contact;numerical model

TV642.3

A

1001-5485(2009)05-0045-05

2008-10-17;

2008-12-17

國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973項(xiàng)目)基金資助項(xiàng)目(2007CB714104);重慶交通大學(xué)省部共建水利水運(yùn)工程教育部重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目(SLK2008B07);重慶交通大學(xué)校內(nèi)青年科學(xué)基金資助項(xiàng)目(2008(G-08))

董玉文(1974-),男,甘肅武威人,博士、講師,主要從事水工結(jié)構(gòu)的教學(xué)與研究,(電話)13647615309(電子信箱)dongyuwen7402 @yahoo.com.cn。