聚丙烯／聚苯乙烯共混体系相态结构的耗散粒子动力学模拟

陆未谷，张玉梅，陈彦模

（东华大学a.材料科学与工程学院；b.纤维材料改性国家重点实验室，上海 201620）

耗散粒子动力学（dissipative particle dynamics，DPD）模拟是对于具有动态和流变性质的简单及复杂的流体的一种计算模拟方法.在DPD模拟中构建一系列粒子（即所谓粗粒化），其相互间的碰撞规则和运动规律遵循牛顿运动定律，这些粒子不是代表单一原子，而是代表整个分子或流体的区域，且计算过程与原子的细节无关.对于小分子而言，粒子构建较为容易，一般为整个分子；对于高聚物，情况变得复杂，DPD模拟采用“珠－簧”模型[1].近年来，DPD模拟已被应用于物理和化学研究的众多领域.

在聚丙烯（PP）纤维实际生产与研究中，其致密的超分子结构成为一个棘手的问题，如纤维可染性和细旦化的集合共混成型加工.因此研究共混改性PP，特别是共混体系相态结构，具有重要的意义.

本文采用DPD模拟对PP／聚苯乙烯（PS）二元共混体系进行系统研究，将温度、PS体积分数、PS相对分子质量以及剪切速率作为考察因素进行分析，为共混体系的相态结构控制提供较为全面的依据和参考.

1 DPD模拟基本原理

DPD模拟中，粒子间有3种成对出现的作用力：保守力（conservative force）、耗散力（dissipative force）、随机力（random force）[2－3].耗散力是指当前模拟尺度下体系消耗掉进入更微观尺度下的能量，随机力则指从更微观尺度进入当前尺度下的能量.因此，用牛顿运动定律计算粒子状态，必须通过引入耗散理论，设置耗散力参数和随机力经验值以保证模拟体系的动量守恒，即温度恒定.文献[4]将DPD模拟与Flory－Huggins平均场理论相结合，建立高聚物体系的原子和分子信息与DPD粒子之间的映射关系.

就相对分子质量确定的高聚物分子而言，分子尺寸会随着构象的变化而变化[5－6].文献[7]认为，应采用高分子链的持久长度所涉及的范围来构建高聚物分子的粒子.因为长度长于持久长度的高聚物，其性质可以用统计学方式描述，就像三维随机游动，可以不考虑构象的影响.由于DPD模拟并不涉及分子细节，因此在这种尺度上构建粒子是合理的.

在DPD模拟中，为了简化单位，通常采用模拟构建的单个粒子的质量和截断半径rc（粒子间相互作用范围）来作为单位量纲，本文也采取这一做法.

2 DPD模拟参数计算

2.1 特征比

对真实的高聚物分子PP和PS进行粗粒化，关键步骤是建立真实分子结构和介观参数之间的联系.本文中用高分子链持久长度中所包含的重复单元数即特征比来建立这一联系，式（1）对真实高聚物分子链进行粗粒化，得到粗粒化模型[8].

其中：NDPD为高分子单链上DPD粒子的个数；Mp为高聚物的相对分子质量；Mm为单元结构的相对分子质量；Cn为特征比.

考虑到计算机的运算能力，选定相对分子质量比较适中的PP和PS样品，其相对分子质量分别为9 000和10 000，特征比分别为6.7和9.9[9].为了运算方便，选择的特征比为10.按照式（1），具体的粗粒化方法[10]以PP为例：相对分子质量为9 000，重复单元的相对分子质量为42，其数均聚合度N≈214，因此NDPD＝214／10≈21.通过此方法得到每条PP分子链上含有的DPD粗粒化粒子数目，之后再用弹簧将相邻DPD粒子链接起来，得到线性PP分子的DPD粗粒化模型.采用同样的方法对PS分子进行粗粒化.

2.2 保守力排斥强度选择

