馬文濤 ，師俊平，李 寧

（1. 西安理工大學(xué) 土木建筑工程學(xué)院，西安 710048；2. 寧夏大學(xué) 數(shù)學(xué)計(jì)算機(jī)學(xué)院，銀川 750021）

1 引 言

接觸摩擦問(wèn)題廣泛存在于機(jī)械、航空航天、土木、水利等工程領(lǐng)域，并且隨著對(duì)工程結(jié)構(gòu)分析精細(xì)程度日益提高的要求，對(duì)該問(wèn)題的準(zhǔn)確描述和求解成為難以回避的問(wèn)題。由于接觸邊界及其接觸狀態(tài)（分離、黏合、滑移）在計(jì)算前均未知，再加之當(dāng)接觸發(fā)生后，需要計(jì)算摩擦力，接觸問(wèn)題本質(zhì)上是復(fù)雜的非線性問(wèn)題，所以接觸摩擦問(wèn)題的求解也被認(rèn)為是固體力學(xué)中極具挑戰(zhàn)性的問(wèn)題之一。目前，有限元是求解接觸問(wèn)題最主要的數(shù)值方法，常用算法包括拉氏乘子法、懲罰函數(shù)法、數(shù)學(xué)規(guī)劃法、互補(bǔ)法等。這些方法都需要引入層狀單元或節(jié)理單元模擬不連續(xù)面的局部接觸特性。但當(dāng)接觸面移動(dòng)時(shí)（如裂紋擴(kuò)展），不得不更新界面單元，網(wǎng)格也隨之需要重新劃分，工作相當(dāng)繁瑣。無(wú)網(wǎng)格法是近年來(lái)發(fā)展速度較快的新型數(shù)值方法，它只需要問(wèn)題域及其邊界上的離散節(jié)點(diǎn)信息，沒(méi)有網(wǎng)格依賴性，非常適合求解流-固耦合、高速碰撞、裂紋擴(kuò)展、局部化等大變形問(wèn)題，受到計(jì)算力學(xué)界和工程界的極大關(guān)注[1]。Belytschko等[2]采用罰函數(shù)法施加接觸條件，基于無(wú)網(wǎng)格Galerkin法（EFGM）模擬了裂紋在壓荷載作用下的擴(kuò)展問(wèn)題。龐作會(huì)等[3]將Goodman單元引入EFGM模擬巖體不連續(xù)面。李臥東等[4]基于EFGM，在接觸面上引入罰參數(shù)，通過(guò)迭代計(jì)算，模擬了彈性體與剛體間的接觸行為。盧波等[5]指出，無(wú)網(wǎng)格法中使用 Goodman單元存在位移模式與假定位移模式不相協(xié)調(diào)的問(wèn)題，通過(guò)將節(jié)理單元的剛度矩陣?yán)奂拥较到y(tǒng)總體剛度矩陣來(lái)解決。以上這些方法在計(jì)算過(guò)程中，均需使用可視準(zhǔn)則。該準(zhǔn)則認(rèn)為，問(wèn)題邊界及任何內(nèi)部不連續(xù)面都是不透明的。也就是說(shuō)，當(dāng)節(jié)點(diǎn)影響域遇到這些邊界和不連續(xù)面時(shí)都將被截?cái)?，那些處于不連續(xù)線另一側(cè)的節(jié)點(diǎn)將會(huì)從影響域中刪除，導(dǎo)致系統(tǒng)方程的稀疏度大大降低，必須使用更大的影響域尺度來(lái)處理。

本文基于單位分解思想，在EFGM的位移模式中增加階躍函數(shù)項(xiàng)和裂尖奇異項(xiàng)，構(gòu)造無(wú)網(wǎng)格不連續(xù)位移場(chǎng)函數(shù)。然后參考文獻(xiàn)[6－7]的作法，將接觸摩擦條件轉(zhuǎn)化為包含懲罰因子的線性互補(bǔ)形式，建立求解摩擦接觸問(wèn)題的不連續(xù)無(wú)網(wǎng)格線性互補(bǔ)模型，利用已經(jīng)成熟的Lemke算法求解。該模型無(wú)需引入可視準(zhǔn)則且不需要迭代計(jì)算，大大提高了求解效率。

2 無(wú)網(wǎng)格不連續(xù)位移模式

標(biāo)準(zhǔn)EFGM的位移場(chǎng)近似函數(shù)可表示為

