朱 偉，劉冬梅，肖繼軍，池旭輝，龐愛(ài)民，肖鶴鳴

(1.南京理工大學(xué) 化工學(xué)院，分子與材料計(jì)算研究所，南京 210094； 2.嘉興學(xué)院 生物與化學(xué)工程學(xué)院，嘉興 314001；3.中國(guó)航天科技集團(tuán)公司四院四十二所，襄陽(yáng) 441003)

?

NEPE推進(jìn)劑/襯層結(jié)構(gòu)-性能MD模擬(Ⅱ)
——復(fù)雜體系組分分子遷移和配方設(shè)計(jì)示例①

朱 偉1，2，劉冬梅1，肖繼軍1，池旭輝3，龐愛(ài)民3，肖鶴鳴1

(1.南京理工大學(xué) 化工學(xué)院，分子與材料計(jì)算研究所，南京 210094； 2.嘉興學(xué)院 生物與化學(xué)工程學(xué)院，嘉興 314001；3.中國(guó)航天科技集團(tuán)公司四院四十二所，襄陽(yáng) 441003)

用分子動(dòng)力學(xué)(MD)方法，對(duì)(PEG/NG/BTTN)/NPBA/HMX/AP/PEG/N-100//HTPB/TDI復(fù)雜的推進(jìn)劑/襯層模型體系進(jìn)行295 K-NVT模擬研究，展示了組分分子的濃度分布和遷移狀況，發(fā)現(xiàn)HMX和NPBA分子有向界面層遷移趨勢(shì)，而AP則呈平均分布態(tài)勢(shì)。以RDX等量取代HMX后所得新配方的MD模擬研究表明，前者拉伸模量(E)、體模量(K)和剪切模量(G)、柯西壓(C12-C44)和K/G值均有明顯下降，表明新配方的剛性、強(qiáng)度和延展性均有下降；新配方中引發(fā)鍵(N—NO2)最大鍵長(zhǎng)(1.528 ?)明顯大于原配方中相應(yīng)值(1.503 ?)，預(yù)示新配方感度增大、安全性將下降；比較RDX、HMX與其他組分之間的結(jié)合能，前者小于后者，預(yù)示新配方的相容性較差。

推進(jìn)劑/襯層界面；組分遷移；力學(xué)性能；感度；相容性；分子動(dòng)力學(xué)模擬

0 引言

相容性、安全性和和力學(xué)性能是推進(jìn)劑、襯層及其界面所需重要性能，直接影響高能固體發(fā)動(dòng)機(jī)的工作效能和壽命長(zhǎng)短。這些性能一般依靠實(shí)驗(yàn)測(cè)定，如能通過(guò)理論模擬和計(jì)算加以定量預(yù)測(cè)或定性評(píng)估，則不僅可節(jié)省人力、物力和財(cái)力，縮短研究周期，而且可指導(dǎo)配方設(shè)計(jì)，是重大基礎(chǔ)研究所迫切需要解決的重要課題。

近年來(lái)興起的分子動(dòng)力學(xué)(MD)模擬，通過(guò)對(duì)龐大復(fù)雜體系結(jié)構(gòu)和性能的模擬計(jì)算，正好滿(mǎn)足為高能復(fù)合材料的配方設(shè)計(jì)提供豐富的信息、規(guī)律和理論指導(dǎo)的需求。通過(guò)MD模擬研究固體推進(jìn)劑和襯層的工作已有一些，研究大都集中在推進(jìn)劑中粘接劑與增塑劑的相容性和力學(xué)性能[1-8]，推進(jìn)劑中可遷移組分的影響作用[9-10]及推進(jìn)劑顆粒填充模型建立[11]。但這些工作主要集中于雙組分體系，個(gè)別多組分體系也只限制在推進(jìn)劑中的粘結(jié)與增塑體系中的簡(jiǎn)單模型[7]。因?qū)嶋H應(yīng)用的固體推進(jìn)劑均為多組分復(fù)雜體系，而對(duì)于推進(jìn)劑/襯層復(fù)雜體系的MD模擬至今未見(jiàn)公開(kāi)報(bào)道。

