李诗尧 于明

1)(中国工程物理研究院研究生院,北京 100088)2)(北京应用物理与计算数学研究所,北京 100094)(2017年11月21日收到;2018年8月20日收到修改稿)

基于固体炸药爆轰过程中化学反应混合区内的固相反应物与气相生成物处于力学平衡状态及热学非平衡状态的事实,提出一种考虑热学非平衡效应的反应流动模型来描述固体炸药的爆轰流动现象.该爆轰流动模型的主要特点是,在反应混合物Euler方程和固相反应物质量守恒方程的基础上,通过附加一套关于固相反应物的组分物理量的流动控制方程来表达固相反应物与气相生成物之间的热学非平衡效应.根据反应混合区内固相反应物与气相生成物这两种化学组分保持各自内能守恒的混合规则,并借助它们具有压力相等的性质以及满足体积分数总和为1的条件,推导获得的附加方程有:固相反应物的内能演化方程、体积分数演化方程及反应混合物的压力演化方程.这样,建立的爆轰模型包括:反应混合物的质量守恒方程、动量守恒方程、总能量守恒方程、压力演化方程,以及固相反应物的质量守恒方程、内能演化方程、体积分数演化方程.对所获得的爆轰模型方程组采用一个时空二阶精度的有限体积法进行数值求解,典型爆轰问题算例结果表明本文提出的固体炸药爆轰模型是合理的.

1 引 言

固体炸药的爆轰是剧烈的化学反应与复杂的激波流体力学高度耦合的非线性动力学现象,其宏观理论模型从粗略到精密不断深入发展着.早在19世纪末及20世纪初,英国工程师Chapman和法国力学家Jouguet各自独立地提出一个简化的理论模型,认为爆轰过程中化学反应是在一个无限薄的间断面上瞬时完成的,化学反应过程无须加以考虑,其作用如同一个外加的能源反映到流体力学方程中来,流体力学的质量、动量和能量守恒方程通过跨过爆轰波阵面的间断面建立,这就是Chapman-Jouguet(CJ)模型[1,2].后来在20世纪40年代至50年代,俄国物理学家Zeldovich、美国数学家von Neumann以及德国力学家Döring分别建立了考虑化学反应过程的理论模型,认为爆轰化学反应具有有限反应速率,爆轰波由惰性的前导激波以及紧随其后的具有一定宽度的化学反应区构成,化学反应区内物质看作是由固相反应物与气相生成物这两种化学组分组成的反应混合物,并且假设固相反应物与气相生成物处于局域热力学平衡状态,即固相反应物与气相生成物具有相同的速度、压力以及温度,这就是Zeldovich-von Neumann-Döring(ZND)模型[3,4].直到目前,CJ模型和ZND模型仍是固体炸药工程中应用最为广泛的两个爆轰理论模型[5,6].

CJ模型的缺点是没有考虑化学反应过程;ZND模型尽管考虑了化学反应过程,但它使用的固体炸药爆轰化学反应混合区内的固相反应物与气相生成物处于局域热力学平衡状态的假设被认为是不合理的.众所周知,ZND模型本身来源于气体爆轰动力学学科,在气体爆轰过程中反应物与生成物均为气体物质,由于气体分子的快速而密集的碰撞运动使得气体爆轰化学反应混合区内的反应物与生成物能够在极其短暂的时间内达到热力学平衡状态,因此气体爆轰过程中假设反应物与生成物处于局域热力学平衡状态是合理的[7,8].然而,固体炸药爆轰现象具有与气体爆轰现象截然不同的动力学图像.固体炸药爆轰过程中,未反应固体炸药在激波冲击作用下起爆时,首先是炸药晶体形成高温高压的“热点”进而分解为细观尺度的颗粒,然后炸药颗粒表面快速燃烧使得固相反应物转变成为气相生成物[9−11],这样爆轰化学反应混合区包含的固相反应物与气相生成物这两种化学组分物质可以看作以接触间断形式进行相互作用,于是固相反应物与气相生成物之间能够达到力学平衡状态,即化学反应过程中固-气两相组分物质具有相同的速度和压力,同时又由于固相反应物与气相生成物的材料物性和微观特性差异极大,气体分子与固体颗粒之间往往需要经过较多次以及较长时间的碰撞作用才能达到热学平衡状态,即化学反应过程中固-气两相物质会具有各自不同的温度或内能,因此,固体炸药爆轰过程中固相反应物与气相生成物应当处于力学平衡状态及热学非平衡状态.这一状况已经得到热力学特征时间尺度估算[12,13]、细观冲击起爆数值模拟[14−16]以及反应过程实验观测[17−19]的证实.

