李華雷,楊 斌,李云龍,,王世杰

(1.汕頭大學(xué) 工學(xué)院土木與環(huán)境工程系, 廣東 汕頭 515063; 2.沈陽(yáng)工業(yè)大學(xué) 機(jī)械工程學(xué)院, 遼寧 沈陽(yáng) 110870)

1 前 言

自1991年碳納米管(CNT)被日本科學(xué)家發(fā)現(xiàn)以來(lái),其優(yōu)異的電學(xué)、力學(xué)和熱學(xué)性能引起了廣泛關(guān)注[1]。超高的楊氏模量和抗拉強(qiáng)度使其成為高性能聚合物基復(fù)合材料的首選增強(qiáng)填料[2-5],在CNT 增強(qiáng)聚合物復(fù)合材料的力學(xué)性能方面,學(xué)者們進(jìn)行了廣泛的實(shí)驗(yàn)研究。汪華鋒等[6]利用催化裂解制備的CNT/環(huán)氧樹(shù)脂復(fù)合材料,力學(xué)性能測(cè)試表明:CNT 可以使環(huán)氧樹(shù)脂基體材料的拉伸及壓縮強(qiáng)度分別提高284%和122%,同時(shí)也可增加環(huán)氧樹(shù)脂材料的韌性。蓬勃等[7]制備了質(zhì)量分?jǐn)?shù)不同的CNT/聚乙烯(PE)復(fù)合材料顯示:當(dāng)填充0.5 wt%的CNT 時(shí),斷裂伸長(zhǎng)率達(dá)到最大值429.2%,再進(jìn)一步增加CNT 填充量,其斷裂伸長(zhǎng)率呈單調(diào)下降趨勢(shì)。由于受到時(shí)間和空間尺度的限制,在納米尺度下材料間的交互過(guò)程和作用機(jī)制(如CNT 和基體之間的界面相互作用)等問(wèn)題的探索在實(shí)驗(yàn)中仍然難以有效開(kāi)展。分子動(dòng)力學(xué)(MD)模擬作為近年來(lái)新興的研究手段,可以在分子層面模擬納米復(fù)合材料基體內(nèi)的交互過(guò)程和作用機(jī)制[8-9],成為彌補(bǔ)實(shí)驗(yàn)研究局限性的有力補(bǔ)充。Frankland等[10]采用MD方法研究了(10, 10)單壁CNT(SWCNT)增強(qiáng)PE 基復(fù)合材料的應(yīng)力-應(yīng)變曲線。結(jié)果表明,長(zhǎng)徑比為10的長(zhǎng)SWCNT 纖維增強(qiáng)的納米復(fù)合材料的剛度顯著增加,對(duì)于長(zhǎng)徑比為4的CNT 增強(qiáng)體,則沒(méi)有觀察到有顯著的增強(qiáng)。Li等[11]采用MD模擬方法研究了具有預(yù)定單邊裂紋的CNT/環(huán)氧樹(shù)脂復(fù)合材料的斷裂行為。結(jié)果表明,引入CNT后,環(huán)氧復(fù)合材料的拉伸強(qiáng)度和斷裂伸長(zhǎng)率分別提高了24.8%和34.3%。Mokashi等[12]對(duì)(10, 10)CNT 增強(qiáng)非晶態(tài)PE聚合物基體的楊氏模量和拉伸強(qiáng)度進(jìn)行了MD模擬研究,納米復(fù)合材料的模擬單元尺寸設(shè)置為4.5 nm×4.5 nm×4.4 nm。模擬結(jié)果表明,CNT/聚合物復(fù)合材料的縱向楊氏模量為82 GPa,是純非晶聚乙烯的25倍左右。