式中：pj為單項(xiàng)式基函數(shù)；m為pj的階次；aij為待定參數(shù)，由離散L2范數(shù)取極小值的條件確定，具體形式為

式中：φI(x)為EFGM的形函數(shù)，具有高階連續(xù)性、單位分解性等。在處理諸如裂紋等不連續(xù)問(wèn)題時(shí)，需考慮由不連續(xù)面或線引起的不連續(xù)性，因此，必須對(duì)式（3）的位移模式進(jìn)行修改。根據(jù)文獻(xiàn)[8－9]的結(jié)論：若形函數(shù)滿足單位分解特性，可根據(jù)特定問(wèn)題解的先驗(yàn)知識(shí)對(duì)近似函數(shù)進(jìn)行擴(kuò)展。因此，考慮線彈性斷裂力學(xué)Westergaard解，對(duì)EFGM近似函數(shù)進(jìn)行擴(kuò)展，得到不連續(xù)近似函數(shù)，具體形式為

式（4）中第1項(xiàng)是標(biāo)準(zhǔn)EFGM近似函數(shù)；第2項(xiàng)反映沿接觸面兩側(cè)的不連續(xù)；第3項(xiàng)則反映接觸面端部奇異性。Ndisc為階躍擴(kuò)展節(jié)點(diǎn)集合，由其影響域被接觸面切割的節(jié)點(diǎn)組成；NTip為裂尖擴(kuò)展節(jié)點(diǎn)集合，由裂尖的影響域內(nèi)包含的節(jié)點(diǎn)組成；αiJ、βiKl為附加自由度；f(x)為水平集函數(shù)，可以隱式表達(dá)不連續(xù)面，具體形式為

S(t)?R2為不連續(xù)界面，當(dāng)節(jié)點(diǎn)位于界面上方時(shí)取正號(hào)，下方時(shí)取負(fù)號(hào)。H(x)為Heaviside階躍函數(shù)，具體形式為

Tl為裂尖擴(kuò)展函數(shù)，表達(dá)式為

其中：r為節(jié)點(diǎn)x到裂尖xTip的距離，

由式（4）可知，計(jì)算點(diǎn)的影響域不受結(jié)構(gòu)內(nèi)部細(xì)節(jié)的影響，不連續(xù)界面引起的位移不連續(xù)性由附加自由度確定，無(wú)需使用可使準(zhǔn)則。

3 接觸面方程

為簡(jiǎn)單起見，考慮二維小變形情況下的面-面接觸。設(shè)ΩA，ΩB為兩個(gè)物體組成的接觸體系，潛在接觸邊界為Γc。接觸面上對(duì)應(yīng)點(diǎn)的位移定義為和法向和切向接觸應(yīng)力分別定義為pn，pτ。則接觸位移用相對(duì)位移描述

通常，接觸力和接觸位移應(yīng)當(dāng)滿足如下兩個(gè)定律：

（1）法向接觸定律

式（9）中3式的意義依次為非穿透條件、法向不受拉條件和互補(bǔ)條件。

（2）切向Coulomb 摩擦定律

式中：μ為摩擦系數(shù)；c為黏聚力。參照理想剛塑性體單向應(yīng)力-應(yīng)變關(guān)系的表述，式（9）、（10）可統(tǒng)一表述為

考慮接觸面條件，標(biāo)準(zhǔn)EFGM的虛功方程修改為

4 離散控制方程

式（17）含有接觸面上的積分，因此，必須對(duì)接觸面進(jìn)行離散。將接觸面劃分為有限個(gè)區(qū)段。每個(gè)區(qū)段兩側(cè)布設(shè)相同數(shù)目的高斯點(diǎn)，組成接觸點(diǎn)對(duì)。對(duì)于某一接觸點(diǎn)對(duì)，由式(4)可求得其相對(duì)位移為

式中：Bstd由式（3）的形函數(shù)導(dǎo)數(shù)組成；Benr由擴(kuò)展部分組成，當(dāng) xI∈NDsic時(shí)，ΨI(x)為階躍擴(kuò)展函數(shù) H(f (x ))-H(f (xJ))；當(dāng) xI∈NTip，ΨI(x)分別為裂尖擴(kuò)展函數(shù)為節(jié)點(diǎn)位移列向量集合，包括節(jié)點(diǎn)常規(guī)自由度和附加自由度；

將式（18）代入式（16）的互補(bǔ)條件，可得

5 數(shù)值算例

