HTPB/Al/AP/RDX推進劑初始燃燒的分子模擬

初慶釗,付小龍,鄭學明,劉金龍,陳東平

(1.北京理工大學 爆炸科學與技術國家重點實驗室,北京 100081;2.西安近代化學研究所,陜西 西安 710065;3.黑龍江北方工具有限公司,黑龍江 牡丹江 157013)

引 言

HTPB推進劑是武器的重要動力來源,其不僅具有較寬的燃速調節范圍且力學性能好,制造工藝簡單,在國內外多種火箭發動機型號中得到應用[1,2]。HTPB推進劑成分包含不同級配的AP、RDX、Al顆粒物以及HTPB黏合劑等物質,是一種非均質復合含能材料。

為了研究HTPB推進劑的燃燒理論模型,Beckstead等[3,4]提出了BDP多火焰模型,將HTPB推進劑燃燒分為3個階段,使用簡化反應機理描述AP/HTPB的反應過程。顏密[5]采用BDP 微尺度燃燒模型對AP/HTPB復合推進劑進行模擬,并進行穩態計算和燃速驗證,研究了工作壓強、AP 濃度和 AP 粒徑這3個參數對組分擴散特性、氣相溫度分布以及火焰位置的影響,為推進劑壓力耦合響應函數的數值計算提供理論基礎。張晏榕等[6]提出了一套氣固耦合的三維空間求解數值框架,使用兩套不同網格描述氣相和固相演化,實現了三維非均質推進劑的燃燒模擬。目前的燃燒模擬研究大多關注AP/HTPB兩種組分,對于HTPB/Al/AP/RDX四組元推進劑,其界面存在復雜的相互作用機制,相關基礎研究缺失,導致難以開展燃燒模擬研究。

推進劑配方直接影響推進劑的物理化學性質,不同組分之間的界面上存在著復雜的相互作用,導致其燃燒、爆炸、安全性能等發生本質改變。國內外學者針對四組元HTPB推進劑的燃燒特性開展了基礎實驗研究。黃蒙等[7]對HTPB/Al/AP/RDX推進劑組分之間的相互作用開展研究,通過熱分解動力學分析,發現在熱分解階段,組分HTPB、AP和RDX之間存在明顯的相互作用。其中, AP會強烈加速HTPB和RDX的分解;反之,HTPB和RDX也會促進AP分解。吳迎春等[8]采用全息技術對HTPB/Al/AP推進劑燃燒過程中的顆粒行為進行可視化分析,發現鋁顆粒燃燒過程中經歷了剝落、微爆、燃燒等過程。盡管前人針對HTPB推進劑已開展了系列實驗研究,對于推進劑組分間影響規律有一定的認識,然而推進劑組分間的相互作用往往發生于微納米尺度,受限于實驗手段和設備,尚未有合理的實驗方法可直接測試其微觀反應機制,這導致HTPB推進劑微觀反應機制缺失,限制了其性能的優化提升。

分子動力學模擬可以深入探究微觀物理化學機制,是研究復雜燃燒過程的重要工具,在復雜體系反應模擬中得到了廣泛的應用[9-13]。其往往依賴于力場或勢函數模型,目前含能材料領域中最著名的力場是van Duin 和Goddard提出的ReaxFF可反應力場[14],該方法通過考慮原子間鍵級來描述化學反應過程,已被應用于研究鋁粉[11]、高能炸藥[15]的燃燒與爆炸機理。在HTPB推進劑中,AP是重要的氧化劑組分。但是,目前還未有適用于AP反應體系的ReaxFF力場參數,這限制了AP熱解與燃燒反應機理的研究。近期,本課題組針對AP單質開發了機器學習勢函數[16],研究了AP熱分解過程中的詳細反應機理,觀察到了AP受熱過程中相態演化、吸/放熱過程以及最終的氣相產物。雖然,在過去20年間分子動力學力場在含能材料領域中取得了重大的成功,但是目前已公開的勢函數模型通常僅能描述推進劑中1～2種組分,這無法滿足推進劑全配方分子模擬的現實需求。推進劑全配方分子模擬一直是世界性的難題,不僅需要描述推進劑單質組分,還需要考慮組分間復雜的相互作用,目前國際上未有成功先例。如若能建立推進劑全配方的勢函數模型,有望開展真實推進劑配方燃燒、力學等物理化學過程的全原子模擬,實現推進劑的理性設計。

