胡榮祖 趙鳳起 高紅旭 姚二崗 張 海 王 堯常象宇 趙宏安
(1西安近代化學(xué)研究所燃燒與爆炸技術(shù)重點實驗室,西安710065;2西北大學(xué)數(shù)學(xué)系數(shù)據(jù)分析和計算化學(xué)研究所,西安710069;3西安交通大學(xué)信息科學(xué)與系統(tǒng)科學(xué)研究所,西安710049;4西北大學(xué)信息科學(xué)與工程學(xué)院,西安710069)
放熱分解反應(yīng)體系熱爆炸的臨界溫升速率(dT/dt)Tb,是評價含能材料(EMs)安定性和安全性的重要參數(shù).在表達(dá)這個參數(shù)方面,胡榮祖等從反應(yīng)進(jìn)度和能量變化的關(guān)系,1熱分解過度到熱爆炸的充分必要條件2和非等溫反應(yīng)的動力學(xué)方程3導(dǎo)出了絕熱及近似絕熱條件下,速率方程為dα/dt=A1exp(-Ea1/RT)(1-α)m+A2exp(-Ea2/RT)αn(1-α)p的表觀經(jīng)驗級數(shù)和m=n=p=1的一級自催化分解反應(yīng)體系(dT/dt)Tb值的估算式;3,4王耘等5導(dǎo)出了絕熱和近似絕熱條件下反應(yīng)體系的dT/dt值的估算式,但對(dT/dt)Tb值的估算式未作研究.在估算這個參數(shù)方面,Hu等6-10估算了氮含量分別為11.92%、11.97%、13.54%、13.86%和14.14%的硝化棉(NC)的一級自催化分解反應(yīng)體系的(dT/dt)Tb值;Eisenreich等11,12估算了NC一級自催化分解反應(yīng)的動力學(xué)參數(shù).在實測這個參數(shù)方面,胡榮祖等13研究了硝仿熱爆炸的爆前加熱、熱分解溫升和時間的關(guān)系.本工作作為文獻(xiàn)3-12的拓展,報道了自催化反應(yīng)速率方程分別為:dα/dt=Aexp(-E/RT)α(1-α),dα/dt=Aexp(-E/RT)(1-α)n(1+Kcatα),dα/dt=Aexp(-E/RT)[αa-(1-α)n],dα/dt=A1exp(-Ea1/RT)(1-α)+A2exp(-Ea2/RT)α(1-α),dα/dt=A1exp(-Ea1/RT)(1-α)m+A2exp(-Ea2/RT)αn(1-α)p,dα/dt=Aexp(-E/RT)·(1-α),dα/dt=Aexp(-E/RT)(1-α)n,dα/dt=A1exp(-Ea1/RT)+A2exp(-Ea2/RT)(1-α)和 dα/dt=A1exp(-Ea1/RT)+A2exp(-Ea2/RT)α(1-α)的(dT/dt)Tb估算式的導(dǎo)出途徑和NC(13.54%N)自催化分解反應(yīng)動力學(xué)參數(shù)和自催化分解轉(zhuǎn)向熱爆炸時(dT/dt)Tb值的估算結(jié)果.
單位時間內(nèi)由于EMs熱分解而放出的熱量q1為
式中,Q為熱分解反應(yīng)的焓(J·mol-1),V為EMs的裝填體積(cm3),d為裝填密度(g·cm-3),M為EMs的摩爾質(zhì)量(g·mol-1),dα/dt為機理函數(shù)為f(α)=α(1-α)的自催化反應(yīng)(Au)速率:
聯(lián)立方程(1)和(2),得
與此同時,單位時間內(nèi)因傳熱由反應(yīng)區(qū)通過器壁向四周環(huán)境散失的熱量q2為
式中,k′為傳熱系數(shù)(J·cm-2·K-1·s-1),S表示藥柱表面積(cm2);Tc為按照線性關(guān)系Tc=T0+βt確定的反應(yīng)器壁和空間溫度;β表示線性加熱速率(K·min-1),T0表示熱分析曲線離開基線的溫度(K).
熱爆炸發(fā)生時,方程(3)變?yōu)?/p>
此處,αb是相應(yīng)于Tb的α值,Tb是EMs的熱爆炸溫度(K).
方程(4)變?yōu)榇颂?Te0是β→0時的onset溫度.
根據(jù)Semenov的熱爆炸理論,2熱分解過渡到熱爆炸的充分必要條件式為
此處(dT/dt)Tb是熱分解轉(zhuǎn)向熱爆炸時EMs中的臨界溫升速率.
方程(4)對T微分,得
方程(14)稱為Au自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得E、A、Te0、Tb和αb,就可從方程(14)得到(dT/dt)Tb的值.
類似地,由機理函數(shù)為 f(α)=(1-α)n(1+Kcatα))速率方程
方程(18)稱為CnB自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得E,A,n,Kcat,Te0,Tb和αb,就可從方程(18)得(dT/dt)Tb值.
由機理函數(shù)為 f(α)=αa-(1-α)n的自催化反應(yīng)(Bna)的速率方程方程(22)為Bna自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得E,A,Te0,n,a,Tb和αb,就可從方程(22)得(dT/dt)Tb值.
由一級自催化分解反應(yīng)速率方程
式中,k1=A1exp(-Ea1/RT),k2=A2exp(-Ea2/RT).
方程(26)稱為一級自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得Ea1,Ea2,A1,A2,Te0,Tb和αb,就可從方程(26)得(dT/dt)Tb值.
由經(jīng)驗級數(shù)自催化分解反應(yīng)速率方程
方程(30)稱為表觀經(jīng)驗級數(shù)自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得Ea1,Ea2,A1,A2,Te0,Tb,αb,m,n和p,就可從方程(30)得(dT/dt)Tb值.
由簡單一級自催化分解反應(yīng)速率方程
方程(34)為簡單一級自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得E、A、Te0和Tb,就可從方程(34)得(dT/dt)Tb值.
由簡單n級自催化分解反應(yīng)速率方程
方程(38)稱為簡單n級自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得E,A,Te0,n,Tb和αb,就可從方程(38)得(dT/dt)Tb值.
由經(jīng)驗級數(shù)m=0,n=0,p=1的自催化分解反應(yīng)速率方程
方程(40)稱表觀經(jīng)驗級數(shù)m=0,n=0,p=1的自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得Ea1,Ea2,A1,A2,Te0,Tb和αb,就可從方程(40)得(dT/dt)Tb值.
由經(jīng)驗級數(shù)m=0,n=1,p=1的自催化分解反應(yīng)速率方程
方程(42)稱為表觀經(jīng)驗級數(shù)m=0,n=1,p=1的自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率表達(dá)式.
一旦從熱流曲線解得Ea1,Ea2,A1,A2,Te0,Tb和αb,就可從方程(42)得(dT/dt)Tb值.
為求方程(2)和(31)中的2參數(shù)(A,E),方程(35)中的3參數(shù)(A,E,n),方程(15)、(19)、(23)、(39)和(41)中的4參數(shù)(分別為A,E,n,Kcat;A,E,a,n;A1,A2,Ea1,Ea2;A1,A2,Ea1,Ea2;A1,A2,Ea1,Ea2)和方程(27)中的7參數(shù)(A1,A2,Ea1,Ea2,m,n,p),進(jìn)行了線性最小二乘法和信賴域方法14的數(shù)值模擬.
對n=1的自催化反應(yīng)速率方程,有
由方程(43)和(44)的左端項與1/T作圖,用最小二乘法從斜率得E,截距得A.
經(jīng)非線性優(yōu)化模型:
對n級自催化反應(yīng)速率方程,有
對CnB速率方程,有
對Bna速率方程,有
對一級自催化分解反應(yīng)速率方程,有
對經(jīng)驗級數(shù)自催化分解反應(yīng)速率方程,有
對經(jīng)驗級數(shù)m=0,n=0,p=1的自催化分解反應(yīng)速率方程,有
對經(jīng)驗級數(shù)m=0,n=1,p=1的自催化分解反應(yīng)速率方程,有
由最小化均方誤差估計方程(45)中的3參數(shù)(A,E,n),方程(46)、(47)、(48)、(50)和(51)中的 4 參數(shù)(A,E,n,Kcat;A,E,a,n;A1,A2,Ea1,Ea2;A1,A2,Ea1,Ea2;A1,A2,Ea1,Ea2)和方程(49)中的7參數(shù)(A1,A2,Ea1,Ea2,m,n,p).
為準(zhǔn)確估算NC(13.54%N)放熱分解反應(yīng)的動力學(xué)參數(shù),進(jìn)行了邊界約束,限制參數(shù)范圍,將方程(45)-(51)分別改寫如下:
對n級自催化反應(yīng)速率方程,有
對CnB速率方程,有
對一級自催化分解反應(yīng)速率方程,有
對經(jīng)驗級數(shù)自催化分解反應(yīng)速率方程,有
對經(jīng)驗級數(shù)m=0,n=0,p=1的自催化分解反應(yīng)速率方程,有
對經(jīng)驗級數(shù)m=0,n=1,p=1的自催化分解反應(yīng)速率方程,有
為敘述簡便,將具有邊界約束的非線性優(yōu)化問題重新描述為
其中,對n級自催化反應(yīng)速率方程,有
對CnB速率方程,有
對Bna速率方程,有
對一級自催化分解反應(yīng)速率,有
u=(1018,1018,106,106)T,i=1,2,…,N
對經(jīng)驗級數(shù)自催化分解反應(yīng)速率方程,有
對經(jīng)驗級數(shù)m=0,n=0,p=1的自催化分解反應(yīng)速率方程,有
對經(jīng)驗級數(shù)m=0,n=1,p=1的自催化分解反應(yīng)速率方程,有
定義矩陣D(x)為D(x)? diag[|vi(x)|-1/2],即,D-2是一對角矩陣且其第i個對角元素為|vi(x)|.考慮非線性方程組
設(shè)xk∈int(Φ),式(59)的牛頓迭代滿足
其中,M(x)?B(x)+C(x),C(x)?D(x)diag[g(x)]Jv(x)D(x),
業(yè)已證明,xk是式(59)的局部最小值點,當(dāng)且僅當(dāng)sk是下式的解
這里Δk為正數(shù),表示信賴域的大小.因此,通過求解子問題(70)來確定xk是否為目標(biāo)函數(shù)的最小值點是合理的,進(jìn)而,當(dāng)信賴域大小Δk足夠大時,局部最小值的鄰域中式(67)的牛頓迭代正好是信賴域子問題(70)的解.
解決具有邊界約束的非線性優(yōu)化問題的信賴域方法模型為:
取x0∈int(Φ),循環(huán)k=0,1,…
(1)計算f(xk),gk,Hk及Ck;
(2)根據(jù)式(70)計算sk且滿足xk+sk∈int(Φ).
(3)計算
(5)修正模型yk,矩陣Dk及信賴域大小Δk.
信賴域大小Δk的修正為:給定0<μ<η<1,γ1<1<γ2及Λl>0,
依據(jù)表1和表2中NC(13.54%N)放熱分解反應(yīng)的T-轉(zhuǎn)化率(α)和β-Te關(guān)系數(shù)據(jù),用方程(14),(18),(22),(26),(30),(34),(38),(40)和(42)分析表1和表2中數(shù)據(jù),得表3中Au,CnB,Bna,一級/表觀經(jīng)驗級數(shù)、簡單一級和n級及表觀經(jīng)驗級數(shù)m=0,n=0,p=1和m=0,n=1,p=1的自催化分解反應(yīng)動力學(xué)參數(shù)和自催化分解反應(yīng)轉(zhuǎn)向熱爆炸時的臨界溫升速率(dT/dt)Tb值.據(jù)此,發(fā)現(xiàn)9個微分方程中Au速率方程[方程(2)]算得的E和A值在EMs分解反應(yīng)動力學(xué)參數(shù)正常范圍[E=80-250 kJ·mol-1,lg(A/s-1)=7-30]17外,簡單一級自催化分解反應(yīng)方程[方程(31)]求E和A的線性相關(guān)系數(shù)為0.9157,遠(yuǎn)離0.98,所得(dT/dt)Tb為負(fù)值,屬不合理,被排除.同時滿足E和A值在正常范圍內(nèi),(dT/dt)Tb為正值的方程有7個,在這7個方程中,以計算值與實驗值的相對誤差為判據(jù),宜做描述NC(13.54%N)自催化分解反應(yīng)過程的方程依次為:方程(27)、(23)、(41)、(39)、(19)、(35)、(15),對應(yīng)誤差依 次 為:0.0206、0.0255、0.0299、0.109、0.1698、0.1932、0.7620,相應(yīng)熱爆炸臨界溫升速率依次為:0.103、0.169、0.189、0.00002、0.615、0.020、0.266 K·s-1.由此認(rèn)為,用相對誤差最小的經(jīng)驗級數(shù)自催化反應(yīng)速率方程描述NC(13.54%N)熱分解過程是可取的,熱分解轉(zhuǎn)向熱爆炸時的臨界溫升速率為0.103 K·s-1是可接受的.
表1 用DSC測得的NC(13.54%N)的熱分解數(shù)據(jù)8Table 1 Thermal decomposition data8of NC(13.54%N)determined by DSC
表2 NC(13.54%N)熱爆炸臨界溫度(Tb)的計算值Table 2 Calculated values of the critical temperature(Tb)of thermal explosion for NC(13.54%N)
表3 用方程(14),(18),(22),(26),(30),(34),(38),(40)和(42)分析表1和表2中數(shù)據(jù)的結(jié)果Table 3 Results for analyzing the data in Tables 1 and 2 by Eqs.(14),(18),(22),(26),(30),34),(38),(40),and(42)
(1)提出了從不同升溫速率條件下的DSC曲線數(shù)據(jù)計算/確定EMs自催化反應(yīng)的動力學(xué)參數(shù)和自催化分解轉(zhuǎn)向熱爆炸時的臨界溫升速率(dT/dt)Tb的方法.
(2)非等溫條件下NC(13.54%N)熱分解過程可用表觀經(jīng)驗級數(shù)自催化分解反應(yīng)動力學(xué)方程描述:
(3)NC(13.54%N)自催化分解反應(yīng)過渡到熱爆炸時的臨界溫升速率值為0.103 K·s-1.
(1) Deng,Y.Chem.J.Chin.Univ.1985,6(7),621.[鄧 郁.高等學(xué)校化學(xué)學(xué)報,1985,6(7),621.]
(2) Semenov,N.N.On Some Problem of Kineties and Reactivity;Ind.,AN SSSR:Moscow,1958;p 421.
(3) Hu,R.Z.;Gao,S.L.;Zhao,F.Q.;Shi,Q.Z.;Zhang,T.L.;Zhang,J.J.Thermal Analysis Kinetics;Science Press:Beijing,2008;pp 322-334.[胡榮祖,高勝利,趙鳳起,史啟楨,張同來,張建軍.熱分析動力學(xué).北京:科學(xué)出版社,2008:322-334.]
(4) Hu,R.Z.;Zhang,H.;Xia,Z.M.;Guo,P.J.;Gao,S.L.;Shi,Q.Z.;Lu,G.E.;Jiang,J.Y.Chin.J.Energ.Mater.2003,11(3),130.[胡榮祖,張 海,夏志明,郭鵬江,高勝利,史啟楨,路桂娥,江勁勇.含能材料,2003,11(3),130.]
(5)Wang,Y.;Feng,C.G.;Zheng,X.Chin.J.Energ.Mater.2000,8(3),119.[王 耘,馮長根,鄭 曉.含能材料,2000,8(3),119.]
(6) Hu,R.Z.;Guo,P.J.;Song,J.R.;Zhang,H.;Xia,Z.M.;Ning,B.K.;Fang,Y.;Shi,Q.Z.;Liu,R.;Luo,G.E.;Jiang,J.Y.Chin.J.Explo.Prop.2003,26(2),53.
(7) Ning,B.K.;Hu,R.Z.;Zhang,H.,Xia,Z.M.;Guo,P.J.;Liu,R.;Luo,G.E.;Jiang,J.Y.Thermochim.Acta 2004,416(1-2),47.doi:10.1016/j.tca.2003.11.029
(8)Hu,R.Z.;Guo,P.J.;Gao,S.L.;Zhang,H.;Xia,Z.M.;Ning,B.K.;Fang,Y.;Shi,Q.Z.;Liu,R.Chinese Journal of Polymer Science 2003,21(3),285.
(9) Guo,P.J.;Hu,R.Z.;Zhang,H.,Xia,Z.M.;Song,J.R.;Gao,S.L.;Ning,B.K;Liu,R.;Lu,G.E.;Jiang,J.Y.Chem.Res.Chin.Univ.2004,20(2),163.
(10)Zhang,H.;Xia,Z.M.;Guo,P.J.;Hu,R.Z.;Gao,S.L.;Ning,B.K.;Fang,Y.;Shi,Q.Z.;Liu,R.Journal of Hazardous Materials 2002,94(3),205.doi:10.1016/S0304-3894(02)00118-8
(11) Eisenreich,N.Beitrag Zur Kinetic Thermischer Zersetzungseaktionen(ThermoanalytischeAuswertung Der Zersetzung Von Nitrocellulose).Ph.D.Dissertation,Technical University of Munich,Munich,1978.
(12) Eisenreich,N.;Pfeil,A.Thermochim.Acta 1983,61(1-2),13.doi:10.1016/0040-6031(83)80300-1
(13) Hu,R.Z.;Song,Q.C.;Xie,J.J.;Liang,Y.J.Explosion and Shock Waves 1987,7(3),217.[胡榮祖,松全才,謝俊杰,梁燕軍.爆炸與沖擊,1987,7(3),217.]
(14)Thomas,F.C.;Li,Y.Y.SIAM Journal on Optimization 1996,6(2),418.doi:10.1137/0806023
(15) Ozawa,T.Bulletin of the Chemical Society of Japan 1965,38(1),1881.
(16) Zhang,T.L.;Hu,R.Z.;Xie,Y.;Li F.P.Thermochim.Acta 1994,244(3),171.
(17)Hu,R.Z.;Yang,Z.Q.;Ling,Y.J.Thermochim.Acta 1988,123,135.doi:10.1016/0040-6031(88)80017-0