PEG/增塑劑共混物相容性的分子動(dòng)力學(xué)模擬和介觀模擬

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

（1.海軍航空大學(xué)岸防兵學(xué)院，山東 煙臺(tái) 264001；2.海軍航空大學(xué)航空基礎(chǔ)學(xué)院，山東 煙臺(tái) 264001）

1 引言

增塑劑用于改善推進(jìn)劑力學(xué)性能，并提供一定的能量，如果增塑劑和粘合劑的相容性不好，將會(huì)導(dǎo)致增塑劑遷移，引起推進(jìn)劑力學(xué)性能變差、內(nèi)部能量流失、界面脫粘、襯層絕熱層不耐燒蝕［1-2］、甚至還會(huì)影響推進(jìn)劑的安定性和感度［3］，所以，增塑劑與其他組分的相容性是評(píng)價(jià)推進(jìn)劑性能的首要指標(biāo)之一。隨著計(jì)算機(jī)軟硬件的發(fā)展，很多研究者采用分子模擬法探究混合體系的相容性［4-5］，此法較實(shí)驗(yàn)法簡(jiǎn)單易行［6］，還可以從微觀角度研究相容性的本質(zhì)，但多數(shù)研究集中于HTPB 推進(jìn)劑［7-9］。

高能固體推進(jìn)劑NEPE常用的增塑劑為硝化甘油（NG），因其感度較高導(dǎo)致推進(jìn)劑的危險(xiǎn)品等級(jí)為1.1級(jí)［10］，并且存在一定程度的增塑劑遷移現(xiàn)象［11］。有研究表明：含能鈍感增塑劑TMETN可與NG起到類似的增塑作用，且TMETN感度更低［12］；法國(guó)已應(yīng)用TMETN降低了雙基推進(jìn)劑的感度［13］；美國(guó)也已將低能硝酸酯增塑劑［14］用于戰(zhàn)術(shù)助推的鈍感NEPE推進(jìn)劑；日本 Yoshio Oyumi［15］發(fā)現(xiàn) TMETN增塑的GAP/AN推進(jìn)劑具有良好力學(xué)性能，TMETN有助于提高燃速。而國(guó)內(nèi)關(guān)于降低NEPE推進(jìn)劑感度的研究還較少，特別是將TMETN應(yīng)用于戰(zhàn)術(shù)NEPE推進(jìn)劑的研究則更少。

本文立足于以上背景，從相容性能的角度研究了TMETN應(yīng)用于NEPE推進(jìn)劑的可行性。采用分子動(dòng)力學(xué)模擬得到溶度參數(shù)和徑向分布函數(shù)，對(duì)聚乙二醇（PEG）與增塑劑NG、BTTN、TMETN的相容性進(jìn)行了預(yù)測(cè)和比較，并通過(guò)結(jié)合能分析相容性的本質(zhì)；然后把微觀結(jié)果轉(zhuǎn)化為介觀模擬的輸入?yún)?shù)，對(duì)比了NG/PEG、BTTN/PEG、TMETN/PEG三種體系的介觀形態(tài)隨時(shí)間演變過(guò)程，為低易損NEPE推進(jìn)劑的研制提供依據(jù)。

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

2.1 建模過(guò)程及分子動(dòng)力學(xué)模擬細(xì)節(jié)

利用Material Studio中Visualizer模塊建立分子量為6002的PEG分子及NG、BTTN、TMETN分子（分子結(jié)構(gòu)式見(jiàn)圖1），利用Forcite模塊對(duì)分子進(jìn)行幾何優(yōu)化，幾何優(yōu)化是為了尋找局部最優(yōu)構(gòu)型，優(yōu)化過(guò)程均選擇Compass力場(chǎng)［16］；然后進(jìn)行若干次退火處理、直到分子的能量穩(wěn)定在最小值，退火可以使分子跨過(guò)局部能壘、尋找全局最優(yōu)構(gòu)型，退火條件為初始溫度為298 K，中間溫度為 598 K［17］，每次退火后自動(dòng)進(jìn)行幾何優(yōu)化，這是在全局最優(yōu)構(gòu)型中進(jìn)一步尋找，從中選取能量最低的一幀構(gòu)型作為輸出的最優(yōu)單鏈構(gòu)型，為下一步建立無(wú)定形結(jié)構(gòu)做好準(zhǔn)備。

圖1 PEG、NG、BTTN、TMETN的分子結(jié)構(gòu)式Fig.1 Chemical structures of PEG、NG、BTTN、TMETN

2.1.1 純物質(zhì)體系模擬細(xì)節(jié)