雖然MD 已經(jīng)廣泛應(yīng)用于納米復(fù)合材料領(lǐng)域,但MD 模擬計(jì)算的力場(chǎng)復(fù)雜,致使工作量巨大,效率較低,使其難以進(jìn)行更大時(shí)間和空間尺度的模擬。為了克服該問(wèn)題,粗?；?CG)方法得到了發(fā)展[13-15]。CG的原理是將一組原子看作一個(gè)整體,忽略內(nèi)部信息,然后映射成對(duì)應(yīng)的粗?；樽?從而保持原子模型的分子細(xì)節(jié)的同時(shí),成倍縮減全原子體系的自由度和粒子數(shù),增大模擬的時(shí)間和空間尺度。目前,已有文獻(xiàn)進(jìn)行了聚合物材料的粗?；Ｐ图捌淞?chǎng)參數(shù)的初步開(kāi)發(fā)工作[16-18]。Marrink 課題組用粗?；椒ǖ玫搅薓ARTINI力場(chǎng)[19],主要用于研究生物大分子的動(dòng)態(tài)行為,比如磷脂分子、蛋白質(zhì)及DNA 等。蔡莊立等[20]根據(jù)特定的映射方案建立了PE 的粗?；Ｐ?并以全原子模型的鍵長(zhǎng)、鍵角分布曲線來(lái)獲取PE 的粗粒化力場(chǎng)參數(shù)。文中驗(yàn)證了粗?；瘏?shù)模擬PE 模型的靜態(tài)結(jié)構(gòu)分布,且粗粒化模型的導(dǎo)熱系數(shù)的模擬值與實(shí)驗(yàn)值符合較好。Langeloth等[21]利用玻爾茲曼反轉(zhuǎn)迭代法得到了環(huán)氧樹(shù)脂的粗?；Ｐ鸵约傲?chǎng)參數(shù),并利用此參數(shù)在更大的尺度下研究了環(huán)氧樹(shù)脂的玻璃化轉(zhuǎn)變溫度與交聯(lián)程度的變化關(guān)系。

然而由于目前的粗?；芯咳鄙賹?duì)粗粒化力場(chǎng)參數(shù)準(zhǔn)確性的驗(yàn)證,粗?；?chǎng)參數(shù)與模型物理性能之間的關(guān)系亦鮮有討論。同時(shí),已有的玻爾茲曼反轉(zhuǎn)迭代法[22-23]等方法操作復(fù)雜,適用范圍有限。本研究使用基于應(yīng)變能守恒定律的恒應(yīng)變方法獲得聚酰亞胺(PI)和CNT 的力場(chǎng)參數(shù),然后通過(guò)不同模擬實(shí)驗(yàn)來(lái)驗(yàn)證粗?；瘏?shù)的準(zhǔn)確性,最后對(duì)粗粒化參數(shù)與物理性能的關(guān)系將進(jìn)行詳細(xì)討論。

2 建模與粗粒化參數(shù)的獲取

粗?；?chǎng)參數(shù)主要由成鍵勢(shì)函數(shù)和非鍵勢(shì)函數(shù)兩部分構(gòu)成。體系總勢(shì)能Etotal可以寫(xiě)成鍵長(zhǎng)能Eb、鍵角能Ea、范德華能Evdw和系統(tǒng)恒定自由能U0之和,表達(dá)式形式如下:

式中組成總勢(shì)能的各個(gè)勢(shì)函數(shù)函數(shù)形式如下:

式中:kd和d0是鍵長(zhǎng)的剛度系數(shù)和平衡鍵長(zhǎng);kθ和θ0是鍵角的剛度系數(shù)和平衡鍵角;D0和r0分別是Lennard—Jones(LJ)勢(shì)函數(shù)的勢(shì)阱參數(shù)和平衡距離。全原子模擬使用COMPASS力場(chǎng)來(lái)描述分子內(nèi)的相互作用力,非鍵相互作用的截?cái)喟霃饺?.2 nm,能量最小化方法使用共軛梯度法。

2.1 粗?；I長(zhǎng)參數(shù)的獲取

