物質點方法與有限元方法在二維接觸問題中的對比

盧嘉錚 姚星宇

【摘 要】針對非線性力學問題，特別是接觸類問題時，有限元法經常暴露出時間成本高，計算結果不容易收斂等先天性缺陷。相比之下，物質點法具有更好的性質，接觸應力由物質點之間的動量守恒直接計算得到，理論上計算開銷更小，同時具備更高的計算精度和準確度。本文采用物質點方法、有限元法和解析方法對赫茲接觸問題進行求解，并將三種方法計算得到的最大接觸應力進行比較。計算分析后發現，相比有限元法，物質點法的計算結果的精確度與準確度略高。本文的分析研究為計算接觸問題提供了新的思路，具備一定工程應用價值。

【關鍵詞】物質點法；有限元法；赫茲接觸

中圖分類號： O241.82；O35 文獻標識碼： A 文章編號： 2095-2457（2018）05-0012-003

【Abstract】For nonlinear mechanical problems，especially the contact problems，the finite element method often exposes the high time cost，and the calculation results are not easy to converge and other congenital defects.In contrast，the material point method has better properties.The contact stress is directly calculated from the momentum conservation between material points，which has less computational cost and higher accuracy.In this paper，the material point method， finite element method and analytical method are used to solve the Hertz contact problem，and the maximum contact stress calculated by the three methods is compared.The calculation results show that compared with the finite element method，the accuracy of the material point method is slightly higher.The analysis and study of this paper provides a new way of thinking for the calculation of contact problems，and has a certain value of engineering application.

【Key words】Material point method；Finite element method；Hertz contact problem

1 前言

接觸問題具有強非線性，提高計算結果的精度和準確度一直以來是學者們研究的重點。20世紀，計算機技術的發展使得有限元等數值方法得以展開拳腳，迅速被應用于求解接觸問題。但是，有限元法中的接觸計算主要采用罰函數法，兩物體之間的接觸狀態未知，在計算的每一個增量步前后，都需要對接觸面進行搜尋，并且約束條件不能被嚴格滿足，因此有限元接觸計算經常出現貫穿、不收斂等問題。

物質點法（Material Point Method， MPM）是由Sulsky和Chen于1994年提出的一種數值方法[1]，其本質是一種采用質點和網格雙重描述的無網格法。物質點法采用質點離散材料區域，通過背景網格計算空間導數和求解動量方程，避免了網格畸變和對流項處理，兼具歐拉和拉格朗日描述的優點，非常適合用于模擬涉及大變形、沖擊和斷裂破碎等問題。但是，MPM中的近似方法不具有克羅內科德爾塔性質，不能解決邊界條件施加的問題[1]。

Arroyo和Ortiz提出了局部最大熵近似[2]（Local Maximum-Entropy Approximation Schemes，LME），該方法擁有“弱”克羅內科德爾塔性質，即在邊界上具有克羅內科德爾塔性質，使得在物質點法中，施加邊界條件變得簡單。Bo Li[3]基于LME提出了最優運輸無網格法（Optimal Transportation Meshfree，OTM），OTM解決了其他無網格法中強制邊界條件施加與伽遼金弱形式數值積分的問題，為強非線性問題的求解提供了一個全新的思路。而OTM本質上也是物質點法的一種。

本文將采用最優運輸無網格法計算經典赫茲接觸問題，對比解析解和有限元解與OTM解。

2 最優運輸無網格法

本文采用的物質點法為最優運輸無網格法，其局部最大熵插值函數具有非常好的性質。局部最大熵插值函數是通過對最大熵插值函數進行寬度限制推導得出的，是一種凸擬合。局部最大熵近似法（LME）具有很多凸近似法的理想性質，并有以下明顯優勢：

1）LME具有弱克羅內科德爾塔性質，第0階和第1階連續性。在頂點處，形函數滿足克羅內科德爾塔性質，且內部點與邊界無關。

2）形函數的局部性實現了有限元形函數和無網格近似法之間的無縫對接。參數決定了形函數的支持寬度。另外，可根據節點位置的不同來調整適應不同的局部度，這使LME適用于流固耦合問題或是極大變形等問題。

3）由于局部度是可調的，所以在構建LME的各向異性形函數和高階近似時可根據不同情況靈活應變，從而使問題變得簡單。

4）形函數的計算非常高效。一方面，LME形函數的計算不需要專門求出N個未知數，而僅是一個無約束最小問題計算結果的附屬產物。另一方面，形函數將按衰減，因此僅有很少的節點對配分函數有貢獻，極大地減少了計算開銷。

基于LME，Bo Li提出了最優運輸無網格法。該方法原理如下：首先對無相互作用流采用最優運輸定理求解，即找到一個最優的運輸映射，使得將初始質量密度運輸到終點質量密度的運輸成本最小。運輸成本可由初始質量密度和終點質量密度之間的Wasserstein距離表示，引入離散的拉格朗日動力學理論，通過對時間和空間的離散作用，得到無相互作用流的離散歐拉-拉格朗日方程。基于無相互作用流的求解方式，考慮采用最優運輸定理求解固體流問題。基于固體流運動方程的變分，通過對時間和空間的離散，得到各項同性彈性固體流的離散歐拉-拉格朗日方程，并得出物質點更新的時間顯式求解步驟：

