自催化反应速率方程的导出途径及其在含能材料热行为研究中的应用

胡荣祖, 姚二岗, 张 海, 王璞玉, 赵宏安, 赵凤起, 高红旭, 罗 阳, 马海霞

(1. 西安近代化学研究所燃烧与爆炸技术重点实验室, 陕西 西安 710065; 2. 西北大学数学系/数据分析和计算化学研究所, 陕西 西安 710069; 3. 西北大学信息科学与工程学院, 陕西 西安 710069; 4. 西北大学化工学院, 陕西 西安 710069)

1 引 言

含能材料(EM)体系自(动)催化反应的动力学行为和动力学三因子,在EM对热抵抗能力的评估,热爆炸临界温度、热爆炸临界温升速率、撞击感度(特性落高)、放热系统热感度、绝热至爆时间、燃烧速度的估算和放热分解反应诱导温度与诱导时间关系的定量描述方面扮演重要角色[1-6]。导出该类反应的速率方程、速率曲线方程和反应进度(α)随时间(t)和温度(T)变化的方程,对描述、定量评估自催化反应的动力学行为,有一定学术意义。本文依据EM体系自催化反应的特性: (1)自催化反应由催化剂生成反应和催化剂催化EM的催化反应组成,体系自催化反应速率则是催化剂生成反应速率和催化反应速率的加和; (2)引导起始催化反应不能从反应进度α=0开始,需要有催化产物; (3)具有催化功能的反应产物使反应经过一段诱导期后才能出现反应加速; (4)反应速率随某一催化反应产物浓度而增长; (5)自催化反应速率有最大值,从α与反应能量变化的关系,导出了经验级数自催化反应的速率方程。由经验级数自催化反应速率方程,导出了13个自催化反应速率的派生式。由自催化反应特性(5)导出了描述自催化反应速率曲线特性[αmax和(dα/dt)max]的方程。通过速率方程的积分处理,导出了α随t和T变化的方程。报道了描述六硝基六氮杂异伍兹烷(HNIW)自催化分解反应的速率方程和硝化棉(NC)(12.82%、12.97%、13.54%、13.61%、13.86%、13.88%、14.14% N)自催化分解反应的动力学参数—催化系数Kcat、速率曲线特性参数和α随t变化的方程。

2 理论和方法

2.1 经验级数自催化反应速率方程的导出途径

EM体系的自(动)催化反应由催化剂(B)的生成反应和催化剂催化EM的催化反应组成,体系自催化反应速率则是催化剂生成反应速率和催化反应速率的加和。

由催化剂生成反应

(1)

和边界条件

(2)

式(1)中,下角标“p”代表产物; “CFR”代表催化剂生成反应。

知反应进度(α或c)与反应能量变化的关系[7]:

(3)

(4)

c=c0(1-α)

(5)

c0-c=c0α

(6)

和催化剂生成速率方程的微分式

(7)

(8)

(9)

令

(10)

则

(11)

式中,c、α、Q和k有通常的含义。

由催化反应

nEM+B→nEMB

(12)

和边界条件

(13)

知方程(3)～(6)和催化反应速率方程的微分式

(14)

(15)

(16)

(17)

(18)

式中

(19)

反应初期,α≪1,式(18)可变为

(20)

α-ndα=k2dt

(21)

方程两边积分

(22)

得

(23)

据此不难看出,引导起始反应不能从α=0开始,需要有催化产物,催化方程(18)应写为

(24)

考虑体系中还独立地进行着从k1(1-α)m速率给出催化产物的反应,因此,自催化方程应写为

(25)

我们称式(25)为经验级数自催化反应速率方程。

2.2 经验级数自催化速率方程派生式的导出途径

由经验级数(m、n、p)速率方程(25)导出的13个派生式见图1。其中,式(25-7)～(25-13)为自催化速率方程; 式(25-1)、(25-2)、(25-3)和(25-4)为催化速率方程; 式(25-2)、(25-3)和(25-4)分别称为第三类微分方程式、第二类微分方程式和第一类微分方程式。

式(25-10)和(25-11)中的Kcat称催化系数,特指自催化方程dα/dt=k1(1-α)m+k2αn(1-α)p中k2与k1的比值:

