交联环氧树脂拉伸性能的分子动力学模拟

杨海霞，王梓霖，孙付宇，高如月，付一政

(1.山西工程科技职业大学信息工程学院,山西晋中 030619； 2.中北大学材料科学与工程学院,太原 030051；3.天津金发新材料有限公司,天津 300308)

环氧树脂是一种低成本,易获得的热固性树脂[1]。具有低密度、高强度、尺寸稳定性好和耐化学腐蚀性等特点[2-3],被广泛应用于粘合剂、汽车、电子等行业[4～6]。但由于其交联密度高,脆性大,耐冲击性和摩擦磨损性能差,限制了其在各领域的应用。因此,提高环氧树脂力学性能具有十分重要的实用价值[7～9]。由于实验方法的局限性,从原子水平对环氧树脂进行实验研究非常困难。分子模拟方法以其独特的优势已成为研究高分子材料结构和性能的有效工具,分子动力学(MD)模拟技术是近年来飞速发展的一种分子模拟方法,其可以从分子或原子层面揭示分子间的相互作用,从微观角度研究材料受力过程中分子、原子的结构变化情况,为材料的分子设计提供理论依据。已有研究者借助MD模拟方法研究了热塑性树脂的力学性能[10-16],但关于热固性树脂的研究却很少[17-19],这主要是由于具有三维交联网络结构的热固性树脂通常是由预聚物与固化剂经交联反应而成,交联结构的构建比较复杂,不能通过简单的分子建模方法得到。课题组前期已经在Materials Studio (MS)软件中通过Perl 语言编写可以生成交联环氧树脂分子模型的脚本,并通过MD模拟预测了交联环氧树脂的玻璃化转变温度和静态力学性能等性质[20]。笔者在以前工作的基础上,进一步采用非平衡态分子动力学模拟方法模拟了不同交联度环氧树脂分子模型的拉伸过程,得到了模型的应力-应变曲线,并从微观角度分析了交联环氧树脂在拉伸过程中微观结构和能量等的变化情况。

1 模型构建与模拟计算

1.1 初始模型构建

首先,在MS 软件中绘制了预聚物4,4-二氨基二苯甲烷四缩水甘油醚(TGDDM)和固化剂4,4-二氨基二苯胺(DDM)的分子模型,其结构式如图1 所示,然后构建了10 个尺寸为5 nm×5 nm×5 nm 的周期性盒子,将TGDDM和DDM按照分子个数为1∶1的比例随机填充到周期性盒子中,最后对10个模型采用Smart Minimizer方法进行结构优化,选择能量最低的分子模型进行500 ps 的恒温恒容(NVT)和500 ps 的恒温恒压(NPT)MD 模拟。最终得到未交联的预聚物与固化剂混合物模型,如图2所示。

图1 TGDDM和DDM的分子结构式

图2 未交联的TGDDM与DDM混合分子模型

1.2 交联模型构建

TGDDM与DDM的交联反应机理如图3所示,反应过程主要分为两步：第一步,环氧树脂上的环氧基团开环后与固化剂上的伯胺基团发生加成反应生成仲胺基团,开环后的氧原子与氢原子形成羟基。第二步,另一个环氧基团开环与仲胺基团发生加成反应生成叔胺基团。根据交联反应机理,通过MS软件中的perl脚本建立了TGDDM与DDM的交联模型,具体交联细节见文献[20]。

图3 反应机理示意图

1.3 拉伸过程模拟

