模擬摩擦接觸問題的新型無網(wǎng)格數(shù)值方法

馬文濤 ，師俊平，李 寧

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

1 引 言

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

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

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

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

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

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

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

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

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

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

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

3 接觸面方程

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

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

（1）法向接觸定律

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

（2）切向Coulomb 摩擦定律

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

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

4 離散控制方程

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

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

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

5 數(shù)值算例

5.1 含節(jié)理平板

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

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

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

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

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

圖3 σ =50 N/m，τ = 5 N/m時節(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)力強度因子計算結(jié)果Fig.5 Computational results of stress intensity factors

5.3 彈性地基中的樁基抗拔問題

本算例分析樁與樁側(cè)土的接觸。考慮對稱性，選取一半?yún)^(qū)域計算，計算模型見圖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。底部邊界固定，兩個側(cè)邊施加水平位移約束，樁頂部施加P =5×102kN/m的拉力作用。

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

圖7為本文方法和有限元中使用接觸單元[11]計算得到的接觸面切向應(yīng)力分布圖。可以看出，與文獻[11]相比，本文方法能更好地體現(xiàn)樁基與地基土之間的滑動摩擦將自上而下地擴展的性質(zhì)。

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

6 結(jié) 論

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

[1]張雄, 劉巖, 馬上. 無網(wǎng)格法的理論及應(yīng)用[J]. 力學(xué)進展, 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]龐作會, 葛修潤, 王水林. 無網(wǎng)格伽遼金法(EFGM)模擬巖體不連續(xù)面[J]. 工程地質(zhì)學(xué)報, 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]李臥東, 陳勝宏. 接觸摩擦問題的數(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ǎng)格法對巖體不連續(xù)面的模擬[J]. 巖石力學(xué)與工程學(xué)報, 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]朱昌銘. 基于虛功原理的彈性接觸問題的線性互補方法[J]. 力學(xué)學(xué)報, 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]余天堂. 摩擦接觸裂紋問題的擴展有限元法[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]. 合肥: 中國科學(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.