中国工程热物理学会第十一届学术会议 分子动力工程热力学及匏源利用编号：0 5 1 1 2 1 学模拟中的变截断半径算法蔡治勇曾丹苓刘娟芳( 重庆夫学动力丁程学院，重庆。4 0 0 0 4 4 )T e l ：0 2 3 6 5 11 2 1 8 4 ；E - m a i l ：d l z e n g c q u e d u c n搐要：在分子动力学模拟中，人们往往根据所研究的对象选取一个固定的截断半径，而不同的作者所推荐的有效截断半径值十分分散。这一方面影响了模拟结果的精度，同时也造成了截断半径选取的困难。本文提出了一种新的变截断半径算法，可以根据对象所处的状态灵活改变截断半径。文中以两相平筏系的非均匀各向异性区的模拟为饲，采用所提出的方法得到了更接近实验的结果，同时也节约了机时。关键词：变截断半径：分子动力学模拟；算法流体的分子动力学模拟近年来得到很火的发展，已证明是一个研究流体物性的重要工具。在分子动力学模拟中，截断半径的选取对流体的模拟结果的影响非常重要，从理论上讲，截断、r 径越火计葬精度越高：但实际上：截断半径的选取受模拟盒尺寸的限制，同时，截断、”径越大，参与作J 【 ；| 的分子数越多，计算速度越慢，计算处理误差越大。传统的模拟方法在模拟过程辛截断半径通常选取为一个固定的数值。由于气一液相平衡体系中，均匀相、非均匀相的分布特性不同，传统的方法中统一的截断半径不可能同时满足均匀相中饱和气体、饱和液体密度以及非均匀相中表面张力的计算精度要求；本文首先采用固定截断半径对L J 流体气液相平衡进行了分子动力学模拟，得到了不同的截断半径下，流体均匀相中饱和气体、饱和液体密度以及非均匀相中表面张力的计算结果。证明了对均匀相、非均匀相的热物理特性计算选取相同的截断半径是不合适的，需取不同的截断半径，才能满足要求。此时对非均匀相特性计算所要求的截断半径癍大= 丁i 二满足均匀相特性计算要求的截断半径。接着，本文采用变截 断半径算法一当分子处于均匀相中截断i F 径l 缸小，分子处于非均匀相中取大的方法， 对L J 流体气液相平衡进行了分子动力学模拟，得到新的流体均匀相中饱和气体、饱和液体密度以及非均匀相中表面张力的计算结果；结果表明：在均匀相中取固定截断半径R = 3 0 4 0 ，在1 I 均匀相中，根据对象所处的状态灵活取截断半径时，均匀相、非均匀相的热物理特性能够同时得到满足。与传统的定截断半径方法相比，变截断半径法计算精度更高计算时间更短。1 分子动力学模拟方法对L J 流体，分子i 和j 间的相互作用势函数表示为 M 州俐1 4 1式中，为能最参数，0 为尺度参数，都是与具体流体性质相关的常数，r I 。为分子i和j 之间的距离。运用L e n n a r d J o n e s 势函数进行分子动力学模拟中采用赢角坐标系O x y z ，在初始时刻流体处于饱和液相，分子以面心立方晶格( F C C ) 均匀分布在边氏为L 的正方体模拟盒中，其质心在原点( 0 ，0 ，0 ) 处。模拟过程中，系统采用正则系综即分子数N 、容积V 和温度T 保持不变，并不断调整体系质心使之处于原点。模拟盒在X 、Y 方向的长度L x ：L y = L ，为了保证足够的气相空间，取在z 方向的长度L z = 3 L ( 见图1 ) 。本模拟在三个方向上都采州周期性边界条件。在计算机分子动力学模拟过程中，为确保模拟体系各物理攀之间的数鼙级较为接近，减小模拟过程中出现的误著将所用劁的各物理量进行无鹫纲化。采用无餐纲物理量进行分子动力学模拟，分子的备物理量利 用L - J 1 2 - 6 势函数中的参数s 和盯，以及体系中分子质量m 对模拟流体备物理簧进行无量纲化处理如下：l|v a p o rl i q u i dv a p o r，|l图1 横拟盘尺寸圉长度、L = L o时间t ：t ( m 孑) ”温度r - l 叮店表面张力Y - Y 叮2 8能量t l * = t l E ；密度p * = p a 3 m压力昨P 仃y 2 计算公式 在进行分子动力学模拟时，分子的运动规律由联立求解各分子的运动方程弘川z 4 r 0 ) 饥争川玉，v得出，对L J 流体，取截断半径为R 时，分子间相互作用力F k J 为： F “) = 一掣?( 3 )本文模拟的晟终结果是计算气液相平衡中均匀相的饱和气体、饱和液体的密度以及非均匀相的表面张力，模拟中将模拟盒沿z 方向平行于x Y 平面平均分为N s 个切片进行数理统计，均匀相的分子数密度：刖= 等界面层中非均匀相的表面张力：1 4 2何 以1 训，弋 陪坩丝勺咖，主 釜，( 七)k = l 州= 壶( 占，。掣嘲。( 5 ，式 为正则系综统计平均值。在求和之中，如果i 和j 两个分子都在第k 切片内，则将其全部计入；如只有其中一个分子位于第k 切片内，则计入其贡献的一、h在本文的模拟中，均将模拟盒沿z 方向等分为N ：1 0 0 个切片，模拟分子数为2 0 4 8 0由于体系达到平衡需要很长的模拟时间，为了得到较为真实的平衡体系，本文取时间步长t = O ，0 0 1 ，达平衡态的模拟艮度为1 0 万分子动力学步，然后以1 2 万分子动力学步进行统计。3 模拟结果与分析 3 1 传统的模拟方法为了得到准确的截断半径对均匀相饱和气体密度、饱和液体密度以及非均匀相表 面张力的影响，本文首先采用传统模拟方法，以L J 流体 r 分子( 参数铲3 4 0 9 x 1 01 0 m ，= 1 6 5 x 1 0 1 1 J ，分子质量m = 6 6 4 x l f f ”k g ) 为模拟对象，初始参数设置为：温度 r = o 7 5 1 3 ，密度衙I ! = 穰7 9 0 6 , 分子数为2 0 4 8 ，取截断半径R - - 2 5 、3 0 、3 5 、4 0 、4 5 、5 - 鼠6 5 等七神情况进行了分子动力学模拟。饱和气体密度和饱和液体的密度以及表面张力模拟结果见表1 ，饱和液体密度和饱和气体的密度及表面张力的模拟结果随截断半径的变化见图2 。豳2 汽相密度、液相密度、表面张 力随截断半径的变化R表1 饱和气体密度和饱和液体的密度Rp LP vY 2 50 7 5 4 70 0 1 2 7 30 5 4 8 63 00 7 8 2 20 0 0 ”8O 6 2 7 33 50 7 9 4 20 0 0 5 1 20 7 1 8 94 00 8 0 3 2O 0 0 4 6 20 1 7 6 9 545O 8 0 6 80 0 0 4 1 30 8 1 2 55 5O 8 1 0 l0 0 0 3 7 50 ，8 4 0 36 50 8 1 3 00 0 0 3 6 70 8 5 2 9从表l 和图2 可以看出：1 在均匀体系中，饱和液体的密度所随截断半径R 的增大而增大；饱和汽体的1 4 , 3密度p ，则随R 的增大而减小。当截断半径较小时，两者的模拟结果均与实验值相差较大，在R = 3 5 。4 O o 时，模拟结栗与实验值( 实验值：研= O 7 9 0 6 ，P 。= O 0 0 4 9 6 7 ) 十分接近，饱和液体的密度误差在】以内，饱和汽体的密度误差在4 以内。2 在非均匀体系中，表面张力Y 的数值随截断半径R 的增加而增加，并逐渐向其实验值逼近；肖截断半径R 取值较小时，模拟结果的误差较大；在R = 3 5 4 O a 处，表面张力的计算误差仍在8 1 4 左右；在S = 5 5 0 处，表面张力的模拟结果( y = O 8 4 0 3 )与实验值( 实验值：y _ o 8 3 5 3 ) 十分接近，但此处饱和汽体的密度的计算误差又明显加大。达到约2 5 。3 不存在一个能同时满足均匀相饱和气体密度、饱和液体密度以及非均匀相表面张力特性计算的最佳的截断F 径；对于均匀相、非均匀相各自的热物理特性计算，需取不同的截断半径，满足非均匀相特性计算要求的截断半径应大于满足均匀相特性计算要求的截断半径。3 2 改进的模拟方法分析上面的模拟结果，我们认为：造成均匀相、非均匀截断半径不相同的主要原 因是因为气一液界面层内分子的不均匀分布造成。一般说来，运甩势函数时，我们能对分子i 所受的来自予它截断半径内的所有分子对它的作用力能准确计算。当截断半径较小时：在均匀相中，来自于截断半径以外的所有分子对这个分子的附加作用力被平衡，故对它的运动轨迹没有影响；丽在非均匀相中，自予截断半径以外的附加作用力则不能被平衡，当求解运动方程时，这个附加作用力必须考虑进去，所以，计算非均匀相的特性时，需要较大的截断半径。 A L o t f i 等早就注意到传统的模拟方法的缺陷，第一次将长程修正( L R C ：l o n g - r a n g ec o r r e c t i o n s ) 的概念引入到分子动力学模拟中；硅碡e c k e 等则在传统的模拟结果的基础上附加一个表面张力部校正值Y “，的方法来计算表面张力； I l d r i j ，T 等从势函数截断半径和作用力函数截断半径不同步出发，提出s T ( s p h e r i c a l l y t r u n c a t e d ) 势函数模型和S T S ( s o N e r i c a l l yt r u n c a t e da n ds h i f t e d ) 势函数模型的概念，指出在进行传统的模拟时，必须在截断半径处加上一个附加作用力。本文提出一种新的计算方法变截断半径的计算方法，该方法充分利用已得到 的计算计算信息，定义饱和气相密度的1 1 5 和饱和液相密度的8 5 构成的区域为气一液界面区，根据分子所处的位置改变截断半径，当分子处于均匀相时，截断半径取表2饱和气体密度和饱和液体的密度以及表面张力模拟结果p tP rVR实验值实验值实验值0 7 9 0 60 0 0 4 9 7o 8 3 5 32 507 8 8 300 0 5 5 30 7 9 8 13 007 9 4 30 0 0 4 9 90 8 2 1 53 507 9 8 20 0 0 4 8 60 8 3 2 54 0O8 0 1 20 0 0 4 8 10 8 4 0 71 4 4一个固定的常数( 通常取值较小) ，当分子处于非均匀相中，截断半径取分子到液相区或汽相区中距离较大者的值，由 于气一液界面区的边界是变化的，同时分子所处的位置也是变化的，因而实现了模拟过程中的变截断半径算法。对前述的模拟体系，本文采用变截断半径算法，在均匀相中取截断半径分别取R = 2 5 、3 0 、3 5 、4 0 四种情况进行了模拟，饱和气体密度和饱和液体密度以及表面张力的模拟结果如表2 。分析表l 、表2 ，可以看出：1 在均匀体系中+ 饱和液体的密度舟随截断? 卜径R 的增大而增大；饱和汽体的 密度A ! J ! I J 随R 的增_ _ = 而减小。其值均比相同截断半径下，采用传统模拟方法更接近实验值。2 在北均匀体系中，表面张力Y 的数值较传统模拟方法有较人的提高，更加接近实验值，说明在汽液界面层内加人截断半径一考虑不平衡附加力的作用对表面张力影响十分明显。3 o 采H j 人截断、r ；| = 的传统模拟方法相比，改进后的模拟方法由于在均匀相计算的分子数较少所以仃省计算时间。4 结论本文提出了汽一液平衡共存的多相系统在进行分子动力学模拟中的变截断、F 径算法，并以A ，为例，将模拟结果与传统的模拟方法和实验值进行了比较，取得了以F 结论：1 采川传统的模拟方法，要满足均匀相、非均匀相的热物理特性的模拟要求，需选取不同的截断、P 径，满足非均匀相特性计算要求的截断半径