一种用分子动力学模拟诱导聚合物链取向结晶以及拉伸力学破坏的方法

摘要

本发明公开一种用分子动力学模拟诱导聚合物链取向结晶以及拉伸力学破坏的方法。建立聚合物的分子模型，设定原子的初始位置、速度和力场。体系总能量包含动能和势能两部分，其中动能和温度有关，而势能主要包括成键和非键相互作用，即键伸缩能、键角弯曲能、二面角扭转能和范德华相互作用能。采用周期性边界条件，依次对体系进行能量优化、NVT和NPT平衡，进行非等温结晶模拟。用GROMACS软件对体系进行单轴恒力拉伸取向结晶；然后用LAMMPS软件恒速拉伸，研究结晶聚合物的应力应变，及对温度和拉伸速率的依赖性。本发明利用计算机模拟方法，克服了实验局限，实时观察到取向结晶和拉伸过程中微晶的演变，结论与实验数据吻合。

说明书

技术领域

本发明属于高分子材料技术领域，更具体地说，涉及一种计算机模拟诱导聚合物链取向结晶以及拉伸力学破坏的方法。

背景技术

半结晶聚合物具有力学性能优异，加工性能良好，性价比高等优点，被广泛应用于各个领域。半结晶聚合物材料要经过加工成型之后才能成为制品，熔体在加工成型时都会不可避免地引入流动场，比如双轴拉伸和薄膜吹塑过程，纤维挤出和塑料成型过程等。流动场的存在能够加速结晶动力学，显著地增强材料的力学性能，并且会诱导生成特殊的晶体形态和新的晶型，其中半结晶聚合物最明显的形态变化是从无规取向的球晶向高度取向的串晶（shish-kebab）的转变，串晶的存在能够大幅度提高半结晶聚合物材料的各项性能。半结晶聚合物的性能主要由产品的结晶形态决定，链结构的多样性和结晶条件的复杂性影响了半结晶聚合物的结晶形态，因此，研究半结晶聚合物的结构-性能之间的关系可以为开发高性能半结晶聚合物产品提供有效的理论指导。由于实验手段和经典理论模型的限制，对拉伸诱导聚合物取向结晶，以及取向结晶的拉伸力学破坏缺少分子尺度的解析，而分子模拟有助于从本质上理解聚合物结晶的微观机理和复杂性，其中分子动力学模拟（MolecularDynamics, MD）是从分子尺度上研究聚合物结晶行为的一个十分重要的研究手段。

分子动力学模拟方法是一种研究凝聚态系统强有力的计算机模拟方法。分子动力学方法中粒子的运动遵循牛顿第二定律，体系的演变接近真实粒子的运动规律，可以同时预测体系的热力学和动力学特征，对经典理论计算和实验研究进行了有力的完善和补充。

分子动力学方法可以通过计算原子间的相互作用势求解体系内每个原子所受的力，用数值方法求解牛顿动力学方程组，模拟预测每个原子在系统中随着时间的推移产生的运动规律。并且可以获取体系重要的参数变化，包括构象变化，配体配合和蛋白质的折叠方式，可以在飞秒的时间尺度下显示所有原子的位置。分子动力学模拟方法通常与各种实验技术相结合，包括X射线结晶学、核磁共振技术和电子顺磁共振等方法。

GROMACS和LAMMPS是一种经典的分子动力学代码，它可对处于液态，固态或气态的粒子团进行模拟计算，在对原子，聚合物，生物，固态（金属，陶瓷，氧化物），颗粒，粗粒度或宏观系统进行建模时，可以向体系施加各种原子间电势（力场）和边界条件。其模拟体系十分广泛，少至几个粒子，多至数百万或数十亿的2D或3D结构，都可以对其进行模拟。

本文采用联合原子模型进行分子动力学模拟，即将半结晶聚合物链中的每个亚甲基都看作一个整体，不计算其内部的相互作用，从而大大简化了计算过程。本发明模拟体系的分子链体系较大，分子链超过盒子边界，所以使用周期性边界条件。整个模拟过程中，盒子的大小会发生改变，但体系的密度始终不变，施加压力后，聚合物处于熔体的状态。在数据处理过程中，以联合原子单元的坐标为准，忽略了周围空位的变化。

发明内容

本发明的目的在于提供一种模拟诱导聚合物链取向结晶以及拉伸力学破坏的方法，进而从微观机理方面解释拉伸诱导聚合物取向结晶，以及拉伸力学破坏中凝聚态结构的演变，克服实验表征无法揭示的现象原因，进一步建立聚合物材料分子链结构—流场控制—结晶行为—宏观性能之间的关联。本发明用分子动力学模拟研究流场下半结晶聚合物的结晶行为和成核过程。

本发明的技术目的通过以下技术方案实现。

一种模拟诱导聚合物链取向结晶以及拉伸力学破坏的方法，建立聚合物模型并模拟其拉伸结晶行为，以及拉伸力学破坏过程，具体按照下述方案进行。

