水系统气泡形成过程的分子动力学模拟

6

（9）
其中，σ 是平衡分离尺寸（势能第一次达到0时分子间的距离）,ε 势阱深度。σ 和 ε 分别是长度和能量的衡量参数。为了节约计算时间，引入截断尺寸D，当分子间距大于 D的时候，这两个分子间的相互作用力就忽略为0。所以，分子间的受力可以用下式来表 示：
24 dij (rij ) F (rij ) rij drij
*基金项目：本项目受到国家自然科学基金面上项目 (NSFCNo.50476049）的资助
Lennard-Jones 分子组成液态氩， 用虚拟常温模型的三层 harmonic 分子表示每一个固体 表面。固体分子和氩分子之间的相互作用势能用 Lennard-Jones 势能表示，结果显示空 穴核在下表面上随机地出现和消失，在许多点能够长大到某一稳定尺寸。但是在汽相非 均相成核的情况下，没有能够成功计算成核速率。 Michael[8] 采用 ST2、LSR 两种水分子相互作用势能函数对水的液-汽界面进行了模 拟，液-汽界面过渡区内的密度变化是平滑的，液气界面过渡区的宽度在密度变化 10%～ 90%的范围内是 0.345nm，与用 X 射线测量的实验值 0.330nm 一致，这个宽度相当于半个 到一个分子直径。径向密度分布随分子簇大小的变化很小，用不同的分子间相互作用势 能模拟得到的径向密度分布几乎是相同的。 Dang[9] 提出了一个刚性四点可极化水分子模型， 将双体势能看成 L-J 势能和 Coulomb 相互作用的和。用此函数对含 550 个水分子、边长 2.6nm，具有完全周期性边界条件的 立方体模拟体系进行了模拟。密度分布表明液-汽界面在 298K 时其厚度为 0.32nm，表面 张力与实验值一致，界面处的分子偶极距接近气相值。 本项工作的前期已经对非极性分子系统（氩）中的微气泡形成过程进行了分子动力 学模拟，获得了微观气泡形成的是个基本过程以及过热度对该过程的影响。为了更接近 于实际的相变过程，本文同样运用平衡态分子动力学方法对含有极性水分子的微正则系 统进行了模拟，验证了极性水分子气泡形成过程的四个阶段，获得温度对微气泡形成过 程的影响以及气泡形成过程产生不稳定的原因。
2 r ij
12 6 , rij D r ij 0, rij >D
（10）
1.3 运动方程有限差分法 对于五阶预估－校正法需要在模拟开始时给出每个粒子的初速度和初始位置。一般 来讲，粒子的初速度和初始位置都是由随机数产生的。这种方法产生的某些粒子位置很 可能非常接近，这样这些粒子间的距离几乎为 0 从而使计算结果溢出。本文采用的方法 是：初速度取 0－1 之间的随机数；以立方面心晶格的方式给出每个分子的初始位置。 1.4 热力学特性参数计算 平衡态下的热力学参数通过统计平均得到，如总能量U、压力P的统计式如下：
12 6 qq 12 R1 , R2 4 OO OO ห้องสมุดไป่ตู้ OO S ( R12 ) i j R R i j 4 0 r ij 12 12
HOH 2 cos 1 (1/ 3) 109.47
摘要：运用平衡态分子动力学理论对含有极性水分子系统的微气泡的形成过程进行模拟。采用五阶预 估－校正有限差分法对每个分子的牛顿运动方程进行求解，该方法能够较好的满足能量守恒特性。通 过计算验证模拟的正确性，统计出了系统的瞬时压力和温度，得到了气泡生长过程中各相分子的分布 形貌，并且分析了温度对气泡生长过程的影响以及不稳定性状态分析。计算得到水系统的势能与现有 的研究成果对比，符合很好，验证了模拟的正确性。极性水分子系统中旗袍的形成基本过程与氩系统 相似，只是温度和无量纲压力比氩系统波动更大，这可能是由于水分子的极性导致的。当过热度大于 14K 时气泡破碎导致界面发生不稳，不同工质发生界面不稳现象所需过热度不同。
，
T k BT ， L L ，
*
*
*
3
m ， P P ， U U 。
* 3 *
所有分子放在体积为 L3 的立方体中，计算中将质量中心定于正方体的中心，这样 便于确定气泡的形貌以及相关的特性参数。同时，六个表面均采用周期性边界条件以保 证模拟空间中的粒子数恒定。 为 了 与 现 有 结 果 对 比 ， 验 证 计 算 选 取 水 分 子 数 为 256 ， σ=0.2791 nm ， =5.85 10-21J ，晶格所占空间由水的密度 1.0g/cm3 导出。 对于两种模型，相同的参数设置如下：各分子起始速度按随机分布取样，然后对速 度进行标度，以确保体系总动量为零，初始位形按面心立方晶格分布；方程求解中采用 周期性边界条件；采用五阶预估-效正进行求解。模拟过程中位能截断采用球形截断法， 截断半径为2.5σ。模拟温度选取T =298 K。模拟过程中选取时间步长为0.5fs ，这是由于 水分子的转动惯量较小，会发生高速旋转，必须采用较小的时间步长以保证运动方程求 解的稳定性。每次模拟的总步数为9000步，体系趋衡后统计各种平衡性质和结构性质。 模拟过程中各参量均无量纲化。 为了验证本文模拟结果的正确性以及与实验值的差距，将势能的计算结果与文献值 和实验值进行了比较，见表 2。从表中可以看出，本文通过两种模型计算的水的势能与 其他研究者的模拟及实验值非常接近， 从而验证了本文计算的正确性， 同时证明了在 L-J 模型中考虑到水的极性而采用了有效参数来代替的这种方法是可行的。


用力采用经典的Lennard-Jones (L-J) 12-6势能函数来表示：


（8） j 个分子间的
距离； ai 和 mi 分别是第i个分子的加速度和质量。 对于惰性气体氩来讲，分子间的作
ij
4 rij

12
r ij
3 U NK BT 2
P
2 结果与讨论
(r )
i j i ij
(11)
N 1 K BT 3 3 L 3L

i j i
(rij ) rij
rij
（12）
2.1 模型验证 对于 L-J 势能函数模型，需要将分子间的偶极、诱导偶极以及多体作用考虑在势能 函数中，所以采用二体位能模型主要原则是将分子参数代换为“有效参数” ，则两水分 子之间的相互作用势可以用式（9）表示。为了避免原子量级的参数间的计算，同时为了 简化运动方程， 模拟中的参数都以长度 σ 、 能量 ε 和质量 m 为基本参数进行无量纲化。 无量纲时间、 温度、 长度、 数密度、 压力、 总能量分别表示如下： t * t / t / m 2
(7)
ε OO
*
电量 e=1.60219 10-19 C
1.2 基本方程 分子动力学模拟主要是基于经典的牛顿第二定律，系统中每个分子的运动方程可以 表示成：
d2r Fi m ia i m i 2i dt 其中， Fi 是作用在第i个分子上的总作用力； rij 是第i个分子和第
关键词：分子动力学；气泡；当地密度
0 前言
对于沸腾机理研究，特别是对于气泡动力学的研究，现在已经形成了经典的气泡动 力学理论，总结了一系列经验公式，得到了广泛的运用。但是当研究者致力于对沸腾的 非线性机理时发现，尽量精细试验系统、提高测试手段或是再精确的计算，都很难得到 明确的沸腾机理。这主要是由于气泡核化现象在时间和空间尺度上非常小，核化过程中 －9 －10 [1] 形成的核子临界直径量级在 10 m，形成临界核子的时间也仅仅在 10 s 的量级 ，想要 在实验过程中达到这个水平几乎是不可能的，这就促使了从微观上对沸腾过程进行更进 一步的研究。近年来运用分子动力学模拟沸腾相关问题的研究逐渐增多。特别是研究的 工质，从简单的惰性气体如氩，逐渐发展到现实中常用的工质，如制冷剂 R134a 以及最 常用的水。 对于均相成核的研究，Yasuoka 和 Matsumoto[2-4] 对 Lennard-Jones 和水分子的核化 过程的非平衡动力学模拟进行了验证， 对氩工质进行了均相成核模拟。 他们在模拟中用 5000 个 Lennard-Jones 分子与 5000 个携带气体分子的“软核”混合，并与保持恒温的 Nose-Hoover 恒温器相连。模拟得到的核化速率是经典核化理论估计的 7 倍。另外采用 TIP4P 势能对水进行了类似的模拟，计算的核化速率是经典核化理论的计算值的一半。 [5-6] Ikeshoji 等 对 Lennard-Jones 分子的核化过程进行了相同的模拟，在实验质谱仪上 进行了大量的观测，采用 65526 个分子进行大尺度模拟，与 Yasuoka 和 Matsumoto 的模 拟时间 3.9ns 相比，他们模拟蒸发过程的时间扩大到 26.4ns(氩)。 在许多热传递问题中，非均相成核实际上比均相成核重要。Maruyama[7] 用 5488 个
SPC 或 SPC/E 的势能表达式为：
(1)
(2)
12 6 qq 12 R1 , R2 4 OO OO OO i j R R 4 0 rij 12 12 i j
HOH 2 cos 1 (1/ 3) 109.47
(3) (4)
TIP4P 的势能函数表达式与 SPC 或 SPC/E 势能相同， 只是 HOH 不同：
HOH 2 cos 1 (1/ 3) 104.52
(5)
CC 势能函数表达式较为复杂，但是相对最准确：
12 R1 , R2
i j
qi q j 4 0 rij
HOH 2 cos 1 (1/ 3) 104.52
上述四种是最常用的水的势能函数，表达式中的参数取值如表 1 所示。 表 1 常用的水的势能函数参数取值 项目 rOH ∠HOH σOO （ 10-21 ） rOM qH * qM 单位 [nm] [0] [nm] [J] [nm] [C] [C] ST2 0.1 109.47 0.31 0.52605 0.08 0.2357e -0.2357e SPC/E 0.1 109.47 0.3166 1.0797 0 0.4238e -0.8476e TIP4P 0.09572 104.52 0.3154 1.0772 0.015 0.52e -1.04e CC 0.09572 104.52 N/A N/A 0.024994 0.18559e -0.37118e