5.1 含節(jié)理平板

如圖1所示，寬為1 m、高為3 m的平板，中間含有一條貫通節(jié)理，底部固定，上部受均布?jí)毫Ζ液颓邢蛄Ζ拥淖饔谩Ｆ桨宀牧蠀?shù)彈性模量為E =2 kPa，泊松比為v = 0.3，節(jié)理面摩擦系數(shù)為μ = 0.3，凝聚力為c = 0。將整個(gè)計(jì)算區(qū)域劃分為21×31個(gè)背景積分網(wǎng)格，網(wǎng)格角點(diǎn)為節(jié)點(diǎn)，采用4×4高斯積分。將節(jié)理面劃分為20個(gè)區(qū)段，每個(gè)區(qū)段采用4個(gè)高斯積分點(diǎn)。首先計(jì)算σ = 50 N/m，τ = 0時(shí)的情況。當(dāng)懲罰因子En= Eτ的取值依次選取為E～1012E，計(jì)算的法向接觸應(yīng)力值與精確解（σ = 50 N/m）均吻合得很好，最大相對(duì)誤差僅為0.77%（En= Eτ= E時(shí)）；圖2給出了沿節(jié)理面上的法向接觸和切向接觸應(yīng)力分布（En= Eτ= 109E時(shí)）?？梢钥闯觯疚挠?jì)算的切向接觸應(yīng)力遠(yuǎn)小于法向接觸應(yīng)力，其分布反對(duì)稱于節(jié)理面。然后計(jì)算σ = 50 N/m，τ = 5的情況。圖3為沿節(jié)理面上的法向接觸和切向接觸應(yīng)力分布。可以看出，切向接觸應(yīng)力為拋物線型分布，法向接觸應(yīng)力為線性分布。以上兩種情況的計(jì)算結(jié)果與文獻(xiàn)[7]十分接近。

5.2 受壓裂紋應(yīng)力強(qiáng)度因子

如圖4所示，裂紋試樣尺寸為50 mm×50 mm，裂紋半長(zhǎng)為a = 5 mm，裂紋中心位于試樣中心。材料參數(shù)取為E = 5 GPa，泊松比v = 0.35，摩擦系數(shù)μ = 0，抗拉強(qiáng)度為σt= 3 MPa。試樣頂部和底部施加均勻的單軸壓縮荷載σ = 1 kPa。節(jié)點(diǎn)離散見圖4。在裂紋面光滑情況下，該問(wèn)題的應(yīng)力強(qiáng)度因子的理論解[10]為

圖1 含節(jié)理方板Fig.1 A rectangular plate with a joint

圖2 σ =50 N/m，τ = 0時(shí)節(jié)理面接觸應(yīng)力Fig.2 Contact stresses on the joint face for σ = 50 N/m, τ=0

圖3 σ =50 N/m，τ = 5 N/m時(shí)節(jié)理面接觸應(yīng)力Fig.3 Contact stresses on the joint face for σ =50 N/m, τ = 5 N/m

圖4 受壓裂紋試樣Fig.4 Compression-load crack specimen

圖5 應(yīng)力強(qiáng)度因子計(jì)算結(jié)果Fig.5 Computational results of stress intensity factors

5.3 彈性地基中的樁基抗拔問(wèn)題

本算例分析樁與樁側(cè)土的接觸?？紤]對(duì)稱性，選取一半?yún)^(qū)域計(jì)算，計(jì)算模型見圖6。土為黏性土，其參數(shù)為Es= 20 MPa，v = 0.25，ρ = 2 g/cm3?；炷翗兜牟牧蠀?shù)為 Ec= 2×104MPa，v = 0.3，ρ=2.5 g/cm3。樁土接觸面間滿足Mohr-Coulomb定律，凝聚力c = 0.05 Pa，摩擦系數(shù)μ = 0.58。底部邊界固定，兩個(gè)側(cè)邊施加水平位移約束，樁頂部施加P =5×102kN/m的拉力作用。

圖6 節(jié)點(diǎn)離散Fig.6 Node discreteness

圖7為本文方法和有限元中使用接觸單元[11]計(jì)算得到的接觸面切向應(yīng)力分布圖?？梢钥闯?，與文獻(xiàn)[11]相比，本文方法能更好地體現(xiàn)樁基與地基土之間的滑動(dòng)摩擦將自上而下地?cái)U(kuò)展的性質(zhì)。

