高速冲击问题的SPH粒子接触算法三维数值计算

张志春，卞 强，展文豪，刘桐宇

(中国航天员科研训练中心人因工程国防科技重点实验室， 北京 100094)

由于在计算大变形问题方面的优势，光滑粒子流体动力学方法(Smoothed Particle Hydrodynamics, SPH)在高速冲击问题的计算中，得到了广泛的应用。肖毅华[1]基于改进PIB搜索法的SPH方法模拟了泰勒杆的冲击问题；于永强[2]使用SPH方法计算了鸟撞复合材料层合板问题；吕东喜[3]使用SPH方法对磨粒冲击工件表面过程进行了数值模拟。

使用SPH方法模拟高速冲击问题，子弹和靶板的接触是计算的关键。传统上采用SPH连续速度算法模拟弹体和靶板的接触。连续速度算法基于SPH粒子求和思想，在求解SPH动量方程时，将支持域内所有的靶板和弹体粒子均考虑在内，实现子弹和靶板的相互作用。由于SPH连续速度算法不能区分不同材料的SPH粒子，当子弹和靶板靠近时，该方法能够描述两者的“排斥”作用；但当子弹和靶板分离时，该方法却使两者相互“吸引”，明显与实际不符合。同时连续速度算法很容易使接触面两侧粒子相互“融合”，不同材料的粒子侵入彼此“领域”，造成接触界面模糊，难以追踪。

为了描述SPH粒子之间的接触作用，Monaghan[4]提出了XSPH算法，包含了相邻粒子的影响作用，使粒子的速度与相邻粒子的平均速度相近。该方法存在如下缺陷：不能计算拉伸力，确保两物体分离；不能计算剪切力。Campbell等[5]提出了一种粒子-粒子罚接触算法，使用Randles和Libersky[6]提出的方法来确定边界，在三维计算中较为复杂。Vignjevic等[7-8]借鉴Monaghan[9]SPH粒子排斥力算法的思想，提出了一种SPH粒子接触算法，通过在SPH粒子上施加接触力的方式实现SPH粒子间的接触作用，取得了较好的效果，文献[10]对该算法进行了应用。

本文基于SPH粒子接触算法，采用C++编写的SPH计算程序，分别对球头钢弹斜冲击钢板以及Arne tool钢制圆柱子弹正冲击Weldox 460 E钢板发生冲塞破坏的过程进行了三维数值计算，通过在子弹和靶板相关粒子上施加接触力实现两者之间的接触。将SPH粒子接触算法的计算结果、LS-DYNA970的计算结果、SPH连续速度算法的计算结果、实验结果进行了对比，验证了SPH粒子接触算法在模拟高速冲击问题中的有效性。

1 考虑材料强度的SPH方程

在具有材料强度的高速冲击问题中，冲击波沿着碰撞体内传播，此时碰撞体具有流体特性。其中运动方程和高压状态方程是控制材料力学行为的关键，并且材料强度只有在这种具有能量传动问题的后期才能显示出其重要性。由SPH离散的具有材料强度的流体动力学控制方程如下[11-12]：

连续性方程：

(1)

动量方程：

(2)

能量方程：

(3)

运动方程：

(4)

其中mi，ρi，vi，ei分别表示SPH粒子i的质量、密度、速度和内能，x表示位置。Wij=W(xi-xj,h)表示SPH光滑核函数，h是光滑长度。Πij表示人工黏度，用于消除非物理振荡。

材料强度体现在应力张量σαβ，由各向同性压力p和黏性剪切应力ταβ组成

σαβ=-pδαβ+ταβ

(5)

其中偏应力率采用Jaumann比率的形式表示为

(6)

(7)

(8)

其SPH离散格式分别表示为

(9)

(10)

2 SPH粒子接触算法

传统的处理SPH接触问题的算法如图1所示，其中小实线圆表示SPH粒子，大虚线圆表示粒子的支持域。当物体A、B靠近时，采用SPH动量方程对A、B内粒子速度进行更新。进行SPH连续速度算法时，粒子i的支持域包含A内的粒子n4,n5,…,n8和B内的粒子n1,n2,n3。Vignjevic等[7-8]提出的SPH无摩擦接触算法通过在SPH粒子上施加接触力的方式实现SPH粒子间的接触作用。如图2所示，当A、B间距达到两倍光滑长度时，在相关SPH粒子上产生接触力(SPH粒子支持域的半径选为两倍光滑长度)；当A、B间距超过2倍光滑长度时，没有接触力产生。首先定义一个接触势函数

