基于改进PIB搜索法的SPH方法在高速冲击模拟中的应用

肖毅华，张浩锋，胡德安，平学成

(1.华东交通大学机电工程学院，南昌330013;2.湖南大学机械与运载工程学院，长沙410082)

基于改进PIB搜索法的SPH方法在高速冲击模拟中的应用

肖毅华1，张浩锋1，胡德安2，平学成1

(1.华东交通大学机电工程学院，南昌330013;2.湖南大学机械与运载工程学院，长沙410082)

SPH方法在高速冲击模拟中具有优势并已得到广泛使用，但其相邻粒子搜索非常耗时，其实际应用受计算效率制约。为了提高计算效率，采用一种改进的Point－In－Box(PIB)搜索法进行相邻粒子搜索。利用基于该相邻粒子搜索法的SPH方法，模拟了泰勒杆、平头弹撞击钢板等一系列高速冲击问题。计算结果表明:SPH方法能有效地模拟冲击过程中的大变形、绝热剪切破坏等现象，模拟结果与实验结果符合良好;同时，与传统相邻粒子搜索法相比，应用改进的PIB搜索法使SPH方法的计算效率得到明显提高。

高速冲击;SPH方法;相邻粒子搜索

SPH(Smoothed Particle Hydrodynamics)方法是一种无网格拉格朗日粒子法，具有良好的大变形计算能力，易于处理断裂、破碎等复杂物理现象，因而在高速冲击模拟中具有优势并得到广泛应用［1－2］。SPH方法无需复杂的数值积分，其场函数的近似可通过相邻粒子进行简单求和而得到，故其计算效率比一般的无网格法高。但是，它需搜索每个粒子的相邻粒子，这是一个很耗时的过程，是制约其计算效率的关键。对于实际高速冲击问题的SPH模拟，往往所需的粒子数比较多，采用高效的相邻粒子搜索技术对于控制计算成本具有重要意义。目前，已有一些比较有效的搜索算法被提出和应用。Monaghan［3］采用链表搜索法实现相邻粒子搜索，该方法的复杂度可以达到O(N)(N为粒子数)，但是它不能适应变光滑长度的情况。Hernquist等［4］提出用树形搜索法进行相邻粒子搜索，该方法对于大规模和变光滑长度问题都具有较高的效率和良好的稳定性，是一种综合性能较好的搜索法。Attaway等［5］应用Point－In－Box(PIB)搜索法完成相邻粒子搜索，该方法比较高效而且非常节省存储空间，但其效率受粒子分布影响，即:当存在某一坐标方向的粒子数很少时，其效率很高，而各坐标方向粒子数都较多时，其效率会明显下降。

为了提高高速冲击模拟的计算效率，本文采用了一种改进的PIB搜索法即条形化PIB搜索法，该方法具有比现有搜索法更优越的性能。应用基于该搜索法的SPH方法，模拟了两个典型的高速冲击算例，通过与实验和其他数值方法的模拟结果对比，验证了计算结果的有效性，同时对比了基于不同相邻粒子搜索法的SPH方法的计算效率。

1 模拟高速冲击过程的轴对称SPH方法

本文仅考虑具有轴对称特性的高速冲击问题，采用二维轴对称SPH方法进行模拟以提高计算效率。轴对称SPH方法的基本方程采用文献［6］中按直接离散方式导出的方程，并引入人工黏性和人工应力来稳定计算和抑制数值断裂。具体的求解格式如下:

式中:t为时间，ρ为密度，m为质量，v为速度，σ为应力，u为单位质量的内能，r、z和θ分别为径向、轴向和环向坐标，i和j表示粒子编号，η=2πρr，W为三次样条核函数，Q为人工黏性，Rfn为人工应力项。人工黏性Q和人工应力项Rfn的具体计算式分别参考文献［2］和文献［7］。

为避免产生虚假的剪切应力和拉伸应力、引起严重的计算误差，本文采用了一种粒子－粒子接触算法显式地处理物体间的接触界面。同时，为了提高接触计算效率，本文的边界粒子识别采用了Randles等［8］提出的方法，即根据物体“颜色”的核近似来确定。

2 条形化PIB搜索法

条形化PIB搜索法是通过改进传统PIB搜索法而得到的。下面先简要介绍传统PIB搜索法的原理。首先对所有粒子各个坐标方向的坐标值分别进行排序(见图1)。当搜索给定矩形区域内的粒子时，根据排序的坐标采用二分法快速查找得到位于x、y两个方向的边界间的粒子，得到两个粒子列表。最后，求出两个粒子列表的交集即为给定矩形区内的粒子。求粒子列表交集的方法为:遍历粒子数最少方向的粒子列表中的每个粒子，判断其在其他方向的坐标排序序号是否在相应方向的最小和最大坐标排序序号之间。例如，假设y方向的粒子列表中的粒子数最少，则检查粒子

