陳思彤，董可海，王 鑫，裴立冠，孔令澤，夏 成

（1. 海軍工程大學兵器工程學院，湖北 武漢 430033；2. 海軍航空大學岸防兵學院，山東 煙臺 264001）

1 引言

硝酸酯增塑聚醚（NEPE）固體推進劑的危險品等級為1.1 級，當受到?jīng)_擊波和撞擊刺激時可能發(fā)生爆轟［1］，其中增塑劑硝化甘油（NG）的敏感性較高、并且存在一定程度的遷移［2］；增塑劑的大量遷移會引起“推進劑?襯層”界面的脫粘、影響導彈的彈道性能、增加發(fā)動機點火后爆炸的可能性等［3］。

目前，關(guān)于組分遷移的研究方法主要有實驗法［4］和分子模擬法［5］。實驗法主要是采用聯(lián)合試件、利用色譜分析或者X 射線光電子能譜（XPS）分析等手段［6-7］，雖然能夠得到最接近于真實條件的遷移量，但存在操作復雜、XPS 對于硝酸酯的檢測不準確等不足［8］；近年來分子模擬法快速發(fā)展［9］，成為預測擴散性能、研究擴散微觀機理的有效手段，具有耗時短、花費低等優(yōu)點［10-11］。

對于新型含能鈍感增塑劑三羥甲基乙烷三硝酸酯（TMETN）［12-13］，法國、日本對其應(yīng)用于雙基推進劑［14］、聚疊氮縮水甘油醚（GAP）推進劑中的性能進行了研究［15］，而將TMETN 應(yīng)用于NEPE 推進劑的研究還比較少。為此，本研究采用分子模擬法比較三種增塑劑TMETN、NG、1，2，4?丁三醇三硝酸酯（BTTN）在NEPE推進劑黏合劑預聚體PEG/N?100 中的擴散性能及機理，并討論溫度和增塑比對TMETN 擴散性能的影響，為鈍感高能推進劑的發(fā)展提供理論依據(jù)。

2 模型構(gòu)建與模擬方法

2.1 分子動力學建模及模擬細節(jié)

真實推進劑的成分和結(jié)構(gòu)都較為復雜，鑒于分子模擬（MD）時間和空間尺度的限制，考慮到固化過程是N?100 中異氰酸酯基與PEG 中羥基發(fā)生反應(yīng)，所以構(gòu)建了高分子預聚體網(wǎng)絡(luò)的一部分（PEG 分子量為6002［16］），如圖1 所示。

圖1 PEG/N?100 聚合物鏈片段Fig.1 Fragments of PEG/N?100 polymer chains

圖2 為優(yōu)化前和優(yōu)化后的高分子鏈構(gòu)型，可見優(yōu)化過程使高分子鏈從直鏈結(jié)構(gòu)變得纏繞卷曲，同時能量也從8188.83 kJ·mol-1降低到-9202.22 kJ·mol-1，因為真實物質(zhì)的能量總是最低的，所以優(yōu)化作用明顯使得分子構(gòu)型更接近于實際［17］。

圖2 優(yōu)化前和優(yōu)化后的預聚物構(gòu)型Fig.2 The configuration of prepolymer before and after opti?mization

PEG 與增塑劑的質(zhì)量比（增塑比）約為2.5，為減小鏈端效應(yīng)［18］，將1 條PEG/N?100 高分子鏈與相應(yīng)數(shù)量的硝酸酯構(gòu)成周期箱，NG、BTTN、TMETN 的分子數(shù)量分別為：65、61、58 個，初始密度按照三種物質(zhì)的體積平均值進行設(shè)定；對周期箱進行smart 方法的幾何優(yōu)化，再進行10 次NPT 系綜的退火處理，退火的初始溫度為298 K、中間溫度為598 K、最終溫度為298 K［17］，每次退火后再進行幾何優(yōu)化，并選取能量較低的構(gòu)型進行1 ns NPT 系綜的密度優(yōu)化，時間步長均為1 fs，然后進行625 ps NVT 系綜動力學計算、其中后400 ps用于采集均方位移隨時間變化曲線，模擬控溫器和控壓器均采用Bederson，有文獻證實Bederson 控溫器對均方位移產(chǎn)生的影響較小［19］。