基于以上分析,本研究針對HTPB/Al/AP/RDX等推進劑關鍵組分開發了一個通用勢函數。該勢函數基于第一性原理計算的數據集,使用深度神經網絡模型進行開發。針對推進劑組分的能量、受力開展了測試驗證。隨后,建立了HTPB四組元推進劑燃面模型,并進行了大規模的分子動力學模擬計算。在模擬計算中,成功觀察到了推進劑燃燒的相變、熱解、傳熱、傳質等微觀過程。在世界范圍內首次實現了HTPB四組元推進劑的燃燒微觀模擬,揭示了其中的界面作用機制,為固體推進劑的理論建模提供了新的思路。

1 模擬方法

1.1 勢函數開發

分子動力學模擬的核心是反應勢函數模型,本研究采用深度勢能框架(DeepMD)中的神經網絡模型來擬合勢函數參數[17]。圖1為機器學習勢函數開發流程。

圖1 機器學習勢函數開發流程Fig.1 Development process of machine learning potential function

在圖1中,訓練數據集由第一性原理MD模擬的軌跡組成。然后,建立了一個神經網絡勢函數(NNP)模型,該模型將原子配位(R)解釋為原子間力(F)和能量(E),從而實現了在大空間和時間尺度上具有第一性原理精度的MD模擬。模型中每個結構的能量由結構內每個原子能量之和表征:

E=∑Ei

(1)

式中:E為結構能量;Ei為原子i的能量,采用公式(2)計算:

Ei=Ewɑi(Ri)

(2)

式中:αi表示原子i的化學元素;wαi表示原子i的神經網絡參數;Ri為原子的環境矩陣。

(3)

式中:rij表示原子i、j間的笛卡爾坐標距離。

隨后建立深度神經網絡對結構的能量與受力進行預測。深度神經網絡包含一個三層(25、50、100節點/層)的濾波(嵌入)網絡和一個三層(240節點/層)的擬合網絡。與經典神經網絡相似,DeepMD方案通過使用反向傳播算法計算損失函數的梯度來訓練模型。NNP訓練共經過4.0×106次迭代,學習率從1.0×10-3到5.0×10-8呈指數衰減。訓練過程中的損失函數L為:

訓練集包括HTPB/Al/AP/RDX等組分的單質以及任意兩兩組分間的界面結構。在300、1000、2000、3000和4000K的溫度下,進行了一組1ps的第一性原理MD模擬,以獲得NVT系綜下的運動軌跡。第一性原理MD模擬使用CP2K[18]直接計算,計算方法采用Goedecker-Teter-Hutter (GTH)贗勢和Perdew-Burke-Ernzerhof (PBE)廣義梯度近似方法處理核心電子,基組考慮了雙ζ高斯基函數加極化(DZVP-MOLOPT)。其他詳細的設置請參考本課題組關于NNP開發的相應工作[19-21]。

1.2 四組元推進劑模型

圖2為構建的HTPB/Al/AP/RDX四組元推進劑裝填模型。

圖2 HTPB/Al/AP/RDX四組元推進劑裝填模型(圖中的 C/H/N/O/Al/Cl 分別使用灰/白/藍/紅/銀/綠色表示)Fig.2 HTPB/Al/AP/RDX four-component propellant model (The C/H/N/O/Al/Cl in the diagram are represented in gray/white/blue/red/silver/green respectively)

