基于恒應(yīng)變方法的碳納米管/聚酰亞胺復(fù)合材料的粗粒化力場開發(fā)

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

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

1 前 言

自1991年碳納米管(CNT)被日本科學(xué)家發(fā)現(xiàn)以來,其優(yōu)異的電學(xué)、力學(xué)和熱學(xué)性能引起了廣泛關(guān)注[1]。超高的楊氏模量和抗拉強度使其成為高性能聚合物基復(fù)合材料的首選增強填料[2-5],在CNT 增強聚合物復(fù)合材料的力學(xué)性能方面,學(xué)者們進行了廣泛的實驗研究。汪華鋒等[6]利用催化裂解制備的CNT/環(huán)氧樹脂復(fù)合材料,力學(xué)性能測試表明:CNT 可以使環(huán)氧樹脂基體材料的拉伸及壓縮強度分別提高284%和122%,同時也可增加環(huán)氧樹脂材料的韌性。蓬勃等[7]制備了質(zhì)量分數(shù)不同的CNT/聚乙烯(PE)復(fù)合材料顯示:當(dāng)填充0.5 wt%的CNT 時,斷裂伸長率達到最大值429.2%,再進一步增加CNT 填充量,其斷裂伸長率呈單調(diào)下降趨勢。由于受到時間和空間尺度的限制,在納米尺度下材料間的交互過程和作用機制(如CNT 和基體之間的界面相互作用)等問題的探索在實驗中仍然難以有效開展。分子動力學(xué)(MD)模擬作為近年來新興的研究手段,可以在分子層面模擬納米復(fù)合材料基體內(nèi)的交互過程和作用機制[8-9],成為彌補實驗研究局限性的有力補充。Frankland等[10]采用MD方法研究了(10, 10)單壁CNT(SWCNT)增強PE 基復(fù)合材料的應(yīng)力-應(yīng)變曲線。結(jié)果表明,長徑比為10的長SWCNT 纖維增強的納米復(fù)合材料的剛度顯著增加,對于長徑比為4的CNT 增強體,則沒有觀察到有顯著的增強。Li等[11]采用MD模擬方法研究了具有預(yù)定單邊裂紋的CNT/環(huán)氧樹脂復(fù)合材料的斷裂行為。結(jié)果表明,引入CNT后,環(huán)氧復(fù)合材料的拉伸強度和斷裂伸長率分別提高了24.8%和34.3%。Mokashi等[12]對(10, 10)CNT 增強非晶態(tài)PE聚合物基體的楊氏模量和拉伸強度進行了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 模擬計算的力場復(fù)雜,致使工作量巨大,效率較低,使其難以進行更大時間和空間尺度的模擬。為了克服該問題,粗?；?CG)方法得到了發(fā)展[13-15]。CG的原理是將一組原子看作一個整體,忽略內(nèi)部信息,然后映射成對應(yīng)的粗粒化珠子,從而保持原子模型的分子細節(jié)的同時,成倍縮減全原子體系的自由度和粒子數(shù),增大模擬的時間和空間尺度。目前,已有文獻進行了聚合物材料的粗?；Ｐ图捌淞鰠?shù)的初步開發(fā)工作[16-18]。Marrink 課題組用粗粒化方法得到了MARTINI力場[19],主要用于研究生物大分子的動態(tài)行為,比如磷脂分子、蛋白質(zhì)及DNA 等。蔡莊立等[20]根據(jù)特定的映射方案建立了PE 的粗?；Ｐ?并以全原子模型的鍵長、鍵角分布曲線來獲取PE 的粗粒化力場參數(shù)。文中驗證了粗粒化參數(shù)模擬PE 模型的靜態(tài)結(jié)構(gòu)分布,且粗?；Ｐ偷膶?dǎo)熱系數(shù)的模擬值與實驗值符合較好。Langeloth等[21]利用玻爾茲曼反轉(zhuǎn)迭代法得到了環(huán)氧樹脂的粗粒化模型以及力場參數(shù),并利用此參數(shù)在更大的尺度下研究了環(huán)氧樹脂的玻璃化轉(zhuǎn)變溫度與交聯(lián)程度的變化關(guān)系。

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

2 建模與粗?；瘏?shù)的獲取

粗粒化力場參數(shù)主要由成鍵勢函數(shù)和非鍵勢函數(shù)兩部分構(gòu)成。體系總勢能Etotal可以寫成鍵長能Eb、鍵角能Ea、范德華能Evdw和系統(tǒng)恒定自由能U0之和,表達式形式如下:

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

式中:kd和d0是鍵長的剛度系數(shù)和平衡鍵長;kθ和θ0是鍵角的剛度系數(shù)和平衡鍵角;D0和r0分別是Lennard—Jones(LJ)勢函數(shù)的勢阱參數(shù)和平衡距離。全原子模擬使用COMPASS力場來描述分子內(nèi)的相互作用力,非鍵相互作用的截斷半徑取1.2 nm,能量最小化方法使用共軛梯度法。

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

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

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

接下來使用恒應(yīng)變方法獲取粗粒化力場參數(shù)。首先改變A-A 珠子對應(yīng)的原子團之間的鍵長,通過分子模擬得到體系的總勢能U,根據(jù)應(yīng)變能守恒,粗粒化體系的總勢能也等于U:

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

