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

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

单、有效的拉伸不稳定性处理方法，其所用的人工应力项 F ij [20] 可表示为：

( )
F ij = f n a R i +R j (7)
ij
{
p i
−ε p i ＜0
R i = ρ i 2 (8)
0 p i ≥0
R j 分别为粒子 i 和 p i 为粒子 n a 为人
式中： f ij 为粒子间距的相关量； R i 和 j 的人工压力参数； i 的压力； ε 和
ε = 0.2 n a = 4 。
工应力项参数，通常取 ，
1.1.3 光滑核函数
在 SPH 方法中，目标粒子的物理状态需要通过光滑核函数进行计算，光滑核函数对于 SPH 方法极
2
为重要。目前，常用的光滑核函数有三次样条核函数和 Wendland C 核函数，Wendland C 核函数较三次
2
样条核函数具有更好的稳定性和计算效率 [21] ，对某一目标粒子有：
® 4
(4R+1)(1−R) R＜1
W C 2 (R,h) = α (9)
0 R≥1
W C 2 为 2 ′ ′ h 之
式中： Wendland C 核函数； R = x− x /h ，为邻居粒子到当前目标粒子距离 x− x 与光滑长度
α 为各光滑函数的归一化系数。
比；
根据式 (9)，当粒子间距离超过光滑长度 h 时，Wendland C 核函数取 0，表示当前目标粒子的球形支
2
h 。该目标粒子球形支持域内的邻居粒子通常由树形搜索法 [14] 确定，随着邻居粒子数增加，
持域半径为
目标粒子物理状态的近似计算精度增大但计算效率降低。为保证各粒子物理状态计算的近似精度一
致，在常用的可变光滑长度 SPH 中，通过调整光滑长度半径内的期望粒子数 n target 来调节光滑长度，使每
个粒子的邻居粒子数保持基本不变 [22] 。
αβ −1
δ αβ δ h ，即可得到椭球形支持域的各向异性核函
此时，通过恒等张量 将光滑长度转化为张量形式
数 [23] 。与球形支持域的光滑函数相比，椭球形支持域的核函数仍采用 n target 作为参数，每个粒子的邻居粒
子数不变，以提高大变形时核函数的近似精度。
1.1.4 粒子分布方法
传统的坐标轴粒子分布法按坐标轴方向均匀分布粒子。多数情况下，学者们只对主要撞击区域的
形变、损伤等状态感兴趣。为此，本文改进了传统的粒子分布法，实现了可对局部区域精细化建模的变
分辨率粒子分布法。
q ，粒子将从最高分辨率所在位置
变分辨率粒子分布法需指定最高分辨率 A res 和分辨率的扩散比例
以球面分布方式逐层向外扩散，每层的间距按扩散比例逐渐扩大，同时每层球面上的粒子数量逐渐减
少。对于变分辨率分布时单层球面上的粒子，其空间位置由菲波那契网格分布给出：
 √
x n = 1−z cos(2πnϕ)+ x 0
2
n

 √
 2
y n =
n
1−z sin(2πnϕ)+y 0
(10)
z n = 2n−1/N −1+z 0


 ( 2 3 )
N = floor 4πR /r
N 为该球面上的总粒子数，floor 为向
式中： n 为当前粒子的序号； x n 、 y n 、 z n 分别为当前粒子的三维位置；
√
下取整函数； R 为当前球面半径； r 为当前层间距； ϕ 为分布均匀性系数，通常取 ϕ = ( 5−1)/2 x 0 、 y 0 、 z 0
；
为最高分辨位置。
1.1.5 时间积分步长
方程积分采用标准二阶龙格-库塔法，并选取自适应时间步长。为保证显式积分的稳定性，自适应
dt CFL 需满足 CFL（Courant-Friedrichs-Lewy）限制：
时间步长
Å ã
h i
dt CFL = λ CFL min (11)
c i
043301-4