一种基于颗粒接触的二维无网格方法 及其在高速冲击模拟中的应用*

冯 春,李世海,刘晓宇

(中国科学院力学研究所，北京 100190)

基于连续介质力学的离散元方法(continumm-based discrete element method, CDEM)通过有限元与离散元的有机结合，在分析材料的变形断裂及断裂后的运动方面具有一定的优势[1-2],但该方法在模拟高速冲击等大变形问题时，会出现网格畸变导致的数值系统能量发散等现象。为解决上述问题，本文中提出了一种基于颗粒接触的二维无网格方法(particle contact-based meshfree method, PCMM)。

由于网格畸变导致的数值系统发散是所有拉格朗日架构下，基于网格的数值计算方法所面临的共同问题。目前的解决方法主要包括生死单元法及自适应网格等两大类。生死单元法是将严重畸变的单元进行钝化处理，从而保证了每一计算时步的网格质量，该方法虽然简单易行，但单元挖空后导致的应力场不连续问题值得探讨。自适应网格[3-4]根据数值系统内局部区域的网格畸变程度，对该区域网格进行调整或加密，从而保证了网格质量，但网格重划分的耗时及新生单元各类场信息的插值问题仍然值得关注。

无网格方法是近年来的研究热点，该方法的核心是通过空间域内一系列离散节点构建连续介质插值函数，由于摆脱了网格的限制，因此可以解决因网格畸变导致的系统能量发散等问题。目前发展的无网格方法有几十种，这些方法的主要区别在于使用不同的插值函数对微分方程进行等效，如光滑粒子法(SPH)、物质点法(MPM)等[5]。

SPH最早用于天体物理方面的研究[6-7]，核心思想是将空间场量的偏导数转化为对设定核函数的偏导数，从而简化了计算。近年来，SPH及其扩展形式在自由表面流动、高速冲击、材料断裂等方面得到了广泛应用[8-10]。尽管如此，该方法在边界处的处理方式、张力不稳定的改进方法及支持域内粒子的快速搜索[11-13]等方面仍然需要进行深入的研究。

MPM的主要思想是空间离散质点系在背景网格中的流动计算。该方法将物质点信息映射到背景网格并进行动量方程求解，获取网格点的速度、位移等信息，而后将网格点运动量映射回物质点，更新物质点的空间坐标[14-15]。该方法发挥了拉格朗日构型(物质点)及欧拉构型(背景网格)的双重优势，但在背景网格尺寸的选取原则及质点穿过背景网格产生的数值振荡等方面仍然需要仔细研究。

本文中在CDEM方法的基础上，结合颗粒离散元中接触对的概念，提出了基于颗粒接触的二维无网格方法(PCMM)，并基于VC++编制了相应的计算程序。

1 PCMM方法的基本原理

1.1 核心思想

采用显式解法计算高速冲击问题，当单元出现畸变后会严重影响计算时步，进而导致系统能量的发散。颗粒离散元中颗粒间的接触信息及拓扑关系复杂而有序，基于颗粒间的接触关系可构建一系列的连续介质单元，通过在单元中引入宏观本构关系即可表述连续介质在冲击载荷作用下的动态响应。当上述连续介质单元的变形逐渐增大，颗粒间的接触关系也将逐渐发生改变；当单元出现严重畸变(如压扁、拉长、扭曲等)时，颗粒间的原有接触关系被打破，基于原有接触关系构建的连续介质单元(畸变单元)自动删除；同时，基于当前颗粒接触关系的单元将自动创建。通过旧单元的删除及新单元的创建，即可解决高速冲击计算中因单元畸变导致的数值系统发散等问题。基于以上想法，文本中提出了PCMM方法，该方法的基本流程如图1所示。

图1 PCMM方法的基本流程Fig.1 Basic steps of PCMM

1.2 三角形单元的创建

三角形单元创建的必备条件包括：

(1)构成三角形单元的3个颗粒必须彼此接触：

dij≤Ri+Rj+δi,j=1,2,3;i≠j

(1)

(2)所构建三角形单元的任何一个内角角度应在30°～150°之间(保证单元不会太奇异)：

30°≤θi≤150°i=1,2,3

(2)

(3)构成的三角形单元3条边的任意一条边不应该小于3个颗粒平均半径的0.5倍：

(3)

式中：dij为颗粒i与j之间的距离，Ri、Rj为颗粒i与j的半径，δ为颗粒i与j之间的接触容差，θi为三角形单元的某个内角。在随机排列的颗粒体系中，为了构建稳定的连续介质单元系统，接触容差δ必须足够大以防止空隙的发生。

