考虑冲击温度影响的物质点法算法

张静静 董龙雷

摘要：针对超高速冲击使结构瞬间产生极高的温度并造成结构软化现象的问题，提出考虑冲击温度影响的改进物质点法算法，采用Taylor杆的撞击过程模拟分析这一现象。将模拟结果与经典物质点法算法对同一算例的计算结果进行对比，结果显示改进物质点法算法计算结果精度可提高0.77%。

关键词：物质点法；塑性变形；冲击温度；超高速冲击

中图分类号：TB302.3；TB115.1

文献标志码：B

0 引 言

越来越多的学者研究有限元模拟大变形过程遇到的网格畸变问题。物质点法是1994年由SULSKY等[1-2]提出的，该算法采用拉格朗日质点和欧拉网格双重描述，是一种完全的拉格朗日质点类法。该方法将连续体离散为一组离散物质点，材料的物理信息，如速度、应力和位置等，存储在每一块材料区域对应的物质点上；欧拉网格仅用于动量方程的求解和空间导数的计算，不携带任何物质信息；在每一个计算时间步中，物质点和计算网格完全固连，不存在相对运动，避免欧拉法因非线性对流项产生的数值计算困难，且易于跟踪物质界面。另外，由于物质点携带物质的所有信息，因此每一步计算结束后，将已变形的背景网格丢弃，在新一步的计算过程中采用新的背景网格，因此物质点法可避免大变形问题中网格畸变产生的数值计算困难，在涉及大变形（如爆炸、冲击等）問题上具有明显优势。

物质点法最初采用显式积分法，适用于载荷作用时间短的瞬态问题，如爆炸、冲击等。很多学者根据各自需求对算法进行改进。目前，物质点法已广泛应用于各类工程模拟研究中。SULSKY等[1-2]模拟Taylor杆和钢球侵彻铝靶问题；HUANG等[3]和MA等[4]采用物质点接触算法模拟中低速冲击侵彻问题。在模拟材料失效方面，CHEN等[5-7]采用物质点法模拟冲击载荷下脆性材料的动态失效问题和材料在局部加热情况下的失效问题。

温度是影响材料特性的重要因素之一。在已有文献中，温度模拟过程大多只考虑大变形对结构造成温度升高的影响。在超高速冲击问题中，如Taylor杆撞击试验等，在冲击绝热压缩和冲击载荷作用下，材料的冲击波耗散效应会引起结构发生剧烈的温度变化，称为冲击温度。然而，经典的物质点法并没有将冲击温度考虑在内，引起模拟结果误差，因此本文在物质点法中考虑冲击温度的影响。

1 物质点法

物质点法将一个连续体离散为一系列物质点，见图1。连续体的密度可近似表达为

在网格的计算方式上，物质点法与有限元法非常相近，两者的区别[8]在于：

（1）有限元法采用高斯积分，将积分转化为被积函数在各个高斯点处的值与该高斯点所代表的体积之积的和。物质点法采用物质点积分，将积分转化为被积函数在各物质点处的值与该物质点所代表的体积之积的和。

（2）有限元法的计算网格始终与物体固连；物质点法的背景网格只在每个时间步内与物体固连，在每个时间步结束时，将已变形的背景网格丢弃，并在下一个时间步内采用新的规则的背景网格进行计算。由于物质点已经携带物体的所有物质信息，因此，在下一个计算时间步内，可将物质点信息映射到新的背景网格上得到网格信息。因此，物质点法在计算大变形的过程中不会出现网格畸变的现象。

物质点法与无网格方法（以光滑粒子流体动力学方法为例）相比，两者都需要对物质点进行离散并且将信息携带在各个物质点上，两者的区别在于：

（1）光滑粒子流体动力学方法中物质点与物质点的联系通过搜索物质点的邻域进行。物质点法在将物质点信息映射到背景网格时，已使各物质点之间产生联系，因此物质点法可避免物质点邻域搜索步骤，提高计算效率。

（2）光滑粒子流体动力学方法通过求解物质点组的动力学方程和跟踪每个物质点的运动轨道求得整个系统的力学行为，即其动量方程的求解也是在物质点组上实现的。物质点法将物质点的质量、速度等参数映射到背景网格上进行动量方程的求解。

开源的物质点法分析结构温度造成的影响因素仅考虑大变形情况，本文针对刨削产生和发展过程中产生高温现象的影响因素，如冲击、摩擦等，提出对标准物质点法的改进，并通过经典的Taylor杆撞击试验予以验证。

2 考虑冲击温度影响的物质点法算法

采用考虑材料的温度软化效应的Johnson-Cook本构模型进行验证，其对应公式为

