张莉，吴开腾

(1.内江师范学院 数学与信息科学学院，四川 内江 641110；2.内江师范学院 四川省高等学校数值仿真重点实验室，四川 内江 641110；3.内江师范学院 四川省数据恢复重点实验室，四川 内江 641110)

0 引言

爆轰是一类特殊的爆炸现象，主要包括起爆、传播和近区效应等基本问题。爆轰波传播是炸药爆轰过程的中间环节，是爆轰近区效应研究的基础。爆轰波传播问题一直是爆轰物理研究领域的重要研究内容，主要涉及到爆轰波阵面及反应区的运动规律、波后爆轰产物的流场分布、炸药与惰性介质或真空界面处爆轰波与介质的相互作用。由于爆轰问题的复杂性，很难用理论分析方法得到解析解，实验研究方法是目前最有效的研究手段，但往往存在实验设备要求过高、实验数据采集有限、实验准备周期较长和存在安全隐患等问题，因此采用数值模拟方法研究爆轰问题具有重要的意义。随着计算机技术的迅速发展，国内外爆轰物理领域工作者进行了大量爆炸力学的数值计算工作，取得了重要的研究成果。Mader等[1]用二维拉格朗日反应流程序模拟了炸药的拐角绕射现象，认为未分解的炸药总量主要依赖于爆轰波拐角前的曲率。Hallquist等[2]用二维有限元分析程序DYNA2D模拟了爆轰波拐角绕射现象，发现反应延迟现象与拐角处的材料有关。Bdzil等[3-5]引入冲击波动力学思想来解决爆轰波传播计算，提出了计算多维爆轰波传播的方法，即爆轰波冲击动力学(DSD)方法。DSD理论指出爆速曲率关系可以将爆轰波阵面传播问题转化为曲率相关的界面运动问题，而Level Set方法是描述运动界面使用最广泛的方法之一，适用于处理曲率相关和拓扑结构变化的界面运动问题。1996年，Aslama等[6]成功地将Level Set方法应用到爆轰波冲击动力学，发展了一种适用于任意复杂二维边界条件的爆轰波传播数值计算方法。符松等[7]将Level Set方法应用到近水面水下爆炸现象的数值模拟，取得了较好的效果。陈森华等[8]研究了贴体坐标系下非理想爆轰波阵面传播计算的Level Set方法，给出了稳定的本质无振荡(ENO)差分格式及边界处理方法，并对圆弧炸药进行了数值计算。文尚刚[9]研究了三维爆轰波传播的Level Set方法，钟敏等[10]采用贴体坐标系中非正交网格的Level Set方法，研究和发展了二维和三维爆轰波传播计算问题。柏劲松等[11]以DSD理论和Level Set方法为基础，利用程序燃烧法实现了波阵面和波后流场的耦合。Wiggins等[12]利用惠更斯原理计算了二维核爆波在Ia型超新星中的传播，并利用流体力学计算了核爆波在Ia型超新星中的传播。Clavin[13]综述了激波前沿三相点的形成及气体爆轰胞状结构的分析研究进展。Ciccarelli等[14]数值模拟了爆轰波在带孔板圆管中传播的三维结构。Avdeev等[15]对气泡爆轰波的结构进行了数值模拟。马天宝等[16]提出了一种弧长算法，计算了气相爆轰波在二维管道中的传播问题，伪弧长算法可以使网格按照一定的形式自适应移动，从而达到在强间断区域自动加密的效果，提高了网格的分辨率。

重新初始化是Level Set方法的一个重要特征，也是Level Set方法数值实现的一个关键问题[17]。重新初始化问题是将随着演化后不再保持符号距离的Level Set函数重新变回符号距离函数，使其Level Set函数的梯度既不能过大、也不能过小，从而避免数值振荡，保证数值计算精度。基于偏微分方程的隐式重新初始化方法数值实现简单，界面正则性好，是经常采用的一种重新初始化方法。但是，基于偏微分方程的隐式方法很难保证Level Set函数在重新初始化过程中界面位置不动。这意味着将Level Set方法用于爆轰波阵面计算时，容易导致爆轰波阵面位置产生偏离。这种偏离在一次计算中可能不会造成太大的影响，但是经过较长时间的演化后，误差的累积会使爆轰波阵面位置产生严重的偏离和失真，计算也就失去了意义。