2.2 介觀建模及模擬細節(jié)

在進行介觀動力學模擬之前，首先要確定各分子所包含的高斯鏈段數(shù)量、各高斯鏈段間的相互作用參數(shù)。高分子預聚物的高斯鏈段數(shù)量有如下公式［20］：

式中，Nmeso為高斯鏈段數(shù)量；n為預聚物的聚合度，此處為136（PEG 分子量為6002）；C∞為全原子高分子鏈的特征比。

由Synthia 模塊得到PEG 的C∞是4.98，故PEG 高斯鏈段的質(zhì)量為219.12（特征比乘以單體的相對分子質(zhì)量），所以1 條PEG 分子鏈等于27 個高斯鏈段、1 個增塑劑分子等于1 個高斯鏈段。

在MesoDyn 中，兩種高斯鏈段的相互作用ν?1εij表示為：

式中，ν-1εij單位為kJ·mol-1；R為氣體常數(shù)，8.314 J·mol-1·K-1；T為溫度，K；χij為Flory?Huggins 參數(shù)：式中，Δδ是兩物質(zhì)的溶度參數(shù)差值，（J·cm-3）1/2；Vr取兩種高斯鏈段摩爾體積的平均值，cm3。

根據(jù)文獻［21］中MD 計算的各物質(zhì)溶度參數(shù)，計算得到兩種高斯鏈段間相互作用參數(shù)ν-1εij、并列于表1。在Meso Dyn 模塊中根據(jù)增塑比設(shè)置珠子的混合比例，三維周期箱尺寸為32.0 nm×32.0 nm×32.0 nm，溫度為298 K，步長為20 ns，總時間為1000 μs。

表1 增塑劑、PEG 的溶度參數(shù)及高斯鏈段間的相互作用參數(shù)Table 1 Solubility parameters of plasticizers and PEG，as well as the interaction parameters among the Gauss chains

3 結(jié)果與討論

3.1 三種增塑劑的擴散性能

首先，應(yīng)闡述“擴散”與“遷移”存在的聯(lián)系與區(qū)別，擴散是發(fā)生在推進劑內(nèi)部，而遷移是指推進劑內(nèi)部的組分擴散至界面處、并進一步擴散到其它層次中（比如襯層），擴散是遷移的前提、可以在一定程度上預測遷移性能，但遷移較擴散更為復雜、還涉及到遷移組分與襯層材料的相互作用等，本節(jié)僅比較TMETN、NG、BTTN 三種增塑劑在預聚體內(nèi)部的擴散性能、不考慮襯層種類對增塑劑擴散的影響。

采用MD 計算擴散系數(shù)的常用方法是Einstein 關(guān)系式［22］：通過MD 計算出均方位移?時間斜率m，即可得到擴散系數(shù)D。

采集后400 ps 增塑劑的均方位移，統(tǒng)計長度采用200 ps，原點間隔采用1 幀（0.01 ps），得到MSD隨時間的變化曲線，如圖3 所示。

圖3 三種增塑劑的均方位移隨時間變化曲線Fig.3 Mean?square displacement（MSD）?t curves of plasti?cizers

