丁銀鳳，梁 昊，丁亞軍，肖忠良

（1. 南京理工大學化工學院，南京 210094；2. 南京理工大學特種能源材料教育部重點實驗室，南京 210094）

1 引言

表面鈍感是實現發(fā)射藥燃燒高漸增性的常用途徑之一，鈍感劑通過擴散作用滲入發(fā)射藥基體，在發(fā)射藥表層呈由表及里下降的趨勢，以降低燃燒初期的氣體生成速率。當前中小口徑武器所使用的鈍感發(fā)射藥，主要采用鄰苯二甲酸二丁酯（DBP）、樟腦等作為鈍感功能組分，取得了較好的彈道效果。但是DBP、樟腦分子量較低，在長時間貯存過程中易遷移，引起鈍感劑濃度和分布的變化，對發(fā)射藥的性能穩(wěn)定性和使用壽命產生負面影響［1］。為降低鈍感劑在發(fā)射藥貯存過程中的遷移速率，提高發(fā)射藥的使用壽命，目前高分子鈍感劑（聚新戊二醇己二酸酯NA）已成功應用到EI 發(fā)射藥（改性單基發(fā)射藥）中，較高的分子量提高了發(fā)射藥的長儲穩(wěn)定 性［2-4］。

但發(fā)射藥貯存過程中，DBP 和NA 的遷移現象仍然存在，國內外主要采用顯微紅外光譜［5］以及激光顯微共聚焦拉曼光譜［6］等實驗手段測定鈍感劑的濃度梯度分布數據，并依據Fick 第二定律的高斯解進行曲線擬合從而得到擴散系數，但存在耗時長、檢測不準確、操作復雜等不足。目前，也有國外文獻報道Bertrand R 通過等溫實驗確定擴散系數與溫度的關系，并考慮發(fā)射藥基體膨脹效應的影響，建立方程，利用AKTS-SML 軟件預測擴散系數，但目前在非穩(wěn)定實驗條件下（例如老化期間的溫度波動、膨脹效應），模擬計算鈍感劑在發(fā)射藥的遷移速率具有一定的局限性，而且對于其擴散的機理也沒有研究［7］。

近年來，分子動力學模擬（MD 模擬）已成功應用于含能材料領域，尤其在力學性能［8］、相容性［9］等方面已經取得一系列研究成果，同時MD 模擬也成為預測擴散性能、研究擴散微觀機理的有效手段，可以得到擴散運動的直觀過程并計算擴散系數［10］。李紅霞［11］等采用分子動力學方法模擬增塑劑在丁羥粘合劑體系中的運動，通過愛因斯坦關系式求得擴散系數，再與采用加速老化實驗方法計算得到的擴散系數對比，模擬值和實驗值基本一致，數量級為10-12m2·s-1。MD 模擬也能夠從物質的微觀信息角度提供理論支持，還具有耗時短、花費低等優(yōu)點。

本研究針對發(fā)射藥貯存過程中鈍感劑的遷移現象影響發(fā)射藥的使用壽命，采用分子動力學方法比較DBP 和NA 兩種不同分子量的鈍感劑在發(fā)射藥中的擴散性能及機理，并探究溫度、硝化甘油（NG）含量對鈍感劑擴散性能的影響，以期為鈍感發(fā)射藥壽命的預估提供理論指導。

2 模擬計算

2.1 物理建模

依據NC（含氮量為13.1%）、NG、DBP 與NA 的化學結構式，利用Materials Studio 軟件的Visualizer 模塊，分別建立分子模型，相對分子量分別為11273、228、277、3000；然后根據配方組成，利用Amorphous Cell 模塊構建含有11 個DBP 分子、5 個NC 分子鏈組成的DBP/NC 體系，兩種組分的質量比為1∶19；另外構建含有1 個NA 分子、5 個NC 分子鏈組成的NA/NC體系，兩種組分的質量比為1∶19；構建DBP/NG/NC、NA/NG/NC 體系，組分之間的質量比見3.1.2 節(jié)的表1。所構建的DBP/NC、NA/NC、DBP/NG/NC、NA/NG/NC體系模型見圖1。