針對接觸問題，OTM與有限元法有著本質上的差異。OTM方法的接觸計算是基于動量守恒，即兩個控制方程，通過物質點間的動量守恒計算直接得到相應物質點處的接觸應力大小，理論上OTM方法的計算開銷較小，同時具備較高精度。另外，由于局部最大熵形函數的優良性質，可直接采用通用有限元軟件對模型進行前處理。

2 赫茲接觸問題研究

兩圓柱或兩圓球之間的接觸是典型的Hertz接觸問題[4]，如圖1所示，兩圓柱體的軸垂直于xy平面，在單位長度上的力P作用下發生接觸。相對接觸區域，兩圓柱的尺寸足夠大，假設接觸面無摩擦，并將兩圓柱看做彈性半空間體，圓柱體材料為無硬化的理想各向同性彈性體，則將該二圓柱接觸問題轉變為二維接觸問題。已知R1=15mm，R2=10mm，彈性模量E=20GPa，泊松比v=0.3。

2.1 OTM解

兩圓柱接觸是平面應變問題，考慮其對稱性，可各取其二分之一作為數值計算的幾何模型。本例的靜態接觸問題，施加位移邊界條件，制定加載參數。在OTM中，加載方向、加載速度和總時間等參數通過submit.sh文件進行設置，邊界條件則需要在OTM算例主程序中，通過C++語句實現。在本例中，加載點為小半圓柱的頂部節點組top_nodes，載荷設置為延Y軸負方向的位移S=1mm。另外，還需在算例的主程序中分別對已分組的節點bottom_nodes和central_nodes設置約束，對bottom_nodes組施加全約束，對central_nodes組施加X方向的位移約束。

主程序編譯通過之后，即可開始計算。本例中，計算共進行了153步，計算結果輸出生成153個.vtu文件，該文件包含了在某一時間步中，所有的計算結果，如平均應力、有效應力、主應變和主應力等。采用ParaView進行計算結果的后處理。t=0ms、t=0.33ms、t=0.67ms和t= 1ms時刻的主應變如下圖2所示。

從上圖可看出，主應變主要發生在半圓的直徑附近，而圓弧附近變形較小，這是由于外力延直徑方向加載，且接觸點也為該直徑的一個端點，所以主要為直徑附近的材料受壓縮變形。下圖3為t=1ms時，接觸區域有效應力分布示意圖。查看最大接觸應力，為：Pmax=5790.9MPa。

2.2 有限元解

采用ANSYS進行計算。由于幾何模型較為簡單，因此選用一階平面單元PLANE182對材料進行離散。本例中選擇點—面接觸方式，大圓柱接觸區域添加目標單元TARGE169，小圓柱接觸區域添加接觸單元CONTA172。有限元模型如下圖4所示：

圖4所示的模型中，固定兩半圓的直徑在X方向的位移，固定大半圓底部節點，在小半圓的頂點施加位移約束，位移延Y軸負方向，位移大小為1mm。邊界條件施加完畢，選用增廣拉格朗日法求解。完成計算后，接觸應力如下圖5所示：提取接觸區域內單元的最大接觸應力，為：Pmax=5929.8MPa。

2.3 赫茲理論解

根據通用赫茲理論，兩圓柱接觸問題的接觸半寬為：

其中，l為圓柱體長度，在本例的解析方法計算中取單位長度1。

聯立（3.1）和（3.2），代入已知的彈性模量，兩圓柱直徑d以及載荷F。計算出本例中的最大接觸應力Pmax=5694.36MPa。

綜上，針對兩圓柱赫茲接觸問題的OTM解、有限元解和解析解如下表2所示，并列出了OTM解和有限元解相對解析解的誤差。

從上表可看出，對于完全相同的網格模型，有限元解的相對誤差是OTM解相對誤差的2.4倍，OTM解的精度比有限元更高，更加逼近于理論解。該結果也證明了，由于計算方法本質上的不同，對于接觸問題，OTM方法比有限元法有著更精準的求解結果。

3 結論

本文針對經典赫茲圓柱接觸問題，將模型進行簡化為平面應變問題。用物質點法（OTM法）有限元法和解析方法分別計算同一接觸問題。對三種不同方法的計算結果進行分析對比，發現以解析解為標準解，有限元解的相對誤差是OTM解的相對誤差的2.4倍，證明了OTM方法在計算接觸問題時相對于有限元法有著較高的精確度與準確度，這是由于OTM方法和有限元法在本質上的差異所決定的。本文的研究工作對工程實際中各類接觸問題的解決方式有一定參考價值，工程技術人員可在處理某些棘手的接觸問題時，考慮采用以OTM等為代表的物質點方法。

【參考文獻】

[1]廉艷平，張帆，劉巖，張雄.物質點法的理論和應用[J].2013，43（2）：237～264.

[2]M.Arroyo and M.Ortiz.Local maximum-entropy approximation schemes：A seamless bridge between finite elements and meshfree methods[J].International Journal for Numerical Methods of Engineering.2006，65：2167～2202.

[3]Bo Li.The optimal transportation method in solid mechanics[D].Pasadena，California：California Institute of Technology，Degree of Doctor of Philosophy，2009.

[4]Wang X C，Chang L M，Cen Z Z.Effective Numerical Methods for Elasto-Plastic Contact Problems with Friction[J]. Acta Mechanica Sinica，1990，6：349～356.