（材料加工工程专业论文）聚合物流体微观和介观尺度下的模拟研究.pdf

摘要 摘要 粘度是聚合物流体的一个重要表征参数，也是聚合物加工成型性能的重要 影响因素。在进行传统的聚合物流动有限元模拟过程中，获得揭示材料物性的 本构方程( 粘度函数) 已经成为了模拟的关键。通过微观模拟或介观模拟的方 法来获得一定条件下的应力应变关系，从而避免封闭本构方程是实现复杂流体 模拟的重要方向。聚合物流体由于分子链的长链结构，尺度跨度比较大，单尺 度模拟往往难以实现。 在本文中，我们采用了将微观的分子动力学模拟和介观的耗散粒子动力学 模拟相结合的方法进行了聚乙烯熔体的非平衡态模拟，考察了在我们建立的模 拟体系中的粘度随剪切速率的变化情况。本文的主要工作如下： 1 )建立了4 5 0 k 、l a t m 下4 0 条链长为9 6 的聚乙烯体系的分子模型，并进行了 分子动力学的模拟。通过径向分布函数、均方回转半径、密度、以及均方 位移来对聚乙烯模型进行表征，获得平衡态下该模型结构和动力学信息。 经过与实验结果的对比分析，发现在所给的条件下，所建立的分子模型及 模拟方法是合理的。 2 )采用基于径向分布函数的方法建立了耗散粒子动力学模拟保守势能，实现 了从分子动力学模拟到耗散粒子动力学模拟的贯通连接。考察了几个关键 模拟参数对模拟结果的影响，将d p d 模拟结果与分子动力学模拟结果进行 对比，选择了合适的参数作为进一步模拟的参数。 3 )建立了平板拖拽流和反转泊萧叶流两种不同的模型对聚合物流体的流动进 行了模拟，通过计算得到了流动过程中的速度分布、剪切应力的分布以及 剪切粘度，并将所得的结果与理论分析的结果进行了对比分析，发现两种 模型的速度分布以及剪切应力的分布都具有与理论分析的结果相同的变化 趋势。 关键词：聚合物流体；分子动力学模拟；耗散粒子动力学模拟；剪切粘度 a b s t r a c t a b s t r a c t v i s c o s i t yi so n eo ft h em o s ti m p o r t a n tc h a r a c t e r i s t i c so ft h ep o l y m e rf l u i d s ，a n di t i sa l s oa k e yi n f l u e n c ef a c t o ro fp o l y m e rm a t e r i a lp r o c e s s i n g g e t t i n gt h ec o n s t i t u t i v e e q u a t i o nt h a tr e v e a l sm a t e r i a lp r o p e r t yh a sb e c o m et h ek e y t oc a r y yo u tp o l y m e r p r o c e s s i n gs i m u l a t i o nw i t ht r a d i t i o n a lf i n i t ee l e m e n ts i m u l a t i o n i th a sb e c o m ea n i m p o r t a n tr e s e r c hd i r e c t i o nt oo b t a i nt h es t r e s s - s t r a i nr e l a t i o n s h i po fc o m p l e xf l u i d u n d e rc e r t a i nc o n d i t i o n st h r o u g ht h em i c r o s c o p i cs i m u l a t i o no rm e s o s c o p i cs i m u l a t i o n m e t h o d ，s oa st oa v o i dt h ec o n s t i t u t i v ee q u a t i o n t h es i m u l a t i o no fp o l y m e rf l u i di so f t e nd i f f i c u l tt oa c h i e v eo n l yt h r o u g hs i n g l e s c a l i n gs i m u l a t i o n ，d u et oi t sl o n gc h a i n so fc h a i ns t r u c t u r ea n db i gs c a l es p a n i nt h i sp a p e r ，w ec o n d u c tan o n e q u i l i b r i u ms i m u l a t i o na b o u tp o l y m e rf l u i d t h r o u g hac o n n e c t i o no fm i c r o s c a l em o l e c u l a rd y n a m i cs i m u l a t i o na n dm e s o s c a l e d i s s i p a t i v ep a r t i c l ed y n a m i c ss i m u l a t i o n a n dw e f o c u so nt h ev i s c o s i t yi no u r s i m u l a t i o nm o l d t h em a i nw o r ki nt h i sp a p e ra rea sf o l l o w s ： 1 ) w es e tu pt h em o l e c u l a rs i m u l a t i o nm o d e lo fp o l y e t h y l e n e ，a n di m p l e m e n tt h e m o l e c u l a rd y n a m i c ss i m u l a t i o nw i t ht h es o f t w a r el a m m p s s t r u c t u r ei n f o r m a t i o n ， t h e r m o d y n a m i cp r o p e r t i e s ，a n dt r a n s p o r tp r o p e r t i e sa r ea n a l y z e di nd e t a i lu n d e r t h ec o n d i t i o no f4 5 0k1a t m ，w i t ht h es y s t e mo f9 6p ec h a i n s ，o n eo fw h i c h c o n t a i n9 6u n i o na t o m s t h er a d i a ld i s t r i b u t i o nf u n c t i o n ，m e a ns q u a r er o t a r y r a d i u s ，d e n s i t y , a n do r i e n t a t i o no fp o l y e t h y l e n ea r eu s e dt oc h a r a c t e r i z et h e m o d e l 2 ) t h r o u g ht h em e t h o dt h a tb a s e do nt h er a d i a ld i s t r i b u t i o nf u n c t i o n ，w es e tu pd p d s i m u l a t i o no fc o n s e r v a t i v ep o t e n t i a le n e r g y , a n dr e a l i z e st h ec o n n e c t i o nb e t w e e n m o l e c u l a rd y n a m i c ss i m u l a t i o na n dd i s s i p a t i o np a r t i c l ed y n a m i c s s i m u l a t i o n s e v e r a lk e yp a r a m e t e r sh a v eb e e ni n v e s t i g a t et oo b t a i nt h ep r o p e r s i m u l a t i o np a r a m e t e r st h r o u g hc o m p a r i n gw i t ht h er e s u l to fm ds i m u l a t i o n 3 ) w es e tu pt w ok i n d so fd i f f e r e n tm o d e l so f p o l y m e rf l u i df l o w ，a n dc a l c u l a t e d t h ev e l o c i t yd i s t r i b u t i o n ，s h e a rs t r e s sd i s t r i b u t i o na n ds h e a rv i s c o s i t yi nt h e i i a b s t r a c t p r o c e s so ff l o w , t h e nt h er e s u l tw i t ht h es i m u l a t i o na r ec o m p a r e dw i t hr e s u l t so f t h et h e o r e t i c a la n a l y s i s ，w h i c hs h o wt h a tb o t ho ft h em o d e l si n d i c a t et h es a m e v a r i a t i o nt e n d e n c yi nv e l o c i t yd i s t r i b u t i o na n ds t r e s sd i s t r i b u t i o n k e yw o r d s ：p o l y m e rf u i d ，m o l e c u l a rd y n a m i cs i m u l m i o n ，d i s s i p a t i o np a r t i c l e d y n a m i c ss i m u l m i o n ，s h e a rv e c o s i t y i i i 第一章绪论 第一章绪论 1 1 计算机模拟技术 计算机模拟技术是随着计算机技术及理论化学的发展而产生的一门新兴的 技术。计算机模拟技术也被称为“计算机实验”【l j ，它将实验化学与建立在第一 性原理上的理论方法联系起来，在模拟计算与实验对比中，将大量实验数据统 一解释，成为规律，深化认识。如今，计算机模拟已被几乎所有的应用科学领 域所采用，如材料科学、生命科学、医学、天体学、气象学以及建筑、工程学 等方面。 随着计算机技术的飞速发展，以及模拟技术的得天独厚的优越性，计算机 模拟技术发展迅猛，计算机模拟技术已经成为继理论和实验之后的第三大科学 研究手段，并与它们成为“三足鼎立”之势【2 训，它们之间相互关联，相互补充， 有效地促进了各自的发展。 计算机模拟技术和实验、理论分析一样都是认识自然世界的手段。在实际 的研究过程中，由于自然界的复杂性，往往一个物理现象被许多个因素影响， 并且这些因素相互关联、制约，而当我们用实验的手段来研究这个问题时，就 需要消去或控制某个或某些因素的影响，然而由于外部条件或实验手段的限制， 无法很好地控制这个影响因素，这样就导致利用实验来解释现象时的局限性。 就以研究聚合物为例，聚合物材料的功能和性质，不仅仅依赖于聚合物分子链 的性质( 聚合物单体种类、聚合度、单体的排列方式等) 、分子链的聚集状态等 微观因素相关，还与宏观的加工过程、外界环境等诸多因素相关，研究这些因 素对材料功能的影响对于开发所需的新材料至关重要。然而，对于实验来说， 各种影响因素的控制受到客观条件的限制，比如说我们很难得到聚合度相同的 聚合物，从而精确地研究聚合度对聚合物性能的影响。另外检测手段也在一定 程度上限制了我们对微观结构的深入研究。计算机模拟技术却在微观结构的模 拟以及消除外部因素的影响方面优势突出，因此计算机模拟是实验研究手段的 一个有效的补充。另一个方面，计算机模拟的模型的是否合适通常需要通过实 验的结果进行验证。 另外，理论分析通常是建立在高度简化的模型之上的，如果利用实验的手 第一章绪论 段来对理论分析的验证通常需要非常苛刻的实验条件，比方说完全真空状态、 绝对零度等理想的实验条件以及观测变化极其快的化学反应过程，在现实中是 无法实现，有时即使能达到非常近似的实验条件也是需要耗费大量的时间和财 力，而利用计算机模拟技术来实现却是非常简单的事情，因此计算机模拟技术 对于理论分析来说也是一个有效的促进手段。相反，计算机模拟技术是建立在 理论研究提供的数学模型之上的，理论研究的进步也会反过来促进模拟技术的 提高。 近年来，计算机在高分子研究中的应用发展很快【5 6 】，建立在严格的量子力 学及统计力学基础上的计算机模拟方法，可以对真实高分子材料的各种性质， 例如，材料的化学性质、物理性质和力学性质，进行深入、系统的分析，得出 结构和这些性质之间的关系，并且总结出内在规律。在某些方面聚合物材料的 计算机模拟以及物理性质的模型化研究己经超前于实验，其结果可以为材料的 设计以及实际应用提供相应的参考【7 。9 1 。另一方面，计算机模拟还可以应用在新 型材料的开发，在内容上既兼顾于各种材料的共性，又兼顾于材料科学所涉及 各个领域，巧妙地把各个学科的广泛性和多学科的相互渗透性联系起来，利用 计算机预测材料的结构和性能，既省时有省力，方便快捷，又降低了成本，可 以更好地指导实验，从而减少了开发的周期，提高了经济效益。另外，计算机 还能结合材料科学中的新方法、新技术，培养了科研工作者的创新意识，辅助 科研工作者们发展新的理论方法，为新方法的创新和发展提供了巨大的潜力和 广阔的空间。 为了模拟不同尺度下的材料的结构和行为，计算机模拟技术按照尺度可以 分为微观尺度的模拟、介观尺度的模拟、以及宏观尺度的模拟。在微观尺度方 面，已经建立起了比较成熟的模拟方法，如分子动力学模拟( m o l e c u l a rd y n a m i c s ， m d ) 方法和蒙特卡罗( m o n t ec a r l o ，m c ) 方法【i 。3 】。此外还有，量子力学( q u a n t u m m e c h a n i c s ，q m ) 方法【1 0 1 、布朗动力学( b r o w n i a nd y n a m i c s ，b d ) 方、法【j 。这 些微观尺度的模拟方法为我们研究微观结构提供了有力的工具。在介观方面， 比较常见的方法有：耗散粒子动力学( d i s s i p a t i v ep a r t i c l ed y n a m i c s ，d p d ) 方 法【1 2 1 4 1 、格子玻尔兹曼方法( l a t t i c eb o l t z m a n nm e t h o d ，l b m ) 1 5 】。介观模拟 方法能够模拟毫秒、微米级尺度的体系，并且能作为微观与宏观模拟联系的纽 带。在宏观方面，有基于流体力学的有限元分析方法，这种传统的方法已经非 常成熟，在聚合物加工中也得到了广泛的实际应用。 2 第一章绪论 1 2 利用多尺度模拟研究聚合物流体的意义 研究聚合物流体的流动行为对于聚合物的加工有着重要的意义，从本质上 讲，加工方法能够改变制品的性能也是由于加工过程中所加的外部作用对材料 的聚集状态或分子结构造成了改变，因此分子结构及其聚集状态是决定性能的 根本，研究微观、介观以及宏观尺度下的结构变化，揭示成型过程中形态结构 生成、演化的规律，及其与成型工艺之间的关系，对于制备性能优良的制品具 有重要的意义。 另外，对聚合物流动的模拟目前主要是用连续介质力学理论，通过推导微 分型或积分型的近似本构方程来描述应力和速度场之间的关系，然后结合流动 控制方程进行求解。然而在一般情况下，近似封闭本构方程是非常困难的，并 且即使得到本构方程，对于微流动，复杂流动等粒子特性明显且有壁面滑移的 情况，n a v i e r - s t o k e s 方程己不再适用【1 6 】。 多尺度贯通的模拟技术，由于其能够从构成体系的基本粒子的运动出发， 来分析结构以及性能，因此能够揭示微观结构与性能的关系。另外，多尺度贯 通的模拟技术对粒子的类型以及流动区域的大小没有限制，因此在研究复杂流 体以及微区域的流动问题方面也有着很大的优势。 1 3 模拟聚合物流体所面临的问题 聚合物流体和简单流体不同，首先从微观上看，聚合物流体是由大量的聚 合物分子组成的，我们知道聚合物分子结构具有多样性，并且构成聚合物流体 的聚合物分子的分子量不均一，常存在较宽的分布，分子结构的多样性以及分 子量的宽分布导致了其呈现出多种构象，体现在模拟上就是巨大的自由度。从 宏观上看，聚合物流体是一种非牛顿流体，存在剪切变稀以及由弹性变形引起 的诸如爬杆效应等特殊的流变现象。 正是由于聚合物流体的这些特征，导致了模拟聚合物流体的流动行为非常 不容易。从建模上来说，若建立全原子模型，当然能出现比较好的结果。但建 立这样的模型来模拟流动行为，在不作大量简化的基础上，基本不可能实现， 这是因为聚合物链的回转半径达到1 0 0 h ，聚合物链的松弛时间最少也达到1 0 。5 s ， 对于分子量稍大的聚合物分子可达到秒级，而全原子模型模拟的时间为l o - i s s ， 要能预测粘度则至少需要几百个松弛时间，那么模拟步数将达到1 0 ”步，而每 3 第一章绪论 一步的计算时间又与原子个数的成几何级数增长，这样的计算量已经远远超过 了计算机的计算能力。因此，利用简单的全原子模型是很难实现模拟聚合物流 体流动行为的。可以看出，聚合物流体的流动行为，从模拟的角度看，是一个 涉及到多尺度的问题，单用一种尺度的模拟方法是很难实现的。研究者们开始 把目光转向了多尺度模拟，多尺度模拟方法从微观结构入手，通过各个尺度的 连接建立起微观与宏观的联系，从而实现从根本上来研究宏观性质。而要研究 这些问题，建立聚合物多尺度耦合的模型是基础。目前，在单一尺度上的研究 比较成熟，多尺度的贯通还没有完全的解决。已有学者尝试了从微观的分子模 拟开始，通过粗粒化方法建立起进行介观模拟的力场，介观模拟能提供剪切粘 度与剪切速率的关系，再把此关系导入到宏观的有限元分析方法中，从而实现 模拟聚合物流体的流动行为，但这种方法还处于探究的阶段。相对而言，从微 观到介观的贯通发展地比较好，从介观到宏观的连接研究还远远不够。目前， 还没有形成一个较为成熟的方法来把各个尺度连通起来。不过多尺度贯通模拟 使从分子设计到材料j j it 全程贯通成为可能，已经成为了当前高分子物理研究 的热点。 1 4 基于粗粒化模型的耗散粒子动力学方法 在高分子多尺度模拟领域，许多研究工作者已经从不同的角度提出了许多 有意义的方法，这些方法主要可以分成两类：建立在结构性质粗粒化上的粗粒 化模拟和建立在非结构性质基础上的粗粒化模拟方法 1 7 - 2 0 】。在粗粒化方法中， 有些是建立在径向分布函数( r d f ) 的基础上的，有些是建立在作用力数据基础 上的。粗粒化方法是实现从微观到介观衔接的手段，聚合物链的粗粒化过程就 是通过整合链上的大量微观细节，将聚合物分子描述成由不同柔性弹簧连接的 粗粒化粒子( 或珠子) 组成的链，并用形式简单的经验势或精确的势表描述珠子间 的相互作用。由于粗粒化模型比原子模型具有较少的自由度和较软的相互作用 势，因而能达到的模拟长度尺度为几个微米，模拟时间尺度为几个微秒。在过去的 二十年里人们采用聚合物粗粒化模型和分子动力学方法研究了高分子链在熔体 中的运动方式【2 m 3 1 ，发现从r o u s e 链到r e p t a t i o n ( 蛇形) 链的动力学转变区域，直接 观察到缠结链在管道中运动的分子图像。 耗散粒子动力学( d p d ) 方法可以看做是“粗粒化”的分子动力学模拟( m d ) 4 第一章绪论 方法，他具有m d 方法以及格子玻尔兹曼( 1 a t t i c eb o l z e m a n nm e t h o dl b m ) 的 特征，与m d 方法相比，在d p d 方法中引入了耗散力和随机力来补偿由于粗粒 化导致的自由度的降低，与l b m 相比，d p d 方法继承了l b m 方法能描述流体 行为的特点，又改进了l b m 所采用的格子模型所带来的在模拟聚合物流体时的 缺陷。 正是由于以上的原因，以粗粒化来实现微观与介观的连接，以耗散粒子动 力学方法来实现流体模拟将是一个非常有意义的尝试。 1 5 利用d p d 研究聚合物流体的进展 利用d p d 研究聚合物流体始于1 9 9 5 年，s c h l i j p e r ，h o o g e r b r u g g e 以及m a n k e 首次利用d p d 方法来研究聚合物【2 4 1 ，他们采用线性刚性弹簧模型对链长为3 3 到5 8 8 9 个粒子的链，采用e u l e r 积分算法，进行了模拟，并将结果与r o u s e 模 型进行了对比，得出标度为y o 5 2 。k o n g 2 5 】等人随后同样采用了刚性弹簧模型 对不同种类的聚合物溶液体系进行了考察，得出标度指数为y 0 6 。k o n g t 2 6 】还 利用同样的模型考察了在稳态剪切下的剪切变稀现象，并首次计算了法向应力 差系数。g r o o t 和w a r r e n 在1 9 9 7 年的文章【2 7 】详细分析了d p d 模拟方法，并将 d p d 模拟结果和f l o r y - h u g g i a s 理论相比较，推导出了模拟参数与实际聚合物体 系的映射关系，这种关系的建立为大尺度的模拟建立了基础。在1 9 9 9 年，g r o o t 【2 8 】 等人利用d p d 方法模拟了嵌段共聚物的相分离，并将d p d 方法和布朗动力学 ( b d ) 方法进行了对比，发现流体动力相互作用在d p d 中是必不可少的。 b o s c h l 2 9 应用d p d 方法对粘弹性流体进行了研究，发现d p d 方法的结果不仅与 n a v i e r - s t o k e s 方程的结果一致，还与宏观模型预测结果相一致。到了2 0 0 0 年， s p e n l e y 3 0 1 采用珠簧模型计算了聚合物的排除体积和自扩散系数，计算结果表明， 自扩散系数和理论分析的结果吻合得很好。 1 6 本文的研究内容及结构 粘度函数是聚合物流体的一个重要表征参数，也是聚合物加工成型性能的 重要影响因素。在进行传统的聚合物流动有限元模拟过程中，获得揭示材料物 性的本构方程( 粘度函数) 已经成为了模拟的关键1 6 】。通过分子模拟或介观模 5 第一章绪论 杂流体模拟的重要方向。本文将以计算聚合物流体的粘度为主要目的。 正如上文所提到的那样，由于聚合物分子的长链结构，导致时间尺度和空 间尺度都较大，要实现计算聚合物流体的粘度采用单一尺度很难实现，因此在 本论文中，采用了微观和介观相耦合的粗粒化的耗散粒子动力学模拟方法来进 行模拟，基于r d f 的粗粒化的方法【3 l 】在保持系统性质的基础上大大扩大了模拟 的尺度，耗散粒子动力学模拟方法在流体动力学方面具有一定的优势。我们希 望通过这样的方法来探究聚合物流体的流动行为。 聚合物流体的流动是一个非平衡态的问题，对于非平衡态的模拟无论是在 实施的方法上还是在理论上都与平衡态的模拟有很大的不同，如何实旌非平衡 态模拟也是本文的一个研究方面。我们分别采用了平板拖曳流模型和反转泊萧 叶流模型来进行非平衡态模拟，并对这两种模型进行了比较。 本文的结构安排如下： 第二章：介绍本文中用到的模拟方法微观的分子动力学模拟、介观的耗散 粒子动力学模拟的方法及其理论，接着，我们还对粗粒化模型的建立方法作了 简要的介绍。 第三章：我们建立了短链的聚乙烯分子模拟模型，采用联合原子力场对该体系 进行了分子动力学模拟，通过模拟我们得到了体系的结构信息( 径向分布函数， 均方回转半径、松弛时间) 、输运性质( 扩散系数) 。 第四章：本章主要介绍贯通微观的分子动力学和介观的耗散粒子动力学方法。 并进行了平衡态的耗散粒子动力学模拟，由于采用不同的耗散粒子动力学模拟 参数对所能模拟的最大粘度有影响，因此我们在这一章里对模拟参数的选择进 行了讨论，并将通过此模拟得到的结构信息与第三章通过分子动力学模拟所得 的结果进行了对比，以证明粗粒化方法的合理性。 第五章：在这一章里我们采用了两种不同的非平衡实施方法，进行非平衡态的 粗粒化耗散粒子动力学模拟，通过模拟来得到聚合物流体的粘度，并对所得的 结果与理论模型预测结果以及理论分析的结果进行了对比。讨论了不同的剪切 速率对粘度的影响。最后，对本文进行了总结和展望。 利用基于r d f 的粗粒化耗散粒子动力学方法对聚合物流体进行模拟还处于 探究的阶段，还有很多需要讨论和研究的问题。但我们也看到，这种方法在研 究聚合物流体方面还是很有前景的。 6 第二章理论基础及模拟方法简介 第二章理论基础及模拟方法简介 本论文主要采用分子动力学和耗散粒子动力学贯通的模拟方法来研究聚乙 烯的流动行为的，因此我们有必要对分子动力学和耗散粒子动力学的模拟的原 理、以及他们相互贯通的方法有一定的了解。另外，考虑到流体的流动行为一 般是在一定的外部作用下( 如重力，压力等) 才发生的，因此模拟流体的流动 行为大多采用非平衡态模拟方法，因此在本章也会对非平衡态的模拟方法做一 些简要的介绍。 2 1分子动力学模拟( m d ) 分子动力学模拟是一种通过计算原子、分子之间的相互作用力随时间的演 化，来得出系统的各种相关物理信息的计算机模拟技术。它采用力场的方法， 通过求解牛顿运动方程，模拟体系随时间的演化过程，得到各个时刻的粒子位 置、速度、作用力等信息，利用统计力学的方法获得体系的宏观性质( 如温度、 压力、能量等) 。分子动力学方法已经被广泛地应用于生物、材料、化学合成等 各个领域，它已经作为一种研究手段，能探究某些实验无法观测或计算的理化 现象。 2 1 1 基本原理 分子动力学的原理即是利用牛顿运动定律来计算原子的运动过程。考虑含 有n 个原子或分子的系统，系统的能量由系统中所有的原子或分子的动能和势 能的总和，其中势能由分子力场( 见2 1 2 节) 进行定义，依照经典力学，系统 中任一原子f 所受的力为势能的梯度： 亏= - v i u = - f 毒+ 了毒靠毒，u q 再由牛顿运动定律可以得到原子f 的加速度a t ，将牛顿运动方程对时间进行 积分： 7 第二章理论基础及模拟方法简介 d 2 一d r ,一 孑2 言钏， ；产v i + 二f f( 2 一2 )y ，=+ 口f f( ) 一一。一0 l 一， ，f2 ，：+ v it + 一口f f ， 上式中，一v ，o 中的上标表示前一时刻，取f ：& ，即通常所说的时间步，通过 上式的计算即能得到经过础的时间间隔后，粒子新的位置和速度，再由新的位 置计算相互作用势、作用力、加速度，如此重复进行就能得到各个时间点下粒 子位置、速度、作用力的信息。对得到的轨迹进行时间平均，即可得到宏观的 物理信息。 2 1 2 分子模拟中的力场 分子力场【2 j 是进行分子动力学模拟的基石，它描述的是原子或分子问的势 能，决定着原子的拓扑结构与运动行为。 一般来说，定义一个力场应该包含以下几部分： ( 1 ) 原子类型列表：原子类型由特定力场定义，它不仅由原子的原子数决定， 通常还包含原子的杂化态信息，有时局部化学环境的信息也包括在内； ( 2 ) 如果在原子类型列表中不包含电荷信息，那么原子电荷可单独列表； ( 3 ) 确定原子类型的规则； ( 4 ) 能量( 势能) 表示的函数形式； ( 5 ) 各函数项的参数； ( 6 ) 有些力场给出产生未知参数的规则； ( 7 ) 有些力场给出标识函数形式和参数的方法。 其中，势函数形式和相应的参数是决定力场有限性的关键因素。 经典力场( 如d r e i d i n g 力场) 采取的函数形式较简单，没有考虑交叉项， 体系的势能由键能和非键能组成，键能一般包括键伸缩、键弯曲、二面角扭转、 离平面振动能，非键能通常包括范德华作用能、库伦作用能以及氢键作用能。 总的势能可以用下式表示： 厶埘甜。也m k w + p - l o n n o n a ，、 = ( e 乙耐+ 五厶幽+ 五，o ，s i 伽+ 瓦嘲) + ( e + 6 + ) 、。7 二代力场如c f f 、p c f f 、c o m p a s s 都包含了交叉项，用于考虑临近原子 8 第二章理论基础及模拟方法简介 对键角、键弯曲等的影响，其设计的目的是为了能精确地计算分子结构、光谱、 晶体性质等。但对于结构简单的分子，经典力场和二代力场的结果差别很小。 在分子模拟中，势函数的选取与动力学的计算关系密切，选择不同的势函 数，就会得到不同的势能面，原子或分子的运动轨迹就有可能不同，进而会影 响到最终的统计结果。 2 1 3 积分算法 所谓运动方程的积分算法就是指根据当前的粒子的位置、速度、受力来计 算经过f 时间间隔后粒子新的位置、速度以及所受的作用力。目前比较常用的 算法有v e r l e t 算法、v e r l e t 1 e a pf r o g 算法、v e l o c i t y - v e r l e t 算法【3 2 。4 1 。在本节我们 将介绍这三种积分算法的原理。 ( 1 ) v e r l e t 算法 其基本原理是对位置坐标函数用t a y l o r 公式在t + a t 处展开得到： ，( t + a t ) = 心) + 丢心) 蚺嘉川) + 刍心) a t 3 + o ( 4 ( 2 4 ) 同理在t a t 处展开得到： 以础) = 心) 一丢m ) 蚺鲁心) 舻一蠡心) a t 3 + o 州a t ( 2 5 ) 将以上两式相加，且根据，( f ) = a ( t ) = f ( t ) m 即可得到位置更新的表达式： ，o + f ) ：2 ，o ) 一厂o 一f ) + 羔尘f 2 + o ( a f ) 4 ( 2 6 ) 新位置的估计误差是( 垃) 4 ，可以看出，新位置的坐标是根据前两步的位置 坐标和前一步的作用力来给出的，而不需要知道速度，t 时刻的速度表达式如下： y “) ：r ( t + a t ) - r ( t - a t )( 2 7 ) 2 a t 这种方法的不足之处是计算t 时刻的值需要知道t + a t 时刻的位置，这给实 际操作带来不便。 ( 2 )v e r l e tl e a p - f r o g 算法 此算法是在v e r l e t 算法的基础上进行改进而获得的。首先定义半个时间步长 的速度为： 9 第二章理论基础及模拟方法简介 1 ，“+ 三f 1 善r ( t + a t ) - r ( t ) 2a t 由此式推导出位置表达式为： 1 r ( t + a t ) = 厂( f ) + v ( t + 寺a t ) a t 二 将此位置表达式代入v e r l e t 算法中的位置更新表达式中， l e a p f r o g 算法的速度更新表达式为： ( 2 8 ) ( 2 9 ) 我们可以导出v e r l e t 地+ 互1a t ) _ v ( r 一扣+ 等缸 ( 2 1 0 ) v e r l e tl e a p f r o g 算法的缺陷在于体系的位置和速度的计算存在不同步性。因 此不能计算某一时刻的总能量。 ( 3 ) v e l o c i t y v e r l e t 算法。 v e l o c i t yv e r l e t 算法的过程是先根据上一步的位置、速度、力来计算下一步 的位置，再由下一步的位置计算作用力，再更新速度。其t + a t 时刻位移和速度 公式如下： ，( f + f ) - - - - - r ( f ) + 加( f ) + 掣a f 2 ( 2 1 1 ) 嵋+ 归印) + 纠盟m + 螋m ) ( 2 1 2 ) z v e l o c i t yv e r l e t 算法克服了v e r l e tl e a p f r o g 算法中体系的位置和速度的计算存在 不同步性的缺陷。 从物理意义上说，以上三种积分算法都把垃时间段内粒子的运动看作是加 速度不变的运动，实际上作用力是会随位置不同而变化的。因此，& 不宜取得 过大，一般取为体系最高频率的1 1 0 左右，只有这样才能满足加速度不变的假 设。 2 1 4 周期性边界条件及最小镜像原理 在通常的模拟过程中，我们希望粒子不会跑出模拟区域，但如果我们没有 在模拟区域的边界处设计一些条件，粒子按照牛顿运动定律运动的话肯定会跑 出模拟区域，使模拟系统的密度发生变化，导致系统失真，这种在边界处设计 的条件我们称之为边界条件。边界条件根据使用的目的不同，分为很多种，常 用的边界条件为周期性边界条件、固定边界条件、自由边界条件。在我们的模 1 0 第二章理论基础及模拟方法简介 拟过程中使用的是周期性边界条件，因此我们将对其进行详细地介绍。 图2 1 为周期性边界条件的示意图，如图示中央的黑色盒子是模拟区域，其 周围的盒子( 周期性镜像系统) 具有和中央盒子相同的运动和排列，当模拟区 域中有任何一个粒子移出盒子( 如2 号粒子) ，则必有一个镜像粒子从相反的方 向进入模拟区域，这样的限制条件就能保证系统的密度保持恒定。 处于边界处的粒子由于有限的尺寸限制，导致所受力与处于模拟区域的中 央位置的粒子不同，而这并非事实，因此计算系统的相互作用力的时候，常采 取最近镜像方法，如图2 2 所示，若计算粒子l 与粒子3 之间的作用力，选取的 是与粒子1 距离最近的镜像，即d 盒中的粒子3 。同样，计算粒子3 与粒子1 的作用力时应取模拟系统中的粒子3 与e 盒中的粒子l 。 由于最近镜像概念的引入，我们需要采用截断半径( c u t o f fr a d i u s ) 的方法 来计算体系中的远程作用力，否则会重复计算相同号码粒子间的作用力，最终 导致不j 下确的结果。截断半径的最大值不能超过模拟盒子的一半，否则会导致 相同号码作用力的重复计算，截断半径的设置至少保证粒子所受力基本符合事 实，并且要在符合事实的前提下尽可能地取小些，这样使得每次参与计算的粒 子数就少，计算的效率就高。一般原子的截断径为l o a 左右。 当当 吕 9 o f l 鞘糊 名罄 厶 圈 l 2 名参名兽 名参 ( )囝0 i ，、 、一， 一工巾 图2 1 周期性边界条件示意图 第二章理论基础及模拟方法简介 名当昌当 a 0 8c 内蓦e i 分 ? 曩 l 1 d 昌昌蓦 n0 h fu 一， 一 第二章理论基础及模拟方法简介 这些力都是对作用力，整个体系随时| 1 自j 的演化过程中满足牛顿运动方程： 堕：m ( 2 1 3 ) j2 m l 厶， d l 朋鲁2z 。善( 够+ 带+ 劈) ( 2 “) 其中，r t 、m 分别是位置和速度矢量，f i ：是i 粒子所受的总的作用力，为了节约 计算时间，这三个作用力都设定在一定范围内起作用，这个距离称为截断半径c ， 这个距离作为体系中唯一的空间尺度标准，将其作为长度约化单位，为了简化， 通常取c = l 。 保守力是作用在粒子中心连线上的软的作用力，经典d p d 方法中，其作用 力可以表示为： 厂a o ( 1 一白) 露 ( r , s ) 劈= ( 2 1 5 ) 。 l 0 ( 吩) 其中，表示f 粒子与粒子之间的排斥强度，弓= 亏一弓，勺= 恳i ，0 = 弓引 保守力和系统的热力学性质密切相关，为了得到能反应真实的体系的热力 学性质，19 9 7 年，g r o o t 将d p d 方法与f l o r y - h u g g i n s 理论耦合在一起【2 7 1 ，从而 建立起真实聚合物体系的原子和分子信息与d p d 粗粒化模型之间的映射关系， 得到了体系的a ，的计算方法，在模拟相变问题时能得到比较好的结果，对于模 拟聚合物熔体的流动行为，近期的文献提出了基于径向分布函数的迭代法来确 定保守力。 涨落耗散理论的引入使得d p d 方法适合于模拟体系的流动行为，其中耗散 力的作用是阻止粒子之间的相对运动，其表达式为： 乒；d = 一胗d ( 吩) ( 疗露) 易 ( 2 1 6 ) 其中，彩“依赖于r 的权重函数，它描述了耗散力随着粒子间距离的增加而衰减 的情况，当芬时，其值为零。v , s 为粒子z 和_ ，之间的相对速度，巧= 谚一巧， y 是摩擦系数。 由于耗散力的阻止作用，会使得系统的动能降低，在d p d 方法中，这个降 低的动能会通过随机力进行补偿，随机力可以表示为： 髫= 俐足( 勺) 岛0 ( 2 1 7 ) 1 3 第二章理论基础及模拟方法简介 其中仃是噪声强度，缈且是依赖于r 的权重函数，岛，是一个平均值为零 ( ( 岛( f ) ) = 0 ) ，且满足： ( o i j ( t ) o i j ( t ) ) = ( 瓯+ 4 ，屯) 万( f t ) ( 2 1 8 ) 的随机数。 与保守力一样，耗散力和随机力作用方向也是沿着粒子对的连线方向，也 是对相互作用，e s p o n a l 和w a 盯翎【2 7 】通过证明指出耗散力和随机力的权重函数可 以随机选择，但为了满足耗散涨落定理，使系统稳定，保证系统满足n 系宗 的统计规律，两权重函数之间、摩擦系数与耗散系数之间必须满足一定的关系： 缈d ( 厂) = - 缈r ( ，) 2 ( 2 1 9 ) o r 2 = 2 y k 占t ( 2 2 0 ) 通常，权重函数可以选择如下形式： 厂( 1 一互) ” ( 勺 ) 0 9 d ( ，- ) = ( 2 2 1 ) l 0 ( 吩) 其中m 为参数，显然m 越小则耗散力越大，传统d p d 模拟中，m 通常取2 。 从以上的运动控制方程可以看出，d p d 有两个特色： ( 1 ) 引入了非常“软”的作用势，这样就使得我们模拟的时间尺度和空间 尺度增大，我们可以用比较大的时间步长，从而减小我们的计算量。 ( 2 ) 耗散力和随机力相互耦合，起到热浴的作用，并且由于是粒子对相互 作用，满足牛顿第三定律，从而整个体系的动量是守恒的，引入的这个热浴， 在满足涨落耗散关系的条件下，使d p d 方法可以正确地表述动量传递，而这一 点对于复杂流体的动力学是非常重要的。 2 2 2 积分方法 d p d 方法也是基于牛顿运动方程的，象m d 模拟一样，采用什么方法对运 动方程进行积分是非常重要的内容，对于m d 方法，其积分方法研究和应用都 比较成熟了。d p d 方法相对于m d 方法其积分方法要更复杂些，这主要是由于 d p d 方法中的耗散力与粒子的相对速度有关，但作用力与与速度的相关性在积 1 4 第二章理论基础及模拟方法简介 分方法上难以实现，并且导致模拟得到的体系的物理量存在错误或偏差。人们 在研究的过程中，提出了一些处理这个问题的方法，最早的是基于m d 模拟所 用的v e l o c i t yv e r l e t 而提出的d p d w 算法，后经g r o o t 和w a r r e n 对其进行改进 形成了g w - w 算法，后者是在目前d p d 模拟中比较常用的方法，我们将对这 两种积分方法作简要介绍。当然，d p d 模拟的积分方法还有s h a r d l o w 提出的s i 方法f 3 7 】，l o w e 提出的l a 方法【3 引。由于我们的模拟不涉及到，所以在此处不作 介绍。 2 2 2 1d p d w 积分法 这种算法的基础是v e l o c i t yv e r l e t ( w 算法) ，和w 算法相比，只是引入 了耗散力对速度的依赖性，这个算法的优点是速度快，不需要任何参数，但对 于不同的时间步长及大的计算步数其稳定性较差。其积分过程如图2 3 所示。首 先根据当前的作用力来更新速度以及位置，再根据新的位置和速度来计算作用 力，再利用作用力对速度进行更新，最后根据更新后的速度来计算耗散力( 耗 散力与速度相关) 。 ( 1 ) e 卜t + 互l 小l ( - c a t + 霞。a t + 霉旯伍) ( 2 ) 亏+ - 亏+ t f ( 3 ) c a l c u l a t e 霞c 弓) ，霉。 弓，t ) ，j ：尺 弓 ( 4 ) e 卜e + 专砉( 霉f a t + 霉。a t + 霉矗压) ( 5 ) c a l c u l a t e 霞。托i 图2 3d p d w 算法 2 2 2 2g w - v v 积分算法 g w - v v 算法是对d p d w 算法的一种改进，这个算法充分考虑了耗散力对 速度的依赖性，其过程如图2 4 所示。首先根据t 时刻的位置、速度、作用力来 计算t + a t 时刻的位置，为了得到t + 血的速度，需要t 时刻的速度、t 时刻的力 以及t + a t 时刻的力，而在计算t + a t 时刻的力的时候由于耗散力对速度的依赖 性，因此需要t + a t 时刻的速度，因此在计算t + a t 时刻的力之前做了一个速度 的预测霹o + a t ) ，最终再进行一次速度的校正。在这个过程中，力只更新了一次 1 5 第二章理论基础及模拟方法简介 ( 计算作用力比较费时) ，因此这种算法并没有增加计算时i 可。在做速度的预测 的时候引入了可调参数兄，g r o o t 和w a r r e n 通过计算比