（物理化学专业论文）二元混合物流体相平衡吉布斯系综模拟.pdf

重庆大学硕士学位论文 中文摘要 摘要 真实流体和复杂模型流体的相行为对化工过程设计和研究具有重要的意义。 流体体系的相行为是流体分子基本的相互作用的体现，计算机模拟是描述体系相 行为的行之有效的方法，因为计算机模拟不受实验条件的限制。当今所关注的相 对简单和计算量不是很大的模拟集中于一个称为“吉布斯系综”方法的技术。尽 管吉布斯系综方法不能提供很高的精度，而且不能应用于一些重要体系的模拟， 但是它现在常用于获得单组分和混合物流体相行为的数据。这种方法应用非常广 泛的原因是它的计算方法简单，很容易描述和编程，并且有直观的物理基础。它 不需要使用者知道太多关于模拟体系的相图和其它相行为的信息。 吉布斯系综m o n t ec a r l o 模拟是在两个微型区域进行的，这两个区域代表模 拟体系的两个相，因为采用周期边界条件，每个区域都被它自身镜像包围着，区 域间没有界面作用。指定温度后，通过区域内粒子的移动以达到内部的平衡，改 变两个区域体积以满足压力的平衡，在两区域间进行粒子交换使各组分的化学势 相等，使模拟体系逐渐达到汽液相平衡。 作者应用吉布斯系综m o n t ec a r l o 方法对极性分子和非线性分子的纯组分体 系和混合物体系进行了模拟计算。 以库仑力作用势能和l - j 势能相结合的模型来表示分子间的相互作用，采用 n v t 吉布斯系综m o n t ec a r l o 方法模拟了二氧化碳、丙烷、正丁烷、正戊烷等纯 组分流体。以e x p 一6 势能和库仑力作用势能相结合的模型来描述二元混合物体 系中分子间的相互作用，采用n p t 吉布斯系综m o n t ec a r l o 方法模拟不同温度下 的二氧化碳和丙烷，丙烷和正戊烷，正丁烷和正戊烷等二元混合物体系。模拟中 分别采用了极性分子和非线性分子模型，根据文献值给定能量参数、直径参数等， 调整分子移动、体积变化和分子交换的比率，经过多次模拟计算得到了上述各体 系的气、液组成和密度等相平衡数据，并作出了汽液平衡相图，模拟结果在很宽 的压力范围内和文献中的实验值吻合较好。 关键词：吉布斯系综，m o n t ec a r l o 模拟，汽液相平衡，纯组分，二元混合物 重庆大学硕士学位论文英文摘要 a b s t r a c t p h a s et r a n s i t i o n so fr e a la n dm o d e l c o m p l e x f l u i d sa r eo f s i g n i f i c a n ts c i e n t i f i ca n d t e c h n o l o g i c a li n t e r e s t f r o mam o d e l i n gp o i n t o fv i e w ，t h ep h a s eb e h a v i o ro fa s i m p l i f i e ds y s t e mc o n t a i nt h ee s s e n t i a li n t e r a c t i o n so fi n t e r e s t c o m p u t e rs i m u l a t i o n s a r ean a t u r a lc h o i c ef o rt h i sp u r p o s e ，a s t h e ya l l o wc o m p l e t ef r e e d o m i ns p e c i f i c a t i o n o ft h ep o t e n t i a lm o d e l t h ef o c u so ft h ep r e s e n tc o n t r i b u t i o ni so na r e l a t i v e l ys i m p l e a n d c o m p u t a t i o n a l l yi n e x p e n s i v e s e to f t e c h n i q u e sk n o w nc o l l e c t i v e l yu n d e r t h en a m e o f ”g i b b se n s e m b l e ”m e t h o d w h i l et h eg i b b se n s e m b l ed o e sn o tn e c e s s a r i l yp r o v i d e a c c u r a t ed a t aa n di sn o ta p p l i c a b l et om a n yi m p o r t a n tc l a s s e so fs y s t e m s ，i ti sn o w c o m m o n l yu s e df o ro b t a i n i n gt h ep h a s eb e h a v i o ro ff l u i d sa n dm i x t u r e s t h em a i n r e a s o nf o rt h i sw i d e s p r e a du s ei sp r o b a b l yt h es i m p l i c i t yo ft h em e t h o d i np a r t i c u l a r ， t h em e t h o di se a s yt od e s c r i b ea n d p r o g r a m ，a n dh a s a l li n t u i t i v ep h y s i c a lb a s i s i ta l s o p l a c e sm i n i m a ld e m a n d so ni t s u s e ri nt e r m so fi n f o r m a t i o no nt h ea p p r o x i m a t e l o c a t i o no ft h ep h a s ed i a g r a mo ro t h e ra s p e c t so ft h eb e h a v i o ro fan e w s y s t e mu n d e r s t u d y t h eg i b b se n s e m b l em e t h o di s t h et e c h n i q u eo fc h o i c ef o rd e t e r m i n i n gt h e p h a s ed i a g r a m so f f l u i dm i x t u r e s g i b b se n s e m b l es i m u l a t i o n sa r ep e r f o r m e di nt w om i c r o s c o p i cr e g i o n sw i t h i nt h e b u l kp h a s e s 。a w a yf r o mt h ei n t e r f a c e e a c hr e g i o ni ss i m u l a t e dw i t h i ns t a n d a r d p e r i o d i cb o u n d a r yc o n d i t i o n sw i t hi m a g e so f i t s e l f t h et h e r m o d y n a m i c r e q u i r e m e n t s f o rp h a s ec o e x i s t e n c ea r et h a te a c hr e g i o ns h o u l db ei ni n t e r n a le q u i l i b r i u m ，a n dt h a t t e m p e r a t u r e ，p r e s s u r ea n dt h ec h e m i c a lp o t e n t i a l so fa l lc o m p o n e n t ss h o u l db et h e s a m ei nt h et w o r e g i o n s s y s t e mt e m p e r a t u r ei nm o n t e c a r l os i m u l a t i o n si ss p e c i f i e d i na d v a n c e t h er e m a i n i n gt h r e ec o n d i t i o n sa r er e s p e c t i v e l ys a t i s f i e db yp e r f o r m i n g t h r e et y p e so fm o n t ec a r l o ”m o v e s ”，d i s p l a c e m e n t so f p a r t i c l e sw i t h i ne a c hr e g i o n ( t o s a r i s f yi n t e m a le q u i l i b r i u m ) ，f l u c t u a t i o n si nt h ev o l u m eo f t h et w or e g i o n s ( t os a r i s f y e q u a l i t yo fp r e s s u r e s ) a n d t r a n s f e r so f p a r t i c l e sb e t w e e nr e g i o n s ( t os a t i s f ye q u a l i t yo f c h e m i c a lp o t e n t i a l so fa l lc o m p o n e n t s ) t h eg i b b se n s e m b l ep r o g r a m sa r eu s e db yt h ea u t h o rt oc a l c u l a t et h ep u r ea n d b i n a r ym i x t u r eo f n o n l i n e a rm o l e c u l a ra n d p o l a r m o l e c u l a r t h ec o n s t a n t - v o l u m eg e m cm e t h o di s a p p l i e d t oc a l c u l a t et h e v a p o r - l i q u i d c o e x i s t e n c ep r o p e r t i e sf o rt h ep u r ec o m p o n e n tc a r b o nd i o x i d e ，p r o p a n e ，n b u t a n ea n d n p e n t a n e t h ec o n s t a n t p r e s s u r e g e m cm e t h o di s a p p l i e d t os i m u l a t e b i n a r y i i 重庆大学硕士学位论文 英文摘要 m i x t u r e so f c a r b o n d i o x i d e p r o p a n e ，p r o p m a e n - b u t a n e ，a n dn - b u t a n e n p e n t a n e s y s t e m s l e r m a r d j o n e sc o m b i n e dw i t hc o u l o m b i cp o t e n t i a lm o d e lw a su t i l i z e dt o d e s c r i b et h ei n t e r a c t i o n sb e t w e e nu n i t e dm o l e c u l ef o r p u r ec o m p o n e n t c a r b o nd i o x i d e p r o p a n e ，n - b u t a n ea n dn p e n t a n es y s t e m s ，w h i l eb u c k i n g h 蹦e x p - 6a n dc o u l o m b i c p o t e n t i a lm o d e l sw e r eu s e df o rb i n a r ym i x t u r e sc a r b o n d i o x i d e p r o p a n e ，p r o p a n e n b u t a n e a n dn b u t a n e n - p e n t a n es y s t e m s t h e nt h el i q u i d v a p o rp h a s ec o e x i s t e n c e c u r v e st h a ta r ep l o t t e da c c o r d i n gt ot h es i m u l a t i o nr e s u l t sr e s p e c t i v e l y d e n s i t i e so fb o t h v a p o rp h a s ea n dl i q u i dp h a s e ，v a p o rp r e s s u r e s f o r p u r e c o m p o n e n t sw e r es i m u l a t e d c o m p o s i t i o n p r e s s u r ed i a g r a m sa n dd e n s i t i e sf o rv a p o r p h a s ea n dl i q u i dp h a s eo f m i x t u r e sw e r ep l o w e d d e n s i t i e sa n dc o m p o s i t i o n sf o rb o t h v a p o rp h n s e a n d l i q u i dp h a s e f r o ms i m u l a t i o ns h o wf a i r a g r e e m e n tw i t ht h e e x p e r i m e n t d a t ao v e raw i d er a n g e o f p r e s s u r e s k e y w o r d s ：g i b b se n s e m b l e ，m o n t e c a r l os i m u l a t i o n ，v a p o r - l i q u i de q u i l i b r i a ， b i n a r ym i x t u r e s ，p u r ec o m p o n e n t s 重庆火学硕士学位论文 绪论 1 绪论 1 1 分子模拟的意义 众所周知，化工热力学数据对于化学工业过程的设计、操作以及优化具有举 足轻重的指导意义。热力学数据一般通过三个途径取得：即实验测定、回归模拟 以及计算机分子模拟。随着计算机技术的飞速发展，计算机分子模拟已经成为存 分子水平上研究流体的一种强有力工具。通过计算机分子模拟，可以较为严格地 从流体的微观相互作用出发，从而预测出流体的宏观性质，特别是在一些极端的 条件：如高温、高压下的宏观性质。这样既可节省大量的实验，又可为工程提供 依据。 采用计算机分子模拟方法，可以得到相当可靠的径向分布函数、热力学性质 以及迁移性质。这为构作与改进各种描述实际现象的理论或模型提供可靠的依 据。计算机分子模拟一方面是模型和理论预测的桥梁，另一方面又是模型和实验 的桥梁，被称为“计算机实验”【1 1 。近2 0 年来，该方法得到了迅速而广泛的发 展。 计算机分子模拟借助分子物理模型和统计力学方法，用计算机进行统计计 算，模拟出宏观物理性质。分子理论应用于研究物理性质包括两步：( 1 ) 体系中 离子或分子间相互作用势的提出；( 2 ) 体系中这些相互作用势的统计力学方程的 求解。第1 步中分子间相互作用势不能用实验方法求解，通常必须用量子力学方 法的从头计算法找出或者从理论与实验数据( 半经验作用势) 的相互联系中得出。 第2 步可用近似法计算或者在计算机上求数值解，求数值解的技术就称为分子模 拟1 2 l 。分子模拟有两种基本方法。一是m o n t ec a r l o 法( 简称m c ) 【3 1 。基本模拟 过程是在一定系综条件下，将系统内粒子进行随机的位移、转动，或粒子在两相 间的转移。根据给定的分子位能函数，进行粒子间内能的加和。采用给定的取样 方法，生成一系列体系的微观粒子随机构型，从而逐渐趋近于平衡时的 b o l t z m a n n 分布。另一种方法是分子动力学方法( m o l e c u l a rd y n a m i c s ，简称m d ) 。 基本的模拟过程是在定的系综及已知分子位能函数条件下，从计算分子间作用 力着手，求解牛顿运动方程，得到体系中各分子微观状态随时问真正的变化，再 将粒子的位置和动量组成的微观状态对时间平均，即可求出体系的压力、能量、 粘度等宏观性质以及组成粒子的空间分布等微观结构。在m d 模拟中，又可根据 是否对体系加一外力场而区分为平衡m d 模拟( e m d ) 和非平衡m d 模拟( n 跏d ) 。 外加力场的n e m d 模拟方法目前尚处于初级阶段。 分子模拟方法不仅可以提供对系统微观结构的直观了解，其结果可用于检验 重庆大学硕士学位论文 绪论 和改进各种理论模型，而且它可以直接提供实验数据，特别是在某些极端条件时 分子模拟方法更显示其优越性。 1 2 分子模拟在化学工程中的应用 许多化工新技术的出现，如膜分离技术、超临界技术、双水相授术、电泳技 术、界面技术、纳米技术等，伴随而生的是一些复杂分子，如高聚物大分子、生 物大分子、两亲分子等。化学工程工作者迫切需要从分予水平来研究系统的微观 结构及其宏观热力学性质及传递性质。从分子水平来研究化工过程及产品的开发 和设计，无疑是2 1 世纪化学工程的一个重要方向n 1 2 1 分子模拟研究流体的扩散性质及超临界状态 分子扩散系数是化工传质过程中最重要的物性参数之一，目前扩散系数数据 的取得主要来自实验测定或经验公式估算。统计力学理论计算还很不完善，仅能 用于稀薄气体，而m d 分子模拟已成为取得传递性质数据的重要手段及建立理 论模型的主要依据。模型流体( 硬球、方阱以及l _ j 流体) 的自扩散系数方程首先 是从搜集大量分子动力学数据拟合得到1 5 j 。后来又在e n s k o g 硬球流体方程基础 上，用m d 模拟数据分别拟合硬球及硬球链流体的校正函数【6 】，从而得到可用于 实际的球形和链状l j 分子( c 原子可达】5 4 个) 的自扩散系数的理论方程。m d 分子模拟方法除了用于建立理论模型外，还可直接测定中性分子、离子和带电胶 体粒子的自扩散系数和互扩散系数。研究表明，在超临界条件下，进行扩散系数 的m d 测定是可靠而易于实现的。超临界流体具有良好的溶解性能和较高的分 子扩散率，因此是一种很好的溶剂。m c 方法可以用于模拟体系在超临界条件下 的热力学特性及吸附平衡i ”。此外，由于超临界流体的高度挥发性，在萃取产物 中溶剂的污染很小，因此可以作为萃取剂。m c 分子模拟已经用于计算n 石蜡 ( c 2 8 h 5 8 和c 3 0 h 6 2 ) 以及高级醇( c 1 6 h 3 3 0 h 、cl 8 h 3 7 0 h 和c 2 0 h a t o h ) 在3 0 8 2 k 时在超临界c 0 2 中的静态性质【8 。c 0 2 简单处理为单点分子，而链分子近似处理 为多点分子，通过引入一个不同点间的分子间参数就可以定量计算在超临界流体 c 0 2 中的溶解性( 固气相平衡) ，n 一石蜡和高级醇的分离性能在超临界c 0 2 中 随着压力的增加而增加，c 0 2 对高级醇的径向分布函数显示c o z 趋于在高分子团 外围聚成束。 1 ，2 ，2 分子模拟用于建立状态方程及反应机理的研究 众所周知，已有众多的立方型状态方程普遍用于化工过程的设计与优化，但 这些方程大多属经验或半经验方程。自2 0 世纪9 0 年代以来，已发展了一类分别 以统计缔合流体理论和硬球链理论为基础的状态方程，可用于链状流体和氢键缔 重庆大学硕士学位论文 绪论 台流体的p v 了计算，现已为国内外普遍关注和应用。另一类是用统计力学平均 球近似( m s a ) 方法求解o m s t e i n z e r n i k e 积分方程而得到的状态方程，可用于 y u k a w a 流体的p 计算p j 。还有一类状态方程是在硬球一阶微扰理论基础上， 直接从l _ j 流体的分子模拟数据回归得到的维里展开多参数状态方程。由此可知， 要获得精度高的具有统计力学基础的状态方程，必须有相当数量在很宽广温度和 密度范围内各种模型流体的分子模拟数据。在表面反应机理和动力学方面， m o n t ec a r l o 方法的研究范围较为广泛。魏庆莉，王贵昌等人对c s t r - b z 复杂化 学反应进行m c 模拟1 1 0 1 。他们在m c 算法中引入流率项，通过对具有流率分又 结构的c s t r - b z 反应体系c n b ( s h o w a l t e r - n o y e s b a r - e i o 模型的随机分析，讨论 了最低反应浓度的意义及流率项连续性的影响。由于c o 在化学工业以及环境 保护方面的重要性，针对c o 的研究非常广泛。研究表明，c o 氧化反应体系具 有复杂的动力学行为，其中包括表面吸附物种的动力学相变、速率振荡和催化剂 表面的结构重排等。曾健青等人【l ”用m o n t ec a r l o 方法研究了表面扩散对c o 氧 化反应一级相变点的影响，将催化剂表面拙象为一个具有理想的二维方格网络， 每一个网络点代表催化剂表面的一个活性中心，模拟中不考虑表面吸附物种的脱 附和吸附态氧原子的扩散，表面c o 的扩散增加了c o 和o z 的接触有利于氧化反 应的进行大大抑制了表面c o 的聚集，从而减缓催化剂的c o 中毒。他们还用 m o n t ec a r l o 方法研究c o 表面催化反应【1 2 1 具有相变和反应速率自振荡的复杂动 力学行为。研究表明反应不仅能产生反应速率上的自振荡现象，而且表明自振荡 与脱附速率无关。郭向云等人把m o n t ec a r l o 方法应用于研究催化反应中的颗粒 形成过程、表面吸脱附过程、表面反应机理、分形表面上的反应动力学u 。 1 2 3 分子模拟研究分子微观结构 分子模拟是研究分子微观结构的强有力手段，如在采用t i p 3 p 及t i p 位能函 数分别对水和甲醇进行m c 分子模拟时，即可区分氢键缔合和非氢键缔合的分子 数，确定分子周围氢键的配位数分布，从而得出：水分子的氢键平均缔合点数为 2 9 3 ，而甲醇分子的氢键平均缔合点数为1 5 8 【i ”。在采用s a f t 状态方程时，即 可参考上述分子模拟结果。与此同时，从分子模拟还可得到体系的静电能、离子 偶极能( 即离子水化能) 、偶极能以及l - j 色散能等数值，这些数值对构筑以微扰 理论和平均球近似理论为基础的电解质分子热力学模型是十分必要的。在液液分 层方面，对l j 流体i ”】、k e e s o m 流体 1 6 】以及高分子格子模型都进行过分子模拟研 究，考察了出现会溶点与分子微观参数间的关系；还发现了在超临界水中苯的亲 水性1 1 7 】。沈良骏用该法研究聚乙烯尾形链构象【1 8 】：他以聚乙烯为例，根据其尾 形链的实际键长、键角，用m o n t ec a r l o 模拟方法生成本链分子，计算尾形链的 均方末端矩及其分量，构象数。此外，m o n t ec a d o 还可以研究聚合物体系的性 重庆大学硕士学位论文 绪论 质f j 9 1 。 1 2 4 分子模拟研究相界面 分子模拟不仅能测定流体主体内分子的近程有序( 即所谓局部组成) ，还能测 定汽液界面、液液界面以及气固、液圃界面的分子分布。对纯流体垂直于汽液界 面的分子密度剖面分布及界面层厚度均可从分子模拟方法得出的双曲正切函数 表示1 2 。虽然也可以采用密度泛函理论与s a f t 状态方程中的微观参数预测纯流 体( 含链状分子和氢键缔合分子) 的表面张力和液液不互溶体系的界面张力。但上 述计算仅能获得链节数沿垂直界面的分布，不能得出链状分子在界面附近的几何 构型。后者必须借助于分子模拟。采用分子模拟可同时从微观和宏观研究固体表 面上的吸附现象。m c 法模拟分子在固液界面的吸附性质【2 ”，可测定吸附层厚度、 分子构型分布、表面覆盖率等一般实验难以获得的信息。汪文川、周健用g i b b s 系综的m c 法模拟甲烷的吸附平衡1 2 2 ：在固定的温度下，对甲烷在1 9 1 n m 的活 性炭孔中的吸附平衡进行模拟，最终得出了该温度下甲烷的吸附等温线。类似的 研究还有c s z 在石墨上的吸附【2 ”。在高压状态下，m c 方法可以对气体的物理吸 附进行模拟【2 “，该研究的目的有两个：得出分子分布函数以及单层覆盖等用实验 方法不易获得的数据；在完全确定的体系下，得出相关的吸附数据。另外，在高 压下，还可以对稀有气体的吸附进行m c 模拟【2 5 1 。 1 2 5 分子模拟研究流体性质 分子模拟方法在流体物性领域的研究时间最长，相关的理论也最为成熟。 其研究内容主要集中在计算流体的化学势、水溶液的活度系数、各种体系的汽一 液相平衡等。吴雄武，时钧1 2 6 】对l j ( l e n n a r d j o n e s ) 流体的化学势；金文正， 汪文j i i t 2 7 】对气体混合物各组分化学势都进行了m c 模拟。郭明学，汪文) l i t 2 8 1 运 用g e m c 方法，对a r 和c 地纯流体汽液共存性质进行了模拟计算，得到各自在 三个温度下汽液两相的饱和密度和汽化潜热。s h y a m a lk n a t h 对支链烷烃的汽液 平衡进行模拟【2 9 】。他针对直链和支链烷烃提出了统一原予力场概念，用于研究 单组分及二元组分烷烃的汽液平衡。j u a nj d ep a b l o 和m a r i a n n eb o n n i n 对碳氢化 合物以及包括c h 4 的二元混合物的相平衡进行模拟【3 0 j 。g r a n ts 和h e f f e l f i n g e r 以及f r a n k v a ns w o l 对圆柱形孔隙中流体的汽液共存情况进行m c 模拟口“。大家 的研究多集中在二元混合物体系的平衡，j ，r e s c i c 和v v l a c h y 对四组分混合物进 行研究【3 2 】。该体系包括两种电解质、一种普通离子型化合物以及中性的物质。 研究表明，由于中性物质的空间作用明显，对化合物的相平衡特性影响显著。从 二维体系变到三维体系后，流体的临界温度增大，相图的形状变宽。但是，在三 维体系里，采用g e m c 方法更容易得到流体的临界点。科学研究需要不断创新， 4 重庆大学硕士学位论文绪论 h l o u c h am 等人就尝试着对四维空间的体系进行模拟【3 3 】。采用g e m c 方法作出 体系的相图，并得出四维流体的临界参数。另外，m c 方法还可以用于研究溶液 的活度系数。金文正等人阐述了无限稀释活度系数的模拟方法1 3 4 , 3 5 】，建立了分子 蜕变偶合参数积分m o n t ec a r l on p t 系综法，模拟无限稀释水溶液的化学势差和 活度系数。 1 3 流体相平衡分子模拟的方法及研究现状 从上世纪五十年代以来，得益于计算机技术的发展，计算机分子模拟也经历 了几个阶段，其研究的体系也越来越多。从最初的单原子体系到双原子体系，再 到碳氢化合物体系，商分子聚合物体系，多元混合物体系，所研究的体系也变得 更加复杂。流体的性质分为两大类：一类是热力学性质即平衡性质；另一类是迁 移性质。流体相平衡性质揭示相平衡时温度、压力、体积、各相组成及其他热力 学函数之间相互关系的规律。它是蒸馏、吸收、结晶、萃取等单元操作的理论基 础，因而对许多工业过程的设计、优化和操作是必不可少的。分子模拟方法在流 体相平衡研究中所起的作用也变得越来越重要，它从分子间的相互作用入手，利 用严格的统计力学规律和先进的计算机技术，对流体的热力学性质和相平衡进行 模拟计算。流体相平衡的分子模拟方法：直接模拟法和间接模拟法。 1 3 1 、直接法 直接法中应用最广泛的是吉布斯蒙特卡罗法( g i b b se n s e m b l em o n t ec a r l o s i m u l m i o nm e t h o d ) 陬”j ，简称g e m c ，由p a n a g i o t o p o u l o s 提出，此方法不考虑 相边界，无需直接计算化学势，模拟计算比间接法简便，因而它广泛用于研究简 单流体体系的相共存性质，如球形l j 流体、y u k a w a 离子体系和一些较复杂体 系。在1 9 9 2 年，d ep a b l o 等人1 3 8 】运用g e m c 模拟了甲烷与正戊烷的相平衡，模 拟结果与实验值只有百分之几的误差。在1 9 9 8 年，a n p i n gl i u 和t h o m a sl b e c k g g l 成功地应用g e m c 模拟了二氧化碳与甲烷以及二氧化碳与乙烷的二元混 合物的汽液相平衡。在2 0 0 1 年，有学者【4 0 1 将吉布斯系综蒙特卡洛法与反应系综 法结合起来，成功模拟了复杂的化学反应体系的汽液相平衡，此方法称为反应吉 布斯系综法( r g e m c ) 。 1 3 2 间接法 流体相平衡的各种间接模拟方法都需要计算一系列状态点的化学势。通过在 相空间中寻找具有相同化学势、相同压力和相同温度的两相来确定相平衡状态。 当模拟复杂分子、固体相或液态晶体相时，间接法具有显著的优越性。在多年的 研究过程中，发展了几种方法：检验粒子法f 4 1 1 ，巨正则系综法( g c m c ) 1 4 孙、 改良抽样法和热力学积分方程法等。 重庆大学硕士学位论文绪论 将此方法应用于聚合物的模拟是一个重要的进展。k u m a r 等人报道了一种基 于检验粒子方程的计算长链分子化学势的方法【4 3 出1 。该方法的特点是在链分子每 增加一个单体，化学势增加一定的值。另一种是基于称之为构型偏倚m c 法的计 算长链分子化学势的方法【4 5 ，4 6 1 。该方法采用在流体中插入一个短链，再逐步加入 其它片段直到形成完整的链分子。 1 4 本文研究的目的和内容 本文依据g e m c 的基本原理，运用f o r t r a n 9 0 结构化程序语言，计算了 纯流体和混合物流体的汽液相平衡数据。通过改变程序输入文件里的各项参数， 控制程序运行中采用的势能模型和分子结构。重点模拟了多种极性分子和非线性 分子纯组分流体及其混合物流体的汽液相平衡，为模拟更复杂的体系积累经验。 本文中引入电荷的相互作用，二氧化碳采用有库仑力存在的三中心线性分子模 型；丙烷、正丁烷和正戊烷采用非线性分子模型。通过模拟不同条件下的二元混 合物体系的相行为，获得相关的汽液相平衡数据。 重庆大学硕士学位论文吉布斯系综, q o n t ec a r l o 模拟的基本理论 2 吉布斯系综m o n t eo a r lo 模拟的基本理论 2 1 抽样方法和周期边界条件 抽样方法包括简单抽样方法，重要性抽样方法，m e t r o p o l i s h 7 1 抽样方法及其 它改良的抽样方法，如伞型抽样方法，构型偏倚抽样方法1 4 8 1 等。 2 1 1 重要性抽样方法和伞型抽样方法简介 重要性抽样是从简单抽样方法发展而来的，其基本观点是，构型的分布函数 是波尔兹曼分布时，在波尔兹曼因子大的地方多抽样，而在其它地方则少抽样。 这样进行非均匀构型空间抽样，以使得权重函数近似正比于波尔兹曼因子。 伞型抽样方法是t o r r i e 和v a l l e a u 【4 9 l 建议的，为了到得到模拟体系0 和体系1 自由能差的精确估算，分别对体系0 和体系1 的空间构型进行抽样，通过用一非 负的权重函数来代替体系的波尔兹曼因子。访问构型空间的一个点的概率正比于 权重函数。伞型抽样效率高，但是要确定的权重函数是未知的，它是通过模拟的 两个体系的波尔兹曼权重方面的信息来构造的。构造一个好的抽样分布需要技巧 和耐心。最近也有一些发展，如v a l l e a u t s 。 5 1 1 建议用一对称的方法计算所需的权 重函数。 2 1 2 m e t r o p o lis 抽样方法 m e t r o p o l i s 等人通过在位形空间构作一个马尔科夫链，使样本出现的概率随 着马尔科夫链的增长逐步趋近于b o l t z m a n n 分布，且随着马尔科夫链的无限增 长，体系中各物理量的平均值接近于平衡分布的系综平均值。 下面就简要介绍一下m e t r o p o l i s 抽样的基本步骤。假定一个体积为v 的矩形 箱中有n 个粒子，势能模型为l k ( ，n ，n 一) ，则新的位形按以下步骤产生： ( 1 ) 从n 个分子中随机地抽取一个分子( 也可轮流) ； ( 2 ) 让这个分子试着移动一步，即将分子质心坐标和方向改变为： 玛=x + 占1 吾1妯=y + 占1 舌2 磊=z + d l 善3 口口= 0 + 6 254 函日= 咖+ 62 5_ = f r 。= i r + 6256 式中：鼠y 、z 空间位置 口、咖、 f r e u l a r 角，代表分子在空间的取向 6 l 、62 一最大允许位移和旋转度 51 ；6 一 o ，1 】上的分布的随机数 ( 3 ) 在新的位形下，计算新的位能u 。； ( 4 ) 如果l l 丑 e x p ( - j bl 1 0 ，于是接受新的位形作为马 重庆大学硕士学位论文 吉布斯系综m o n g ec a r l o 模拟的基本理论 尔科夫链的一个状态点； ( 5 ) 如果户。，则以正比于e x p 一卢( 己k l 珀几率接受新位形。为此将 e x p 【一卢( 。一l 铀】与一个【o ，l 】上的随机数进行比较，如果前者大于后者， 刚接受新的位形；反之，则将原位形荐次作为马尔科夫链的新的状态点。 m e t r o p o l i s 抽样方法就是从某一个初始构型( 比如将n 个粒子按面心立方晶 格放置) 出发，不断重复上述步骤产生一系列构型，最后作必要的统计。为了获 得满意的模拟结果，一般需要数百万乃至数千万个构型，且参数占1 和6 2 还要选 得合理。参数的选择和接受率( “试移”被接受的次数进行加和，再除以总移动 次数，这一比值称为接受率) 相关联。 2 1 3 周期边界条件 原子和分子体系的m c 及m d 模拟的目的是提供一个宏观样本的物性信息。 但当今的计算机中通常还只能处理几百至几百万这样数目的自由度，绝大多数模 拟只检测几百至几千个粒子的体系的结构及热力学性质，这个数目与热力学极限 相去甚远。粒子处于界面的概率很大，如1 0 0 0 个粒子的简单立方晶体构型的体 系，约4 9 的粒子位于界面。为了消除这种表面效应，需要引入周期边界条件。 如图2 1 所示，“中心盒子”为n 个模拟分子组成的标本盒子，其周围盒子是虚 拟的映像，并且像盒子是无限延伸的。 0o o 。o o oo o oo o o o 。 。 。oo o o o 口o o 。o 。 。 o 。o o o o oo o 。 or 屯 怠，。 d 。o o o 0 o。一。 。o o l 口o 。 二。， 三 uo o oo o r 。c f z | i oo o o0 o o o o o o o 。 o oo o oo o o 。 。 o o o o o o o o 0 o 。o o 。o o 。o o 图2 1 周期性边界条件的平面示意图 f i g 2 1l c h n o g r a p h y0 fp e r l e d i eb o u n d a r i e s 重庆大学硕士学位论文吉布斯系综m o n t ec a r o 模托f 的基本理论 2 2 相互作用的截断 在周期边界条件下，由于不同盒子之间允许存在相互作用，因此，每一个分 子与其它分子的相互作用延续至无穷远，势能的计算也变得复杂起来。但是，对 于文中所应用的l j 势能模型而畜，势能随分子问距离的增大衰减得很快，因此， 在计算势能时可以近似地截断到适当的距离，而远程分子之间的相互作用就可以 忽略不计。下面就简要介绍一下常用的两种近似方法，即“最近映像”和“势能 球形截断”。 如图2 1 所示：每一个分子必有而且只有n - 1 个近邻分子不重复地落入以该 分子为中心的盒子( 不是原来的盒子，但尺寸相同) 中。这n - 1 个分子可能是中 。i i , 箱中的，也可能是像盒子中的，它们是中心分子的紧邻分子，称为最近映像( 分 子) 。在最近映像近似中，中心分子与其余分子的相互作用能只能考虑到紧邻分 子，也就是说将分子间力的作用范围截断到模拟盒表面。 在采用最近映像规则后，两两相互作用势仍然包括0 5 n ( n - i ) 项。对于2 5 6 个粒子的体系，仍然是很大的，为了进一步减少计算量，采用球形势能截断。将 截断距离r c 确定为一个常数( 在程序中取模拟盒边长的一半) ，作用范围则是以 中心分子为圆心、r c 为半径的圆球，称为截断球。在该条件下，与中心分子发生 相互作用的对象仅限于球内分子。此时，与中- i i , 分子发生作用的分子数就明显小 于n - 1 个。不难看出，使用截断方法会使计算的热力学性质产生一定的偏差。为 此，可以引入长程校正方法来解决这一问题。 2 3 长程校正 7 对模拟结果进行长程校正可以提高模拟的准确度。 e 删= e c + e 腓= e c + 2 刀 p 【r 2 v ( r ) d r ( 2 1 ) ( p v ) 。= ( p v ) 。+ ( p v ) 。= ( p v ) 。+ 玎n p r r 2 c o ( r ) d r ( 2 2 ) 肛血= c + a t t c = p c + 4 印【r 2 v ( r ) d r ( 2 - 3 ) 式中点钿，俨功舢，_ “。，分别表示模拟所要得到的体系的总体值，e - ( p 垓，_ u 。 是体系在丘处截断计算得到的值，e z x c ，俨。阳，。即为长程校正项a 对于 恒容吉布斯系综和恒压吉布斯系综，由于模拟盒的体积和里面的粒子数可以改 变，它们对于瞬间的能量值和压力在是有很大的影响，从而使得模拟体系的密度 波动很大，所以进行长程校正是必不可少。 蘑庆大学硕士学位论文吉布斯系综m o n t ec a r l 0 模拟的基本理论 2 ，4 吉布斯系综原理 2 4 1 系综原理 统计力学旨在根据物质的微观状态的动力学行为来解释宏观状态的性质，其 基本原理就是系综原理【5 2 , 5 3 1 。热力学状态是由一系列参数确定的，热力学性质可 由状态方程和热力学基本方程得出。虽然象扩散系数d 这样的量与体系的微观 结构和动力学性质有关，但它也是状态函数，其值完全由少数变量( 如、p 、d 决定，而不是由瞬间的力学状态比如原予位盟、动量等决定。在多维相空间中某 一点a ，其某一瞬态性质m 的值记为撇) ，有理由设宏观观测值m o k 是各瞬时 朋似，的平均值： 10 a 2m ) t i m e2 ( m ( 以f ) ) ) 2 憋j m ( a ( o ) d t ( 2 4 ) 0 牛顿运动方程仅是常微分方程，根据所需精度用计算机解此类方程是完全可 行的。只需要在一个较长的有限时间内取平均值就行，这就是分子模拟所采用的 方法。事实上，解运动方程是一步一步进行的，步长占，= 形，f 为有限步数，此 时，可将方程写为： 1i m = ( 聊。= m ( 4 f ) )( 2 5 ) l 由分子动力学方法进行的时间平均不是传统统计力学意义上的热力学方法。 大量分予体系的 d 矗晔) 汁分复杂，要进行时间平均有一定的困难。为此，g i b b s 在1 9 0 1 年提出了一种巧妙的方法：他将大量宏观状态相同的系统集合在一起， 而每一个系统各处在它所辗转经历的某一个微观状态中，这样的系统称为标本系 统，由此而构成的标本系统的集合就叫做统计系综( 简称为系综) 。因为标本系 统的数目n 很大( n 的数量级通常为1 0 2 4 ) ，足以代表所有的微观状态，所以对 于给定的热力学系统，某一性质的时间平均值就等于它的系综平均值。通常将系 综看成是相空间所有点a 的集合，各点按概率密度p 分布。它由指定的宏观参 数n 、u 丁确定，因此采用符号佛w 来标记，一般用j d 缸表示。每一点代表任一 特定瞬间的系综，每一体系随时间变化，遵循力学运动方程而与其它系综无关。 如果幽。似康示某平衡系综，即，( m 锣卢0 ，则体系的变化显得相当特 殊：某一体系离开状态搿力到另一状态川计p 。就有另一体系从a ，- 口状态进入 来代替它。如果通过相空间的所有点中只有一个轨迹p 。( 砷为非零，那么每一体 系最终在所有点都可预见，这称为“各态遍历”。然而现在关心的是平均性质， 可以对系综的每个成员的平均来代替时间平均： 重庆大学硕士学僚论文吉布斯系综m o n t ec a r l o 模拟的基本理论 m = (