式（4）中：第一项A+Bεnp为应变的表达式，A为材料名义屈服强度，B和n反映应力硬化的影响程度；第二项和第三项分别反映应变率和温度的影响程度。在应变和应变率均非常小且材料温度为室温的情况下，式（4）可简化为σy=A。当应力或应变率增大时，σy增加；当材料温度接近熔融温度时，温度影响项1-（T*）m接近于0，该现象导致材料屈服强度趋于0，结构极易发生变形。因此，当材料温度趋于熔融温度时，结构出现温度软化现象，结构屈服强度非常小，容易发生变形。以Taylor杆撞击试验[9]为例对该现象进行验证。

Taylor杆为铜杆，以190 m/s的速度向刚性墙撞击，铜杆初始长度L0=25.40 mm。在物质点法模拟过程中，物质点间距为0.38 mm，Taylor杆共离散为21 172个物质点。模拟时间80 μs，此时动能为0，基本不再变化，不考虑任何能量产生或损耗因素。Taylor杆的材料参数见表1。在该试验中，Taylor杆的最终长度L=16.20 mm，底端直径D=13.50 mm，距离底部0.2L0处直径W=10.10 mm。

2.1 温度软化效应验证

模拟材料温度为30、300、1 000和1 500 K的Taylor杆撞击，在80 μs时Taylor杆沿Y轴中心位置的切面见图2。在不同初始温度条件下，Taylor杆的长度L、底部直径D、距离底部0.2L0位置处的直径W的计算结果见表2。

圖2和表2表明，在相同的初始条件和计算参数下，随着材料温度的升高，结构的变形增大，说明随着温度的升高，材料硬度降低，结构越容易发生变形。刨削的产生和发展过程伴随着大量的热量产生，因此在对其过程进行数值分析时必须考虑温度，主要可以从塑性变形和冲击2个方面考虑物质点法中温度的影响。

2.2 塑性功

低应变率和高应变率变形过程的处理往往是不同的：低应变率现象常常处理为等温过程；高应变率现象常常处理为绝热过程。绝热过程中结构变形所做的功常转换为热能进而导致结构温度升高。因此，当材料发生高应变率变形时，大部分塑性功会转换为热能。

以Taylor杆试验为例，分析材料初始温度为30 K、考虑塑性变形做功的影响下物质点法算法的计算结果。Taylor杆沿Y轴中心位置的切面在80 μs时的结果见图3，最终杆长L=16.09 mm，底部直径D=12.89 mm，距离底部0.2L0位置处的直径W=9.28 mm。与试验结果相比，计算结果误差为4.44%。

2.3 冲击温度

冲击温度是材料在冲击压缩下的重要物理参数。[10]从原理上讲，冲击温度可通过材料的热力学函数、守恒方程和冲击绝热线计算获得。

根据文献统计，冲击温度的计算大致可分为3种近似方法：（1）利用Mie-Grüneisen状态方程和Hugoniot曲线计算；（2）利用三项式状态方程计算；（3）利用以等熵线为参考线的Mie-Grüneisen状态方程计算。[11-16]汤文辉等[17]对冲击温度近似算法的研究结果表明：方法（1）涉及的参数少，没有考虑电子热运动和固液相变的影响；方法（2）考虑的因素最全面，既考虑电子热运动的贡献，又考虑固液相变的影响，因此涉及的参数较多；方法（3）与等熵线有关，不考虑电子热运动和固液相变的影响，其可靠性取决于等熵方程的选取。从计算结果与试验结果对比分析看，方法（2）在很宽的压强范围内与试验结果最符合，而方法（1）仅仅在固相区与试验结果符合。这是由于方法（1）未考虑固液相变的影响，因此不适用于材料发生熔融后的状态区域。本文在研究过程中不考虑熔融现象，数值仿真过程仅考虑固相结构，由于方法（1）涉及的参数少，因此将Mie-Grüneisen状态方程和Hugoniot曲线计算冲击温度的方法添加到物质点法算法中。

由热力学关系可知

同样以Taylor杆试验为例，分析材料初始温度为30 K、考虑冲击温度下物质点法算法的计算结果。

Taylor杆沿y轴中心位置的切面在80 μs时的结果见图4，最终杆长L=15.88 mm，底部直径D=13.56 mm，距离底部0.2L0位置处的直径W=9.25 mm。与试验结果相比，计算结果误差为3.61%。

综合考虑塑性功和冲击温度对结构计算结果的影响，Taylor杆沿y轴中心位置的切面在80 μs时的计算结果见图5，最终杆长L=16.00 mm，底部直径D=13.40 mm，距离底部0.2L0位置处的直径W=9.26 mm。与试验结果相比，计算结果误差为3.31%。