(11)

其中ncont表示SPH粒子i支持域内属于不同体的粒子数，如图2中，对于粒子i，ncont=3，即计算粒子i的接触力时其支持域仅包含物体B内的粒子n1,n2,n3。当xA和xB属于同一体时，SPH核函数W(xA-xB,h)=0。Δdavg表示粒子间距的平均值。K和n表示用户定义参数，K的选取与有限元接触中的罚刚度类似，与材料性质和冲击速度相关。该接触势函数具有如下特点：在物体内部为0；通常为正值；随粒子间距的减小而增大。对该接触势函数的梯度进行SPH格式离散，得到施加在粒子上的接触力

▽xiW(rij)

(12)

接触力的方向由接触势函数的梯度确定。将式(12)代入式(2)，得到含有接触力的SPH动量方程

(13)

采用连续速度算法处理SPH的接触问题，当物体A、B分离时，由于粒子求和的作用，A、B粒子会产生相互“吸引”，与实际情况不符合。而接触力算法能够客观的反应A、B之间的接触作用，较为真实的描述A、B的靠近和分离。该粒子接触算法与SPH的求解格式保持一致，可以充分利用SPH的临近搜索列表。只要当两物体接触界面粒子间距达到两倍光滑长度，就在相关粒子上施加接触力，当接触界面粒子间距超过两倍光滑长度时，没有接触力，因此不需要定义接触粒子。该算法和有限元接触算法不同，在每个时间步不需要定义接触界面和界面法线方向，便于进行三维程序设计。

3 球头钢弹斜冲击钢板

本节采用SPH粒子接触算法对球头钢弹斜冲击钢板进行了三维数值计算。子弹和靶板的材料参数为：密度ρ=7 830 kg/m3，弹性模量E=2.07×1011Pa，泊松比v=0.3。子弹长0.078 m，直径0.018 m，靶板尺寸为0.09 m×0.09 m×0.009 m，子弹和靶板的SPH粒子分别为2 740和2 700，粒子间距均为0.003 m。子弹初速v0=300 m/s，沿Z轴正向，子弹与靶板的初始夹角为60°。为了进行对比验证，采用LS-DYNA970对相同问题进行了计算，子弹和靶板均采用六面体单元离散，单元数分别为1 824和2 700。计算时间均为100×10-6s。

在t=70×10-6s时，SPH接触算法和LS-DYNA970数值计算结果如图3所示。在冲击载荷作用下，靶板发生较大的变形，而子弹变形相对较小。子弹和靶板沿Z轴的速度时程曲线如图4所示，表明两种算法的计算结果吻合较好。两种方法的误差主要有两方面的原因：由于离散方式的不同，在材料密度相同的前提下，两种方法中子弹和靶板的质量存在一定差异；SPH粒子接触算法使子弹和靶板发生接触的时间比LS-DYNA970要早。通过减小SPH粒子的初始间距，可以一定程度上减小这两种原因造成的误差。通过球头钢弹斜冲击钢板的数值计算表明：SPH粒子接触算法适用于斜冲击和接触界面形状复杂的问题。

4 平头钢弹正冲击钢板

本节采用SPH粒子接触算法对圆柱形钢弹正冲击钢板发生冲塞破坏的过程进行了三维数值计算，计算模型的尺寸及离散方式如图5所示。子弹和靶板均由SPH粒子离散，粒子总数57 728。子弹直径0.02 m，高0.08 m，粒子间距为10-3m。靶板尺寸0.28 m×0.28 m×0.012 m，中心0.04 m×0.04 m×0.012 m区域内粒子间距10-3m，为减小计算量，粒子间距从中心正碰区向边缘逐渐增大。子弹与靶板初始间距2×10-3m，靶板四周采用固支边界条件，计算总时间160×10-6s。采用Monaghan[9]提出的人工应力法来消除拉伸不稳定性造成的数值断裂，采用强洪夫[13-14]提出的完全变光滑长度法来消除光变光滑长度带来的影响。