步骤1，建立聚合物联合原子模型，得到聚合物初始坐标、键长和键角参数,也就是势函数的选取，通常采用的势函数是LJ势函数，不同的物质状态可以用来描述不同的势函数。模型的势函数确定之后，就可以根据物理学规律求取模拟中的守恒量。对于初始条件的设定，需要知道粒子的初始位置和速度。实际上也不需要精确选择待求系统的初始条件，由于模拟时间足够长，系统就会忘记初始条件，但是合理的初始条件可以加速系统趋向平衡的过程，以便获得好的精度，初始速度从玻尔兹曼分布随机抽样得到。

步骤2，对模型设置合适的对应力场，其中体系势能包括如下四个方面，键伸缩能（

（1）键伸缩能

分子中互相成键的原子形成的化学键的键长并非是固定不变的，而是在其平衡位置附近呈现小幅度的震荡，描述这种作用的势能称为键伸缩能

其中，

（2）键角弯曲能

分子中连续键接的三个原子形成键角，同键伸缩项一样，其并非维持在固定值，而是在其平衡值附近小幅度震荡，描述此种作用的势能为键角弯曲能

其中，

（3）二面角扭转能

分子中连续成键的四个原子形成二面角，描述此种内扭转的势能为二面角扭转能

其中，

（4）范德华相互作用能

在分子力场中，如果i、j两原子属于同一个分子，但其之间被三个或者三个以上的化学键所分隔开，或分属于两个不同的分子，那么i、j原子间存在范德华相互作用能

其中，

步骤3，对模型设定相应的周期性边界条件，且所有的原子运动方程采用Verlet速度算法，用标准Nosé-Hoover系综保持体系的温度稳定，采用标准的Parrinello-Rahman系综保持体系的压力恒定。

步骤4，趋于平衡计算，为使系统趋于平衡，在整个过程中需要增加或者从系统中移除能量，直至能够持续给出确定的能量值，此时，系统已经达到平衡，我们称达到平衡的这段时间为弛豫时间。在模拟过程中，时间步长的选择非常重要，适当时间步长的选取既可以减小误差又可以大幅度缩短系统模拟的弛豫时间。

具体方案：将模型导入GROMACS中进行体系能量最小化，在固定体积以及T* = 8.0温度下平衡（NVT），再在T* = 8.0，P=0.1 MPa下平衡（NPT），最后温度T* = 8.0降至3.0并固定拉力进行拉伸结晶，通过这四个阶段确保体系接近真实的温度和压力，最终输出所有原子的坐标，能量等数据，进行宏观物理量的计算以及结晶行为分析，并得到相应的结晶曲线（如图1）和成核过程的变化情况（如图3）。

步骤5，将降温结晶的模型导入LAMMPS中，在NPT系综中以固定的速率进行拉伸模拟，最后输出拉伸过程中原子坐标、应力、应变和能量等数据，观察模型的构象变化并对拉伸行为进行分析，得到应力~应变曲线和能量演化曲线。改变温度和速率对模型拉伸，则可以分析应力应变行为和能量演化对温度及拉伸速率的依赖性。

本发明的计算机模拟是在实验基础上通过基本原理建立起的一套模型与算法，从而计算出合理的分子结构与分子行为，运用分子模拟，更加清晰地从分子尺度描述和解释在实际流场下聚向结晶行为，以及拉伸过程中凝聚态结构的演变，对实际操作有一定的指导意义。

附图说明

图1是分子量分布为20（MWD20）的40000C的聚乙烯链（PE），由T* = 8.0降至T* =3.0的结晶度，有序度，能量和反式二面角占比关于温度/运行时间的关系图。

图2是是分子量分布为20（MWD20）的40000C的聚乙烯链（PE），由T* = 8.0降至T* =3.0的拉伸诱导结晶图。

图3是分子量分布为1、2.5、10、15和20的40000C的聚乙烯链（PE），由T* = 8.0降至T* = 3.0的结晶度曲线。

图4是分子量分布为20（MWD20）的40000C的聚乙烯链（PE）在不同拉力下（k

图5是具有共聚单体为己烯（-B）的支链分布在长链上的不同支链含量（2C/1000C、4C/1000C、6C/1000C、8C/1000C）的双峰聚乙烯模型的结晶度变化曲线。

图6是具有共聚单体己烯（-B）的支链分布在长链和短链上的不同支链含量（0、2C/1000C、4C/1000C、6C/1000C、8C/1000C）的双峰聚乙烯模型的结晶度对比图。

图7是具有共聚单体己烯（-B）的支链分布在长链和短链上的不同支链含量（0、2C/1000C、4C/1000C、6C/1000C、8C/1000C）的双峰聚乙烯模型的反式二面角占比对比图。

图8是20条主链链长为1000C的聚乙烯链（PE），在T* = 1.0，10

图9是20条主链链长为1000C的聚乙烯链（PE），在T* = 1.0，在10

