（高分子化学与物理专业论文）多链聚乙烯单轴拉伸形变过程的monte+carlo模拟.pdf

复且大学硕士论文摘要 摘要 高分子材料具有独特的力学性能，其形变时的应力应变曲线相对于无机和金属 材料而言，有很大的特殊性。随着高分了材料在生产和牛活中的广泛应用，这种力 学性能一直是理论和实验中的研究热点。 计算机模拟是研究这种形变过程的常用方法。其优点在于可以通过模型化研究， 揭示物质宏观现象与其微观结构的相互关系。分子动力学( m d ) 和m o n t ec a r l o ( m c ) 是常用的两种研究高分子材料的计算机模拟方法。前者注重考察模型随时间的演变 过程，适合于研究物质的动态过程：后者是基于一种随机概念的方法，注重考查体 系的构象状态，适合于计算物质的统计平均属性。由于m d 方法存在体系结构的演 变效率问题，即需要较长的时间使体系完全脱离初始结构，对计算机硬件和研究时 间的周期提出了考验；m c 方法则避免了这种问题，因为基于m c 方法设计的结构 运动规律，可以人为的从大尺度上演变体系结构，具有较高的效率。 本文中我们将用m c 方法，模拟高分子材料线形聚乙烯( p e ) 的形变过程。 通篇主要围绕三个问题讨论：( 1 ) 分子模型的内应力计算；( 2 ) 弹性分了键模型和 刚性分子键模型；( 3 ) 高级m c 算法在模拟过程中的应用。结合以上三个问题，考 查了p e 在单轴拉伸形变时的应力应变关系，及其温度和拉伸速度的依赖性等问题。 我们希望通过此研究，开拓性的把复杂但有效的m c 算法应用于模拟高分子模型的 形变过程，同时系统性的考查此算法应用于弹性分子键模型和刚性分子键模型的区 别，正确的计算不同模型下体系的内应力，为将来应用于复杂高分子体系奠定基础。 为合理的解决这三个问题，我们的研究遵循由简入难的原则，从简单体系开始，逐 渐增加体系复杂度，最终研究p e 体系。文章组织结构主要分为六个章节： 第一章是绪论，我们将对此领域已有的研究成果做一个简短的回顾，同时提出 对本文工作的目标。 第二章中我们讨论m o n t ec a r l o 模拟中分子模型体系的内应力计算方法。 第三章中，我们将用最基本的m c 算法模拟单原予分了体系的平衡过程。单原 子分子体系由单原子分子组成，分子间具有范德华力相互作用。其结构的简单性， 可以让我们初步认识m c 算法如何应用于分子体系，同时我们可以学习如何用经典 的维利定理( v i r i a lt h e o r e m ) 计算简单体系的内应力。模拟得到的平衡态量与己发 表的同体系m d 研究结果相比较，从而判断m c 算法和内应力计算的正确性。 第四章中，我们将体系复杂化至双原子分子。伴随着分子键的引入，我们着重 考察了弹性分子键模型和刚性分子键模型下的内应力计算。在正则系综( n p t ) 下， 第1 页 复旦大学硕士论文摘要 计算得到的体系内应力，应该与体系外部的应力相平衡，据此可以判断应力计算的 正确性。这样，我们初步了解了结构限制条件的存在与否，对内应力计算带来的影 响。 第五章将讨论线形p e 模型的建立、平衡和形变过程。我们将会仔细介绍所应 用的高级m c 算法，讨论弹性分了键模型和刚性分子间模型下的内应力计算方法。 然后着重考查单轴拉伸形变过程的应力应变关系，及其温度和拉伸速率的依赖性。 同时将试图解析形变过程中的微观结构转化因素，力图将宏观形变过程和微观结构 转变聪系起来。 第六章是结论和对未来工作的展望。 关键词：m o n t ec a r l o 模拟、聚乙烯、拉伸形变、应力应变关系 第1 i 页 复旦大学硕士论文摘要 a b s t r a c t e x i s t i n gm o n t ec a r l o ( m c ) a l g o r i t h m sf o rt h es i m u l a t i o no f u n i t e d a t o mm o d e l so f p o l y e t h y l e n e ( p e ) h a v eb e e ne x t e n d e dw i t hp r e s s u r et e n s o rc a l c u l a t i o n sa n du n i a x i a l d e f o r m a t i o n i no r d e rt ou n d e r s t a n dh o wt oc a l c u l a t ep r e s s u r ei nm c ，t w o l o w - m o l e c u l a r - w e i g h ts y s t e m s ，am o n a t o m i cl e n n a r d j o n e sl i q u i da n das y s t e mo f d i a t o m i cm o l e c u l e s ，h a v eb e e ns t u d i e df i r s t e s p e c i a l l yt h ei n f l u e n c eo fc o n s t r a i n t s ( f i x e d b o n dl e n g t h s ) i si n v e s t i g a t e di nt h es e c o n dc a s e r e s u l t ss h o wa g r e e m e n tb e t w e e nc a l c u l a t e di n t e r n a lp r e s s u r ea n di m p o s e de x t e r n a l p r e s s u r ei nt h en p t e n s e m b l eb o t hf o rt h em o n a t o m i cl e n n a r d j o n e sl i q u i da n df o rt h e s y s t e mo fd i a t o m i cm o l e c u l e s ( w i t ha n dw i t h o u tc o n s t r a i n t s ) f i n a l l yt h ep r e s s u r e c a l c u l a t i o ni ss u c c e s s f u l l yi m p l e m e n t e df o rp em o d e l s ( w i t ha n dw i t h o u tc o n s t r a i n t s ) u p o nu n i a x i a ld e f o r m a t i o nt h es t r e s sr e s p o n s et ot h es t r a i ni si n i t i a l l yl i n e a re l a s t i c ， s u b s e q u e n t l ya st h es t r a i ni n c r e a s e sf u r t h e ry i e l d i n gi so b s e r v e d ，a n df i n a l l ys t r a i n h a r d e n i n gi sd e v e l o p e d t h es i m u l a t e dy o u n g sm o d u l u sa n dp o i s s o nr a t i ot a k er e a l i s t i c v a l u e s f u r t h e r m o r et h et e m p e r a t u r e - a n ds t r a i n - r a t e d e p e n d e n c i e so fs t r e s s s t r a i nc h i v e s h a v eb e e ni n v e s t i g a t e da n dt h er e s u l t sa l ei nq u a l i t a t i v ea g r e e m e n tw i t ht h ee x p e r i m e n t a l o b s e r v a t i o n s c h a i nc o n f o r m a t i o na n de n e r g ya n ds t r e s sp a r t i t i o n i n gw i t hi n c r e a s i n g s t r a i na r ef o l l o w e di nd e t a i l d u r i n gt h ed e f o r m a t i o nt h ec h a i n sa d o p tm o r ee x t e n d e d c o n f o r m a t i o n sa n dt h ef r a c t i o no f d y a d si nt h et r a n ss t a t ei n c r e a s e s i nt h ee l a s t i cr e g i o n m e c h a n i c a lw o r kd o n eo nt h es a m p l e ，i sp r i m a r i l ys t o r e da sn o n b o n d e di n t e m a le n e r g n w h e r e a sf r o mt h ey i e l dr e g i o no n w a r d st h ei n t r a c h a i nc o n t r i b u t i o n ss t a r tt op l a yar o l e k e y w o r d s ：m o n t e c a r l o s i m u l a t i o n ，p o l y e t h y l e n e ，u n i a x i a ld e f o r m a t i o n ， s t r e s s s t r a i nb e h a v i or 弟i i l 页 论文独创性声明 木论文是我个人在导师指导卜进行的研究工作及取得的研究成果。论文中除 了特别加以标注和致谢的地方外，刁i 包含其他人或其它机构已经发表或撰写过的 研究成果。其他例志对本研究的启发和所做的贞献均已在论文中作了明确的声明 并表示厂谢意。 作者签名：秘 论文使用授权声明 本人完全了解复旦大学有关保留、使用学位论文的规定，即：学校有权保留 送交论文的复印件，允许论文被查阅和借阅；学校可以公布论文的全部或部分内 容，可以采用影印、缩印或其它复制于段保存论文。保密的论文在解密后遵守此 规定。 未p 作者签名：! 塑!导师签期：三竺：至重! ! o 复旦大学硕士论文绪论 1 1 引言 第一章绪论 高分子材料通常可以分成塑料、橡胶和纤维三大类，在我们的日常牛活中都是 不可或缺的。由于其具有的独特的化学和物理性质，他们已经被应用到各个领域。 从化学的角度说，高分子是一种长链状分予，具有很大的分子量，通常一个高分子 是由几千乃至上万个单体分了组成。2 0 世纪产牛了很多理论，用于揭示实验中得到 的高分子材料性能的本质，但是仍然存在着一些开放性的问题。 略 也 z 一 三 11 112 d r a wr a t i o 卜】 图1 ：三种常见高分子材料拉伸时的应力应变关系【l 】。 其中一个开放性问题就是高分子材料特有的应力应变关系。人们从实验中观察 到典型的玻璃态高分子材料，比如聚苯乙烯( p s ) 和聚碳酸酯( p c ) ，在拉伸形变 过程中表现出完全不同的力学性能 1 2 ：p s 是一种脆性的材料，在发牛很小的应变 时就会断裂；而p c 是韧性的，他可以支持比较大应变量( 如图l 所示) ；同时，经 过深入研究，人们发现这些特性还依赖于对实验样品的热处理和力学处理过程。高 分子材料典型的应力应变曲线是由四个区域组成：弹性形变区、应力屈服( 随着应 变增加，应力达到一个极大值) 、应变软化区( 应力随应变增加而减小) 和应变强化 区( 随着应变增加，应力又重新开始增加) 。人们试图从实验和理论角度来研究高分 子材料的这种现象，比如b o y c e 等口，4 1 发现屈服应力随着拉伸速率的增加、温度的降 低以及压力的增加而增加，而且应变强化行为强烈的依赖于形变的方式；m e l i e k 等 1 5 1 发现应变强化行为与形变后期分子链的交联或者缠绕网状结构相关。但是应变强 化模量与交联密度的精确关系，以及应变软化想象的物理根源，仍然没有被很好的 解释 5 1 。 第1 页共7 6 页 复旦大学硕士论文 绪论 计算机模拟提供了一条从分予层面上研究高分了材料形变过程的途径。他把理 论和实验结果联系起来，通过模拟结果和实验数据的对比，考查理论的正确性，或 者基于合理的理论，来预测不同实验条件下的实验结果。使用这种技术，人们可以 模拟各种极端实验条件，比如超低温、超高压等，其预测的结果对实验有很重要的 指导性；或者对于不能肯定的理论模型，通过模拟结果与实验数据的对比，测试理 论的正确性，比如要考查一个物质分子间相互作用理论，可以建立一个模型，基于 这种相互作用，模拟得到物质的诸如密度等属性，看其是否与实验数据相符合从 而判断此相互作用理论是否完善。常用的分子模拟技术有两个：分子动力学 ( m o l e c u l a r d y n a m i c s ，简称m d ) 和m o n t ec a r l o ( 简称m c ) 。 1 2 分子动力学简介 a l d e r 和w a i n w r i g h t 6 最早于1 9 5 6 年使用分子动力学模拟刚性球集合的动态过 程。分子动力学基于牛顿运动定律，解决体系随时间的演变规律，适合用来研究平 衡属性或者动态过程。他与常规实验过程非常相似，即制备一个样品，通过( 计算 机) 实验来得到体系某个属性随时间的演变，如果体系达到平衡，对数据平均可以 得到该属性的时间平均值。模拟时间越长，平均得到的数据就越可靠。 要进行分子动力学模拟，首先需要对研究休系进行初始化，即指定体系粒了的 初始位置和初始速度。对于初始位置的选择，通常以靠近平衡时的位置为原则，保 证粒子间没有明显的交叠或导致高能量状态的相对位置。以单原了分了模型为例， 通常初始位置被设置在立方点阵格子上。粒子的初始速度则以一个平均速度值为基 础，从他附近的一个平均分布中选择。 接着需要解决的问题是分子动力学模拟中最耗时的部分计算作用在每个粒 子上的力。以成对相互作用为例，由个粒子构成的体系，如果我们简单的考虑一 个粒子与其他所有粒子的相互作用，我们需要计算f n 一1 1 2 个相对距离，使用 这种方法计算作用力所需要的时间正比于2 。当然，存在一些算法技巧，可以使计 算时间缩减到正比于，在后文中我们会详细讲述。要计算一个力的x 分量，可以 简单的使用下述公式： 正( ，) ：一掣：一例掣1 ，( 1 - 1 ) “ 7 八o r ， 其中，和u f r l 分别代表粒子的位置和体系的势能。 第三个问题是通过积分牛顿运动方程，确定体系粒子在各个时间点的位置。需 第2 页共7 6 页 复旦人学硕士论文 绪论 要注意的是牛顿运动方程是具有时间可逆性的，而积分算法同样应该具有时间可逆 性，即将来和过去两个时间点的相空间坐标体现在算法中应该是对称的，否则，如 果从将来的一个状态点反推回来，将不能得到与正向运动相同的轨迹。最简单而且 最好的算法之一称为v e r l e t 算法，它的基本思想是对时间点t 处的粒子坐标做泰勒展 开： ，( m r ) ；帅v ( 协辔f 2 + 争+ 。( ) ， 其中“是分子动力学模拟中的时间步长，v ( f ) 为f 时间时粒子的速度，m 为粒子质 量，o ( a t 4 ) 是关于a t 的四阶无穷小量。那么，同样也有 rt - n t ) _ r ( f ) - 雌) f + 等一争+ 。( “4 ) ， 加和上述两式，我们得到 r ( m r ) + ，( r ) _ 2 r ( r ) + 掣+ 。( ) ， 也就是 ，( m r ) ：2 ，( r ) 一，( t ) + 掣 ( 1 - 2 ) 此式的误差仅仅是“4 数量级，合理选择a t ，会使计算变得相当精确。注意， 这里并没有用到速度来计算新的位置坐标。但是速度参数也是可以由运动轨迹中的 坐标计算出来，在此就不详细讲述了。同样比较常用的还有l e a p f r o g 算法1 7 j 及其他 高阶算法【6 1 ，目的都是提高体系的精确性和运算速度。 通过分子动力学模拟，可以得到体系某些量随时间的演化，对某一时间段里的 演化量进行平均，就可以得到体系在这一时间段的时间平均( t u n ea v e r a g e ) ，这也 正是真实实验中观测体系某个量的方法。通过运动轨迹，人们可以计算体系的许多 平均属性，比如平均能量、平均密度，或者一些动态量，比如与扩散效应相关的粒 子均方移动距离随时间的关系等。 1 3 高分子中的分子动力学模拟 分子动力学方法在高分子研究中的应用十分广泛。通过建立高分子材料的分子 模型，结合已知的分子问相互作用，可以模拟模型中所有分予在一定时间里的微观 第3 页共7 6 页 复旦大学硕士论文绪论 运动。这种分子模型可以完全按照高分子真实结构，构造分子链上的所有原子( c 、 h 、n 等) ，通过此方法获得的模拟结果包含的信息非常充足，但是这种模型的一个 明显的缺点是需要引入大量的相互作用，增加了模型的复杂度；另一种常用的分予 模型称为联合原予模型( u n i t e d a t o mm o d e l ) 1 7 ，它把分子粗粒化为原子基团之间 直接的相互作用。已聚乙烯( p e ) 为例，通常把一个c h 2 链节粗粒化为一个相互作 用的基本单元，而整条p e 链则由许多这样的基本单元通过分子键相连。很明显， 人为设计的分了模型，它的相互作用参数与实验测得的参数有差别，不同高分子模 型的分子相瓦作用参数可以参见f r e n k e l 等人所著的 u n d e r s t a n d i n gm o l e c u l a r s i m u l a t i o n s ) ) 一书 6 】o 分子动力学模拟大量的被应用于研究各相同性的高分子模型，如聚乙烯的静态 s - l l l 和动态f 1 2 | ”l 的性质；同时也被应用到其他复杂高分子，如聚苯乙烯的微观动力 学和玻璃化转变1 1 4 - 1 8 】。这些研究工作从热力学和动力学的角度考查高分了体系的微 观运动。同时，也有大量的工作投入到研究高分子材料的形变过程。最早用分子动 力学方法研究高分子的形变过程是b r o w n 和c l a r k e l l 9 1 于1 9 9 1 年研究聚乙烯单链体 系的单轴拉伸形变过程。他们采用单链聚乙烯模型，链长1 0 0 0 ，同时应用了不同的 拉伸速率，最终得到于实验值相近的弹性形变、应力屈服和塑性流动等性质。此后， 分子动力学模拟被更深入的应用在聚乙烯1 2 0 2 3 】、聚苯乙烯( p s ) 2 4 , 2 5 】、聚碳酸酯( p c ) 1 2 5 , 2 6 1 、聚甲基丙烯酸甲酯( p e t ) 2 7 1 等多种高分子模型的形变研究中。通过模拟德 到的体系各项热力学或者动力学属性可以力学【2 8 1 、核磁共振【2 明以及中子散射 3 0 , 3 1 1 等实验得到的数据对比，从而更深层的揭示高分子材料的所具有的独特性质。 分子动力学模拟包含许多微小的时间步长，最终得到体系运动的轨迹。而模拟 所需要的计算机时间与时间步长的设定值和体系的大小密切相关。通常的时间步长 在飞秒( 1 0 。1 5s ) 数量级，这样，在有限的计算机模拟时间范围里，分子动力学可 以模拟从几个飞秒到上百纳秒( 1 f f 9s ) 的运动过程。很显然，这种比较短的时间范 围，使得模拟一些较慢的体系运动过程，变得非常没有效率。比如，高分子链的管 道运动( r e p t a t i o n ) 3 2 j 的驰豫时间与高分子的链长度n 成指数关系( n 3 4 ) ，对于长 度很长的高分子链，蛇行运动的驰豫时间已经超过分子动力学模拟的时间范围，超 过驰豫时间范围的蛇行运动则不能用分子动力学很好的模拟出来。而m o n t ec a r l o 方法则很好的解决了这个问题，因为m o n t ec a r l o 方法一般不引入时间概念，体系 的运动由构象空间( c o n f i g u r a t i o n a ls p a c e ) 里的体系结构互相转化来实现。这种转 化的方法就是构造合适的m o n t ec a r l o 运动( m o v e ) ，而不同的m o n t ec a r l o 运动可 以对应不同的驰豫时间范围，从来很好的模拟各种时间范围内的分子运动。 第4 页共7 6 页 复旦大学硕士论文绪论 1 4m o n t ec a r l o 模拟简介 m o n t ec a r l o 方法是一种基于随机过程的方法，要了解这种方法，我们先要了解 什么叫做重要性抽样( i m p o r t a n c es a m p l i n g ) 。根据经典统计热力学，知道体系的配 分函数，就可以计算一个可观察量a 的系综平均( e n s e m b l ea v e r a g e ) 。假设a 是关 于体系构象函数一的函数，其中n 是体系的粒子组成单元数( 粒子数) ，则( 4 ) 可 以写成： ( 爿) = f d r 丽e x p 可e - f l u 面( r u 研) a ( 厂r n ) ( 1 - 3 ) 其中卢是波尔兹曼因子，它于体系的温度r 成反比卢= 1 ，k f ，u 是体系的构象能， 它也是一个只于体系构象函数r ”相关的函数。设想，如果知道体系的所有构象，就 可以计算每个构象对应的能量，从而得到配分函数，最后通过式( 1 3 ) 计算a 的构象 平均，也就是系综平均。但是，这种思路清晰的方法却严重依赖体系的大小v ，因 为对于多元体系，的数值一般较大，对应这种体系的构象数量通常是非常庞大， 计算每个构象变得很不现实。 要解决这个问题，就要引入重要性抽样的概念。对于多元体系，不同构象所对 应的能量不尽相同，高能量的构象出现的几率较低，对应这类构象的构象数较少， 反之亦然；由于配分函数中波尔兹曼分布的特殊性( 指数关系) ，使得多数高能态构 象对应的波尔兹曼权重几乎等于零，而系综平均仅仅由少数低能态构象主导。因此， 只计算对应波尔兹曼权重较大的构象，忽略小权重构象，对重要性的构象进行平均， 就可以近似的计算体系的系综平均，这就是重要性抽样的基本思路。 m o n t ec a r l o 方法就是基于重要性抽样概念发展而来的一种方法。m e t r o p o l i s 等 人i j 圳最早于1 9 5 3 年创造m o n t ec a r l o 重要性抽样算法，称为m e t r o p o l i ss c h e m e ，它 是最简单的m o n t ec a r l o 算法，但是却代表了最基本的m o n t ec a r l o 方法的所有元素。 现在让我们来讨论如何获得那些对应高波尔兹曼权重的构象。首先建立一个体系利 它的一个确定的构象r ”，用0 ( o l d ) 表示，根据它可以计算出一个非零的波尔兹曼 权重e x p 卜p u ( o ) i ：接着，对这个构象施加一个非常小的尝试性扰动，使体系演 变到一个新的构象r ，用n ( n e w ) 表示，它对应的波尔兹曼权重为e x p i p u ( n ) i ： 将新构象的波尔兹曼权重和旧的做比较，来确定是否要接受这个扰动。如何来确定 这个接受几率石( o 专n ) 呢? 我们可以假设对这两个构象做许多次重复但不相同的 m o n t e c a r l o 试验，通过计数获得处于0 构象和n 构象的总构象数m ( 0 1 和m ( n ) ，可 第5 页共7 6 页 复旦大学硕士论文 以想象，m 是和体系的构象密度p 成正比的。综合考虑体系的构象密度和转化几率， 可以得到如下的细致平衡原理( d e t a i l e db a l a n c ec o n d i t i o n ) ： p ( o ) 万( o n ) = p ( n ) 万( n o ) ， ( 1 4 ) 万( o 斗n ) 其实代表着一种结构的转化控制矩阵，这种转化过程通常分为两个步骤， 第一步是口( o n ) ，它对应一种潜在的m a r k o v 进程嘲：觥a c c ( o - - ) n ) ， 它用来确定是否接受这次尝试，所以有如下关系 石( o _ n ) = 口( o 斗n ) a c c ( o 斗n ) ( 1 - 5 ) 在较为简单的m o n t ec a r l o 算法中，口( o 斗1 3 ) 通常是一个对称矩阵，即 a ( o 畸n ) = 口( n 一0 ) ，但是在众多复杂的m o n t ec a r l o 算法中，这个矩阵是不对称 的，在此就不详述了。如果口( o - - n ) 是对称的，式( 1 - 4 ) 可以改写成 p ( o ) a o o ( o 呻n ) = p ( n ) a c c ( n 寸o ) ( 1 6 ) 进一步变形可以得到 竺! ( 1 2 尘 a c c ( n o )鬻一一m ) 叫。) ( 1 - ，) 注意，代表几率的数值总是小于等于l 的。很多不同的a c c ( o _ n 1 形式都满足式 ( 1 7 ) ，在此我们给出m e t r o p o l i s 形式 。( 。斗n 1 _ je x p p ( n ) - u ( 。) u ( n ) u ( 。) ( 1 - 8 ) 、7 i l i f u ( n ) u ( o ) 知道了接受尝试性运动的几率a c c ( o 斗n ) ，它是一个( 0 ，1 1 区间内的实数，在m o n t e c a r l o 模拟上具体接受运动的方法是调用一个区间为( 0 ，1 1 的伪随机数( 伪随机数由 随机数生成器程序产生，具体请参考口4 1 ) 与接受几率相比，来确定是否接受这次尝 试性运动。可见，m o n t ec a r l o 方法是一种基于随机概念上的方法。 至此，我们简要介绍了m o n t ec a r l o 模拟运动的接受几率问题，其次需要解决 的问题就是如何设计一个合适的m o n t ec a r l o 运动，这种运动代表着体系构象的转 变，它可以是具有物理意义的运动，比如键角的弯曲振动；也可以是没有物理意义 的运动，比如随机调换体系中的无关的两个结构单元的位置。设计的目的是为了有 效的使体系在构象空间里演化。一般说来，新出现一种运动，判断其是否正确的依 据是看它是否满足细致平衡原理1 ，不满足细致平衡原理的运动，会使体系向错误的 1 满足细致平衡原理是用来判断新m o n t e c a r l o 运动是否正确的充分条件，但不是必要条什，请参阅 m a n o u s i o u t h a k i s ，vi ；d e e m ，m ws t r i c td e t a i l e db a l a n c ei su n n e c e s s a r yi nm o n t ec a r l os i m u l a t i o n 明j o u r n a lo f 第6 页共7 6 页 复岜大学硕士论文绪论 方向演化，最终得到错误的结论。要用m o n t ec a r l o 方法计算体系的系综平均，所 采用的m o n t ec a r l o 运动要能够使体系达到各态历经( e r g o d l c i t y ) ，这样才能避免体 系因为运动设计没有效率而停留在一些亚稳态上，导致错误的系综平均结果。 这里我们给出一个最简单的m o n t ec a r l o 运动粒子跳跃( h o p ) 。每次粒了 跳跃只操纵一个粒子发牛位移，假设该粒子的原始坐标是r f 0 1 ，则尝试性运动之后 粒了的新坐标r f n l 可以用下式表示 r ( n ) = r ( o ) + ( 2 r a n f 1 ) - a t e ， ( 1 9 ) 其中r a n f 是一个区间为( 0 ，1 1 的伪随机数，是每次跳跃允许的最大距离，e 是一个 随机的单位向量。可见，粒子的跳跃范围是在一个以原始位置为球心，a 为半径的 球体内。值的选取原则是在获得大尺度运动和获得高的接受几率之间取得平衡， 大的值可以使粒子活动范围增加，但是会降低跳跃的接受几枣；而小的值能够 带来高的接受几率，但却减少了粒子的活动范围，是体系演化变得没有效率。通常 人们选择一个合适的值，使得跳跃运动的接受几率维持在5 0 左右。在选择随机 取向的单位向量e 时，人们通常使用球极坐标，但是通过两个随机的方向角获得的 单位向量并不是真正的随机取向向量，因为这种方法并不能带来向量端点在单位球 表面的随机分布。正确的随机取向单位向量的获得，请参考【6 1 。在通过跳跃运动获 得一个新的坐标之后，计算两个状态下对应的体系能量，然后可以通过上面介绍的 m e t r o p o l i ss c h e m e 来判断这次运动的接受几率，实现m o n t ec a r l o 运动。 1 5 高分子中的m o n t ec a r l o 模拟 高分子中的m o n t ec a r l o 模拟最初是从研究高分子单链的构象统计开始的，通 过模拟高分了链在真空或者溶剂中的构象变化，可以得到高分了链的均方末端距或 者均方旋转半径等口”。最早的m o n t ec a r l o 算法是建立在高分了的格了模型( l a t t i c e m o d e l ) 上的，以最简单的二维格子模型为例，高分子的链节被分布在一个正交的格 了交点上，而交点之间的连线代表着高分子的化学键。通过研究高分子链节在这些 格子上的移动，可以模拟高分子链的运动，或者计算高分子链的结构特征。早期用 格子模型研究高分子链构象的工作有w a l l 等人【3 6 】和r o s e n b l u t h 等人m 1 ，后者是 c o n f i g u r a t i o n a l b i a sm o n t ec a r l o 方法的开端。随着理论和计算机技术的发展，人们 发展简单的格子模型已经不能满足模拟工作的需要，所以有许多高级的格子算法诞 c h e m i c a l p h y s i c s , 1 9 9 01 1 0 f 6 ) ：2 7 5 3 2 7 5 6 第7 页共7 6 页 复咀大学硕十论文绪论 牛，比如c a r m e s i n 等人于1 9 8 8 年提出的键长涨落模型( b o n df l u c t u a t i o nm o d e l ) 3 8 1 ， 把m o n t ec a r l o 模拟的格子模型引入了一个新的阶段。 诚然，格子模型取得了许多重要的成果，它至今仍然是人们常用的一种m o n t e c a r l o 模型，但是它本身的局限性是很显而易见的。格点的引入使得整个体系的构象 数大量减少，以最简单的二维格子模型为倒，从一个节点出发构造一根化学键，它 的取向只有四个方向可以选择，同时键长的变化也受到不同程度的阻碍。这与真实 物质的运动是截然不同的，真实的高分子结构处于一个连续的三维空间中，它所能 达到的构象数量大大超过格子模型可以模拟的数量，所以说格了模型只是一种近似 的模拟方法，从统计意义上趋近于真实结果。于是，人们想到了使用连续模型( o f f l a t t i c em o d e l ) ，这种建立在连续空间中的模型，从一个相对真实的基础上研究高分 子链，是比较合适的。v a c a t e l l o 等人【3 9 于1 9 7 9 年第一次用m o n t ec a r l o 方法在连续 空间里研究c 3 0 熔体的特征，同时提出了高分子链的管道运动模型( r e p t a t i o n ) 。 连续模型增加了体系的自由度，使模型变得更加真实可信，但是也一定程度上也增 加模型的复杂性和可操作性。 随着高分子研究的进展，人们不能满足于短链长高分子模型研究的结果开始 从单链体系拓展到多链体系，从短链长体系发展到长链长体系，从研究单链构象转 化为研究多链相互作用等等。但这些新的思想无疑给计算机模拟的效率带来的挑战， 在当时允许的计算机硬件条件基础上，人们不得不对已有的m o n t ec a r l o 运动进行 改进或者提出新的高效的算法。t h e o d o r o u 等人f 4o 4 1 1 提出了高效的链内联合旋转运 动( c o n c e r t e dr o t a t i o n ) 和尾端桥接运动( e n d b r i d g i n g ) ，最初用来研究聚乙烯多链 体系的熔体特征。这些新型运动从大尺度上改变链的整体构象，使得体系在构象空 间里能够高效的演化，尤其是尾端桥接运动，它改变高分子多链体系的连通性，从 某种意义上说不具备物理意义，但其大尺度的改变体系构象增加了统计平均数据的 可靠性。这些运动被用在我们的研究之中，具体细节详见后文。 至此可见，m o n t ec a r l o 方法相对于分子动力学方法而言，在演化体系结构的方 法上存在着差异：分子动力学依靠的是解决牛顿运动方程，从而获得体系的运动轨 迹；而m o n t ec a r l o 方法是依靠m o n t ec a r l o 运动搜索体系的构象空间里的低能态构 象，通过构象平均算出体系的一些平均量。分子动力学很难解决的长链高分子的演 化效率问题，在m o n t ec a r l o 模拟中却可以通过设计高效m o n t ec a r l o 运动得以改善。 分子动力学顾名思义，可以用来研究高分子的动力学特征，在高分了材料形变方面 有大量的研究( 见前文所述) ；而m o n t ec a r l o 方法并没有引入时间的概念，它是一 个用来搜索构象空间的静态方法，所以m o n t ec a r l o 方法用于高分予材料形变过程 第8 页共7 6 页 复旦大学硕十论文 的研究并不是很多。杨玉良教授在其高分子科学中的m o n t e c a r l o 方法一书中 也曾指出，只要m o n t ec a r l o 算法驱使的体系状态转移过程与实际物理过程十分逼 近，那么m o n t ec a r l o 动力学可以被认为是真实反应动力学的客观描写【3 “。w i t t k o p 等人【4 2 】用二维和三维格子模型研究了单链高分子链体系的形变，他们发现在不同的 外力条件下模拟得到链末端距与理论估计值非常吻合。o s p i n a 等人1 4 3 4 4 悃二维格子 模型模拟了线形聚乙烯的形变初期过程，发现其应力应变关系与实验【4 5j 得到的数据 是相似的。h o l z l 等人【4 6 ，4 ”用键长涨落模型考查了多链高分子体系形变时候的应力 应变关系受链长、链在拉伸短接枝密度以及玻璃化转变温度等因素的影响。b o y c e 等人【3 】用三维连续模型模拟了聚乙烯网络的掩伸和压缩形变，获得了和实验非常相 似的应力应变关系，同时也阐述了它的温度和拉伸速率依赖性。 1 6 本文研究的目的和内容 分子动力学方法用于高分了材料形变的研究已经取得很大的成功，但它在解决 一些驰豫时间较长的过程，比如长链高分子体系的管道运动时，有着本质的缺陷： 而m o n t ec a r l o 方法在运动的设计上具有很强的灵活性，使得其演化体系的效率高 于分子动力学方法，因此，本文试图用m o n t ec a r l o 方法研究高分子的形变动力学。 然而，无论是格了模型还是连续模型，上述已经发表的m o n t ec a r l o 方法用于 高分子材料形变的研究都是基于简单的粒了跳跃运动，这种运动不能很有效的使体 系在构象空间里演化，即并没有非常好的利用m o n t ec a r l o 方法的优越性。使用这 种运动获得模拟结果，要么需要较长的模拟时间，要么不能使体系发生大的构象改 变，比如模拟大应变形变，所以通过这种方法获得的结果，与用分子动力学方法， 并不具备很大的独特性。 t h e o d o r o u 等人【4 0 4 1 】开发的新型m o n t ec a r l o 运动，即联合旋转运动和尾端桥接 运动，在研究高分子熔体的平衡过程上取得了巨大的成功。这些运动比分子动力学 方法或者其他低效m o n t ec a r l o 运动，能更有效的使高分子多链体系演化至平衡态； 甚至对于分子动力学无法企及的超长链高分子( 单链平均聚合度为5 0 0 ) 体系，使 用尾端桥接运动也能使其有效的脱离初始结构，让这样的体系达到各态历经成为可 能。 因此，本文研究的出发点是将这些高效的m o n t ec a r l o 运动用于线形多链聚乙 烯连续模型的形变，研究其形变时候的应力应变规律和结构变化因素。在此之前， 有两个主要的问题需要解决：第一，为了提高模拟的效率，许多研究中都对高分子 第9 页7 6 页 复旦大学硕士论文 体系施加了限制条件，比如使用刚性键长和固定键角，这样可以减少体系的自由度， 降低模型的复杂度；同时，因为键长振动和键角振动一般属于高频振动，他们通常 会使计算结果发生很大的涨落，所以使用固定键长键角可以减少输出数据的涨落幅 度。因此，本文对聚乙烯模型没有限制条件和只施加刚性键长限制条件下的形变过 程进行研究，来比较两者的区别。第二，虽然刚性键长的引入简化了模拟过程，却 使体系内应力的计算变得复杂，因为这些限制条件对体系总应力的贡献的计算相对 比较复杂，所以我们先考查了不同限制条件下，体系内应力的计算方法。 第二章中我们将讨论不同限制条件下体系内应力的计算方法，给出后面计算中 将要用到的具体内应力表达形式。从考虑多链聚乙烯连续模型的复杂度出发，我们 循序渐进的把工作分为接下来三章：第三章中我们先研究了一个最简单的单原子体 系的平衡过程，从中我们考查简单m o n t e c a r l o 跳跃运动在简单体系中的表现，同 时将检查基本的内应力计算方法：第四章中我们将建立双原子分子模型，由于分子 键的引入，我们可以考查刚性键长和弹性键长两种条件下，体系内应力的计算，尤 其注重判断刚性键长限制条件下计算的正确性：第五章中我们将前面获得的正确结 果推广到聚乙烯体系，我们将系统的阐述从分了模型的建立、平衡，知道最终形变 的全部过程，最后给出我们得到的初步结果。最后，第六章中我们给出本文的结论 和对将来工作的展望。 第1 0 页共7 6 页 复吐大学硕士论文 参考文献 1 】m e l i c k , h gh ；g o v a e r t ，l e ；m e i j e r , h e h o nt h eo r i g i no f s t r a i nh a r d e n i n g i ng l a s s yp o l y m e r s j 】p o l y m e r ，2 0 0 3 ，4 4 ( 8 ) ：2 4 9 3 2 5 0 2 ( 2 】m e l i c k , h gh ；g o v a e r t ，l e ；m e i j e r , h e h l o c a l i s a t i o np h e n o m e n a i ng l a s s y p o l y m e r s ：i n f l u e n c eo f t h e r m a la n dm e c h a n i c a lh i s t o r y 【j p o l y m e r ，2 0 0 3 ，4 4 ( 1 2 ) ： 3 5 7 9 3 5 9 1 3 】c h u i ，c ；b o y c e ，m c m o n t ec a r l om o d e l i n go fa m o r p h o u sp o l y m e rd e f o r m a t i o n ： e v o l u t i o no fs t r e s sw i t hs t r a i n 【j 】m a c r o m o l e c u l e s ，1 9 9 9 ，3 2 ( 11 ) ：3 7 9 5 3 8 0 8 【4 a 仃u d a ，e m ；b o y c e ，m c e v o l u t i o no fp l a s t i ca n i s o t r o p yi na m o r p h o u s p o l y m e r sd u r i n gf i n i t es t r a i n i n g 【j i n t e r n a t i o n a l j o u r n a lo f p l a s t i c i t y ，1 9 9 3 ，9 ( 6 ) ： 6 9 7 7 2 0 5 】m e l i c k , h gh d e f o r m a t i o na n df a i l u r eo f p o l y m e rg l a s s e s 【d 】e i n d h o v e n ： t e c h n i s c h eu n i v e r s i t e i te i n d h o v e n ，2 0 0 2 6 】f r e n k e l ，d ；s m i t ，b u n d e r s t a n d i n gm o l e c u l a rs i m u l a t i o n f i o ma l g o r i t