楊海霞，王梓霖，孫付宇，高如月，付一政

(1.山西工程科技職業(yè)大學(xué)信息工程學(xué)院,山西晉中 030619； 2.中北大學(xué)材料科學(xué)與工程學(xué)院,太原 030051；3.天津金發(fā)新材料有限公司,天津 300308)

環(huán)氧樹脂是一種低成本,易獲得的熱固性樹脂[1]。具有低密度、高強(qiáng)度、尺寸穩(wěn)定性好和耐化學(xué)腐蝕性等特點[2-3],被廣泛應(yīng)用于粘合劑、汽車、電子等行業(yè)[4～6]。但由于其交聯(lián)密度高,脆性大,耐沖擊性和摩擦磨損性能差,限制了其在各領(lǐng)域的應(yīng)用。因此,提高環(huán)氧樹脂力學(xué)性能具有十分重要的實用價值[7～9]。由于實驗方法的局限性,從原子水平對環(huán)氧樹脂進(jìn)行實驗研究非常困難。分子模擬方法以其獨特的優(yōu)勢已成為研究高分子材料結(jié)構(gòu)和性能的有效工具,分子動力學(xué)(MD)模擬技術(shù)是近年來飛速發(fā)展的一種分子模擬方法,其可以從分子或原子層面揭示分子間的相互作用,從微觀角度研究材料受力過程中分子、原子的結(jié)構(gòu)變化情況,為材料的分子設(shè)計提供理論依據(jù)。已有研究者借助MD模擬方法研究了熱塑性樹脂的力學(xué)性能[10-16],但關(guān)于熱固性樹脂的研究卻很少[17-19],這主要是由于具有三維交聯(lián)網(wǎng)絡(luò)結(jié)構(gòu)的熱固性樹脂通常是由預(yù)聚物與固化劑經(jīng)交聯(lián)反應(yīng)而成,交聯(lián)結(jié)構(gòu)的構(gòu)建比較復(fù)雜,不能通過簡單的分子建模方法得到。課題組前期已經(jīng)在Materials Studio (MS)軟件中通過Perl 語言編寫可以生成交聯(lián)環(huán)氧樹脂分子模型的腳本,并通過MD模擬預(yù)測了交聯(lián)環(huán)氧樹脂的玻璃化轉(zhuǎn)變溫度和靜態(tài)力學(xué)性能等性質(zhì)[20]。筆者在以前工作的基礎(chǔ)上,進(jìn)一步采用非平衡態(tài)分子動力學(xué)模擬方法模擬了不同交聯(lián)度環(huán)氧樹脂分子模型的拉伸過程,得到了模型的應(yīng)力-應(yīng)變曲線,并從微觀角度分析了交聯(lián)環(huán)氧樹脂在拉伸過程中微觀結(jié)構(gòu)和能量等的變化情況。

1 模型構(gòu)建與模擬計算

1.1 初始模型構(gòu)建

首先,在MS 軟件中繪制了預(yù)聚物4,4-二氨基二苯甲烷四縮水甘油醚(TGDDM)和固化劑4,4-二氨基二苯胺(DDM)的分子模型,其結(jié)構(gòu)式如圖1 所示,然后構(gòu)建了10 個尺寸為5 nm×5 nm×5 nm 的周期性盒子,將TGDDM和DDM按照分子個數(shù)為1∶1的比例隨機(jī)填充到周期性盒子中,最后對10個模型采用Smart Minimizer方法進(jìn)行結(jié)構(gòu)優(yōu)化,選擇能量最低的分子模型進(jìn)行500 ps 的恒溫恒容(NVT)和500 ps 的恒溫恒壓(NPT)MD 模擬。最終得到未交聯(lián)的預(yù)聚物與固化劑混合物模型,如圖2所示。

圖1 TGDDM和DDM的分子結(jié)構(gòu)式

圖2 未交聯(lián)的TGDDM與DDM混合分子模型

1.2 交聯(lián)模型構(gòu)建

TGDDM與DDM的交聯(lián)反應(yīng)機(jī)理如圖3所示,反應(yīng)過程主要分為兩步：第一步,環(huán)氧樹脂上的環(huán)氧基團(tuán)開環(huán)后與固化劑上的伯胺基團(tuán)發(fā)生加成反應(yīng)生成仲胺基團(tuán),開環(huán)后的氧原子與氫原子形成羥基。第二步,另一個環(huán)氧基團(tuán)開環(huán)與仲胺基團(tuán)發(fā)生加成反應(yīng)生成叔胺基團(tuán)。根據(jù)交聯(lián)反應(yīng)機(jī)理,通過MS軟件中的perl腳本建立了TGDDM與DDM的交聯(lián)模型,具體交聯(lián)細(xì)節(jié)見文獻(xiàn)[20]。

圖3 反應(yīng)機(jī)理示意圖

1.3 拉伸過程模擬