其中鋁粉、AP、RDX均為顆粒物,它們之間使用HTPB進行填充。鋁粉直徑5nm,包含1nm厚度的氧化層;AP直徑為20nm,RDX顆粒直徑為5nm,兩者均通過切割晶體結構獲得。計算域大小為2nm×25nm×100nm,共計126165個原子。計算域的x和y方向設為周期性邊界,z方向上設為反射邊界。此外,推進劑模型下方5?內原子均被固定住。為了模擬推進劑表面的燃燒過程,對模型的上表面區域(z>23.5nm)施加4000K高溫熱浴的正則系綜,采用Nose-Hoover恒溫器,阻尼因子設為20fs。其他區域采用微正則系綜(NVE),初始溫度設為500K,采用速度Verlet法對運動方程進行積分。分子動力學模擬計算使用LAMMPS軟件[22],時間步長為0.2fs,共模擬40ps,在NVIDIA V100 GPU上運行約24h。

2 結果與討論

2.1 機器學習勢函數驗證

用推進劑組分的第一性原理計算數據庫對勢函數進行測試,以評價勢函數的準確性,結果見圖3和表1。

數據庫中的分子結構包括HTPB四組元推進劑(Al/RDX/AP/HTPB)的各個單質與兩兩界面構型。所有模型均使用CP2K在300～4000K下進行采樣,每個溫度區間各有1000個構型。

如圖3所示,最終的四組元模型對于推進劑單質組分預測精度很高,對于不同組分組成的界面結構,誤差略高于單質體系。表1中展示了每個體系的能量與受力誤差。總體上,所有組分的MAE(平均絕對誤差)均小于40meV,受力誤差小于0.5eV/?,這表明四組元勢函數可較好地描述推進劑單質組分和界面上的化學反應過程。

圖3 單組分和兩兩界面的機器學習勢函數的能量與受力精度測試Fig.3 Energy and force accuracy testing of machine learning potential for single-component and interfaces

表1 機器學習勢函數對推進劑單質和界面組分的能量與受力精度Table 1 Prediction accuracy of the machine learning potential function for energies and forces on single and interfacial components of solid propellants

2.2 推進劑復合模型燃燒結構演化

采用新開發的勢函數模型對HTPB推進劑的燃燒過程進行全配方分子動力學模擬。對推進劑表面1.5nm區域內的原子施加4000K高溫,其余部分溫度為500K,以模擬推進劑燃面燃燒模型,燃燒過程如圖4所示,初始時刻(t=2ps)推進劑表面的HTPB與RDX先發生分解反應,釋放出氣體小分子。隨后在t=14ps時,HTPB進一步分解,氣體小分子產物向燃面上方擴散,AP顆粒上方部分也發生分解。而鋁粉反應較弱,在HTPB分解氣體的推動下逐漸脫離燃面。在t=22ps時,AP上方形成錐形分解產物區,同時AP/HTPB界面上也觀察到明顯的分層現象。當t=34ps時,推進劑整體上密度顯著降低,完全擴散至計算域中,生成氣體分解產物。

圖4 推進劑燃燒過程中結構演化Fig.4 Structural evolution during propellant combustion

進一步分析分子動力學模擬軌跡,可發現推進劑上方形成了擴散火焰結構。推進劑擴散火焰的局部放大圖如圖5所示,通過觀察原子軌跡觀察到此時RDX在高溫區加速的條件下發生了快速分解,生成了大量NO2、H2O、N2O等產物。而AP顆粒則分解產生NO、HCl、O2等氧化性氣體組分,這些組分可作為氧化劑與HTPB分解產生的碳氫小分子進行燃燒反應,主要反應區位于AP與HTPB的交界面處。鋁粉雖然同樣暴露在高溫區,但是此時并沒有直接參與燃燒反應,而是隨著AP和HTPB分解產生的小分子氣流向上運動。通過分子動力學模擬,可定量地對整個推進劑模型進行時間分辨的三維重建,進而獲得推進劑燃燒的微觀機理。