粗?；谝徊绞歉鶕?jù)設(shè)計(jì)好的映射方案將原子團(tuán)映射成粗?；樽印Ｈ鐖D1所示,組成PI的每個(gè)單體被映射為兩個(gè)粗?；樽覣(C8H3NO2)和一個(gè)粗粒化珠子B(C6H4),將(5, 5)的CNT 中每50個(gè)碳原子映射成1個(gè)粗?；樽覥(C50),以原子團(tuán)的質(zhì)量中心為粗?；樽拥闹行?。A 珠子的相對(duì)原子質(zhì)量為145.21,半徑為2.61 ?;B 珠子的相對(duì)原子質(zhì)量為70.09,半徑為1.94 ?;C 珠子的相對(duì)原子質(zhì)量為600.55,半徑為3.96 ?。

圖1 PI和CNT 的粗?；成浞桨窮ig.1 Coarse-grained mapping scheme of polyimide (PI)and carbon nanotube (CNT)

接下來(lái)使用恒應(yīng)變方法獲取粗?；?chǎng)參數(shù)。首先改變A-A 珠子對(duì)應(yīng)的原子團(tuán)之間的鍵長(zhǎng),通過(guò)分子模擬得到體系的總勢(shì)能U,根據(jù)應(yīng)變能守恒,粗?；w系的總勢(shì)能也等于U:

得到d0=6.78 ?,圖2(a)展示了U與Δd的變化曲線。kd由式(6)求出:

圖2 原子團(tuán)勢(shì)能與鍵長(zhǎng)增量的變化 (a)A-A 原子團(tuán);(b)A-B原子團(tuán);(c)C-C原子團(tuán)Fig.2 Variation of potential energy versus deformation of bond between (a) A-A group; (b) A-B group; (c) C-C group

得出kd=400.56 kcal·mol-1·?-2。

同樣地,改變A-B珠子對(duì)應(yīng)的原子團(tuán)之間的鍵長(zhǎng),得出U與Δd的變化曲線。如圖2(b)所示,根據(jù)式(5)、(6),依次得到d0=4.80 ?,kd=340.78 kcal·mol-1·?-2。改變C-C 珠子對(duì)應(yīng)的原子團(tuán)之間的鍵長(zhǎng),得出U與Δd的變化曲線,如圖2(c)所示。根據(jù)式(5)、(6)得出:d0=6.12 ?,kd=4 542.0 kcal·mol-1·?-2。

2.2 粗粒化鍵角參數(shù)的獲取

改變A-A-B珠子對(duì)應(yīng)的原子團(tuán)的角度,通過(guò)分子模擬得到體系的總勢(shì)能U,根據(jù)應(yīng)變能守恒,粗?；w系的總勢(shì)能也等于U:

得到θ0=158.8°,整理出U與Δθ的變化曲線,如圖3(a)所示。kθ由式(8)求出:

圖3 原子團(tuán)勢(shì)能與鍵角增量的變化 (a)A-A-B原子團(tuán);(b)A-B-A 原子團(tuán);(c)C-C-C原子團(tuán)Fig.3 Variation of potential energy versus the bond angle’s deformation of (a) A-A-B group; (b) A-B-A group; (c) C-C-C group

得出kθ=980.80 kcal·mol-1·rad-2。

同樣地,改變A-B-A 珠子對(duì)應(yīng)的原子團(tuán)的角度,得到U與Δθ的變化曲線,如圖3(b)所示。根據(jù)式(7)、(8)得出θ0=179.4°,kθ=670.21 kcal·mol-1·rad-2。最后,改變C-C-C珠子對(duì)應(yīng)的原子團(tuán)的角度,得出U與Δθ的變化曲線,如圖3(c)所示。根據(jù)式(7)、(8)得出θ0=180.0°,kθ=66 880.0 kcal·mol-1·rad-2。

2.3 粗?；兜氯A參數(shù)的獲取

改變A-A 珠子對(duì)應(yīng)的原子團(tuán)間的距離r,得到內(nèi)聚能U隨r的變化曲線,如圖4所示。U是體系范德華能與單個(gè)原子團(tuán)的范德華能之間的差值。U的最小值作為D0,U的最低點(diǎn)對(duì)應(yīng)的橫坐標(biāo)為r0,得出A 珠子的D0=7.56 kcal·mol-1,r0=4.01 ?。以同樣地步驟獲取其余珠子的D0和r0,整理成表1。

