《蒙特卡洛方法在高分子材料中的应用ppt课件.ppt》由会员分享，可在线阅读，更多相关《蒙特卡洛方法在高分子材料中的应用ppt课件.ppt（56页珍藏版）》请在三一办公上搜索。
1、第六章 高分子科学中的Monte Carlo 方法,Monte Carlo方法一个十分独特的名字,Monte-Carlo, Monaco,Monte Carlo原为地中海沿岸Monaco的一个城市的地名，气候温和，景色怡人，人口不到一万，是世界闻名的大赌场。将Monte Carlo作为一种计算方法的命名固然已经赋予了新的内容。然而，顾名思义， Monte Carlo方法的随机抽样特征在它的命名上得到了反映。,MC方法的发展归功于核武器早期工作期间Los Alamos（美国国家实验室中子散射研究中心）的一批科学家。von Neumann, Metropolis, Ulam和Kahn等人在电子计算
2、机上对中子行为进行随机抽样模拟，通过对大量中子行为的观察推断出所要求算的参数。Los Alamos小组的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了MC在各种问题中的应用。学术界一般将Metropolis和Ulam在1949年发表的论文作为Monte Carlo方法诞生的标志。,6.1 Monte Carlo方法的基本思想,Monte Carlo方法在数学上称其为随机模拟(random simulation)方法、随机抽样(random sampling)技术或统计试验(statistical testing)方法它的最基本思想是：为了求解数学、物理及化学等问题，建立一个概率模型或随机过程
3、，使它的参数等于问题的解；当所解的问题本身属随机性问题时，则可采用直接模拟法，即根据实际物理情况的概率法则来构造Monte Carlo模型；然后通过对模型或过程的观察抽样试验来计算所求参数的统计特征，最后给出所求解的近似值。在高分子科学中的Monte Carlo模拟主要采用直接模拟方法。,Monte Carlo方法的突出特点是，它的解是由试验得到的，而不是计算出来的。其程序结构简单，解题时受问题条件限制的影响较小，具有广泛的适应性。但不能解决精确度要求很高的问题。蒙特卡洛方法需要大量的随机数，计算量很大，人工计算需耗费大量的时间，利用计算机可大大减少计算时间，增加试验次数以提高计算精度，因此，
4、蒙特卡洛方法的广泛应用与计算机技术的发展是不可分割的。,设所要求的量x是随机变量的数学期望E()，那么用Monte Carlo方法来近似确定x的方法是对进行N次重复抽样，产生相互独立的值的序列l， 2， N，并计算其算术平均值：,根据Kolmogorov的大数定理则有：,即当N充分大时， 成立的概率等于1，亦即可以用 作为所求量x的估算值。,例6-1 用统计试验方法求圆周率考虑边长为1的正方形，以其一角为圆心和边长为半径，在正方形内画一条14圆弧，如图所示。,在正方形内等概率地产生n个随机点(xi，yi)，i = l，2，3，n，设n个随机点中有k个点落在四分之一圆弧内，显然，当n 时有以下关
5、系成立：,因而，圆周率的估值为：,判断随机点(xi，yi)是否位于圆内的判别式为：,用一对(0，1)随机数Ul，U2分别模拟随机变量的取值xi和yi，当 时，则计数器k值增1。这个判别式就是蒙特卡洛方法的概率模型。当试验次数n足够大时，所得的估值的精度也随之提高。,例6-2. 蒲丰氏问题,Comte de Buffon (1707-1788) French Needle experiment, 1777,Buffon投针问题：平面上画很多平行线，间距为a。向此平面投掷长为l（ la）的针, 求此针与任一平行线相交的概率p。,可以证明,求出值,其中为投计次数，n为针与平行线相交次数。这就是古典概
6、率论中著名的蒲丰氏问题。,一些人进行了实验，其结果列于下表 ：,6.2 Monte Carlo方法与高分子科学,Monte Carlo模拟与高分子科学结下了不解之缘是由于高分子科学本身的特点所决定的，因为在高分子科学中存在着大量可供进行Monte Carlo直接模拟的随机性问题。如：由于聚合反应本身的随机性特点，高分子系综内各个成员之间存在着与其生成机理密切相关的特定分布，即体系中所生成的高分子链并非具有相同的分子量，而是存在着所谓的分子量分布问题；在多元聚合中，多元共聚物不仅具有分子量分布，而且导致了不同种单元在高分子链上的排列问题，即所谓的序列分布；在多官能团的聚合反应中的支化和凝胶化问题
7、；高分子链的热降解和辐射降解等等，无一不是随机性问题。,Monte Carlo方法在现代高分子科学中的应用主要具有以下特征:,由于高分子凝聚态物理的发展，高分子体系的Monte Carlo研究从对单链的研究转向对高浓度多链体系的研究。由静态平衡态问题向动态和非平衡态问题发展也是当前高分子Monte Carlo模拟的重要特征。高分子链的分子运动学，尤其是高浓度多链体系的分子运动问题是当前研究的重要方面。人们对共混和嵌段共聚物的界面、高分子和液晶的界面、高分子链的吸附、晶态和非晶态的界面性质和相互扩散问题开展了Monte Carlo模拟研究。 高分子Monte Carlo方法的新算法也是值得研究的
8、。,6.3 随机数与伪随机数,产生均匀分布随机数的方法可以采用物理方法和数学方法。最简单的产生随机数的物理方法是掷骰子游戏；采用电学噪声的变化也可产生随机数。但物理方法产生随机数的“费用”很高，且速度慢。因此，实际应用的随机数一般均在计算机上采用数学方法来产生。用数学方法产生的随机数一般均采用某种确定性的表达式来实现，因此其并非真正的随机，故通常称其为“伪随机数”。用数学方法产生伪随机数的优点是因为它借助于迭代公式，所以特别适合于计算机。而且其产生的速度快、费用低。目前，多数的计算机均附带有“随机数发生器”。,用数学迭代方法产生的随机数存在两个问题：,1、递推公式和初始值a1、a2、ak确定后
9、，整个随机数序列便被唯一确定下来。即任意一个随机数被前面的随机数唯一确定了，不满足随机数相互独立的要求。2、既然随机数序列是用递推公式确定的，而在计算机上所能表示的0，1上的数又是有限多的，因此这样的随机数序列就不可能不出现重复地无限继续下去。这种随机数序列出现周期性的循环现象是与随机数的要求相矛盾的。,对第一个问题不能从本质上改变，但只要递推公式选得好随机数的相互独立性是可近似满足；第二个问题，则不是本质的，因为用Monte Carlo方法解任何问题时，所用随机数个数总是有限的，只要保证不超过伪随机数序列出现循环现象的长度即可。,用数学迭代方法产生随机数均存在周期现象，随着迭代过程的不同，其
10、效果也各不相同。一般满足下列要求的产生方法才可被认为是好的： (1)随机性和统计独立性要好； (2)容易在计算机上实现； (3)省时，存贮量小； (4)伪随机数的周期长。,乘同余法,乘同余法由Lehmer首先提出。由于采用乘同余法具有在计算机上容易实现、快速等优点，因此乘同余法已被广泛采用。乘同余法的迭代公式为，,一个简单的例子,上面的例子中，第一个随机数生成器的周期长度是 10，而后两个的周期长度只有它的一半。我们自然希望随机数的周期越长越好，这样得到的分布就更接近于真实的均匀分布。,表：乘同余法的参数及周期,Monte Carlo方法的核心就是随机数的使用，因此计算机模拟结果的优劣将强烈地
11、依赖于伪随机数的质量。,伪随机数的均匀性伪随机数的独立性,对于已经产生的随机数质量的检验主要是：,伪随机数的均匀性检验可用xn的矩来判别，均匀性好的随机数序列在N时应满足下列要求：,一阶矩,二阶矩,三阶矩,四阶矩,伪随机数独立性检验一般采用2检验。,随机变量的抽样：前面讨论了0，1均匀分布的伪随机数的产生，然而在实际应用中概率分布的形式是多种多样的。,并满足：,产生0，1随机数r，如果条件满足，则认为事件Ai发生。,一、从随机事件中抽样：假设随机事件的出现概率分别为Pi (i1，2，n)。为了对随机事件Ai进行抽样，首先需构造累积概率：,例6-3. 掷骰子点数的抽样,掷骰子点数X=n的概率为：
12、选取随机数，如 则在等概率的情况下，可使用如下更简单的方法：其中 表示取整数。,二、连续型分布的抽样：连续型分布的一般形式如下：这里f(t)为分布的概率密度函数。如果分布函数的反函数存在，则连续型分布的一般抽样方法是通过其反函数直接抽样：,这里r是0，1均匀分布的随机数，F-1为F(x)的反函数。,在a,b上均匀分布的分布函数为：,例6-4. 在a,b上均匀分布的抽样,其抽样方法为：,这里r是0,1区间均匀分布的随机数。,Monte Carlo方法的估值精度与试验次数N的平方根成反比，若精度提高10倍，则试验次数N要增加100倍。收敛速度慢是蒙特卡洛方法的主要缺点。蒙特卡洛方法的精度估算有概率
13、性质，它并不断言精度一定好于 ，而只是表明，所算精度以接近于1的概率不超过某一界限，这是蒙特卡洛方法与其它确定性误差计算的根本区别之处。,例6-5：中子扩散问题,原子核反应堆的壁是铅制的，对中子起屏蔽作用。中子从反应堆内侧进入壁内与铅原子发生碰撞。求出穿透铅壁中子数的百分比，被吸收入铅壁中子数的百分比，以及重新返回反应堆中子数的百分比。,解：设壁厚为常量3d，中子是垂直进入壁内的，并设每个中子在壁内每次走过d（平均自由程）才与铅原子碰撞，碰撞后以随机的方向弹射，再走过d的距离，和第二个铅原子碰撞，如此继续下去。最后，有三种情况（1）中子穿透铅壁；（2）被铅壁吸收（假定经过8次碰撞后，没有穿透或
14、返回，则认为被吸收；（3）重新返回反应堆。,现在研究对中子运动的模拟：假设一个中子在壁内处于与壁内侧距离为x的位置上与铅原子碰撞，然后以角的方向弹射，那么是0，2之间的均匀分布的随机数。中子经过弹射后，与壁内侧的距离x变为：x + d cos(2 y)若（1）x 3d 则中子穿透铅壁 （2）x 0 则中子返回反应堆 （3）0 x 3d 则继续下一次碰撞，重复这个过程直至中子脱离铅壁或8次碰撞后被吸收为止。,对5000个中子进行模拟的结果为：穿透 26.3%；吸收 22% ；返回 51.7%,6.4 Monte Carlo方法在聚合物研究中的应用示例,聚乙烯分子结构的模拟共聚物序列分布的Mont
15、e Carlo算法高分子无规行走(random walks)链的模拟研究高浓度多链体系动力学的“空格扩散算法”,聚乙烯分子结构的模拟,聚乙烯分子是由重复的 (CH2 )单体组成的长分子链，由于碳原子有四个共价键，其空间构型如图所示。,Flory根据统计力学理论，导出柔性分子链的非晶态结构取无规线团的构象，各高分子链之间可以相互贯通，它们可以缠结。该无规线团模型在70年代利用中子小角散射技术得到了证实。,聚乙烯分子的空间构型在平面上的投影，可以近似地看成如图所示的结构。模拟程序先定义八个方向，并给出每个方向对应的数值，如图。当分子链段方向为3时，其后面分子链的可能取向方向为2、3或4，它们在聚乙
16、烯中是等概率的。至于下面分子链向哪个方向运动，可由计算机产生的随机数来决定。这样就可模拟出聚乙烯的无规线团状分子结构。,共聚物序列分布的Monte Carlo算法,共聚反应的Monte Carlo研究开展得较早，所涉及的主要问题是组成和序列分布问题，其主要目的是通过共聚产物的序列分布来获得单体的活性比和鉴别不同的反应机理。共聚反应的Monte Carlo算法比较简单，因此我们只是简要地介绍其基本算法。,具有末端效应两元共聚反应：末端效应是指只有端点上的单体单元对聚合反应的速率常数有影响。对于两元共聚反应的四种增长反应可记为：,由此还可定义活性比，,假定，各速率常数与链长无关，而且引发和终止过程
17、可忽略(一般当高分子的链长很长时均可认为引发和终止过程的影响可忽略)，则由M1*到M1*的转变概率为：,这里M1表示投料浓度，而由M1*转变为M2*的转变概率为：,相应地有：,Monte Carlo模拟程序可由如下步骤构成:,(1) 设增长链的第一个单元为Mi(i1，2)，根据M1、M2、M1* 、M2* 的浓度可计算活性比r1，r2和转变概率pij；(2) 产生一个单位区间内均匀分布的随机数。(3) 因pi1+ pi2l，故若 pi1则在增长链上加上一个Ml单体，并认为其生成了M1* 。若认为单体Ml和M2的浓度在增长过程中一直保持恒定，则转回步骤(2)继续进行模拟。但若认为单体浓度是可变的
18、，则由于Mi单体消耗了一个分子故必须重新计算浓度Mi(il，2)，然后再回到步骤(1)继续进行模拟。,(4) 若上式不满足，即pi1，则表明发生M1*到M2*的增长反应。因此，在增长链上加上一个M2单体，并认为增长链的端基己转变为M2* 。对于恒定浓度的情况，转回步骤(2)，而对于非恒定浓度的情况，则计算变化后的浓度再转回步骤(1)；(5) 上述步骤一直重复，直至达到所需的链长或所需的单体转化率。在模拟过程中可统计各感兴趣的量，如链上Ml和M2单体的组成和序列分布等。必须指出，上述过程只模拟了一根链的情况。为了获得较高的统计精度，可重复多条链后进行平均。,采用上述Monte Carlo算法，M
19、otoc等模拟了丙烯酸甲酯(M1)和氯丙烯(M2)的共聚反应。他们设r1=0.08，r2=5.1，模拟所得的共聚物组成和三元组的百分数与实验值的比较见表。结果表明，对于该体系可近似地认为只存在末端效应。,Monte Carlo模拟结果与实验结果的比较,高分子无规行走(random walks)链的模拟,无规行走(random walks，简称RW)链模型用来研究柔性高分子链在稀溶液中的大尺度性质。RW链即所谓的自由连接链(freely jointed chain，简称FJC)。其基本特征是：链中两相邻键的夹角(键角)可任意选择，每个键的内旋转角也可任意取值，链中非直接键接的链单元与链单元之间不
20、存在任何相互作用。,以计算均方末端距R2为例，RW链的简单抽样法计算过程如下： (1)键向量r1的起点放在坐标原点，并令k1； (2)对于其后的键向量通过产生在半径为a(通常令a1)的球面上均匀分布的随机点，并以其作为rk向量的终点以及rk+1向量的起点； (3)计算末端距向量RkRk-1+ rk； (4)如果k=n，则将令Rk=R，并求R2，如果kn，则将k+l替代k，并返回到(2)。对于RW链的最基本特征是： RW n高分子物理： RW = n l2,格子链模型,格子链(lattice chain)模型的基本做法是将空间离散化，即链单元只能取空间中某些人为规定的格点(lattice sit
21、e)。显然，格子链在细节上与真实链有较大的差别，但高分子链的许多统计性质(大尺度行为)并不依赖于链模型的细节。,方格子模型是把空间离散化为一个立方点阵，即链单元的空间坐标只能在这个点阵空间所定义的格点上取值。为了便于计算，通常格子的边长取为1。因其空间维数的不同，人们给予方格子链以不同的名称。在两维空间里，人们一般称其为方格子链(square lattice chain)；在三维空间中，称其为立方格子链(cubic lattice chain)。,一、RW链的抽样：,采用直接抽样法生成RW链的方法十分简单，以两维方格子链为例，可由如下几个步骤构成：,(1) 将第一个链节固定在坐标原点上，并设格
22、子的边 长为1； (2) 产生(0，3)整数序列的随机数； (3) 由随机数的数值按预先规定的规则，决定下一个链节所达的格点。设有一链单元的坐标为(x，y)，则由伪随机数的具体数值来决定下一个链单元的坐标位置；,当随机数为0时，x+1x ；当随机数为l时， y+1y ；当随机数为2时，x-1x ；当随机数为3时，y-1 y ；,(4) 重复(2)，(3)直至所需链长n，记录所需的结果，诸如最后一个链单元的坐标位置rn等；(5) 重复(1)(4)直至达到所需的高分子链的分子数(样本容量)M，由所生成的链的样本，可计算链构象统计的特征量，如均方末端距。,这里， 为样本中第l条链的末端距平方，,研究
23、高浓度多链体系动力学的“空格扩散算法”,为了使得格子链模型能更为有效地推广到研究高浓度直至熔体”的高分子体系的动力学问题，陆建明和杨玉良在Larson等提出的键长涨落模型的基础上提出了适合于研究高浓度多链体系动力学的“空格扩散算法”。,按照键长涨落模型，模型链的键长允许取两个数值，即方格子的边长(一般取为1)和格子的对角线(边长的2倍)。由于键长的可涨落性，因此每个格点的配位数分别为8(两维)和18(三维)。,下图给出了典型的微松弛模式和禁阻运动模式。考虑到主要想模拟高浓度多链体系，因此采用空格作为算法的运动主体，具体算法可归结为以下几个步骤。,在具有周期边界条件的元胞中规则地按所需浓度排入所
24、需链长的链，设链长为n，链的总数为N条，而少量的空格在排布中尽可能分布均匀。本文元胞大小取为LXLY=4444，高分子链长取为n=21，高分子数为N=88。链所占的格子分数(浓度)为0.9545 周期边界条件：即某个高分子链单元离开元胞的某一边界时，则在元胞中与该边界相对的边界上有一链单元进入元胞。,由于数据统计的需要和节省计算耗时，不仅要记录第l条链的第m链单元的坐标位置M(l,m)，而且还需知道某一坐标点P(i,j)为哪条链的第几个单元所占据。为此采用编码技术，以二维情况为例，即定义： M(l，m)i 100+ j P(i，j)l 100+ mM(l，m)和P(i，j)的关系是显而易见的。
25、设P(10，10)1014，其表明坐标点(10，10)被编号为10的链上的第14号单元所占据。由此可知M(10，14)1010。,随机地选择一个空格，并在空格周围的配位格点上随机地选择一个链单元，并尝试与空格交换位置；若交换后与该链单元连接的键的键长不大于2，并且不出现禁阻情况，则运动接受，若这个交换不能接受，则放弃交换，但体系再参加一次统计，这样的交换直至体系达到平衡。,(3) 对平衡的体系仍然执行步骤(2)，并统计相应的静态和动态性质。在动态问题的模拟中，时间的取法采用一个Monte Carlo cycle作为时间单位，更具体地说即在(2)中的尝试次数(空格选择次数)等于nN时，时间取为1。,二维高分子动态模拟,表中的1，2，3三点是初始单纯形的三次实验，其中1点为最坏点，2点为次坏点，3点为最好点（转化率最高）。舍去最坏点3。求出重心点后，求出1点的反射点4点。因此，第四次实验结果所得转化率为97%，明显高于2，3点。所以将反射点4点进行扩展到5点(取为1.618)，其转化率为81.5%，此转化率小于3点，扩展失败。因而2，3，4三点构成新的单纯形。依次调优到第6点，其结果为97.5%，第六次实验即得到最佳实验条件。,
