兵器装备工程学报

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

分类:主编推荐 发布时间:2018-12-13 19:42 访问量:2272

分享到:
ü  引用本文

 ____________________________________________________________________________________________________________________________________

Citation format :ZHANG Zhichun, BIAN Qiang, ZHAN Wenhao, et al.High Velocity Impact Simulation Using SPH Contact Algorithm[J].Journal of Ordnance Equipment Engineering,2018,39(9):1-6.

本文引用格式 : 张志春,卞强,展文豪,等.高速冲击问题的SPH粒子接触算法三维数值计算[J].兵器装备工程学报,2018,39(9):1-6.

 ____________________________________________________________________________________________________________________________________


作者简介 : 张志春(1982—),男,博士,助理研究员,主要从事计算力学研究,E-mail:zzchun323@163.com。

___________________________________________


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

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

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

摘要 : 光滑粒子流体动力学方法模拟高速冲击问题时通常采用连续速度算法来计算子弹和靶板的接触,但当子弹和靶板分离时,该算法存在缺陷;为解决高速冲击中的接触问题,基于SPH粒子接触算法,分别对球头钢弹斜冲击钢板和平头钢弹正冲击钢板的过程进行了三维数值计算;通过在SPH粒子上施加接触力的方式计算子弹和靶板发生接触时的接触作用;将SPH粒子接触算法、LS-DYNA970算法、SPH连续速度算法的计算结果和实验结果进行了对比,表明SPH粒子接触算法能够适用于正冲击、斜冲击、复杂接触面的问题。

关键词 : 光滑粒子流体动力学方法;接触算法;斜冲击;正冲击


High Velocity Impact Simulation Using SPH ContactAlgorithm

ZHANG Zhichun, BIAN Qiang, ZHAN Wenhao, LIU Tongyu

(National Key Laboratory of Human Factors Engineering, China Astronauts Research and Training Center, Beijing 100094, China)

Abstract : Contact is often treated with continuity velocity algorithm in Smoothed Particle Hydrodynamics (SPH) for high velocity impact, which inevitably introduces inaccuracies into the calculation when simulating the separation between bullet and target. The SPH contact algorithm was employed to solve the contact in high velocity impact in this paper. Contact forces were enforced on SPH particles once contact occurs, and the identification of the contact surface and surface normal were not required. The oblique impact of target plate against hemispherical-nosed bullet and the normal impact of target plate against blunt-nosed bullet were simulated in 3D with this SPH contact algorithm. The numerical simulation results of the SPH contact algorithm, LS-DYNA970, the continuity velocity algorithm and the experimental observations were compared, showing that the SPH contact algorithm can treat oblique impact, normal impact, and problems with complicated contact surface.

Key words : Smoothed Particle Hydrodynamics; contact algorithm; oblique impact; normal impact


由于在计算大变形问题方面的优势,光滑粒子流体动力学方法(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)

其中 , ρ , , 分别表示SPH粒子 的质量、密度、速度和内能, 表示位置。 ij )表示SPH光滑核函数, 是光滑长度。 Π ij 表示人工黏度,用于消除非物理振荡。

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

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

(5)

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

(6)

是剪切模量, δ αβ 是Kronecker张量,应变率张量  和扭转率张量 αβ 分别为

(7)

(8)

其SPH离散格式分别表示为

(9)

(10)

2 SPH粒子接触算法

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

(11)

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

▽ ij )

(12)

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

(13)

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

图1 SPH连续速度算法

图2 SPH粒子接触力算法

3 球头钢弹斜冲击钢板

本节采用SPH粒子接触算法对球头钢弹斜冲击钢板进行了三维数值计算。子弹和靶板的材料参数为:密度 ρ =7 830 kg/m ,弹性模量 =2.07×10 11 Pa,泊松比 =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。子弹初速 =300 m/s,沿 轴正向,子弹与靶板的初始夹角为60°。为了进行对比验证,采用LS-DYNA970对相同问题进行了计算,子弹和靶板均采用六面体单元离散,单元数分别为1 824和2 700。计算时间均为100×10 -6 s。

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

图3 数值计算结果

图4 子弹和靶板的速度时程曲线

4 平头钢弹正冲击钢板

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

图5 计算模型离散方式及尺寸(单位:m)

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

(14)

其中 是材料常数, 是累积损伤塑性应变,是累积塑性应变,  是用户定义参考应变率。无量纲温度


T*

(15)

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

(16)

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

(17)

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

=[ exp( σ 

(18)

其中 ~ 为材料常数, σ σ σ eq 为应力三轴度, σ =( σ σ σ )/3为平均正应力。子弹和靶板的材料参数分别见表1、表2。

表1 子弹材料参数

表2 靶板材料参数

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

图6 冲击过程数值模拟

图7 冲击后的子弹和靶塞

图8 冲击后的靶板横截面

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

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

5 结论

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

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

参考文献:

[1] 肖毅华,张浩锋,胡德安,等.基于改进PIB搜索法的SPH方法在高速冲击模拟中的应用[J].振动与冲击,2015,34(23):88-92.

[2] 于永强,李成,铁瑛.基于SPH方法的鸟撞复合材料层合板数值分析[J].玻璃钢/复合材料,2017(5):48-52.

[3] 吕东喜,黄燕华,唐永健,等.基于SPH算法的磨粒冲击工件表面过程数值模拟[J].振动与冲击,2013,32(7):169-174.

[4] MONAGHAN J J.On the Problem of Penetration in Particle Methods[J].J Comput Phys,1989,82(11):1-15.

[5] CAMPBELL J,VIGNJEVIC R,LIBERSKY L.A Contact Algorithm for Smoothed Particle Hydrodynamics[J].Comput Method Appl M,2000,184:49-65.

[6] RANDLES P W,LIBERSKY L D.Smoothed Particle Hydrodynamics:Some Recent Improvements and Applications[J].Comput Method Appl M,1996,139:375-408.

[7] VIGNJEVIC R,DE VUYST T,CAMPBELL J C.A Frictionless Contact Algorithm for Meshless Methods[C]//International Conference on Computational & Experimental Engineering and Sciences.Miami,USA ICCES,2007,3(2):107-112.

[8] VIGNJEVIC R,DE VUYST T,CAMPBELL J C.A Frictionless Contact Algorithm for Meshless Methods[J].Comput Model Eng Sci,2006(13):35-48.

[9] MONAGHAN J J.SPH without a Tensile Instability[J].J Comput Phys,2000,159:290-311.

[10] ZHANG Z C,QIANG H F,GAO W R.Coupling of Smoothed Particle Hydrodynamics and Finite Element Method for Impact Dynamics Simulation[J].Engineering Structures,2011,33(1):255-264.

[11] LIU G R,LIU M B.Smoothed Particle Hydrodynamics:A Meshfree Particle Method[M].World Scientific,2004.

[12] MONAGHAN J J.Smoothed Particle Hydrodynamics[J].Rep Prog Phys,2005,68:1703-1759.

[13] 强洪夫,高巍然.修正变光滑长度SPH方法及其应用[J].解放军理工大学学报(自然科学版),2007,8(5):419-424.

[14] 强洪夫,高巍然.完全变光滑长度SPH法及其实现[J].计算物理,2008,25(5):569-575.

[15] B?RVIK T, LANGSETH M, HOPPERSTAD O S.Ballistic Penetration of Steel Plates[J].Int J Impact Eng,1999,22:855-886.

[16] JOHNSON G R,COOK W H.Fracture Characteristics of Three Metals Subjected to Various Strains,Strain Rates,Temperatures,and Pressures[J].Eng Fract Mech,1985,21:31-48.