综合以上结果，采用物质点法对高速冲击过程进行模拟须考虑冲击温度对结构温度分布的影响。

3 结 论

针对物质点法算法中温度的影响展开研究，在物质点法中同时考虑塑性功、冲击温度和摩擦对温度的影响，对算法进行改进。以Taylor杆撞击试验为例，温度对结构计算影响很大，随着温度的升高，金属出现温度软化效应，且冲击温度造成的温升影响远远高于大变形造成的温升影响，验证典型物质点法算法中考虑冲击温度等因素的必要性。这一改进可为刨削机理的研究打下基础。

参考文献：

[1] SULSKY D， CHEN Z， SCHREYER H L. A particle method for history-dependent materials[J]. Computer Methods in Applied Mechanics and Engineering， 1994， 118（1-2）： 179-196. DOI： 10.1016/0045-7825（94）90112-0.

[2] SULSKY D， ZHOU S J， SCHREYER H L. Application of a particle-in-cell method to solid mechanics[J]. Computer Physics Communications， 1995， 87（1-2）： 236-252. DOI： 10.1016/0010-4655（94）00170-7.

[3] HUANG P， ZHANG X， MA S， et al. Contact algorithms for the material point method in impact and penetration simulation[J]. International Journal for Numerical Methods in Engineering， 2011， 85（4）： 498-517. DOI： 10.1002/nme.2981.

[4] MA Z， ZHANG X， HUANG P. An object-oriented MPM framework for simulation of large deformation and contact of numerous grains[J]. Computer Modeling in Engineering and Sciences， 2010， 55（1）： 61-87. DOI： 10.3970/cmes.2010.055.061.

[5] CHEN Z， HU W， SHEN L， et al. An evaluation of the MPM for simulating dynamic failure with damage diffusion[J]. Engineering Fracture Mechanics， 2002， 69（17）： 1873-1890. DOI： 10.1016/S0013-7944（02）00066-8.

[6] CHEN Z， FENG R， XIN X， et al. A computational model for impact failure with shear-induced dilatancy[J]. International Journal for Numerical Methods in Engineering， 2003， 56（14）： 1979-1997. DOI： 10.1002/nme.651.

[7] CHEN Z， GAN Y， CHEN J. A coupled thermo-mechanical model for simulating the material failure evolution due to localized heating[J]. Computer Modeling in Engineering and Sciences， 2008， 26（2）： 123-137. DOI：10.3970/cmes.2008.026.123.

[8] 張雄， 廉艳平， 刘岩， 等. 物质点法[M]. 北京： 清华大学出版社， 2013： 41-42.

[9] JOHNSON G R， HOLMQUIST T J. Evaluation of cylinder： Impact test data for constitutive model constants[J]. Journal of Applied Physics， 1988， 64（8）： 3901-3910. DOI： 10.1063/1.341344.

[10] WILLIAMS Q， JEANLOZ R， BASS J， et al. Melting curve of iron to 250 Gigapascals： A constraint on temperature at Earths center[J]. Science， 1987， 236（4798）： 181-183. DOI： 10.1126/science.236.4798.181.

[11] BROWN J M， MCQUEEN R G. Phase transitions， Grüneisen parameter， and elasticity for shocked iron between 77 GPa and 400 GPa[J]. Journal of Geophysical Research： Solid Earth， 1986， 91（B7）： 7485-7494. DOI： 10.1029/JB091iB07p07485.

[12] RICE M， MCQUEEN R G， WALSH J. Compression of solids by strong shock waves[J]. Solid State Physics， 1958， 6： 1-63. DOI： 10.1016/S0081-1947（08）60724-9.

[13] MCQUEEN RG， MARSH S P. The equation of state of solids from shock wave studies[C]∥Proceedings of High Velocity Impact Phenomena. New York： Elsevier， 1970： 293.

[14] LYZENGA G， AHRENS T J. Multiwavelength optical pyrometer for shock compression experiments[J]. Review of Scientific Instruments， 1979， 50（11）： 1421-1424. DOI： 10.1063/1.1135731.

[15] MILLER G H， AHRENS T J， STOLPER E M. The equation of state of molybdenum at 1 400 ℃[J]. Journal of Applied Physics， 1988， 63（9）： 4469-4475. DOI： 10.1063/1.341124.

[16] AHRENS T J. Shock-induced optical radiation from solids[C]// Proceedings of 2nd International Symposium on Intense Dynamic Loading and its Effects. Chengdu： Sichuan University Press， 1994.

[17] 汤文辉， 张若棋， 胡金彪， 等. 冲击温度的近似计算方法[J]. 力学进展， 1998， 28（4）： 479-487.

（编辑 武晓英）