PVDF单轴拉伸的第一性原理和分子动力学模拟

曾凡林，孙 毅，刘一志

（哈尔滨工业大学航天科学与力学系，150001 哈尔滨）

PVDF单轴拉伸的第一性原理和分子动力学模拟

曾凡林，孙 毅，刘一志

（哈尔滨工业大学航天科学与力学系，150001 哈尔滨）

为得到聚偏二氟乙烯(PVDF)的大变形拉伸特性，应用第一性原理方法模拟了β相PVDF分子链的在单轴拉伸下的能量及构型变化，得到了β相分子链拉伸过程中发生断裂时的临界键长、临界断裂能量、临界拉力以及分子链弹性系数的变化规律.接着用一种新颖的方法构造了无定形PVDF的立方元胞模型并用于单轴拉伸的分子动力学模拟.模拟中分别考虑链不被拉断和链允许被拉断两种情况.第1种情况下，拉伸应变率分别为1010，109和108s－1，得到了应力应变关系和构型变化;第2种情况下，拉伸应变率分别为109s－1，5×108s－1和准静态，得到了应力应变关系和断裂前后的构型.同时分析了一些关键的形变信息和影响模拟结果的一些关键因素.模拟结果表明，不仅对于理解PVDF的拉伸变形机制具有较重要意义，而且对于其他无定形材料的拉伸模拟具有参考价值.

PVDF;单轴拉伸;第一性原理;分子动力学;断裂

作为一种被广泛应用的压电聚合物，聚偏二氟乙烯(Poly(vinylidene difluoride)，PVDF)还应用于对热稳定性、耐化学腐蚀性和耐磨性有特殊要求的领域.因为有与其他聚合物相比更卓越的力学性能以及突出的压电、热电性能，PVDF一直以来都引起广大研究者的注意［1－2］.PVDF 的压电性能主要源自其分子结构中的β相和γ相.相比于其他的压电材料例如 PZT，虽然 PVDF的力——电耦合效应低于PZT，但是它的可塑性要高出10倍左右［3－4］;此外作为传感材料更重要的一个特征参数PVDF的g要比PZT高20倍左右，而且它能被制成厚度＜10 μm的薄膜，振动质量极小，具有比PZT大得多的阻尼值，因此具有产生极短波长所需的动态特性.正因为如此，PVDF正被尝试应用于各种航天智能系统结构中，例如可充气式展开天线等，有望在各类航天器智能结构中得到应用［5－6］.

单轴拉伸是制备β相PVDF的传统且有效的方法.拉伸过程中，应变高达400%，使得PVDF的分子链由蜷曲的无定形态转变为直链态.然后在强电场的作用下利用CH2和CF2基团的极性，使得分子链中原先呈交错状态的CH2和CF2的相对位置发生改变，所有的H原子位于一侧，而所有的F原子则转向另一侧，使材料产生宏观的整体极性，从而具有压电性能.实验结果表明，如果分子链之间交联度较小的话，较小的拉应变即会使链和链发生分离从而产生链间断裂.因此在制备无定形的PVDF时，尽可能多地在分子链间形成相互交联的网状结构是拉伸法成功的关键.但这样会产生另外一个问题，即是相互交联的分子链中不可避免地会存在死缠接.一旦这种情况出现，即使较小的拉伸应变也会使得分子链断裂，从而在断裂处形成空位或银纹，极化时非常容易被击穿，使得材料失效.因此研究无定形的PVDF分子链间的死缠接以及拉伸过程中的分子链断裂情况，对于PVDF的制备和应用是非常有意义的.

因为传统的实验手段很难观察到PVDF内部的分子链结构，更难以表征拉伸过程中链的解缠接和断裂过程，因此本文采用计算机模拟的方法.分子链的断裂需要了解单个链在拉伸过程的构型信息和能量信息从而得到临界断裂值，本文采用第一性原理来模拟单链的拉伸来进行;链间的解缠结或断裂过程则通过大规模的模拟来进行，本文采用分子动力学模拟方法.第一性原理对单链的模拟结果直接应用于大规模的分子动力学模拟中，两者结合成功地对PVDF的大变形拉伸过程进行了模拟.该方法对于其他材料的大变形分子动力学模拟同样具有参考价值.