本文選擇常用粘接劑聚乙二醇(PEG)、增塑劑硝化甘油(NG)和1，2，4-丁三醇三硝酸酯(BTTN)、聚合物鍵合劑(NPBA)、氧化劑奧克托今(HMX)和高氯酸銨(AP)以及推進(jìn)劑和襯層固化體系PEG(聚乙二醇)/N-100(異氰酸酯)//HTPB(端羥基聚丁二烯)/TDI(甲苯二異氰酸酯)混合體系，即以推進(jìn)劑/襯層(PEG/NG/BTTN)/NPBA/HMX/AP/PEG/N-100//HTPB/TDI復(fù)雜體系為研究對(duì)象，搭建合理模型進(jìn)行MD模擬，首次展示典型組分分子向界面遷移的有趣圖像，并以RDX等量取代HMX，比較2種配方的相容性、力學(xué)性能和安全性(感度)，提供MD模擬用于配方設(shè)計(jì)的范例。

1 模型構(gòu)建和模擬細(xì)節(jié)

對(duì)于(PEG/NG/BTTN)/NPBA/HMX/AP/PEG/N-100//HTPB/TDI推進(jìn)劑/襯層十分復(fù)雜的混合體系，按一定配比和添加次序，由簡(jiǎn)到繁加以構(gòu)建。

1.1 (PEG/NG/BTTN)/NPBA/HMX/AP混合體系模型搭建與MD模擬

在(PEG/NG/BTTN)/NPBA/HMX/AP體系中，共含10 015個(gè)原子，依據(jù)其配比和先后次序要求加以混合。即按照實(shí)驗(yàn)要求，先將NG和BTTN以1∶1的比例混合搭好模型，然后與PEG混合組成粘合劑；再將構(gòu)建好的粘合劑體系與NPBA、HMX和AP分子，均勻地放入邊長(zhǎng)為200 ?×200 ?×200 ?的周期箱中。在修正和擴(kuò)展的MPCFF力場(chǎng)[12]和NVT系綜下進(jìn)行20 ps MD模擬，以使體系達(dá)到平衡；然后，緩慢縮小周期箱體積，同時(shí)進(jìn)行MD模擬使達(dá)到平衡；重復(fù)此過(guò)程，直到體系密度接近其理論值。取其最終結(jié)構(gòu)，溫度設(shè)定為295 K，再進(jìn)行200 ps MD周期性模擬，直至獲得平衡結(jié)構(gòu)。圖1(a)和(b)分別示出(PEG/NG/BTTN)/NPBA/HMX/AP混合體系的初始結(jié)構(gòu)和平衡結(jié)構(gòu)。

(a)初始結(jié)構(gòu) (b)平衡結(jié)構(gòu)圖1 (PEG/NG/BTTN)/NPBA/HMX/AP 混合體系初始結(jié)構(gòu)和平衡結(jié)構(gòu)

1.2 PEG/N-100和HTPB/TDI固化體系模型搭建與MD模擬

將2條端基飽和鏈節(jié)數(shù)為20的PEG高分子鏈與N-100構(gòu)成的混合模型(如圖2所示)，先進(jìn)行MM優(yōu)化，然后在295 K-NVT系綜下，進(jìn)行200 ps的MD模擬，所得結(jié)構(gòu)用作研究固化反應(yīng)的初始模型。對(duì)此模型，先在323 K-NPT系綜下進(jìn)行100 ps的MD模擬，然后使用熱力學(xué)方法、手動(dòng)連接鄰近的兩分子—NCO與—OH活性基團(tuán)，完成反應(yīng)[13]。對(duì)最后所得混合體系，在295 K-NVT系綜下，再進(jìn)行200 ps MD模擬。體系中含430個(gè)原子，PEG和N-100的質(zhì)量分?jǐn)?shù)分別為64.1%和35.9%。

圖2 N-100和PEG混合模型

仿照PEG/N-100固化體系搭建方式與模擬方法，HTPB/TDI固化體系也進(jìn)行類(lèi)似的搭建與模擬。該固化體系中含有5條鏈節(jié)數(shù)為20端基飽和的HTPB分子鏈和5個(gè)TDI分子，共含1 145個(gè)原子，HTPB和TDI的質(zhì)量分?jǐn)?shù)分別為86.8%和13.2%。圖3給出了HTPB/TDI混合模型。