在DPD模拟中，最为重要的参数是DPD粒子之间的保守力排斥强度（同类型粒子间的保守力排斥强度为aii、不同粒子间的保守力排斥强度为aij，其中aii与模拟体系的数密度和Flory-Huggins相互作用参数有关）.为了有效关联DPD粒子体系的压缩系数与高分子熔体的压缩系数，且能合理地描述流体的密度变化，aii的选取应遵循式（2）[2]：

其中：kBT反映出温度信息（kB为波尔兹曼常数，T为绝对温度）；ρ为数密度.在模拟中，选取kBT＝1，ρ＝3，计算可得同类型粒子间的保守力排斥强度aii＝25.不同粒子间的保守力排斥强度aij可使用Flory－Huggins相互作用参数χ计算.

对于PP／PS的二元共混体系，通过 Material Studio软件的Blend模块进行计算，拟合直线如图1所示.

图1 PP／PS共混体系的Flory-Huggins参数χ与 温度倒数1／T之间的关系Fig.1 Interaction parameterχas a function of inverse temperature for a series of PP／PS blending system

由图1可得，拟合直线为

依据式（3）可求得450K时，χ为0.541.但需要注意，通过式（3）求得的是丙烯单体和苯乙烯单体之间的Flory-Huggins相互作用参数，为了得到粗粒化之后的DPD粒子与粒子之间的相互作用参数，还需要在单体间相互作用参数的基础上乘以特征比Cn.如前所述，选取的特征比为10，因此，DPD粒子与粒子之间的Flory-Huggins相互作用参数应为5.41.根据式（4）可以求得不同种类DPD粒子之间的保守力排斥强度[4]：

不同温度下的PP和PS的保守力排斥强度如表1所示.

表1 PP和PS分子参数Table 1 PP and PS molecular parameters

2.3 其他模拟参数的选择

本 文 采 用 正 则 系 综 （canonical ensemble，NVT，表示具有确定的粒子数N、体积V、温度T）作为宏观约束条件，并应用周期性边界条件进行模拟.计算盒子的大小为20×20×20r3c，体系数密度为3，因此体系包含DPD粒子24 000个.为了便于模拟，截断半径rc、质量m以及温度kBT的单位都设置为l，弹簧常数C为4.0，耗散因子、噪声强度均选择DPD方法的经典值（4.5和3）.由于所采用的积分方法为velocity-verlet算法，为了使在该积分方法中能较好地控制温度，所选用的时间步长为0.05.每一次模拟均从均匀共混状态开始，运行至少1×104步直至相态结构达到平衡.

3 PP／PS相态结构的形成

图2所示为PS相在共混体系中的聚并过程，通过与文献[11]中的扫描电镜图样（如图3所示）的对比，发现随着动力学运行时间的增加分散相颗粒的体积都会因聚并而逐渐变大，但是，由于模拟中看到的分散相情况与扫描电镜照片[11]所反映出的现实中共混体系的结构在空间尺度上的差异，模拟的有效性和真实性需要验证.

图2 PS体积分数为80%，kBT＝1，每1 000步模拟截取的图形Fig.2 PS volume fraction is 80%，kBT＝1，the figure of each 1 000steps in simulation

图3 PP／PS（20／80）不同聚并时刻的电镜图片[11]Fig.3 PP／PS（20／80）electron microscopic images of different coalescence moments

在DPD模拟中看到的分散相情况与扫描电镜照片有共同点并非偶然.由于计算机性能的局限性，模拟的体系不能很大，而现实中的共混体系无论是平衡时间或是微相尺寸都要比模拟的体系大好几个数量级.另外计算机模拟时的初始状态两相共混得非常均匀，分子链之间的距离很短，随后开始的聚并过程就很迅速；而现实中的共混不可能完全均匀，PP和PS在相对较大的尺度中依然是两相分离的；PS颗粒之间的距离相对较远，热扰动的效果不明显，碰到一起就需要很长的时间，如果没有剪切诱导，完全静态地进行聚并则时间更长.

分子模拟过程提供了一个能在小尺度下直观考察相行为的机会.因为，无论是宏观聚并还是微观聚并，本质上是一样的，都是分子在一定的初始热扰动的基础上开始运动而产生的结果.