1 PVDF单链拉伸的第一性原理模拟

1.1 模拟模型及方法

本文依据文献［7］构造了β相PVDF的分子链，如图1所示.模型采用周期性边界条件以利于实施链长方向的轴向加载.为了减少周期性边界所产生的镜像分子链对当前分子链的影响，本文将模拟元胞的垂直于链长方向(x和y方向)的尺寸设置成远远大于截断半径的1.5 nm.为了使得变形只发生在链长方向，将元胞垂直于链长方向(x和y方向)的尺寸固定.为了减少计算量，模型中只取到两个分子链的重复单元，整个模型含有12个原子，其优化后的x、y、z方向的尺寸分别为1.5 nm(固定)、1.5 nm(固定)和 0.509 9 nm.

图1 β相PVDF分子链的构型及模拟模型

链长方向的单轴拉伸过程的第一性原理模拟采用Accelrys公司的 Materials Studio软件中的Castep模块来进行.通过改变元胞的z向尺寸来实现加载过程，每步加载应变为1%，加载步数不限，直至分子链发生断裂.加载过程中x和y方向的尺寸保持不变，具体的模型及加载信息参考文献［8］.本文在计算中采用了Vanderbilt超软赝势和 LDA 的 Ceperley-Alder形式的函数［9］，布里渊区采用2×2×5的 k点取样(其收敛性设为1.0 e－8eV/原子)来进行数值积分.元胞的能量计算中收敛标准设置为1.0 eV/原子.整个模拟过程参考文献［10］方法.

考虑到模型只包含单个分子链并且分子链之间的相互影响忽略不计，此外分子链只在链长方向受到约束而x和y方向形变自由，因此拉伸过程中能量的变化可看作全部由链长方向的变形所贡献.这样分子链所受到的拉力F可表示为F=－ (∂ΔE/∂εzz)/d0.其中:d0为元胞的 z向原始尺寸;εzz为当前的拉伸应变;ΔE为当前拉伸加载后相对于初始结构的能量变化.同样，C-C键的弹性系数k可表示为:k= － (∂ΔF/∂εzz)/d0.其中F 为计算得到的分子链所受到的拉力.

1.2 模拟结果

PVDF单分子链在链长方向的单轴拉伸下，能量变化如图2所示.曲线显示系统能量随着应变的增加而不断升高，当应变εzz=0.3时，能量发生突降，预示着分子链可能已经被拉断，这一点可被随后的构型观测结果证实.考虑到本文的模型选取较小，能量的变化主要来自于C-C之间的共价键作用，由此得到β相PVDF分子链上C-C键的临界断裂能约为1.14 eV.

图2 单轴拉伸载荷下PVDF单分子链的能量变化曲线

当εzz=0.29(第29步)和εzz=0.3(第30步)时，分子链构型如图3(a)、(b)所示.显然，拉伸到第30步时(εzz=0.3)分子结构已经失稳，相邻两个C-C键显示巨大差异(键长分别为0.292 6、0.131 3 nm)，键长较大的那个 C-C 键(0.292 6 nm)已经被拉断.而在此前的拉伸步(εzz=0.29)，两个相邻C-C键的键长接近(分别为0.182 5 、0.183 2 nm)，表明此时结构稳定，分子链并未被拉断.构型分析的结果与前面能量变化的结果一致，当应变εzz达到0.3时，β相PVDF单分子链将会被拉断，其临界断裂键长可认为是0.183 2 nm.

图3 PVDF分子链在应变为0.29和0.30时的构型

拉伸过程中PVDF单分子链的拉力变化曲线和C-C键的弹性系数k的变化分别如图4(a)、(b)所示.拉力变化曲线显示，当能量发生突变时(εzz=0.3，对应分子链被拉断)拉力也发生了突变.而在εzz=0.28时拉力达最大值，为6.75 nN，显示PVDF单链在断裂之前有一个非常微小的屈服过程，屈服拉力约为6.75 nN.图4(b)的曲线显示只有当应变非常小时(εzz＜0.03，对应着弹性理论里面的小变形)，弹性系数k才可近似认为是常数，而当应变增大时，k迅速降低，在断裂之前趋于零.普通的势函数描述C-C之间的弹性能时并未考虑这一因素，而是处理成常数，由此可见普通的势函数(分子力场)只适用于弹性范围内的小变形模拟，对于大变形的模拟则必须考虑这一因素.