值得指出的是,文献[13]在估算固体炸药爆轰过程的热力学特征时间尺度的基础上推导出一种考虑热学非平衡的爆轰反应流动模型,本文借鉴了这种热学非平衡思想,并且在如下三方面做了改进:1)文献[13]组分体积分数通过物质等熵混合假设得出,本文通过物质内能演化得出,更符合物理实际情况;2)文献[13]对反应混合物定义了混合状态方程来封闭反应流动控制方程组,本文将混合物压力直接作为一个自变量而无须定义混合物状态方程,这样的处理方式有更好的适应性;3)正是定义了混合物状态方程使得进入混合区的激波往往产生附加的体积压缩,为克服这一问题文献[13]在数值求解过程中进行了压力松弛和重新初始化处理,本文没有定义混合物状态方程因而无须进行压力松弛和重新初始化处理,从而减小了误差.

由此,基于固体炸药爆轰过程中化学反应混合区内的固相反应物与气相生成物处于力学平衡状态及热学非平衡状态的事实,本文提出一种考虑热学非平衡效应的反应流动模型来描述固体炸药的爆轰流动现象.建立的爆轰模型包括:化学反应混合物的质量守恒方程、动量守恒方程、总能量守恒方程、压力演化方程,以及固相反应物的质量守恒方程、内能演化方程、体积分数演化方程.典型爆轰问题算例结果表明本文提出的固体炸药爆轰模型是合理的.

2 固体炸药爆轰模型方程的推导

首先给出物理量的定义.考虑一个控制体包含有固相反应物与气相生成物两种组分物质,固相反应物的物理量标记下标“s”,气相反应物的物理量标记下标“g”,则有组分密度ρs和ρg、组分内能es和eg、组分压力ps和pg、组分速度us和ug、组分体积分数fs和fg、组分温度Ts和Tg、组分比容τs和τg,且有fs+fg=1.对该控制体内的混合物密度ρ、混合物内能e、混合物比容τ有如下定义:

如果控制体内固相反应物与气相生成物处于力学平衡状态,则对该控制体的混合物压力p及速度u有如下关系式:

不考虑各类耗散因素及外力做功情况,则在力学平衡状态下固体炸药爆轰反应流动方程组可用如下反应混合物的Euler方程与固相反应物的质量守恒方程来表达[20,21]:

(3a)式为反应混合物的质量守恒方程,(3b)式为反应混合物的动量守恒方程,(3c)式为反应混合物的总能量守恒方程,(3d)式为固相反应物的质量守恒方程.

对于(3a)—(3d)式,如果再补充反应混合区内的固相反应物与气相生成物按照等温等压(热学平衡及力学平衡)的规则进行混合,即有以下各式:

则(3a)—(3d)式满足偏微分方程组封闭条件因而能够得到定解,其中考虑了固体炸药发生化学反应过程中单位质量的固相反应物转变成为气相生成物时释放了数量为q的热能.(3a)—(3d)式连同(4)式一起就构成了前述的ZND模型.

显然,在热学非平衡状态下,固相反应物与气相生成物之间的等温条件已经不再满足,必须寻找另外的混合规则以使得(3a)—(3d)式能够满足偏微分方程组封闭条件从而得到定解.根据爆轰反应流动的物理机理,可以认为反应混合区内的固相反应物与气相生成物按照各自化学组分的内能保持守恒这一规则进行混合,以此为基础可以给出与固相反应物相关的组分物理量的流动控制方程.为此,由(3a),(3b)和(3c)式可以把反应混合物的总能量守恒方程转化为反应混合物的物质导数形式的内能方程:

将(1)式代入到(5)式,有

(6)式也可化为

(7)式两个括号中式子分别表示固相反应物和气相生成物的内能变化,再考虑到发生化学反应时固相反应物转变成为气相生成物会释放热量,以及由于固相反应物与气相生成物之间存在温度差而导致的热量扩散,可设单位质量的固相反应物向气相生成物进行了总数为q的热量交换,则可以将(7)式分裂为固相反应物和气相生成物的内能演化方程:

(8a)和(8b)式说明,尽管固相反应物和气相生成物进行了热量交换,但它们的总能量是保持不变的.

针对固相反应物流动过程,(8a)式借助混合物质量守恒方程可以化成如下形式:

(9)式中出现了固相反应物的体积分数变化这一变量.为了求出固相反应物体积分数的控制方程,需借助化学反应混合区内固相反应物与气相生成物具有压力相等的性质以及满足体积分数总和为1的条件.

为简化推导过程可不失一般性地认为,固相反应物具有ps(ρs,es)形式的状态方程,气相生成物具有pg(ρg,eg)形式的状态方程,则在固相反应物与气相生成物压力相等即

的条件下,可以得到如下物质导数形式的固相反应物的压力控制方程:

由(3d)式可得

由(3d)和(9)式可得

将(11)和(12)式代入(10)式,可得

(13)式也可以转化成物质导数形式的固相反应物的体积分数控制方程:

采用与(14)式相同的思路,可以得到物质导数形式的气相生成物的体积分数控制方程:

如果定义化学组分等熵声速:

则(14)和(15)式变成

在固相反应物与气相生成物满足体积分数总和为1的条件下,即有fs+fg=1,则由(16)和(17)式可以获得反应混合物的压力演化方程:

将(18)式代入(16)式,可以获得固相反应物的体积分数演化方程:

将(19)式代入到(9)式中,可以获得固相反应物的内能守恒方程:

综上推导过程,本文所获得的考虑热学非平衡效应的固体炸药爆轰模型的反应流动控制方程组由(3a)—(3d)式、(18)—(20)式组成,具体包括:反应混合物的质量守恒方程、动量守恒方程、总能量守恒方程、压力演化方程,以及固相反应物的质量守恒方程、内能演化方程、体积分数演化方程.

3 爆轰模型方程的数值方法与典型爆轰问题算例

3.1 爆轰模型方程的数值方法

为了便于数值求解,获得的固体炸药爆轰反应流动控制方程组(3a)—(3d),(18)—(20)式在二维直角坐标系下可以写成如下矢量-矩阵形式的偏微分方程组:

式中状态变量q=[ρ,ρu,ρv,ρE,ρsfs,ρsesfs,fs,p]T;矩阵A(q)与矩阵B(q)是具有类似表达形式的Jacobi矩阵,其中

源项s(q)=[0,0,0,0,−ρr,−ρr(q−pϕ),−ρrϕ,−ρrφ]T.这里,u,v为速度分量;

Jacobi矩阵A(q)的8个特征值和特征向量如下.

右特征向量为

相应的左特征向量为

Strang分裂方法[22]被用来求解这个带源项的非守恒型双曲律方程组(21),求解的方程组分裂为双曲律方程组求解步与常微分方程组求解步.在常微分方程组求解步中,根据固体炸药爆轰通常具有快反应过程和慢反应过程的特点[9−11],可以将其化学反应率分解成快反应率和慢反应率两部分,即有r=rF+rS,相应地(21)式的源项也分解成与快反应率和慢反应率相关的两部分,由此(21)式可以分裂为

其中,

在双曲律方程组求解步(22a)中,一个具有波传播特征的二阶精度Godunov型格式被采用[23],该格式能够很好地处理双曲律方程组中的非守恒型对流项:

这里A±=RΛ±L,L和R是Jacobi矩阵A的左右特征向量构成的矩阵,Λ±是A的非负与非正特征值构成的对角矩阵,且∆qi,j=qi,j−qi−1,j,˜α=φ(θ)L∆q,其中φ(θ)是限制函数;˜Gi,j为横向通量修正.B±有与之类似的性质.