圖3 HTPB和TDI的混合模型

1.3 (PEG/NG/BTTN)/NPBA/HMX/AP/PEG/N-100//HTPB/TDI模型搭建與模擬

將得到的(PEG/NG/BTTN)/NPBA/HMX/AP平衡結(jié)構(gòu)與PEG/N-100固化平衡結(jié)構(gòu)，構(gòu)成一個(gè)整體模型，共含10 445個(gè)原子。對(duì)此模型做周期性MD模擬，直到達(dá)到熱力學(xué)平衡。再將獲得的此含固化反應(yīng)的推進(jìn)劑模型結(jié)構(gòu)，與HTPB/TDI襯層固化平衡結(jié)構(gòu)構(gòu)成一個(gè)整體模型，共含11 690個(gè)原子，類(lèi)似地完成MD模擬，最終獲得如圖4所示推進(jìn)劑/襯層平衡結(jié)構(gòu)。

圖4 (PEG/NG/BTTN)/NPBA/HMX/AP/PEG/ N-100//HTPB/TDI體系的平衡結(jié)構(gòu)

對(duì)于以RDX等量取代HMX，模型搭建與模擬方法同HMX一樣，為節(jié)省篇幅，不再贅述。圖5給出了(PEG/NG/BTTN)/NPBA/RDX/AP/PEG/N-100//HTPB/TDI混合體系的平衡結(jié)構(gòu)。

圖5 (PEG/NG/BTTN)/NPBA/RDX/AP/PEG/N-100// HTPB/TDI體系的平衡結(jié)構(gòu)

2 結(jié)果與討論

2.1 組分分子的濃度分布和遷移狀況

從圖4多組分體系的平衡結(jié)構(gòu)可見(jiàn)，在推進(jìn)劑/襯層界面處，亦即在(PEG/NG/BTTN)/NPBA/HMX/AP/PEG/N-100與HTPB/TDI交界層，明顯可見(jiàn)HMX分子的富集，并大致可推斷HMX分子有向界面層遷移的趨勢(shì)。為仔細(xì)探究組分分子的運(yùn)動(dòng)遷移狀況，作為例子以下給出HMX、NPBA和AP分子在OA(X)、OB(Y)和OC(Z)方向上的濃度分布曲線(xiàn)，如圖6所示。

因?yàn)橥七M(jìn)劑/襯層之界面層在OC(Z)方向上，故可細(xì)致考察HMX、NPBA和AP分子在OC(Z)方向上的濃度分布圖。從圖6(c)可見(jiàn)，HMX和NPBA分子在6 ?以?xún)?nèi)分布較少，隨著向界面層接近，HMX和NPBA分子的濃度迅速增加，在30 ?處，達(dá)到頂峰，而后又緩慢下降，直至減小為零。從HMX和NPBA分子的整體濃度分布圖，可清晰看到HMX和NPBA分子有向界面層遷移的趨勢(shì)。這與實(shí)驗(yàn)研究HMX有向界面遷移的結(jié)果一致[14]。而AP分子則基本上呈較平均分布的態(tài)勢(shì)，并無(wú)明顯遷移情況的方向性。

2.2 結(jié)合能和相容性比較

先前曾建議以組分之間的相互作用能(結(jié)合能)表征其相容性的優(yōu)劣[6]。對(duì)于(PEG/NG/BTTN)/NPBA/HMX/AP/PEG/N-100//HTPB/TDI復(fù)雜體系(記為配方I)，可將(PEG/NG/BTTN)/NPBA/AP/PEG/N-100//HTPB/TDI視作一個(gè)整體，記作子體系1，其能量用E1表示；將HMX記為子體系2，其能量用E2表示。對(duì)于含RDX混合體系(記為配方II)也作類(lèi)似處理。

(a)OA(X)方向

(a)OB(Y)方向

(c)OC(Z)方向圖6 HMX、NPBA和AP分子在OA(X)、OB(Y)、OC(Z)方向的濃度分布

定義結(jié)合能Ebind是相互作用能Einter的負(fù)值，即Ebind=-Einter。則配方I中HMX與子體系之間的相互作用能為

Einter=ETotal-E(PEG/NG/BTTN)/NPBA/AP/PEG/N-100//HTPB/TDI-EHMX

