0 引言
氟气具有强烈的氧化性,可用于许多需要高反应能的工业过程,氟气被逐渐应用在推进系统中,氟气与氢气相遇会发生爆炸反应,学者们对氢氟反应机理做了很多研究,然而氢氟爆炸反应难以控制,实验危险性大,导致应用研究推进缓慢,氢氟爆轰特性尚不清楚。
氢气与氟气混合就可以发生爆炸,其爆炸极限在20世纪六七十年代[1]被研究过,Levy等[2]通过流动反应器研究了氢氟的动力学行为,Orkin等[3]通过实验测得了激发后的氢气与分子氟的反应速率,Matsugi等[4]和Werner等[5]通过实验与分子动力学模拟研究了氢氟链式反应机理。Yuri[6]通过流动反应器确定了氟原子与硫化氢的反应速率系数。Zhang[7]通过数值计算了氟在氢中的流动及爆炸特性。
在计算炸药爆轰参数的方法中,应用最广的Beckker-Kistiakowsk-Wilson(BKW)状态方程已经有五六十年的发展历史,其计算程序比较适合于凝聚类炸药,但对于计算液体和气体燃料存在严重缺陷[8]。Jacobs-Cowperthwaite-Zwisler(JCZ)状态方程是基于自由体积理论的液体模型;Jones-Wilkins-Lee(JWL)状态方程同时兼顾了固体模型和液体模型中冷能冷压贡献和分子热运动贡献,然而这些经验模型含有多个可调参数,在应用上局限性大,经验参数经过反复调整仍与实际应用环境会有较大误差。
吴雄基于维里理论和相似论建立起来的VLW状态方程,其中所有参数都能从分子相互作用势参数[8]中导出,具有更强的理论基础和更广的适用范围。近年已发展为从常压到几十万大气压广阔压力范围均适用的通用状态方程,其炸药密度范围适用于从0.005 g/cm3到3.5 g/cm3,且克服了BKW状态方程的爆温计算值较实验值偏低的问题。
吴雄[9-10]利用VLW状态方程计算了CHNO、CNO、HNO、HN、CHNO及CHNF等各种类型炸药的爆轰参数,与BKW和LJD状态方程的计算值进行了比较,表明VLW状态方程计算值更符合实验值。张金龙等[11]通过VLW状态方程理论计算了RDX爆轰产物组成和爆热值。Yue[12-14]应用VLW状态方程计算了几种典型的工业炸药爆轰参数和叠氮化铜微雷管的性能。保福成等[15]利用3种状态方程计算了高通量含能分子的性能,并进行了评估。Zhao等[16]利用VLW状态方程计算了3种新型烈性炸药的爆轰性能,并与传统炸药进行了比较。
基于VLW状态方程,结合爆轰理论方程组,利用Matlab编写计算程序,计算氢氟爆轰参数;利用分子动力学模拟拟合出了氢氟反应的速率与温度的关系,用于Fluent计算软件定义氢氟反应,模拟了氢氟爆轰并进行验证,对研究实验危险系数大的化学品的爆轰性能具有一定的借鉴意义。
1 VLW爆轰产物状态方程
1.1 维里理论方程
众所周知,由统计力学理论推导出来的维里(VIRIAL)方程,完美地解决了爆轰产物状态方程问题:
(1)
式(1)中:B、C、D为第二、第三、第四维里系数;V为气态爆轰产物的摩尔体积;R为气体常数。
(2)
(3)
式(2)中,η12和dω1等符号的定义见文献[8]。该文献指出:“计算这些高级维里系数是很困难的,除了用最简单的刚球模型外,计算非常麻烦。第四维里系数尚未计算出来”。文献[8]中只给出了其大致猜想曲线。
这就是为什么理论完美无缺的维里方程,自1940年Mayer用统计力学理论推导出以来,至今没有得到实际应用的原因所在。
1.2 VLW状态方程
统计物理指出:高温下,高级维里系数与第二维里系数之间存在某种幂函数关系。寻找这个关系,将计算难度很大的高级维里系数均通过第二维里系数来计算解决。
文献[8]已经给出了第二、第三维里系数如下理论关系式:
B=b0B*(T*)
(4)
(5)
(6)
(7)
(8)
其中:
式(5)和式(7)中:B*(T*)、C*(T*)分别为无量纲第二、第三维里系数;N为亚佛加德罗常数;κ为玻尔兹曼常数; T*为无量纲温度;σ、ε为Lennard-Jones势参数,用来模拟2个电中性的分子或原子间相互作用势能的一个比较简单的数学模型,最早由数学家Lennard-Jones在1924年提出,由于其解析形式简单而被广泛使用,特别是用来描述气体分子间相互作用尤为精确。文献只是将以表格形式给出来。
由式(5)和式(8)可知:当T*>1时,这2个级数收敛很快,当T*很大时(爆轰条件下,T*>20),可以只取第1项而忽略后面各项,此时:
于是:
(9)
钱学森指出:高温下,高级维里系数随温度的变化不大,它们都趋于某一正数,其数量级各为b0的某次幂,高温时,C(T)和
同数量级幂,D(T)和
同数量级,其余类推。由高级维里系数与第二维里之间存在某种幂函数关系,式(1)可以写成:
其中,
上式收敛很快,前3项在方程中起主要作用,后面诸项都是小项,对计算结果影响甚小,取前5项已足够精确。因此得到:
(10)
式(10)中,n≥3, m≥5。
式(10)称为VLW爆轰产物状态方程,VL代表VIRIAL理论,W代表吴雄(Wu)[10]。
联立VLW状态方程和爆轰理论方程组可解除爆轰性能参数。
1.3 理论计算结果
基于VLW爆轰产物状态方程和爆轰理论方程组在Matlab中编写计算程序,解出氢氟爆轰的性能参数如表1所示。
表1 氢氟爆轰参数理论计算结果
Table 1 Hydrogen and fluorine detonation parameters
Mole of H2∶F2D/(m·s-1)p/MPaT/K2∶17 94410.225 7471.5∶18 03512.306 1611∶18 04415.316 4791∶1.57 41114.785 9161∶26 90113.695 173
从计算结果可以看出,当氢气氟气的摩尔比为1∶1时爆轰性能最好,以下的数值模拟计算均在氢氟摩尔比1∶1的情况下进行。
2 氢氟反应参数的分子动力学模拟
2.1 分子动力学方法
通过基于密度泛函理论的第一性原理计算来获取氢与氟之间的反应参数,第一性原理依托于量子力学,利用电荷、质量、普朗克常量和原子序数来求解薛定谔方程,从而获得基态电子结构与能量,这个过程里面没有使用任何经验参数。
系综的概念是在采用统计方法,描述热力学系统的统计规律时引入的,它表示相同宏观状态下,分别处于不同微观状态的相互独立,并且无相互作用的热力学体系的集合。本文中的计算采用了正则系综(NVT)。正则系综表示系统与外界不存在物质交换并且体积不变。为了保证系统温度恒定,将其置于温度恒为T的热浴中,因此,系统中的原子数N、体积V和温度T均保持不变。
为了实现正则系综,需要对系统的温度和压力进行调节,需要选择正确的温和控压法。本文中仿真所选择的温控方法是Nose-Hoover法。
为了得到氢氟反应的反应速率随时间的变化关系,通过分子动力学计算氢氟在不同温度下的反应速率,利用阿伦尼乌斯公式拟合出反应参数指前因子A和活化能Ea。
2.2 氢氟反应分子动力学模型与计算结果
在MS软件中建立H2分子与F2分子结构和10×10×10的立方体晶胞。将10个H2分子和10个F2分子随机放置在盒子中,得到的晶胞计算模型如图1所示。在COMPASS力场、298 K下进行优化,获得H2/F2分子混合物在此力场下最稳定的结构体系。将此结构作为进行高温反应计算的初始模型。
图1 氢氟反应晶胞模型结构
Fig.1 Model structure of crystal cell
采用CP/2K软件对H2/F2分子混合物结构进行能量最小化计算,获得0 K下的稳定结构,之后在10 K下进行弛豫5 ps,获得10 K下的稳定构型,最后在1 000、1 500、2 000、2 500、3 000 K高温下进行化合反应。
不同温度下的反应速率通过拟合势能曲线得到,如图2所示。势能曲线拟合参数如表2所示。
图2 不同温度下氢氟反应的势能曲线(曲线)及相应的拟合曲线(虚线)
Fig.2 Potential energy curves(full line) and corresponding fitting curves(dotted line) at different temperatures
表2 不同温度下氢氟反应势能曲线拟合参数
Table 2 Potential energy curve fitting parameters
T/KUmax/(a.u.)tmax/psU∞/(a.u.)ΔU/(a.u.)k/ps-11 000-493.4461.407-494.8911.444 22.0421 500-493.4010.254-494.8431.442 31.0972 000-493.3480.290-495.0111.662 91.2162 500-493.4140.013-495.2991.884 60.6753 000-493.3690.243-495.5642.194 50.742
U=U∞+ΔUexp[-k(t-tmax)]
(11)
式(11)中:U∞为产物的势能;ΔU=Umax-U∞,Umax为势能的最大值; tmax为势能达到最大值所用的时间;k为该温度下的反应速率。
势能曲线拟合参数如表2所示。
拟合出的指前因子A为4.36×1011 s-1,活化能Ea为-12.79 kJ/mol。
3 氢氟爆轰的FLUENT数值计算
3.1 FLUENT控制方程
质量守恒方程:在Fluent中,质量守恒方程的一般形式如式(12)所示。无论是对于可压缩流还是不可压缩流来说,这个方程都是需要满足的。
(12)
式(12)中: ρ为密度;t为时间;
为速度矢量;Sm为用户定义的源项,即仿真过程中加入到连续相中的质量,其中包括但不限于液滴的汽化。
动量守恒方程:该定律可以表述为,在给定的流体系统,其动量的时间变化率等于作用于其上的外力总和。其形式如式(13)所示。
(13)
式(13)中: p为静压;
为应力张量;
和
分别为重力体力和外部体力;为依赖于模型的源项。
能量守恒方程:
(14)
式(14)中:Keff为有效导热系数; h为可压缩气体的显焓;
为物质j′的扩散通量; Sh为化学反应能量源项; xi为空间坐标在i方向上的分量; ui为速度矢量在i方向上的分量。
标准k-ε模型:该模型是一个比较经典且十分常用的双方程涡粘模型,由湍流动能运输方程(k方程)与湍流耗散率运输方程(ε方程)组成,其湍流动能运输方程如式(15)所示,湍流耗散率方程如式(16)所示。
Gk+Gb-ρε-YM+Sk
(15)
(16)
式(15)和式(16)中: Gk为由平均速度梯度产生的湍流动能; Gb为由浮力产生的湍流动能; YM为可压缩湍流脉动膨胀对总耗散率的影响; σk为k方程的普朗特数; σε为ε方程的普朗特数; Sk与Sε为源项; μt为湍流粘性系数; C1ε、C2ε、C3ε和σε都是常数,由实验确定。
气相燃烧采用了有限速率/涡耗散模型,该模型同时对阿伦尼乌斯反应速率RA和湍流混合速率RE进行计算,并取两者中的较小值作为最终计算反应速率
R=min{RA,RE}
(17)
涡耗散模型认为反应速率由大涡混合尺度k/ε控制,只要湍流出现(k/ε>1),反应即可进行。通过涡耗散模型计算的湍流反应速率RE取燃料、氧化剂和产物三者反应速率中的最小值
RE=min{Rfu,Rox,Rp}
(18)
式(18)中: fu为燃料; ox为氧化剂; p为反应产物。
3.2 氢氟爆轰的FLUENT计算模型
计算模型如图3所示。结论分子动力学模拟出的参数对氢氟反应的详细化学机理[5]进行简化,得到106步的基元反应,在FLUENT中进行模拟计算。计算模型为φ100×1 000的管道,对摩尔比1∶1的氢氟预混气体进行2 000 K高温点火,操作压力为1个大气压。
图3 氢氟爆轰计算模型
Fig.3 The calculation model of hydrogrn/fluorine detonation
3.3 网格收敛性验证
本文中利用Ansys Mesh进行网格划分,采用四边形网格,通过控制网格单元大小来控制网格疏密程度。网格参数如表3所示。其中,element quality代表了网格质量,其值越接近1,代表网格质量越好;aspect ratio代表网格的正交比,其值越接近1越好,一般不要超过5;Skewness代表网格的倾斜度,其值越接近0越好,一般不要超过0.8。
表3 不同尺度网格参数
Table 3 Parameters of different scale grids
默认单元大小网格质量网格正交比网格倾斜度网格数量0.3 mm0.997 931.003 31.961 9×10-3556 6000.5 mm0.998 251.001 21.588 7×10-3200 0070.7 mm0.996 871.007 53.104 7×10-3101 6581.0 mm0.999 081.000 65.212 4×10-450 0001.5 mm0.998 271.010 61.381 1×10-322 011
在管壁上50 mm设置监测点,在同样的初始条件下进行计算,参数如表4所示。计算3 000个时间步长,得到温度时间曲线如图4所示。
图4 不同网格密度温度仿真结果
Fig.4 Temperature simulation results at different mesh densities
表4 计算初始参数设置
Table 4 The settings of the initial parameters
参数数值操作压力/Pa101 325表压/Pa0温度/K298壁面导热率/(W·(m2·k)-1)58.2高温点火半径/mm2点火温度/K2 000时间步长/s1×10-8
从图4中可以看出,不同的网格密度计算出的温度波动趋势和温度峰值基本相同,区别在于温度的上升沿时间和上升速率,不同网格密度计算出的参数相对偏差如表5所示。
表5 计算出的参数与相对偏差
Table 5 Calculated parameters and relative deviations
温度峰值/K相对误差/%上升沿开始时间(time steps)相对误差/%0.3 mm6 809-396-0.5 mm6 7081.474226.570.7 mm7 1364.8147219.191.0 mm7 0223.1358146.721.5 mm6 2298.5270878.79
从表5中可以看出,0.5 mm的网格相较于0.3 mm网格计算出的参数相对偏差在7%之内;0.7、1和1.5 mm网格计算出的峰值温度相较于0.3 mm的网格相对偏差小于9 %;温度上升沿偏差较大,随着网格的变大,上升沿开始时间偏差逐渐增大,网格不合适;0.3 mm的网格数约为0.5 mm网格数的2倍,占用过多计算计算资源,综上,选用0.5 mm的计算网格,既保证了计算准确性,又不占用过多计算资源。
3.4 氢氟爆轰的FLUENT计算结果
计算得到的爆轰波传播结果如图5所示。从图5中可以看出,预混气体在被点火激发之后,反应向外扩散,压力逐渐升高,爆炸波在22 μs时到达管壁处;在26 μs时整个管道外围气体被激发反应,开始出现从边缘向轴线方向的反应,在29 μs时能明显的看到管壁处压力升高;在43 μs时爆炸波传播到了管道轴线处,此时轴线上的超压达到峰值;随后爆炸波沿径向传播,其中在52 μs时可以看到爆炸波向管壁传播;在63 μs时爆炸波传播到管壁,此时轴线上压力很小,管壁处压力很高;随后爆炸波向轴线方向传播,在69 μs时传播到轴线上,此时轴线压力又达到一个峰值,出现了压力的波动现象。
图5 管道内爆轰波传播云图
Fig.5 Cloud diagram of detonation wave propagation in pipeline
管道轴线上不同距离处计算出的温度(峰值温度、稳定温度)与超压(峰值超压、稳定超压)结果如图6所示。
图6 管道内轴线上不同距离处的温度与超压结果
Fig.6 Temperature and overpressure results at different distances along the inner axis of the pipeline
从图6中可以看出,在管道轴线上不同距离处的峰值温度都在6 800 K左右;最后轴线不同距离处的温度都稳定在5 000 K。
从图6中可以看出,在管道轴线上不同距离处的峰值超压虽有波动,但大都集中在46 MPa左右;最后轴线不同距离处的超压都稳定在2 MPa左右。
4 结论
1) 基于VLW爆轰产物状态方程,与爆轰基础理论分析了氢气氟气在标准大气压下的预混状态爆轰性能参数,发现随着氟气摩尔占比的提高,氢氟爆轰性能出现先增加后减小的趋势;当氢气氟气摩尔比1∶1时爆轰性能最好,此时爆速为8 044 m/s,爆压为15.31 MPa,爆温为6 479 K。
2) 通过materials studio与CP/2K进行分子动力学模拟研究了氢与氟的反应,完善现有氢氟反应机理并通过FLUENT流体动力学仿真软件对摩尔比1∶1时的氢氟爆轰进行了仿真计算,因爆轰波在管道内多次叠加,所以计算出爆压高达46 MPa,高于理论爆压值的15.31 MPa;管道内计算出的温度可达爆温6 800 K,与理论爆温误差4.6 %,计算结果与理论计算结果较好的验证。
3) 氢氟爆轰在本文中管道中的行为表现为:爆炸冲击波形成后波阵面上超压逐渐增大,在22 μs传播到管壁上发生反射;在29 μs时管道内气体整体被激发,反应从管壁向轴线处传播;随后爆轰波在管道内传播,最后管道内的温度场与压力场趋于稳定。
[1] KAPRALOVA G A,MARGOLINA E M,CHAIKIN A M.The effect of fluorine,hydrogen,or inert gases on the upper self-ignition limit inthe reaction of fluorine with hydrogen[J].Kinet.Catal,1969(10):23-27.
[2] LEVY J B,COPELAND B K W.The kinetics of the hydrogen fluorine reaction.ii.the oxygen-inhibited reaction[J].Phys.Chem,1965,69:408-416.
[3] ORKIN V L,CHAIKIN A M.Determination of the rate constant for the reaction of vibrationally excited hydrogen with molecular fluorine[J].Kinet.Khim.Reakts,1980,6:63-66.
[4] AKIRA M,HIROUMI S,KENTARO T,et al.Chain reaction mechanism in hydrogen/fluorine combustion[J].The Journal of Physical Chemistry A,2013,117(51):14042-14047.
[5] WERNER H J,KNOWLES P J,KNIZIA G,et al.Molpro:A general-purpose quantum chemistry programpackage[J].WIREs Comput.Mol.Sci.,2012(2):242-253.
[6] YURI BEDJANIAN.Rate coefficients of the reactions of fluorine atoms with h2s and sh over the temperature range 220-960 K[J].Molecules,2022,27(23):8365.
[7] ZHANG Y.Flow process and energy release of hydrogen in fluorine[J].International Journal of Hydrogen Energy,2022,48(5):2044-2054.
[8] 钱学森.物理力学讲义(新世纪版)[M].上海:上海交通大学出版社,2007:253-254.QIAN Xuesen.Notes on physical mechanics[M].Shanghai:Shanghai Jiao Tong University Press,2007:253-254.
[9] 吴雄.应用VLW状态方程计算炸药的爆轰参数[J].兵工学报,1985(3):11-20.WU Xiong.Detonation properties of condensed explosives compacted with the vlw equation of state[J].Acta Armamentarii,1985(3):11-20.
[10] 吴雄.VLW爆轰产物状态方程的发展及应用[J].火炸药学报,2021,44(1):1-7.WU Xiong.Development and application of vlw equation of state for detonation products[J].Chinese Journal of Explosives &Propellants,2021,44(1):1-7.
[11] 张金龙,郭子如,杜明燃,等.RDX的爆炸产物组成和爆热的计算与分析[J].煤矿爆破,2019,37(4):24-27.ZHANG Jinlong,GUO Ziru,DU Mingran,et al.Calculation and analysis of explosive product composition and heat of RDX[J].Coal Mine Blasting,2019,37(4):24-27.
[12] YUE P.A new high order virial equation of state and its application in the chapman-jouguet parameters calculation of explosives[J].Propellants Explosives Pyrotechnics,2021,371.
[13] YUE P.Performance investigation of a copper azide micro detonator using experiments and simulations[J].Propellants,Explosives,Pyrotechnics,2021,343.
[14] YUE P.A new simplified virial equation of state for high temperature and high pressure gas[J].Aip Advances,2022,12(1):015119.
[15] 保福成,彭汝芳,张朝阳.面向高通量含能分子设计筛选的三种生成焓计算方法评估[J].含能材料,2022,30(7):726-735.BAO Fucheng,PENG Rufang,ZHANG Chaoyang.Evaluation of three heat of formation calculation methods for high-throughput energetic molecule design and screening[J].Chinese Journal of Energetic Materials,2022,30(7):726-735.
[16] ZHAO Y.Evaluation of detonation performance and working capacity of explosives by optimized vlw eos[J].Combustion and Flame,2022,235.