在常微分方程组求解步(22b)中,一个具有显-隐离散特性的二阶精度Runge-Kutta格式被采用[24],显式离散被用来处理非刚性源项,隐式离散被用来处理刚性源项,该格式能够很好地处理爆轰化学反应过程中出现的快慢反应差异:

3.2 典型爆轰问题算例

对固体炸药PBX9404的爆轰流动情况进行数值考察.该炸药的主要特征参数如下(本节均采用cm-g-µs单位制):ρ0=1.842,pVN=0.563,uVN=0.347,ρVN=3.041,DCJ=0.880,pCJ=0.370,uCJ=0.229,ρCJ=2.488,q=0.102.炸药爆轰化学反应率选用Lee-Tarver模型[25]:r=I(ρs/ρ−1−a)mλx+Gλy(1−λ)zpn,其中λ=ρsfs/ρ为炸药爆轰固相反应物的质量分数,这里I,G,a,m,n,x,y,z均为与炸药相关的常数;固相反应物和气相生成物选择应用广泛的Jones-Wilkins-Lee形式状态方程[26]:

这里Ak,Bk,R1k,R2k,ωk均为与炸药相关的常数.

3.2.1 平面一维爆轰问题

考虑长度为4.0的一维炸药爆轰起爆与爆轰传播问题,左端以CJ条件起爆,在如下几种网格分点条件下:5 cells/mm,10 cells/mm,30 cells/mm,50 cells/mm,100 cells/mm,获得爆轰波传播过程中几个典型时刻反应混合物的压力、密度、速度分布以及固相反应物的质量分数分布,如图1—图3所示(文中仅展示5 cells/mm,30 cells/mm,100 cells/mm三种情况),图中每条曲线对应的时刻为t=0.4,0.8,1.2,1.6,2.0,2.4,2.8,3.2,3.6,4.0.从计算结果可以看出:

1)随着网格变细,计算获得的压力von Neumann尖峰值变大,爆轰传播速度变大,达到定常状态的时间变短,化学反应区变窄;这五种网格分点条件下的压力尖峰值分别约为0.461,0.503,0.543,0.555,0.561;爆轰传播速度分别约为0.810,0.840,0.865,0.870,0.875;定常时间分别约为1.2,0.8,0.5,0.3,0.2;化学反应区宽度分别约为0.0175,0.0115,0.0091,0.0087,0.0085;爆轰压力尖峰值、传播速度、定常时间、化学反应区宽度随网格尺度的变化如图4所示;

图1 5 cells/mm条件下爆轰波传播过程物理量变化趋势 (a)压力;(b)速度;(c)密度;(d)固相反应物质量分数Fig.1.Physical variables of detonation wave under 5 cells/mm:(a)Pressure;(b)velocity;(c)density;(d)mass fraction.

图2 30 cells/mm条件下爆轰波传播过程物理量变化趋势 (a)压力;(b)速度;(c)密度;(d)固相反应物质量分数Fig.2.Physical variables of detonation wave under 30 cells/mm:(a)Pressure;(b)velocity;(c)density;(d)mass fraction.

图3 100 cells/mm条件下爆轰波传播过程物理量变化趋势 (a)压力;(b)速度;(c)密度;(d)固相反应物质量分数Fig.3.Physical variables of detonation wave under 100 cells/mm:(a)Pressure;(b)velocity;(c)density;(d)mass fraction.

图4 不同网格尺度下爆轰波传播过程物理量变化 (a)尖峰压力;(b)传播速度;(c)定常时间;(d)反应区宽度Fig.4.Physical variables of detonation wave under different mesh scales:(a)Peak pressure;(b)propagation velocity;(c)stationary time;(d)width of reaction area.

图5 柱面散心爆轰波爆轰波传播过程物理量变化趋势 (a)压力;(b)径向速度;(c)密度;(d)固相反应物质量分数Fig.5.Physical variables of cylindrically divergent detonation wave:(a)Pressure;(b)radial velocity;(c)density;(d)mass fraction.

