[ 文章编号 ] 10012246X (2009) 0420548205
高速碰撞 SP H 方法模拟中的初始光滑长度
和粒子间距
徐金中 , 汤文辉
( 国防科学技术大学理学院技术物理研究所 , 湖南 长沙 410073)
[ 摘 要 ] 采用光滑粒子流体动力学方法对高速碰撞问题作数值模拟 , 分析初始光滑长度和非一致粒子间距对
计算结果的影响 , 提出修正光滑长度法 , 并引入 XS PH 速度纠错公式. 数值结果表明 , 初始光滑长度越大 , 弹丸 的刚性越小 , h 0 的合理取值范围应为 d 0 < h 0 ≤115 d 0 ; 粒子间距越小 , 材料的刚性越小. 非一致粒子间距和均 匀粒子分布的计算结果吻合的较好.
[ 关键词 ]
高速碰撞 ; 光滑粒子流体动力学 ; 粒子间距 ; 数值模拟
[ 中图分类号 ]
O385
[ 文献标识码 ]
A
0 引言
光滑粒子流体动力学 ( s m oothed particle hydrodynamics ，简称 SPH ) 方法最初是由 Monaghan 于 1977 年提出 来的一种纯 Lagrange 粒子方法. 该方法的基本思想是将整个流场的物质离散为一系列具有质量 、速度和能量 的粒子 ，然后通过核函数估值 ，求得流场中不同位置在不同时刻的各种动力学量. 由于该方法避免了拉格朗
日网格方法中的网格缠绕和扭曲问题 ，可以广泛应用于大变形计算分析中. 1993 年 Libersky 等1
率先将材料 强度效应引入 SPH 方法 ，成功进行了高速碰撞数值模拟. R andles 2 在此基础上引入了守恒光滑解和边界条 件 ，得到了超高速正碰撞和斜碰撞的合理图像.
在采用 SPH 方法进行高速碰撞模拟时 ，初始光滑长度的选取对计算结果有重要影响. 初始光滑长度太 小 ，可能出现数值断裂等问题 ；初始光滑长度太大会导致计算量的急剧增加 ，核函数影响域内的粒子数太多 , 可能产生压缩不稳定性3 . 一般情况下 ，为了保证计算精度和数值稳定性 ，整个计算区域的粒子间距应该保 持一致. 但是 ，如果研究的问题是三维高速或超高速碰撞 ，一致的粒子间距必然会产生大量的粒子 ，这就对计 算能力提出了更高的要求. 理想的情况是 ，对于感兴趣的区域 ，粒子间距尽量要小 ，而对计算结果影响不大的 区域的粒子间距应相对增大 ，以提高计算效率. 但是 ，这样的粒子布置对计算结果有何影响还不得而知.
本文自编 SPH 程序 ，对高速碰撞问题进行模拟 ，分析初始光滑长度对计算结果的影响. 比较一致粒子间 距和非一致粒子间距的数值结果 ，分析产生差异的原因.
1 SP H 基本方程
考虑材料的弹塑性效应 ，全应力张量空间中的 SPH 插值公式可表示为
1 ,4
d ρi
d t = m i
u ij ·
i W ij , (1) α
d u i
Πij
= - ∑m (2) W ij ,β , d t j d E i
α
= ∑m j ( u i - ij W ij ,β + H i ,
(3)
d t
j
其中 ρ为密度 , m 为粒子质量 ,σ为粒子应力 , u 为粒子速度 , E 为比内能 ,Π 为人工粘性 , H 为人工热流 , W
[ 收稿日期 ] 2007 - 11 - 19 ; [ 修回日期 ] 2008 - 09 - 11 [ 基金项目 ] 国家自然科学基金
(10672180 ,10802097) 资助项目 [ 作者简介 ] 徐金中
(1977 - ) ，男 ，山西原平 ，博士生 ，从事超高速碰撞动力学方面的研究.
为核函数. 下标i 为粒子编号，下标j 为i 粒子的近邻粒子编号.
2核函数
核函数是SPH 算法的一个重要组成部分，其形式多种多样，如B2spline 核函数、G auss 核函数、二次核函数及指数核函数等. B2spline 核函数，其形式为4
≤|≤| z |
z |
z |
< 1 ,
≤2 ,(4) > 2 ,
3
式中z = | x i - x j | Πh h i j ，分别对应一维、二维和三维问题.
中的光滑长度应该随材料密度和空间维数的变化而变化，它有多种形式
(5)其中u 为粒子速度
3 数值模拟及分析
311 不同初始光滑长度下的高速碰撞结果
计算模型为长径比LΠD = 1215 的钨合金长杆弹( L= 510 cm) 高速撞击4185 cm 厚的硬钢靶板，碰撞速度为1168 km·s- 1 . 材料参数如表1 所示. 弹靶粒子间距均为d = 014 mm ，初始光滑长度h 分别取d ,115 d ,
0 0 0 0
2 d
0 ,3 d
, t = 70μs 时计算终止.
表1 材料参数
T a b le 1M aterial para m eters
ρ0Πg·c m - 3c0Πcm·μs - 1γ0
材料s Y/ G Pa G/ GPa 1716
我帮妈妈做家务
7185
01385
01450
闻官军收河南河北的意思
1158
2117
1144
1149
112
1145
钨合金
硬钢
120
77
图1 给出了四种初始光滑长度下的弹靶变形. 从图中可以看出，当h0 = d0 时，靶板中部存在明显的数值
断裂，靶板背面存在数值不稳定性；当h0 = 115 d0 时，数值断裂消失，靶板背面粒子分布更加有序; h0 = 210 d0得到的计算结果与h0 = 115 d0 基本相似; 而当初始光滑长度增大到310 d0 时，弹丸粒子的空间分布变得无序，弹靶界面不再清晰光滑，亦即数值不稳定性又重新产生.
凉拌腐竹的家常做法图1 四种初始光滑长度下的弹靶变形图
Fig11 Deformations o f projectile and targ et w ith four initial sm oothing leng ths
弹丸头部和尾部位置随时间的变化如图2 所示，图中的实验数据见文7 . 从图中可以看出, h0 分别取
d0 ,115 d0 和2 d0 时，计算结果差别不大，而h0 = 3 d0 得到的结果与前三者有明显差别. 弹丸头部位置表明了靶板的变形程度，说明初始光滑长度越大，靶板的变形越小. 头部和尾部位置坐标的变化说明初始光滑长度
叠信封计算物理第26 卷550
越大，弹丸变形越大, 即弹丸的刚性越小. 导致这一
结果的原因可能是由于核函数的影响域内包含了太
多粒子而导致压缩失稳. 如果取h0 = d0 ，在应用可变
光滑长度之后，随着光滑长度和粒子体积的改变，可
能出现h n < d n 的情况, 这就必然导致数值不稳定
性. 从图1 ,2 可以看出, h0 分别取115 d0 和2 d0 得到
的计算结果相差不大，但是h0 = 210 d0 需要花费更多
的计算时间. 因此，在模拟碰撞动力学问题时，为了得
到稳定的数值结果, h0 的合理取值范围应为d0 < h0
≤115 d0 ，这一点可在相关文献中得到证实.
312 非一致粒子间距的高速碰撞结果及分析
伏腊为提高SPH 方法的计算效率, 采用非一致的粒
子间距进行粒子的初始布置. 计算模型同前. 粒子间
图2 弹丸头部和尾部位置随时间的变化
Fig12 No and tail positions versus time
距的选取分两种情况: ①弹丸粒子间距d P0 = 012 mm , 靶板粒子间距d T0 = 014 mm ; ②弹丸粒子间距d P0 = 014 mm ，靶板粒子间距d T0 = 012 mm. 这里定义弹丸靶板粒子间距的比值R PT = d P0 Πd T0 . 图3 为两种粒子间距比值下的弹靶变形. 可以看出，当R PT = 015 ，弹丸的变形较大，剩余长度较小，靶板的变形较小；当R PT = 210 , 靶板的变形较大，弹丸的变形较小，呈现出一定的刚性. 图4 为两种粒子间距比值下弹丸头部和尾部位置随时间的变化. 可以看出，非一致粒子间距对计算结果有较大影响. 粒子间距越小，相应的碰撞体的刚度越小, 反之亦反.
图3 两种粒子间距比值下的弹靶变形图
Fig13 Deformations o f projectile and
targ et at tw o R PT
图4 两种粒子间距比值下弹丸头部和
尾部位置随时间的变化
Fig14 No and tail positions versus time at d ifferent R P T 对于任意粒子i ，其影响域是以该粒子质心为中心，以2 h i 为半径的圆形区域. 那么，如果粒子间距不同，其初始光滑长度就会不同，进而导致粒子的影响域大小不同，从而造成搜索到近邻粒子总质量的不同. 这种近邻粒子总质量的改变，势必造成粒子压力的变化，而弹丸和靶板的变形是与粒子压力紧密相关的. 一般来说，对于同一体(如弹丸) ，若其本身具有一致的粒子间距，那么间距的大小对计算结果影响不大. 因此，造成上面计算结果的原因实际上是因为弹丸和靶板界面两侧粒子分布的不对称性.
不同粒子间距下界面粒子压力随时间的变化如图5 所示. 由图可见，随着R P T的增大，界面弹丸粒子压力逐渐减小，靶板粒子压力逐渐增大，这意味着弹丸的变形是随着R PT 的增大而逐渐减小的，而靶板的变形是随着R P T的增大而逐渐增大的，这就表明粒子间距越小碰撞体的刚度越小.
似组词语从上面界面压力的分析结果看，粒子间距越小，其压力越大，一种可行的方法是对界面粒子的影响域作一定的调整，使得其压力减小，材料刚性增大. 为此，本文对界面粒子的初始光滑长度进行修正:
简笔画车子
图 5 不同粒子间距下界面粒子压力随时间的变化
Fig 15 Pressure o f interfacial particles versus time w ith d ifferent R P T
αPT h i ,
(2 - αPT ) h ,
i ∈ V 1 , i ∈ V 2 ,
(6) h i =
2
αPT
= 01777 + 01270 R P T - 01047 R PT , (7)
其中 V 1 , V 2 分别代表弹丸和靶板.
采用以上的表达式对算例作重新计算. 图 6 给出了修正后的弹丸头部和尾部位置随时间的变化. 从中可以看出 ，修正后的结果得到了明显改善 ，表明本文提出的修正方法可以较好地处理非一致粒子间距物质界面
问题 ，有助于进一步提高 SPH 方法的计算效率. 还可以看出 ，虽然采用修正的光滑长度法可以明显改善计算
结果 ，但是并不能完全消除非一致粒子间距对计算结果带来的影响 ，头部位置坐标的计算结果仍需进一步改善. 为此 ，在使用修正光滑长度的基础上 ，本文尝试引入粒子的 XSPH 速度纠错公式8 ,9 对数值结果进一步修正. XSPH 速度纠错实际上是对粒子 i 与它的近邻粒子进行速度平均 ，其形式如下 :
m i
α
α
α α
u i = u i + β ∑ ( u j - u i ) W ij , (8)
ρ ij
j 这里 ,ρ ij = 015 (ρi +ρj ) ,
β 的取值范围为 0～1 ，一般取 015. XSPH 速度纠错公式是为了消除拉伸不稳定性 ，但它在一定程度上抑制了粒子的运动 ，从而使得速度计
算结果与真实情况不符 ，因此 ，式 (8) 仅用于对本文的修正结果进行微调 ，且仅应用于界面粒子. 对于 R P T =
015 和 210 两种情况 ,β 分别取 01000 5 和 01000 3 . 在引入修正光滑长度和调整系数后的 XSPH 速度纠错公式之后 ，对前面的算例再次计算 ，弹丸头部和尾部位置随时间的变化如图 7 所示. 由图可见 ，弹丸头部位置的修正
客户经理结果比较理想 ，而尾部位置的修正结果不是很理想. 因此 ，对于 XSPH 速度纠错公式的使用仍需进一步探讨.
图 6 修正后的弹丸头部和尾部位置随时间的变化
Fig 16 M od ified no and tail positions versus time
图 7 引入 XS PH 后的弹丸头部和尾部位置随时间的变化
Fig 17 No and tail positions versus time w ith XS PH
计 算 物 理 第 26 卷
552 4 结论
1) 采用 SPH 方法对不同初始光滑长度下的高速碰撞结果比较分析 ，结果表明 ，初始光滑长度越大 ，弹丸
的刚性越小 , h
0 的合理取值范围应为 d 0 < h 0 ≤115 d 0 . 2) 对非一致粒子间距的高速碰撞问题进行模拟 ，计算表明 ，弹丸靶板粒子间距比值 R P T 越小 ，弹丸的变
形越大 ，即粒子间距越小 ，材料的刚性越小.
3) 对导致非一致粒子间距和均匀粒子布置的计算结果产生明显差异的原因进行分析 ，提出修正的光滑
长度法 ，将该方法引入数值模拟 ，得到了很好的结果. 因此 ，本文提出的修正光滑长度法是处理非一致粒子间
距高速碰撞问题的一种有效方法.
4) 通过进一步引入 XSPH 速度纠错公式并对其中的系数进行调整 ，使得弹丸头部位置的修正结果比较
理想 ，但尾部位置的修正结果不够理想 ,XSPH 速度纠错公式的使用问题仍需进一步研究.
[ 参 考 文 献 ]
Lib ersky L D , Petschek A G . Hig h strain Lag rang ian hydrod ynamicsJ . Journal of C
om p u tational Ph ysics , 1993 , 109 :67 - 71. Benz W , Asp hang E. S imu lation of brittle s olids using sm oothed particle hydrod ynamicsJ . C omp u ter Physics C ommunication , 1995 ,
87 :253 - 265.
S wegle J W , Attaw ay S W , Heinstein M W , et al . An analysis of sm oothed particle hydrod ynamics R . Albuqu erque : S an d ia National Lab oratories , S AND93 - 2513 , 1994.
Petschek A G , Lib ersky L D. C ylindrical sm oothed particle hydrod ynamicsJ . Journal of C omp u tational Physics , 1993 , 109 : 76 - 80.
B ø
rve S , Omang M , Tru ln J . R egu larized sm oothed particle hydrod ynamics w ith improved mu lt 2i resolu tion hand lingJ . Journal of C omp u tational Ph ysics , 2005 , 208 :345 - 367.
Benz W. Nu merical m od eling of nonlinear s
tellar p u lsation : prob lems and prospectsR . K luw er Acad emic , B oston , MA , 1990 : 268 - 269.
And erson J r C E , H oler V , Walker J D , S tilp A J . Tim e 2resolved penetration of long rods into steel targ etsJ . In t J I mpact Eng ng ,
1995 ,16 :1～18.
M onag han J J . S PH w ithou t a tensile instab ilityJ . Journal of C
om p u tational Ph ysics , 2000 , 159 :290 - 311. M onag han J J . On the prob lem of penetration in par ticle m ethodsJ . Journal of C om p u tational Ph ysics , 1989 , 82 :1 - 15.
1 2 ] ] 3 ] 4 ]
5 ]
6 ]
7 ]
8
9 ] ] I nitial Smoothing Length and Sp ace Bet w een Particles in SP H Method f o r
N umerical Simu lation of H igh 2speed Imp acts
X U Jinzhong , TANG Wenhu i
( Institute o f Te chnical P hysics , C ollege o f Science , N a tional Univ ersity o f D efen Technology , Changsha 410073 , C hina )
Abstract :
Hig h 2speed impact prob lems are simu lated w ith sm oothed particle hydrod ynamics method. Initial sm oothing leng th and
noncon forming space b etw een particles in simu lation are analyzed. A m od ified sm oothing leng t
h approach is p u t forw ard. With increa of
initial sm oothing leng th , projectile rig id ity d ecreas g radu ally. A rational rang e of initial sm oothing leng th is d 0 < h 0 ≤1. 5 d 0 . W ith d ecrea
of space b etw een particles , rigid ity of material d ecreas. With XS PH correction on velocity of a particle , resu lts b y noncon forming space b etw een particles ag ree w ell w ith tho fr om id entical particle d istribu tion.
K ey w ords : hig h 2speed impacts ; sm oothed particle hydrod ynamics ; space b etw een particles ; nu merical simu lation
R eceived date : 2007 - 11 - 19 ; R evid date : 2008 - 09 - 11
