全文HTML
--> --> -->Einstein的布朗动力学理论指出颗粒的随机运动是周围流体分子热运动的结果. 颗粒在周围流体分子的碰撞下产生随机运动的同时, 又受到流体分子的阻力, 由碰撞所产生的随机力和流体分子对颗粒的摩擦阻力通过涨落耗散定理约束. 颗粒的热涨落和摩擦力的关系为




图 1 球形颗粒在不同边界条件下的平动摩擦系数, 颗粒向右运动 (a)黏性边界条件; (b)光滑边界条件Figure1. Translational friction coefficients of spherical particles under different boundary conditions and the particle move to the right: (a) Translational friction coefficient for stick boundary; (b) translational friction coefficient for slip boundary.
本文采用分子动力学模拟方法建立了单个Cu纳米颗粒在水溶液中扩散的全原子模型, 在原子水平上研究流体分子与纳米颗粒之间的相互作用, 在计算能力能够承受的范围内分析比较了4种不同半径的纳米颗粒在水溶液中的扩散特性. 采用单颗粒追踪算法[17-19]并结合最小二乘法对颗粒的平转动扩散系数进行拟合, 解决了扩散系数需要对大量颗粒系综平均的问题. 将摩擦系数的计算结果与黏性边界和滑移边界条件下的结果进行比较, 讨论了颗粒表面的溶剂吸附对颗粒水动力学半径的影响. 胶体系统的粗粒化模拟过程中, 摩擦系数的选择直接决定布朗力和黏滞力的大小, 本文的结果可以为摩擦因子的选择提供支持, 而不必在黏性边界与滑移边界条件之间做出选择, 为布朗动力学和耗散动力学模型的完善提供理论依据.
2.1.全原子模型的构建
本文建立了不同大小的单个Cu纳米颗粒在纯水溶液中扩散的全原子模型(图2(a)), Cu颗粒半径分别为0.5, 1.0, 1.5和2.0 nm, 颗粒中原子按照FCC的晶格结构排列. 初始时刻颗粒位于模拟域的中心, 模拟域的大小和总原子数目需要根据颗粒的大小进行调整, 如表1所列. 边界条件均采用周期性边界条件, 从而保证颗粒的运动为自由布朗运动. 模拟开始时颗粒的初始速度符合麦克斯韦分布, 模拟时间步长为1 fs, 系统温度均为298 K. 由于溶剂动力学黏度对颗粒摩擦系数的计算有直接的影响, 而不同水分子模型的动力学黏度相差较大, 本文中水分子的构建选用与实验符合较好的刚性TIP4P/2005模型[20]. 除了两个氢原子和一个氧原子, 该模型在H—O—H的键角平分线上增加了一个无质量的作用位点M, 水分子结构如图2(b)所示, 其中键长b = 0.9572 ?, 键角θ = 104.52°, O和H的电荷分别为–1.1128和+0.5564 eV, 水分子的半径d = 2 ?, OM = 0.1546 ?. 水分子的键长和键角采用Shake算法固定, 因此键合势为零, 长程静电力采用PPPM算法处理[21]. 原子间的非键合相互作用包括范德瓦耳斯相互作用和静电相互作用, 计算公式如下:
图 2 (a) Cu-H2O纳米流体的物理模型; (b) TIP4P/2005水分子结构示意图Figure2. (a) Physical model of Cu-H2O nanofluid; (b) schematic diagram of TIP4P/2005 water.
| 颗粒半 径/nm | Cu原 子数 | H2O分 子数 | 总原 子数 | 模拟域大小 Lx × Ly × Lz /?3 |
| 0.5 | 43 | 971 | 2956 | 30 × 30 × 30 |
| 1.0 | 351 | 7380 | 23841 | 60 × 60 × 60 |
| 1.5 | 1199 | 7459 | 23576 | 60 × 60 × 60 |
| 2.0 | 2928 | 14381 | 45972 | 75 × 75 × 75 |
表1各工况下模拟域大小与原子数目
Table1.The size of the simulation domain and number of atoms of different working conditions.
颗粒与水分子的相互作用只考虑Cu原子与氧原子的范德瓦耳斯作用, Cu-H2O纳米流体相互作用的势能参数如表2所列, 计算过程不考虑氢原子和氧原子及铜原子的范德瓦耳斯作用. 颗粒中原子间的相互作用采用适用于金属和合金的多体势(EAM)[22], 这种EAM势在描述金属原子间的相互作用时比对势更精确, 晶体的总势能表示为
| 参数 | 取值 |
| εo—o/eV | 0.00803 |
| σo—o/? | 3.1589 |
| εCu—o /eV | 0.052 |
| σCu—o /? | 2.743 |
表2Cu-H2O纳米流体势能参数
Table2.Potential energy parameter of Cu-H2O nanofluid.
系统能量最小化后, 在NPT系综下弛豫平衡, 采用Nose-Hoover方法控温控压. 目标温度和压力分别为298 K和 1 bar (1 bar = 105 Pa), 时间步长为1 fs, 经过20 ps的弛豫后, 压力和温度趋于稳定. 体系稳定后, 在NVE系综下使颗粒做自由布朗运动, 总步数为6000 万步, 对应总时长为60 ns. 上述纳米流体模型的建立以及模型验证工作均采用LAMMPS软件包完成.
2
2.2.水分子模型的验证
为了得到稳定的水分子体系, 需要确保系统在模拟过程中有足够的时间达到平衡. 选取随机分布的4096个水分子, 大小为49.6 ? × 49.6 ? × 49.6 ?的模拟域, 对弛豫过程中体系的温度和总能进行监测, 温度和总能随平衡时间的变化如图3(a)和图3(b)所示. 可以看出, 体系温度和总能量在300 fs左右时已经达到平衡, 总体积和密度亦趋于稳定. 由于水分子的动力学黏度对摩擦系数的计算至关重要, 本文采用green-kubo公式(5)式计算水的黏度, 该方法将流体的剪切黏度和压力张量中非对角元素的自相关函数联系起来, 通过计算压力张量的非主对角元素pxy, pxz, pyz的平均值得到TIP4P/2005刚性水分子模型的剪切黏度[23],
图 3 水分子模型的验证 (a) NPT系综下体系的温度随时间变化; (b) NPT系综下体系的总能量随时间变化Figure3. Validation of water molecular model: (a) Curve of temperature over time; (b) total energy of the system over time.
NPT系综下平衡之后, 在NVT系综下计算体系的黏度, 经过2 ns的模拟, 得到水分子的黏度为0.825 × 10–3 Pa·s, 该黏度与实验结果的误差约为7%[24], 其他水分子模型(TIP3P, TIP4P, SPE/C)误差与实验结果在18%—64%之间[23], 因此TIP4P/2005水分子模型的动力黏度已明显优于其他水分子模型.
3.1.单颗粒布朗运动平动扩散
在统计力学中, 均方位移是指颗粒随时间运动的距离相对于参考位置的度量, Einstein在布朗运动理论中指出颗粒移动距离的平方的均值与时间呈线性关系:

由于布朗运动的方向具有随机性, 它的轨迹是一条由在空间随机分布的点依次连接所构成的曲线, 该曲线处处连续且处处不可导[25]. 随着颗粒尺寸的增加, 体系内的总原子数因模拟域的增大而成倍增长, 受到计算资源的限制, 所有体系中只放置一个颗粒. 颗粒在总模拟时长内轨迹如图4所示, 蓝色代表轨迹的起点, 红色表示轨迹的终点, 模拟温度T = 298 K, 颗粒半径a = 1 nm. 从图4可以看出, 颗粒的运动方向并没有发生明显的变化, 但是轨迹同样满足处处连续且处处不可导的特征, 这是因为颗粒在运动过程中并未与其他颗粒发生碰撞. 对于单颗粒布朗运动轨迹, 采用传统的计算方法得到的扩散系数随计算时长的增加而增加, 为了准确地计算单个颗粒的均方位移, 本文采用单颗粒追踪方法并结合最小二乘法对颗粒的平动扩散系数进行拟合. 下面以半径a = 1 nm的颗粒为例说明扩散系数的拟合过程.
图 4 单颗粒布朗运动轨迹Figure4. Brownian motion trajectory of a single particle.
1)确定最佳关联时间. 颗粒的轨迹为含有N = 60000个数据点的连续曲线, 曲线的每一段的时间间隔为1 ps. 将该轨迹按不同的关联时间(时间间隔)Nt分成Ns = N/Nt份, 对每一段轨迹, 根据下式计算颗粒的均方位移:
根据所求均方位移进而确定扩散系数, 然后对每一个关联时间下的Ns个扩散系数进行直方图分析. 如图5所示, 当关联时间较小时, 扩散系数呈正偏态分布, 关联时间很大时, 扩散系数呈负偏态分布, 一个呈正态分布的关联时间即是最佳关联时间[26].
图 5 半径a = 1 nm的Cu颗粒在不同关联时间下扩散系数分布Figure5. Diffusion coefficient distribution of Cu particles with radius a =1 nm at different correlation time.
2)确定最佳拟合点. 从(7)式可以看出, MSD曲线的第一个点代表每一段位移的平均, 而最后一个数据点没有被平均. 计算扩散系数时统计点过多或者过少都会对计算结果带来不利的影响, 因此本文计算了每一段轨迹采用不同拟合点时扩散系数D的变异系数(二阶矩/一阶矩), 变异系数最低的拟合点数目即为最佳拟合点[27]. 图6为关联时间t = 1363 ps时得到的变异系数随拟合点变化的曲线, 可以看出, 拟合点数目在20左右时变异系数最低.
图 6 半径a = 1 nm的Cu纳米颗粒在采用不同拟合点拟合时扩散系数的变异系数, 关联时间Nt = 1363 psFigure6. Variation coefficient of the diffusion coefficient of Cu nanoparticles with radius a =1 nm when different fitting points were used for fitting, correlation time Nt = 1363 ps.
将最佳拟合点下的扩散系数取平均可以得到半径a = 1 nm的Cu颗粒的扩散系数, 采用同样的方法可求得a = 0.5, 1.5, 2.0 nm时颗粒的扩散系数. 图7(a)表示4种颗粒的平动扩散系数, 红色曲线和蓝色曲线分别表示黏性边界条件和滑移边界条件下的预测结果, 从图7(a)可以看出, 平动扩散系数均在预测值之间. 平动扩散的摩擦因子可由