子弹材料为Arne tool钢，冲击过程中变形较小，采用线弹性模型。靶板材料为Weldox 460 E钢，冲击中发生冲塞型破坏，变形较大，为描述靶板材料的屈服应力及损伤演化，采用BØrvik等[15]提出的修正Johnson-Cook强度模型和Gruneisen状态方程。该模型中将材料的屈服强度表示为损伤变量、等效应变、等效应变率和温度的函数

(14)

(15)

T0是室温，Tm是材料熔点。假设为绝热条件，靶板材料温度的升高是由于冲击过程中的塑性功转换为热量造成的，温度的变化计算公式为

(16)

Cp是材料比热，α是比例常数，Wp是塑性功。D为损伤变量，D=0表示材料没有损伤，D=1表示材料完全失效。实际上材料出现宏观裂纹时，损伤变量的临界值小于1，失效准则可描述为D=DC≤1，其中DC是损伤变量临界值。损伤变量D是累积塑性应变p的函数，可描述如下：

(17)

其中pd是损伤阈值，pf是断裂等效塑性应变，与材料的应力三轴度、应变率和温度相关。Johnson和Cook在文献[16]提出的剪切损伤演化模型中将pf描述如下

(18)

其中D1～D5为材料常数，σ*=σm/σeq为应力三轴度，σm=(σx+σy+σz)/3为平均正应力。子弹和靶板的材料参数分别见表1、表2。

表1 子弹材料参数

表2 靶板材料参数

图6显示了子弹初速为285.4 m/s时的数值计算结果，而相关的实验结果如图7、图8所示[15]。在t=7×10-6s，子弹与靶板开始接触，子弹正前方区域靶板开始加速，子弹头部周围的靶板出现剪切区。靶板被挤压变形，背部出现凸起。由于子弹-靶板接触界面的塑性应变，剪切区内粒子出现损伤。靶板变形继续增大，粒子损伤迅速演化，部分粒子损伤达到临界值，开始失效，形成裂纹并逐渐向靶板背部扩展。在冲击后期，靶板的失效模式包含靶板内部的剪切断裂和靶板背部的拉伸损伤。最终形成靶塞，完全从靶板脱离。从中可以看出子弹、靶塞和靶板冲塞型破坏的计算结果均与实验吻合较好。子弹变形很小，表明数值计算中子弹选择线弹性模型可行。

表3给出了SPH粒子接触算法、SPH连续速度算法和实验结果[15]，其中v0表示子弹初速，vpr表示子弹余速，vtr表示靶塞速度，Δvpr=(vpr计算值-vpr实验值)/vpr实验值，Δvtr=(vtr计算值-vtr实验值)/vtr实验值。对比显示SPH粒子接触算法的计算精度要明显高于SPH连续速度算法，表明该方法能够更为真实的描述子弹和靶板的接触问题。并且计算精度随子弹初速的增大而提高，表明SPH粒子接触算法对于高速冲击领域具有较好的适用性。

表3 数值模拟结果与实验结果

5 结论

由于在处理大变形方面具有明显的优势，SPH方法被成功地应用于高速冲击问题的数值计算。为了计算不同物体之间SPH粒子的接触作用，本文将SPH粒子接触算法应用于高速冲击。该算法是通过在SPH粒子上施加接触力的方式实现，不需要定义接触界面以及接触界面的法线方向，便于实现三维编程。接触力的计算与接触物体的相对运动方向、接触面的几何形状无关，因此该接触算法对正冲击、斜冲击、复杂接触界面都有较好的适用性。算例中与其他数值方法以及实验结果的对比显示了该SPH粒子接触算法的有效性。

将SPH粒子接触算法应用于高速冲击问题，当子弹和靶板间距达到两倍光滑长度时，在相关SPH粒子上施加接触力。因此接触作用产生的时机比实际情况早，产生一定的误差。为了减小该误差，在前处理建模时，要加密接触部位的SPH粒子设置，使接触作用的计算与实际情况更加相符。