圖2 原子團勢能與鍵長增量的變化 (a)A-A 原子團;(b)A-B原子團;(c)C-C原子團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珠子對應(yīng)的原子團之間的鍵長,得出U與Δd的變化曲線。如圖2(b)所示,根據(jù)式(5)、(6),依次得到d0=4.80 ?,kd=340.78 kcal·mol-1·?-2。改變C-C 珠子對應(yīng)的原子團之間的鍵長,得出U與Δd的變化曲線,如圖2(c)所示。根據(jù)式(5)、(6)得出:d0=6.12 ?,kd=4 542.0 kcal·mol-1·?-2。

2.2 粗?；I角參數(shù)的獲取

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

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

圖3 原子團勢能與鍵角增量的變化 (a)A-A-B原子團;(b)A-B-A 原子團;(c)C-C-C原子團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 珠子對應(yīng)的原子團的角度,得到U與Δθ的變化曲線,如圖3(b)所示。根據(jù)式(7)、(8)得出θ0=179.4°,kθ=670.21 kcal·mol-1·rad-2。最后,改變C-C-C珠子對應(yīng)的原子團的角度,得出U與Δθ的變化曲線,如圖3(c)所示。根據(jù)式(7)、(8)得出θ0=180.0°,kθ=66 880.0 kcal·mol-1·rad-2。

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

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

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

圖4 A-A 原子團的勢能與距離的變化Fig.4 Variation of energy versus distance between A-A group

利用得到的粗?；瘏?shù)對CNT/PI復(fù)合材料的粗粒化模型進行充分弛豫,然后將弛豫后的粗?；Ｐ偷慕Y(jié)構(gòu)分布曲線,包括鍵長、鍵角以及徑向分布曲線,與全原子模型的進行對比,并根據(jù)分布曲線的誤差修正粗?；鰠?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ù)的驗證和探討

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

如圖5(a)所示,一條PI分子鏈含有四個單體,一條PI珠子鏈含有8個A 珠子,4個B珠子,分子鏈的長度為66.17 ?,對兩條鏈進行拉伸,應(yīng)變率為8%,得到圖6(a)中的應(yīng)力-應(yīng)變曲線。對曲線進行擬合得到全原子與粗?；疨I鏈的楊氏模量分別為16.82 和17.18 GPa,相對誤差為2.1%。圖5(b)展示的是(5,5)的CNT(C500)和含有10個C 珠子的CNT 的周期性模型,從圖6(b)的應(yīng)力-應(yīng)變曲線得出全原子與粗?；疌NT 的楊氏模量分別為183.0和178.0 GPa,相對誤差為2.7%。兩組曲線表明粗?；鰠?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 驗證CNT/PI復(fù)合材料的物理性質(zhì)

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

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

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

接著利用經(jīng)過充分弛豫的全原子和粗?；Ｐ瓦M行CNT 的拔出模擬,以此驗證CNT 與PI基體間的界面相互作用。如圖8(a)展示了CNT 的拔出模型,將CNT 沿著z軸方向以0.000 1 ?/fs的恒定速度拔出,輸出拔出能,并繪制出拔出能隨拔出距離變化的曲線,如圖8(b)所示。CNT 的拔出模擬是通過在CNT端部系一根彈簧,讓彈簧以固定速度移動,實現(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)周期波動的原因是 CNT的粗?；Ｐ椭械?C 珠子之間相距較遠,所以不能被近似看成連續(xù)一維材料,導(dǎo)致拔出能在 C珠子離開基體過程中具有明顯的周期性變化。從圖8(b)中可以看出,粗?；Ｐ偷陌纬瞿芮€保持在全原子模型的拔出能曲線上下波動,該模擬結(jié)果證明粗?；鰠?shù)可以還原出CNT和PI基體的界面相互作用。

最后將納米復(fù)合物沿著z軸方向拉伸,應(yīng)變率為2%,全原子和粗粒化模型的應(yīng)力-應(yīng)變曲線如圖9b所示。全原子與粗?；Ｐ偷臈钍夏Ａ糠謩e為8.13和8.16 GPa,相對誤差為3.6%,證明粗?；鰠?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ù)對復(fù)合材料物理性能的影響探討

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

3.3.1r0對密度的影響 式(4)表明,D0與分子間相互作用有關(guān),r0與分子間的距離有關(guān)。增加A、B和C珠子的r0,Δr0依次取值0、0.3、0.6和1.0 ?,然后在NPT 系綜下對模型動力學(xué)弛豫,繪制不同力場參數(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 下模型密度隨時間的變化Fig.10 Variation of density versus r0

3.3.2D0對界面相互作用的影響 在CNT 的拔出模擬中,拔出力要克服CNT 與PI基體間的相互作用力做功?？梢酝ㄟ^不同D0下CNT的拔出能與拔出距離的變化曲線,比較界面相互作用強度。將C珠子與A、B珠子間的D0分別增加0、2、4和8 kcal·mol-1,然后進行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對楊氏模量的影響 探究D0與CNT/PI復(fù)合材料的楊氏模量的關(guān)系。改變D0,ΔD0分別取0.0、2.0、4.0及6.0 kcal·mol-1,模型進行拉伸模擬,應(yīng)變率為3%,整理出不D0下模型的應(yīng)力應(yīng)變曲線,如圖12所示。由圖12和表5可以得出結(jié)論:在彈性拉伸階段,D0越大,分子間相互作用力越大,納米復(fù)合材料的楊氏模量越大。

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

圖12 不同勢阱參數(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 的粗?；鰠?shù)。

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

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

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

5.對模型進行拉伸模擬,得出全原子與粗?；Ｐ偷臈钍夏Ａ康南鄬φ`差低于5%,證明粗?；鰠?shù)可以還原復(fù)合物的彈性拉伸性質(zhì)。

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