為減小尺寸效應(yīng)，確保無(wú)定形分子模型中有1000個(gè)以上原子［7］，PEG、NG、BTTN、TMETN 無(wú)定形分子模型中包含的分子數(shù)量及物質(zhì)的密度見(jiàn)表1。

表1 PEG、NG、BTTN、TMETN 4種純物質(zhì)無(wú)定形分子模型參數(shù)Table 1 The parameters of pure substance（PEG、NG、BTTN、TMETN）in amorphous cells

每種物質(zhì)均采用Amorphous cell模塊構(gòu)建10個(gè)構(gòu)型，從中篩選能量最低者，進(jìn)行幾何優(yōu)化和若干次退火處理，為保證純物質(zhì)密度保持為設(shè)定值，采用微正則系綜（NVE）［18］進(jìn)行退火，退火后進(jìn)行幾何優(yōu)化，從若干禎中篩選能量最低者進(jìn)行分子動(dòng)力學(xué)模擬［17］。

動(dòng)力學(xué)模擬用Andersen控溫法［19］，范德華和靜電力的計(jì)算分別用 Atom-based［20］和 Ewald 方法，非鍵截?cái)喟霃綖?0.95 nm，時(shí)間步長(zhǎng) 1 fs［7］，其它為默認(rèn)，先進(jìn)行300p正則系綜（NVT）的動(dòng)力學(xué)松弛計(jì)算，再進(jìn)行400 ps的NVE系綜動(dòng)力學(xué)計(jì)算，后50 ps已經(jīng)平衡、用于采集溶度參數(shù)及徑向分布函數(shù)數(shù)據(jù)；此外，對(duì)PEG進(jìn)行了1 ns的模擬，1 ns得到的溶度參數(shù)20.809與400 ps的結(jié)果20.941比較接近，證實(shí)了400 ps模擬時(shí)間確實(shí)已足夠長(zhǎng)。

2.1.2 混合體系模擬細(xì)節(jié)

在保證周期箱內(nèi)原子數(shù)量大于1000的基礎(chǔ)上，為了便于比較增塑劑分子與PEG分子間相互作用，將三種增塑劑分子的數(shù)量均取50個(gè)，周期箱的初始密度按體積比例進(jìn)行加和［7］，NG/PEG、BTTN/PEG、TMETN/PEG三種混合體系的具體參數(shù)設(shè)置如表2，無(wú)定形分子模型（圖2）的構(gòu)造和優(yōu)化與2.1.1節(jié)純物質(zhì)的類似。

表2 混合物無(wú)定形分子模型參數(shù)Table 2 The parameters of blends in amorphous cells

無(wú)定形分子模型要選取重復(fù)出現(xiàn)且能量較低的構(gòu)型，退火后進(jìn)行400 ps NVT系綜的動(dòng)力學(xué)松弛，再采用Berendsen控壓方法［7］在等溫等壓系綜中（NPT）進(jìn)行400 ps的動(dòng)力學(xué)計(jì)算，后50 ps用于分析性能。為了證明模擬的時(shí)間足夠長(zhǎng)，又對(duì)NG/PEG體系進(jìn)行了1 ns NPT系綜計(jì)算，不同模擬時(shí)間的徑向分布函數(shù)結(jié)果如圖3，因徑向分布函數(shù)只是分析曲線的相對(duì)位置高低、進(jìn)而判斷分子間作用的相對(duì)大小，所以400 ps的結(jié)果與1 ns的結(jié)果很相近、可以用于性能的分析。

2.2 介觀建模過(guò)程及模擬細(xì)節(jié)

介觀模擬需要確定：每種物質(zhì)的粗?；咚规湺螖?shù)量，及不同粗?；湺伍g的相互作用參數(shù)。

高分子預(yù)聚物的高斯鏈段數(shù)量采用式（1）［8］：

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

采用Synthia模塊計(jì)算預(yù)聚物PEG的C∞為4.98，粗?；樽拥馁|(zhì)量為219.12（用特征比乘以單體分子質(zhì)量），所以需采用27個(gè)粗粒化珠子表示分子量為6002的PEG分子，根據(jù)各粗?；樽拥馁|(zhì)量相同，采用一個(gè)珠子代替一個(gè)增塑劑分子。

此外，不同珠子間的相互作用參數(shù)為：

圖2 NG/PEG、BTTN/PEG、TMETN/PEG混合物的無(wú)定形分子模型Fig.2 Amorphous cells of the NG/PEG、BTTN/PEG、TMETN/PEG

