基于分子动力学模拟的天然橡胶黏弹性材料力学行为

徐业守 徐赵东 郭迎庆 黄兴淮 戴 军 盖盼盼

(1东南大学混凝土及预应力混凝土结构教育部重点实验室, 南京 211189)(2东南大学土木工程学院, 南京 211189)(3南京林业大学机械电子工程学院, 南京 210037)(4南京东瑞减震控制工程有限公司, 南京 210033)

黏弹性阻尼减震在实际工程中应用广泛,是一种比较成熟的消能减震技术.黏弹性阻尼器的减震性能和耗能能力主要取决于材料的力学性能[1].目前,学者们已在黏弹性材料和阻尼器研究方面取得了大量的研究成果[2-3].

Shi等[4]研究了乙烯-乙酸乙酯橡胶与三元乙丙橡胶共混(EVM/EPDM)和乙烯-乙酸乙酯橡胶与丁腈橡胶共混(EVM/NBR)2种共混体系的黏弹特性.吴波等[5]采用试验方法研究了环境温度、加载频率、激励幅值和滞回加载圈数等对阻尼器性能的影响.Xu等[6]对不同黏弹性材料种类、加载幅值和频率下阻尼器的动态力学性能和减震能力进行了研究.黏弹性材料和阻尼器的力学性能和耗能能力受环境温度、激励幅值和加载频率影响较大,选择合适的力学模型来描述材料和阻尼器的力学性能尤为重要.常规力学模型包括Kelvin模型、Maxwell模型、广义黏弹性模型和标准线性固体模型等.Tsai等[7]提出了一种可以表征外界动态加载位移幅值、频率和环境温度对阻尼器力学性能和耗能能力影响的力学模型.徐赵东等[8]考虑微观分子结构对黏弹性材料宏观力学性能和耗能性能的影响,提出了黏弹性材料的微观分子链结构力学模型.

然而,针对黏弹性材料和阻尼器力学性能的研究相对较少.目前,关于基体橡胶材料超弹特性的研究方法主要分为2种:① 基于连续介质力学的唯象理论,如Neo-Hookean模型、Yeoh模型、Mooney-Rivlin模型和Ogden模型等;② 从橡胶分子网络结构出发,结合统计力学从微观层面对橡胶类材料力学行为进行理论分析[9].

鉴于微观分子结构的复杂性,一些微观尺度的黏弹性材料试验存在试验室内操作难度大、成本高等问题.而基于计算机技术的分子动力学模拟则提供了一种新的研究方法.Kitamura等[10]采用分子动力学模拟方法对聚丙烯橡胶在单轴拉伸、单轴压缩、等双轴拉伸和等双轴压缩等加载条件下的力学行为进行了研究.Payal等[11]基于全原子分子动力学模拟研究了聚丁二烯橡胶的拉伸变形,分析了拉伸过程中材料的破坏过程.

常用的分子动力学模拟软件包括LAMMPS、MS和VMD等.本文以天然橡胶为例,采用Dreiding力场,分别基于MS和LAMMPS软件进行材料的微结构建模和单轴拉伸加载条件下的分子动力学模拟.考虑分子链聚合度参数p、环境温度T和加载应变率R的影响,对黏弹性材料的力学性能进行分析.

1 模型信息

本文以天然橡胶为例,对黏弹性材料的弹性力学行为进行建模分析.模型的建立主要借助MS软件来实现.

1.1 模型建立

天然橡胶是以聚异戊二烯为主要成分的高分子材料.借助MS软件中Visualizer 模块,建立异戊二烯分子单体,生成不同聚合度的天然橡胶分子链,对单根分子链进行能量最小化处理.采用Amorphous Cell模块建立天然橡胶的周期模型,密度为0.92 g/cm3[12].聚合度为40时天然橡胶模型建模过程如图1所示.采用msi2lmp.exe工具将已建好的天然橡胶周期模型转化为LAMMPS可读取的模型数据,以备后续模拟运算.

图1 天然橡胶黏弹性材料周期模型的建立

1.2 力场信息

Dreiding力场可以较好地模拟大分子聚集体的结构及各项性质,常用于高分子材料的模拟计算[13].研究黏弹性材料力学性能的过程中,采用全原子模型,并选用Dreiding力场进行计算.参照文献[13],Dreiding力场中键长拉伸势能为