由表2 可見，NG、BTTN 在PEG/N?100 中擴散系數(shù)的數(shù)量級為10-11、TMETN 擴散系數(shù)的數(shù)量級為10-12，將本次模擬的結(jié)果與類似的文獻模擬結(jié)果進行比較：李紅霞［23］采用MD 計算了DOS在HTPB 中的擴散系數(shù)為2×10-11，王曉［14］采用MD 計算NG 在PET/MDI、TDI、IPDI 中的擴散系數(shù)分別為6.5×10-13、1.25×10-12、1.11×10-11；本文研究的預聚體、固化劑與文獻中的稍所不同，但模擬結(jié)果的數(shù)量級在同一范圍內(nèi)，說明本模擬與文獻結(jié)果的吻合性較好。由表2 可見，三種增塑劑在預聚體PEG/N?100 中擴散系數(shù)的大小順序為：NG＞BTTN＞TMETN，結(jié)果表明NG 最容易在黏合劑基體中發(fā)生擴散，BTTN 次之，TMETN 最難發(fā)生擴散。

表2 三種增塑劑均方位移?時間曲線的斜率及其對應(yīng)的擴散系數(shù)Table 2 Slope of MSD?t curves and the diffusion coefficient of three plasticizers

3.2 三種增塑劑的擴散機理

據(jù)文獻［14，24-25］，影響增塑劑擴散的微觀因素主要有：增塑劑與預聚體的分子間相互作用、體系的自由體積分數(shù)和增塑劑分子的尺寸，下面從這三個角度分析增塑劑在PEG/N?100 預聚體中的微觀擴散機理。

3.2.1 增塑劑與預聚體的分子間相互作用

徑向分布函數(shù)g（r）可以反映體系的相互作用類型（峰值位置）及大?。ㄇ€高低）［26］，在混合物的徑向分布函數(shù)曲線中，異類分子間的相互作用曲線高于同類分子間的相互作用曲線，則體系相容性較好［27，29］。

由圖4 可見：TMETN/PEG/N?100 體系的相容性較好，而在另外兩種體系中，增塑劑與增塑劑分子間的相互作用大于增塑劑與PEG/N?100 分子間的相互作用（NG 所在的體系更為明顯），所以，BTTN 與PEG/N?100 的相容性較差、NG 與PEG/N?100 的相容性最差；體系相容性的順序與擴散系數(shù)的順序在邏輯上保持一致，即體系的相容性越好，增塑劑越難發(fā)生擴散。

圖4 NG/PEG、BTTN/PEG、TMETN/PEG 混合物的分子間徑向分布函數(shù)Fig.4 Intermolecular radial distribution functions for NG/PEG，BTTN/PEG and TMETN/PEG

相容性就是復雜的分子間相互作用，為了進一步比較分析各體系內(nèi)分子間相互作用的差異，將混合物中分子間徑向分布函數(shù)按照增塑劑?PEG、PEG?PEG、增塑劑?增塑劑的作用類型進行處理，結(jié)果如圖5 所示。

在三個體系中，PEG 的分子數(shù)量相同，NG、BTTN、TMETN 的分子數(shù)量依次遞減，由圖5a 可見，雖然體系中TMETN 的分子數(shù)量最少，但該體系中PEG分子附近出現(xiàn)TMETN 的概率卻最大，故PEG 與TMETN 的結(jié)合能力最強，同理可發(fā)現(xiàn)，BTTN 與PEG的結(jié)合能力居中、NG 與PEG 的結(jié)合能力最弱；由圖5b 可見，NG/PEG 體系中PEG 的自聚集能力最強，BTTN/PEG 體系中次之、TMETN/PEG 體系中最弱；由圖5c 可見，因三個體系中增塑劑分子數(shù)量的大小順序與該徑向分布函數(shù)曲線的高低一致，故無法判定“NG附近出現(xiàn)NG 概率最大”是由NG 的自聚集能力強引起的、還是由NG 分子數(shù)量多引起的。綜上，PEG 與TMETN 的結(jié)合能力最強、該體系中PEG 的自聚集能力最弱，這也是TMETN 不易發(fā)生擴散的微觀原因。

圖5 三個體系中增塑劑?PEG、PEG?PEG、增塑劑?增塑劑的分子間徑向分布函數(shù)Fig.5 Intermolecular radial distribution functions of plasticizer?PEG，PEG?PEG and plasticizer?plasticizer in three systems