那么点i在给定的矩形区域内，否则不在其中。在式(5)和式(6)中，I和R是记录坐标排序结果的数组，数组元素Ix(i)记录x坐标排序序号为i的粒子号;数组元素Rx(i)记录粒子i的x坐标排序的序号，Iy和Ry针对y坐标设置，功能分别与Ix和Rx相同;和分别为构造粒子列表时通过二分法得到的第1个和最后1个位于x方向边界间的粒子的x坐标排序序号，和对应于y方向。

图1 传统PIB搜索法示意图Fig.1 Schematic diagram of traditional PIB search algorithm

确定粒子列表交集是PIB搜索法中最耗时的部分。对于条状分布粒子集，始终有一个方向的粒子列表中的粒子数很少，按式(5)和式(6)确定粒子列表交集的计算量很小，故PIB搜索法对于此类分布的粒子集搜索效率很高。条形化PIB搜索法就是利用这一特点而提出的。

图2 条形化PIB搜索法示意图Fig.2 Schematic diagram of traditional PIB search algorithm

该算法的基本原理见图2。它把任意分布的粒子集划分为一些条形分布的子集。划分子集时，先计算出所有粒子各方向的最大、最小坐标值，得到一个恰好包含粒子集的矩形;再以该矩形尺寸较大的方向作为条形方向，并确定一个条形尺寸Δs，将该矩形划分为若干条状区域，每一个条状区域内的粒子为一个子集。这些子集沿坐标方向顺序编号，且每个子集内的粒子坐标值按各个方向独立排序。当搜索某个给定矩形区域内的粒子时，首先根据该矩形域的两条边界位置坐标，计算出其所在的两个子集的编号，它们为可能包含位于该矩形域内的粒子的首、末两个子集。然后，利用与传统PIB搜索法完全一样的粒子列表构造方法和粒子列表交集确定方法，得到这两个子集中位于该矩形域内的粒子;首、末两个子集之间的粒子集也包含位于该矩形域内的粒子，对于这些子集，只需利用传统PIB搜索法中的粒子列表构造方法确定沿条形方向的粒子列表，此粒子列表中的粒子都是位于该矩形域内的粒子。最后，合并在每个子集中搜索到的粒子即得到该矩形域内的所有粒子(见图2)。在本文中，条形尺寸大小按式(7)确定

式中:hmax和hmin分别为所有粒子的最大和最小光滑长度。

应用条形化PIB搜索法在SPH方法中进行相邻粒子搜索时，首先以粒子i为中心构造边长为4hi的正方形，然后按前述过程搜索得到该正方形内的粒子。如果搜索到的粒子j满足相邻粒子条件则将粒子i和粒子j记录为一个相邻粒子对。

条形化PIB搜索法能有效地克服传统PIB搜索法对粒子分布形式敏感的缺陷，同时具有很高的搜索效率。下节将通过算例分析与传统搜索法中综合性能较好的树形搜索法进行对比，以验证其优越性。

3 算例分析

3.1 泰勒杆实验模拟

本算例模拟House等［9］开展的泰勒杆实验，即直径为7.595 mm、长度为37.97 mm的4340钢圆柱杆撞击刚性壁面，模拟了三种不同初始冲击速度的实验情况，分别为181 m/s、224 m/s和270 m/s。模拟结果与实验结果［9］以及Batra和Zhang［10］的模拟结果进行了对比。计算模型中，圆柱杆用均匀分布的粒子离散，径向粒子数为20，轴向粒子数为200，总共有4 000个粒子;刚性壁面通过镜像虚粒子方法进行处理;4340钢的材料特性采用Johnson－Cook模型描述，模型相关参数参考文献［10］。

图3为三种不同初始冲击速度下杆在60μs时的变形。该时刻杆的动能已经很小，杆的冲击端基本不再产生塑形变形，只有杆的尾部仍沿轴向做微小幅度的振荡。以此时的杆件形状作为SPH模拟的杆的最终形状，得到三种情况下杆的最终长度和冲击面直径，结果见表1。由该表可见，三种工况下SPH模拟结果与实验结果均符合良好，冲击面直径的最大绝对误差仅0.3 mm，最终长度的最大绝对误差也仅0.8 mm。与MSPH相比，SPH在冲击面直径方面的计算精度要略高，而在最终长度方面的计算精度稍低。

