环状聚合物及其对应的线性链熔体在启动剪切场下流变特性的分子动力学模拟研究*

杨俊升 黄多辉

1)(中国科学技术大学,国家同步辐射实验室,合肥 230026)

2)(宜宾学院,计算物理四川省高等学校重点实验室,宜宾 644007)

1 引 言

环状聚合物是自然界比较常见的分子形态,如质粒、基因组、肌动蛋白和多糖等分子结构[1−3].关于环状聚合物体系的研究一直以来也是高分子物理研究中的热点和难点之一[4−7].主要原因是环状分子链与其相同分子量的线性链比较没有链端,进而在其加工和实验研究中就会表现出不同的应力响应和流变特性[8−10].目前实验对于这一特殊的过冲响应的分子机理还没有做出更清晰的分子机理解释,实验上也很难从分子层面观察链段的演化信息,所以急需在分子层面对于这一应力过冲现象给出一个定性和定量的解释.

目前,关于非串联的环状聚合物熔体的研究主要集中在其构象演变和流变性质.关于环状聚合物构象的动力学研究,有些人认为环状链是一个双折叠的线性链[5],另一种观点认为环状聚合物在熔体中是网状折叠分布的[11,12],而且应力是通过环与环之间的轮廓滑移松弛的.模型不一样,所带来的对应力松弛的解释也就不同,目前这个问题仍然没有得到解决.关于非串联的环状链熔体的流变学方面: 实验结果显示环状链熔体的零剪切黏度明显要比其线性链对应的黏度小,而且环状聚合物熔体系对应的最大过冲应变与剪切速率之间并不满足其线性链在高剪切速率下的1/3幂指数的标度关系,而且其对应的线性前驱体在剪切流场下具有明显的应变过冲现象.最近的实验研究还发现环状分子熔体在剪切场下表现出了更弱的剪切变稀行为[8,9,13,14].故而,环状聚合物结构和流变学研究中表现出来的差异性表明针对环状聚合物熔体还需深入进行分子机理的研究.

目前关于环状分子体系的研究,实验很难测得高分子加工过程中分子链段的信息变化.分子动力学(moleculardynamics,MD)和非平衡分子动力学(non-equilibrium molecular dynamics,NEMD)模拟由于其精确的原子模型和小的时间尺度跟踪能力,可以作为解决上述拉伸和剪切流变实验问题的有效工具[15,16],从而为环状聚合物熔体在不同流场强度下的分子过程提供有价值的分子图像和理论解释.本文主要以环状聚合物及其对应的线性链模型体系为研究对象,通过非平衡分子动力学模拟方法探索其在剪切流场下所对应的结构和流变特性.

本文通过非平衡分子动力学方法对线性链和环状聚合物体系在不同剪切流场下应力响应的分子机理做了系统的研究.虽然理想情况下需要一个全原子化的模型,但理解缠结聚合物熔体的宏观动力学所需的体系的大小也是相当具有挑战性的一个问题,特别是考虑到聚合物运动的时间和长度尺度时.因此,本文采用粗粒化模型(coarse-grained,CG)方法来获取环状链体系应力过冲对应的链拉伸和取向分子信息.动力学模拟结果不仅仅是再现了实验观察到的环状链体系的流变特性,通过分子层面的观察与分析,给出了环状聚合物在不同剪切流场强度下应力响应所对应的分子结构与取向分子机理.

2 计算方法与模型

2.1 计算方法

本文采用分子动力学方法对缠结的聚合物熔体进行模拟,模拟软件和平台主要是在宜宾学院计算物理实验室自行搭建集群上搭建的LAMMPS平台上进行的.对于所研究的体系,粒子之间的非键相互作用采用Lennard-Jones(LJ)势来描述[16]:

其中,rij为第i个珠子到第j个粒子或单体之间的距离,s是粒子或单体半径.两个粒子相互作用的截止距离为21/6s,对应于粒子之间的纯斥力.在当前聚合物流变学研究中,吸引相互作用不应发挥重要作用,在此不予考虑.粒子之间的键合势采用有限伸展非线性弹簧模型(FENE)来描述:

其中,弹性常数是kb=30ε/σ2,最大键长为R0=1.5σ.选择较大的弹簧常数值,以减小网络变形时键的拉伸效应.在模拟中,链上两个珠子之间的实际平均键长是0.97s.

在MD模拟中,压力张量P也称为应力张量σαβ,或者仅仅是压力.对角线单元称为拉应力,非对角线单元称为剪应力.因为压强P是一个热力学量应力,σαβ是一个力学量,前者定义为时间平均值,后者可即时计算.它们都是张量,满足下面的关系:

在微观系统中,压力是根据维里方程作用在所有原子上的力之和计算出来的.对于加力场,这个方程通常写成:

其中,N为计算域中存在的珠数,kB为玻尔兹曼常数,V和T为计算域中的体积和温度,fij为珠i和j之间的力.

2.2 模型构建

为了研究线性链与它所对应的首尾相接的环状链体系所的流变特性,本文构建了长度N=400的粗粒化线性链和环状链分子体系,体系中包含了100根链,如图1所示.为了比较方便,在数据分析过程中将它们分别标记为linear和ring.环状链结构在熔体中非链接或者非串联的.根据Likhtman 和Larson[17]的工作,知道分子链结构体系对应的缠结长度为Ne=50,并且计算出线性链和环状链对应的松弛时间分别为1.45 × 105t和0.38 × 105t.初始结构建完之后,体系分别在T=1.0下松弛1 × 106t.两种体系的初始密度都为0.85σ−3.当体系经历了充分的松弛之后,再给体系施加剪切流场来计算其对应的流变特性.

图1 线性链与其对应的环状链对应的结构Fig.1.The structure of linear and ring chains.

3 计算结果与讨论

3.1 线性链及其环状链熔体在不同流场强度下的应力及黏性响应

图2(a)和图2(c)示出了线性链与其同分子量的环状链体系对应在不同流场强度下的应力-应变曲线.对于线性链,随着应变的增加,剪切应力演化在各个应变速率下都存在一个过冲应变点,而且随着剪切速率的增加过冲应变也在不断地增加.应力过冲和软化之后,线性链熔体对应的应力响应会进入一个平台区,也就是稳态.但与其对应的环状链体系比较,会发现当应变速率小于1 × 10–4t–1时,剪切应力没有出现明显的过冲现象,应力屈服之后直接进入的是一个平台区.当剪切速率大于1 × 10–4t–1时,才出现了微弱的过冲现象.大应变速率下的应力响应与线性链的响应一致.继续分析线性链和环状链最大应力的大小,不难发现线性链最大应力比其环状链对应的最大应力大很多.本文模拟结果与实验所观察到的结果也是一致的[8],证明我们选取的模型及模拟过程的合理性.同时,这个差异性也引起了我们的关注,但是实验上关于这种现象的分子机理很难检测到,也没有给出明确的分子解释.所以,很有必要在分子层面研究线性链与其对应的环状链在不同剪切应变速率下的分子链段的拉伸和取向.黏性和剪切应力满足以下公式[18,19]:

图2 (a),(b)线性链体系对应的应力随着应变及剪切黏性随着时间的演化过程；(c),(d)环状链体系对应的应力随着应变及剪切黏性随着时间的演化过程Fig.2.(a),(b)Stress-strain and nonlinear startup shear viscosity as function of strain and time for linear polymers,respectively;(c),(d)stress-strain and nonlinear startup shear viscosity as function of strain and time for ring polymers,respectively.

其中,σxy是剪切应力,是剪切速率.根据剪切应力和黏性的关系,图2(b)和图2(d)给出了线性链和环状链的黏性随着时间的演化过程.稳态峰值后的瞬时黏度降低反映了链的收缩(表现为部分解缠结或翻滚).随着剪切速率越小,在稳态时对应的黏性越大；随着剪切速率增大,黏性也会随着减小,表现出了剪切变稀行为.

