基于恒應變方法的碳納米管/聚酰亞胺復合材料的粗?；鲩_發

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

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

1 前 言

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

雖然MD 已經廣泛應用于納米復合材料領域,但MD 模擬計算的力場復雜,致使工作量巨大,效率較低,使其難以進行更大時間和空間尺度的模擬。為了克服該問題,粗粒化(CG)方法得到了發展[13-15]。CG的原理是將一組原子看作一個整體,忽略內部信息,然后映射成對應的粗?；樽?從而保持原子模型的分子細節的同時,成倍縮減全原子體系的自由度和粒子數,增大模擬的時間和空間尺度。目前,已有文獻進行了聚合物材料的粗粒化模型及其力場參數的初步開發工作[16-18]。Marrink 課題組用粗粒化方法得到了MARTINI力場[19],主要用于研究生物大分子的動態行為,比如磷脂分子、蛋白質及DNA 等。蔡莊立等[20]根據特定的映射方案建立了PE 的粗?；Ｐ?并以全原子模型的鍵長、鍵角分布曲線來獲取PE 的粗?；鰠?。文中驗證了粗?；瘏的MPE 模型的靜態結構分布,且粗?；Ｐ偷膶嵯禂档哪M值與實驗值符合較好。Langeloth等[21]利用玻爾茲曼反轉迭代法得到了環氧樹脂的粗?；Ｐ鸵约傲鰠?并利用此參數在更大的尺度下研究了環氧樹脂的玻璃化轉變溫度與交聯程度的變化關系。

然而由于目前的粗?；芯咳鄙賹Υ至；鰠禍蚀_性的驗證,粗?；鰠蹬c模型物理性能之間的關系亦鮮有討論。同時,已有的玻爾茲曼反轉迭代法[22-23]等方法操作復雜,適用范圍有限。本研究使用基于應變能守恒定律的恒應變方法獲得聚酰亞胺(PI)和CNT 的力場參數,然后通過不同模擬實驗來驗證粗?；瘏档臏蚀_性,最后對粗?；瘏蹬c物理性能的關系將進行詳細討論。

2 建模與粗?；瘏档墨@取

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

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

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

2.1 粗粒化鍵長參數的獲取

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

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

接下來使用恒應變方法獲取粗粒化力場參數。首先改變A-A 珠子對應的原子團之間的鍵長,通過分子模擬得到體系的總勢能U,根據應變能守恒,粗?；w系的總勢能也等于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珠子對應的原子團之間的鍵長,得出U與Δd的變化曲線。如圖2(b)所示,根據式(5)、(6),依次得到d0=4.80 ?,kd=340.78 kcal·mol-1·?-2。改變C-C 珠子對應的原子團之間的鍵長,得出U與Δd的變化曲線,如圖2(c)所示。根據式(5)、(6)得出:d0=6.12 ?,kd=4 542.0 kcal·mol-1·?-2。

2.2 粗?；I角參數的獲取

改變A-A-B珠子對應的原子團的角度,通過分子模擬得到體系的總勢能U,根據應變能守恒,粗?；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 珠子對應的原子團的角度,得到U與Δθ的變化曲線,如圖3(b)所示。根據式(7)、(8)得出θ0=179.4°,kθ=670.21 kcal·mol-1·rad-2。最后,改變C-C-C珠子對應的原子團的角度,得出U與Δθ的變化曲線,如圖3(c)所示。根據式(7)、(8)得出θ0=180.0°,kθ=66 880.0 kcal·mol-1·rad-2。

2.3 粗?；兜氯A參數的獲取

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

表1 A、B、C珠子的勢阱參數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

利用得到的粗粒化參數對CNT/PI復合材料的粗?；Ｐ瓦M行充分弛豫,然后將弛豫后的粗粒化模型的結構分布曲線,包括鍵長、鍵角以及徑向分布曲線,與全原子模型的進行對比,并根據分布曲線的誤差修正粗粒化力場參數,最后得出調整后的參數,如表2所示。

表2 修正后的粗?；瘏礣able 2 Modified coarse-grained parameters

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

3 粗粒化參數的驗證和探討

3.1 驗證PI分子鏈和CNT的拉伸性質