在DSD理论的基础上，为了避免由于误差累积导致的爆轰波阵面偏离现象，本文提出一种固定界面Level Set方法，其基本思想是在重新初始化相邻两次迭代过程中，以固定界面为基本原则，构建保持界面不动应满足的等式条件，进而推导光滑参数满足的计算公式，主要改进在于将传统Level Set方法中固定的光滑参数修正为变化的光滑参数。

1 基于DSD理论的Level Set方程

根据爆轰冲击动力学理论，爆轰波阵面传播问题可以转化为计算界面运动速度与界面当地曲率相关的界面运动问题，计算所需要的参数由实验给出，且计算波阵面位置时不涉及具体状态方程形式，具有较强的通用性。在计算区域Ω中，假设爆轰波阵面为一强间断面，构造一个Level Set函数φ(x,t)，使得在任意时刻，爆轰波阵面Γ(t)恰是φ(x,t)的零等值面，即

Γ(t)={x∈Ω|φ(x,t)=0}.

(1)

取φ(x,0)为x点到爆轰波阵面Γ(0)的符号距离函数，当φ=0时(1)式表示爆轰波阵面坐标满足的方程，当φ>0时(1)式表示炸药未反应区，当φ<0时(1)式表示爆轰产物区和化学反应。爆轰波阵面传播时，Level Set函数φ(x,t)应当满足

(2)

假设爆速曲率关系Dn(κ)近似呈线性关系，即

Dn(κ0)=DCJ(1-ακ0)，

(3)

式中：DCJ为CJ爆速；α为与炸药有关的参数；κ0为φ=0时波阵面的中值曲率，即

(4)

式中：φx、φy、φz分别为函数φ(x,t)对x、y、z的1阶偏导数；φxx、φyy、φzz分别为函数φ(x,t)对x、y、z的2阶偏导数；φxy、φxz、φyz分别为函数φ(x,t)对x、y，x、z，以及y、z的混合2阶偏导数。

同样地，流场中其他部分(φ>0和φ<0)的Dn(κ)可类似处理。若所在点的Level Set曲面中值的曲率为κ，则将κ替换为κ0即可得到所在点的Dn(κ)。

事实上，(2)式是一个Hamilton-Jacobi类型方程，只要给定适当的初始条件和边界条件，通过求解(2)式就可以捕捉界面的运动和确定界面的位置，从而掌握爆轰波阵面的传播规律。

2 固定界面Level Set方法

Level Set函数φ(x,t)的初值满足符号距离函数，可以克服由Level Set函数梯度过大或过小所造成的数值振荡，然而在进行了几个时间步求解后，由于数值方法的内在效应，Level Set函数φ(x,t)将不再满足符号距离函数，因此需要对Level Set函数φ(x,t)重新初始化，即求解初值问题(5)式的稳态解来实现：

(5)

式中：φt为φ(x,t)对t的1阶偏导数；φ0为初始Level Set函数值；sign(φ0)为符号函数，即

(6)

ε为光滑参数。

在Level Set方法重新初始化过程中，最理想的情况是既使Level Set函数重新变成符号距离函数的同时，又保证界面在这过程中不能偏离原来位置。文献[18]指出，为了减少在界面附近产生的数值振荡，光滑参数一般取ε=Δx(Δx为空间步长)，但这种取法并不能避免界面位置的偏离。

以一维重新初始化问题为例，即(5)式退化为

(7)

(8)

对(8)式进行一系列恒等变形后发现，光滑参数ε与固定界面位置之间存在数量关系[19]，即

(9)

显然，有

(10)

因此，有

(11)

3 控制方程

爆轰产物流场内动力学方程组为

(12)

p=f(ρ,e).

(13)

4 爆轰波阵面与爆轰产物流场的耦合计算