表1 DBP、NA 在雙基發(fā)射藥體系中的擴散系數Table 1 Diffusion coefficients of DBP and NA in double base gun propellant

圖1 DBP/NC、NA/NC、DBP/NG/NC、NA/NG/NC 體系模型Fig.1 The system models of DBP/NC，NA/NC，DBP/NG/NC and NA/NG/NC

2.2 MD 模擬細節(jié)

采用Amorphous Cell 構建DBP/NC、NA/NC 共混體系時，應用周期性邊界條件，即以立方元胞為中心，周圍有26 個相鄰的鏡像立方元胞，以達到利用較少分子模擬宏觀性質的目的。利用Forcite 模塊對建立的體系模型先進行幾何優(yōu)化（Geometry Optimization），幾何優(yōu)化是為了尋找局部最優(yōu)構型，優(yōu)化過程均采用Smart Minization 方法，選用COMPASS 立場和Fine 精度，分別用Atom-based［12］和Ewald［13］方法求范德華作用和靜電作用；其次利用Anneal 進行若干次退火處理，退火可以使分子跨過局部能壘、尋找全局最優(yōu)構型，退火的初始溫度為300 K、中間溫度為800 K、最終溫度為300 K，每次退火后自動進行幾何優(yōu)化，這是在全局最優(yōu)構型中進一步尋找，從中選取能量最低的一幀構型作為輸出的最優(yōu)構型。模擬過程中位能均采用球形截斷法，非鍵截斷半徑（cutoff distance）取1.5 nm，樣條寬度（spline width）取0.1 nm，緩沖寬度（buffer width）取0.05 nm，時 間 步 長 為1fs，采 用Andersen 控溫方法［14］和Berendsen 控壓方法［15］。

利用Forcite 模塊的Dynamic 進行動力學模擬，先運行50 ps 的正則系綜NVT（系統(tǒng)的粒子數N、體積V和溫度T恒定）進行弛豫，然后選擇恒溫恒壓系綜NPT（系統(tǒng)的粒子數N、壓強p和溫度T恒定），在1 個標準大氣壓下，分別在5，25，65，85 ℃下進行1000 ps 的MD 模擬，保存全運動軌跡，后400 ps 體系已達到溫度和能量平衡（溫度和能量隨時間的變化率小于5%），選擇后400ps 期間的軌跡進行分析均方位移、徑向分布函數、自由體積。用同樣方法構建DBP/10NG/NC、DBP/15NG/NC、NA/10NG/NC、NA/15NG/NC 體系，并采用同樣的參數進行能量優(yōu)化、MD 模擬。

2.2.1 擴散系數

擴散系數是表征分子擴散遷移能力的重要參數，指的是沿擴散方向在單位時間每單位濃度梯度的條件下，垂直通過單位面積所擴散某物質的質量或摩爾數。分析MD 模擬過程中所得到的粒子的運動軌跡，可以得到均方位移曲線（MSD曲線）。均方位移（MSD）反映模擬體系中分子的空間位置對初始位置的偏移程度，MSD與時間（t）的關系式如方程（1）：

式中，rt表示t時刻分子的坐標；r0表示分子的初始坐標。

通過MSD曲線擬合計算出均方位移-時間斜率（m，m2·s-1），將 其 代 入 計 算 分 子 擴 散 系 數D的Einstein-Smoluchowski［16］方程，如方程（2）：

2.2.2 徑向分布函數

徑向分布函數（RDF）是反映材料微觀結構的特征物理量，表示在離所研究原子距離為r處其他原子出現的概率，概率越大，說明這幾種原子間的相互作用力越強。通常情況分子間作用力包括氫鍵和范德華力，氫鍵長度為0.21～0.31 nm，強范德華力相互作用鍵長范圍為0.31～0.50 nm，弱范德華力相互作用鍵長大于0.50 nm。從徑向分布函數的峰值位置r可以判斷相互作用的類型，從峰值g（r）的高低進行推斷作用力的強弱［17］。g（r）計算式如方程（3）為：式中，NAB表示與中心粒子A 距離r到r+dr處粒子B 數目；ρ表示粒子B 的平均數目密度。

2.2.3 自由體積分數