图4 拉伸过程中PVDF单分子链的拉力变化和C-C键的弹性系数变化曲线

1.3 讨论

由于第一性原理模拟的计算量大，本文的模拟所采用的模型较小，模拟所选取的精度也不是非常高.不过考虑到研究β相分子链拉伸过程中的断裂性能并用于分子动力学的模拟中，这样的计算精度是可以接受的.对单链的模拟结果表明，β相的PVDF分子链被拉伸时最终会发生链的断裂，而在这个过程中，系统能量和C-C共价键弹性系数k的变化都是非线性的.因为高分子材料的分子动力学模拟所用的分子力场如COMPASS、PCFF、MM3、UFF等都只能模拟小的弹性变形，无法描述大变形甚至断裂过程，本文所得到的这些拉伸结果及影响参数都是在利用分子动力学模拟来研究块体材料大变形拉伸性质时需要考虑并注意的问题.

2 无定形PVDF的元胞模型的构建

在以往研究中给出了高分子材料的单分子链模型和包含多个无定形分子链的元胞模型的构建方法［11－12］.为了使构建的模型更加与实际的分子链构型相符，本文采用了一种较为新颖的方法来进行.该方法主要采用Materials Studio软件中的Perl script脚本语言完成.其基本步骤及实现过程为:1)将80个聚合度为200的PVDF分子链按照1.60 g·cm－3的密度装入长方体元胞中建立具有周期性边界条件的初始长方体模型(其长宽高a、b、c分别为 a=9.335 nm、b=9.335 5 nm、c=12 nm).在这一过程中，分别将每个分子链端部的C原子设置成“活性”原子，意味着在随后的过程中，这些较短的链之间有机会相互连接而形成长链.2)在NVT系统下对第1步建立的初始模型进行分子动力学模拟.在模拟过程中，每个“活性”C原子不断地在其截断半径内(1 nm)寻找其他分子链上的“活性”原子，一旦发现有满足条件的原子，两个“活性”原子发生相互连接，同时删掉“活性”标记，变为“稳定”原子，这一过程贯穿整个动力学模拟过程，直到该“活性”原子完成连接或自始至终也无法找到可连接的原子.为了避免同一链上的两个“活性”原子相互连接而形成环链，此处又对每个原子加了一个额外标记用以记录分子链信息，那些只有来自不同的分子链上的“活性”原子才能发生反应，并且一旦发生反应后，新形成的链上的其他“活性”原子重新进行标记用以表明它们位于同一分子链上.3)检查生成的含有不同长度的分子链的模型的完整性，对不饱和C原子进行加H处理并对模型进行彻底的NVT模拟进行能量优化和弛豫后，得到的模型作为单轴拉伸的分子动力学模拟模型.本文利用这种方法最终建立的模型包含的分子链数目为30(一共是96 060个原子)，最大的聚合度为2 200(意味着该分子链由11个初始分子链连接而成)，最小的200(该分子链未和其他分子链相互连接)，更多的分子链聚合度为3～8(整体符合正态分布).这样的模型与实际情况更加吻合，因为在实际材料中，PVDF的链长是不尽相同的，从聚合度为数万到数千的都可能存在，此处建立的模型在一定程度上可以描述分子链间的这种差异.图5显示了初始模型和经过上述方法处理过的模型，可明显见到分子链数目减少且变得更长.

图5 “活性”原子发生反应前和发生反应后的无定形PVDF的元胞模型

3 无定形PVDF单轴拉伸的分子动力学模拟