对框图1中各方程的参数: 方程(25-2)、(25-3)、(25-4)和(25-6)中的2参数(A2、Ea2;A2、Ea2;A2、Ea2;A2、Ea2),方程(25-5)和(25-11)中的3参数(A2、Ea2、p;A1、Ea1、Kcat),方程(25-1)、(25-7)、(25-8)、(25-10)、(25-12)和(25-13)中的4参数(分别为A2、Ea2、n,p;A1、A2、Ea1、Ea2;A1、A2、Ea1、Ea2;A1、Ea1、m、Kcat;A1、A2、Ea1、Ea2;A1、A2、Ea1、Ea2),方程(25-9)中的5参数(A1、A2、Ea1、Ea2、m)和方程(25)中的7参数(A1、A2、Ea1、Ea2、m、n、p)可用线性最小二乘法和信赖域方法求解[4-6,8],也可将方程(25-13)改写为

和

的形式,通过解伯努利方程(Bernoulli方程),得4参数(A1、A2、Ea1、Ea2)[9]。

这里,y是未反应物质的分数; 1-y是已发生反应物质的分数;T为温度;t为时间;k1(T)和k2(T)是一阶自催化方程在温度T的速率常数;k1(T)=A1exp(-E1/RT),k2(T)=A2exp(-E2/RT); 为指前因子;E为活化能; dt=dT/β,β是线性升温速率。

图1经验级数自催化反应速率方程和派生式的导出过程框图

Fig.1Block diagram of the derivation process of empiric-order autocatalytic reaction rate equation and its thirteen-derived formulate

2.3 自催化反应速率曲线特性—αmax和(dα/dt)max—表达式的导出途径

由方程(25)两边对t求导

k2p(α+α0)n(1-α)p-1]

(26)

和d2α/dt2=0,速率达最大值,得

k2n(αmax+α0)n-1(1-αmax)p

=k1m(1-αmax)m-1+k2p(αmax+α0)n(1-αmax)p-1

(27)

由m=0,得

(28)

n(1-αmax)=p(αmax+α0)

(29)

n-nαmax=pαmax+pα0

(30)

(n+p)αmax=n-pα0

(31)

(32)

将αmax表达式(32)代入方程(25-1),得

(33)

当α0≪1时

(34)

(35)

对n=2,p=1的派生式(25-2),由式(32)和(33)知

(36)

(37)

当α0≪1时

(38)

(39)

对n=1,p=2的派生式(25-3),由式(32)和(33)知

(40)

(41)

当α0≪1时

(42)

(43)

对n=p=1的派生式(25-4)

(44)

(45)

当α0≪1时

(46)

(47)

对m=0,n=p=1的派生式(25-7),由式(27)知

k2(1-αmax)=k2(αmax+α0)

(48)

1-αmax=αmax+α0

1-α0=2αmax

(49)

(50)

当α0≪1时

(51)

(52)

对m=n=0,p=1的派生式(25-8),由式(27)知

k2=0

(53)

(54)

对m=n=p=1的派生式(25-13),由式(27)知

k2(1-αmax)=k1+k2(αmax+α0)

(55)

k2-k2αmax=k1+k2αmax+k2α0

2k2αmax=k2-k1-k2α0

(56)

=(1-αmax)[k1+k2(αmax+α0)]

(57)

令Kcat=k2/k1,则

(58)

当α0≪1时

(59)

(60)

2.4 自催化反应的α随t变化方程的导出途径

对方程(25-13),由

(61)

知

(62)

得

(63)

(64)

k1+k2(α+α0)=(1-α)(k1+k2α0)e[k1+k2(1+α0)]t

(65)

k1+k2α0+k2α= (k1+k2α0)e[k1+k2(1+α0)]t-

(k1+k2α0)αe[k1+k2(1+α0)]t

(66)

k2α+(k1+k2α0)αe[k1+k2(1+α0)]t

=(k1+k2α0)e[k1+k2(1+α0)]t-(k1+k2α0)

(67)

(68)

(69)

(70)

对α0=0.0001≈0的方程(25-13),由

(71)

知

(72)

得

(73)

(74)

k1+k2α=k1(1-α)e(k1+k2)t

(75)

k2α+k1αe(k1+k2)t=k1e(k1+k2)t-k1

(76)

(77)

(78)

对方程(25-8),由

(79)

知

(80)

(81)

得

(82)

(83)

(84)

k1+k2(1-α)=(k1+k2)e-k2t

(85)

(86)

对方程(25-11),由

(87)

知

(88)

得

(89)

(90)

1+Kcatα+Kcatα0=(1-α)(1+Kcatα0)ek1(1+Kcat+Kcatα0)t

(91)

1+Kcatα+Kcatα0= (1+Kcatα0)ek1(1+Kcat+Kcatα0)t-

