Page 105 - 《爆炸与冲击》2026年第01期

第 46 卷 陈 丁，等： 非药式水下爆炸冲击波加载的PD-SPH建模与分析 第 1 期

 ∑ ∑ ∑ ∑ −1
V b W ab V b W ab x ba V b W ab y ba V b W ab z ba
 
b∈Ω a b∈Ω a b∈Ω a b∈Ω a
     
∑ ∑ ∑ ∑
h ab   W ab
 V b ∇ x W ab V b ∇ x W ab x ba V b ∇ x W ab y ba V b ∇ x W ab z ba 
     
 g x,ab   b∈Ω a b∈Ω a b∈Ω a b∈Ω a   ∇ x W ab 
  ≈  ∑ ∑ ∑ ∑    (16)
 g y,ab     ∇ y W ab 
V b ∇ y W ab V b ∇ y W ab x ba V b ∇ y W ab y ba V b ∇ y W ab z ba 

 
g z,ab  b∈Ω a b∈Ω a b∈Ω a b∈Ω a  ∇ z W ab
∑ ∑ ∑ ∑
 
V b ∇ z W ab V b ∇ z W ab x ba V b ∇ z W ab y ba V b ∇ z W ab z ba
b∈Ω a b∈Ω a b∈Ω a b∈Ω a
x = [x,y,z] T 分别为三维坐标的 3 z ba 分别为粒子所对应的方向位置差，核函数
式中： 个方向， x ba 、 y ba 和
s s s 为核函数梯度。注
W ab = W(x a − x b ,l ) 为带有光滑长度 l 的预定义的核函数， ∇ x W ab = ∂W(x− x b ,l )/∂x| x=x a
a a a
意到，辅助核近似函数不要求整个计算域的光滑长度恒定。
2 PD-SPH 耦合策略耦合

引入虚粒子耦合的方案实现 PD-SPH 的流
固耦合算法 [20, 34] 。将 PD 粒子视为 SPH 的边界 SPH particle (a, b, c) PD particle (A, B, C)
虚粒子，参与 SPH 的控制方程计算。如图 1 所
Ω A
示 ， 一 共 有 3 类 支 持 域 （ 阴 影 区 域 ） ， 分 别 为
S P H 粒 子 a 支 持 域 Ω a ， P D 粒 子 A 支 持 域 Ω A ，

A|A′
A ′ Ω A ′ 。PD A 的流固耦合力
虚粒子 支持域 粒子
A ′ A ′
与虚粒子 受到的流固耦合力一致，虚粒子
Ω A′
的除应力外的物理量由 PD 粒子 A 提供。根据
式 (2) 计算得到 PD 粒子 A 的运动。根据 SPH 控 A|A′
制方程的内力计算公式 (12)，得到 SPH 粒子 a
a
与 PD 粒子 A 或其虚粒子 A ′ 之间的流固耦合力公
Ω a
式为：
∑ ∑
fsi
fsi
f = f ′ = V a V A ′ p A ′ g A ′ a − V a V A ′ p a g aA ′ (17)
A A
A ′ ∈Ω a a∈Ω A ′
图 1 PD-SPH 耦合模型示意图
其中 Fig. 1 Schematic diagram of PD-SPH coupled model
∑
p A ′ = V c p c h A ′ c (18)
c∈Ω A ′
σ A ′ 为虚粒子参与 计算所需要的插值应力，仅包含 粒子的插值。
式中： SPH SPH
3 GPU 并行框架平台介绍
根据 SPH-PD 模型计算公式，采用预测-校正时间积分方法 [31,34] ，多 GPU 计算的时间积分方案及
SPH-PD 计 算 流 程 如 图 2 所 示 。 在 初 始 化 步 骤 中 ， 根 据 用 户 自 定 义 设 置 初 始 边 界 条 件 ， 并 预 先 定 义
GPU 数量和区域分解方式。区域分解通过欧拉格式进行划分，当粒子从一个区域运动到另一个区域时，
则将粒子物理量进行信息交换。该过程完成后，将每个 GPU 分配到对应的子区域，并根据初始步骤建
立 GPU 间数据通信机制。因 SPH/PD 均为非局部计算思想，相互作用力仅与邻近粒子相关。只需要将
靠近邻近线程区域的粒子信息进行通信，采用数据包形式进行传递。计算耗时主要集中于图中加粗实
线红框内的相互作用计算部分。通过子区域分解的形式，将主要的计算量分成若干部分，每个部分分别
在单独 GPU 上进行计算，具有较高的并行效率。

011107-5