利用无定形PVDF的元胞模型构建的模拟模型，采用LAMMPS软件和PCFF分子力场对无定形PVDF的单轴拉伸进行了分子动力学模拟，模拟的温度为300 K，最大拉伸应变设为1，常应变率模拟.模拟中模型沿着z轴(c方向)被拉伸，两端各占模型长度10%的部分被固定，3个方向都采取周期性边界条件.模拟分两种情况:第1种在拉伸过程不考虑分子链断裂(传统模拟)，模拟采用的应变率分别为1010、109和108s－1;第2种情况考虑了链的断裂情况，其应变率分别为109s－1，5×108s－1和准静态(模拟应变率为零).后一种情况下，依据PVDF单链拉伸的第一性原理模拟结果，设置C-C键的临界断裂键长为0.184 nm.实际的模拟中，依据模拟结果，实际得到的最大应变从0.22～1.00不等.因为较低的应变率下模拟需要耗费非常长的时间，本文只进行了最低应变率为108s－1的动态模拟.

图6显示了不同应变率下不同情况的模拟中所得到的应力－应变曲线.当不考虑C-C键的断裂时，应变为 1010，109，108s－1的应力 － 应变曲线非常相似，应力能够随着应变的增加而达到非常高的值，而这与实际结果是不相吻合的.事实上在如此高的应变率下，PVDF很容易被拉断.只有在准静态拉伸下才可能达到非常高的应变.而当考虑C-C键的断裂因素时，应变率为109s－1的应力－应变曲线显示当应变达到较低的0.21时，模型即被拉断，应力急剧下降，这表明模型实际是在冲击作用下发生的断裂，109s－1的应变率实际对应着很高的冲击速度;当应变率降为5×108s－1时，实际结果相差不大，发生断裂时的应变仅仅有少许的增加.这些结果表明动态模拟PVDF的拉伸过程，因为计算速度的限制，所选择的应变率太高而难以模拟实际拉伸过程.为此本文采取了准静态模拟，并考虑键的断裂因素.从准静态所对应的应力－应变曲线上，可以发现当应变较小时，其曲线形状与其他情况下得到的结果是一致的，但是当应变在0.15附近时，应力有一个小的突降，但整个曲线走势变化不大，表明在这个应变下模型内部可能发生了某个分子链的断裂，但是直到应变达到0.85时，从应力的变化上看，整个模型都未发生断裂.这样的结果与实际情况是相吻合的.由此可见，考虑键断裂因素的准静态模拟是模拟大变形的有效途径.

进一步地，当不考虑键的断裂时应变率为109s－1下和考虑键断裂时准静态下的部分模型构型如图7所示.可明显看出，当应变率为109s－1，ε=0.34时，模型中出现明显的微孔洞，而当ε=0.85时，模型中微孔洞发展成了比较大的空洞，而且银纹也已连接成大的微孔洞，模型本质上已经破坏，但因为链不断裂，所以整个模型被拉成了两部分，中间由几个分子链保持连接.而当考虑键断裂因素，在准静态加载条件下，即使当ε=0.85，模型中也仅见少量微孔洞，而未见整体断裂迹象，但是通过对每个分子链的构型进行分析，确实发现了有3个分子链已被拉断，这也证明了上述应力－应变分析的结果.由此可见在大变形的拉伸过程中，PVDF分子链确实会发生死缠接现象而被拉断，但整体体现得并不明显，即使对每个分子链都采取了较低的断裂准则(键长为0.184 nm)，整体模型在较大应变下仍未发现明显的损伤破坏.

图6 无定形PVDF在考虑键断裂和不考虑键断裂时在不同应变率下的应力－应变响应

图7 PVDF的构型变化

4 结论

1)β相PVDF分子链在链长方向被持续拉伸会发生断裂，断裂时有一个微小的屈服过程，分子链的弹性系数只有当应变很小时才可以近似看作是常数.

2)分子链被拉断时存在临界拉力、临界键长及临界断裂能，这些参数可作为大变形模拟时分子链的断裂准则.

3)构造大规模无定形材料模拟模型的方法对于构造与实际构型相吻合的无定形材料的模拟模型具有很好的实用性.

4)应变率是在拉伸分子动力学模拟时必须考虑的重要因素，并且只有当考虑键的断裂因素时才能够较真实地模拟材料的大变形行为，PVDF在大变形的拉伸下，存在少量分子链发生死缠接而被拉断的现象，但整体体现得并不明显.

［1］YEE W A，KOTAKI M，LIU Ye，et al.Morphology，polymorphism behavior and molecular orientation of electrospun poly(vinylidene fluoride)fibers［J］.Polymer，2007，48(2):512－521.