(1+Kcatα0)ek1(1+Kcat+Kcatα0)tα

(92)

Kcatα+(1+Kcatα0)ek1(1+Kcat+Kcatα0)tα

=(1+Kcatα0)ek1(1+Kcat+Kcatα0)t-(1+Kcatα0)

(93)

(94)

对α0=0.0001≈0的方程(25-11),由

(95)

知

(96)

得

(97)

(98)

1+Kcatα=(1-α)ek1(1+Kcat)t

(99)

1+Kcatα=ek1(1+Kcat)t-ek1(1+Kcat)tα

(100)

Kcatα+ek1(1+Kcat)tα=ek1(1+Kcat)t-1

(101)

(102)

对方程(25-4),由

(103)

知

(104)

得

(105)

(106)

α+α0=α0(1-α)ek2(α0+1)t

(107)

[1+α0ek2(α0+1)t]α=α0[ek2(α0+1)t-1]

(108)

(109)

反应初期,α0≪1

α=α0[ek2(α0+1)t-1]≈α0ek2(α0+1)t

(110)

显示反应随时间呈指数变化的规律。

对α0≪1,α0=0.0001≈0的方程(25-4),由

(111)

知

(112)

得

(113)

(114)

9999α=(1-α)ek2t

(115)

(9999+ek2t)α=ek2t

(116)

(117)

对方程(25-2),由

(118)

知

(119)

得

(120)

对α0≪1,α0=0.0001≈0的方程(25-2),由

(121)

知

(122)

得

(123)

(124)

(125)

对方程(25-3),由

(126)

知

(127)

得

(128)

(129)

对α0≪1,α0=0.0001≈0的方程(25-3),由

(130)

知

(131)

得

(132)

(133)

(134)

对方程(25-5),由

(135)

知

(136)

得

(137)

(138)

(139)

对方程(25-6),由

(140)

知

(141)

得

(142)

ln(1-α)=-k2t

(143)

α=1-e-k2t

(144)

2.5 自催化反应的α随T变化方程的导出途径

对第Ⅰ类动力学方程[11]

(145)

由Frank-Kameneskii近似式[10]

(146)

知

(147)

由第三类微分方程式(25-2),知

f(α)=(α+α0)2(1-α)

(148)

(149)

由式(147)和式(149)联立,得

(150)

由第二类微分方程式方程(25-3),知

f(α)=(α+α0)(1-α)2

(151)

(152)

由式(147)和式(152)联立,得

(153)

由第一类微分方程式方程(25-4),知

f(α)=(α+α0)(1-α)

(154)

(155)

由式(147)和式(155)联立,得

(156)

(157)

(158)

(159)

(160)

对第Ⅱ类动力学方程[11]

(161)

由

(162)

知

(163)

由式(149)和(163)式联立,得

(164)

由式(152)和式(163)联立,得

(165)

由式(155)和式(163)联立,得

(166)

(167)

(168)

(169)

(170)

由Harcourt-Esson方程的积分式(171)[12]

(171)

与式(149)联立,得

(172)

由式(152)和式(171)联立,得

(173)

由式(155)和式(171)联立,得

(174)

(175)

(176)

(177)

(178)

由Berthelot方程的积分式(179)[12]

(179)

与式(149)式联立,得

(180)

由式(152)和式(179)联立,得

(181)

由式(155)和式(179)联立,得

(182)

(183)

(184)

(185)

(186)

由积分式(187)[11]

(187)

与式(149)联立,得

(188)

由式(152)和式(187)联立,得

(189)

由式(155)和式(187)联立,得

(190)

(191)

(192)

(193)

(194)

由积分式(195)

(195)

与式(149)联立,得

(196)

由式(152)和式(195)联立,得

(197)

由式(155)和式(195)联立,得

(198)

(199)

(200)

(201)

3 计算实例

(202)

的非线性优化模型

(203)

s.t. 1012.5s-1≤A≤1013.5s-1, 150000 J·mol-1≤E≤165000 J·mol-1, 0.9≤n≤1.1, 10≤Kcat≤15

由线性最小二乘法和信赖域方法[4-6, 8]得模型中的4参数(A、E、n、Kcat)。

从计算结果可知,Δδ很小,表明非等温条件下,用CnB速率方程

dα/dt=1013.46exp(-160230/RT)(1-α)0.90(1+15α)

描述HNIW的热分解过程是可取的。若视n=0.90≈1,

表1由DSC测得HNIW的数据[13]

Table1Data of HNIW determined by DSC[13]

