自催化反應速率方程的導出途徑及其在含能材料熱行為研究中的應用

胡榮祖, 姚二崗, 張 海, 王璞玉, 趙宏安, 趙鳳起, 高紅旭, 羅 陽, 馬海霞

(1. 西安近代化學研究所燃燒與爆炸技術重點實驗室, 陜西 西安 710065; 2. 西北大學數學系/數據分析和計算化學研究所, 陜西 西安 710069; 3. 西北大學信息科學與工程學院, 陜西 西安 710069; 4. 西北大學化工學院, 陜西 西安 710069)

1 引 言

含能材料(EM)體系自(動)催化反應的動力學行為和動力學三因子,在EM對熱抵抗能力的評估,熱爆炸臨界溫度、熱爆炸臨界溫升速率、撞擊感度(特性落高)、放熱系統(tǒng)熱感度、絕熱至爆時間、燃燒速度的估算和放熱分解反應誘導溫度與誘導時間關系的定量描述方面扮演重要角色[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)可變?yōu)?/p>

(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是已發(fā)生反應物質的分數;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)

顯示反應隨時間呈指數變化的規(guī)律。

對α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)

的非線性優(yōu)化模型

(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,表明: 催化分解反應速率大于催化劑生成反應速率,催化產物存在下的催化分解反應易于發(fā)生。

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]. 高等學?；瘜W學報, 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.