這里ETotal(ET)、E(PEG/NG/BTTN)/NPBA/AP/PEG/N-100//HTPB/TDI(E1)和EHMX(E2)分別為混合體系、子體系(PEG/NG/BTTN)/NPBA/AP/PEG/N-100//HTPB/TDI和HMX的平均單點(diǎn)能。配方Ⅱ與此類(lèi)似。計(jì)算結(jié)果見(jiàn)表1。

表1 2種不同配方體系的結(jié)合能

結(jié)合能愈大。表明形成的多組分混合體系越穩(wěn)定，反映了各組分間的彼此相容性越好。所以，曾建議以結(jié)合能Ebind的大小來(lái)預(yù)示和比較高能體系的相容性[6-7，15]。從表1可見(jiàn)，體系Ⅰ的ET(306 573.3 kJ/mol)大于體系Ⅱ的ET(289 438.6 kJ/mol)，即含HMX的體系Ⅰ整體的熱力學(xué)穩(wěn)定性較含RDX的混合體系Ⅱ穩(wěn)定，且HMX與子體系1的Ebind(223 319.2 kJ/mol)明顯大于RDX與其他子體系1的Ebind(202 257.5 kJ/mol)，即含HMX配方的子體系之間的相容性?xún)?yōu)于含RDX配方的子體系之間的相容性。

2.3 力學(xué)性能比較

表2給出了基于MD模擬原子運(yùn)動(dòng)軌跡由靜態(tài)力學(xué)分析方法[16]，求得的含HMX配方體系Ⅰ和含RDX配方體系Ⅱ的力學(xué)性能。

表2 2種不同配方的力學(xué)性能

從表2可見(jiàn)，氧化劑為HMX的配方Ⅰ的彈性系數(shù)、工程模量及柯西壓(C12-C44)和K/G值等力學(xué)性能參數(shù)均高于含RDX的配方Ⅱ，表明前者的剛性和強(qiáng)度均較后者為強(qiáng)，且延展性也以前者為優(yōu)。這就是說(shuō)，若以RDX取代HMX，則推進(jìn)劑/襯層體系(界面)的剛性、強(qiáng)度和延展性將減小，力學(xué)性能將變差。

2.4 安全性能比較

如何利用MD模擬結(jié)果關(guān)聯(lián)感度特性(安全性)，建議重點(diǎn)關(guān)注其中易爆燃組分，因?yàn)橐妆冀M分的引發(fā)鍵是分子中最弱的化學(xué)鍵，在外界刺激下，將優(yōu)先斷裂、引發(fā)熱解或起爆。所以，研究易爆燃組分的引發(fā)鍵很有必要。

表3給出了混合體系Ⅰ和混合體系Ⅱ中引發(fā)鍵的平均鍵長(zhǎng)(Lave)和最大鍵長(zhǎng)(Lmax)。為便于比較，表3中還列出純NG、BTTN、HMX和RDX晶體的MD模擬所得Lave和Lmax。從表3可見(jiàn)，在混合體系中，由于存在多組分之間的相互作用，它們的Lave和Lmax與純晶體中的相應(yīng)值比較，均有所增大。例如，在含HMX混合體系中，其引發(fā)鍵(N—NO2)最大鍵長(zhǎng)(Lmax)從純晶體的1.486 ?增加到1.503 ?；而在含RDX的混合體系中，其引發(fā)鍵(N—NO2)的最大鍵長(zhǎng)(Lmax)則由純晶體的1.512 ?增加到1.528 ?。重要的是比較2種配方，NG和BTTN的O—NO2鍵長(zhǎng)幾無(wú)變化；但HMX引發(fā)鍵(N—N)的最大鍵長(zhǎng)Lmax(1.503 ?)遠(yuǎn)小于RDX的相應(yīng)值(1.528 ?)。這就是說(shuō)，以RDX取代HMX后，新配方中引發(fā)鍵(N—NO2)的Lmax遠(yuǎn)大于原配方的相應(yīng)值。以建議的引發(fā)鍵Lmax判據(jù)[17-20]，可預(yù)測(cè)新配方的熱和撞擊等感度將大于原配方，這主要?dú)w因于RDX比HMX的感度大[21-22]。