datapointαiTi/K2K·min-15K·min-110K·min-120K·min-110473.20483.05491.58498.5620.025488.29499.59508.08516.0230.050490.63502.22510.46518.7840.075491.96503.83512.03520.4950.100492.97504.98513.19521.7760.125493.75505.88514.13522.8270.150494.45506.67514.93523.6980.175495.06507.35515.63524.4790.200495.61507.96516.26525.16100.225496.13508.52516.84525.79110.250496.66509.04517.37526.38120.275497.09509.54517.86526.93130.300497.49510.02518.34527.45140.325497.84510.48518.79527.94150.350498.19510.91519.23528.42160.375498.51511.32519.65528.87170.400498.85511.71520.06529.29180.425499.17512.09520.45529.70190.450499.50512.46520.84530.10200.475499.82512.83521.22530.48210.500500.14513.19521.59530.85220.525500.43513.53521.95531.20230.550500.71513.88522.31531.56240.575501.00514.20522.65531.89250.600501.29514.53522.99532.23260.625501.59514.86523.31532.55270.650501.88515.17523.64532.87280.675502.14515.48523.94533.17290.700502.42515.78524.25533.47

则可用方程(94)和(102)描述α随t的变化。

类似地,由描述NC(12.82% N)自催化分解反应的速率方程[3]

dα/dt= 1016.4exp(-178000/RT)(1-α)+

1017.0exp(-174000/RT)α(1-α)

可得:Tmax=479.55 K,αmax=0.4539,Kcat=10.86。

由描述NC(12.97% N)自催化分解反应的速率方程[5]

dα/dt= 1016.00exp(-174520/RT)(1-α)+

1016.00exp(-163510/RT)α(1-α)

可得:Tmax=481.19 K,αmax=0.4680,Kcat=15.68。

由描述NC(13.54% N)自催化分解反应的速率方程[4]

dα/dt= 1015.82exp(-170020/RT)(1-α)1.11+

1015.82exp(-157140/RT)α1.51(1-α)2.51

可得:Tmax=482.25 K,αmax=0.3411,Kcat=24.84。

由描述NC(13.61% N)自催化分解反应的速率方程[3]

dα/dt= 1016.5exp(-184700/RT)(1-α)+

1016.9exp(-174700/RT)α(1-α)

可得:Tmax=478.75 K,αmax=0.4838,Kcat=30.98。

由描述NC(13.86% N)自催化分解反应的速率方程[14]

dα/dt= 1016.30exp(-181860/RT)(1-α)+

1016.70exp(-173050/RT)α(1-α)

可得:Tmax=477.00 K,αmax=0.4783,Kcat=23.16。

由描述NC(13.88% N)自催化分解反应的速率方程[3]

dα/dt= 1016.40exp(-181860/RT)(1-α)+

1016.70exp(-171730/RT)α(1-α)

可得:Tmax=476.95 K,αmax=0.4805,Kcat=25.67。

由描述NC(14.14% N)自催化分解反应的速率方程[6]

dα/dt= 1015.76exp(-170800/RT)(1-α)0.95+

1015.76exp(-159100/RT)α1.81(1-α)1.16

可得:Tmax=489.61 K,αmax=0.5756,Kcat=17.71。

上述计算结果表明,非等温条件下,NC(13.54%、14.14% N)的热分解过程可用经验级数自催化反应速率方程(25)描述。NC(12.82%、12.97%、13.61%、13.86%、13.88% N)的热分解过程可用一级(m=1、n=1、p=1)自催化分解反应速率方程(25-13)描述。对一级自催化分解反应,可用方程(56)和(59)描述αmax和Kcat的关系,可用方程(57)和(58)描述(dα/dt)max和αmax、Kcat、E1/RTmax的关系,可用方程(70)和(77)描述α随的t变化。E2大于E1,表明: 催化分解反应速率大于催化剂生成反应速率,催化产物存在下的催化分解反应易于发生。

4 结 论

(1) 导出了经验级数自催化反应速率方程和13个派生式。提出了描述自催化分解反应速率曲线特性的方程,以及反应进度随时间和温度变化的方程。

(2) 描述HNIW分解过程的速率方程为自催化的n级反应(CnB)速率方程:

(3) 提出了正文中描述NC(12.82%、12.97%、13.54%、13.61%、13.86%、13.88%、14.14% N)自催化分解过程的反应动力学参数—催化系数Kcat、速率曲线特性参数和α随t变化的方程。

参考文献:

[1] Eisenreich, N. Beitrag zur kinetic thermischer zersetzungseaktionen (Thermoanalytische auswertungder zersetzung von nitrocellulose)[D]. Munich : Technical University of Munich, 1978.

[2] Eisenreich N, Pfeil A. Non-linear least-squares fit of non-isothermal thermoanalytical curves. Reinvestigation of the kinetics of the autocatalytic decomposition of nitrated cellulose[J].ThermochimActa, 1983, 61(1-2): 13-21.

[3] 宁斌科, 杨正权, 刘蓉,等. 硝化棉一级自催化分解反应动力学参数数值模拟[J]. 含能材料, 1999, 7(4): 162-165.

NING Bin-ke, YANG Zheng-quan, LIU Rong, et al. Numerical simulation of kinetic parameters of the first-order autocatalytic decomposition of NC[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 1999, 7(4): 162-165.

[4] 胡荣祖, 赵凤起, 高红旭, 等. 从DSC曲线数据计算/确定含能材料自催化分解反应动力学参数和热爆炸临界温升速率的方法[J]. 物理化学学报, 2013, 29(8): 1623-1631.

HU Rong-zu, ZHAO Feng-qi, GAO Hong-xu, et al. A method of computing/determining the kinetic parameters of the autocatalytic decomposing reaction and the critical rate of temperature rise of thermal explosion of energetic materials from data of DSC curves at different heating rate[J].ActaPhys-ChimSin, 2013, 29(8): 1623-1631.

[5] WANG Hui-ya, ZHANG Hai, HU Rong-zu, et al. Estimation of the critical rate of temperature rise for thermal explosion of apparent empiric-order autocatalytic decomposing reaction of nitrocellulose (12.97% N) using non-isothermal DSC[J].JouranalofThermalAnalysisandCalorimetry, 2014, 115: 1099-1110.

[6] 姚二岗, 胡荣祖, 赵凤起, 等. 用非等温DSC估算硝化棉(14.14% N)的CnB和表观经验级数自催化分解反应热爆炸临界温升速率[J]. 火炸药学报, 2013, 36(5): 72-76.

YAO Er-gang, HU Rong-zu, ZHAO Feng-qi, et al. Estimation of the critical rate of temperature rise for thermal explosion ofCnBand apparent empiric-order autocatalytic decomposing reaction of nitrocellulose (14.14% N) using non-isothermal DSC[J].ChineseJournalofExplosives&Propellants, 2013, 36(5): 72-76.

[7] 邓郁. 热动力学的研究(I): 热动力学的理论和方法[J]. 高等学校化学学报, 1985, 6(7): 621-626.

DENG Yu. Studies on thermokinetics (I): Theory and method of thermokinetics[J].ChemJChineseU, 1985, 6(7): 621-626.

[8] Coleman T F, Li Yuying. An interior trust region approach for nonlinear minimization subject to bounds[J].SIAMJournalonOptimization, 1996, 6(2): 418-445.

[9] 胡荣祖, 史启祯. 热分析动力学[M]. 北京: 科学出版社, 2001: 180-191.

HU Rong-zu, SHI Qi-zheng. Thermal Analysis Kinetics[M]. Beijing: Science Press, 2001: 180-191.

[10] Frank-Kameneshii D. A. Diffusion and Heat Exchange in Chemical Kinetics[M]. Princeton: Princeton University Press, 1955.

[11] 胡荣祖, 高胜利, 赵凤起, 等. 热分析动力学[M]. 北京: 科学出版社, 2008: 20-26.

HU Rong-zu, GAO Sheng-li, ZHAO Feng-qi, et al. Thermal Analysis Kinetics[M]. Beijing: Science Press, 2008: 20-26

[12] 胡荣祖, 赵凤起, 高红旭, 等. 量热学基础与应用[M]. 北京: 科学出版社, 2011: 342-351.

HU Rong-zu, ZHAO Feng-qi, GAO Hong-xu, et al. Fundamentals and Application of Calorimetry[M]. Beijing: Science Press, 2011: 342-351.

[13] GAO Hong-xu, ZHAO Feng-qi, HU Rong-zu, et al. Estimation of the kinetic parameters of thermal decomposition reaction and thermal safety evaluation on hexanitrohexaazaisowurtzitane(HNIW)[J].DefenceTchnology, 2014, 10: 28-33.

[14] GUO Peng-jiang, HU Rong-zu, NING Bing-ke, et al. Kinetics of the first order autocatalytic decomposition reaction of nitrocellulose (13.86% N)[J].ChinJChem, 2004, 22: 19-23.