4 PP／PS相态结构的影响因素分析

4.1 模拟参数选择

静态聚并过程模拟（剪切速率为0）：kBT＝1.0，1.1，1.2，1.3；PP体积分数＝0.1，0.2，0.3，0.5；PS分子链粒子数NDPD＝5，10，20，50；其他模拟参数参照2.3节选取.

剪切聚并过程模拟（剪切速率为x轴方向上）：kBT＝1.2；剪切速率r′＝0.2，0.4，0.6；PS体积分数＝0.3；PS分子链粒子数NDPD＝10；其他模拟参数参照2.3节选取.

4.2 结果讨论

4.2.1 温度对共混体系相形态的影响

不同温度条件下的模拟结果如图4所示，图中的颗粒为PS分散相.由图4可以看出，温度越高，则PS相的分散性越好，颗粒尺寸越小.因为温度的升高加快了高分子及其链段运动速度，同时降低了两相组分的黏度，降低了体系内界面的表面张力.

图4 不同温度下PS分散相的最终结构Fig.4 The PS final structure at different temperatures

4.2.2 剪切速率对共混体系相行为的影响

图5所示为不同剪切速率下共混体系动力学模拟的情况.从图5中可以看出，随着剪切速率的增大，PS分散相的颗粒明显变小，同时沿剪切方向，分散相被拉长，说明剪切破坏了分散相的聚并过程.

图5 不同剪切速率下PS分散相的最终结构Fig.5 The PS final structure at different shear rates

另外，分析分散相的形成过程发现，随着剪切速率的增加，稳定的分散相颗粒形成的时间提前了.图6列出了图5所示模拟过程中前3帧过程的PS相态结构.从图6可以看出，初始结构相同，剪切速率为0.6的体系在3帧后基本形成最终结构，而没有剪切的体系，则很难看出最终结构情况，说明剪切诱导聚并的发生.

图6 不同剪切速率下PS分散相的演变过程Fig.6 The PS evolution process of the dispersed phase at different shear rates

在外界力场作用下，分相颗粒的分散过程和聚并过程共同决定了共混体系的相态结构.条件不同则这两个过程对最后形成的相态结构的影响程度也不同，对此可用Elemendorp相图（图7所示）[12]来研究.从图7可知，相图包含了颗粒破裂曲线和颗粒聚并曲线，在全图中仅在区域Ⅱ内，相态结构的演变仅由聚并过程控制.

在前述模拟过程中，不同剪切速率下颗粒的平均尺寸和粒间距离都变大，这一现象与静态的聚并过程相似；不同的是，最终稳定的分散相颗粒形成的速率明显大于静态聚并的情况，说明剪切作用的存在能够有效地加快颗粒间的聚并.

4.2.3 PS体积分数对共混体系相行为的影响

图7 Elemendorp相态图[12]Fig.7 Elemendorp phase schemes

通过观察PS相的聚并过程，分析PS体积分数对共混体系相态结构的影响，不同PS体积分数下共混体系的相态结构如图8所示.试验开始默认PS均匀分布在PP相中，随着时间的推移，粒子在运动过程中发生黏结和聚并，达到平衡时，形成相态结构.随着PS体积分数的增加，PS相的尺寸变大，且分布变宽，聚并的倾向明显，PS分散相的数目减少.由于体积分数增加，PS相之间接触的机会增加，由于PS与PP之间的相容性不好，所以都朝着各自组分的方向前进，最终均形成连续相.

图8 不同PS体积分数下（kBT＝1.3）共混体系进行模拟后最终的相态结构Fig.8 The PS final phase morphologies structure at different PS volume fraction（kBT＝1.3）

4.2.4 PS相对分子质量对共混体系相行为的影响

相对分子质量对共混体系的影响较为特殊.由于研究体系发生变化，特别是PS的分子结构发生变化，分子间的Flory-Huggins相互作用参数也发生了变化，这会导致保守力排斥强度aij发生变化，模拟的结果可能不是很精确，因此对图9所示的不同相对分子质量PS的共混体系最终相态结构的图像变化趋势描述可能存在着误差.