由于分子間徑向分布函數(shù)曲線是對各原子間相互作用進行平均而得到的，所以可能抹平原子間的氫鍵作用，氫鍵作用的范圍為2.6～3.1 nm，3.1 nm 以后為范德華作用的范圍，針對于此，對可能存在氫鍵的增塑劑與PEG/N?100 原子對進行徑向分布函數(shù)分析，最終發(fā)現(xiàn)PEG 中的H 原子與增塑劑中的硝基O 原子存在氫鍵作用。

三個體系中氫鍵作用的結(jié)果見圖6 和表3。由圖6和表3可見，三個體系均存在氫鍵作用，且TMETN_O 與PEG/N100_H 的氫鍵作用最強，峰值比BTTN 體系中的高0.096、比NG 體系中的高0.1456，說明原子間氫鍵作用越強、增塑劑越難發(fā)生擴散。

圖6 增塑劑與PEG 間O—H 原子的徑向分布函數(shù)Fig.6 Radial distribution functions between the oxygen atom of plasticizer and hydrogen atom of PEG

表3 氫鍵作用的位置及峰值Table 3 Location and peak of hydrogen bond interaction

增塑劑的擴散會導致體系的混亂程度增加、引起相分離，這將對裝藥的性能產(chǎn)生較大威脅，而相分離的時間和特征區(qū)尺寸都超出了MD 的范圍、需借助介觀動力學模擬（Meso Dyn）來研究。根據(jù)體系的增塑比為2.5，參數(shù)設(shè)置如表1，最終得到三個體系的密度分布如圖7 所示，顯示顏色的地方代表著密度高于該物質(zhì)平均密度的位置和水平。

從圖7可見，TMETN/PEG體系中的兩相分布較均勻，而BTTN/PEG 和NG/PEG 體系發(fā)生了相聚集、且后者更為嚴重，這也說明TMETN在PEG相中發(fā)生擴散的能力最弱，而BTTN和NG發(fā)生擴散近而聚集的能力更強；此外還發(fā)現(xiàn)，PEG易向周期箱四周聚集、增塑劑則是向中間聚集。

圖7 三個體系中的濃度分布Fig.7 Concentration distribution in three systems

3.2.2 體系的自由體積分數(shù)

體系內(nèi)的自由體積為高分子鏈的運動和小分子的擴散提供空間，MD 中自由體積表示為：

其中，Vf是自由體積，Vsp是周期箱體積，Vw是范德華表面占有體積，單位均為?3。

自由體積是指除去PEG/N?100 分子和增塑劑分子共同占有體積后的剩余體積，采用Atom Volumes& Surfaces 工具對三種增塑體系的自由體積進行計算，探針半徑為1.4 ?，因為各體系的總體積大小不一樣，為了能夠比較三種增塑劑的相對運動空間，定義自由體積分數(shù)FFVSim.為自由體積與各周期箱體積的比值：

體系的自由體積分數(shù)越大，提供給增塑劑分子運動的相對空間就越大，擴散也就越容易發(fā)生。由表4可見，體系內(nèi)自由體積分數(shù)的大小順序為TMETN/PEG/N?100＞NG/PEG/N?100＞BTTN/PEG/N?100，可見TMETN 具有最大的相對運動空間，但在“自由體積分數(shù)”和“TMETN 與PEG/N?100 的強分子間作用”的綜合影響下，TMETN 的擴散系數(shù)還是最小的；此外，NG/PEG/N?100 體系中的自由體積分數(shù)居中、BTTN/PEG/N?100 體系的最小，這說明NG 的相對運動空間較BTTN 的 大，又 由3.2.1 節(jié) 知“NG 與PEG/N?100 的相容性也較差”，這兩種因素的綜合作用導致了NG 的擴散系數(shù)最大，BTTN 次之。

