偏二甲肼與二氧化氮?dú)庀喾磻?yīng)的數(shù)值模擬研究

詹 翔，崔村燕，周宵燈，辛騰達(dá)，韓向陽

（航天工程大學(xué)，北京，101416）

0 引 言

目前，中國在航天領(lǐng)域使用的大型火箭所采用的液體推進(jìn)劑主要是肼類燃料（偏二甲肼UDMH）和硝基氧化劑（四氧化二氮 NTO），具有很強(qiáng)的毒性、腐蝕性、吸濕性和易燃、易爆、易揮發(fā)等特點[1]。由于推進(jìn)劑在貯存過程中采取了一系列防護(hù)措施，發(fā)生大規(guī)模泄漏的可能性較小，容易出現(xiàn)的泄漏問題主要是滲漏或滴漏。滲漏或滴漏過程中，偏二甲肼主要以蒸汽形式擴(kuò)散，四氧化二氮則分解為二氧化氮擴(kuò)散。當(dāng)兩者相遇時發(fā)生劇烈燃燒、爆炸，需要溫度和濃度條件如下[2]：

a）25 ℃以下，9%氣相偏二甲肼與 100%氣相二氧化氮；

b）50 ℃以下，61%氣相偏二甲肼蒸氣與 35%氣相二氧化氮接觸；

c）70 ℃以上，兩種推進(jìn)劑蒸汽接觸。

由此可見，偏二甲肼與二氧化氮同時泄漏后發(fā)生反應(yīng)不會立即爆炸，但兩者發(fā)生緩慢化學(xué)反應(yīng)所釋放的熱量有可能觸發(fā)爆炸。本文通過數(shù)值模擬的方式研究常規(guī)推進(jìn)劑間的氣相反應(yīng)規(guī)律，為相關(guān)事故處理措施提供依據(jù)。

1 反應(yīng)機(jī)理

Ernesto等[3]通過試驗研究給出了偏二甲肼與二氧化氮在無光照、室溫（15 ℃）、氣相條件下的反應(yīng)過程。所有反應(yīng)物和生成物均被傅里葉變換紅外光譜儀進(jìn)行檢測。試驗中各組分濃度隨時間的變化如表1所示。

表1 試驗中相關(guān)物質(zhì)在各時刻的濃度Tab.1 Concentration of the Related Substances at All Times in the Test

續(xù)表1

Ernesto等人根據(jù)試驗產(chǎn)物給出如下反應(yīng)路徑：

由試驗產(chǎn)物反應(yīng)可以看出，在室溫條件下偏二甲肼與四氧化二氮分解產(chǎn)生的二氧化氮?dú)庀喾磻?yīng)較為緩和，二氧化氮奪取偏二甲肼上與氮原子相連的氫原子，生成中間產(chǎn)物（B）以及亞硝酸（HONO）；在無光照和室溫條件下，中間產(chǎn)物（B）在二氧化氮的作用下進(jìn)一步脫去氮上的氫原子生成（H），此時中間體（H）上的氮原子相連接的氫全部被二氧化氮奪走，氮原子上存在兩個未配對的電子。當(dāng)兩個（H）相遇時，有未配對電子的氮原子就會配對形成氮氮雙鍵，形成最終產(chǎn)物（I）1，1，4，4-四甲基-2-四氮烯。結(jié)合反應(yīng)物消耗速率和生成物合成速率，給出如下反應(yīng)方程式：

其整體反應(yīng)速率系數(shù)由肼的衰變率確定為2.3×10-17cm3/（mol·s）。該試驗在無光照、室溫 15 ℃條件下進(jìn)行，是絕大多數(shù)液態(tài)推進(jìn)劑貯存的環(huán)境條件。因此該反應(yīng)機(jī)理可以代表偏二甲肼與四氧化二氮推進(jìn)劑泄漏時的氣相反應(yīng)過程。

基于FLUENT軟件模擬化學(xué)反應(yīng)時，放熱量是根據(jù)生成物與反應(yīng)物的生成焓之差計算得到的。1，1，4，4-四甲基-2-四氮烯（TMTZ）是一種含氮化合物，分子結(jié)構(gòu)如圖 2所示，常見于偏二甲肼和二氧化氮的反應(yīng)產(chǎn)物中[4～6]。為獲得 TMTZ分子的生成焓，采用量子化學(xué)的方法對TMTZ分子的生成焓進(jìn)行計算。

圖2 TMTZ分子結(jié)構(gòu)示意Fig.2 TMTZ Molecular Structure Diagram

2 TMTZ分子生成焓計算

2.1 計算方法