表3 2個(gè)配方中引發(fā)鍵的平均鍵長(zhǎng)(Lave)和最大鍵長(zhǎng)(Lmax)

3 結(jié)論

(1)搭建了(PEG/NG/BTTN)/NPBA/HMX/AP/PEG/N-100//HTPB/TDI復(fù)雜的推進(jìn)劑/襯層模型(配方I)，并以RDX等量取代HMX，搭建相應(yīng)的配方II模型。首次對(duì)2種復(fù)雜的配方體系進(jìn)行了系統(tǒng)的MD比較研究。

(2)展示了配方I中組分分子的濃度分布和遷移狀況，發(fā)現(xiàn)HMX和NPBA組分的分子有向界面層遷移的趨勢(shì)，而AP分子則呈平均分布態(tài)勢(shì)。

(3)配方II子體系之間的結(jié)合能比配方I小，表明相容性不如配方I；配方II工程模量減小，表明剛性和強(qiáng)度下降；柯西壓(C12-C44)和K/G值亦減小，表明延展性亦減??；配方II中引發(fā)鍵N—NO2的Lmax增大，表明感度將增大，安全性下降?？梢?jiàn)，含RDX的配方II沒(méi)有含HMX的配方I好。本文首次提供了MD對(duì)復(fù)雜高能體系配方設(shè)計(jì)的指導(dǎo)作用。

[1] 姚維尚,李倩,譚惠民.NEPE推進(jìn)劑粘合劑性能的分子模擬研究[J].含能材料,2007,15(6):650-655.

[2] 楊月誠(chéng),焦東明,強(qiáng)洪夫,等.HTPB推進(jìn)劑組分溶度參數(shù)的分子模擬研究[J].含能材料,2008,16(2):191-195.

[3] 付一政,劉亞青,蘭艷花.端羥基聚丁二烯/增塑劑共混物相容性的分子動(dòng)力學(xué)模擬[J].物理化學(xué)學(xué)報(bào),2009,25(7):1267-1272.

[4] 趙貴哲,馮益柏,付一政,等.端羥基聚丁二烯/增塑劑共混物相容性的分子動(dòng)力學(xué)模擬和介觀模擬[J].化學(xué)學(xué)報(bào),2009,67(19):2233-2238.

[5] 焦東明,楊月誠(chéng),強(qiáng)洪夫,等.HTPB固體推進(jìn)劑增塑劑選取分子模擬研究[J].化學(xué)研究與應(yīng)用,2009,21(6):805-809.

[6] Xu X J,Xiao J J,Huang H,et al.Molecular dynamics simulations on the structures and properties ofε-CL-20-based PBXs[J].Sci.China Ser.B,2007,50(6):737-745.

[7] 于艷春,朱偉,肖繼軍,等.四組分高能體系結(jié)合能和力學(xué)性能的分子動(dòng)力學(xué)模擬[J].化學(xué)學(xué)報(bào),2010,68(12):1181-1187.

[8] 付一政,胡雙啟,蘭艷花,等.HTPB/增塑劑玻璃化轉(zhuǎn)變溫度及力學(xué)性能的分子動(dòng)力學(xué)模擬[J].化學(xué)學(xué)報(bào),2010,68(8):809-813.

[9] 李紅霞,強(qiáng)洪夫,武文明.丁羥推進(jìn)劑黏結(jié)體系中增塑劑遷移的分子模擬[J].火炸藥學(xué)報(bào),2008,31(5):74-78.

[10] 李紅霞,強(qiáng)洪夫,王廣,等.基于MD方法的增塑劑擴(kuò)散行為的模擬研究[J].含能材料,2009,17(1):36-41.

[11] 史佩,李高春,王玉峰,等.復(fù)合推進(jìn)劑顆粒填充模型的分子動(dòng)力學(xué)模擬方法[J].計(jì)算機(jī)與應(yīng)用化學(xué),2007,24(5):665-668.

[12] Zhu W,Wang X J,Xiao J J,et al.Molecular dynamics simulations of AP/HMX composite with a modified force field[J].Journal of Hazardous Materials,2009,167(1):810-816.

[13] 李紅霞,強(qiáng)洪夫,武文明.HTPB與TDI固化的分子模擬研究[J].固體火箭技術(shù),2008,31(6):602-606.

