物质点方法与有限元方法在二维接触问题中的对比

卢嘉铮 姚星宇

【摘 要】针对非线性力学问题，特别是接触类问题时，有限元法经常暴露出时间成本高，计算结果不容易收敛等先天性缺陷。相比之下，物质点法具有更好的性质，接触应力由物质点之间的动量守恒直接计算得到，理论上计算开销更小，同时具备更高的计算精度和准确度。本文采用物质点方法、有限元法和解析方法对赫兹接触问题进行求解，并将三种方法计算得到的最大接触应力进行比较。计算分析后发现，相比有限元法，物质点法的计算结果的精确度与准确度略高。本文的分析研究为计算接触问题提供了新的思路，具备一定工程应用价值。

【关键词】物质点法；有限元法；赫兹接触

中图分类号： 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.