3.2 线性链及其环状链熔体在不同流场强度下的流变特性

为了进一步给出线性和环状链体系的流变特性的差异,从图2中分别提取了线性和环状体系在过冲点和稳态所对应的流变特性随着流场强度的变化.为了分析和比较方便,这里将流场的强度转换成维森伯格数(WiR=R),其中τR是线性链和环状链对应的Rouse松弛时间.线性链的松弛时间是其对应的环状链松弛时间的4倍(τR.ring=τR.linear/4)[8,20].根据Likhtman 和Larson[17]的工作,可计算出线性链和环状链对应的松弛时间分别为1.45 × 105t和0.38 × 105t.根据松弛时间τR和剪切速率,就能给出相应的剪切速率下其对应的维森伯格数WiR.

图3(a)示出了稳态黏性随着时间的变化.随着WiR的增加,线性和环状链体系都表现出了剪切变稀的行为,但是环状链体系所表现出来的剪切变稀行为明显要比线性链的更弱.线性链和环状链体系所应的黏性变化与WiR之间所满足的标度关系与实验测量出来的标度关系也是吻合的,这也从侧面证实我们构建的模型是合理的[8,13],可以进一步探索其流变特向所对应的结构变化.

图3(b)示出了体系在剪切过程中最大黏性和稳态黏性的比值(ηmax/ηsteady)随着WiR的变化,其反映的是稳态有效链的形变量.事实上,当WiR<1时,对于线性链和环状链,ηmax/ηsteady都是趋近于1的,在这个区间链实际上有足够的时间松弛.当WiR>1 时,随着剪切速率的增加ηmax/ηsteady也随着增加,并且线性链和环状链体系与时间的剪切强度表现出了不同的标度关系,0.08明显也要比线性链的0.13小.进一步证实了稳态流场下线性链的线团要比环状链的线团容易发生形变.

图3(c)示出了最大应力所对应的峰应变(γmax)随着WiR的变化过程.对于线性链,应力过冲响应中所展现的标度关系在早期的流变学实验中也有发现,而且这些标度关系目前也被其他理论模拟证实.当WiR<1 时,应力过冲主要是由于聚合物链段的取向过冲所引起的,而且最大过冲应变主要处于应变为1.8—2.3之间[8,14,21,22].当WiR>1 时,应力过冲主要是由聚合物链段的取向和拉伸的协同作用决定的,而且γmax与WiR之间满足指数为1/3的标度关系.对于环状链体系,分析结果显示γmax与WiR之间满足指数为0.1,也明显小于其线性链对应的标度指数.图3(d)给出了线性和环状链体系所对应的最大应力σmax与WiR之间的关系,我们发现相同的WiR下最大应力与体系的构象没有关系,随着剪切速率增加线性链和环状链体系都满足指数为0.5的标度关系.

从流变特性的变化,可以认为双折环的构象是其变形能力较弱的主要原因.换言之,与其对应的线性链相比,环型聚合物由于结构紧凑,在剪切流作用下改变构象的自由度较小.此外,环状链的链段在剪切流作用下的协同运动与线性链段的运动不同.这是都是实验所无法观察到的,需要通过分子模拟进一步验证.

3.3 稳态线性及其环状分子整链构象与取向在不同流场强度下的演化

为了进一步探索线性链及其对应的环状链达到稳态以后的分子构象和取向的变化,本文计算了回转半径(Rg)和取向角(tan2q)的随着流场强度的变化,如图4所示分别表示整链的均方回转半径和初始结构的均方回转半径)代表分子链的形变程度.从图4(a)中不难发现,小WiR下线性链体系更容易形变,但是当WiR>20时,在相同的WiR下,两者的形变也基本相同,并且都满足的幂指数关系.

图4(b)给出了线性链和环状链体系与流场方向取向角的变化.环状分子和线性分子的取向采用以下的方程来描述[13]:

其中,Gxy,Gxx和Gxy分别代表分子链构象在不同方向上的回转张量平方.从图4(b)中可以看到线性链和环状链的取向随着WiR的增加而减小,但是线性链明显要比环状分子更容易发生取向,并且线性链和环状链的平均取向角与流场强度WiR分别满足的标度关系.

