模擬摩擦接觸問題的新型無網格數值方法

馬文濤 ，師俊平，李 寧

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

1 引 言

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

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

2 無網格不連續位移模式

標準EFGM的位移場近似函數可表示為

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

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

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

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

Tl為裂尖擴展函數，表達式為

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

由式（4）可知，計算點的影響域不受結構內部細節的影響，不連續界面引起的位移不連續性由附加自由度確定，無需使用可使準則。

3 接觸面方程

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

通常，接觸力和接觸位移應當滿足如下兩個定律：

（1）法向接觸定律

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

（2）切向Coulomb 摩擦定律

式中：μ為摩擦系數；c為黏聚力。參照理想剛塑性體單向應力-應變關系的表述，式（9）、（10）可統一表述為

考慮接觸面條件，標準EFGM的虛功方程修改為

4 離散控制方程

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

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

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

5 數值算例

5.1 含節理平板

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

5.2 受壓裂紋應力強度因子

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

圖1 含節理方板Fig.1 A rectangular plate with a joint

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

圖3 σ =50 N/m，τ = 5 N/m時節理面接觸應力Fig.3 Contact stresses on the joint face for σ =50 N/m, τ = 5 N/m

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

圖5 應力強度因子計算結果Fig.5 Computational results of stress intensity factors

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

本算例分析樁與樁側土的接觸。考慮對稱性，選取一半區域計算，計算模型見圖6。土為黏性土，其參數為Es= 20 MPa，v = 0.25，ρ = 2 g/cm3。混凝土樁的材料參數為 Ec= 2×104MPa，v = 0.3，ρ=2.5 g/cm3。樁土接觸面間滿足Mohr-Coulomb定律，凝聚力c = 0.05 Pa，摩擦系數μ = 0.58。底部邊界固定，兩個側邊施加水平位移約束，樁頂部施加P =5×102kN/m的拉力作用。

圖6 節點離散Fig.6 Node discreteness

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

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

6 結 論

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

[1]張雄, 劉巖, 馬上. 無網格法的理論及應用[J]. 力學進展, 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]龐作會, 葛修潤, 王水林. 無網格伽遼金法(EFGM)模擬巖體不連續面[J]. 工程地質學報, 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]李臥東, 陳勝宏. 接觸摩擦問題的數值模擬[J]. 巖土力學, 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]盧波, 丁秀麗, 鄔愛清. 無網格法對巖體不連續面的模擬[J]. 巖石力學與工程學報, 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]. 力學學報, 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]. 工程力學, 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]李世愚, 和泰名, 尹祥礎. 巖石斷裂力學導論[M]. 合肥: 中國科學技術大學出版社, 2010.

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