表1 A、B、C珠子的勢(shì)阱參數(shù)D0 與平衡距離r0Table 1 Well depth (D0) and equilibrium distance(r0) of beads A, B and C

圖4 A-A 原子團(tuán)的勢(shì)能與距離的變化Fig.4 Variation of energy versus distance between A-A group

利用得到的粗粒化參數(shù)對(duì)CNT/PI復(fù)合材料的粗?；Ｐ瓦M(jìn)行充分弛豫,然后將弛豫后的粗粒化模型的結(jié)構(gòu)分布曲線,包括鍵長(zhǎng)、鍵角以及徑向分布曲線,與全原子模型的進(jìn)行對(duì)比,并根據(jù)分布曲線的誤差修正粗?；?chǎng)參數(shù),最后得出調(diào)整后的參數(shù),如表2所示。

表2 修正后的粗?；瘏?shù)Table 2 Modified coarse-grained parameters

表3 PI鏈與CNT的楊氏模量Table 3 Young's modulus of PI chain and CNT

3 粗粒化參數(shù)的驗(yàn)證和探討

3.1 驗(yàn)證PI分子鏈和CNT的拉伸性質(zhì)

如圖5(a)所示,一條PI分子鏈含有四個(gè)單體,一條PI珠子鏈含有8個(gè)A 珠子,4個(gè)B珠子,分子鏈的長(zhǎng)度為66.17 ?,對(duì)兩條鏈進(jìn)行拉伸,應(yīng)變率為8%,得到圖6(a)中的應(yīng)力-應(yīng)變曲線。對(duì)曲線進(jìn)行擬合得到全原子與粗?；疨I鏈的楊氏模量分別為16.82 和17.18 GPa,相對(duì)誤差為2.1%。圖5(b)展示的是(5,5)的CNT(C500)和含有10個(gè)C 珠子的CNT 的周期性模型,從圖6(b)的應(yīng)力-應(yīng)變曲線得出全原子與粗?；疌NT 的楊氏模量分別為183.0和178.0 GPa,相對(duì)誤差為2.7%。兩組曲線表明粗?；?chǎng)參數(shù)能準(zhǔn)確還原PI與CNT 的拉伸性能。

圖5 全原子及粗?；Ｐ?(a)PI鏈;(b)CNTFig.5 All-atomic and coarse-grained model of (a) PI; (b) CNT

圖6 全原子與粗?；膽?yīng)力應(yīng)變曲線(a)PI;(b)CNTFig.6 Stress-strain curves between all-atomic and coarse-grained model of (a)PI chain; (b)CNT

3.2 驗(yàn)證CNT/PI復(fù)合材料的物理性質(zhì)

利用Accelrys Materials Studio 8.0 軟件的Amorphous Cell Packing模塊將含有4 個(gè)單體的PI分子鏈填充到盒子中,并將一根(5, 5)的CNT 置于盒子中間,建造一個(gè)大小為6.0 nm×6.0 nm×6.0 nm的CNT/PI復(fù)合材料的周期性模型,全原子模型共包含16 650個(gè)原子。通過(guò)全原子模型映射而成的粗粒化模型包含2 580個(gè)粗?；樽?模型的自由度相對(duì)于全原子的減少了約8倍。

首先對(duì)模型進(jìn)行充分的結(jié)構(gòu)幾何優(yōu)化,然后在NPT 系綜下對(duì)模型進(jìn)行動(dòng)力學(xué)弛豫,弛豫溫度設(shè)為298 K,壓強(qiáng)為101 kPa,全原子和粗?；Ｐ偷臅r(shí)間步長(zhǎng)分別設(shè)置為1 fs和10 fs,弛豫總時(shí)長(zhǎng)為1.5 ns。全原子模型弛豫耗時(shí)約18 h,粗?；Ｐ秃臅r(shí)約1.5 h,粗?；竽M效率可提升約10倍。圖7展示了全原子和粗?；Ｐ偷拿芏入S時(shí)間的變化曲線,平衡密度分別為1.185和1.184 g/cm3,相對(duì)誤差為0.8%。模擬結(jié)果表明粗?；?chǎng)參數(shù)可以準(zhǔn)確模擬出納米復(fù)合材料的密度。

