Processing math: 100%
  • ISSN 1001-1455  CN 51-1148/O3
  • EI、Scopus、CA、JST、EBSCO、DOAJ收录
  • 力学类中文核心期刊
  • 中国科技核心期刊、CSCD统计源期刊

含空穴炸药硝基甲烷冲击转爆轰过程的数值模拟

肖敏 王成 杨同会

肖敏, 王成, 杨同会. 含空穴炸药硝基甲烷冲击转爆轰过程的数值模拟[J]. 爆炸与冲击. doi: 10.11883/bzycj-2024-0333
引用本文: 肖敏, 王成, 杨同会. 含空穴炸药硝基甲烷冲击转爆轰过程的数值模拟[J]. 爆炸与冲击. doi: 10.11883/bzycj-2024-0333
XIAO Min, WANG Cheng, YANG Tonghui. Numerical simulation for shock to detonation process of explosive nitromethane containing cavities[J]. Explosion And Shock Waves. doi: 10.11883/bzycj-2024-0333
Citation: XIAO Min, WANG Cheng, YANG Tonghui. Numerical simulation for shock to detonation process of explosive nitromethane containing cavities[J]. Explosion And Shock Waves. doi: 10.11883/bzycj-2024-0333

含空穴炸药硝基甲烷冲击转爆轰过程的数值模拟

doi: 10.11883/bzycj-2024-0333
基金项目: 国家自然科学基金(12102052,12221002);中国石油大学(北京)科研基金(2462023YJRC008);北京应用物理与计算数学研究所计算物理重点实验室基金(6142A05QN23005)
详细信息
    作者简介:

    肖 敏(1992- ),女,博士,讲师,xiaomin@cup.edu.cn

    通讯作者:

    王 成(1972- ),男,博士,教授,wangcheng@bit.edu.cn

  • 中图分类号: O381