2)30 cells/mm分点条件下von Neumann压力值和爆轰传播速度已接近理论值.

3.2.2 柱面一维散心爆轰问题

考虑半径为4.0的圆柱状炸药爆轰起爆与爆轰传播问题,中心线以CJ条件起爆,遵守上节得出的30 cells/mm网格分点条件.计算获得爆轰波传播过程中几个典型时刻反应混合物的压力、密度、径向速度分布以及固相反应物的质量分数分布,如图5所示,图中每条曲线对应的时刻为t=0.05,0.1,0.2,0.4,0.8,1.2,1.6,2.0,2.4,2.8,3.2,3.6,4.0.从计算结果可以看出:

1)炸药爆轰波传播过程中爆速和爆压均随波面半径的增加而增加,但压力传播3.0µs后,速度和密度传播2.0µs后增加不显著;

2)散心爆轰传播过程中最大压力值约为0.527,低于平面爆轰最大压力值;反应区宽度约为0.013,比平面爆轰反应区宽50%;这是由于散心爆轰波在传播过程中的面积效应使得化学反应变慢,因而传播压力更低,化学反应区更宽.

3.2.3 柱面散心爆轰的反射问题

考察无限厚炸药表面上有相距2l长度的两线状起爆器A和B同时起爆,如图6所示,产生两列相同的柱面散心爆轰波a,这两列散心爆轰波传播到对称面Ox处产生碰撞爆轰波b,接着继续向前传播形成正规反射波c,之后进一步向前传播形成爆轰波Mach反射d.由于在相互碰撞作用点处压力相对增高,所以在对称面附近会引起传播速度增加,初始落后的Mach杆会逐渐追赶上散心入射波,最后两散心爆轰波阵面趋于一致而形成平面爆轰波.

图6 柱面散心爆轰相互作用Fig.6.Interaction of two cylindrically divergent detonation waves.

图7 柱面散心爆轰波反射压力图Fig.7.Pressure contours of two cylindrically divergent detonation waves.

采用文献[27]中的算例,两起爆点距离为5.08 cm,计算域取下半平面并设四边为固壁条件,仍采用30 cells/mm网格分点条件,得到不同时刻爆轰流场的压力等值线(如图7(a)—(f)所示).其中图7(a)为炸药起爆形成拟定常散心爆轰波,图7(b)为两散心爆轰波在对称面处发生对碰,图7(c)为两爆轰波在碰撞角度较小条件下发生的正规反射,图7(d)为爆轰波的正规反射向Mach反射过渡状态,图7(e)为两爆轰波在碰撞角度较大条件下发生的Mach反射,图7(f)显示整个爆轰波阵面趋于一致而形成平面爆轰波.该数值模拟结果与文献[27]结果符合.

当两爆轰波碰撞作用的夹角大于100◦时,出现Mach反射现象,在爆轰波相互作用点附近形成三波点A(图8),其中AB为进行入射的爆轰波阵面,AC为Mach杆,AD为反射的激波.Mach杆后的流场压力比入射爆轰波后的流场压力高约25%—30%,这种高压状态驱动爆轰波使得Mach杆运动速度比散心爆轰波运动速度快约11%—13%,因此Mach杆会逐渐赶上散心爆轰波.

图8 柱面散心爆轰波Mach反射压力图Fig.8.Mach reflected pressure for two cylindrically divergent detonation waves.

4 结 论

本文提出一种考虑热学非平衡效应的反应流动模型来描述固体炸药的爆轰流动现象,该爆轰模型通过附加固相反应物的内能演化方程、体积组分演化方程以及化学反应混合物的压力演化方程来表达固相反应物与气相生成物之间的热学非平衡效应.从固体炸药的平面一维爆轰、柱面一维散心爆轰以及柱面散心爆轰波反射等几个典型算例来看,该模型能够正确预测固体炸药爆轰的起爆及传播过程中的关键特征,如von Neumann压力、爆轰波传播速度、化学反应区宽度、爆轰波反射等,数值模拟结果证明本文提出的固体炸药爆轰的理论模型是合理的.