體系內的自由體積為分子運動和擴散提供空間，MD 中自由體積表示如方程（4）：

式中，Vf表示自由體積，Vsp表示周期箱體積，Vo表示分子占據體積，單位均為nm3。

自由體積是指除去NC 分子和鈍感劑分子共同占據體積后體系的剩余體積，采用Atom Volumes &Surfaces 工具，選取0～0.2 nm 的探針半徑，對DBP/NC、NA/NC 體系的自由體積進行計算。因為各體系的總體積大小不一樣，為了能夠比較這DBP、NA 的相對運動空間，定義自由體積分數FFV［18］為自由體積與各周期箱體積的比值，如方程（5）：

3 結果與討論

3.1 擴散系數

3.1.1 DBP、NA 在純NC 體系中的擴散系數

為了探究溫度對DBP 和NA 在NC 體系中擴散性能的影響，主要取決于發(fā)射藥鈍感、儲存、使用的溫度環(huán)境，比如鈍感溫度65～85 ℃，冬天低溫0 ℃左右，常溫 是25 ℃，所 以DBP/NC、NA/NC 體 系 選 取5、25、65、85 ℃四個代表溫度進行MD 模擬，模擬得到的DBP、NA 的MSD曲線，如圖2 所示。

圖2 DBP 和NA 在NC 中不同溫度下的均方位移Fig.2 MSD curves of DBP and NA in NC at different temperatures

將圖2 的MSD曲線進行擬合得到的斜率代入方程（2）求得擴散系數D，結果如圖3 所示。從圖3可以看出，在相同溫度下，DBP 的擴散系數一直高于NA，表明高分子NA 的遷移率低于小分子DBP，這是由于鈍感劑分子的大小會影響擴散，分子越小，越容易在聚合物孔洞中跳躍，2 種鈍感劑分子的大小順序為NA＞DBP，這與擴散系數大小的順序在邏輯上是一致的，即DBP 分子越小，越容易擴散。此外，隨著溫度升高，DBP、NA 的擴散系數也隨著增加，并且DBP 的MSD值隨溫度的升高的增幅高于NA，尤其65 ℃、85 ℃高溫下，DBP 的擴散系數約是NA 的2 倍，比 如85 ℃時DBP 與NA 的 擴 散 系 數分別為3.42×10-11m2·s-1和1.11×10-11m2·s-1。如圖4所示，采用指數方程擬合擴散系數D與溫度t的關系，t越 高，D的 增 長 趨 勢 也 越 大，t無 窮 小 時，D接 近于0。因此，與DBP 相比，NA 具有較好的抗遷移特性，該特性在高溫時更為顯著。

圖3 不同溫度下DBP 和NA 在NC 中的擴散系數Fig.3 Diffusion coefficients of DBP and NA in NC at different temperatures

3.1.2 DBP、NA 在雙基發(fā)射藥體系中的擴散系數

NG 既是發(fā)射藥中的能量組分，又對NC 起到增塑的作用，提高加工性能，但NG 的存在也會對鈍感劑的擴散遷移產生影響。通常中小口徑武器一般用含有10%NG 左右的發(fā)射藥，所以選取10%，15% NG 含量的雙基發(fā)射藥體系研究NG 含量對鈍感劑擴散的影響具有實際意義。外加5%DBP 和5%NA 在雙基發(fā)射藥體系中的物理模型，分別為DBP/10NG/NC、 DBP/15NG/NC、 NA/10NG/NC、NA/15NG/NC 體系，由于溫度較多，僅選取25 ℃作為代表進行MD 模擬，得到的MSD曲線如圖5 所示。將圖5 均方位移隨時間的變化曲線，經擬合得到各直線斜率，代入方程（2）計算，得到的擴散系數結果如表1 所示。從表1 中可以看出，在相同溫度下，DBP、NA 的擴散系數均隨著NG 含量增加而增大，在含有15%NG 的雙基發(fā)射藥體系中運動最快，擴散系 數 分 別 為6.96×10-12m2·s-1和9.33×10-12m2·s-1。因此，提高NG 含量將會加速鈍感劑的擴散遷移，從而降低鈍感發(fā)射藥的使用壽命。