图3 不同冲速度下杆的变形Fig.3 Deformations of bars at different impact velocities

图4对比了基于三种不同相邻粒子搜索方法的SPH模拟泰勒杆实验时一个计算循环内各部分的平均计算时间。由图4可知，相邻粒子搜索是SPH模拟中最耗时的部分。采用树形搜索法时，相邻粒子搜索时间占总时间的比例达66%。采用PIB搜索法时，相邻粒子搜索时间减少约2/3，其占总时间的比例降为42%。采用条形化PIB搜索法时，相邻粒子搜索时间略少于基于PIB搜索法的搜索时间，其占总时间的比例为41%。在本算例中，基于条形化PIB搜索法的SPH效率最高，略高于基于PIB搜索法的SPH，远高于基于树形搜索法的SPH，与其相比节省了大概45%的总计算时间。PIB搜索法在本算例中也表现出很高的计算效率，这主要由于泰勒杆沿径向的粒子数目较少(仅为20)。

表1 泰勒杆实验模拟结果Tab.1 Simulation results of Taylor tests

图4 泰勒杆模拟中一个计算循环内各部分的平均计算时间Fig.4 Average computational time for each part of a calculation cycle inTaylor test simulation

3.2 平头弹撞击钢板模拟

本算例模拟Børvik等［11］的平头弹撞击金属板冲塞破坏实验，即直径为20 mm、长度为80mm的Arne工具钢平头弹正撞击直径为500 mm、厚度为8 mm的Weldox 460E钢板。模拟工况的初始冲击速度为173.7m/s。图5为模拟该实验的计算模型。模型共使用104 047个粒子，弹体和靶板中心2倍弹体半径以内的区域采用均匀粒子离散，粒子间距为0.1 mm;靶板外围区域采用沿径向向外粒子间距按等比例逐渐增大的粒子离散，当粒子间距达到约1 mm时再采用恒定粒子间距离散，保证厚度方向有足够多的粒子。弹体材料采用双线性硬化模型，靶板材料采用耦合黏塑性和延性损伤的本构模型［12］。

图5 平头弹撞击金属板的计算模型Fig.5 Computationalmodel for steel plate impacted by blunt projectile

图6为SPH模拟得到的3个典型时刻的弹体和靶板的变形图像。在整个过程中，弹体变形很小。靶板发生绝热剪切变形，产生冲塞破坏。它在弹体半径附近发生了大的塑性变形并产生裂纹，裂纹沿厚度方向发展，最终穿透靶板的整个厚度，形成一个半径与弹体半径大致相等的圆柱塞块而被弹体冲出。上述模拟得到的现象与实验观察结果基本相符。图7给出了弹体的速度时间历程曲线。计算的弹体速度变化与实验结果趋势大体相同，在冲击过程前期计算的弹体速度比实验结果衰减稍快，但最终的弹体剩余速度相差很小。上述结果说明:SPH方法有效地模拟绝热剪切失效过程、预测弹体的剩余速度。

图6 不同时刻弹体和靶板的变形Fig.6 Deformations of projectile and target plate at different time

图8对比了基于三种不同相邻粒子搜索法的SPH模拟平头弹撞击钢板时一个计算循环中各部分的平均计算时间。在本算例中，粒子总数多且沿径向和轴向粒子数都较多，同时光滑长度在空间上变化也大(1个数量级)。由图8可知，相邻粒子搜索仍是最耗时的计算部分。采用树形搜索法时，相邻粒子搜索时间占总时间的比例为60%。采用PIB搜索法时，相邻粒子搜索时间比树形搜索法减少45%，占总时间的比例为44%，虽然搜索效率仍有所提高，但提高幅度比上一算例有大幅下降，这主要是受到粒子分布和总数的影响。而采用条形化PIB搜索法时，相邻粒子搜索时间仍比树形搜索法减少约2/3，且比PIB搜索法减少了41%，占总时间的比例为32%。在本算例中，基于条形化PIB搜索法的SPH的计算效率明显高于基于PIB搜索法和树形搜索法的SPH，与其相比分别节省了20%和40%左右的总计算时间。

图7 弹体的速度变化历程Fig.7 Velocity history of projectile

图8 平头弹撞击钢板模拟中一个计算循环内各部分的平均计算时间Fig.8 Average computational time for each part of a calculation cycle in simulation of plate impacted by blunt projectile

4 结论