生成焓對于高能物質(zhì)而言是十分重要的一項，是評價物質(zhì)安全性和能量水平的標(biāo)準(zhǔn)之一[7]。目前量子化學(xué)方法能十分精確地計算高能化合物的能量，并且采用理論計算可以避免重復(fù)大規(guī)模試驗帶來的經(jīng)濟(jì)負(fù)擔(dān)以及試驗過程中可能發(fā)生的毒害和爆炸危險。通過量子化學(xué)的方法計算得到 TMTZ的生成焓既經(jīng)濟(jì)又安全。駱艷嬌等[8]對幾種不同多氮含能材料的生成焓進(jìn)行了計算，同時對比了幾種不同的密度泛函方法，研究過程表明：對于多氮含能材料，雜化密度泛函方法（DFT-B3LYP）計算出的生成焓可靠度高、占用資源少，計算速度快。本文采用的密度泛函理論方法計算生成焓，存在計算結(jié)果不直接給出生成焓，只能得到分子的一些熱力學(xué)參數(shù)，要借助輔助方法才能求得目標(biāo)化合物的生成焓。常用的輔助方法包括原子化反應(yīng)法[9～11]、高斯理論法等。

高斯理論法（Gaussian theory）是目前使用較為廣泛的理論計算化合物能量的方法之一。1989年John等人發(fā)表 Gaussian-1(簡稱 G1)方法以來，已經(jīng)有 G2、G3和G4等4代方法。在缺少試驗值的情況下，其計算結(jié)果有時也被當(dāng)作標(biāo)準(zhǔn)數(shù)值使用[12]。若分子體系中的非氫原子超過7個，高斯理論方法計算耗時將大大增加。而TMTZ分子中的非氫原子為7個，計算量將成指數(shù)級增加，因此不適宜使用高斯理論法。

基于上述分析，本文采用原子化反應(yīng)法對熱力學(xué)參數(shù)進(jìn)行處理以獲得TMTZ分子的生成焓。

2.2 生成焓計算

選用含有電子相關(guān)效應(yīng)的B3LYP雜化密度泛函關(guān)聯(lián) 6-31G(d)基組對 TMTZ分子結(jié)構(gòu)進(jìn)行幾何優(yōu)化計算，獲得了能量低、穩(wěn)定的構(gòu)型，結(jié)果如圖3所示。

圖3 TMTZ最穩(wěn)定分子構(gòu)型示意Fig.3 Stable Molecular Configuration of TMTZ

在穩(wěn)定構(gòu)型的基礎(chǔ)上進(jìn)行頻率計算，獲得焓的熱校正量Hcorr和零點矯正能εZPE；分別在B3LYP雜化密度泛函關(guān)聯(lián)6-31G(d)基組下計算TMTZ分子中H、C、N的單個原子總能量，并計算分子的總能量。將計算結(jié)果代入原子化反應(yīng)公式[13，14]，有：

式中 M為分子；X為M分子的組成原子；Δf(X，0K)為各組成原子生成焓；ε0(X)為各組成原子總能量之和；[ε0( M) + εZPE(M)]為M分子的單點能加上其零點矯正能；(298K )- Δ HΘ(0K)為M分子的生成焓溫度修正值且 f m

等于 Hcorr-εZPE、[( 298K)-(0K)]各原子的焓值溫度修正值之和。

表2 幾種氣態(tài)原子在0K時的生成焓值和生成焓矯正值[15]Tab.2 Enthalpy and Enthalpy of Formation of Atoms at 0K

基于量子化學(xué)方法計算出H、C、N原子的單點能ε0( X)、分子的焓值差( 298K)-(0K)、單點能ε0(X)以及零點能εZPE(M)之和，結(jié)果如表3所示。

表3 TMTZ分子計算結(jié)果Tab.3 Calculation Results of TMTZ

密度泛函理論是目前應(yīng)用較為廣泛的量子化學(xué)方法，采用B3LYP雜化密度泛函計算多氮含能物質(zhì)生成焓的誤差平均為[8]18.84 kJ/mol。對于TMTZ分子來說，理論誤差約占計算值的5%，占比較小。該計算結(jié)果具有一定的可靠度，可用來計算化學(xué)反應(yīng)的放熱量，分析反應(yīng)過程中溫度的變化趨勢。

3 數(shù)值模擬過程

3.1 計算域和基本設(shè)置

氣相反應(yīng)試驗在1.2 m×2.4 m×2.4 m的箱式密閉容器中進(jìn)行[3]。計算域根據(jù)試驗采用2.4 m×2.4 m的2D模型，四周采用壁面條件。反應(yīng)條件為15 ℃、一個標(biāo)準(zhǔn)大氣壓。假設(shè)反應(yīng)物通入到容器中后立即完全混合，在數(shù)值模擬初始化時直接根據(jù)表1設(shè)置反應(yīng)物的預(yù)混濃度。反應(yīng)的背景填充物為空氣。

3.2 數(shù)值模擬與試驗結(jié)果對比

設(shè)置監(jiān)測面記錄反應(yīng)物、生成物在計算域中的平均濃度變化，計算結(jié)果如圖4所示。從圖4中可以看出，反應(yīng)速率根據(jù)反應(yīng)物濃度經(jīng)歷了一個由快到慢的過程，開始時，由于反應(yīng)物濃度較大，反應(yīng)速率較快，當(dāng)反應(yīng)進(jìn)行到5 min左右時，反應(yīng)物的濃度減少導(dǎo)致反應(yīng)速率下降，生成物濃度幾乎不再變化。根據(jù)計算結(jié)果，生成物、反應(yīng)物的增減關(guān)系嚴(yán)格遵循式（1）給出的化學(xué)計量數(shù)。