圖5 DBP、NA 在雙基發(fā)射藥體系中的均方位移Fig.5 MSD curves of DBP and NA in double base gun propellant

3.2 徑向分布函數

為 了 直 觀 地 探 討NC 與DBP、NC 與NA、NC 與NG 之間的相互作用，選取有可能形成氫鍵作用的NC（O、H、N）與DBP（O、H）、NA（O、H）、NG（H）原子對進行了徑向分布函數分析。根據原子的力場類型，NC硝基中的O 原子、H 原子、N 原子分別標記為o12、o2n、h1、n3o，NA、DBP 的O 原子和H 原子都標記為o1=和h1，NG 的H 原子標記為h1，如圖6 所示。

圖6 NC、DBP、NA、NG 的原子分布Fig.6 Atomic distrubutions of NC，DBP，NA and NG

3.2.1 DBP、NA 在純NC 體系的徑向分布函數

根據氫鍵形成的原理，可以推測NC 分子中的h1、o12、o2n、n3o 原子可能與DBP、NA 分子的o1=原子形成氫鍵，5 種原子對的徑向分布函數如圖7 所示。

從 圖7a 可 看 出，h1（NC）-o1=（DBP）原 子 對 在0.265 nm 附近出現峰，g（r）值為2.27；o2n（NC）-o1=（DBP）原子對在0.315 nm 附近出現峰，g（r）值為2.8；由圖7b 看出，h1（NC）-o1=（NA）原子對在0.275 nm附近出現峰，g（r）值為1.96；o2n（NC）-o1=（NA）原子對在0.355 nm 附近出現峰，g（r）值為1.41。表明NC與DBP、NC 與NA 之間存在較強的分子間作用力。

根據圖7a 和圖7b 的峰值位置r和g（r）的強度，與h1（NC）-o1=（NA）相比，h1（NC）-o1=（DBP）的g（r）峰值更大，表明NC 與DBP 之間的氫鍵作用比NC 與NA 強。h1（NC）-o1=（DBP）、o12（NC）-h1（DBP）、n3o（NC）-o1=（DBP）、h1（NC）-o1=（NA）、o12（NC）-h1（NA）原子對的g（r）最大時所對應的r值均在0.21～0.31 nm 范圍內，表明它們的相互作用力均屬于氫鍵作用力，其它原子對的g（r）最大時的r值均大于0.31 nm，表明它們的相互作用力均屬于范德華力。結果表明，NC 與DBP之間的氫鍵和范德華力均比NC 與NA 之間的大，意味著NC 與DBP 分子間的作用力比NC 與NA 之間的強，表明DBP 更易受到NC 的束縛，從而不利于DBP 擴散。雖然分子間作用力具有阻礙或者減弱擴散的作用，但是DBP 分子量小更容易擴散，通過比較3.1.1 節(jié)的擴散系數的結論可知：兩者共同的作用結果是分子量影響更大。

進一步分析在不同溫度下（5、25、65、85 ℃）NC與DBP、NA 之間的作用力，由于篇幅有限，因此僅選取NC（h1）與DBP（o1=）、NA（o1=）形成氫鍵作用力最大的徑向分布函數進行分析，見圖8。

如圖8 所示，隨著溫度的升高，NC 與DBP、NA 的原子間氫鍵作用的峰值變低，意味著NC 與DBP、NA的分子間作用力減弱，這說明高溫使NC 與鈍感劑之間的相互作用變弱，因此DBP、NA 擴散速度隨著溫度的升高而加快。

圖7 DBP/NC、NA/NC 模型的徑向分布函數Fig.7 Radial distribution functions of DBP/NC and NA/NC model

圖8 DBP/NC、NA/NC 體系在不同溫度下的徑向分布函數Fig.8 Radial distribution functions of DBP/NC and NA/NC model at different temperatures

3.2.2 DBP、NA 在雙基發(fā)射藥體系中的徑向分布函數