图9 不同相对分子质量PS的（kBT＝1.2）共混体系进行模拟后最终的相态结构Fig.9 The PS final phase morphologies structure at different PS molecular weight（kBT＝1.2）

从图9可以看到，相对分子质量较小的PS相形成的分散相尺寸较均匀，相对分子质量较大的PS相形成的分散相尺寸差异较大.这是由于在分散相聚并过程中，相对分子质量越大的PS分子整体上受到的力就越多，而且各个方向都有，使得其所受合力较小，所以运动速度较慢，因此相对分子质量大的PS相发生聚并的概率就小，但是由于越长的高分子链之间的缠结点就越多，相互结合越紧密，分子间的作用力越大，表现出黏度很大，这导致其与PP相的黏度差异变大，因此PS相一旦发生聚并就很难分开，因而相对分子质量大的PS相会形成尺寸差异较大的分散相颗粒.

5 结 语

本文依据DPD模拟的原理，通过Material Studio 5.0软件的Blend模块计算了不同温度下Flory－Huggins相互作用参数χ及DPD粒子之间的排斥参数.通过对PP／PS二元共混体系的相态模拟，研究了温度、剪切速率、PS体积分数以及PS相对分子质量对PP／PS共混体系的影响，得到下述结论.

（1）共混体系温度越高，分散相PS相的分散性越好、颗粒尺寸越小.

（2）剪切速率增大，分散相的颗粒明显变小，分散性变好，同时通过对聚并过程的观察，了解到剪切又能诱导聚并发生.

（3）随着PS体积分数的增加，分散相尺寸逐渐增加，最终成为连续相.随着PS相对分子质量的增加，PS的黏度增加，造成PP与PS的黏度比变大，形成尺寸差异较大的分散相颗粒.

参 考 文 献

[1]SCHLIJPER A G，HOOGERBRUGGE P J，MANKE C W.Computer-simulation of dilute polymer-solutions with the dissipative particle dynamics method[J].Journal of Rheology，1995，39（3）：567-579.

[2]ESPANOL P，WARREN P.Statistical-mechanics of dissipative particle dynamics [J]. Europhysics Letters，1995，30（4）：191-196.

[3]朱孑系.液体自扩散系数与温度的关系[J].化工学报，1959（2）：172-175.

[4]GROOT R D，WARREN P B.Dissipative particle dynamics：Bridging the gap between atomistic and mesoscopic simulation[J].Journal of Chemical Physics，1997，107（11）：4423-4435.

[5]MA P W，LIU W C，WOO C H，et al.Large-scale molecular dynamics simulation of magnetic properties of amorphous iron under pressure[J].Journal of Applied Physics，2007，101（7）：073908-9.

[6]RUBINSTEIN M.Polymer physics of the perciliary layer[J].Pediatric Pulmonology，2007，30（sl）：112-113.

[7]GROOT R D，MADDEN T J.Dynamic simulation of diblock copolymer microphase separation [J].Journal of Chemical Physics，1998，108（20）：8713-8724.

[8]GAI J G，LI H L，SCHRAUWEN C，et al.Dissipative particle dynamics study on the phase morphologies of the ultrahigh molecular weight polyethylene／polypropylene／poly（ethylene glycol）blends[J].Polymer，2009，50（1）：336-346.

[9]BRANDUP J，IMMERGUT E H，GRUIKE E A.Polymer handbook[M].4th ed.New York：John Wiley and Sons，1999.

[10]苑世领，刘成卜，徐桂英，等.聚合物PVP与表面活性剂AOT相互作用的介观模拟[J].高等学校化学学报，2003，24（6）：1048-1051.

[11]李云岩.聚丙烯／聚苯乙烯共混物相结构形成演变及其动力学研究[D].天津：天津大学材料科学与工程学院，2007.

[12]ELMENDORP J J，VANDERVEGT A K.A study on polymer blending microrheology 4：The influence of coalescence on blend morphology origination [J].Polymer Engineering and Science，1986，26（19）：1332-1338.