表4 體系的自由體積分數(shù)Table 4 Fractional free volume of systems

3.2.3 增塑劑分子的尺寸

增塑劑分子的大小也會影響擴散性能，分子越小，越容易在聚合物孔洞中跳躍；由增塑劑分子結(jié)構(gòu)可知，三種增塑劑分子的大小順序為：TMETN＞BTTN＞NG，這與擴散系數(shù)的順序在邏輯上是一致的，即NG 分子最小、最容易在聚合物大分子孔洞中發(fā)生跳躍，BTTN次之，TMETN 最難。

綜上，從“分子間相互作用”、“體系的自由體積分數(shù)”、“增塑劑分子的尺寸”三個方面討論了增塑劑擴散性能的微觀機理，在TMETN/PEG/N?100 體系中，體系相容性最好、TMETN 分子尺寸最大，這兩方面都加大了擴散的難度，但同時該體系的自由體積分數(shù)也最大，這方面會促進擴散的發(fā)生，最終在三種因素的綜合作用下，TMETN 的擴散系數(shù)還是最小的；在NG/PEG/N?100 體系中，體系的相容性最差、NG 的分子尺寸最小、體系的自由體積分數(shù)較BTTN/PEG/N?100 的大，這三種因素均導致NG 的擴散系數(shù)比BTTN 的大。

3.3 溫度對TMETN 擴散性能的影響

3.3.1 溫度對體系相容性的影響

根據(jù)推進劑儲存、使用及固化等條件，比如北方冬天低溫可達0 ℃，常溫是25 ℃，固化溫度大約為50 ℃，所以選取0、25、40、50 ℃，分析TMETN/PEG/N?100 體系中TMETN 的擴散系數(shù)，采集的MSD 隨時間變化曲線如圖8，并以圖8 為根據(jù)、經(jīng)過擬合得到擴散系數(shù)，將擴散系數(shù)列于表5。不同溫度下擴散系數(shù)的變化如圖9 所示。

圖8 不同溫度下TMETN 均方位移隨時間的變化曲線Fig.8 Mean?square displacement（MSD）?t curves of plasti?cizers at various temperatures

由表5 和圖9 可見，隨著溫度升高，擴散系數(shù)先緩慢增加、后劇烈增加，即溫度越高、擴散系數(shù)增加的速率越大，這與文獻［23］中的規(guī)律一致、也與高溫加速老化的規(guī)律一致；常溫298 K 時擴散系數(shù)為低溫273 K的1.2 倍左右；固化溫度323 K 時的擴散系數(shù)大約為常溫298 K 時擴散系數(shù)的2 倍，這解釋了為何增塑劑在固化過程中已發(fā)生較大的遷移［30］。

表5 不同溫度下增塑劑TMETN 的均方位移?時間曲線的斜率及對應(yīng)的擴散系數(shù)Table 5 The slopes of MSD?t curves and the diffusion coeffi?cient of TMETN at various temperatures

圖9 不同溫度下擴散系數(shù)的變化Fig.9 The diffusion coefficient at various temperatures

進一步分析不同溫度下TMETN 中硝基O 原子與PEG 中H 原子間的氫鍵作用，為了便于比較各溫度下氫鍵的峰值和位置，將273，298，313，323 K 條件下的徑向分布函數(shù)曲線均置于圖10 中。不同溫度下氫鍵作用的位置及峰值見表6。

圖10 不同溫度下TMETN 與PEG 間H—O 原子的徑向分布函數(shù)曲線Fig.10 Curves for radial distribution functions between plas?ticizer（O）and PEG（H）at various temperatures

由 圖10 和 表6 可 見，隨 著 溫 度 升 高，TMETN 與PEG 原子間氫鍵作用峰的位置后移、峰值變低，這說明高溫使增塑劑與黏合劑間的相互作用變?nèi)酰@導致了TMETN 更容易發(fā)生擴散。