圖4 各物質(zhì)摩爾濃度隨時間變化曲線Fig.4 Molar Concentration of Each Component with Time

各組分試驗值和數(shù)值模擬結(jié)果對比如圖5所示。

續(xù)圖5

從圖5中可以看出，數(shù)值模擬結(jié)果與試驗值在變化趨勢上基本保存一致。為了定量地研究數(shù)值模擬結(jié)果與試驗結(jié)果的差異，采用Hanna等人提出的誤差估計方法評價的模擬效果[16，17]。該方法的誤差分析包括F（平均分?jǐn)?shù)偏差）和M（幾何平均偏差）、V（幾何平均方差）、N（標(biāo)準(zhǔn)化均方差）和FAC2（模擬值與觀測值之比），分別表示為

式中 Xo為試驗觀測值；Xp為模擬值。

該誤差統(tǒng)計分析法可以較為全面地比較模擬結(jié)果與試驗值的差異[18]。通過以上分析方法對數(shù)值模擬結(jié)果進(jìn)行統(tǒng)計，結(jié)果如表4所示。

表4 各組分?jǐn)?shù)值模擬結(jié)算誤差分析Tab.4 Numerical Simulation’s Deviation of Each Component

對于可靠性較高的數(shù)值模擬結(jié)果，其誤差分析得到的各參數(shù)應(yīng)滿足[19]：-0.67＜F＜0.67，N＜4，F(xiàn)AC2＞0.5。由表4可知，數(shù)值模擬結(jié)果的誤差均在合理的區(qū)間內(nèi)，數(shù)值仿真較好地復(fù)現(xiàn)了試驗過程。從圖5對比中分析，對于反應(yīng)物來說，數(shù)值模擬計算的結(jié)果總是略高于試驗值，說明在實際反應(yīng)過程中，還有少量的UDMH和NO2轉(zhuǎn)換成了其它物質(zhì)。由文獻(xiàn)[20]的分析可知，生成物可能是二甲基亞硝銨(CH3CH3N-N=O)、二甲胺（CH3CH3N-H）、甲醇（CH3OH）、二甲基甲酰胺（CH3CH3N-CH=O）和二甲基肼硝酸鹽（CH3CH3N-NH3NO3）等甲基衍生物。從偏二甲肼的濃度誤差來看，該轉(zhuǎn)化量比較少，在誤差可接受的范圍內(nèi)。結(jié)合生成物濃度情況，HONO的濃度值總是小于試驗值，說明生成物HONO還有其它合成來源，可能是NO2與空氣中的組分發(fā)生了一系列反應(yīng)生成的額外的亞硝酸HONO。

設(shè)置監(jiān)測點對反應(yīng)區(qū)域溫度進(jìn)行檢測，取整個區(qū)域的平均溫度，結(jié)果以288 K為零點，如圖6所示。

圖6 溫度隨時間變化曲線Fig.6 Temperature Along with the Time Variation

由圖 6可知，在試驗給出的反應(yīng)物濃度下，反應(yīng)放熱使氣體升溫的效果不明顯。數(shù)值模擬結(jié)果顯示：溫度增長趨勢與反應(yīng)速率增長趨勢一致，反應(yīng)釋放的熱量使整個區(qū)域溫度上升約0.072 K。但是在實際泄漏的情況下，由于偏二甲肼與二氧化氮均比空氣重，兩推進(jìn)劑蒸汽受到重力的影響在地勢較低的地方聚集。當(dāng)積累的濃度較高時，反應(yīng)釋放的大量熱會使推進(jìn)劑間的反應(yīng)類型由冷反應(yīng)階段（該反應(yīng)階段默認(rèn)反應(yīng)溫度在室溫 15 ℃左右）轉(zhuǎn)變?yōu)闊岱纸怆A段[21]。在熱分解階段，一些大分子物質(zhì)可能受到高溫的影響發(fā)生熱分解反應(yīng)放出大量熱，進(jìn)一步加速熱分解反應(yīng)，使推進(jìn)劑蒸汽發(fā)生劇烈燃燒、爆炸。

4 結(jié) 論

a）采用DFT-B3LYP方法計算出TMTZ分子的熱力學(xué)數(shù)據(jù)，經(jīng)原子化反應(yīng)法處理后得到標(biāo)準(zhǔn)生成焓，為-369.15 kJ/mol；

b）數(shù)值模擬結(jié)果與試驗值基本一致，反應(yīng)速率增長趨勢相同，數(shù)值仿真能較好地反映生成物濃度變化；

c）計算結(jié)果與實驗結(jié)果存在偏差的原因是：少量偏二甲肼與二氧化氮反應(yīng)生成了其它甲基衍生物，同時二氧化氮與空氣中的組分反應(yīng)生成了額外的亞硝酸；

d）反應(yīng)釋放的熱量使整個區(qū)域溫度上升了約0.072 K。