图 7 纳米颗粒的平动扩散特性 (a)不同尺寸Cu纳米颗粒的平动扩散系数; (b)不同尺寸Cu纳米颗粒的平动摩擦因子Figure7. Translational diffusion characteristics of nanoparticles: (a) Translational diffusion coefficients of Cu nanoparticles with different sizes; (b) translational friction factors of Cu nanoparticles of different sizes.
2
3.2.单颗粒布朗运动的旋转扩散
与平动扩散均方位移的定义类似, 颗粒的均方旋转角度可按下式定义:


图 8 颗粒旋转扩散系数计算的物理模型Figure8. Physical model for particle rotation diffusion coefficient calculation.
不同于平动扩散系数, 旋转扩散系数的拟合对关联时间的选择不敏感, 因此可以直接对颗粒的MSR曲线进行拟合. 颗粒的均方旋转角度如图9所示, MSR曲线随时间单调递增, 对该曲线进行线性拟合即得到颗粒的旋转扩散系数. 图10(a)表示颗粒的旋转扩散系数随颗粒尺寸的变化, 红色曲线代表由Stokes-Einstein关系得到的颗粒在黏性边界条件下的扩散系数, 从图10(a)可以看出, 随颗粒尺寸的增加, 扩散系数逐渐趋近于黏性边界调下的结果. 图10(b)给出了旋转摩擦因子的大小, 旋转摩擦因子β可由

图 9 不同半径的Cu纳米颗粒的均方旋转角度Figure9. Mean-square rotation angles of Cu nanoparticles with different sizes.
图 10 颗粒旋转扩散特性 (a)颗粒的转动扩散系数; (b) 颗粒转动扩散摩擦因子Figure10. Rotational diffusion characteristics of particles: (a) Rotational diffusion coefficients of particle; (b) rotational friction factors of particle.
2
3.3.纳米颗粒表面溶剂分子的吸附特性
在滑移边界条件假设下, 纳米颗粒运动过程中表面不会携带流体分子. 然而考虑到纳米颗粒的粗糙度和固-液相互作用, 颗粒周围往往会吸附溶剂分子, 这部分溶剂分子与纳米颗粒一起做随机运动, 因此吸附的水分子层会对颗粒的有效水动力学半径产生影响. 本文统计了纳米颗粒周围水分子的分布[28], 如图11(a)所示, 通过统计距离纳米颗粒表面为r, 厚度为dr的球壳内的原子数目, 得到了水分子的径向分布函数(RDF), RDF可用下式来计算:
图 11 纳米颗粒周围水分子的RDF (a) RDF计算的物理模型; (b)不同尺寸的颗粒周围水分子的RDF曲线Figure11. Radial distribution function (RDF) of water molecules around nanoparticles: (a) Physical model for RDF calculation; (b) RDF curves of water molecules around particles of different sizes.



图11(b)为纳米颗粒周围原子数目的径向分布函数. 可以看出, 不同尺寸的纳米颗粒周围原子数分布均在1.7 ? (略小于水分子半径)左右时达到峰值, 峰值越大, 代表溶剂分子的吸附程度越紧密, 曲线中第一个峰的高度随着颗粒尺寸的增加依次降低; 第二个峰值出现在距离表面4.3 ?处(略大于一个水分子的直径), 且峰的高度没有明显差别, 之后水分子数密度逐渐趋于标准水分子密度. 因此, 在距离颗粒表面一个水分子直径的区域内水分子排列呈现结构化, 导致颗粒周围水分子的密度大于水分子的平均密度, 固体和水分子之间悬空的氢键是导致水分子在颗粒表面产生吸附的原因[12], 氢键存在于水分子中, 但水分子不能和金属原子之间形成氢键, 因此水分子趋于向颗粒表面靠近. 水分子在颗粒表面的吸附增大了颗粒的有效半径, 而且颗粒尺寸越小, 表面水分子层相对于颗粒的尺寸越大, 对颗粒水动力学半径的影响越大, 这是导致在转动摩擦因子的计算中, 尺寸a = 0.5和1.0 nm的颗粒摩擦因子向非滑移边界移动的原因.
