分子动力学模拟研究NR/BR力学性能和界面相互作用

江 浩，岳 红，刘 倩

（西北工业大学理学院，陕西 西安710129）

分子动力学模拟研究NR/BR力学性能和界面相互作用

江 浩，岳 红，刘 倩

（西北工业大学理学院，陕西 西安710129）

运用分子动力学（MD）方法，模拟研究了天然橡胶（NR）/顺丁橡胶（BR）混合胶体的力学性能，从微观上解释了NR和BR共混后性能得到改善的原因，并通过界面结合能和径向分布函数分析揭示了混合物组分之间的相互作用的本质。结果表明，NR/BR的平均结合能为21.35kJ/mol，除去程序自行产生的能量校正值（－3.64kJ/mol），结合能在数值上近似等于非键能（17.91kJ/mol），表明NR/BR组分界面间的相互作用主要由范德华力构成。

分子动力学；力学性能；界面结合能；径向分布函数

0 前言

橡胶自问世以来，随着工业生产的不断发展和对橡胶需求的日益增加，人们对橡胶的改性研究一直没有间断过。目前，对橡胶改性已经有了诸多改性手段，其中橡胶共混并用就是常用改性方法之一。长期以来，橡胶并用的配方设计主要依靠实验，缺乏基础理论指导，需要进行大量的实验研究，耗费大量的人力、物力和财力。与实验研究相比，分子模拟是一种更直接的研究方法，可以通过分子模拟技术先于实验预测和了解材料的性能，缩短新材料的研制周期，大大降低实验的成本和危险性。在橡胶工业中，NR和BR是应用最为广泛的并用体系，以其制备的轮胎胎面胶、胎侧胶以及橡胶输送带的覆盖胶等都具有良好的技术效果，其加工、形态和性能都已被广泛研究。通过分析NR/BR的共混原理有利于推进橡胶并用的改性手段。近年来，对共混物力学性能与界面间相互作用进行分子模拟已有相关研究［1－4］，但利用分子模拟研究 NR和BR共混物的微观机理，至今尚未有这方面的报导。

本文运用MD方法，使用Material Studio 4.0程序和 Compass力场［5－6］分别模拟计算BR和 NR/BR的原子运动轨迹，并比较了其力学性能（弹性系数、柯西压、模量和泊松比等）。求解出混合物组分间的界面相互作用（结合能），并以径向分布函数揭示了相互作用的本质，从微观尺度上解释了NR/BR共混后力学性能得到改善的原因。我们知道，通过分子模拟来解释NR和BR共混物的一切实验现象是不可能的，尤其是在如此复杂的界面相互作用的情形下有很多因素（如混合胶填料、添加剂的影响等）都不能考虑在内。所以，我们只希望得到能解释或预测特定性能的结果，就已经达到模拟的目的了。这些研究为本课题组后续研究“杜仲胶/NR/BR共混胶的结构与性能”和多种橡胶并用的分子模拟建立了一定的基础数据，所得结果可用于定性分析，具有一定的指导意义。

1 模型构建与方法

NR主要成分是顺式聚异戊二烯，BR主要成分是聚丁二烯，这两种聚合物分子式如图1所示。

图1 顺式聚异戊二烯和顺式聚丁二烯分子结构图Fig.1 Molecular structure of cis－polyisoprene and cis－polybutadiene

运用Materials Studio 4.0软件包中Visualizer模块分别构建聚合度为20的NR分子模型和聚合度为15的BR的分子模型，然后依据在298K和1个标准大气压下NR、BR相应的密度（NR相对密度为0.96g/cm3，BR相对密度为0.91g/cm3），利用 Amorphous cell模块分别搭建NR、BR纯物质和NR/BR共混物（质量比为2∶1［7］）三维周期性边界条件下的无定形高分子模型，如图2、3所示。

图2 BR与NR三维无定形结构模型Fig.2 Three－dimensional structural model of amorphous for BR and NR