式中，ν-1εij為 MesoDyn 相互作用參數(shù) kJ·mol-1；χij為Flory-Huggins參數(shù)（無(wú)量綱）；R為氣體常數(shù)8.314 J·（mol·K）-1，T取 298 K。

在確定分子粗粒化珠子數(shù)量及相互作用參數(shù)后，在Meso Dyn模塊中進(jìn)行模擬：根據(jù)增塑劑與粘合劑預(yù)聚物的質(zhì)量比（增塑比）為3∶1，設(shè)置珠子的混合比例，三維周期箱尺寸為32.0 nm×32.0 nm×32.0 nm，溫度為298 K，步長(zhǎng)為20 ns，總時(shí)間為1000 μs，其它為默認(rèn)。

3 結(jié)果與討論

3.1 MD模擬結(jié)果分析

3.1.1 相容性的溶度參數(shù)判斷

圖3 不同模擬時(shí)長(zhǎng)的NG/PEG混合物中的分子間徑向分布函數(shù)Fig.3 Intermolecular radial distribution function for NG/PEG with different simulation time

根據(jù)高分子溶液理論［21］，恒溫和恒壓下的溶解過(guò)程自發(fā)進(jìn)行的前提是：Gibbs自由能的變化小于零。兩組分的溶度參數(shù)差值越小，越容易滿足此條件、相容性越好，這就是物質(zhì)的“相似相溶”原理，也可以通過(guò)幾種物質(zhì)的極性相近、分子間作用力類型和大小相近、各自的徑向分布函數(shù)曲線相近，來(lái)評(píng)價(jià)相容性的優(yōu)劣。

采用分子動(dòng)力學(xué)模擬計(jì)算得到PEG和增塑劑的溶度參數(shù)結(jié)果見(jiàn)表3，δ為本次結(jié)算結(jié)果，為了便于比較，同時(shí)將文獻(xiàn)［22-23］的分子動(dòng)力學(xué)計(jì)算結(jié)果δ′和δ"、文獻(xiàn)［24］采用 Hoftyzer-Krevelen（1976）方法的估算值δ′"均列于表3。由表3可見(jiàn)，本次模擬（δ）與文獻(xiàn)模擬值（δ′）、（δ"）相近，但與估算值（δ′"）有一定差距，分析原因?yàn)樗捎玫挠?jì)算方法不同，而本模擬是在相同力場(chǎng)和參數(shù)設(shè)置下進(jìn)行的，所以可以考察溶度參數(shù)的相對(duì)值［4］；并且不同方法得到物質(zhì)溶度參數(shù)的大小順序是相同的，即通過(guò)計(jì)算Δδ判斷共混物相容性優(yōu)劣的結(jié)果是一致的，由表4可見(jiàn)，PEG與三種增塑劑的相容性優(yōu)劣順序?yàn)椋篢METN/PEG＞BTTN/PEG＞NG/PEG。

四種物質(zhì)的內(nèi)聚能密度和溶度參數(shù)的分量結(jié)果見(jiàn)表5。分析表5結(jié)果發(fā)現(xiàn)，在分子間相互作用中，每種物質(zhì)內(nèi)聚能密度的范德華分量（CEDv）均大于靜電分量（CEDe），即范德華作用均大于靜電作用、占據(jù)主要地位，特別是對(duì)于聚乙二醇而言；因?yàn)榉菢O性分子之間只存在范德華作用，所以范德華力代表著非極性作用、而靜電力代表著極性作用［22］，由表5可見(jiàn)，PEG的極性較弱，TMETN與PEG的極性更相近、相容性更好。

表3 溶度參數(shù)模擬值與文獻(xiàn)值Table 3 Simulated values and literature values of solubility parameters

表4 3種增塑劑與PEG的溶度參數(shù)的差值Table 4 The difference of solubility parameters between plasticizers and PEG （J·cm-3）1/2

表5 內(nèi)聚能密度及其范德華分量、靜電分量（CED=CEDv+CEDe）Table 5 Cohesive energy density and the van der Waals、electrostatic components（CED=CEDv+CEDe）J·m-3

此外，運(yùn)用公式（3）可以將表4所示的溶度參數(shù)差值轉(zhuǎn)換為Flory-Huggins參數(shù)（χ）：

式中，Δδ為兩種物質(zhì)溶度參數(shù)的差值，（J·cm-3）1/2；Vr為參比體積，取兩種珠子摩爾體積的平均值，cm3；R為氣體常數(shù) 8.314 J·（mol·K）-1；T為 298 K。三種增塑劑與PEG的Flory-Huggins參數(shù)結(jié)果列于表6，再利用公式（2）將其轉(zhuǎn)換為介觀動(dòng)力學(xué)模擬的輸入?yún)?shù)（ν-1εij），結(jié)果列于表 7。