［2］WANG Mian，SHI Jiahua，PRAMODA K P，et al.Microstructure，crystallization and dynamic mechanical behaviour of poly(vinylidene fluoride)composites containing poly(methyl methacrylate)-grafted multiwalled carbon nanotubes［J］.Nanotechnology，2007，18(23)，235701.

［3］黄维恒，闻建勋.铁电压电高分子材料［M］//闻建勋，王贤姗.高技术有机高分子材料进展.北京:化学工业出版社，1994:526－553.

［4］陈贻瑞，王健.基础材料与新材料［M］.天津:天津大学出版社，1994.

［5］TUZZOLINO A J，ECONOMOU T E，CLARK B C，et al.Dust measurements in the coma of comet 81P/Wild 2 by the dust flux monitor instrument［J］.Science，2004，304(5678):1776－1780.

［6］McDONNELL J A M，BURCHELL M J，GREEN S F，et al.The stardust dust flux monitor［J］.Advances in Space Research，2000，25(2):335－338.

［7］KARASAWA N，GODDARD W A.Force-fields，structures，and properties of poly(vinylidene fluoride)crystals［J］.Macromolecules，1992，25(26):7268 －7281.

［8］曾凡林，孙毅.PVDF单链拉伸的第一性原理模拟［C］//第七届海峡两岸工程力学研讨会论文摘要集.哈尔滨:中国力学学会，2011:111－116.

［9］CEPERLEY D M，ALDER B J.Ground-state of the electron-gas by a stochastic method［J］.Physical Review Letters，1980，45(7):566 －569.

［10］FUJINAMI Akinori，OGATA Shigenobu，SHIBUTANI Yoji.Ab initio study of the tensile behavior of single polyimidemolecularchain ［J］. Polymer， 2004，45(26):9023－9028.

［11］曾凡林，孙毅，周玉.有机材料的分子模拟模型(I):单分子链的构建［J］.哈尔滨工业大学学报，2009，41(10):95－99.

［12］曾凡林，孙毅，周玉.有机材料的分子模拟模型(II):模拟元胞的构建［J］.哈尔滨工业大学学报，2010，42(2):281－287.

Ab initio and molecular dynamics study of uniaxial tensile behaviors of PVDF

ZENG Fan-lin，SUN Yi，LIU Yi-zhi

(Dept.of Astronautic Science and Mechanics，Harbin Institute of Technology，150001 Harbin，China)

To obtain the behaviors of Poly(vinylidene difluoride)(PVDF)under large tensile deformations，this paper investigates the tensile behavior of a single molecular chain of rod-like PVDF with β phase structure at zero temperature using a Density Functional Theory(DFT)calculation.We obtain the critical bond length，critical fracture energy，critical tensile force and the elastic constant when the β phase PVDF chain is stretched to break.Then a novel method is proposed to build the cell model of the amorphous polymer，which is used for the molecular dynamics(MD)simulations of uniaxial tensile of PVDF.During the simulations，two cases that the chains will never break and are permitted to break are considered.In the first case，the simulations with three different strain rates 1010s－1，109s－1，and 108s－1were carried out and the strain rate dependence of the stress-strain relationship，the structure deformation were analyzed as well.In the second case，the simulations with a strain rate of 109s－1，5 ×108s－1and a quasistatic strain loading were performed and the stressstrain relationship，the cell structures before and after the breaking were researched.Some key factors to influence the simulation results and some important deformation information were analyzed.This is helpful to elucidate the inherent deformation mechanisms during the tensile procedure.

PVDF;uniaxial tensile;DFT;molecular dynamics;fracture

O631

A

0367－6234(2012)09－0046－05

2011－08－26.

国家自然科学基金资助项目(11102053);哈尔滨工业大学科研创新基金资助项目(HIT.NSRIF.2010070);哈尔滨市科技创新人才研究专项资金资助项目(2012RFQXG001).

曾凡林(1978—)，男，副教授，硕士生导师;

孙 毅(1961—)，男，教授，博士生导师.

曾凡林，zengfanlin@hit.edu.cn.

(编辑 张 红)