介绍了一种基于条形化PIB搜索法的轴对称SPH方法，应用其成功模拟了泰勒杆和平头弹撞击钢板两个问题。结果表明，SPH方法可以较准确地模拟泰勒杆变形和钢板的绝热剪切破坏，能有效地预测弹体的剩余速度。文中对计算效率进行了对比分析，结果表明:基于条形化PIB搜索法的SPH有很高的计算效率，与基于树形搜索法的SPH相比可以大大减少相邻粒子搜索时间，使整体计算效率得到明显提高，与基于PIB搜索法的SPH相比能克服其效率受粒子分布影响的缺陷，在计算效率方面具有良好的稳定性。此外，基于条形化PIB搜索法的SPH对于大规模和光滑长度空间不均的问题具有良好的适应性。

［1］Libersky L D，Petschek A G，Carney TC，etal.High strain Lagrangian hydrodynamics:a three－dimensional SPH code for dynamic material response［J］.Journal of Computational Physics，1993，109(1):67－75.

［2］Johnson G R，Beissel S R.SPH for high velocity impact computations［J］.Computer Methods in Applied Mechanics and Engineering，1996，139(1/2/3/4):347－373.

［3］Monaghan J J.Particle methods for hydrodynamics［J］.Computer Physics Reports，1985，3(2):71－124.

［4］Hernquist L，Katz N.TREESPH－A unification of SPH with the hierarchical tree method［J］.The Astrophysical Journal Supplement Series，1989，70:419－446.

［5］Attaway SW，Heinstein M W，Mello F J，et al.Coupling of smooth particle hydrodynamics with PRONTO［R］.SAND93－1781C.Albuquerque:Sandia National Laboratories，1993.

［6］杨刚.光滑粒子法的改进及其若干典型应用［D］.长沙:湖南大学，2011.

［7］Xiao Y H，Hu D A，Han X，et al.Simulation of normal perforation of aluminum plates using axisymmetric smoothed particle hydrodynamics with contact algorithm［J］.International Journal of Computational Methods，2013，10(3)1350039 (21pages).

［8］Randles PW，Libersky L D.Smoothed particle hydrodynamics:some recent improvements and applications［J］.Computer Methods in Applied Mechanics and Engineering，1996，139(1/2/3/4):375－408.

［9］House JW，Lewis JC，Gillis P P，et al.Estimation of flow stress under high rate plastic deformation［J］.International Journal of Impact Engineering，1995，16:189－200.

［10］Batra R C，Zhang G M.Modified smoothed particle hydrodynamics(MSPH)basis functions for meshless methods，and their application to axisymmetric Taylor impact test［J］.Journal of Computational Physics，2008，227:1962－1981.

［11］Børvik T，Hopperstad O S，Berstad T，et al.Numerical simulation of plugging failure in ballistic penetration［J］.International Journal of solids and structures，2001，38(34/35): 6241－6264.

［12］Børvik T，Hopperstad O S，Berstad T，et al.A computationalmodel of viscoplasticity and ductile damage for impact and penetration［J］.European Journal of Mechanics－A/Solids，2001，20(5):685－712.

Application of SPH w ith modified PIB search algorithm to high-velocity impact simulation

XIAO Yi－hua1，ZHANG Hao－feng1，HU De－an2，PING Xue－cheng1
(1.College of Mechanical and Electrical Engineering，East China Jiaotong University，Nanchang 330013，China; 2.College of Mechanical and Vehicle Engineering，Hunan University，Changsha 410082，China)

SPH method hasmany advantages and iswidely used in high－velocity impact simulation，but its practical applications are still restricted by its low computational efficiency due to the time－consuming neighboring particle search (NPS).To improve its efficiency，a modified point－in－box(PIB)search algorithm was employed to determine neighboring particles.By using the SPH with modified PIB search algorithm，a series of impact problems，including Taylor tests and steel plate impacted by blunt projectile，were simulated.The calculated results indicated that SPH can effectively simulate large deformations and adiabatic shear failure during the impact processes，and give results in good agreementwith experiments.The results also showed that the modified PIB search algorithm can make the SPH method achieve amuch better efficiency than a traditional search algorithm.

high－velocity impact;SPH method;neighboring particle search

O353.4;O242.1

A

10.13465/j.cnki.jvs.2015.23.016

国家自然科学基金项目(11302077);江西省教育厅科学技术研究项目(GJJ14398)

2014－09－12修改稿收到日期:2014－11－06

肖毅华男，博士，讲师，1984年生

平学成男，博士，教授，1975年生