图3 不同流场强度下线性链及其环状链的归一化稳态黏性 ηsteady/η0(a),最大黏性和稳态黏性比值 ηmax/ηsteady(b),峰应变γmax(c)和最大应力 σmax(d)随着 WiR 的变化Fig.3.(a)Evolution of steady viscosity normalized with the zero-shear viscosity ηsteady/η0 ,(b)maximum viscosity scaled with steady viscosity ηmax/ηsteady ,(c)peak strains γmax and(d)peak shear stress σmax as a function of WiR.

3.4 线性及其环状分子链段构象和取向分布的演化过程

为了进一步探究不同剪切流场下,环状链与其对应的线性链分子链段信息的变化,将线性和环状分子按缠结长度Ne=50分成小段来考察链段的分布和拉伸比分别对应的是开始剪切之前链段的平均末端距和不同应变下链段的平均末端距)及链段沿着流场方向取向角分布及平均值链段矢量与流场方向的平均夹角)的变化.

图5示出了不同WiR下链段的长度分布概率及链段平均拉伸比随着应变的变化过程.在弱流场下,线性链与其环状链对应的链段的分布基本上是高斯分布,平均拉伸比保持在1附近.但是当WiR增加到10以后,会发现线性链的链段的平均拉伸比会明显增加,而且到达最大值后会有一个明显的拉伸过冲现象,暗含链段有一个明显的回弹.相比于线性链,环状链体系的链段拉伸比在强流场下没有出现明显的过冲现象,这与目前看到的应力没有明显的过冲现象也是一致的,从而可以推断环状链体系之所以没有明显的过冲现象是由于环状链自身弱的拉伸比及剪切过程中链段没有明显的回弹导致的.最近Kremer等[1]的结论也证实了环状高分子链在小应变速率下构象不发生明显的拉伸,从侧面也验证了我们的结论[6].不仅仅是环状高分子链分子存在这样的现象,环状的DNA分子与其线性链比较的实验也发现了类似的结果.图6给出了链段取向角的分布变化.对于线性链和环状分子,在低剪切速率下,链段的取向都没有发生明显的变化,而且θ分布都很宽,但是随着流场强度的增加,θ的分布也明显变窄,表明分子链段存在明显的取向.当体系达到稳态之后,线性链体系和环状分子体系θ并没有减小.从而可以推断: 线性链体系表现出明显的过冲主要是由于链段的拉伸回弹导致的.环状链体系之所以没有表现出明显的过冲现象是由于其分子的链段到达稳态之后并没有出现明显的回弹.

图4 (a)线性链和环状链体系的 ＜＞/＜＞ 随着 WiR 的变化；(b)线性链和环状链熔体的取向 ＜tan2θ＞ 随着 WiR 的变化Fig.4.The ＜＞/＜＞(a)and tan2θ(b)as a function of WiR for linear and ring polymers.

图5 不同 WiR 下线性与环状分子链段长度分布随着应变的演化Fig.5.Evolution of segmental length distribution of linear and ring polymers under the different WiR.

4 结 论

本文通过非平衡分子动力学方法研究了环状分子及其对应的线性链分子体系在启动剪切流场下的结构和流变特性.动力力学模拟结果显示: 线性链与其对应的环状分子比较在低维森伯格数下具有更明显的应变过冲现象.这主要是由于线性分子链段在剪切场下发生了明显的回弹、但是环状链体系并没有出现明显的回弹所引起的.本文还提取了线性链和环状链体系在过冲点和稳态的流变数据(黏性、过冲应变、过冲应力及最大黏性与稳态黏性的比值),并给出了线性链和环状链体系流变特性分别满足的标度关系,为进一步的实验研究提供了理论参考.

图6 不同 WiR 下线性与环状分子链段沿着剪切方向角分布随着应变的演化Fig.6.Evolution of segmental angle distribution for linear and ring polymers under the different WiR.