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

盧嘉錚 姚星宇

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

【關(guān)鍵詞】物質(zhì)點法；有限元法；赫茲接觸

中圖分類號： 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 前言

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

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

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

本文將采用最優(yōu)運輸無網(wǎng)格法計算經(jīng)典赫茲接觸問題，對比解析解和有限元解與OTM解。

2 最優(yōu)運輸無網(wǎng)格法

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

1）LME具有弱克羅內(nèi)科德爾塔性質(zhì)，第0階和第1階連續(xù)性。在頂點處，形函數(shù)滿足克羅內(nèi)科德爾塔性質(zhì)，且內(nèi)部點與邊界無關(guān)。

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

3）由于局部度是可調(diào)的，所以在構(gòu)建LME的各向異性形函數(shù)和高階近似時可根據(jù)不同情況靈活應(yīng)變，從而使問題變得簡單。

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

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

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

2 赫茲接觸問題研究

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

2.1 OTM解

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

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

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

2.2 有限元解

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

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

2.3 赫茲理論解

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

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

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

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

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

3 結(jié)論

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

【參考文獻】

[1]廉艷平，張帆，劉巖，張雄.物質(zhì)點法的理論和應(yīng)用[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.