拉伸过程模拟通过美国Sandia 国家实验室开发的软件(LAMMPS）进行。首先将得到的交联环氧树脂模型转化为LAMMPS 软件的输入文件。然后对交联模型进行降温处理,具体策略如下：①2 000 步几何优化；②温度 600 K,200 ps 的 NVT 系综和 200 ps 的 NPT 系综的 MD 模拟；③以 30 K/200 ps 的降温速率,进行温度从600 K 降至300 K 的NPT 系综的 MD 模拟,模拟时间为 2 ns；④300 K,500 ps的NPT系综的MD模拟。

最后在NPT系综下对模型进行单轴拉伸模拟,在z方向以恒定的应变速率对模型施加应变,并保证体系在横向(x和y方向)的压强稳定在标准大气压。在整个模拟过程中使用COMPASS 力场,控温和控压方法分别采用Andersen 和Berendsen 方法,压强设定为0.1 MPa,时间步长为1 fs。Atom-based和Ewald方法计算范德华(vdW)和静电作用。

2 结果与讨论

2.1 不同交联度环氧树脂的力学性能

笔者通过MS软件中Perl脚本从交联反应截断半径0.55 nm 开始交联,到1.25 nm 结束,建立一系列不同交联度的交联环氧树脂分子模型,交联过程中交联密度、模型的密度和交联反应截断半径的关系列于表1,最终得到的最大交联度为90.21%。随后将表1 中的模型导入到LAMMPS 软件中,在300 K下采用deform命令对交联模型进行拉伸模拟。

表1 交联度、密度与截断半径的关系

图4 为交联度为90.21%模型在应变速率为5×108s-1时的应力-应变曲线。由图4 可见,在整个拉伸过程中,垂直与拉伸方向(图 2 中的x和y方向)上的应力σx和σy基本为0。当应变＜0.01 时,沿着拉伸方向(图2 中的z方向)的应力σz随应变的增加而线性增加,为弹性区,超过0.01 后材料进入屈服阶段,屈服强度可达110 MPa,在屈服应变之后,应力随应变的增加而下降,这与真实的实验现象一致[21]。

图4 交联度为90.21%模型在不同方向的应力-应变图

但当应变超过0.012 后,应力值又持续增长,这主要是由于交联环氧树脂的断裂伸长率较低,实验中当应变超过一定值后,交联分子链会发生断裂,但由于受到本文分子模拟中使用力场的限制,无法模拟分子链的断裂过程,所以应力值会持续增长。另外,由于实验测量时应变率一般为1×10-1s-1,远小于模拟时的应变速率(5×108s-1),由于“应变硬化”导致模拟值较实验值偏高,这与文献[22]等的研究结果一致。

图5为不同交联反应截断半径下得到的交联模型的应力-应变图,由图5和表1可见,交联反应截断半径越大,得到的模型交联度越大。各个模型的应力-应变图变化趋势比较相似,拉伸强度随交联度的增大而增大。

图5 不同交联反应截断半径下环氧树脂模型的应力-应变图

2.2 拉伸过程中结构变化情况

为了进一步观察拉伸过程中交联环氧树脂模型微观结构的变化,截取了不同应变时的结构快照如图6所示。从图6可知,由于交联结构不完整,模型会出现空洞,交联度越小,空洞的数目越多,形状越大；随着应变不断增加,模型沿z方向的长度不断增大,空洞也随之变大、变长。但由于模拟的应变较小(＜0.06)和三维网状交联结构的限制,模型的结构和空洞的形状变化并不十分明显。

图6 不同交联度环氧树脂模型在不同应变时的快照图

2.3 拉伸过程中的质量密度变化情况

为了进一步探索拉伸过程中交联环氧树脂模型微观结构变化情况,分析了交联度为90.21%环氧树脂模型在拉伸不同时刻即初始阶段(应变为0)、弹性阶段(应变为 0.005)、屈服阶段(应变为 0.010)、屈服过后应力下降阶段(应变为0.015)和应力上升阶段(应变分别为0.020,0.025,0.030)时,交联环氧树脂分子模型沿z方向的质量密度分布曲线,如图7所示。从图7可以看出,由于交联结构的限制,在整个拉伸过程中,模型沿z方向的质量密度不会发生较大的变化,与2.2分析结果一致。

图7 不同应变下交联环氧树脂模型z方向的质量密度分布图

2.4 交联环氧树脂的泊松比

通过计算横向应变与纵向应变的比值可以计算材料的泊松比。图8 为交联度为90.21%的环氧树脂模型在拉伸过程中横向应变与纵向应变的变化情况,在这里横向应变的值为x和y方向压缩应变的平均值,纵向应变为z方向的拉伸应变。通过线性回归求出泊松比为0.30,实验中典型热固性树脂在玻璃态下的泊松比为0.30～0.46,模拟值在实验范围内[23]。

图8 横向应变与纵向应变的关系

2.5 拉伸过程中各能量项分析

进一步分析了拉伸过程中各能量项的变化情况如图9所示,为了便于比较,所有能量值均为相对于初始模型能量的差值(Ei-E0),结合图4 的应力-应变图可以看出,在弹性区域键能、角能和二面角能随应变的增加而增加,在屈服阶段达到最大值,这说明在弹性区,模型的结构在拉力的作用下发生了小幅度的变形和扭转,达到屈服阶段后,结构变化趋于均匀一致；超过屈服应力后,由于交联结构的限制,使得交联模型的键,角和二面角不会发生明显变化,各能量项基本处于波动平衡状态。而vdW能的变化趋势与应力-应变图变化趋势大体一致,说明拉伸过程中,体系应变主要引起了vdW相互作用能的改变。

图9 应变过程中各能量项的变化情况

3 结论

通过分子动力学模拟方法,利用Perl 脚本通过Materials Studio分子模拟软件构建了不同交联度的环氧树脂模型,然后对交联模型进行单轴拉伸,研究了交联模型的拉伸性能,从拉伸过程中结构和能量的变化对拉伸机理进行了分析,主要结论如下：

(1)不同交联度的环氧树脂的应力-应变趋势比较相似,拉伸强度随交联度的增大而增大,当交联度为90.21%时,拉伸强度为110 MPa,泊松比为0.30,模拟值较实验值偏高,主要是由于模拟的应变速率较实验值高。

(2)由于三维交联网状结构的限制,拉伸过程中模型沿z方向的质量密度分布和微观结构变化并不十分明显。

(3)在拉伸过程中,在弹性区域键能、角能和二面角能随应变的增加而增加,超过屈服应力后,由于交联结构的限制,使得交联模型的键,角和二面角不会发生明显变化。体系应变主要引起了vdW相互作用的改变。