[14] 吳豐軍,彭松,池旭輝,等.NEPE推進(jìn)劑/襯層粘接界面XPS表征[J].固體火箭技術(shù),2009,32(2):192-196.

[15] Xiao J J,Huang H,Li J S,et al.A molecular dynamics study of interface interactions and mechanical properties of HMX-based PBXs with PEG and HTPB[J].Journal of Molecular Structure:THEOCHEM,2008,851:242-248.

[16] Swenson R J.Comments on virial theorems for bounded systems[J].Am.J.Phys.,1983,51:940-942.

[17] 朱偉,肖繼軍,鄭劍,等.高能混合物的感度理論判別——不同配比和不同溫度AP/HMX的MD研究[J].化學(xué)學(xué)報(bào),2008,66(23):2592-2596.

[18] Xiao J J,Zhao L,Zhu W,et al.Molecular dynamics study on the relationships of modeling,structural and energy properties with sensitivity for RDX-based PBXs[J].Sci.China Ser.B,2012,55:2587-2594.

[19] Xiao J J,Wang W R,Chen J,et al.Study on structure,sensitivity and mechanical properties of HMX and HMX-based PBXs with molecular dynamics simulation[J].Comput.Theor.Chem.,2012,999:21-27.

[20] Xiao J J,Li S Y,Chen J,et al.Molecular dynamics study on the correlation between structure and sensitivity for defective RDX crystals and their PBXs[J].J.Mol.Model,2013,19:803-809.

[21] 董海山,周芬芬.高能炸藥及其相關(guān)性能[M].北京:科學(xué)出版社,1989:252-253,266-267.

[22] 肖繼軍,朱衛(wèi)華,朱偉,等.高能材料分子動(dòng)力學(xué)[M].北京:科學(xué)出版社,2013:120-129.

(編輯：劉紅利)

Molecular dynamics simulation of the structures and properties of NEPE propellant/liner(Ⅱ) —Demonstrating component molecule migrating and formulation for complex system

ZHU Wei1，2，LIU Dong-mei1，XIAO Ji-jun1，CHI Xu-hui3，PANG Ai-min3，XIAO He-ming1

(1.Molecules and Materials Computation Institute,Nanjing University of Science and Technology,Nanjing 210094,China； 2.College of Biological,Chemical Sciences and Engineering,Jiaxing University,Jiaxing 314001,China； 3.The 42nd Institute of the Fourth Academy of CASC,Xiangyang 441003,China)

Molecular dynamics(MD)simulation was applied to propellant/liner model(PEG/NG/BTTN)/NPBA/HMX/AP/PEG/N-100//HTPB/TDI at 295 K with NVT ensemble.Component concentration distribution and migration behavior were presented and analyzed.The simulation results show that HMX and NPBA have a tendency to migrate to interface,while AP tends to have a homogeneous distribution.Moreover,a new model built by replacing HMX with RDX equivalently was studied by the molecular dynamic simulation.It is found that tensile modulus(E),bulk modulus(K),shear modulus(G)and Cauchy pressure(C12-C44)of the new formula decreas apparently.This implies that its stiffness,strength and ductibility are all decreased;The maximum bond length of the N—NO2trigger bonds of the new formula(1.528 ?) is larger than that of the original formula(1.503 ?),indicating that the sensitivity increases and safety decreases of the new formula.The results shows that the binding energies of HMX and other constituent is larger than RDX and other constituent indicate that compatibility of HMX and other constituent,which is better than RDX and other constituent.

propellant/liner interface;constituent migration;mechanical property;sensitivity;compatibility;molecular dynamics simulation

2013-07-31；

2014-01-16。

國(guó)防973項(xiàng)目(613142010102；613142010302)。

朱偉(1977—)，男，博士/副教授，研究方向?yàn)橛?jì)算化學(xué)及其應(yīng)用。E-mail：zhuwmail@126.com

肖繼軍(1964—)，男，博士生導(dǎo)師，主要從事高分子物理和分子材料學(xué)研究。E-mail：xiao_jijun@163.com

V512

A

1006-2793(2014)05-0678-06

10.7673/j.issn.1006-2793.2014.05.018