图3 NR/BR三维无定形结构模型Fig.3 Three－dimensional structural model of amorphous of NR/BR blends

为了使体系更好的越过势能面上局部极小值之间的位垒，我们首先利用Smart Minimization方法对所构建的无定形分子模型进行结构优化，之后做退火模拟，使体系从200K开始逐渐升温，间隔为50K，选择NPT系综（恒定分子数、恒定压力以及恒定温度的系综），在体系压强为0.101MPa时对每个温度进行50ps的MD模拟，直至温度达到600K时为止；然后再以50K为间隔，将温度逐渐降至200K；进行2次循环，体系在经历了高温下的多次松弛，空间结构被逐渐优化，这一过程使建造的模型中产生的局部不合理结构基本消除。之后进行MD模拟计算分析力学性能等，具体模拟细节见表1.

表1 MD模拟细节Tab.1 Details of MD simulation

力学性能通过静态模型分析获得。关于应力－应变关系的广义虎克定律如式（1）所示。

式中 Cij——6×6弹性系数矩阵

σi——应力

εj——应变

由于应变能的存在，且是对称的，因此描述任意材料的应力－应变行为最多只需21个系数。在原子水平计算中，静态模型用维里方程计算内应力张量σ如式（2）所示。

式中 mi——原子质量

vi——原子速度

V0——无形变时体系的体积

施加于体系的应力使体系内粒子相对位置发生变化。对平行六面体（如模拟所取各边长分别为a、b、c的周期箱），若以列向量a0、b0、c0表征参比状态，以向量a、b、c表征形变状态，则应变张量为：

式中 h0——列向量a0，b0，c0组成的矩阵

A——度量张量

h——a，b，c组成的矩阵

因此通过计算拉伸和剪切形变的斜率，即可求得弹性系数矩阵，然后依据弹性力学原理，即可求得模量和泊松比等各向同性力学参数。各向同性材料的拉伸模量（E）、泊松比（v）、体模量（K）和剪切模量（G）可以用拉梅系数λ和μ 表示为式（5）～式（8）［10］。

2 结果与讨论

2.1 力学性能分析

拉伸模量、体模量和剪切模量均可以用来衡量材料的刚性强弱，即在外力作用下抵抗变形的能力。泊松比是材料力学和弹性力学中的名词，定义为材料受拉伸或压缩力时，材料会发生变形，而其横向变形量与纵向变形量的比值，是一无量纲的物理量。柯西压（C12－C44）通常可用来衡量体系的延展性，即材料发生形变而不产生裂缝的能力。一般来讲，一个具有较好延展性的物质，其（C12－C44）为正值，反之则为负值［11］。根据静态力学分析，对 MD模拟所得体系的平衡运动轨迹分析得到BR和NR/BR的弹性系数、模量和泊松比等，如表2所示。

表2 BR与NR/BR的弹性系数和有效各向同性力学性能Tab.2 Elastic constants and effective isotropic moduli for BR and NR/BR

对比表2中BR和NR/BR的弹性系数，发现后者的C11、C22和C33的数值相差不大，且C12、C13和C23组以及C44、C55和C66组内数值彼此更加接近，表明NR/BR更加趋近于各向同性弹性体［12］。另外，BR的拉伸模量、体模量和剪切模量值均大于混合物NR/BR的模量值，表明BR刚性较强；NR/BR共混后，体系模量下降，表明混合物刚性减弱，柔性增强，较易于发生弹性形变。BR的柯西压值大于NR/BR，即BR的延展性强于NR/BR，表明共混后胶料强度增大。这些都表明共混后体系的力学性能有较为显著的改善。

2.2 径向分布函数分析

