专利名称：分子动力学模拟中壁面边界的模拟方法
技术领域：
本发明涉及一种分子动力学模拟中的方法，尤其是模拟粒子在壁面边界条件下的 运动轨迹的方法。
背景技术：
现在越来越多的科学家用分子动力学方法模拟平行平板中流体和纳米通道内流 体流动的动力学特性。随着研究的深入，科学家们开始对非平直纳米通道进行研究和采用 分子动力学方法对此通道进行模拟，如=Xi-Jun Fan等人采用非平衡分子动力学方法模拟 了简单流体在周期性喷管中的流动特性；为研究粗糙度对微通道内体流动及其边界滑移性 质的影响，曹柄阳等人研究了氩气在锯齿状钼通道内的流动J. Castillo-Tejas等人分别 对线性聚合链和牛顿流体在4 1 4收缩-扩展通道中的流动过程进行了模拟。在使用 分子动力学方法模拟两平板间和纳米通道内流体的流动特性时，以前的文献中对虚拟热壁 面(Virtualthermal wall)情况下壁面边界条件的模拟有两种方法，一种是Allen，M. P. and D. J. Tildesley，在文献Computer Simulation of Liquids.中提出的全反射法，全反射法 将壁面看成光滑镜面，粒子在壁面上发生镜面反射，采用全反射法模拟粒子在壁面边界条 件下的运动轨迹参见图1(a)，当粒子从A点运动到B’时，由于虚拟壁面的存在，粒子不会 在B’的位置出现，粒子会在壁面发生全反射后运动到B点，B与B’点关于壁面对称，与壁 面平行坐标方向的速度矢量保持不变，与壁面垂直坐标方向的速度矢量大小保持不变，方 向反向；虽然全反射法可以真实地反应出粒子在发生发射后的位置，但是全反射法不能模 拟出粒子在定壁温和壁面不光滑的条件下的运动轨迹，并且如果在非平直纳米通道中，若 采用全反射法来模拟粒子的壁面边界条件下的轨迹，则在计算速度矢量方面也存在一定的 X隹度；另一禾中是 Raraport，D. C.在文献 The art of molecular dynamics simulation.中 提出了随机反射法，随机反射法与物理上的漫发射对应，采用随机反射法模拟粒子在壁面 边界条件下的运动轨迹参见图1(b)，粒子反射到B’向壁面做垂线的垂点B，速度大小由壁 面温度确定，方向随机。虽然随机反射法可以模拟出壁面不光滑以及定壁温的物理条件，但 是随机反射法由于粒子反射后的位置在壁面垂直方向的单一性，因此在模拟粒子运动轨迹 的过程中存在一定的不真实性，不能真实反映出粒子反射后的位置。因此，如何真实模拟粒子在壁面边界条件下的运动状态，成为了本领域内的研究 人员急需解决的技术问题。
发明内容
针对现有技术存在的上述不足，本发明提供一种在分子动力学模拟中，当粒子处 于壁面边界条件下时，可综合全反射法和随机反射法优点，真实体现粒子在壁面边界条件 下的运动状态的壁面边界模拟方法。为了实现上述发明目的，本发明采取以下技术方案分子动力学模拟中壁面边界 的模拟方法，其特征在于，所述模拟方法包含以下步骤
A)定义分子动力学模拟中粒子的数据结构，粒子在三维模型中由六个信息来表示 粒子在分子动力学模拟通道中的运动状态，所述六个信息分别是X方向坐标、y方向坐标、Z 方向坐标、χ方向速度、y方向速度和ζ方向速度；其中χ方向坐标为壁面横坐标，y方向坐 标为壁面纵坐标，ζ方向为与XY平面垂直的方向；B)确定粒子在通道中的当前位置A(X(1，y0, z0, Utl，V0, W(1)，其中X(1表示粒子在χ方 向上的当前位置坐标，ytl表示粒子在y方向上的当前位置坐标，Ztl表示粒子在ζ方向上的 当前位置坐标，Utl表示粒子在χ方向上的当前位置速度，Vtl表示粒子在y方向上的当前位 置速度，W0表示粒子在ζ方向上的当前位置速度；C)根据粒子当前位置及粒子受到其它粒子的作用力信息，确定粒子的运动方向， 并根据粒子当前位置及受力信息判断粒子下一步是否会运动到壁面边界以外位置；如果粒 子没有运动到壁面边界以外位置，则粒子会沿着确定的运动方向运动，到达下一步指定位 置B’(x’，y’，Z’，U’，V’，w’)，其中χ’为粒子在χ方向上的指定位置坐标，y’表示粒子在 y方向上的指定位置坐标，ζ’表示粒子在ζ方向上的指定位置坐标，u’表示粒子在χ方向 上的指定位置速度，V’表示粒子在y方向上的指定位置速度，W'表示粒子在Z方向上的指 定位置速度；D)若粒子运动到壁面边界以外位置，粒子将不能到达下一步指定的位置B’，粒子 根据壁面的温度，经过壁面反射到达下一步反射位置B (X，y, z, U, ν, w)，当χ方向存在壁 面，则其中 y = y’，ζ = ζ’，χ = 2* 壁面 χ 方向的坐标-χ，，u=a* ^JkeTfm ,v=a* ^kBT/m，
Vf= a* Hm，粒子的速度方向随机；当y方向存在壁面，则其中χ = χ’，ζ = y = 2* 壁面y方向的坐标_y’，VV> ^Tjm，w= a* ,]kBT/m，粒子的速度方向随机； 当ζ方向存在壁面，则其中χ = x’，y = y’，z = 2*壁面ζ方向的坐标-ζ，，u=a* ^kBT/m， v=a*」kBT/m ,w= a* ^kBT/m，粒子的速度方向随机。上述公式中，a为-1 1之间的随机数(该随机数的取得采用文献The art of moleculardynamics simulation(D· C· Rapaport，Second Edition，Cambridge University Press, 2004)中的程序进行计算，在此不作详细描述)，kB为波尔兹曼常数，T为当前温度，m 为质量，χ表示粒子在χ方向上的反射位置坐标，y表示粒子在y方向上的反射位置坐标，ζ 表示粒子在ζ方向上的反射位置坐标，u为粒子在χ方向上的反射位置速度，ν为粒子y方 向上的反射位置速度，w为粒子在ζ方向上的反射位置速度。本发明分子动力学模拟中壁面边界的模拟方法，其主要原理是在镜面反射基础上 考虑了壁面粗糙度对粒子速度的影响，因此，本发明方法又可以称为半反射法。通常情况下，统计系统的宏观量，可以通过粒子的运动状态，即在本发明方法通过 粒子分别在X，y，Z方向上的坐标和粒子分别在X，y，Z方向的速度来反映粒子的运动状态， 当粒子运动状态发生变化时，χ, ι, ζ方向上的坐标和粒子分别在χ，y，ζ方向的速度都会发 生变化。当采用半反射方法对粒子处于壁面边界条件下的运动轨迹进行模拟时，其既考虑 粒子的随机性，即粒子在遇到壁面进行反射后，其反射位置速度(u，ν, w)会根据壁面的温 度有关，而其反射位置x，y，ζ方向坐标(x，y，z)则与粒子未碰到壁面的预定的下一步指定 位置的χ，ι, ζ方向坐标(χ’，y’，ζ’ )有关。因此本方法在模拟粒子下一步位置时，不仅考 虑到壁面的温度，也考虑到了粒子经过壁面发生反射后，其位置必定与粒子的未经过反射的位置有关系。本发明产生的有益效果(1)本发明方法不仅可以在分子动力学模拟中体现恒壁温条件和壁面不光滑的物 理条件，而且也可以在分子动力学模拟中真实反映粒子反射后的位置。其解决了随机反射 法中不能真实反应粒子反射后位置的问题。(2)本发明方法还可以在不规则纳米通道中，简化速度矢量计算方法，从而便于计 算。解决了全反射边界法不能模拟恒壁温条件以及壁面不光滑条件下粒子的真实运动状态 的问题和在不规则纳米通道中，全反射法不能简化速度矢量计算方法的问题。(3)本发明方法可以用于模拟粒子在各种不同通道中运动到虚拟壁面边界外之后 的运动情况，从而为研究粒子在不同通道中运动状态，提供了技术支撑。
图1是分子动力学模拟中的粒子运动轨迹示意图；其中(a)为采用全反射法模拟 粒子在壁面边界条件下的运动轨迹示意图；其中(b)为采用随机反射法模拟粒子在壁面边 界条件下的运动轨迹示意图；图2是分子动力学模拟中的粒子在两种边界条件下反射示意图其中(a)半反射 粒子轨迹运动示意图；(b)随机反射粒子轨迹运动示意图；图3是分子动力学模拟中的变截面通道示意图；其中(a)为三维变截面通道示意 图；(b)为三维变截面通道正视示意图；图4是分子动力学模拟中粒子数密度分布随到壁面距离的比较图；图5是本发明方法模拟变截面纳米通道中势能随时间变化曲线图。
具体实施例方式下面结合附图和具体实施方式
对本发明作进一步详细地说明。本文中提出一种综合了全反射法和随机反射法优点的新的壁面边界条件处理方 法——半反射法。当粒子A点运动到B’时，如图1(a)，由于壁面的存在，粒子经过壁面反射 后到达B点，坐标位置与全反射法粒子反射后相同，但是速度大小由壁面温度确定，方向随 机，即速度矢量的确定方法与随机反射法相同。在模拟非平直微纳米通道时，比如非45° 倾斜通道，如果壁面采用全反射边界条件处理法，粒子反射后，速度矢量的确定有一定的难 度。随机反射法相对于全反射法隐含了两个条件(1)壁面恒温；(2)壁面不光滑，与真实 情况更加接近。但是随机反射法在模拟纳米通道时有一定的不真实性。如图2所示，χ方 向为流动方向，当粒子运动到区域A时，如果采用随机反射法，见图2 (b)，无论粒子B处在A 中的任意位置经过坐标调节后都出现在B’，不能真实反应出粒子运动的坐标位置。但如果 采用半反射边界条件模拟方法，如图2(a)，粒子B出现在区域A中的任意位置，采用半反射 法经过坐标调节后可以出现在B对应任意B’点，从而与粒子的真实运动位置更为接近。实施例1本实施例通过在分子动力学模拟过程中，分别采用随机反射法、全反射法和半反 射三种壁面边界条件模拟方法来模拟平行平板中的氩流体流动过程，从而进一步说明本发 明方法与现有的随机反射法和全反射法的不同之处。具体实施如下
采用氩原子模拟NVT系中，用速度修正法控制系统温度，原子间相互作用采 用12-6LJ势能模型，运动方程采用Verlet数值积分法，时间步长为lfs，模拟时间为 lOOOOOfs，模拟尺寸为6nmX6nmXl. 5nm，粒子数为648，密度为lg/cm3.在本实施例中，将 坐标原点定义为0点，χ和y方向采用周期性边界条件，当ζ方向(Z方向是指三维通道中 与xy平面垂直的方向，粒子只可能超出ζ方向的壁面边界，因为χ和y方向采用周期性边 界条件，所以在这个实施例中都是无限延伸的，没有壁面)为0时，xy平面为下壁面，当ζ 为1.5nm时，xy平面为上壁面，下面分别采用三种不同的壁面边界条件模拟方法来获取粒 子的运动状态随机反射法χ= χ，；y = y，；z = 0 ；M= a*^jk13T/m ；V= a*^kBT/m ；w= a*^jkeT/m全反射法x= x，;y = y，;z = 2* 壁面 ζ 方向的坐标 _z，;u = u，;v = v，;w = i’半反射法x= X，；y = y，；z = 2* 壁面 Z 方向的坐标-Z，-,U= a*^keTIm ； ν= a*yjkBT/m ；w= a*yJkBT/m上述公式中，a为-1 1之间的随机数，kB为波尔兹曼常数，T为当前温度，m为 质量，X表示粒子在X方向坐标上的反射位置坐标，y表示粒子在y方向坐标上的反射位置 坐标，ζ表示粒子在ζ方向坐标上的反射位置坐标，u为粒子在χ方向坐标上的反射位置速 度，ν为粒子在y方向坐标上的反射位置速度，w为粒子在ζ方向坐标上的反射位置速度。 本文中的“*”为乘号，“_”为减号。从采用随机反射法、全反射法和半反射法的模拟在特定条件下平行平板中的氩 流体流动中碰壁的情况来看，当采用随机反射法模拟时，参见图1(b)粒子A(X(I，y0, z0, u0, v0, w0)在运动过程中超出模拟区域边界时，粒子不会沿着既定线路到达下一步指定位 置B’(χ’，y’，ζ’，U’，V’，W’)，而会沿着粒子反射到B’向壁面做垂线的垂点B，χ = χ’ ；y =y’ ；ζ =壁面ζ方向的坐标；速度大小由壁面温度确定，《=口*彳‘7>2 ,ν= a*^kJjm、 W= a*^bTIm ；方向随机，；当采全反射法模拟式，参见图1(a)，粒子A (X(1，％，、，％，％，％) 在运动过程中超出模拟区域边界时，粒子不会沿着既定线路到达下一步指定位置B’(χ’， y，，ζ，，U，，V，，W，)，粒子会在壁面发生全反射后运动到B (x, y, ζ, u, ν, w)点，B (χ, y, ζ, u, v,w)与B’(X’，y’，Z’，U’，V’，w’)点关于壁面对称，与壁面平行坐标方向的速度矢量保持 不变，与壁面垂直坐标方向的速度大小保持不变，速度方向反向；当采用半反射法模拟时， 粒子A(X(I，y0, z0, u0, v0, w0)在运动过程中超出模拟区域边界时，粒子不会沿着既定线路到 达下一步指定位置B’(χ’，y’，ζ’，U’，V’，W’)，而会在虚拟的壁面边界上发生反射，到达反 射位置B (x, y, z,u, ν, w)，其反射后位置B既与粒子未发生反射时的指定位置B’(χ’，y’， z’，u’，v’，W’)有关，也与粒子发生反射时的壁面的温度有关，即χ = χ’ ；y = y’ ；z = 2* 壁面ζ方向的坐标-ζ， ’u= a*Hm -,ν= a*^keTlm ；w= a*^kBT/m，粒子反射后的速度 方向随机，从而得到粒子反射后的位置。由于粒子反射后的速度与三个方向的速度分量u， ν和w有关，因此，粒子反射后的速度大小为λΑ/2 +V2 + V ,而速度的方向为三个速度分量的 速度方向的合成，而速度分量是有正有负的，所以合成的速度方向也是各个角度都有的，从 而粒子反射后的速度方向是随机的。本说明书中，壁面χ方向的坐标、壁面y方向的坐标和 壁面ζ方向的坐标是指壁面在三维笛卡尔坐标系的坐标，由于在实际应用过程中三维坐标系的原点可以为0点，也可以不为0点，因此，壁面本身在三维坐标系中的坐标与三维坐标 系原点的设置有关。2)模拟结束后，统计了三种壁面边界条件模拟方法情况下粒子数密度分布的情 况，如附图4所示。比较结果说明采用半反射边界条件模拟法和采用全反射边界条件模拟法模拟得 到的结果100%吻合，具有很强的实用性。从而表明本发明方法在模拟粒子在壁面边界条件 下的运动状态时，不仅可以体现恒壁温条件和壁面不光滑的物理条件，而且也可以在分子 动力学模拟中真实反映粒子反射后的位置。半反射边界法与随机反射边界法和全反射边界法异同见表1和表2 表1半反射边界法与随机反射边界法和全反射边界法特性比较表
权利要求
分子动力学模拟中壁面边界的模拟方法，其特征在于，所述模拟方法包含以下步骤A)定义分子动力学模拟中粒子的数据结构，粒子在三维模型中由六个信息来表示粒子在分子动力学模拟通道中的运动状态，所述六个信息分别是x方向坐标、y方向坐标、z方向坐标、x方向速度、y方向速度和z方向速度；其中x方向坐标为壁面横坐标，y方向坐标为壁面纵坐标，z方向为与xy平面垂直的方向；B)确定粒子在通道中的当前位置A(x0，y0，z0，u0，v0，w0)，其中x0表示粒子在x方向上的当前位置坐标，y0表示粒子在y方向上的当前位置坐标，z0表示粒子在z方向上的当前位置坐标，u0表示粒子在x方向上的当前位置速度，v0表示粒子在y方向上的当前位置速度，w0表示粒子在z方向上的当前位置速度；C)根据粒子当前位置及粒子受到其它粒子的作用力信息，确定粒子的运动方向，并根据粒子当前位置及受力信息判断粒子下一步是否会运动到壁面边界以外位置；如果粒子没有运动到壁面边界以外位置，则粒子会沿着确定的运动方向运动，到达下一步指定位置B’(x’，y’，z’，u’，v’，w’)，其中x’为粒子在x方向上的指定位置坐标，y’表示粒子在y方向上的指定位置坐标，z’表示粒子在z方向上的指定位置坐标，u’表示粒子在x方向上的指定位置速度，v’表示粒子在y方向上的指定位置速度，w’表示粒子在z方向上的指定位置速度；D)若粒子运动到壁面边界以外位置，粒子将不能到达下一步指定的位置B’，粒子根据壁面的温度，经过壁面反射到达下一步反射位置B(x，y，z，u，v，w)，当x方向存在壁面，则其中y＝y’，z＝z’，x＝2*壁面x方向的坐标 x’，粒子的速度方向随机；当y方向存在壁面，则其中x＝x’，z＝z’，y＝2*壁面y方向的坐标 y’，粒子的速度方向随机；当z方向存在壁面，则其中x＝x’，y＝y’，z＝2*壁面z方向的坐标 z’，粒子的速度方向随机；上述公式中，a为 1～1之间的随机数，kB为波尔兹曼常数，T为当前温度，m为质量，x表示粒子在x方向上的反射位置坐标，y表示粒子在y方向上的反射位置坐标，z表示粒子在z方向上的反射位置坐标，u为粒子在x方向上的反射位置速度，v为粒子在y方向上的反射位置速度，w为粒子在z方向上的反射位置速度。FSA00000292915600011.tif,FSA00000292915600012.tif,FSA00000292915600013.tif,FSA00000292915600014.tif,FSA00000292915600015.tif,FSA00000292915600016.tif,FSA00000292915600017.tif,FSA00000292915600018.tif,FSA00000292915600019.tif
全文摘要
本发明公开提出一种分子动力学模拟中壁面边界的模拟方法，该方法综合了全反射法和随机反射法优点，能够真实反应出粒子碰壁后的运动状态，也称为半反射法，其实现步骤如下1)定义分子动力学模拟中粒子的数据结构；2)确定粒子在通道中的当前位置A(x0，y0，z0，u0，v0，w0)；3)根据粒子当前位置及粒子受到其它粒子的作用力信息，确定粒子的运动方向，并根据粒子当前位置及受力信息判断粒子下一步是否会运动到壁面边界以外位置；4)若粒子运动到壁面边界以外位置，经过壁面反射到达下一步反射位置B(x，y，z，u，v，w)，当x方向存在壁面，则其中y＝y’，z＝z’，x＝2*壁面x方向的坐标-x’，粒子的速度方向随机。本发明在模拟中可以体现出壁面恒温和壁面不光滑的物理条件；可以更加真实地反应出粒子反射后的位置；也可应用在不规则的纳米通道中，简化速度矢量的计算方法。
文档编号G06F17/50GK101944151SQ201010299639
公开日2011年1月12日 申请日期2010年9月30日 优先权日2010年9月30日
发明者丁玉栋, 冯洁, 叶丁丁, 吴睿, 廖强, 朱恂, 李俊, 王宏, 王永忠 申请人:重庆大学