为提升颗粒接触对的检索效率，本文中采用了二维空间分割法(staticcell)，计算复杂度为O(N)。该方法将二维空间域按照颗粒的最大半径进行正交格子划分，并将每个颗粒按照其质心坐标映射到对应的格子内(见图2)；接触寻找时循环所有颗粒，根据当前颗粒的质心坐标确定颗粒所在格子的位置(图2中颗粒A的位置为x=3,y=3),在颗粒A周围的9个格子内(包含自身所在格子)寻找潜在接触颗粒。

如潜在颗粒i到颗粒A的距离小于设定容差(见图3)，将颗粒i增加至颗粒A的接触链表中。完成颗粒A的邻居检索后，计算邻居颗粒到颗粒A的方位角，并将邻居颗粒按照方位角进行排序。排序过程中，如果相邻2个邻居颗粒之间的方位角之差小于30°，将距离颗粒A较远的那个接触颗粒从颗粒A的接触链表中删除(防止产生低质量单元)。完成了颗粒A邻居颗粒的排序及局部删除后，重新循环颗粒A的接触链表，判断颗粒A，第i个及i+ 1个邻居颗粒组成的三角形是否满足式(1)～(3)所示条件，如果满足条件且该三角形单元不在颗粒A的单元链表中(新单元)，创建三角形单元并增加至颗粒A的单元链表中。

图2 用于邻居搜索的二维空间盒Fig.2 2D bin array for neighbour searching

图3 基于方位角的邻居排序Fig.3 Neighbour sequence according to position angle

1.3 三角形单元的删除

经过一段时间的迭代计算后，单元将出现大变形，如果不进行任何修正继续计算，将出现由于单元畸变导致的系统发散等问题。考虑到所有场量(如密度、速度、温度、应力、应变等)均位于颗粒上，且单元畸变后某一个方向尺寸非常小(单元体积可忽略)，将畸变后的单元删除不会引起系统总质量及总能量的改变。因此，提出如下的三角形单元删除原则(满足以下3个条件的任何一个，该三角形单元即被删除)：

(1)组成该三角形单元的3个颗粒不再彼此接触；

(2)该三角形单元的任何一个内角小于30°或者大于150°；

(3)该三角形单元的任何一条边长小于3个颗粒平均半径的0.5倍。

1.4 单元内力及节点运动的计算

三角形单元创建完毕后，即可采用经典连续介质理论进行单元变形力的求解。为了模拟高速冲击问题，引入流体弹塑性模型，并采用增量的计算方式，基本计算流程为：

(4)

(6)根据物态方程计算压力(体应力)p，本文中采用Mie-Grüeisen物态方程。

(8)计算数值阻尼应力：采用人工q阻尼消除强间断带来的计算稳定性问题

(5)

(10)根据弹性力学中斜截面公式，将应力转换为节点力

(6)

式中：N表示第N个节点，k为节点N所在棱的序号(三角形中每个节点隶属于2个棱)。

(11)计算颗粒所受合外力Fi，根据牛顿第二运动定律计算颗粒的加速度、速度、位移：

(7)

(12)回到步骤(1)进行下一时步速度梯度的求解。

2 流体弹塑性模型的引入

流体弹塑性模型是描述材料在爆炸冲击作用下动态响应的经典模型，该模型中体应力由物态方程控制，偏应力由胡克定律控制并受塑性准则修正。本文中采用Mie-Grüeisen物态方程表征材料的流体特征，采用胡克定律表征材料的弹性特征，采用Johnson-Cook模型表征材料的塑性特征。

2.1 Mie-Grüeisen物态方程的引入

以冲击Hugoniot态为参考态的Grüeisen物态方程可表述为：

(8)

式中：μ=ρ/ρ0-1，ρ0为初始密度，ρ为当前密度(ρ=(V0/V)ρ0)，c0、λ为高压下冲击波速度与粒子速度拟合函数的系数(D=c0+λu，D冲击波速度，u粒子速度)，γ0、a为量纲一材料参数。

上述状态方程中还有一个变量E，它是系统内能，可用下式求得，

(9)

2.2 Johnson-Cook模型的引入

Johnson-Cook模型(JC模型)是Johnson和Cook于1983年提出的用于高应变率和高温下材料大变形的本构模型：

(10)

当达到上述屈服应力后，需对6个偏应力分量进行修正。由于JC模型只定义了屈服条件，没有定义流动法则，因此实际计算时，当应力超过上述屈服强度，让其等于此屈服强度。

对于当前温度T的处理有2种方式：一种方式认为数值计算过程中当前温度T保持不变，这种方式在参数输入时设定初值即可。还有一种认为当前温度T在数值计算过程中会因为塑性功而改变。