从(2)式中容易看出，方程不包含其他流体动力学参数，仅仅包含1个未知量φ，因此一种简便可行的方法是，将(2)式独立在只包含整个炸药的区域内求解。波后流场的流体动力学方程(12)式和(13)式中，如果反应进程变量λ是关于时间和坐标的已知函数，则波后流场的流体动力学方程组也可以独立求解，可以根据包括惰性介质的整个流场特征来选取坐标系及其差分格式。这两部分独立求解结果的耦合问题，可以通过程序燃烧法[11]来实现。程序燃烧法是根据已知的波阵面位置构造1个人工反应速率代替真实反应速率。为了避免对时间的积分，以化学反应分数作为反应进程变量，其数值由以下反应速率函数给出：

(14)

式中：d为所在点到波阵面的符号距离，d<0表示点在反应区内，d=0表示点在爆轰波阵面上，d>0表示点在未反应区炸药内；dc>0为人工反应速率确定的反应区宽度(一般为经验值)；对于给定的x、x1、x2，函数Λ满足性质

(15)

因此，利用(14)式便可计算反应进程量λ.

5 数值算例

空间采用5阶加权无本质振荡格式离散，时间采用3阶全变分下降龙格- 库塔格式离散，求解(2)式和(12)式。为了叙述方便，本节将传统的Level Set方法简称为LS方法，将固定界面Level Set方法简称为FLS方法。

5.1 剪切速度场

剪切速度场是测试流体介质发生拓扑变形时运动界面捕捉的经典算例。考虑剪切速度分布为

(16)

图1 t=2.0 s时刻剪切流场界面位置Fig.1 Interface location of shear flow field for t=2.0 s

计算区域为[0,1]×[0,1]，初始界面为圆心在(0.5,0.3)处、半径为0.2的圆周，计算网格为100×100. 若流场在t=2.0 s时刻后反转剪切速度，则理论上在t=4.0 s时刻流体界面应返回到初始时刻界面状态。从图1中可以看出，FLS方法能很好地捕捉到狭窄、细长的流带，而LS方法无法捕捉到狭长、细长的流带，在流体头部和尾部均出现严重的体积损失。从图2中可以看出，采用LS方法得到的返回界面已经严重失真，出现了较严重的体积损失，而采用FLS方法得到的返回界面非常接近初始圆周。这是因为剪切速度场在整个流场分布固定，LS方法和FLS方法都没有反作用于速度场，在界面变形过程中，FLS方法在重新初始化过程中能够很好地固定界面位置，使得在其剪切流场中边界能够很好地返回初始状态。

图2 t=4.0 s时刻剪切流场界面位置Fig.2 Interface location of shear flow field for t=4.0 s

若记剪切运动前界面所围的面积为S0、剪切运动后界面所围的面积为S，则定义运动前、后界面面积的比值为err，即err=S/S0，表1给出LS方法和FLS方法计算的剪切前、后界面所围面积的对比结果。从表1中可以看出，随着网格数的增加，传统LS方法和FLS方法的质量守恒性都会提高，但在同等网格条件时，FLS方法质量守恒性的优势更明显。

表1 剪切流场界面面积比值对比结果

5.2 中心一点起爆

假定爆轰波在炸药中是理想传播的，即以不变的理想爆速DCJ在炸药中进行稳定传播。正方形炸药的边长为16 mm，起爆点在炸药的中心，坐标为(8 mm,8 mm,8 mm)，设炸药的理想爆速为8.0 km/s.

图3～图5分别给出了LS方法和FLS方法在不同时刻爆轰波的传播情况。从图3～图5中可以看出，爆轰波以球形向四周发散，与文献[9]计算结果是一致的，表明LS方法和FLS方法计算的正确性。在此基础上，图6分别给出了两种方法在不同时刻的爆轰波阵面，其中截面为y=8 mm的Oxz平面。从图6中也可以看出爆轰波阵面的发展演化过程。根据中心一点起爆的爆速曲率关系，如果忽略爆速曲率的高阶项，则爆轰波在没有达到炸药边界的情况下是能够得到解析式[20](理想爆轰的CJ值)的，这可以作为检验LS方法和FLS方法计算精度的参考依据。为了简化，考虑一维情形，设初始时刻的爆轰波阵面半径为r0，根据爆轰曲率关系，有