3.1.2 相容性的徑向分布函數(shù)判斷

圖4為四種純物質(zhì)的分子內(nèi)、分子間徑向分布函數(shù)，分子內(nèi)徑向分布函數(shù)可體現(xiàn)物質(zhì)的有序度信息［25］：圖 4a中NG-NG、BTTN-BTTN兩條曲線進(jìn)行了向上平移10個(gè)單位處理，目的在于方便觀察，可見(jiàn)PEG與三種增塑劑均屬于無(wú)定性結(jié)構(gòu)（僅在0.3 nm以前出現(xiàn)峰）［26］。

表6 增塑劑/PEG的Flory-Huggins參數(shù)Table 6 Flory-Huggins parameters of plasticizer/PEG

表7 混合體系粗粒化珠子的相互作用參數(shù)Table 7 Exclusion parameters between different beads

圖4 純物質(zhì)分子內(nèi)和分子間徑向分布函數(shù)Fig.4 Intramolecular and intermolecular radial distribution functions for pure substance

分子間徑向分布函數(shù)可以反映分子間相互作用的類型（峰值位置）及大?。ㄇ€高低）［27］，氫鍵作用范圍為0.26～0.31 nm，范德華作用范圍為0.31～0.50 nm。由圖4b可見(jiàn)，四種純物質(zhì)內(nèi)的分子間作用均為范德華作用，此外，根據(jù)“純物質(zhì)的分子間徑向分布函數(shù)越接近，其混合后的相容性越好［25］”，也可得出體系相容性優(yōu)劣順序?yàn)椋篢METN/PEG＞BTTN/PEG＞NG/PEG；通過(guò)分析三種增塑劑分子結(jié)構(gòu)發(fā)現(xiàn)，分子內(nèi)均含有3個(gè)硝基，相比于NG分子，BTTN分子中僅多了一個(gè)亞甲基，根據(jù)“結(jié)構(gòu)決定性質(zhì)”，可以說(shuō)明亞甲基削弱了BTTN分子間的相互作用、降低了溶度參數(shù)值，而TMETN分子是一個(gè)中心碳原子周圍連接四部分官能團(tuán)，其類四面體結(jié)構(gòu)導(dǎo)致分子間作用力最小、溶度參數(shù)與PEG最接近，相容性最好。

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

進(jìn)一步分析混合物分子間的徑向分布函數(shù)，由圖5可見(jiàn)，增塑劑與PEG的相互作用方式均為VdW作用；有文獻(xiàn)［28-30］記載，在混合物中，兩種分子間的徑向分布函數(shù)越是高于單一物質(zhì)自身的分子間徑向分布函數(shù)，混合物的相容性越好。在NG/PEG體系中，NG與NG分子間作用明顯大于NG與PEG的分子間作用，所以體系容易發(fā)生相分離、相容性最差；在BTTN/PEG體系中，當(dāng)r＜0.4 nm時(shí)，BTTN與BTTN的分子間作用略大于BTTN與PEG的分子間作用，當(dāng)r＞0.4 nm時(shí)，前者明顯大于后者，所以BTTN/PEG體系相容性較差；在TMETN/PEG體系中，TMETN與PEG分子間作用絕大多數(shù)大于TMETN與TMETN分子間作用，說(shuō)明TMETN與PEG最難發(fā)生相分離、相容性最好，這與溶度參數(shù)的判斷結(jié)果一致。

3.1.3 相容性的本質(zhì)分析

相容的本質(zhì)就是組分間的相互作用，結(jié)合能可以表示分子間相互作用的類型和大?。?1］，其表達(dá)式為：

式中，EPEG/plasticizer為混合體系的總能量，Eplasticizer和EPEG分別為增塑劑和PEG的平均單點(diǎn)能。

混合體系內(nèi)分子間結(jié)合能的結(jié)果見(jiàn)表8。由表8可見(jiàn)，在增塑劑/PEG體系中，兩種分子間的相互作用僅來(lái)源于非鍵作用，在NG/PEG體系中范德華力約占61.64%，靜電力約占38.36%，在BTTN/PEG體系中范德華力約占69.11%，靜電力約占30.89%，在TMETN/PEG體系中范德華力約占72.67%，靜電力約占27.33%，這說(shuō)明混合物分子相容的本質(zhì)為：“分子間主要通過(guò)范德華作用相結(jié)合，并且，增塑劑與PEG的相容性越好，分子間通過(guò)非極性作用（范德華作用）相結(jié)合的占比越大”。