图10是20条主链链长为1000C的聚乙烯链（PE），在T* = 1.0及T* = 3.0，以10

具体实施方式

下面结合具体实施例进一步说明本发明的技术方案。

实例1。

建立分子量分布为20（MWD=20）的40000C的聚乙烯链（PE），依次进行体系能量最小化，在固定体积以及T* = 8.0温度下平衡（NVT），再在T* = 8.0，P=1bar下平衡（NPT），最后温度T* = 8.0降至T* = 3.0并固定拉力的单轴拉伸结晶，输出整个过程的原子坐标并进行统计分析。

如图1和图2，结合结晶度和反式二面角含量的变化趋势以及直观的成核结晶行为，整个拉伸过程可以分为三个阶段，I：shish前驱体的形成，结晶度几乎稳定在0左右；II：从前驱体到shish核的演化，随着拉伸过程的进行，连续的反式二面角的比例和

实例2。

建立五种具有不同分子量分布（MWD）的聚乙烯链模型，分别为1,2.5,10,15,20的40000C的聚乙烯链（PE），依次进行体系能量最小化，在固定体积以及T* = 8.0温度下平衡（NVT），再在T* = 8.0，P=1bar下平衡（NPT），最后温度T* = 8.0降至T* = 3.0并固定拉力的单轴拉伸结晶，输出整个过程的原子坐标并进行统计分析。如图3，进一步证实了结晶行为，并得到分子量分布对结晶行为的影响。

实例3。

建立分子量分布为20（MWD=20）的40000C的聚乙烯链（PE），依次进行体系能量最小化，在固定体积以及T* = 8.0温度下平衡（NVT），再在T* = 8.0，P=1bar下平衡（NPT），最后温度T* = 8.0降至T* = 3.0，用不同的固定拉力（k

实例4。

建立一根长链（主链链长为10000 C）和10根短链（每根短链的主链链长均为1000C），在长链主链上设计了共聚单体为己烯（-B）的不同支链含量的模型。支链含量逐渐由2C/1000C增加到4C/1000C、6C/1000C、8C/1000C（每1000个主链C上插入8个支链C）。该模型依次进行体系能量最小化，在固定体积以及T* = 8.0温度下平衡（NVT），再在T* = 8.0，P=1bar下平衡（NPT），最后温度T* = 8.0降至T* = 3.0并固定拉力的单轴拉伸结晶，输出整个过程的原子坐标并进行统计分析。

如图5，结晶度曲线继续验证了结晶行为，且支链含量增加对聚合物结晶有阻碍作用。5 ns之前，结晶度增加较缓慢，说明在拉力的作用下体系只产生了较少的反式富集片段，结晶度较小。约5 ns之后，大量的shish核出现并且快速生长为大的晶区，体系的结晶度迅速上升至较大值，随后晶区的生长速度减缓，结晶速率略有降低，结晶度达到最大值。

实例5。

建立一根长链（主链链长为10000 C）和10根短链（每根短链的主链链长均为1000C），在长链主链上设计了共聚单体为已烯（-B）的不同支链分布的模型。一是支链全部均一的分布在长链上（L

如图6和图7，当支链为己基时，支链分布在短链上的聚乙烯链模型（L/10S

当支链分布在短链上时，长链更易产生反式富集片段，这些片段易于连接在一起形成晶核，含有支链的短链较难进入长链成核的生长前沿，形成规整片晶，系带分子浓度相对较低，结晶度较高。当支链分布在长链上时，易产生反式富集片段的长链由于支链的存在相对难以形成晶核，短链易进入长链成核的生长前沿，不易形成规整的片晶结构，系带分子浓度较高，结晶度较低。

实例6。

建立20根主链链长为1000C的聚乙烯链（PE），依次进行体系能量最小化，在固定体积以及T* = 8.0温度下平衡（NVT），再在T* = 8.0，P=1bar下平衡（NPT），最后温度T* = 8.0降至T* = 1.0降温结晶。结晶的聚乙烯链在T* = 1.0，零压条件下（NPT）以10

如图8，显示了聚乙烯链在T* = 1.0下以10

实例7。

建立20根主链链长为1000C的聚乙烯链（PE），依次进行体系能量最小化，在固定体积以及T* = 8.0温度下平衡（NVT），再在T* = 8.0，P=1bar下平衡（NPT），最后温度T* = 8.0降至T* = 3.0和T* = 1.0分别降温结晶。结晶的聚乙烯链在T* = 1.0和T* = 3.0，零压条件下（NPT）分别以10

如图9，聚乙烯体系以10

如图10，显示了由20条链组成的聚乙烯在T* = 1.0, T* = 3.0温度下以10

不同温度和变形速率下的应力应变行为表明，随着应变速率的增加，应力-应变响应与温度降低时的情况相似。

以上对本发明做了示例性的描述，应该说明的是，在不脱离本发明的核心的情况下，任何简单的变形、修改或者其他领域技术人员能够不花费创造性劳动的同等替换均落入本发明的保护范围。