对于第2种方式，需根据塑性功与温度的等效关系进行当前温度T的求解，每时步单位体积内的塑性功增量可由下式得出，

(11)

每时步的温度增量可由下式得出，

(12)

式中：c为材料的比热，α为功热转化率。当前时步的温度值可表示为T=Troom+∑ΔT。

3 高速冲击算例

3.1 弹性杆撞击

根据文献[9]，建立直径1 mm、高20 mm的杆件，并用2 211个规则排布的颗粒进行离散，颗粒半径为0.05 mm(见图4)。弹性杆密度7 850 kg/m3，弹性模量206 GPa，泊松比0.3。弹性杆以50 m/s的速度向基板撞击，观察杆件正中间测点轴向速度随时间的变化(见图5)。由图5可得，数值解与理论解基本一致，由此表明了PCMM方法在求解弹性动力问题方面的准确性。

图4 弹性杆撞击数值模型Fig.4 Numerical model of elastic bar impacting rigid wall

图5 杆件中部测点轴向速度曲线Fig.5 Axial velocity curve on the middle point of bar

3.2 泰勒杆

根据文献[16]，建立直径10 mm、高70 mm的OFHC铜杆，撞击速度165 m/s。杆件采用17 901个规则排布的颗粒进行离散，颗粒半径为0.1 mm，并采用流体弹塑性模型进行撞击分析。取A=89.63 MPa，B=45MPa，n=0.31，C=0.025，m=1.09；Mie-Grüeisen物态方程中取ρ=8 940 kg/m3，c0=3 940m/s，λ=1.49，γ=2.02，a=1.5；胡克定律中剪切模量G取60.74 GPa。不同时刻杆件的变形如图6所示，最终状态下杆件各位置的半径如图7所示。由图6～7可得，PCMM计算获得的杆件运动过程及最终形态与文献[16]的实验结果基本吻合，表明了PCMM方法在模拟高速冲击破坏效应时的准确性及合理性。

图6 不同时刻杆件的变形Fig.6 Deformation of Taylor bar at different times

图7 最终状态下杆件各位置的半径Fig.7 Bar radius at different position in final state

3.3 碎片云

碎片云数值计算模型中，弹丸直径2 mm，冲击速度5 km/s，靶板厚1 mm，靶板高16 mm。弹丸及靶板采用8 491个随机排布的颗粒进行离散，颗粒半径0.015～0.040 mm。弹丸及靶板的材料均为铝，Johnson-Cook模型中，取A=324 MPa，B=114 MPa，n=0.42，C=0.016，m=1.34，初始温度298 K，融化温度877 K，比热容875 J/(K·kg)，功热转化率0.9；Mie-Grüeisen物态方程中取ρ=2 703 kg/m3，c0=5 350m/s，λ=1.34，γ=1.97，a=1.5；胡克定律中剪切模量G取26.7 GPa。弹丸撞击靶板后2.68 μs时刻碎片云的空间形态如图8所示，该结果与文献[17-19]给出的数值及实验结果基本一致，表明了PCMM方法在模拟超高速碰撞问题中的准确性及合理性。

图8 2.68 μs 时碎片云的形态Fig.8 Shape of debris clouds at 2.68 μs

图9 不同时刻子弹、靶板系统的水平速度云图Fig.9 Horizontal velocity contour of bullet target system at different times

3.4 子弹入射靶板

子弹入射靶板模型中，子弹材料为钨，入射速度1 140 m/s，长2.84 cm，直径8 mm；靶板材料为铝，板高20 cm，宽5 cm。对子弹及靶板采用11 428个随机排布的颗粒进行离散，颗粒半径0.3～0.7 mm。靶板及子弹均采用流体弹塑性模型，靶板的材料参数与3.3节碎片云算例的参数一致。子弹的材料参数如下：Johnson-Cook模型中，取A=1.2 GPa，B=1.03 GPa，n=0.019，C=0.034，m=0.4，初始温度298 K，融化温度1 723 K，比热容134 J/(K·kg)，功热转化率0.9；Mie-Grüeisen物态方程中取ρ=19.22 g/cm3，c0=4.02km/s，λ=1.24，γ=1.67，a=1.3；胡克定律中剪切模量G取134.8GPa。子弹撞击靶板后10、32及65μs时,子弹、靶板系统水平速度云图(按照颗粒组成的连续介质单元显示)如图9所示。图9给出的子弹对靶板的侵彻作用与实际较为接近，证明了PCMM方法的合理性。

4 结 论

