胡榮祖, 姚二崗, 張 海, 王璞玉, 趙宏安, 趙鳳起, 高紅旭, 羅 陽(yáng), 馬海霞
(1. 西安近代化學(xué)研究所燃燒與爆炸技術(shù)重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710065; 2. 西北大學(xué)數(shù)學(xué)系/數(shù)據(jù)分析和計(jì)算化學(xué)研究所, 陜西 西安 710069; 3. 西北大學(xué)信息科學(xué)與工程學(xué)院, 陜西 西安 710069; 4. 西北大學(xué)化工學(xué)院, 陜西 西安 710069)
含能材料(EM)體系自(動(dòng))催化反應(yīng)的動(dòng)力學(xué)行為和動(dòng)力學(xué)三因子,在EM對(duì)熱抵抗能力的評(píng)估,熱爆炸臨界溫度、熱爆炸臨界溫升速率、撞擊感度(特性落高)、放熱系統(tǒng)熱感度、絕熱至爆時(shí)間、燃燒速度的估算和放熱分解反應(yīng)誘導(dǎo)溫度與誘導(dǎo)時(shí)間關(guān)系的定量描述方面扮演重要角色[1-6]。導(dǎo)出該類反應(yīng)的速率方程、速率曲線方程和反應(yīng)進(jìn)度(α)隨時(shí)間(t)和溫度(T)變化的方程,對(duì)描述、定量評(píng)估自催化反應(yīng)的動(dòng)力學(xué)行為,有一定學(xué)術(shù)意義。本文依據(jù)EM體系自催化反應(yīng)的特性: (1)自催化反應(yīng)由催化劑生成反應(yīng)和催化劑催化EM的催化反應(yīng)組成,體系自催化反應(yīng)速率則是催化劑生成反應(yīng)速率和催化反應(yīng)速率的加和; (2)引導(dǎo)起始催化反應(yīng)不能從反應(yīng)進(jìn)度α=0開始,需要有催化產(chǎn)物; (3)具有催化功能的反應(yīng)產(chǎn)物使反應(yīng)經(jīng)過一段誘導(dǎo)期后才能出現(xiàn)反應(yīng)加速; (4)反應(yīng)速率隨某一催化反應(yīng)產(chǎn)物濃度而增長(zhǎng); (5)自催化反應(yīng)速率有最大值,從α與反應(yīng)能量變化的關(guān)系,導(dǎo)出了經(jīng)驗(yàn)級(jí)數(shù)自催化反應(yīng)的速率方程。由經(jīng)驗(yàn)級(jí)數(shù)自催化反應(yīng)速率方程,導(dǎo)出了13個(gè)自催化反應(yīng)速率的派生式。由自催化反應(yīng)特性(5)導(dǎo)出了描述自催化反應(yīng)速率曲線特性[αmax和(dα/dt)max]的方程。通過速率方程的積分處理,導(dǎo)出了α隨t和T變化的方程。報(bào)道了描述六硝基六氮雜異伍茲烷(HNIW)自催化分解反應(yīng)的速率方程和硝化棉(NC)(12.82%、12.97%、13.54%、13.61%、13.86%、13.88%、14.14% N)自催化分解反應(yīng)的動(dòng)力學(xué)參數(shù)—催化系數(shù)Kcat、速率曲線特性參數(shù)和α隨t變化的方程。
EM體系的自(動(dòng))催化反應(yīng)由催化劑(B)的生成反應(yīng)和催化劑催化EM的催化反應(yīng)組成,體系自催化反應(yīng)速率則是催化劑生成反應(yīng)速率和催化反應(yīng)速率的加和。
由催化劑生成反應(yīng)
(1)
和邊界條件
(2)
式(1)中,下角標(biāo)“p”代表產(chǎn)物; “CFR”代表催化劑生成反應(yīng)。
知反應(yīng)進(jìn)度(α或c)與反應(yīng)能量變化的關(guān)系[7]:
(3)
(4)
c=c0(1-α)
(5)
c0-c=c0α
(6)
和催化劑生成速率方程的微分式
(7)
(8)
(9)
令
(10)
則
(11)
式中,c、α、Q和k有通常的含義。
由催化反應(yīng)
nEM+B→nEMB
(12)
和邊界條件
(13)
知方程(3)~(6)和催化反應(yīng)速率方程的微分式
(14)
(15)
(16)
(17)
(18)
式中
(19)
反應(yīng)初期,α?1,式(18)可變?yōu)?/p>
(20)
α-ndα=k2dt
(21)
方程兩邊積分
(22)
得
(23)
據(jù)此不難看出,引導(dǎo)起始反應(yīng)不能從α=0開始,需要有催化產(chǎn)物,催化方程(18)應(yīng)寫為
(24)
考慮體系中還獨(dú)立地進(jìn)行著從k1(1-α)m速率給出催化產(chǎn)物的反應(yīng),因此,自催化方程應(yīng)寫為
(25)
我們稱式(25)為經(jīng)驗(yàn)級(jí)數(shù)自催化反應(yīng)速率方程。
由經(jīng)驗(yàn)級(jí)數(shù)(m、n、p)速率方程(25)導(dǎo)出的13個(gè)派生式見圖1。其中,式(25-7)~(25-13)為自催化速率方程; 式(25-1)、(25-2)、(25-3)和(25-4)為催化速率方程; 式(25-2)、(25-3)和(25-4)分別稱為第三類微分方程式、第二類微分方程式和第一類微分方程式。
式(25-10)和(25-11)中的Kcat稱催化系數(shù),特指自催化方程dα/dt=k1(1-α)m+k2αn(1-α)p中k2與k1的比值:
對(duì)框圖1中各方程的參數(shù): 方程(25-2)、(25-3)、(25-4)和(25-6)中的2參數(shù)(A2、Ea2;A2、Ea2;A2、Ea2;A2、Ea2),方程(25-5)和(25-11)中的3參數(shù)(A2、Ea2、p;A1、Ea1、Kcat),方程(25-1)、(25-7)、(25-8)、(25-10)、(25-12)和(25-13)中的4參數(shù)(分別為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參數(shù)(A1、A2、Ea1、Ea2、m)和方程(25)中的7參數(shù)(A1、A2、Ea1、Ea2、m、n、p)可用線性最小二乘法和信賴域方法求解[4-6,8],也可將方程(25-13)改寫為
和
的形式,通過解伯努利方程(Bernoulli方程),得4參數(shù)(A1、A2、Ea1、Ea2)[9]。
這里,y是未反應(yīng)物質(zhì)的分?jǐn)?shù); 1-y是已發(fā)生反應(yīng)物質(zhì)的分?jǐn)?shù);T為溫度;t為時(shí)間;k1(T)和k2(T)是一階自催化方程在溫度T的速率常數(shù);k1(T)=A1exp(-E1/RT),k2(T)=A2exp(-E2/RT); 為指前因子;E為活化能; dt=dT/β,β是線性升溫速率。
圖1經(jīng)驗(yàn)級(jí)數(shù)自催化反應(yīng)速率方程和派生式的導(dǎo)出過程框圖
Fig.1Block diagram of the derivation process of empiric-order autocatalytic reaction rate equation and its thirteen-derived formulate
由方程(25)兩邊對(duì)t求導(dǎo)
k2p(α+α0)n(1-α)p-1]
(26)
和d2α/dt2=0,速率達(dá)最大值,得
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表達(dá)式(32)代入方程(25-1),得
(33)
當(dāng)α0?1時(shí)
(34)
(35)
對(duì)n=2,p=1的派生式(25-2),由式(32)和(33)知
(36)
(37)
當(dāng)α0?1時(shí)
(38)
(39)
對(duì)n=1,p=2的派生式(25-3),由式(32)和(33)知
(40)
(41)
當(dāng)α0?1時(shí)
(42)
(43)
對(duì)n=p=1的派生式(25-4)
(44)
(45)
當(dāng)α0?1時(shí)
(46)
(47)
對(duì)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)
當(dāng)α0?1時(shí)
(51)
(52)
對(duì)m=n=0,p=1的派生式(25-8),由式(27)知
k2=0
(53)
(54)
對(duì)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)
當(dāng)α0?1時(shí)
(59)
(60)
對(duì)方程(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)
對(duì)α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)
對(duì)方程(25-8),由
(79)
知
(80)
(81)
得
(82)
(83)
(84)
k1+k2(1-α)=(k1+k2)e-k2t
(85)
(86)
對(duì)方程(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)
對(duì)α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)
對(duì)方程(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)
反應(yīng)初期,α0?1
α=α0[ek2(α0+1)t-1]≈α0ek2(α0+1)t
(110)
顯示反應(yīng)隨時(shí)間呈指數(shù)變化的規(guī)律。
對(duì)α0?1,α0=0.0001≈0的方程(25-4),由
(111)
知
(112)
得
(113)
(114)
9999α=(1-α)ek2t
(115)
(9999+ek2t)α=ek2t
(116)
(117)
對(duì)方程(25-2),由
(118)
知
(119)
得
(120)
對(duì)α0?1,α0=0.0001≈0的方程(25-2),由
(121)
知
(122)
得
(123)
(124)
(125)
對(duì)方程(25-3),由
(126)
知
(127)
得
(128)
(129)
對(duì)α0?1,α0=0.0001≈0的方程(25-3),由
(130)
知
(131)
得
(132)
(133)
(134)
對(duì)方程(25-5),由
(135)
知
(136)
得
(137)
(138)
(139)
對(duì)方程(25-6),由
(140)
知
(141)
得
(142)
ln(1-α)=-k2t
(143)
α=1-e-k2t
(144)
對(duì)第Ⅰ類動(dòng)力學(xué)方程[11]
(145)
由Frank-Kameneskii近似式[10]
(146)
知
(147)
由第三類微分方程式(25-2),知
f(α)=(α+α0)2(1-α)
(148)
(149)
由式(147)和式(149)聯(lián)立,得
(150)
由第二類微分方程式方程(25-3),知
f(α)=(α+α0)(1-α)2
(151)
(152)
由式(147)和式(152)聯(lián)立,得
(153)
由第一類微分方程式方程(25-4),知
f(α)=(α+α0)(1-α)
(154)
(155)
由式(147)和式(155)聯(lián)立,得
(156)
(157)
(158)
(159)
(160)
對(duì)第Ⅱ類動(dòng)力學(xué)方程[11]
(161)
由
(162)
知
(163)
由式(149)和(163)式聯(lián)立,得
(164)
由式(152)和式(163)聯(lián)立,得
(165)
由式(155)和式(163)聯(lián)立,得
(166)
(167)
(168)
(169)
(170)
由Harcourt-Esson方程的積分式(171)[12]
(171)
與式(149)聯(lián)立,得
(172)
由式(152)和式(171)聯(lián)立,得
(173)
由式(155)和式(171)聯(lián)立,得
(174)
(175)
(176)
(177)
(178)
由Berthelot方程的積分式(179)[12]
(179)
與式(149)式聯(lián)立,得
(180)
由式(152)和式(179)聯(lián)立,得
(181)
由式(155)和式(179)聯(lián)立,得
(182)
(183)
(184)
(185)
(186)
由積分式(187)[11]
(187)
與式(149)聯(lián)立,得
(188)
由式(152)和式(187)聯(lián)立,得
(189)
由式(155)和式(187)聯(lián)立,得
(190)
(191)
(192)
(193)
(194)
由積分式(195)
(195)
與式(149)聯(lián)立,得
(196)
由式(152)和式(195)聯(lián)立,得
(197)
由式(155)和式(195)聯(lián)立,得
(198)
(199)
(200)
(201)
(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參數(shù)(A、E、n、Kcat)。
從計(jì)算結(jié)果可知,Δδ很小,表明非等溫條件下,用CnB速率方程
dα/dt=1013.46exp(-160230/RT)(1-α)0.90(1+15α)
描述HNIW的熱分解過程是可取的。若視n=0.90≈1,
表1由DSC測(cè)得HNIW的數(shù)據(jù)[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)自催化分解反應(yīng)的速率方程[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)自催化分解反應(yīng)的速率方程[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)自催化分解反應(yīng)的速率方程[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)自催化分解反應(yīng)的速率方程[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)自催化分解反應(yīng)的速率方程[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)自催化分解反應(yīng)的速率方程[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)自催化分解反應(yīng)的速率方程[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。
上述計(jì)算結(jié)果表明,非等溫條件下,NC(13.54%、14.14% N)的熱分解過程可用經(jīng)驗(yàn)級(jí)數(shù)自催化反應(yīng)速率方程(25)描述。NC(12.82%、12.97%、13.61%、13.86%、13.88% N)的熱分解過程可用一級(jí)(m=1、n=1、p=1)自催化分解反應(yīng)速率方程(25-13)描述。對(duì)一級(jí)自催化分解反應(yīng),可用方程(56)和(59)描述αmax和Kcat的關(guān)系,可用方程(57)和(58)描述(dα/dt)max和αmax、Kcat、E1/RTmax的關(guān)系,可用方程(70)和(77)描述α隨的t變化。E2大于E1,表明: 催化分解反應(yīng)速率大于催化劑生成反應(yīng)速率,催化產(chǎn)物存在下的催化分解反應(yīng)易于發(fā)生。
(1) 導(dǎo)出了經(jīng)驗(yàn)級(jí)數(shù)自催化反應(yīng)速率方程和13個(gè)派生式。提出了描述自催化分解反應(yīng)速率曲線特性的方程,以及反應(yīng)進(jìn)度隨時(shí)間和溫度變化的方程。
(2) 描述HNIW分解過程的速率方程為自催化的n級(jí)反應(yīng)(CnB)速率方程:
(3) 提出了正文中描述NC(12.82%、12.97%、13.54%、13.61%、13.86%、13.88%、14.14% N)自催化分解過程的反應(yīng)動(dòng)力學(xué)參數(shù)—催化系數(shù)Kcat、速率曲線特性參數(shù)和α隨t變化的方程。
參考文獻(xiàn):
[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] 寧斌科, 楊正權(quán), 劉蓉,等. 硝化棉一級(jí)自催化分解反應(yīng)動(dòng)力學(xué)參數(shù)數(shù)值模擬[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曲線數(shù)據(jù)計(jì)算/確定含能材料自催化分解反應(yīng)動(dòng)力學(xué)參數(shù)和熱爆炸臨界溫升速率的方法[J]. 物理化學(xué)學(xué)報(bào), 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īng)驗(yàn)級(jí)數(shù)自催化分解反應(yīng)熱爆炸臨界溫升速率[J]. 火炸藥學(xué)報(bào), 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] 鄧郁. 熱動(dòng)力學(xué)的研究(I): 熱動(dòng)力學(xué)的理論和方法[J]. 高等學(xué)?;瘜W(xué)學(xué)報(bào), 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] 胡榮祖, 史啟禎. 熱分析動(dòng)力學(xué)[M]. 北京: 科學(xué)出版社, 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] 胡榮祖, 高勝利, 趙鳳起, 等. 熱分析動(dòng)力學(xué)[M]. 北京: 科學(xué)出版社, 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] 胡榮祖, 趙鳳起, 高紅旭, 等. 量熱學(xué)基礎(chǔ)與應(yīng)用[M]. 北京: 科學(xué)出版社, 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.