圖7 CNT/PI復(fù)合材料的密度變化Fig.7 Variation of the density of CNT/PI composites

接著利用經(jīng)過(guò)充分弛豫的全原子和粗?；Ｐ瓦M(jìn)行CNT 的拔出模擬,以此驗(yàn)證CNT 與PI基體間的界面相互作用。如圖8(a)展示了CNT 的拔出模型,將CNT 沿著z軸方向以0.000 1 ?/fs的恒定速度拔出,輸出拔出能,并繪制出拔出能隨拔出距離變化的曲線,如圖8(b)所示。CNT 的拔出模擬是通過(guò)在CNT端部系一根彈簧,讓彈簧以固定速度移動(dòng),實(shí)現(xiàn)CNT的拔出,拔出力等于彈簧的拉力,拔出能等于拔出力與CNT 位移的乘積。

圖8 全原子與粗?；疌NT/PI復(fù)合材料的(a)CNT 拔出模型;(b)CNT 拔出能與拔出位移的曲線Fig.8 (a) Schematic illustration of CNT pull-out process; (b) pull-out energy response versus pull-out distance

圖8中粗?；Ｐ偷陌纬瞿艹尸F(xiàn)周期波動(dòng)的原因是 CNT的粗?；Ｐ椭械?C 珠子之間相距較遠(yuǎn),所以不能被近似看成連續(xù)一維材料,導(dǎo)致拔出能在 C珠子離開(kāi)基體過(guò)程中具有明顯的周期性變化。從圖8(b)中可以看出,粗?；Ｐ偷陌纬瞿芮€保持在全原子模型的拔出能曲線上下波動(dòng),該模擬結(jié)果證明粗?；?chǎng)參數(shù)可以還原出CNT和PI基體的界面相互作用。

最后將納米復(fù)合物沿著z軸方向拉伸,應(yīng)變率為2%,全原子和粗?；Ｐ偷膽?yīng)力-應(yīng)變曲線如圖9b所示。全原子與粗?；Ｐ偷臈钍夏Ａ糠謩e為8.13和8.16 GPa,相對(duì)誤差為3.6%,證明粗?；?chǎng)參數(shù)可以準(zhǔn)確還原出納米復(fù)合材料的楊氏模量。

圖9 全原子與粗?；疌NT/PI復(fù)合材料的(a)拉伸模型;(b)應(yīng)力應(yīng)變曲線Fig.9 CNT/PI composites between all-atomic and coarse-grained (a) the schematic illustration of stretching process; (b) stress-strain curves

3.3 粗粒化參數(shù)對(duì)復(fù)合材料物理性能的影響探討

探討粗?；瘏?shù)對(duì)模型的影響有助于深入理解粗?；黜?xiàng)參數(shù)對(duì)模型的物理性質(zhì)的影響效果和影響程度,便于后續(xù)對(duì)參數(shù)的修正,使粗粒化模擬更加準(zhǔn)確、適用性更強(qiáng)。該部分主要探討范德華參數(shù)中的r0和D0對(duì)CNT/PI復(fù)合材料的物理性能的影響,內(nèi)容包括密度、CNT 與PI基體的相互作用以及小應(yīng)變下的拉伸性能。小應(yīng)變下,影響拉伸性能的主要的分子間的相互作用力,鍵長(zhǎng)鍵角的參數(shù)主要與大應(yīng)變下的拉伸性能有關(guān)。

