Page 9 - 《爆炸与冲击》2026年第2期

第 46 卷 寿列枫，等： 大尺度复杂环境下的强爆炸冲击波传播数值模拟技术研究 第 2 期

dU
= S(U) (19)
dt
由于爆炸冲击波传播速度极高，介质与外界的热传导、辐射等能量交换时间远大于冲击波的作用时
µ → ∞ 。不失一般性，以 种介质
间，一般假设爆炸冲击波传播处于绝热过程，可认为压力松弛系数满足 2
情形为例，式 (19) 可进一步简化为非线性代数方程：
(0)
(0)
(0)
∗
e (p ,v )−e i (p ,v )+ p (v −v ) = 0 (20)
∗
∗
∗
i i i I i i
p = p ∗ e 和e i 分别表示利用压力、密度和第
∗
式中：上标“0”和“*”分别表示刚性松弛前、后的值，取 I ， i
i 相的状态方程计算得到的比内能。
∑
为了使式 (20) 闭合，基于各相的体积分数应满足的饱和条件 α i = 1 ，可得
i
∑
∗ ∗
(α i ρ i ) v = 1 (21)
i
i
∗
(α i ρ i ) = (α i ρ i ) (0) ，式 (20) 可以转换成以压力 p 为单一变量的非线性代数
由于在刚性松弛阶段满足
方程：
∑
(0) ∗
(α i ρ i ) v (p) = 1 (22)
i
i
v (p) 由相应的状态方程来获得。
∗
式中： i
对式 (22) 进行迭代求解即可获得刚性松弛后的平衡压力 p，以及每一相的体积分数 α ∗ i 和密度 ρ ∗ i ，定
义函数：
∑
f(α 1 ) = (α i ρ i ) v (p)−1 (23)
(0) ∗
i
i
采用 Newton-Raphason 方法进行迭代求解：
N
(0)
(n)
f(p ) ∑ (α i ρ i ) ∂ρ i
(n)
(n)
′
p (n+1) = p − , f (p ) = − (24)
f (p ) (ρ i ) ∂p
(n) 2
(n)
′
i=1
式 中 ： ∂ρ i /∂p 为 Hugoniot 曲 线 上 各 组 分 密 度 对 压 力 的 偏 导 数 ， 即 Hugoniot 曲 线 斜 率 的 倒 数 。 对 于
h(ρ) 为密度的函数（状态方程中的非冷凝项），Hugoniot 曲线的
p(ρ,e) = Γρe+h(ρ) 这一类特殊的状态方程，
斜率可表示为：
(0) ′ (n) (0) (0) (n) 0
∂p 2ρ i h (ρ i )+2Γρ i e i +Γ(p + p )
= (25)
(n)
(0)
∂ρ H 2ρ i −Γ(ρ i −ρ i )
(0)
p ∗ α ∗ ρ ∗ 。最后利用这些
当迭代得到收敛解 后，每一相的体积分数和密度也会得到相应的收敛值 i 和 i
收敛值来修正每一组分的守恒量，即可完成整个时间步的数值求解。
2.2 网格自适应和并行计算技术
2.2.1 网格自适应技术
为提高数值计算的整体效率，网格自适应技术如今被广泛地应用于数值模拟中。本文采用 h-网格
自适应方法，利用局部加密或稀疏网格的方法以提高数值模拟的效率。基于一个后验误差估计子作为
网格自适应的指示子，指示子的拟一致分布性质使得误差在所有自由度上几乎是等分布的，最优收敛解
能够获得。特征型指示子针在实际问题中的特定区域，用于识别数值解自身的变化和间断，并利用相应
的特征识别量试探获得。
作为数值计算前剖分的背景网格一般都是共形的，即任意 2 个拥有公共边（棱或面）的单元都共享
一整个边（棱或面），但是在网格自适应过程中，随着网格的加密和稀疏，相邻的 2 个单元都有可能出现
不同的网格层级，导致在单元边界上出现悬点。在共形网格的基础上，需进一步引入半正则化网格的概
念（如图 2(a) 所示），其基本要求为：
021001-6