径向分布函数g（r）是反映材料微观结构的特征物理量，表示在一个分子周围距离为r的地方出现另一个分子的概率密度相对于随机分布概率密度的比值。分子内径向分布函数可以提供模拟体系有序度的信息，一般出现大于0.3nm的峰表明分子链长程有序，属于结晶体系；而出现小于0.3nm的峰则表明分子链短程无序，属于无定形结构。分子间径向分布函数可以揭示非键原子间相互作用方式和本质，氢键作用范围为0.26～0.31nm，范德华作用范围为0.31～0.5nm。我们对整个平衡轨迹进行径向分布函数g（r）分析，做出NR，BR与NR/BR分子内与分子间径向分布函数，如图4和图5。

图4 BR和NR分子内及分子间径向分布函数Fig.4 Plots for pair correlation function computed for atoms of carbon of BR and NR

图5 NR/BR分子内及分子间径向分布函数Fig.5 Plots for pair correlation function Computedfor atoms of carbon of NR/BR blends

从图4可见，分子内径向分布函数中0.3～0.7nm范围内两种纯物质均没有出现强的峰，说明BR、NR的分子结构属于无定形结构；分子间径向分布函数则揭示，BR、NR分子间主要作用方式为范德华力作用。二者的径向分布函数曲线非常接近，根据相似相容的原理，可以说明其共混体系的相容性较好。从图5可见，分子内径向分布函数中0.3～0.7nm范围内两种纯物质均没有出现强的峰，说明BR/NR混合物属于无定形聚合物；分子间径向分布曲线显示NR/BR混合物内氢键作用不大，而形成了较强的范德华力。表明混合物内各组分间的相互作用主要由范德华力构成。

2.3 界面相互作用和结合能

界面相互作用的本质是各组分分子间的相互作用，而各组分间的相互作用直接影响其力学强度等性能。分子间相互作用能或分子间结合能可定量地表征分子间相互作用的强度［13－15］。为细致考察混合物组分间的相互作用，我们以稳定构型下各体系的总能量进行分子间相互作用能计算，NR与BR的平均相互作用能（Einter）可表示为：Einter＝Etotal－（EBR＋ENR），其中Etotal为NR/BR的平均总能量，EBR和ENR分别为BR和NR的平均单点能，定义结合能（Ebind）为相互作用能的负值，即Ebind＝－Einter。

经MD模拟计算，所得结果如表3。

表3 NR/BR界面的相互作用 kJ/molTab.1 Interface interaction of NR/BR blends kJ/mol

由表3可见，NR/BR 的平均结合能为21.35kJ/mol，若不考虑MS程序自行产生的能量校正值（－3.64kJ/mol），则结合能在数值上近似等于非键能（17.91kJ/mol），即为范德华力（17.71kJ/mol）和静电力（0.20kJ/mol）之和。其中只有微弱的静电力，范德华力起主导作用，表明NR/BR组分界面间的相互作用主要由范德华力构成，这与径向分布函数分析结果相一致。

3 结论

（1）与BR相比较，BR与NR共混可显著提高生胶强度，其动态拉伸疲劳性能极佳，其结论与实验结果一致。表明可以通过分子模拟来计算橡胶共混后的特定性能；

（2）比较NR与BR纯物质的分子内、分子间径向分布函数发现，二者分子链短程无序，属于无定形结构，同时揭示了NR、BR非键原子间主要作用方式为范德华力作用。分析NR/BR共混胶径向分布函数发现混合物内各组分间的相互作用主要由范德华力构成；

（3）在NR/BR共混胶中，界面结合能主要由范德华力构成。

［1］ Abou－Rachid H，Lussier L S，Ringuette S.On the Correlation Between Miscibility and Solubility Properties of Energetic Plasticizers/Polymer Blends： Modeling and Simulation Studies［J］.Propellants Explosives Pyrotechnics，2008，33（4）：301－310.

［2］ Karl G Wagner，Martin Maus，Andreas Kornherr.Glass Transition Temperature of a Cationic Polymethacrylate Dependent on the Plasticizer Content－simulation vs Experiment［J］.Chemical Physics Letters，2005，406（1/2/3）：90－94.