拉伸過程模擬通過美國Sandia 國家實驗室開發(fā)的軟件(LAMMPS）進(jìn)行。首先將得到的交聯(lián)環(huán)氧樹脂模型轉(zhuǎn)化為LAMMPS 軟件的輸入文件。然后對交聯(lián)模型進(jìn)行降溫處理,具體策略如下：①2 000 步幾何優(yōu)化；②溫度 600 K,200 ps 的 NVT 系綜和 200 ps 的 NPT 系綜的 MD 模擬；③以 30 K/200 ps 的降溫速率,進(jìn)行溫度從600 K 降至300 K 的NPT 系綜的 MD 模擬,模擬時間為 2 ns；④300 K,500 ps的NPT系綜的MD模擬。

最后在NPT系綜下對模型進(jìn)行單軸拉伸模擬,在z方向以恒定的應(yīng)變速率對模型施加應(yīng)變,并保證體系在橫向(x和y方向)的壓強(qiáng)穩(wěn)定在標(biāo)準(zhǔn)大氣壓。在整個模擬過程中使用COMPASS 力場,控溫和控壓方法分別采用Andersen 和Berendsen 方法,壓強(qiáng)設(shè)定為0.1 MPa,時間步長為1 fs。Atom-based和Ewald方法計算范德華(vdW)和靜電作用。

2 結(jié)果與討論

2.1 不同交聯(lián)度環(huán)氧樹脂的力學(xué)性能

筆者通過MS軟件中Perl腳本從交聯(lián)反應(yīng)截斷半徑0.55 nm 開始交聯(lián),到1.25 nm 結(jié)束,建立一系列不同交聯(lián)度的交聯(lián)環(huán)氧樹脂分子模型,交聯(lián)過程中交聯(lián)密度、模型的密度和交聯(lián)反應(yīng)截斷半徑的關(guān)系列于表1,最終得到的最大交聯(lián)度為90.21%。隨后將表1 中的模型導(dǎo)入到LAMMPS 軟件中,在300 K下采用deform命令對交聯(lián)模型進(jìn)行拉伸模擬。

表1 交聯(lián)度、密度與截斷半徑的關(guān)系

圖4 為交聯(lián)度為90.21%模型在應(yīng)變速率為5×108s-1時的應(yīng)力-應(yīng)變曲線。由圖4 可見,在整個拉伸過程中,垂直與拉伸方向(圖 2 中的x和y方向)上的應(yīng)力σx和σy基本為0。當(dāng)應(yīng)變＜0.01 時,沿著拉伸方向(圖2 中的z方向)的應(yīng)力σz隨應(yīng)變的增加而線性增加,為彈性區(qū),超過0.01 后材料進(jìn)入屈服階段,屈服強(qiáng)度可達(dá)110 MPa,在屈服應(yīng)變之后,應(yīng)力隨應(yīng)變的增加而下降,這與真實的實驗現(xiàn)象一致[21]。

圖4 交聯(lián)度為90.21%模型在不同方向的應(yīng)力-應(yīng)變圖

但當(dāng)應(yīng)變超過0.012 后,應(yīng)力值又持續(xù)增長,這主要是由于交聯(lián)環(huán)氧樹脂的斷裂伸長率較低,實驗中當(dāng)應(yīng)變超過一定值后,交聯(lián)分子鏈會發(fā)生斷裂,但由于受到本文分子模擬中使用力場的限制,無法模擬分子鏈的斷裂過程,所以應(yīng)力值會持續(xù)增長。另外,由于實驗測量時應(yīng)變率一般為1×10-1s-1,遠(yuǎn)小于模擬時的應(yīng)變速率(5×108s-1),由于“應(yīng)變硬化”導(dǎo)致模擬值較實驗值偏高,這與文獻(xiàn)[22]等的研究結(jié)果一致。

圖5為不同交聯(lián)反應(yīng)截斷半徑下得到的交聯(lián)模型的應(yīng)力-應(yīng)變圖,由圖5和表1可見,交聯(lián)反應(yīng)截斷半徑越大,得到的模型交聯(lián)度越大。各個模型的應(yīng)力-應(yīng)變圖變化趨勢比較相似,拉伸強(qiáng)度隨交聯(lián)度的增大而增大。

圖5 不同交聯(lián)反應(yīng)截斷半徑下環(huán)氧樹脂模型的應(yīng)力-應(yīng)變圖

2.2 拉伸過程中結(jié)構(gòu)變化情況

為了進(jìn)一步觀察拉伸過程中交聯(lián)環(huán)氧樹脂模型微觀結(jié)構(gòu)的變化,截取了不同應(yīng)變時的結(jié)構(gòu)快照如圖6所示。從圖6可知,由于交聯(lián)結(jié)構(gòu)不完整,模型會出現(xiàn)空洞,交聯(lián)度越小,空洞的數(shù)目越多,形狀越大；隨著應(yīng)變不斷增加,模型沿z方向的長度不斷增大,空洞也隨之變大、變長。但由于模擬的應(yīng)變較小(＜0.06)和三維網(wǎng)狀交聯(lián)結(jié)構(gòu)的限制,模型的結(jié)構(gòu)和空洞的形狀變化并不十分明顯。

圖6 不同交聯(lián)度環(huán)氧樹脂模型在不同應(yīng)變時的快照圖

2.3 拉伸過程中的質(zhì)量密度變化情況

為了進(jìn)一步探索拉伸過程中交聯(lián)環(huán)氧樹脂模型微觀結(jié)構(gòu)變化情況,分析了交聯(lián)度為90.21%環(huán)氧樹脂模型在拉伸不同時刻即初始階段(應(yīng)變?yōu)?)、彈性階段(應(yīng)變?yōu)?0.005)、屈服階段(應(yīng)變?yōu)?0.010)、屈服過后應(yīng)力下降階段(應(yīng)變?yōu)?.015)和應(yīng)力上升階段(應(yīng)變分別為0.020,0.025,0.030)時,交聯(lián)環(huán)氧樹脂分子模型沿z方向的質(zhì)量密度分布曲線,如圖7所示。從圖7可以看出,由于交聯(lián)結(jié)構(gòu)的限制,在整個拉伸過程中,模型沿z方向的質(zhì)量密度不會發(fā)生較大的變化,與2.2分析結(jié)果一致。

