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

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

通过积分形式构造平衡方程。PD 粒子 A 的控制方程为：
w
d˙ u A ′ ′ ′
ρ A = {T (x A ,t)⟨x − x A ⟩−T (x ,t)⟨x A − x ⟩}dV x ′ + b(x A ,t) (1)
dt H A
ρ A 分别为 粒子 H A ∈ {0＜|ξ|＜δ} 为 粒子
式中： x A 、 ˙ u A 和 PD A 的坐标、速度和密度， b(x A ,t) 为体力项， PD A
⟨ ⟩
′ ′ ′ ′
的近场范围， ξ = x − x A 为相对位置， T (x A ,t) 和 T (x ,t) 分别为坐标 x A 和 x 的力矢量状态， T (x A ,t) x − x A
为粒子 A 的键层次力矢量状态中邻近粒子相互作用力。
采用非常规态型近场力学（non-ordinary state-based peridynamics, NOSB-PD）方法，引入经典连续介质
力学中的本构模型，其粒子离散形式如下：

∑ (
d˙ u A ) fsi
= T ⟨ξ AB ⟩−T ⟨ξ BA ⟩ V A V B +V A b A + f (2)
m A A B A
dt
B∈Ω A
其中：
−1
T ⟨ξ AB ⟩ = ω⟨ξ AB ⟩ P A K ξ AB (3)
A
A
P A 分别为 粒子 的质量、体积和第一类 fsi 为流固耦合作用
式中： m A 、 V A 和 PD A Piola-Kirchhoff 应力， f A
ξ AB = x A − x B 为初始构型下的 PD 粒子 A 和 B K A 的具体离散形式为：
力， 的相对位置， K A 为形状张量。
∑Ä ä
K A = ω⟨ξ AB ⟩ξ AB ⊗ξ AB V B (4)
B∈Ω A
在 NOSB-PD 中，变形梯度张量是通过对其领域内所有键的变形进行加权平均得到的。有必要引入
额外的健力，即增加粒子之间的内部约束，以抑制数值不稳定性的发生。引入人工刚度形式，修正后的
˜ σ A 如下：
应力
˜ σ A = σ A +Π AB I (5)
2
Π AB = −γ Ⅱ ρ c ϕ AB (6)
AB AB
η AB ·ξ AB
ϕ AB = 2 (7)
ξ AB ·ξ AB +(0.1δ AB )
ρ = (ρ A +ρ B )/2 c AB =
式中： I 为二阶单位张量， γ Ⅱ 为人工刚度系数， η AB = u A −u B 为当前构型的位移差， ，
AB
δ AB = (δ A +δ B )/2 分别为 PD 粒子 A 和 B 的平均密度、声速和近场范围。
(c A +c B )/2 和
固体考虑应变硬化、温度软化以及率效应，采用 Johnson-Cook 损伤模型 [27] ，其表达式为：
ñ ô
Å ã Å ã M
( N ) ˙ ε p T −T 0
σ = (1− D) A+ Bε p 1+C ln 1− (8)
˙ ε p,0 T m −T 0
ε p 为等效塑性
式中： A JC 为屈服应力， B JC 和 N JC 为应变强化系数， C JC 为率效应系数， M JC 为温度软化系数，
T m 分别为参考温度和熔点温度，D 为 JC 模型损
应变， ˙ ε p 为等效塑性应变率， ˙ ε p,0 为参考塑性应变率， T 0 和
伤。JC 模型损伤率 [28] 表述为：

0 ε p ≤D th
˙
D = ˙ ε p (9)
D c ε p ＞D th
ε f − D th
Å ãÅ ã D 4 Å ã
pD 3 ˙ ε p T −T 0
ε f = D 1 + D 2 e − σ 1+ 1+ D 5 (10)
˙ ε p,0 T m −T 0
D c 为临界损伤。根据文献 [29]，
式中：p 为静水压力， D 1 、 D 2 、 D 3 、 D 4 和 D 5 均为损伤参数， D th 为损伤阈值，
本文所采用的靶体材料的 JC 模型参数（包括弹性模量 E 及泊松比 ν 参数）如表 1 所示。

011107-3