Eb(rm)=0.5kb(rm-rm0)2

(1)

键角转动势能为

Ea(θ)=0.5kθ(θ-θ0)2

(2)

二面角转动势能为

Ed(φa)=0.5Vφ{1-cos[nφ(φa-φa0)]}

(3)

非键结势能为

(4)

则系统的总势能可以表述为

Et=Eb+Ea+Ed+En

(5)

式中,kb和kθ分别为键长和键角势函数的刚度参数;rm为2个原子间的距离;θ为相邻2个化学键的夹角;rm0和θ0分别为相邻2个化学键间的平衡键长和平衡角度;φa为二面角;Vφ为二面角发生旋转的阻力常数;nφ为周期常数;φa0为二面角的平衡角度;σr为零势能间距;εr为势能最低点与零势能之差;rc为截断半径.

天然橡胶周期模型基本信息包括原子类型、键长、键角、二面角、非键结势能类型等.键长、键角势能的计算均采用简谐振动形式,见式(1)和式(2).二面角势能的表达式如式(3)所示,非键结势能采用L-J/cut 类型进行定义.

同样,结合天然橡胶的分子构型简式,参考文献[13],可以得到分子链结构中键角转动势能和二面角转动势能的类型分别为8和11种.天然橡胶分子链结构中非键结势能类型及其参数见表 1.

表1 天然橡胶非键结势能参数取值

2 分子动力模拟

分子动力学模拟过程中,采用LAMMPS编程并进行计算.根据从MS中导入的模型数据,建立天然橡胶的不定型态周期模型.在Dreiding力场下,首先采用NVE系综运算,对模型运行1×104步动态模拟,温度为500 K,步长为0.5 fs,使模型得到优化和平衡;其次,运行1×104步NPT系综动态模拟,时间步长为1 fs,使模型系统更加合理;然后,在NPT系综下运行5×104步动态模拟,时间步长为0.5 fs,使系统温度达到目标温度;进而采用NPT系综,在目标温度下运行5×104步动态模拟,时间步长为0.5 fs,使系统在目标温度下达到稳定状态;最后,对模型进行单轴拉伸性能模拟.图2为天然橡胶黏弹性材料在拉伸加载前后的状态对比.模拟温度为T=100 K,分子链的聚合度参数为p=40,拉伸应变加载率为R=109s-1.由图可知,材料沿拉伸方向伸长,垂直于拉伸方向的尺寸逐渐减小,在低温与高加载率条件下材料表现出塑性破坏特性.图2(b)中可以清楚地观察到材料分子链间相互分离,说明断裂破坏已经发生.由此表明,该模型可以准确模拟天然橡胶材料的力学行为.

(a) 加载前

2.1 应力-应变分析

天然橡胶黏弹性材料的标准拉伸曲线见图3.材料分子链聚合度参数p=40,加载率R=108s-1,环境温度T=250 K.由图可知,在温度较低、加载率较高的情况下,材料表现出明显的塑性特性,经历了弹性、屈服、软化、增强以及最终的颈缩与断裂阶段.初始时刻,材料表现为线性弹性;到达屈服点后,出现应变软化现象;应力下降到局部最小值后出现转折并开始逐渐增大;应力增加到一定程度后进入破坏阶段,截面面积逐步减小,直至发生断裂.本文重点分析分子链聚合度参数p、环境温度T和加载率R对黏弹性材料应力-应变行为的影响.

图3 黏弹性材料标准应力-应变曲线

图4给出了材料的应力-应变曲线随天然橡胶分子链聚合度参数的变化规律.由图可知,相同密度下,材料的应力峰值随着聚合度参数p的增加而增加.以图4(a)为例,p=10，20，40时,应力峰值分别为71.7、79.8、99.8 MPa.此外，p=20,40时,分子链聚合度较高,初始状态下物理缠结和分子间作用力较强,加载过程中应力峰值较高,较早发生塑性破坏；p=10时分子链聚合度较低,物理缠结作用影响较小,材料应力发展平缓.由图4(b)可知,在T=400 K的高温条件下,p=10，20，40时,材料均未发生明显的塑性变形.

(a) R=108 s-1,T=250 K