圖5 推進劑擴散火焰微結構局部放大圖(t=34ps)Fig.5 Local magnification of diffusion flame microstructure (t=34ps)

2.3 溫度演化

此外,推進劑燃燒過程中的溫度演化過程也可被完整地從分子模擬中重建出來。圖6展示了不同時刻推進劑溫度演化過程。

圖6 推進劑燃燒過程中溫度演化Fig.6 Temperature evolution during propellant combustion

初始時刻表面為4000K高溫,其余部分為500K。隨著熱量的傳遞,推進劑上方溫度快速升高。當t=2ps時,推進劑上方維持3000K左右的高溫,RDX和AP顆粒由于快速分解反應的吸熱效應,其所在區域溫度顯著降低。在AP周圍與HTPB交界處,也觀察到了顯著的高溫區,表明兩者在界面上發生了小分子產物之間的放熱反應。隨著整體燃燒過程的推進,高溫區逐漸從AP/HTPB交界處向兩側擴散,對AP/HTPB進行加熱,促進兩者的熱分解過程。由于鋁顆粒具有較好的熱導率,其溫度通過熱傳導,一直保持在較高的溫度(約3000K)與擴散火焰中的溫度基本保持一致。相對應的,AP顆粒內部的溫度在模擬過程中緩慢的從約500K升高至約1000K。

2.4 推進劑燃燒中組分演化

由圖7可知,除去部分熱點區域,HClO4生成的位置均勻分布在整個AP顆粒內部,而并非是從AP/HTPB界面處向顆粒內部擴散。從時間演化上看,在初始時刻HClO4濃度較高(約1.5～2mol/cm3),隨著反應的進行,HClO4濃度逐漸降至約0.5mol/cm3。這說明HClO4是整個燃燒過程中的一種中間產物,在后續反應中被快速消耗。

AP經過初始分解后,通過進一步的演化最終產生O2和Cl2,見圖8,它們將作為四組元HTPB推進劑中最重要的兩種氧化劑參與推進劑的燃燒過程。

圖8 O2/Cl2組分演化(圖中濃度為兩組分總濃度)Fig.8 Evolution of O2/Cl2 components(The concentration in the figure represents the total concentration of two components)

如圖8所示,和HClO4相比,O2/Cl2的生成具有顯著的空間依賴性,在AP/HTPB交界面處濃度較高,說明界面處AP與HTPB間的相互作用實際促進了O2/Cl2產生。這也為實驗中觀察到的HTPB對AP分解的促進作用提供了理論解釋[7]。在14ps時,AP/HTPB界面處氧化劑濃度進一步升高,同時在推進劑上方產生了微弱的尖錐狀氣相產物區域,該部分O2和Cl2的濃度要比推進劑下方區域的濃度低2～3倍。這主要由于尖錐狀區域的溫度顯著高于推進劑下方區域(圖6),導致HTPB分解產物和O2/Cl2在該區域發生了快速的氧化反應,釋放出了大量的能量。同時,這一尖錐狀區域與圖5所示的推進劑擴散火焰微結構一致,說明擴散火焰中主要發生的是HTPB分解的小分子碳氫化合物與O2/Cl2之間的氣相反應。

在推進劑中也添加了少量的RDX顆粒。在此前研究中發現RDX熱解的主要途徑是通過N—N鍵均裂發生的,產生NO2和RDR分子[21]。隨后NO2可通過進一步反應轉化為NO。圖9展示了NO/NO2的空間分布,可觀察到在2ps時RDX的位置產生了大量的NO/NO2分子,同時在AP顆粒周圍也有較少的分布。在推進劑組分的DSC實驗中,前人觀察到RDX的DSC峰溫為242.7℃,低于AP的峰溫296.7℃,這表明RDX更容易發生分解,與本研究的模擬結果一致。從濃度分布上看,RDX分解產生的NO/NO2分子濃度較高(大于0.4mol/cm3),而AP分解產生的NO/NO2分子濃度相對較低(小于0.15mol/cm3)。這些NO/NO2分子都將作為氧化劑與HTPB分解產物發生燃燒反應,最終轉化為N2和H2O等物質。