表6 不同溫度下氫鍵作用的位置及峰值Table 6 Location and peak of hydrogen bond interaction at various temperatures

3.3.2 溫度對體系自由體積分數(shù)的影響

由表7 可見，隨著溫度升高，TMETN/PEG/N?100體系的自由體積分數(shù)變大，說明高溫使聚合物的活動性增強、分子間距離增加，這給增塑劑提供了更大的活動空間，在宏觀上表現(xiàn)為聚合物的溶脹度提高［23］；并且溫度的升高也會增加小分子的動能，即小分子的跳躍能力變強，這均使增塑劑更易發(fā)生擴散。

表7 不同溫度下體系的自由體積分數(shù)Table 7 Fractional free volume for systems at various tem?peratures

3.4 增塑比對TMETN 擴散性能的影響

根據(jù)NEPE 推進劑中增塑劑的含量，研究了在常溫298 K 條件下，增塑比為2.5、2.8、3 時TMETN 的擴散 性 能［31］，得 到MSD?t曲 線 擬 合 斜 率 及 擴 散 系數(shù)（表8）。

表8 增塑劑TMETN 的均方位移?時間曲線斜率及對應(yīng)的擴散系數(shù)Table 8 Slopes of MSD?t curves and the diffusion coeffi?cients of TMETN

由表8 可見，隨著增塑比的增加，TMETN 擴散系數(shù)在減小，為了觀察增塑比對體系形態(tài)分布的影響，采用MesoDyn 進行了1000 μs 的介觀動力學模擬，得到TMETN 等密度圖及有序度參數(shù)。由圖11 可見，隨著增塑比的增加，組分發(fā)生同相歸并的程度變小。由圖12 可見，TMETN 的有序度參數(shù)減小、體系的熵增減小，所以體系相容性變好，這是TMETN 的擴散系數(shù)減小的一個原因；此外有研究表明［23］，隨著增塑劑含量的增加，大量的增塑劑會導致高分子網(wǎng)絡(luò)收縮、使得增塑劑在孔隙中發(fā)生有效跳躍的難度變大、從而使擴散系數(shù)減小，這一點在達到其共容增塑極限后更為明顯。

圖11 不同增塑比條件下TMETN 的等密度圖Fig.11 The isosurface of the density fields for TMETN at vari?ous plasticization ratios

圖12 不同增塑比條件下TMETN 的有序度參數(shù)Fig.12 The order parameters for TMETN at various plastici?zation ratios

4 結(jié)論

（1）增塑劑在PEG/N?100 中的擴散系數(shù)大小順序為：NG＞BTTN＞TMETN，這說明含能鈍感增塑劑TMETN 具有較強地抗擴散能力，可以進一步研究將其代替或部分代替NG、應(yīng)用于NEPE 推進劑中。

（2）從增塑劑與黏合劑的分子間相互作用、體系的自由體積分數(shù)和增塑劑分子的尺寸三個方面分析了增塑劑擴散性能的微觀機理，在三個增塑體系中，TMETN 與PEG/N?100 的分子間相互作用最強、體系的自由體積分數(shù)最大、TMETN 分子的尺寸最大，各因素綜合作用導致了TMETN 的擴散系數(shù)最小。

（3）隨著溫度升高，TMETN 擴散系數(shù)先緩慢增加、后劇烈增加，這與高溫加速老化的規(guī)律一致；固化溫度323 K 時的擴散系數(shù)為常溫298 K 時的2 倍左右，這在一定程度上解釋了“增塑劑在固化過程中已發(fā)生較大遷移”的現(xiàn)象；從微觀角度來解釋擴散系數(shù)增大的原因為：高溫導致TMETN 與PEG 的原子間氫鍵作用減弱、體系的自由體積分數(shù)變大。

（4）隨著增塑比增加，TMETN 擴散系數(shù)減小，介觀動力學研究表明“體系的有序度參數(shù)減小、體系相容性變好”是其中的原因之一。