(17)

式中：r为爆轰波阵面半径；r(0)为初始时刻爆轰波阵面半径；β为与反应区长度和外壳材料性质相关的系数。

图3 t=0.35 μs时爆轰波的计算结果Fig.3 Calculated results of detonation wave for t=0.35 μs

图5 t=1.05 μs时爆轰波的计算结果Fig.5 Calculated results of detonation wave for t=1.05 μs

图6 在不同时刻爆轰波阵面传播的计算结果Fig.6 Calculated results of detonation wavefront propagation at different times

求解(17)式可得到爆轰波阵面半径r(t)的解析式：

(18)

图7 中心一点起爆爆轰波阵面半径的对比结果Fig.7 Comparison results of radii of detonation wavefront

图7给出了两种方法计算的爆轰波阵面半径与CJ理论值的对比情况。从图7中可以看出，LS方法与CJ理论值的差距较大，FLS方法的计算结果比较接近于CJ理论值。随着时间的演化，LS方法的计算结果与CJ理论值的偏离程度越来越大，而FLS方法的计算结果始终保持接近于CJ理论值。

如表2所示，将一个计算时间(t=1.4 μs)分成4个相同的计算周期(Δt=0.35 μs)，观察在每一段计算周期内LS方法和FLS方法计算的爆轰波阵面半径与CJ理论值的绝对误差之和变化情况。从表2中可以看出，随着时间的演化，LS方法的绝对误差之和越来越大，而FLS方法的绝对误差之和变化不大，与图7的结论是一致的。这是因为LS方法不重视重新初始化过程中的界面偏离，误差的累积导致爆轰波阵面产生了偏离，而FLS方法在重新初始化过程中，通过固定界面有效地抑制了爆轰波阵面偏离。

表2 爆轰波阵面半径绝对误差之和对比结果

5.3 圆弧炸药平面起爆

设圆弧片材料为JB9014炸药，其内径、外径分别为70 mm和100 mm，圆心角为直角，主装药内、外界面均为5 mm厚的钢外壳封闭。爆轰曲率关系见表3，空间步长为0.5 mm，时间步长为0.01 μs，dc=4 mm，流体动力学计算取穿透边界条件，其他部分取刚性壁边界条件。

表3 炸药的爆速曲率关系

图8 爆轰波阵面的对比结果Fig.8 Comparison results of detonation wavefront

图8给出了LS方法、FLS方法计算的爆轰波阵面与实验结果[11]的对比情况。从图8中可以看出，由于约束条件处理的不同，LS方法和FLS方法的计算结果比实验结果[11]总体上要小些，LS方法计算结果与实验结果的平均相对误差大约为9.625 7%，FLS方法计算结果与实验结果的平均相对误差大约为2.792 3%. 因此，FLS方法比LS方法的计算结果与实验结果更加吻合。

6 结论

本文针对误差累积导致的爆轰波阵面偏离问题，根据DSD理论，在传统Level Set方法的基础上，提出了一种固定界面Level Set方法，其基本思想是在重新初始化相邻两次迭代过程中，以固定界面为基本原则，构建保持界面不动应满足的等式条件，进而推导出光滑参数满足的计算公式，主要改进在于将传统Level Set方法中固定的光滑参数修正为变化的光滑参数。

为了验证固定界面Level Set方法在描述爆轰波阵面传播计算的有效性和适用性，主要讨论了剪切速度场、中心一点起爆和圆弧炸药平面起爆的验证模型。剪切速度场数值模拟结果表明，固定界面Level Set方法在重新初始化过程中能够很好地固定流场界面位置，能够较为准确地捕捉流体介质发生拓扑变形时的运动界面，使其流场边界能够很好地返回初始状态。中心一点起爆和圆弧炸药平面起爆数值模拟结果表明，与传统Level Set方法相比，固定界面Level Set方法能够有效地避免由于误差累积导致的爆轰波阵面偏离现象，使其计算结果更接近于CJ理论值和实验结果，能够更准确地捕捉到合理的爆轰波阵面。