3.2 Meso Dyn模擬結(jié)果分析

分子動(dòng)力學(xué)模擬反映的是微觀角度的信息，而高分子預(yù)聚物的相分離特征區(qū)尺寸和時(shí)間都遠(yuǎn)大于其模擬范圍，下面進(jìn)一步借助介觀模擬，忽略原子尺度的信息，得到三種體系的等密度圖在1000 μs內(nèi)的變化規(guī)律，如圖6所示，為了便于觀察，圖中只顯示了增塑劑相。

由圖6可見(jiàn)，在模擬初期，三種增塑劑都均勻分散在PEG相中，隨著時(shí)間延長(zhǎng)，NG/PEG體系中NG相開(kāi)始聚集，最終形成較大的顆粒相分散于體系之中；相比之下，BTTN/PEG體系相分離發(fā)生時(shí)間更晚、程度更??；而TMETN/PEG體系中PEG僅發(fā)生輕微同相歸并，未出現(xiàn)明顯相分離。

表8 混合體系分子間的結(jié)合能及其源項(xiàng)Table 8 The binding energy between molecules and the components kJ·mol-1

圖6 體系中增塑劑相NG、BTTN、TMETN的等密度圖隨時(shí)間變化過(guò)程Fig.6 The isosurface of the density fields for plasticizer NG，BTTN和TMETN with time evolution in blends

體系內(nèi)硝酸酯的有序度參數(shù)隨時(shí)間變化曲線如圖 7所示。由圖7可見(jiàn)，在模擬初期NG的有序度參數(shù)迅速增大，50 μs后增大速率變慢，說(shuō)明體系的熵已明顯變大、發(fā)生了相分離，在300 μs后有序度趨于不變，說(shuō)明體系處于新的穩(wěn)定狀態(tài)；BTTN有序度參數(shù)增大的速率和程度均低于NG，在200 μs后處于新的穩(wěn)定狀態(tài)，說(shuō)明BTTN發(fā)生相分離的程度低于NG；相比于前兩者，TMETN的有序度參數(shù)基本穩(wěn)定在0附近，體系熵變化很小、基本處于穩(wěn)定狀態(tài)。從體系自由能密度變化也可以看出：NG/PEG體系的自由能密度迅速下降；BTTN/PEG體系的自由能密度稍有下降；而TMETN/PEG體系的自由能密度基本穩(wěn)定，可見(jiàn)TMETN/PEG體系的相容性最好，其次為BTTN/PEG，最差為NG/PEG體系。

圖7 有序度參數(shù)和自由能密度隨模擬時(shí)間變化圖Fig.7 Order parameter and free energy density for plasticizer with time evolution

4 結(jié)論

通過(guò)對(duì)PEG/NG，PEG/BTTN，PEG/TMETN共混體系的分子動(dòng)力學(xué)和介觀動(dòng)力學(xué)模擬，主要得到以下結(jié)論：

（1）分子動(dòng)力學(xué)模擬PEG和三種增塑劑的溶度參數(shù)、純物質(zhì)及混合物的分子間徑向分布函數(shù)，結(jié)果均表明體系相容性的優(yōu)劣順序?yàn)椋篢METN/PEG＞BTTN/PEG＞NG/PEG，討論其微觀原因?yàn)椋築TTN中的亞甲基以及TMETN的分子結(jié)構(gòu)削弱了各自的分子間相互作用，減小了與PEG的溶度參數(shù)差值。

（2）對(duì)混合體系的分子間徑向分布函數(shù)和結(jié)合能分析得到：PEG與三種增塑劑相容本質(zhì)為“通過(guò)非鍵作用相結(jié)合，其中范德華作用占主要比重”，此外，增塑劑的極性越接近于PEG，體系通過(guò)范德華作用相結(jié)合的比重越大。

（3）增塑劑/PEG體系介觀形態(tài)演變、等密度圖和有序度參數(shù)圖表明：TMETN/PEG體系未發(fā)生相分離、熵變最小，混溶性最好；BTTN/PEG體系發(fā)生較輕微相分離、混溶性較好；而NG/PEG體系發(fā)生相分離、混溶性較差。

（4）MD和MesoDyn模擬結(jié)果表明：含能鈍感增塑劑TMETN與粘合劑預(yù)聚物PEG的相容性優(yōu)于增塑劑NG和BTTN，同時(shí)TMETN的感度遠(yuǎn)低于NG，可以考慮將其代替或部分代替NG、用于改善NEPE推進(jìn)劑的易損性，研制低易損戰(zhàn)術(shù)武器裝備。