Numerical simulation for shock to detonation process of explosive nitromethane containing cavities

  • 摘要: 为研究冲击波加载下含空穴液体炸药硝基甲烷的起爆过程,提出了一种基于水平集方法的欧拉多介质计算方法。采用反应欧拉方程组作为控制方程,通过水平集方法追踪化学反应混合物与空穴之间的界面。为提高计算方法的鲁棒性,在界面附近的计算单元内应用修正的虚拟流体方法,将多介质问题转化成单介质问题。对于这两种流体,均采用高阶加权本质无振荡(weighted essentially non-oscillatory, WENO)有限差分方法计算单元边界的数值通量,使得模拟结果具有高可靠性。由于Jones-Wilkins-Lee(JWL)状态方程与理想气体状态方程形式差别很大,爆轰产物的质量分数又直接影响了化学反应区内守恒变量和原始变量的相互转化过程,难以给出爆炸混合物状态方程的显式表达形式,因此发展了一种能够解决以上难题的虚拟流体状态预测方法。通过求解涉及化学反应的复杂多介质黎曼问题,获得界面两侧虚拟流体的变量状态。对不同强度冲击波加载下硝基甲烷与空穴的相互作用问题开展了数值模拟,数值结果可以说明:提出的计算方法能够捕捉到空穴压缩、塌陷、闭合以及消失后的流体动力学全过程。
  • 炸药在浇铸、压装过程中很难避免在内部形成空穴。根据热点理论,非均质炸药在冲击波作用下,空穴受压塌陷产生局部高温高压区域,当热点区域足够大或者高温持续时间足够长时,会引起炸药局部发生点火,实现爆轰[1],点火时间会明显短于点燃纯净炸药所需的时间。因此,考察炸药关于空穴的敏感性具有重要的研究意义和应用价值。

    科研人员采用数值模拟方法对不同结构冲击波-孔洞的相互作用问题进行过一些研究。Mader[2]模拟了液体炸药硝基甲烷内的空穴坍塌过程,这是一项早期研究成果,受到了分辨率的限制。Lauer等[3]分别用Tait状态方程和理想气体状态方程描述了水和空气这2种介质,通过处理欧拉方程组模拟了冲击波载荷下空穴阵列在水中的变形过程,探讨了水平空穴阵列的坍塌机制。Ozlem等[4]采用Godunov型计算方法,重点关注了冲击波与空穴相互作用形成的高压区域。Kapila等[5]延续了Ozlem等[4]的工作,引入关于反应率、未反应炸药和爆轰产物的平流方程,模拟了嵌入椭圆形空穴的HMX型固体炸药在冲击波加载下的爆轰过程。Kapahi等[6]研究了强瞬态载荷下微米尺度多孔高能材料的响应过程。Michael等[7]通过整合多组分方程以及欧拉方程,采用MUSCL-Hancock方法讨论了冲击波诱发空穴的坍塌过程。Xiang等[8]采用刚性状态方程描述双流体,通过计算Godunov型数值通量模拟了含空穴圆柱形水柱与平面冲击波的相互作用过程。Betney等[9]利用拉格朗日超曲面捕捉材料之间的界面,讨论了冲击波与铸入空穴水凝胶的相互作用。Sun等[10]利用陡度可调谐波(steepness-adjustable harmonic, SAH)技术进行重构,考察了空穴列对非均匀含能材料LX-17冲击起爆过程的影响。上述研究有助于理解空穴的坍塌机制,然而对于空穴压缩、塌陷以及消失后爆轰波传播全过程的数值模拟依然面临着困难和挑战。可压缩流体多组分问题常使用2种计算方法,一种是多相流方法[11-12],另一种是求解欧拉方程组的增广形式[13-14],欧拉方程组的增广形式包含质量守恒、动量守恒、能量守恒方程以及描述混合物某种性质的方程,比如化学组分的质量分数演化方程。Shyue[15-17]采用欧拉方程组作为控制方程,讨论了针对刚性状态方程、Van der Waals状态方程和Mie-Grüneisen状态方程的多组分模型,这是非常经典的研究成果,并没有考虑化学反应。而本文采用反应欧拉方程组作为控制方程,计算区域被看作由两部分组成,一部分是炸药与爆轰产物的混合物,另一部分是不与爆炸混合物融合的惰性成分空气。

    对比多相流方法的混合界面模型,多介质界面追踪方法将界面的厚度看成零,能够更精准地捕捉到界面的变形过程,但是需要面临不同介质间的界面难题。常见的处理方法有任意拉格朗日欧拉(arbitrary Lagrange Euler, ALE)方法[18]、流体体积(volume of fluid, VOF)方法[19-20]、界面追踪(front tracking)方法[21]以及水平集(level set, LS)方法[22-23]。由于水平集方法能够自然地描述界面的拓扑变化,并且能通过高阶方法处理水平集函数的控制方程,因此采用水平集/虚拟流体方法追踪化学反应混合物与空穴间的界面。传统虚拟流体方法应用起来很简便,并且容易推广至高维空间,但是不适用于求解强冲击波或者高速射流冲击等复杂问题。刘铁钢等[24]率先引入求解黎曼问题的思想,陆续发展了各种虚拟流体的变体方法,比如修正的虚拟流体方法[25]、真实的虚拟流体方法[26]以及实用的虚拟流体方法[27]等,更适用于处理含强激波或强间断的计算难题。在修正的虚拟流体方法中,需要在界面附近的计算单元内构造并求解多介质黎曼问题给出接触间断的近似状态,再根据预测值定义虚拟流体状态。该方法已经成功地应用于可压多介质流、气体化学反应流、水下爆炸[28-30]等多种问题的数值模拟。然而,当化学反应涉及复杂形式的状态方程,只能应用数值迭代方法实现守恒变量和原始变量的相互转化,通过求解多介质黎曼问题定义虚拟流体的状态成为一种难题。本文中,运用一种能够处理复杂多介质界面问题的欧拉计算方法,实现了对炸药内空穴压缩、塌陷以及消失后爆轰波传播全过程的数值模拟。

    考察如下形式的可压缩反应欧拉方程组:

    Ut+F(U)x+G(U)y=S(U) (1)

    其中:

    U=(ρρuρvρEρλ) F(U)=(ρuρu2+pρuv(ρE+p)uρuλ) G(U)=(ρvρuvρv2+p(ρE+p)vρvλ) S(U)=(0000ρR)

    式中:ρuvpEλ分别为密度、x方向的运动速度、y方向的运动速度、压力、比总能以及爆轰产物的质量分数,R为化学反应率函数。比总能E满足:

    E=e+12(u2+v2)+(1λ)Q (2)

    式中:e为比内能,Q为热量参数。

    空气与炸药和爆轰产物的混合物被认为是不会融合的,并且空气不参与化学反应。通过水平集方法追踪化学反应混合物与空气的界面时,计算区域被分割为2个子区域。其中水平集函数大于零的区域定义为化学反应混合物所在区域;水平集函数小于零的区域定义为空气所在区域。这2种流体的质量、动量、总能量和爆轰产物的总质量分数关于反应欧拉方程组(1)分别进行更新,其中对于水平集函数小于零的介质即空气,方程组(1)中的化学反应率函数R取为0,控制方程退化为无反应欧拉方程组。

    对于水平集函数大于零的介质即炸药与爆轰产物的混合物,化学反应率函数取为简化的点火增长模型,它是依赖压力项的单步反应模型[7]

    R=G1(1λ)c(λ+λ0)dpy (3)

    式中:G1cdyλ0均为常数。

    为封闭方程组(1),需要给出状态方程的表达形式。硝基甲烷是一种液体炸药,其内部比较均匀,微介观尺度下不存在裂纹、晶错等结构缺陷。因此,未反应炸药与爆轰产物可以用Jones-Wilkins-Lee状态方程描述[7]

    p=A(1ωρR1ρ0)exp(R1ρ0ρ)+B(1ωρR2ρ0)exp(R2ρ0ρ)+ωρeρ0 (4)

    式中:ωρ0ABR1R2为常数。用式(4)分别描述未反应炸药和爆轰产物时需要不同的状态方程参数。

    空气满足理想气体状态方程:

    p=(γ1)ρe (5)

    式中:γ为比热比。

    通过水平集方法确定化学反应混合物与空气的界面,水平集函数的控制方程为对流方程[31]

    ϕt+uϕx+vϕy=0 (6)

    式中:ϕ为水平集函数。零水平集标志着界面所在处。

    使用式(6)更新下一时刻的水平集函数后,还需要根据下式:

    ϕt+S(ϕ)(ϕx2+ϕy21)=0 (7)

    对水平集函数进行重新初始化,维持距离函数的性质,式中?、?、?、和?分别为??????S为符号函数。

    1999年,虚拟流体方法由Fedkiw等[32]首次提出,通过熵外推技术处理虚拟流体的密度。结合水平集和虚拟流体方法将多介质问题转化为单介质问题,有助于提高算法的鲁棒性,从而能够运用高精度数值方法分别处理各流体。冲击波加载下的空穴塌陷以及爆轰波传播过程是含强激波和强间断的计算难题,修正的虚拟流体方法考虑了界面附近各介质的物质属性,因此采用修正的虚拟流体方法处理多介质界面。在界面附近的计算单元内引入黎曼问题,通过求解接触间断两侧的物理量预测虚拟流体的状态。

    首先通过延拓方程,将当地流体的变量外推到虚拟流体所在计算单元:

    Vτ±NV=0 (8)

    式中:V为需要构造的虚拟流体变量;τ为人工时间步长;N为单位法方向向量,法方向定义为负水平集指向正水平集。

    在二维空间,需要基于真实流体和虚拟流体的状态沿着界面法方向构造黎曼问题。通常,各介质的状态方程是确定的,求解黎曼问题并没有特别困难。然而,由于本文中需要处理的未反应炸药和爆轰产物具有不同的状态方程参数,不得不通过一些数值迭代方法实现守恒变量和原始变量的相互转化,这使得构造黎曼解成为一种难题。JWL状态方程与理想气体状态方程形式差别很大,特别是JWL状态方程形式较复杂,HLLC近似黎曼解法器仅通过Rankine-Hugoniot关系式就可以计算出接触间断区域的物理变量,因此基于HLLC黎曼解法器发展了一种涉及化学反应率的虚拟流体变量求解方法。

    接下来将介绍HLLC黎曼解法器,假设黎曼解结构中由密度、速度、压力组成的原始变量从左至右依次用ULULURUR表示,那么如图1所示,黎曼解满足[33]

    图  1  近似黎曼解结构示意图
    Figure  1.  Schematic diagram of approximate Riemann solution structure
    ˜U(x,t)={ULx/tSLULSLx/tSURSx/tSRURx/tSR (9)

    黎曼解中3个波两侧的变量状态均满足Rankine-Hugoniot关系式:

    FL=FL+SL(ULUL) (10)
    FR=FL+S(URUL) (11)
    FR=FR+SR(URUR) (12)

    中间的波为接触间断,还满足式(13):

    {pL=pRuL=uR=S (13)

    综合式(10)~(13),可知接触间断两侧变量为:

    S=pRpL+ρLuL(SLuL)ρRuR(SRuR)ρL(SLuL)ρR(SRuR) (14)
    ρK=ρK(SKuKSKS) K=L, R (15)
    EK=EK+(SuK)[S+pKρK(SKuK)] K=L, R (16)

    其中最慢波速和最快波速能够通过式(17)~(18)近似地求得:

    SL=min(uLcL, uRcR) (17)
    SR=max(uL+cL, uR+cR) (18)

    在本文中,虚拟流体变量状态的求解过程可以分为如下4个步骤。

    (1)计算接触间断区域的变量:

    p=pL+pR2uRuL2ρL+ρR2cL+cR2u=pRpL+ρLuL(SLuL)ρRuR(SRuR)ρL(SLuL)ρR(SRuR)
    ρK=ρK(SKuKSKS)EK=EK+(SuK)[S+pKρK(SKuK)]K=L, R

    (2)令λK=λK (K=L, R),根据守恒变量ρKρKuKρKEK计算原始变量pK。对于爆炸混合物,这个过程需要联立以下6个方程,炸药与爆轰产物各自的状态方程、温压平衡条件、混合内能和混合比容关系式,得到未反应炸药比容与爆炸产物比容满足的非线性方程,再通过一些求解非线性方程的数值方法如牛顿迭代法处理该方程,计算出混合压力。

    (3)重置黎曼问题初始条件:

    ρK=ρKuK=upK=pKK=L, R

    (4)如果|pLpR|大于可允许误差,重新进入步骤1;否则跳出循环,ρK (K=L, R)puλK(K=L, R)即为预测的虚拟流体状态。

    方程组(1)的半离散格式为:

    (Ut)i,j=Fi1/2,jFi+1/2,jΔx+Gi,j1/2Gi,j+1/2Δy+Si,j (19)

    式中:Fi1/2,jFi+1/2,j为沿着x方向的半节点处数值通量,Gi,j1/2Gi,j+1/2为沿着y方向的半节点处数值通量,Si,j为节点处源项的数值近似。

    采用三阶TVD Runge-Kutta方法进行时间离散,应用五阶WENO-LF有限差分方法[34]求解数值通量Fi±1/2,jGi,j±1/2。以x方向为例,对WENO-LF有限差分方法进行简要描述。将节点上的通量Fi,j进行Lax-Friedrichs分裂:

    F+i,j=12(F(Ui,j)+αUi,j) (20)
    Fi,j=12(F(Ui,j)αUi,j) (21)

    式中:α=max(|λ(Ui,j)|)λ(Ui,j)为雅可比矩阵A=FU的特征值。假设矩阵LARA分别为矩阵A的左、右特征矩阵,令W±i,j=LAF±i,j,先将数值通量投影到局部特征空间,对W±i+1/2,j分别进行重构,然后再将其变换回物理空间ˆF±i+1/2,j=RAˆW±i+1/2,j

    将反应欧拉方程组进行特征分解,计算获得雅可比矩阵为:

    FU=[01000u2+pρ2u+p(ρu)p(ρv)p(ρE)p(ρλ)uvvu00upρu(E+pρ)E+pρ+up(ρu)up(ρv)u(1+p(ρE))up(ρλ)uλλ00u] (22)

    由此可知,求解雅可比矩阵的左右特征向量以及特征值时,不可避免地需要计算出压力关于守恒量的偏导数。在化学反应区域,由于同时存在炸药和爆轰产物2种物质,无法直接给出统一形式的状态方程。因此,设置温压平衡条件:

    {Tr=Tp=Tpr=pp=p (23)

    式中:下标r、p分别表示炸药与爆轰产物的对应变量。混合内能和混合比容关系式满足:

    {e=(1λ)er+λepτ=(1λ)τr+λτp (24)

    基于方程组(23)~(24)以及炸药与爆轰产物各自的状态方程,能够通过求解线性方程组计算出2种比容τrτp关于各守恒量的偏导数,再获得压力关于各守恒量的偏导数[35]

    硝基甲烷及其爆轰产物的状态方程和反应率函数的相关参数参考文献[7]。空穴直径为1 mm,比热比为1.4,初始密度、x方向与y方向的运动速度和压力分别为1.205 kg/m3、0 m/s、0 m/s和101.3 kPa。

    在空穴左侧设置入射冲击波,考察了4、6和8 GPa这3种压力加载下的空穴塌陷过程。计算区域左边界为入流边界条件,其余边界均为出流边界条件。计算网格尺寸为Δxy=1 mm/90。

    图210展示了4 GPa冲击压力下的计算结果,其中图23为空穴压缩以及塌陷过程中不同时刻的密度和压力分布。由于这段时间内爆轰产物的质量分数一直很小,因此不作展示。如图2所示,冲击波与空穴开始相互作用后,空穴的前表面受压缩后向右移动发生形变,空穴对冲击波产生稀疏作用,使得左侧的压力和密度下降。当前表面与后表面接近时,空穴周围的密度和压力开始上升。随着空穴被进一步压缩发生塌陷,会出现局部高温、高压、高密度区域,如图3所示,热点附近区域的压力达到6 GPa以上。

    图  2  4 GPa冲击压力下不同时刻的密度分布
    Figure  2.  Density distributions under the impact pressure of 4 GPa at different times
    图  3  4 GPa冲击压力下不同时刻的压力分布
    Figure  3.  Pressure distributions under the impact pressure of 4 GPa at different times
    图  4  4 GPa冲击压力下不同时刻的密度分布
    Figure  4.  Density distributions under the impact pressure of 4 GPa at different times
    图  5  4 GPa冲击压力下不同时刻的压力分布
    Figure  5.  Pressure distributions under the impact pressure of 4 GPa at different times
    图  6  4 GPa冲击压力下不同时刻的爆轰产物质量分数分布
    Figure  6.  Detonation product mass fraction distributions under the impact pressure of 4 GPa at different times
    图  7  4 GPa冲击压力下0、0.17和0.37 μs时的中轴线密度、压力和速度分布
    Figure  7.  Density, pressure and velocity distributions on the central axis under the impact pressure of 4 GPaat 0, 0.17, and 0.37 μs
    图  8  4 GPa冲击压力下0.74、1.35和2.20 μs时的中轴线密度、压力和速度分布
    Figure  8.  Density, pressure and velocity distributions on the central axis under the impact pressure of 4 GPa at 0.74, 1.35, and 2.20 μs
    图  9  不同网格尺寸下0.37 μs时的密度分布
    Figure  9.  Density distributions with different mesh sizes at 0.37 μs

    图46展示的是空穴塌缩后3个时刻的密度、压力和爆轰产物质量分数分布。热点附近区域的压力更高,因此化学反应速率更大,在热点附近形成新的冲击波后,向四周进行扩散。但是,新产生的冲击波强度并不足以引发爆轰。2.20 μs时冲击波波阵面处压力峰值约为5 GPa,远低于CJ爆轰压力,最高的爆轰产物质量分数也低于0.03,可以认为几乎没有发生化学反应。

    图78展示的是不同时刻中轴线处的密度、压力和x方向运动速度分布。从计算结果可以看出:起初受到稀疏波的影响,0.17 μs时的密度和压力从左至右下降。随后空穴进一步受到压缩发生塌缩,0.37 μs时在7.5~7.6 mm处形成能量汇聚的高密度、高压力区域,压力峰值约为7.5 GPa,速度峰值约为3.6 km/s,此时的压力峰值明显低于CJ爆轰压力,因此化学反应强度较低。空穴塌陷引起的冲击波往四周扩散,冲击波波阵面处的压力峰值逐渐下降,2.20 μs时压力峰值约为4.8 GPa。

    图910展示了不同网格尺寸下0.37 μs时的计算结果。显然,网格数量越大,计算区域内流场分布以及气泡闭合与热点产生过程刻画得越清晰。如图10所示,将网格尺寸为Δxy=1 mm/180的计算结果当作参考解,当网格尺寸为Δxy=1 mm/22.5时,捕捉不到准确的空穴塌陷过程。随着计算网格的加密,中轴线处的密度和压力分布曲线与参考解越吻合。特别是流场分布较平缓的区域,网格尺寸为Δxy=1 mm/90与Δxy=1 mm/180计算出的密度和压力分布曲线几乎重合。数值结果显示:Δxy=1 mm/22.5, Δxy=1 mm/45, Δxy=1 mm/90与Δxy=1 mm/180计算得到的密度峰值分别相差0.377、0.092和0.0285 g/cm3,压力峰值分别相差4.856、1.781和0.85 GPa,能够说明本文算法的有效性。

    图  10  0.37 μs时的中轴线密度和压力分布
    Figure  10.  Density and pressure distributions on the central axis at 0.37 μs

    图1117展示了6 GPa冲击压力下7个时刻的计算结果,其中图1112为空穴压缩以及塌陷过程中不同时刻的密度和压力分布。空穴被压缩发生塌陷,会出现高温、高压、高密度的热点区域,局部区域压力达到9 GPa以上。

    图  11  6 GPa冲击压力下不同时刻的密度分布
    Figure  11.  Density distributions under the impact pressure of 6 GPa at different times
    图  12  6 GPa冲击压力下不同时刻的压力分布
    Figure  12.  Pressure distributions under the impact pressure of 6 GPa at different times
    图  13  6 GPa冲击压力下不同时刻的密度分布
    Figure  13.  Density distributions under the impact pressure of 6 GPa at different times
    图  14  6 GPa压力下不同时刻的压力分布
    Figure  14.  Pressure distributions under the impact pressure of 6 GPa at different times
    图  15  6 GPa冲击压力下不同时刻的爆轰产物质量分数分布
    Figure  15.  Product mass fraction distributions under the impact pressure of 6 GPa at different times
    图  16  6 GPa冲击压力下0、0.15和0.30 μs时的中轴线密度、压力和速度分布
    Figure  16.  Density, pressure and velocity distributions on the central axis under the impact pressure of 6 GPa at 0, 0.15, and 0.3 μs
    图  17  6 GPa冲击压力下0.70、1.11、1.94 和2.49 μs时的中轴线密度、压力和速度分布
    Figure  17.  Density, pressure and velocity distributions on the central axis under the impact pressure of 6 GPa at 0.70, 1.11, 1.94, and 2.49 μs

    图1315展示的是空穴塌缩后4个时刻的密度、压力和爆轰产物质量分数分布。能够看出:在热点附近化学反应速率显著提高,形成新的冲击波向四周传播。在冲击波作用下高压区域内的炸药发生化学反应,快速生成新的爆轰产物。如图1415所示,随着时间演化,冲击波波阵面处压力上升,1.94 μs时压力峰值为10 GPa,2.49 μs时压力峰值增长到远高于CJ爆轰压力的22 GPa。高压区域爆轰产物的质量分数持续上升,2.49 μs时2.7~16.4 mm区域内的炸药几乎全部转化为爆轰产物。

    图1617展示的是不同时刻中轴线密度、压力和x方向运动速度的分布。起初随着空穴被压缩,炸药受到稀疏作用,密度以及压力从左至右下降。空穴进一步受到压缩发生塌陷,0.30 μs时在10.5~10.6 mm处形成能量汇聚的高密度、高压力区域,压力峰值约为10 GPa,速度峰值约为4.4 km/s,热点附近区域内的化学反应速率大幅提升,产生新的冲击波。冲击波压力峰值起先略微下降,然后随着化学反应的进行不断提高。1.94 μs时新产生的冲击波传播到5 mm处,该处压力值为10.5 GPa,2.49 μs时传播到3 mm处,该处压力值为21.7 GPa。

    图1824展示了8 GPa冲击压力下7个时刻的计算结果,其中图1819为空穴压缩以及塌陷过程中不同时刻的密度和压力分布。冲击波与空穴相互作用后,空穴的前表面受压缩后发生形变。随着空穴被进一步压缩,热点区域压力达到11 GPa以上。

    图  18  8 GPa冲击压力下不同时刻的密度分布
    Figure  18.  Density distributions under the impact pressure of 8 GPa at different times
    图  19  8 GPa冲击压力下不同时刻的压力分布
    Figure  19.  Pressure distributions under the impact pressure of 8 GPa at different times
    图  20  8 GPa冲击压力下不同时刻的密度分布
    Figure  20.  Density distributions under the impact pressure of 8 GPa at different times
    图  21  8 GPa冲击压力下不同时刻的压力分布
    Figure  21.  Pressure distributions under the impact pressure of 8 GPa at different times
    图  22  8 GPa冲击压力下不同时刻的爆轰产物质量分数分布
    Figure  22.  Detonation product mass fraction distributions under the impact pressure of 8 GPa at different times
    图  23  8 GPa冲击压力下0、0.14和0.27 μs时的中轴线密度、压力和速度分布
    Figure  23.  Density, pressure and velocity distributions on the central axis under the impact pressure of 8 GPa at 0, 0.14 and 0.27 μs
    图  24  8 GPa冲击压力下0.58、0.78、0.90和1.50 μs时的中轴线密度、压力和速度分布
    Figure  24.  Density, pressure and velocity distributions on the central axis under the impact pressure of 8 GPa at 0.58, 0.78, 0.90, and 1.5 μs

    图2022展示的是空穴塌缩后4个时刻的密度、压力和爆轰产物质量分数分布。空穴塌缩后,热点附近的压力和温度升高,化学反应速率迅速升高。随后,由热点引发新的冲击波向四周传播,冲击波波后高压区的炸药发生剧烈反应,如图2122所示,0.78 μs时压力峰值高达25 GPa,部分区域内的炸药已完全反应,转化为爆轰产物。0.90 μs时压力峰值升至27 GPa,在加载冲击波以及新产生冲击波的作用下,高压区的剩余炸药继续剧烈反应,转化为爆轰产物。1.50 μs时14 mm左侧区域的炸药几乎全部发生反应。

    图2324展示的是不同时刻的中轴线密度、压力和x方向运动速度分布。0.27 μs时7.50~7.67 mm处为能量汇聚的高密度、高压力区,速度峰值约为4.9 km/s。冲击波波阵面往四周传播,0.78 μs时在5.50 mm处冲击压力已经达到25 GPa,高于约14 GPa的CJ爆轰压力,随后压力峰值一直高于CJ爆轰压力。

    通过数值模拟3种不同强度冲击波加载下的液体炸药硝基甲烷与气泡相互作用过程,可以得到结论:在热点即空穴塌缩处能够实现化学反应的加速,并且化学反应的剧烈程度随着冲击波压力的增加而增大。4 GPa冲击压力加载下,热点产生后爆轰产物的质量分数增长非常缓慢。而6 GPa冲击压力加载下,炸药内空穴塌缩形成热点,能够实现爆轰,引起化学反应区域的扩散。这与Turley等的实验结果[36]显示:冲击波强度不够高时,化学反应区域会停止增长,而当冲击压力超过5 GPa时,在空穴塌缩处发生爆轰,爆轰产物区域会逐渐向四周扩大吻合。

    提出了一种处理炸药和爆轰产物的混合物以及空穴间界面问题的虚拟流体方法,并且采用基于水平集技术的多介质界面追踪方法对炸药内空穴压缩、塌陷等介观尺度下动力学过程开展了数值模拟。根据模拟结果发现:在冲击波载荷下含空穴液体炸药硝基甲烷比纯净炸药的化学反应速率显著提高。4 GPa压力加载下,由于空穴塌陷形成新冲击波的强度不足以引发爆轰,生成的爆轰产物非常少。当冲击压力升至6 GPa时,空穴塌缩、闭合处炸药的化学反应速率大幅提高,成功实现爆轰,反应区内不断产生爆轰产物。含空穴硝基甲烷在8 GPa压力加载下能更快实现爆轰。

  • 图  1  近似黎曼解结构示意图

    Figure  1.  Schematic diagram of approximate Riemann solution structure

    图  2  4 GPa冲击压力下不同时刻的密度分布

    Figure  2.  Density distributions under the impact pressure of 4 GPa at different times

    图  3  4 GPa冲击压力下不同时刻的压力分布

    Figure  3.  Pressure distributions under the impact pressure of 4 GPa at different times

    图  4  4 GPa冲击压力下不同时刻的密度分布

    Figure  4.  Density distributions under the impact pressure of 4 GPa at different times

    图  5  4 GPa冲击压力下不同时刻的压力分布

    Figure  5.  Pressure distributions under the impact pressure of 4 GPa at different times

    图  6  4 GPa冲击压力下不同时刻的爆轰产物质量分数分布

    Figure  6.  Detonation product mass fraction distributions under the impact pressure of 4 GPa at different times

    图  7  4 GPa冲击压力下0、0.17和0.37 μs时的中轴线密度、压力和速度分布

    Figure  7.  Density, pressure and velocity distributions on the central axis under the impact pressure of 4 GPaat 0, 0.17, and 0.37 μs

    图  8  4 GPa冲击压力下0.74、1.35和2.20 μs时的中轴线密度、压力和速度分布

    Figure  8.  Density, pressure and velocity distributions on the central axis under the impact pressure of 4 GPa at 0.74, 1.35, and 2.20 μs

    图  9  不同网格尺寸下0.37 μs时的密度分布

    Figure  9.  Density distributions with different mesh sizes at 0.37 μs

    图  10  0.37 μs时的中轴线密度和压力分布

    Figure  10.  Density and pressure distributions on the central axis at 0.37 μs

    图  11  6 GPa冲击压力下不同时刻的密度分布

    Figure  11.  Density distributions under the impact pressure of 6 GPa at different times

    图  12  6 GPa冲击压力下不同时刻的压力分布

    Figure  12.  Pressure distributions under the impact pressure of 6 GPa at different times

    图  13  6 GPa冲击压力下不同时刻的密度分布

    Figure  13.  Density distributions under the impact pressure of 6 GPa at different times

    图  14  6 GPa压力下不同时刻的压力分布

    Figure  14.  Pressure distributions under the impact pressure of 6 GPa at different times

    图  15  6 GPa冲击压力下不同时刻的爆轰产物质量分数分布

    Figure  15.  Product mass fraction distributions under the impact pressure of 6 GPa at different times

    图  16  6 GPa冲击压力下0、0.15和0.30 μs时的中轴线密度、压力和速度分布

    Figure  16.  Density, pressure and velocity distributions on the central axis under the impact pressure of 6 GPa at 0, 0.15, and 0.3 μs

    图  17  6 GPa冲击压力下0.70、1.11、1.94 和2.49 μs时的中轴线密度、压力和速度分布

    Figure  17.  Density, pressure and velocity distributions on the central axis under the impact pressure of 6 GPa at 0.70, 1.11, 1.94, and 2.49 μs

    图  18  8 GPa冲击压力下不同时刻的密度分布

    Figure  18.  Density distributions under the impact pressure of 8 GPa at different times

    图  19  8 GPa冲击压力下不同时刻的压力分布

    Figure  19.  Pressure distributions under the impact pressure of 8 GPa at different times

    图  20  8 GPa冲击压力下不同时刻的密度分布

    Figure  20.  Density distributions under the impact pressure of 8 GPa at different times

    图  21  8 GPa冲击压力下不同时刻的压力分布

    Figure  21.  Pressure distributions under the impact pressure of 8 GPa at different times

    图  22  8 GPa冲击压力下不同时刻的爆轰产物质量分数分布

    Figure  22.  Detonation product mass fraction distributions under the impact pressure of 8 GPa at different times

    图  23  8 GPa冲击压力下0、0.14和0.27 μs时的中轴线密度、压力和速度分布

    Figure  23.  Density, pressure and velocity distributions on the central axis under the impact pressure of 8 GPa at 0, 0.14 and 0.27 μs

    图  24  8 GPa冲击压力下0.58、0.78、0.90和1.50 μs时的中轴线密度、压力和速度分布

    Figure  24.  Density, pressure and velocity distributions on the central axis under the impact pressure of 8 GPa at 0.58, 0.78, 0.90, and 1.5 μs

  • [1] 彭亚晶, 叶玉清. 含能材料起爆过程“热点”理论研究进展 [J]. 化学通报, 2015, 78(8): 693–701. DOI: 10.14159/j.cnki.0441-3776.2015.08.004.

    PENG Y J, YE Y Q. Research progress of ‘hot-spot’ theory in energetic materials initiation [J]. Chemistry, 2015, 78(8): 693–701. DOI: 10.14159/j.cnki.0441-3776.2015.08.004.
    [2] MADER C L. Numerical modeling of detonations [M]. Berkeley: University of California Press, 1979.
    [3] LAUER E, HU X Y, HICKEL S, et al. Numerical investigation of collapsing cavity arrays [J]. Physics of Fluids, 2012, 24(5): 052104. DOI: 10.1063/1.4719142.
    [4] OZLEM M, SCHWENDEMAN D W, KAPILA A K, et al. A numerical study of shock-induced cavity collapse [J]. Shock Waves, 2012, 22(2): 89–117. DOI: 10.1007/s00193-011-0352-9.
    [5] KAPILA A K, SCHWENDEMAN D W, GAMBINO J R, et al. A numerical study of the dynamics of detonation initiated by cavity collapse [J]. Shock Waves, 2015, 25(6): 545–572. DOI: 10.1007/s00193-015-0597-9.
    [6] KAPAHI A, UDAYKUMAR H S. Dynamics of void collapse in shocked energetic materials: physics of void-void interactions [J]. Shock Waves, 2013, 23(6): 537–558. DOI: 10.1007/s00193-013-0439-6.
    [7] MICHAEL L, NIKIFORAKIS N. A hybrid formulation for the numerical simulation of condensed phase explosives [J]. Journal of Computational Physics, 2016, 316: 193–217. DOI: 10.1016/j.jcp.2016.04.017.
    [8] XIANG G M, WANG B. Numerical study of a planar shock interacting with a cylindrical water column embedded with an air cavity [J]. Journal of Fluid Mechanics, 2017, 825: 825–852. DOI: 10.1017/jfm.2017.403.
    [9] BETNEY M R, ANDERSON P A, DOYLE H, et al. Numerical and experimental study of shock-driven cavity collapse [C]//9th International Symposium on Cavitation (CAV2015). Lausanne: IOP Publishing, 2015: 012011. DOI: 10.1088/1742-6596/656/1/012011.
    [10] SUN J, YANG P F, MENG B Q, et al. Effects of micrometer-scale cavities on the shock-to-detonation transition in a heterogeneous LX-17 energetic material [J]. Physics of Fluids, 2023, 35(12): 126117. DOI: 10.1063/5.0174851.
    [11] 于明. 固体炸药爆轰与惰性介质相互作用的一种扩散界面模型 [J]. 爆炸与冲击, 2020, 40(10): 104202. DOI: 10.11883/bzycj-2019-0435.

    YU M. An improved diffuse interface model for the numerical simulation of interaction between solid explosive detonation and inert media [J]. Explosion and Shock Waves, 2020, 40(10): 104202. DOI: 10.11883/bzycj-2019-0435.
    [12] SCHMIDMAYER K, BRYNGELSON S H, COLONIUS T. An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics [J]. Journal of Computational Physics, 2020, 402: 109080. DOI: 10.1016/J.JCP.2019.109080.
    [13] XIAO M, NI G X, WANG C, et al. Front capturing by level set method for the reactive Euler equations [J]. International Journal for Numerical Methods in Fluids, 2021, 93(8): 2723–2743. DOI: 10.1002/fld.4995.
    [14] KLOPSCH R, GARAN N, BACH E, et al. Parametric influence on rotating detonation combustion: insights from fast reactive Euler simulations [J]. AIAA Journal, 2024, 62(1): 127–139. DOI: 10.2514/1.J063193.
    [15] SHYUE K M. An efficient shock-capturing algorithm for compressible multicomponent problems [J]. Journal of Computational Physics, 1998, 142(1): 208–242. DOI: 10.1006/jcph.1998.5930.
    [16] SHYUE K M. A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state [J]. Journal of Computational Physics, 1999, 156(1): 43–88. DOI: 10.1006/jcph.1999.6349.
    [17] SHYUE K M. A fluid-mixture type algorithm for compressible multicomponent flow with Mie-Grüneisen equation of state [J]. Journal of Computational Physics, 2001, 171(2): 678–707. DOI: 10.1006/jcph.2001.6801.
    [18] ZHANG F, CHENG J. A bound-preserving and positivity-preserving high-order arbitrary Lagrangian-Eulerian discontinuous Galerkin method for compressible multi-medium flows [J]. SIAM Journal on Scientific Computing, 2024, 46(3): B254–B279. DOI: 10.1137/23M1588810.
    [19] HIRT C W, NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries [J]. Journal of Computational Physics, 1981, 39(1): 201–225. DOI: 10.1016/0021-9991(81)90145-5.
    [20] MULBAH C, KANG C, MAO N, et al. A review of VOF methods for simulating bubble dynamics [J]. Progress in Nuclear Energy, 2022, 154: 104478. DOI: 10.1016/j.pnucene.2022.104478.
    [21] TRYGGVASON G, BUNNER B, ESMAEELI A, et al. A front-tracking method for the computations of multiphase flow [J]. Journal of Computational Physics, 2001, 169(2): 708–759. DOI: 10.1006/jcph.2001.6726.
    [22] GIBOU F, FEDKIW R, OSHER S. A review of level-set methods and some recent applications [J]. Journal of Computational Physics, 2018, 353: 82–109. DOI: 10.1016/j.jcp.2017.10.006.
    [23] 姚成宝, 王宏亮, 浦锡锋, 等. 空中强爆炸冲击波地面反射规律数值模拟研究 [J]. 爆炸与冲击, 2019, 39(11): 112201. DOI: 10.11883/bzycj-2018-0287.

    YAO C B, WANG H L, PU X F, et al. Numerical simulation of intense blast wave reflected on rigid ground [J]. Explosion and Shock Waves, 2019, 39(11): 112201. DOI: 10.11883/bzycj-2018-0287.
    [24] 刘铁钢, 许亮. 模拟多介质界面问题的虚拟流体方法综述 [J]. 气体物理, 2019, 4(2): 1–16. DOI: 10.19527/j.cnki.2096-1642.0746.

    LIU T G, XU L. A review of ghost fluid methods for multi-medium interface simulation [J]. Physics of Gases, 2019, 4(2): 1–16. DOI: 10.19527/j.cnki.2096-1642.0746.
    [25] LIU T G, KHOO B C, YEO K S. Ghost fluid method for strong shock impacting on material interface [J]. Journal of Computational Physics, 2003, 190(2): 651–681. DOI: 10.1016/S0021-9991(03)00301-2.
    [26] WANG C W, LIU T G, KHOO B C. A real ghost fluid method for the simulation of multimedium compressible flow [J]. SIAM Journal on Scientific Computing, 2006, 28(1): 278–302. DOI: 10.1137/030601363.
    [27] XU L, FENG C L, LIU T G. Practical techniques in ghost fluid method for compressible multi-medium flows [J]. Communications in Computational Physics, 2016, 20(3): 619–659. DOI: 10.4208/cicp.190315.290316a.
    [28] HUO Z X, JIA Z P. A GRP-based tangential effects preserving, high resolution and efficient ghost fluid method for the simulation of two-dimensional multi-medium compressible flows [J]. Computers and Fluids, 2024, 276: 106261. DOI: 10.1016/j.compfluid.2024.106261.
    [29] XU L, YANG W B, LIU T G. An interface treatment for two-material multi-species flows involving thermally perfect gases with chemical reactions [J]. Journal of Computational Physics, 2022, 448: 110707. DOI: 10.1016/j.jcp.2021.110707.
    [30] ZHAO Z T, RONG J L, ZHANG S X. A numerical study of underwater explosions based on the ghost fluid method [J]. Ocean Engineering, 2022, 247: 109796. DOI: 10.1016/j.oceaneng.2021.109796.
    [31] OSHER S, SETHIAN J A. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations [J]. Journal of Computational Physics, 1988, 79(1): 12–49. DOI: 10.1016/0021-9991(88)90002-2.
    [32] FEDKIW R P, ASLAM T, MERRIMAN B, et al. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method) [J]. Journal of Computational Physics, 1999, 152(2): 457–492. DOI: 10.1006/jcph.1999.6236.
    [33] TORO E F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction [M]. 3rd ed. Heidelberg: Springer Science & Business Media, 2009.
    [34] JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes [J]. Journal of Computational Physics, 1996, 126(1): 202–228. DOI: 10.1006/jcph.1996.0130.
    [35] WANG C, LIU X Q. High resolution numerical simulation of detonation diffraction of condensed explosives [J]. International Journal of Computational Methods, 2015, 12(2): 1550005. DOI: 10.1142/S021987621550005X.
    [36] TURLEY W D, LA LONE B M, MANCE J G, et al. Experimental observations of shock-wave-induced bubble collapse and hot-spot formation in nitromethane liquid explosive [J]. Journal of Applied Physics, 2021, 129(14): 145102. DOI: 10.1063/5.0039414.
  • 加载中
图(24)
计量
  • 文章访问数:  123
  • HTML全文浏览量:  28
  • PDF下载量:  38
  • 被引次数: 0
出版历程
  • 收稿日期:  2024-09-09
  • 修回日期:  2025-04-10
  • 网络出版日期:  2025-04-16

目录

/

返回文章
返回