圖7 接觸面上的剪應(yīng)力分布Fig.7 Shear stress distributions along contact surface

6 結(jié) 論

本文首先利用標(biāo)準(zhǔn)的無(wú)網(wǎng)格Galerkin法的單位分解性質(zhì)，通過(guò)在位移模式中嵌入不連續(xù)項(xiàng)表達(dá)由接觸面引起的位移不連續(xù)和接觸面端點(diǎn)的應(yīng)力奇異，構(gòu)造了不連續(xù)無(wú)網(wǎng)格位移場(chǎng)函數(shù)。然后結(jié)合摩擦接觸定律的線性互補(bǔ)描述，給出一種新型的求解摩擦接觸問(wèn)題的無(wú)網(wǎng)格線性互補(bǔ)方法。幾個(gè)數(shù)值算例驗(yàn)證了本文方法的可行性和有效性。巖石、混凝土等準(zhǔn)脆性材料的破壞形式大多屬于壓剪性裂紋開裂破壞。這類破壞形式涉及接觸非線性和裂紋斷裂擴(kuò)展兩方面結(jié)合的問(wèn)題，十分復(fù)雜，本文為無(wú)網(wǎng)格方法分析壓剪性裂紋擴(kuò)展問(wèn)題提供了有力的支持，具有良好的應(yīng)用前景。

[1]張雄, 劉巖, 馬上. 無(wú)網(wǎng)格法的理論及應(yīng)用[J]. 力學(xué)進(jìn)展, 2009, 39(1): 1－36.ZHANG Xiong, LIU Yan, MA Shang. Meshfree methods and their applications[J]. Advances in Mechanics, 2009,39(1): 1－36.

[2]BELYTSCHKO T, FLEMING M. Smoothing, enrichment and contact in the element-free Galerkin method[J].Computers and Structures, 1999, 71: 173－195.

[3]龐作會(huì), 葛修潤(rùn), 王水林. 無(wú)網(wǎng)格伽遼金法(EFGM)模擬巖體不連續(xù)面[J]. 工程地質(zhì)學(xué)報(bào), 2000, 8(3):364－368.PANG Zuo-hui, GE Xiu-run, WANG Shui-lin. Simulation discontinuity with element-free Galerkin method(EFGM)[J]. Journal of Engineering Geology, 2000, 8(3):364－368.

[4]李臥東, 陳勝宏. 接觸摩擦問(wèn)題的數(shù)值模擬[J]. 巖土力學(xué), 2003, 24(3): 385－388.LI Wo-dong, CHEN Shen-hong. Numerical modeling for frictional contact problems[J]. Rock and Soil Mechanics,2003, 24(3): 385－388.

[5]盧波, 丁秀麗, 鄔愛清. 無(wú)網(wǎng)格法對(duì)巖體不連續(xù)面的模擬[J]. 巖石力學(xué)與工程學(xué)報(bào), 2008, 27(10): 2108－2116.LU Bo, DING Xiu-li, WU Ai-qing. Modeling of rock discontinuity with meshless method[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(10):2108－2116.

[6]朱昌銘. 基于虛功原理的彈性接觸問(wèn)題的線性互補(bǔ)方法[J]. 力學(xué)學(xué)報(bào), 1995, 27(2): 189－197.ZHU Chang-ming. A linear complementarity method for elastic contact problems based on the principle of virtual work[J]. Acta Mechanica Sinica, 1995, 27(2): 189－197.

[7]余天堂. 摩擦接觸裂紋問(wèn)題的擴(kuò)展有限元法[J]. 工程力學(xué), 2010, 27(4): 84－89.YU Tian-tang. An extended finite element method for modeling crack problems with frictional contact[J].Engineering Mechanics, 2010, 27(4): 84－89.

[8]BABUSKA T, MELENK J M. The partition of unity method[J]. International Journal for Numerical Method in Engineering, 1997, 40: 727－758.

[9]VENTURA A, XU J X, BELYTSCHKO T. A vector level set method and new discontinuity approximations for crack growth by EFG[J]. International Journal for Numerical Methods in Engineering, 2002, 54: 923－944.

[10]李世愚, 和泰名, 尹祥礎(chǔ). 巖石斷裂力學(xué)導(dǎo)論[M]. 合肥: 中國(guó)科學(xué)技術(shù)大學(xué)出版社, 2010.

[11]LEI X. Contact friction analysis with a simple interface element[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190: 1955－1965.