圖7 不同應(yīng)變下交聯(lián)環(huán)氧樹脂模型z方向的質(zhì)量密度分布圖

2.4 交聯(lián)環(huán)氧樹脂的泊松比

通過計算橫向應(yīng)變與縱向應(yīng)變的比值可以計算材料的泊松比。圖8 為交聯(lián)度為90.21%的環(huán)氧樹脂模型在拉伸過程中橫向應(yīng)變與縱向應(yīng)變的變化情況,在這里橫向應(yīng)變的值為x和y方向壓縮應(yīng)變的平均值,縱向應(yīng)變?yōu)閦方向的拉伸應(yīng)變。通過線性回歸求出泊松比為0.30,實驗中典型熱固性樹脂在玻璃態(tài)下的泊松比為0.30～0.46,模擬值在實驗范圍內(nèi)[23]。

圖8 橫向應(yīng)變與縱向應(yīng)變的關(guān)系

2.5 拉伸過程中各能量項分析

進(jìn)一步分析了拉伸過程中各能量項的變化情況如圖9所示,為了便于比較,所有能量值均為相對于初始模型能量的差值(Ei-E0),結(jié)合圖4 的應(yīng)力-應(yīng)變圖可以看出,在彈性區(qū)域鍵能、角能和二面角能隨應(yīng)變的增加而增加,在屈服階段達(dá)到最大值,這說明在彈性區(qū),模型的結(jié)構(gòu)在拉力的作用下發(fā)生了小幅度的變形和扭轉(zhuǎn),達(dá)到屈服階段后,結(jié)構(gòu)變化趨于均勻一致；超過屈服應(yīng)力后,由于交聯(lián)結(jié)構(gòu)的限制,使得交聯(lián)模型的鍵,角和二面角不會發(fā)生明顯變化,各能量項基本處于波動平衡狀態(tài)。而vdW能的變化趨勢與應(yīng)力-應(yīng)變圖變化趨勢大體一致,說明拉伸過程中,體系應(yīng)變主要引起了vdW相互作用能的改變。

圖9 應(yīng)變過程中各能量項的變化情況

3 結(jié)論

通過分子動力學(xué)模擬方法,利用Perl 腳本通過Materials Studio分子模擬軟件構(gòu)建了不同交聯(lián)度的環(huán)氧樹脂模型,然后對交聯(lián)模型進(jìn)行單軸拉伸,研究了交聯(lián)模型的拉伸性能,從拉伸過程中結(jié)構(gòu)和能量的變化對拉伸機(jī)理進(jìn)行了分析,主要結(jié)論如下：

(1)不同交聯(lián)度的環(huán)氧樹脂的應(yīng)力-應(yīng)變趨勢比較相似,拉伸強(qiáng)度隨交聯(lián)度的增大而增大,當(dāng)交聯(lián)度為90.21%時,拉伸強(qiáng)度為110 MPa,泊松比為0.30,模擬值較實驗值偏高,主要是由于模擬的應(yīng)變速率較實驗值高。

(2)由于三維交聯(lián)網(wǎng)狀結(jié)構(gòu)的限制,拉伸過程中模型沿z方向的質(zhì)量密度分布和微觀結(jié)構(gòu)變化并不十分明顯。

(3)在拉伸過程中,在彈性區(qū)域鍵能、角能和二面角能隨應(yīng)變的增加而增加,超過屈服應(yīng)力后,由于交聯(lián)結(jié)構(gòu)的限制,使得交聯(lián)模型的鍵,角和二面角不會發(fā)生明顯變化。體系應(yīng)變主要引起了vdW相互作用的改變。