Page 87 - 《爆炸与冲击》2026年第4期

第 46 卷 刘晏东，等： 超高速撞击玄武岩材料的Riemann-SPH仿真参数分析与验证 第 4 期

1 仿真模型

1.1 SPH 仿真模型
1.1.1 Riemann-SPH 求解算法
对于具有材料强度的高速冲击问题，材料行为类似于流体，此时流体动力学的控制方程（governing
equations）宜选用拉格朗日描述下的 Navier-Stokes （NS）方程。NS 方程包含连续性方程（质量守恒方
程）、动量守恒方程、能量守恒方程和运动方程。类似于一般的 SPH 方法，Riemann-SPH 采用具体的物
理状态处理公式近似求解 NS 方程，并使用黎曼求解器引入粒子间接触界面的速度 v ∗ 以描述冲击波面与
不同粒子间的相互作用。在保证能量、动量守恒的前提下，Riemann-SPH 方法能更准确地描述接触界
面，有利于解决在密度不连续处的人工表面张力问题 ：
[6]
v i +v j ( ) p
v ∗ = b ij + 1−b ij e v ∗ +v ij (1)
2

®
if same material
1− D i − D j
b i j = (2)
0 otherwise
p
b ij 为界面耦合系数，v 和 i v 分别为粒子 i 和 v ij 为使用 HLLC
j
式中： j 的速度， e v ∗ 为各自波速加权的平均速度，
[6] D j 分别
（Harten-Lax-van Leer contact）黎曼求解器求解局部接触面黎曼问题得到的压力解附加项 ， D i 和
为粒子 i 和 j 的的损伤程度。
Riemann-SPH 方法将 SPH 方法中的粒子间相对速度替换为粒子间接触界面的速度 v ∗ ，得到求解公式：
 ∑
dρ i m j
 ‹
 = 2ρ i (v i −v ∗ )·∇W i
 dt ρ j

 j

 αβ αβ

 dv i ∑ σ i +σ j

 = ‹
 m j ·∇W i
dt ρ i ρ j
j (3)
αβ
 ∑
 de i σ i

 = −2 m j (v i −v ∗ )·∇W i
‹

 dt
 ρ i ρ j
 j


 dx i

= v i
dt
ρ j 分别为粒子 i 和 m j 分别为粒子 i 和 αβ αβ 分别为粒子 i 和 j 的应
式中： ρ i 和 j 的密度， m i 和 j 的质量， σ i 和 σ j
e i 为粒子 x i 为粒子 ∇W i 为粒子 i 的修正核函数梯度。
力张量， i 的能量， i 的位置， ‹
式 (3) 求解方程中，连续性方程采用了基于黎曼求解器的连续性密度法，粒子材料的总应力张量 σ αβ
可表示为：
)
(
αβ
αβ
σ = (1− D)τ − p e +Π ij + F ij δ αβ (4)
®
(1− D) p p＜0
p e = (5)
p p≥0
τ αβ δ αβ p 为压力，受压为正。
式中： D 为损伤程度； 为剪切应力； 为恒等张量；
Monaghan 型人工黏度项 Π i j [19] 可表示为：

ρ i ρ j ( )
 −α Π c ij ν ij +β Π ν 2 ij ν ij · x ij ＜0
Π ij = 2ρ ij (6)

0 ν ij · x ij ≥0
ρ ν i j 为粒子 i 相对粒子 x ij 为粒子 i 和粒子
式中： 为平均密度， j 的速度， j 的相对位移， c ij 为平均声速， ν ij
i j
α Π ≈ 1 β Π = 1 ∼ 2 。
为相对速度 ν i j 的标量， α Π 为线性系数， β Π 为二次系数。在固体撞击中，常取 ，
1.1.2 人工应力法
拉伸不稳定性是 SPH 方法中的相关问题，它是指与应力状态相关的不稳定性，包括压缩应力下的不
稳定性。在固体大变形仿真中出现的拉伸不稳定性，会导致结果错误或计算失败。人工应力法是最简
043301-3