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

盧嘉錚 姚星宇

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

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

中圖分類號(hào)： O241.82；O35 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 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 前言

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

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

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

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

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

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

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

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

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

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

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

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

2 赫茲接觸問(wèn)題研究

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

2.1 OTM解

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

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

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

2.2 有限元解

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

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

2.3 赫茲理論解

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

其中，l為圓柱體長(zhǎng)度，在本例的解析方法計(jì)算中取單位長(zhǎng)度1。

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

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

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

3 結(jié)論

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

【參考文獻(xiàn)】

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