《空泡流现象的非平衡分子动力学模拟和机理研究》由会员分享，可在线阅读，更多相关《空泡流现象的非平衡分子动力学模拟和机理研究（121页珍藏版）》请在金锄头文库上搜索。
1、浙江大学理学院 博士学位论文 空泡流现象的非平衡分子动力学模拟和机理研究 姓名：龚博致 申请学位级别：博士 专业：物理化学 指导教师：张秉坚 20081027 浙江人学博l j 学位论文 摘要 流体中物体超高速运动的微观机理是现代科技最前沿的研究领域之一，无论 在军事还是在民用技术中都具有极高的应用价值，特别是超空泡现象的研究，对 水下航行体的减阻效应和超高速鱼雷的设计至关重要，有着重要的科学价值。因 此本工作得到了国家自然科学基金项目超高速流体的汽液相变与界面现象研 究( 2 0 4 7 3 0 7 3 ) 的资助。 本论文运用非平衡分子动力学模拟方法，研究了微观体系条件下流体分子与 固体壁
2、面高速相对运动过程中流体流粒子的分布规律、状态性质、界面现象和摩 擦情况等，特别是对开放体系超空泡现象的计算机模拟的实现，对超空泡流动中 分子运动规律的观察和微观机理的探讨，为真实环境中物体超高速运动的研究提 供了重要的基础数据。 为了逐步深入研究不同体系不同条件下的流动，我们模拟了从简单窄孔到开 放体系，从简单分子到甲烷和水等真实流体，从单系统到复合系统等模型，应用 N E M D 方法计算了流体压力、密度、速度、界面张力、摩擦系数等物理化学性 质，通过对流场细节的观察，提出了一系列有关高速流动现象的形成机理，并与 宏观流体力学理论进行了关联。 为了验证分子动力学模拟程序和所建分子模型，本工
3、作首先利用所建立的 N E M D 模型，代入平衡流动条件，模拟了宏量条件下甲烷的P V T 关系，模拟结 果与实验测量值或文献模拟值都能较好地吻合，说明本工作的模拟方法和设置基 本正确。随后，对微孔无阻碍流进行了非平衡模拟，获得了符合宏观水动力学的 模拟结果，验证了模拟程序与模型在非平衡流动模拟中的可行性。 在上述研究的基础上，通过设定阻碍物模型，模拟计算了受限于窄缝孔中甲 烷流体在1 6 0 K 温度、不同力场作用以及不同阻碍物条件下的阻碍流流场细节。 模拟结果发现在高速条件下阻碍流会在流经阻碍物时形成高密区和低密区，同时 引发相应的压力变化。流速对低密区形成有重要影响；阻碍物前端形态会影
4、响流 体两相形成和自 f 端压力梯度，因此钝头阻碍物与其他流线形阻碍物相比将受到更 大的压差阻力，但对两相界面表面张力影响不大；增加阻碍物长度或外力场强度 会增大两相间表面张力。同时流动中观察到的摆动波现象符合R a y l e i g h T a y l o r 理论，表明该模型可以作为数值模拟的补充，研究宏观现象的微观机理。 浙江人学博I ：学位论文 为了将流体模型应用到更有实际意义的超空泡流动模拟中，本工作引入了新 的开放体系模型。通过将E M D 模拟体系与N E M D 模拟体系结合，模拟系统由 硅酸盐狭缝孔扩展到开放体系，壁面限制作用被消除，外加力场改为恒流速流动， 使得流体更为贴
5、近真实流动状态。通过引入局部空化数a 概念，模拟体系被划分 为多个独立方格，以观察超空泡的形成过程以及不同因素对超空泡现象的影响。 模拟结果表明，微观条件下，空泡流的传统判据空化数t r r c f 2 11 ) 这样L e n n a r d J o n e s 势能在截断半径处为零，由此消除了非连续性。 在计算作用力的时候，要对L e n n a r d J o n e s 作用势能求梯度。将式( 2 9 ) 以对比 单位表示如下： 浙江人学博l ：学位论文 , r , 4 l i r a 1 2 例 亿 式中，带星号的变量都表示对比单位。对上式求导，得到势能梯度： 警十2 ( 妒删( 去
6、) 亿 根据作用力和势能函数梯度的关系，可以求得各个方向上作用力的分量，如工方 向上有： 删2 侧警) _ 等( 专一o 5 专) 亿 Y 、z 方向的作用力分量可用相同的办法得到。 2 2 3M i s h i n 镶嵌原子势 2 0 世纪8 0 年代初。镶嵌原子法 3 7 ( E m b e d d e dA t o mM e t h o d ，E A M ) ， F i n n i s S i n c l a i r 方法 8 ，9 ( F S ) ，有效介质理论 1 0 ( E f f e c t i v eM e d i u mT h e o r y , E M T ) 以及G l
7、u e 势 1 1 等多体势能模型相继提出。它们考虑了原子在区域背景电子云密 度环境下的结合能，将系统总能量分解为如下统一的形式： = F ( 肛) + 去中( ) ( 2 1 5 ) 肛= 厂( ) ( 2 1 6 ) 其中肛为第i 个原子处的电子云密度，( 乃) 为原子在距离勺出产生的电子云 密度，镶嵌能函数F 表示为了将原子f 嵌入电子云密度肛的位置所需要的嵌入能， 函数鳏用中心对势形式。勺为原子i 、J 间距离。式中1 2 表示结合能为两原子 共有，避免重复计算影响体系总能量。镶嵌能函数F 的形式与背景电子云来源无 关，具有一般形式，可用于直接计算不同分子之间的势能。 目前基于E A
8、M 等多体势框架已发展了适合多种材料的原子势函数。每种原 子势函数都有一定的适用范围，如M i s h i n 于2 0 0 1 发展的C u 镶嵌原子势 6 】和2 0 0 3 年提出的T i A l 系统镶嵌原子势 7 】。这些方法虽然实现比传统U 势能模型复杂， 但计算量并没有显著提升。 2 2 4 不同种分子间相互作用势能 浙江人学博一l j 学位论文 对于不同种L e n n a r d J o n e s 流体分子间的相互作用，可以采用参数混合的办 法： 4 z = 丢( 磊- + 嘎：) ( 2 1 7 ) q ：= 石忑 求得相应的相互作用L e n n a r d J o n
9、 e s 势能参数。 在计算作用力时，将同种分子间的作用势能和异种分子问的作用势H , - 匕E 加和： 办= 确。+ 办： 兜= 办，+ 兜： ( 2 1 8 ) ( 2 1 9 ) 在这个过程中，要特别注意单位的统一。 由于L e n n a r d J o n e s 势能模型实现简单，计算速度快捷，且有足够的实验数 据进行理论支持，本工作中简单粒子模型均采用U 势能模型。 2 3N E M D 统计热力学实现 2 3 1 解牛顿运动方程 牛顿运动方程求解中，需要利用有限积分法( F i n i t eD i f f e r e n c eM e t h o d ) ，将模 拟时间段分为
10、固定间隔为万f 的小段，从粒子间交互作用力的向量总和，计算出 在t 时问每一个粒子所受的合力。 V e r l e t 算法 1 2 假设在该时间间隔内作用力恒定，则可以利用式( 2 2 ) 求得粒子的加速度 在分子动力学模拟中，常用的解运动方程的方法是V e r l e t 算法。这个算法不 仅是最简单的，通常也是最好的 1 3 。它是通过对t 时刻某粒子的坐标和下 一时刻坐标的差分的泰勒展开获得的 1 4 】： 心+ A 归心) 州必H 掣+ 等j = ：+ o ( 血a ) 2 mj ! f 2 2 0 ) 式中o ( 6 t q ) 表示截断误差的数量级。相似地，可以得到： 心地一叫倘
11、+ 等以等+ O ( A t 4 m ) ( 2 2 1 ) Zj ! ，1 、 浙江人学博l j 学位论义 上面两式相加，可得： r ( f + ，) 十，| O 一，) ：2 ，_ O ) + 翌2 f 2 + 2 0 ( A t 4 ) m 省略高次项，得到： 衍出川一( f 一f ) + 等 根据轨迹的牛顿力学知识可以得到速度的表达式： V ( r ) = r ( t + A 1 t ) - 万r ( 一t - A t ) V e r l e t 算法实现的程序片断如下： S U M V 2 = 0 ： F o r ( I = 0 ；I T E M P E = S U M V 2 3
12、N P R ： 计算温度 注意到该算法中速度并不是直接由上一时刻的速度得到的，也没有用速 度来计算新时刻的位置。速度的这一表达式的精度只有& 2 左右。为了使速 度的计算更加精确，可以对V e r l e t 算法进行改进，例如“蛙跳法”( L e a p - F r o g ) ： r O + A t ) = r ( t ) + v ( t + A t 2 ) A t ( 2 2 5 ) 浙江人学博I j 学位论文 1 ，( f + A t 2 ) ：1 ，( f _ A t 2 ) + 盟f m ( 2 2 6 ) 这罩，定义丛垒：口( f ) 为加速度。这样，在计算的时候存储的数据就是
13、m 当前时刻的位置r ( t ) 和加速度口( f ) ，首先计算的是速度项，当前时刻的速度 可以用下面的式子来计算： V ( 归j 1 ( ，( f + A t 2 ) + V ( f 一垃2 ) ) ( 2 2 8 ) 从该算法可以看出，位置的计算总是比速度的计算提前一个时刻，这样 新时刻加速度的可以先求出来，但是动能和势能就不能同时计算出来。为了 能够同时得到某一时刻粒子的位置、速度和加速度，可以采用另一种改进的 V e r l e t 算法，如下面方程所示： 心地一州归+ 扣逾2 ( 2 2 9 ) V ( H 缸2 ) = 州+ 吉) 岔 y ( H = y ( H a t 2 )
14、+ 吉础+ & ) 出 ( 2 3 0 ) ( 2 3 1 ) 这样一来速度就可以直接从上一时刻的对应值得到，其精度可以达到 & 4 。还有一种更精确的改进的V e r l e t 算法，称作B e e m a n 算法，如下面所示： ，( f + & ) = ，( f ) + V ( f ) + ；口( f ) f 2 一否1 口。一垃) 出2 ( 2 3 2 ) V ( f + 岔) = V ( ，) + 三口( f ) 垃一丢口( ，一出) 缸 ( 2 3 3 ) B e e m a n 算法能够更好地给出能量的计算结果，但是该算法的形式稍微复 杂一些，可能会需要更多的计算时间。 牛顿动
15、力学方程本身是时间可逆的，因此解方程的方法最好也要能够符 合这一条件。V e r l e t 算法理论上符合时间可逆要求，但由于计算机精度的截 断影响，实际的计算模拟中不可能真J 下实现时间可逆。同时运动方程的积分 浙江人学博I j 学位论义 算法最好是能在长时f n 范围内重现粒子的运动轨迹，但分子动力学模拟体系 中，粒子在相空间穿过的轨迹十分依赖于初始设置。即使初始设置相差十分 微小，随着模拟时间的推移，轨迹也会出现显著的差异。在积分算法中，无 论计算精度有多高，这种差异总是存在的，即所谓的“李雅普诺夫不稳定性”。 不过对于分子模拟的来说，这种差异并不会对结果带来大的影响，因为模拟 值都是
16、在模拟时间上的统计平均值，只要模拟时间足够长，这种差异可以被 修正。 预测校正法 预测校正法是另一种解牛顿运动方程的积分算法。它通过三个步骤实现： ( a ) 预测。根据t 时刻i 粒子的坐标位置，速度，加速度，通过泰勒展丌式计 算t + d t 时刻相应物理量的值。如下： ，( ) = r ( r ) + V ( ，) 班+ 扣) 肌扣) 矿+ 去c ( ，) 舢( 2 3 4 ) V ( t + d O = V ( ，) + 口( ，) 历+ 丢6 ( r ) 魂2 + 吉c ( z ) 疵3 + ( 2 3 5 ) 口( ，+ 出) = 口( f ) + 6 ( z ) 衍+ 吉c ( f ) 加2 + ( 2 3 6 ) b ( t + d O = + c ( f ) 坊+ ( 2 3 7 ) ( b )