3.3.1r0對(duì)密度的影響 式(4)表明,D0與分子間相互作用有關(guān),r0與分子間的距離有關(guān)。增加A、B和C珠子的r0,Δr0依次取值0、0.3、0.6和1.0 ?,然后在NPT 系綜下對(duì)模型動(dòng)力學(xué)弛豫,繪制不同力場(chǎng)參數(shù)下模型的密度隨時(shí)間變化的曲線,如圖10所示。隨著Δr0從0增加到1.0 ?,模型的平衡密度從1.18減小到0.84 g/cm3,結(jié)果表明:r0增加,模型密度減小。

Δr0/? Density/(g·cm-3)0.0 1.18 0.3 1.07 0.6 0.95 1.0 0.84

圖10 不同r0 下模型密度隨時(shí)間的變化Fig.10 Variation of density versus r0

3.3.2D0對(duì)界面相互作用的影響 在CNT 的拔出模擬中,拔出力要克服CNT 與PI基體間的相互作用力做功?？梢酝ㄟ^(guò)不同D0下CNT的拔出能與拔出距離的變化曲線,比較界面相互作用強(qiáng)度。將C珠子與A、B珠子間的D0分別增加0、2、4和8 kcal·mol-1,然后進(jìn)行CNT的拔出模擬,并繪制不同D0下CNT的拔出能隨位移變化的曲線。如圖11所示,結(jié)果表明:CNT與PI間的D0越大,界面相互作用越大。

圖11 不同D0 下CNT 拔出能隨位移變化的曲線Fig.11 Curves of CNT pull-out energy versus displacement with different well depth parameters

3.3.3D0對(duì)楊氏模量的影響 探究D0與CNT/PI復(fù)合材料的楊氏模量的關(guān)系。改變D0,ΔD0分別取0.0、2.0、4.0及6.0 kcal·mol-1,模型進(jìn)行拉伸模擬,應(yīng)變率為3%,整理出不D0下模型的應(yīng)力應(yīng)變曲線,如圖12所示。由圖12和表5可以得出結(jié)論:在彈性拉伸階段,D0越大,分子間相互作用力越大,納米復(fù)合材料的楊氏模量越大。

表5 不同勢(shì)阱參數(shù)與對(duì)應(yīng)的復(fù)合物的楊氏模量Table 5 Different well depth parameters and corresponding Young's modulus of nanocomposite

圖12 不同勢(shì)阱參數(shù)下模型的應(yīng)力-應(yīng)變曲線Fig.12 Stress-strain curves of the model with different well depth parameters

4 結(jié) 論

1.根據(jù)相應(yīng)的映射方案,利用恒應(yīng)變方法可以完整地得出PI與CNT 的粗?；?chǎng)參數(shù)。

2.PI分子鏈以及CNT 的拉伸模擬中,粗?；Ｐ团c全原子模型的楊氏模量的相對(duì)誤差低于5%,證明粗?；?chǎng)參數(shù)能準(zhǔn)確還原單鏈的拉伸性質(zhì)。

3.粗?；Ｐ涂梢猿杀稖p少原子數(shù),并使得模擬效率提高約10倍。

4.在CNT 拔出模擬中,全原子與粗?；Ｐ偷腃NT 拔出能與拔出距離的變化曲線基本吻合,證明粗粒化力場(chǎng)參數(shù)可以準(zhǔn)確還原CNT 與PI基體間的界面相互作用。

5.對(duì)模型進(jìn)行拉伸模擬,得出全原子與粗?；Ｐ偷臈钍夏Ａ康南鄬?duì)誤差低于5%,證明粗粒化力場(chǎng)參數(shù)可以還原復(fù)合物的彈性拉伸性質(zhì)。

6.在關(guān)于粗?；?chǎng)參數(shù)與材料物理性能的關(guān)系的研究中得出結(jié)論:分子間的平衡距離r0增加,模型體積增大,密度減小;勢(shì)阱參數(shù)D0增加,分子間作用力增加,CNT 與PI基體的界面相互作用增大,楊氏模量增加。該結(jié)論可為未來(lái)利用粗?；椒ㄌ剿鰿NT表面缺陷、不同改性等因素對(duì)復(fù)合基體性能影響方面提供理論支撐。