由于篇幅有限，因此僅選取形成氫鍵作用力最大的NC（h1、o12）與DBP（o1=）、NA（o1=）、NG（h1）原子對進行徑向分布函數分析，進一步探究DBP、NA 在不同NG 含量的雙基發(fā)射藥體系擴散能力不同的原因。由 圖9 可 見，h1（NC）-o1=（NA）、h（NC）-o1=（DBP）的峰值g（r）均隨著體系中的NG 含量增加而減小，這 是 由 于NG 的h1 原 子 與NC 的-ONO2基 團 的o12 原 子 形 成 氫 鍵，使 得NC 與DBP、NA 之 間 的 作 用力相對微弱很多。如圖10 所示，不同NG 含量的雙基發(fā)射藥體系的h1（NG）-o12（NC）的峰值g（r）均隨著NG 的含量增多而增強。由此推測，增塑劑NG 的添加減弱了NC 與DBP、NA 之間的作用力，因此DBP、NA的運動更活躍，擴散能力增強。

3.3 自由體積分數

通過分子動力學模擬分析，分別得到了5，25 ℃和65 ℃下DBP/NC、NA/NC 的自由體積分布，結果如圖11 所示。

Voccupied（Vo）是NC 和鈍感劑DBP、NA 共同所占的體積，Vfree（Vf）表示分子間的間隙，即自由體積。如圖11 所示，NC/DBP 的Vf均大于相同溫度下的NA/NC。此外，DBP/NC、NA/NC 的Vf值均隨溫度的升高而顯著增加，而Vo值均隨溫度的升高而略有減小。

圖9 DBP、NA 在雙基發(fā)射藥體系的徑向分布函數Fig.9 Radial distribution functions of DBP and NA in double base propellant

圖10 NG（h1）與NC（o12）在不同體系中的徑向分布函數Fig.10 Radial distribution functions of NG（h1）and NC（o12）in different systems

圖12 DBP/NA 、NA/NC 在不同溫度下的自由體積分數Fig.12 Fractional free volume of DBP/NA and NA/NC at different temperatures

體系的自由體積分數越大，提供給鈍感劑運動的有效活動空間就越大，擴散也就越容易發(fā)生。由圖12可以看出，相同溫度下，體系內自由體積分數的大小順序為DBP/NC＞NA/NC，可見DBP 的相對運動空間比NA 的大，這與擴散系數的順序在邏輯上是一致的；此外，隨著溫度升高，DBP/NC、NA/NC 體系的自由體積分數也在增大，意味著分子的有效活動空間隨著溫度的升高而增大，這是由于高溫使NC 的活動性增強、分子間距離增加，從而增加了體系的自由體積分數，這給DBP、NA 提供了更大的活動空間，這使DBP、NA 更易發(fā)生擴散。

4 結論

對DBP、NA 在NC 體系中擴散進行了分子動力學研究。研究結果表明：

（1）在相同溫度下，鈍感劑在發(fā)射藥中的擴散系數大小順序為：DBP＞NA，這說明NA 具有較強的抗遷移能力；此外，在25 ℃下，DBP與NA 的擴散系數分別為1.13×10-11m2·s-1和5.13×10-12m2·s-1，但在高溫85 ℃下，DBP 與NA 的擴散系數分別為3.42×10-11m2·s-1和1.11×10-11m2·s-1，因此NA 相較于DBP 具有的抗遷移特性在高溫時更為凸顯。

（2）從NC 與DBP、NA 的分子間相互作用、體系的自由體積分數2 個方面分析了鈍感劑擴散性能的微觀機理。從微觀角度解釋擴散系數隨著溫度的升高而增大的原因為：高溫導致NC 與DBP、NA 的分子間作用力減弱，DBP、NA 的運動變得更加活躍，同時體系的自由體積分數也在增大，增大了分子運動的有效活動空間，這均是DBP、NA 擴散能力增強的原因。

（3）在25 ℃下，DBP、NA 的擴散系數均隨著NG含量的增加而增加，當雙基發(fā)射藥體系中NG 含量分別 為0%、10%、15%，DBP 的 擴 散 系 數 分 別 為7.03×10-12，8.05×10-12，1.41×10-11m2·s-1；NA 的擴散系數分別為4.05×10-12，6.96×10-12，9.33×10-12m2·s-1。這是由于NG 含量的增加使得NC 與DBP、NA 之間的分子作用力均減小，DBP、NA 的運動變得更加活躍，擴散能力增強。