如圖5(a)所示,一條PI分子鏈含有四個單體,一條PI珠子鏈含有8個A 珠子,4個B珠子,分子鏈的長度為66.17 ?,對兩條鏈進行拉伸,應變率為8%,得到圖6(a)中的應力-應變曲線。對曲線進行擬合得到全原子與粗粒化PI鏈的楊氏模量分別為16.82 和17.18 GPa,相對誤差為2.1%。圖5(b)展示的是(5,5)的CNT(C500)和含有10個C 珠子的CNT 的周期性模型,從圖6(b)的應力-應變曲線得出全原子與粗粒化CNT 的楊氏模量分別為183.0和178.0 GPa,相對誤差為2.7%。兩組曲線表明粗?；鰠的軠蚀_還原PI與CNT 的拉伸性能。

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

圖6 全原子與粗?；膽兦€(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復合材料的物理性質

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

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

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

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

圖8 全原子與粗?；疌NT/PI復合材料的(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周期波動的原因是 CNT的粗?；Ｐ椭械?C 珠子之間相距較遠,所以不能被近似看成連續一維材料,導致拔出能在 C珠子離開基體過程中具有明顯的周期性變化。從圖8(b)中可以看出,粗粒化模型的拔出能曲線保持在全原子模型的拔出能曲線上下波動,該模擬結果證明粗粒化力場參數可以還原出CNT和PI基體的界面相互作用。

最后將納米復合物沿著z軸方向拉伸,應變率為2%,全原子和粗粒化模型的應力-應變曲線如圖9b所示。全原子與粗粒化模型的楊氏模量分別為8.13和8.16 GPa,相對誤差為3.6%,證明粗粒化力場參數可以準確還原出納米復合材料的楊氏模量。

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

3.3 粗粒化參數對復合材料物理性能的影響探討

探討粗?；瘏祵δＰ偷挠绊懹兄谏钊肜斫獯至；黜梾祵δＰ偷奈锢硇再|的影響效果和影響程度,便于后續對參數的修正,使粗?；M更加準確、適用性更強。該部分主要探討范德華參數中的r0和D0對CNT/PI復合材料的物理性能的影響,內容包括密度、CNT 與PI基體的相互作用以及小應變下的拉伸性能。小應變下,影響拉伸性能的主要的分子間的相互作用力,鍵長鍵角的參數主要與大應變下的拉伸性能有關。

3.3.1r0對密度的影響 式(4)表明,D0與分子間相互作用有關,r0與分子間的距離有關。增加A、B和C珠子的r0,Δr0依次取值0、0.3、0.6和1.0 ?,然后在NPT 系綜下對模型動力學弛豫,繪制不同力場參數下模型的密度隨時間變化的曲線,如圖10所示。隨著Δr0從0增加到1.0 ?,模型的平衡密度從1.18減小到0.84 g/cm3,結果表明: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所示,結果表明: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復合材料的楊氏模量的關系。改變D0,ΔD0分別取0.0、2.0、4.0及6.0 kcal·mol-1,模型進行拉伸模擬,應變率為3%,整理出不D0下模型的應力應變曲線,如圖12所示。由圖12和表5可以得出結論:在彈性拉伸階段,D0越大,分子間相互作用力越大,納米復合材料的楊氏模量越大。

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

圖12 不同勢阱參數下模型的應力-應變曲線Fig.12 Stress-strain curves of the model with different well depth parameters

4 結 論

1.根據相應的映射方案,利用恒應變方法可以完整地得出PI與CNT 的粗粒化力場參數。

2.PI分子鏈以及CNT 的拉伸模擬中,粗?；Ｐ团c全原子模型的楊氏模量的相對誤差低于5%,證明粗?；鰠的軠蚀_還原單鏈的拉伸性質。

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

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

5.對模型進行拉伸模擬,得出全原子與粗?；Ｐ偷臈钍夏Ａ康南鄬φ`差低于5%,證明粗?；鰠悼梢赃€原復合物的彈性拉伸性質。

6.在關于粗?；鰠蹬c材料物理性能的關系的研究中得出結論:分子間的平衡距離r0增加,模型體積增大,密度減小;勢阱參數D0增加,分子間作用力增加,CNT 與PI基體的界面相互作用增大,楊氏模量增加。該結論可為未來利用粗粒化方法探索CNT表面缺陷、不同改性等因素對復合基體性能影響方面提供理論支撐。