［3］ Jaidann M，Abou－Rachid H，Lafleur－Lambert X.Mode－ling and Measurement of Glass Transition Temperatures of Energetic And Inert Systems［J］.Polymer Engineering ＆Science，2008，48（6）：1141－1150.

［4］ Morita H，Tanaka K，Kajiyama T.Study of the Glass Transition Temperature of Polymer Surface by Coarse－Grained Molecular Dynamics Simulation［J］.Macromolecules，2006，39（18）：6233－6237.

［5］ H Sun.COMPASS：an Abinitio Force Field Optimized for Condensed Phase Applications Overview with Details on Alkane and Benzene Compounds［J］.Phys Chem，1998，102（38）：7338－7364.

［6］ S W Bunte，H Sun.Molecular Modeling of Energetic Materials：The Parameterization and Validation of Nitrate Esters in the COMPASS Force Field ［J］.Phys Chem B，2000，104（11）：2477－2489.

［7］ 姚钟尧，寇利涛.NR/BR并用胶的共混方法与微观结构和强度之间的相互关系［J］.世界橡胶工业，2007，34（7/8）：16－24.Yao Zhongyao，Kou Litao.NR/BR Rubber′s Blends Method and the Relationship Between Microstructure and Strength［J］.World Rubber Industry，2007，34（7/8）：16－24.

［8］ Andersen H C.Molecular Dynamics Simulations at Constant Pressure and/or Temperature［J］，Chem Phys，1980，72（4）：2384－2393.

［9］ Parrinello M，Rahman A.Polymorphic Transitions in Sin－gle Crystals：A New Molecular Dynamics Method［J］.J App Phys，1981，52（12）：7182－7190.

［10］ 张慧军，岳 红，刘 倩.形状记忆高分子材料性能评价的分子模拟研究［J］.材料导报，2010，24（9）：116－119.Zhang Huijun，Yue Hong，Liu Qian.Performance Evaluation of Shape Memory Polymer by Molecular Simulation［J］.Materials Review，2010，24（9）：116－119.

［11］ 夏 露.高能材料结构和性能的分子动力学模拟［D］.苏州：苏州大学化学化工学院，2008.

［12］ Xu Z L.Elasticity Mechanics［M］.Beijing：Higher Education Press，1990：328－376.

［13］ Radovan T，Alessandro C，Marco F.Computer Simulation of Polypropylene/Organoclay Nanocomposites：Characterization of Atomic Scale Structure and Prediction of Binding Energy［J］.Polymer，2004，45（23）：8075－8083.

［14］ Paolo C，Giulio S，Sabrina P.Many－scale Molecular Simulation for ABS－MMT Nano Composites：Up Grading of Industrial Scraps ［J］.Microporous and Mesoporous Materials，2008，107（1/2）：169－179.

［15］ Subhra M，Golok B，Nando K.Mechanical and Dynamic Mechanical Properties of Miscible Blends of Epoxidized Natural Rubber and Poly（ethylene－co－acrylic acid）［J］.Polymer，1996，37（24）：5387－5394.

Molecular Dynamics Simulation of Mechanical Properties and Surface Interaction for NR/BR Blends

JIANG Hao，YUE Hong，LIU Qian
（School of Science，Northwestern Polytechnical University，Xi’an 710129，China）

Molecular dynamics simulation was applied to investigate the mechanical properties of NR/BR blends，and the reasons of their superior performance were interpreted in microstructure.Finally，by calculating binding energy and pair correlation function，the interface interaction and its essence of the blends were elucidated.The results showed that the average binding energy of NR/BR blends was 21.35kJ/mol，the number of average binding energy was approximately the nonbond energy（17.91kJ/mol）except energy calibration （－3.64kJ/mol）produced by program，which meant that interaction between NR and BR was formed by Vander Waals force.

molecular dynamics；mechanical property；interface binding energy；pair correlation function

TQ332

B

1001－9278（2012）05－0064－05

2011－12－12

联系人，russelljiang＠163.com

（本文编辑：刘本刚）