(b) R=109 s-1,T=400 K

图5给出了材料的应力-应变曲线随系统温度的变化规律.由图可知,随着温度的增加,材料的应力峰值明显下降.以图5(b)为例,T=100，250，400 K时,应力最大值分别为134.9、103.9、83.3 MPa.天然橡胶玻璃化转变温度在200～213 K之间,故图5(a)中R=109s-1时,100和250 K条件下,材料进入塑性阶段；而在400 K时,材料表现出一定的黏弹性特性.由图5(b)可知,R=1010s-1时不同温度下材料均较早屈服.

图6为材料应力-应变曲线随加载率的变化规律.由图可知,随着加载率增加,材料应力峰值呈现增大趋势.图6(a)中,R=108，109，1010s-1时,应力幅值分别为122.5、122.7、134.9 MPa.由于系统温度较低,聚合度较高,图6(a)和(b)中天然橡胶材料在R=108，109，1010s-1时均较早达到屈服点,较低应变下即产生塑性变形.

(a) p=40,R=109 s-1

(b) p=40,R=1010 s-1

(a) p=40,T=100 K

(b) p=40,T=250 K

2.2 能量分析

图7给出了单轴拉伸过程中天然橡胶材料系统的非键结势能增量ΔEn、键长拉伸势能增量ΔEb、键角转动势能增量ΔEa、二面角转动势能ΔEd以及系统总势能增量ΔEt随应变的变化规律.

(a) p=40,R=109 s-1,T=100 K

(b) p=20,R=108 s-1,T=250 K

由图7(a)可知,在材料的拉伸过程中,系统总势能的增加主要归因于非键结势能,键长拉伸势能、键角和二面角转动势能对总势能的变化贡献不大.系统非键结势能在弹性阶段随应变增加而线性增大,到达屈服阶段后趋于稳定,在颈缩破坏阶段开始下降.由图7(b)可知,在弹性和屈服阶段,非键结势能影响占据主导作用,进入应变软化阶段后,非键结势能下降明显,对系统总势能的贡献与键长拉伸势能、键角和二面角转动势能相当.

2.3 自由体积分析

系统自由体积变化对材料性能影响显著.研究中认为系统内原子体系的体积不发生变化,采用周期模型总体积的变化来等效表征自由体积的改变,分别研究了分子链聚合度、系统温度和应变率变化对天然橡胶黏弹性材料自由体积改变的影响.

由图8(a)可知,p=10时,由于分子链较短,自由体积增量随应变增加变化不大.p=20，40时,分子间物理缠结和范德华作用随应变的增加而减弱,自由体积增量增加明显.相同应变下,自由体积增量随参数p增加而增加,这是因为p值越大,分子链越长,变形前材料物理缠结和范德华作用越显著,拉伸加载后物理缠结和分子作用力受到的破坏越多,对自由体积变化的影响也越大.

(a) R=108 s-1,T=250 K

(b) p=40,R=109 s-1

(c) p=40,T=250 K

由图8(b)可知,T=400 K时,系统温度高于玻璃化转变温度,自由体积增量随应变变化不明显.T=100，250 K时,自由体积增量均随应变增加而显著增大.同时,系统自由体积增量随温度的升高而降低.

由图8(c)可知,不同加载率下,系统自由体积增量均随应变的增加而增加.相同应变下,较高的应变加载率具有更高的自由体积增量.

系统自由体积增量随应变变化曲线的斜率在颈缩破坏起始点附近突然增大,究其原因在于颈缩破坏阶段材料颈缩截面面积不断减小,材料逐步发生破坏,导致黏弹性材料中孔隙体积大量增加.

3 结论

1) 分子动力学模拟方法可以较好地分析材料微结构和加载条件对黏弹性材料力学性能的影响.

2) 黏弹性材料的应力和弹性模量随分子链聚合度和加载率的增加而增加,随温度升高而明显降低.

3)较大加载率和低温下,黏弹性材料表现出塑性破坏的特点.玻璃化转变温度附近一定温域内,材料表现出黏弹性特性.

4) 在拉伸加载过程中,非键结势能对系统总势能的变化起主导作用.系统自由体积增量随分子链聚合度和加载率的增加而增加,随温度的升高而降低.