圖9 分解產物 NO/NO2分子的時空演化(圖中濃度為兩組分總濃度)Fig.9 The spatiotemporal evolution of NO/NO2 molecules(The concentration in the figure represents the total concentration of two components)

推進劑燃燒過程中,HTPB可分解產生一系列碳氫小分子,這些小分子可與AP分解產物形成擴散火焰,生成CO、CO2、H2O等燃燒產物。因為AP中不含C元素,所以通過CO/CO2的空間分布可具體確定擴散火焰發展的過程和具體形貌。圖10展示了CO/CO2的時空演化過程,初始時刻(2ps)CO/CO2分子主要分布在AP/HTPB交界處,隨著反應的進行,在推進劑AP顆粒上方逐漸形成了錐形的燃燒產物區。產物不斷向上擴散,最終充滿整個計算域。然而,在這個過程中,RDX顆粒附近并沒有看到較多CO/CO2。如果比較RDX顆粒的空間位置(圖9),能更清楚地發現RDX分解產物中沒有CO/CO2。前期關于RDX的熱分解過程研究發現RDX中CO/CO2的產生存在一定的滯后[21]。在2500K下觀察到RDX先分解產生NO、NO2、HNO2等物質,在10ps左右后才觀察到CO/CO2的生成。在推進劑燃燒中,推進劑表面的RDX先產生氮氧化物,RDX生成CO/CO2主要發生在氣相反應區。

圖10 CO/CO2組分演化(圖中濃度為兩組分總濃度)Fig.10 Evolution of CO/CO2 components(The concentration in the figure represents the total concentration of two components)

圖11 H2O組分演化Fig.11 Evolution of H2O components

相對于HTPB/AP/RDX,鋁粉與各個組分之間的反應相對較慢(圖12)。隨著燃燒的進行,可觀察到鋁粉顆粒快速向遠離燃面的方向運動,但是反應程度有限。從組分分布來看,鋁顆粒周圍是3000K高溫的H2O/CO2/O2/Cl2等氧化性氣體,但是在本工作中,未發現這些氧化性氣體與鋁粉之間的反應。另外,從顆粒溫度來看(圖6),鋁顆粒由于良好的熱傳導能力,在燃燒初期就已經達到了較高的溫度,但是由于外層的氧化層存在,其反應活性被嚴重制約。即使本工作中為了加速潛在的鋁粉與高溫氣體之間的反應,使用了1nm的超薄氧化層,仍然未發現鋁粉參與反應的過程。在推進劑燃燒過程中,AP、HTPB等組分反應產生了大量氣體,局部壓力升高,推動鋁顆粒離開推進劑表面,這一現象也在近期的推進劑燃燒實驗的全息測量中得到了證實[8]。

圖12 Al組分演化過程Fig.12 Evolution of Al components

3 結 論

(1)基于第一性原理計算的訓練集,開發了一個通用勢函數,可以開展HTPB/Al/AP/RDX推進劑關鍵組分的大規模分子動力學模擬,以研究推進劑的熱解和燃燒。

(2)建立了HTPB四組元推進劑燃面模型,并開展分子動力學模擬研究推進劑的燃燒過程,實現了對溫度、組分等信息的統計分析,發現了AP/HTPB間相互作用機制,AP分解產生的氧化劑在界面處與HTPB反應放熱,進而促進了AP的分解過程。

(3)觀察到了推進劑燃燒過程中鋁粉的剝離過程,推進劑燃燒過程中產生大量氣體,推動鋁粉離開推進劑表面,而鋁粉本身反應比較微弱。

(4)證明了分子動力學模擬能夠在原子尺度上實現時間分辨的三維重建,進而獲得推進劑燃燒的微觀機理,為固體推進劑的理論研究提供了新的工具。