PCMM方法通过颗粒的接触关系实现了连续介质单元的快速创建及删除，解决了传统CDEM方法中因网格畸变导致的系统能量发散等问题。在PCMM中引入流体弹塑性模型后，可精确计算高速冲击下材料的动态响应及变形破坏特征。弹性杆撞击、泰勒杆、碎片云及子弹入射靶板等算例，从不同的角度证明了PCMM方法计算高速冲击问题的精确性及合理性。

当然，PCMM方法也存在许多值得进一步探讨的地方，如颗粒质量的等效方式、单元删除与重建的条件、单元重建中力学信息的继承与更新、颗粒接触的快速检索方法、三维PCMM的实现等。

[1] Li S H, Zhao M H, Wang Y N, et al.A continuum-based discrete element method for continuous deformation and failure process[C]∥WCCM VI in Conjunction with APCOM’04.Beijing, 2004:77.

[2] Wang Y N, Zhao M H, Li S H, et al.Stochastic structural model of rock and soil aggregates by continumm-based discrete element method[J].Science in China Series E-Engineering & Materials Science, 2005,48(suppl):95-106.

[3] Littlefield D L.The use ofr-adaptivity with local, intermittent remesh for modeling hypervelocity impact and penetration[J].International Journal of Impact Engineering, 2001,26(1):433-442.

[5] 张雄,刘岩,马上.无网格法的理论及应用[J].力学进展,2009,39(1):1-36.Zhang Xiong, Liu Yan, Ma Shang.Meshfree methods and their applications[J].Advances in Mechanics, 2009,39(1):1-36.

[6] Lucy L B.A numerical approach to the testing of the fission hypothesis[J].The Astronomical Journal, 1977,82:1013-1024.

[7] Gingold R A, Monaghan J J.Smoothed particle hydrodynamics-theory and application to non-spherical stars[J].Monthly Notices of the Royal Astronomical Society, 1977,181:375-389.

[8] Antuono M, Colagrossi A, Marrone S, et al.Free-surface flows solved by means of SPH schemes with numerical diffusive terms[J].Computer Physics Communications, 2010,181(3):532-549.

[9] 肖毅华,胡德安,韩旭,等.一种自适应轴对称FEM-SPH耦合算法及其在高速冲击模拟中的应用[J].爆炸与冲击,2012,32(4):384-392.Xiao Yi-hua, Hu De-an, Han Xu, et al.An adaptive axisymmetric FEM-SPH coupling algorithm and its application to high velocity impact simulation[J].Explosion and Shock Waves, 2012,32(4):384-392.

[10] Ma G W, Wang X J, Ren F.Numerical simulation of compressive failure of heterogeneous rock-like materials using SPH method[J].International Journal of Rock Mechanics and Mining Sciences, 2011,48(3):353-363.

[11] Chen J K, Beraun J E, Carney T C.A corrective smoothed particle method for boundary value problems in heat conduction[J].International Journal for Numerical Methods in Engineering, 1999,46(2):231-252.

[12] Sigalotti L D G, López H.Adaptive kernel estimation and SPH tensile instability[J].Computers & Mathematics with Applications, 2008,55(1):23-50.

[13] Ma S, Zhang X, Qiu X M.Comparison study of MPM and SPH in modeling hypervelocity impact problems[J].International Journal of Impact Engineering, 2009,36(2):272-282.

[14] Zhang X, Sze K Y, Ma S.An explicit material point finite element method for hyper-velocity impact[J].International Journal for Numerical Methods in Engineering, 2006,66(4):689-706.

[15] Ambati R, Pan X F, Yuan H, et al.Application of material point methods for cutting process simulations[J].Computational Materials Science, 2012,57:102-110.

[16] 吕剑,何颖波,田常津,等.泰勒杆实验对材料动态本构参数的确认和优化确定[J].爆炸与冲击,2006,26(4):339-344.Lü Jian, He Ying-bo, Tian Chang-jing, et al.Validation and optimization of dynamic constitutive model constants with Taylor test[J].Explosion and Shock Waves, 2006,26(4):339-344.

[17] Beissel S R, Gerlach C A, Johnson G R.A quantitative analysis of computed hypervelocity debris clouds[J].International Journal of Impact Engineering, 2008,35(12):1410-1418.

[18] Huang J, Ma Z, Ren L S, et al.A new engineering model of debris cloud produced by hypervelocity impact[J].International Journal of Impact Engineering, 2013,56:32-39.

[19] Chi R Q, Pang B J, Guan G S, et al.Analysis of debris clouds produced by impact of aluminum spheres with aluminum sheets[J].International Journal of Impact Engineering, 2008,35(12):1465-1472.