Page 130 - 《爆炸与冲击》2025年第12期

第 45 卷 郭景琪，等： 马赫反射波系在平面重/轻界面的入射加载 第 12 期

了提高求解精度，该求解器还引入了自适应网格加密算法，可以实现在激波与界面附近的动态网格加

密。blastFoam 求解器已经在单相双马赫反射、含颗粒爆炸等场景得到验证和应用 [26-27] 。此外针对扰动
激波冲击界面问题，Li 等 [28] 在之前的工作中已经成功使用该求解器完成了不同工况下界面 RM 失稳的
数值模拟研究。该求解器所采用的控制方程为非定常、可压缩、考虑黏性的 Navier-Stokes (NS) 方程：
∂ρ ∂(ρu i )


 + = 0

 ∂t
 ∂x i
 ( )

 ∂(ρu i )

 ∂ ρu i u j + pδ ij ∂σ ij
 + =

 ∂t ∂x j ∂x j
[ ] ( ) (1)
∂(ρe t ) ∂ (ρe t + p)u j ∂q j ∂ u i σ ij

 + = − +
 ∂t ∂x j ∂x j ∂x j



 ( )

∂(ρY s ) ∂ ρY s u j

 + = 0

∂t ∂x j
式 中 ： ρ 为 密 度 ； t 为 时 间 ； u 为 i 流 体 在 x 和 y 方 向 上 的 速 度 ， i=1,2； x 为 i 空 间 坐 标 ； p 为 压 力 ； δ 为
j
i
Kronecker 符号，i=j 时 δ = 1，i ≠ j 时 δ = 0；e 为流体的比能量；q =−k(∂T/∂x ) 为 j 方向上的热通量，k 为热
t
j
j
ij
ij
传导系数，T 为温度；Y 为组分 s 的质量分数；σ 为黏性应力张量。σ 的表达式为：
j
j
s
i
i
Å ã
∂u i ∂u j 2 ∂u i
σ i j = µ + − δ ij (2)
∂x j ∂x i 3 ∂x i
式中：μ 为流体的动力黏性系数。由 Sutherland 法则 [29] 可知，μ 与温度 T 相关：
Å ã 3
T 2 T s +T 0
µ = µ 0 (3)
T s T +T 0
−5
式中：T =273.15 K；μ 为 0 T 下的流体动力黏性系数，μ (N )=1.76 × 10 kg/(m·s)，μ (SF )=3.15 × 10 kg/(m·s)；
−5
0
6
0
0
0
2
T 为 s Sutherland 常数，T (N )=120.01 K，T (SF )=150.00 K。
s
s
6
2
此外，本研究采用较低马赫数的激波进行模拟，因此不必考虑波后流场的非完全气体效应。数值模
拟中所有气体均为理想气体，相应的气体状态方程为：
R
p = ρ T (4)
W
式中：R 为理想气体常数，W 为摩尔质量。
本文数值模拟计算域如图 1 所示，采用绕射刚体圆柱的方式生成马赫反射波系。流场整体尺寸为
L =150 mm，L =100 mm，上下边界为固壁面，左右边界为开口。激波初始位置为 x =5 mm，激波马赫数
y
x
0
M=1.8。刚体圆柱圆心位置为 x =15 mm（y 方向位置为 50 mm），直径 d=6 mm。界面初始位置为 x =
1 2
55 mm，界面为无扰动的平面 SF /N 界面。激
2
6
波前流场的初始温度为 293.15 K，初始压力为 L x
101 325 Pa。重、轻气体的初始参数如表 1 所示，
Incident Solid wall y
表中 γ 为比热容比，a 为声速。 x 0 shock wave
O x
网格划分方法采用结构网格，选择 0.500、 Interface
0.250、0.125 和 0.100 mm 等 4 种不同尺寸的初始 Cylinder
网格对研究问题进计算。提取 4 种网格尺寸下 L y Inlet d SF 6 N 2 Outlet
t=150 μs 时沿 x=40.5 mm 方向的密度分布，绘制
x 1
密度曲线如图 2 所示。结果显示，随着网格尺寸
x 2
的减小，密度曲线逐渐收敛。初始网格尺寸为 Solid wall
0.500 和 0.250 mm 的情况下，曲线两侧出现振
荡 ， 计 算 结 果 的 收 敛 性 较 差 ； 初 始 网 格 尺 寸 为 图 1 计算域示意图
0.125 和 0.100 mm 的情况下，计算结果收敛性较 Fig. 1 Schematic diagram of the computational domain
123201-3