不均匀环境中二聚体微纳米马达的动力学性质

金英俊，董康杰，温正城

(杭州电子科技大学理学院，浙江 杭州 310018)

0 引 言

自驱动微纳米马达是一种尺寸在微纳米量级的微小器件，可以将其他形式的能量转化为动能，实现自主运动[1]。目前，自驱动微纳米马达的探索尚处于早期阶段，譬如Paxton等[2]研究以H2O2为燃料的Au/Pt双元金属棒马达在H2O2中的运动；Gao等[3]先后制备了以H2O2为驱动力的Al/Ga和Mg微米马达。马达动力学性质的研究往往在较简单的均匀环境中进行，主要关注的是马达的运动速度、驱动机理、功能应用等方面，燃料时空动力学特性引起马达运动性质变化的相关研究甚少。

众所周知，生物系统中的环境绝大部分是不均匀耗散体系，展示复杂的时空动力学现象。比如，在不均匀的环境中，细菌往往趋向有营养物质区域，远离有毒性物质区域，展示趋利避害的重要性质。随着相关研究的深入，不均匀环境对马达动力学的影响引起研究者关注[4-5]。本文通过构造燃料分布不均匀环境，运用混合的分子动力学和多粒子碰撞动力学(Molecular Dynamics-Multi-Particle Collision Dynamics，MD-MPC)方法，探究自驱动二聚体微纳米马达的动力学性质，以及马达间相互作用对其形成团簇的影响。

1 系统模型与仿真

本文主要针对自驱动二聚体微纳米马达在不均匀环境中的动力学性质和群集效应展开研究。建立尺寸为V=Lx×Ly×Lz=100×100×20的三维系统进行仿真模拟。系统中，随机放置大量的燃料粒子F以及不具有化学活性的惰性溶液粒子S，粒子的总数目NT=NF+NS，其中NF为燃料粒子F的数目，NS为惰性粒子S的数目。同时，系统中放入ND个二聚体微纳米马达，每个二聚体微纳米马达由催化球C和非催化球N组成。当随机运动的溶剂粒子碰到马达的催化球C时，C表面发生不可逆的催化反应，F+C→P+C，将燃料粒子F转化为产物粒子P；而在非催化球N处，没有催化反应发生。二聚体微纳米马达的构型如图1所示，主要包括催化球C和非催化球N。

图1 二聚体马达构型

将溶液中的F粒子、S粒子以及反应生成的P粒子作为质量均为m的质点。为了限制马达在Z方向上的运动，采用9-3兰纳-琼斯(Lennard-Jones，LJ)排斥势进行描述，具体表达式如下：

(1)

式中，r为2个粒子之间的距离，εW表示马达与上下墙壁之间相互作用的能量参数，σW表示马达与上下墙壁之间相互作用的距离参数。令σW=LZ/2，其中LZ表示系统在Z方向上的尺寸，由于这个势的存在，马达只能在XY平面上做准二维运动[6]。因为壁面的存在，系统中的溶液点粒子不断通过回弹碰撞与坐标Z=0和Z=LZ处的壁相互作用，与壁面碰撞后，小球的速度不变，并沿碰撞前的镜面对称方向弹回。对于含有多达1023量级的微观粒子的宏观系统，通过对整个微观系统的演化来获得宏观物理量并不容易。为了获得较为准确的物理量，在系统的X方向和Y方向上，采用周期性边界条件来正确描述二聚体马达的动力学性质。

对于溶液中的粒子与二聚体微马达之间的相互作用，采用12-6形式的LJ排斥势进行描述，具体表达式如下：

(2)

式中，λ为C或N，表示马达的C球或N球，μ为F，S或P，表示溶液中的粒子，rc=21/6σ为马达的截断距离，θ(x)为Heaviside函数。当r≤rc时，可以通过式(2)计算得出粒子与马达之间的作用势；当r>rc时，该作用势被截断[7]。除了产物粒子P对N球的能量参数不同外，其他所有溶液粒子与C球和N球的相互作用能量参数相同，即εCF=εCS=εNF=εNS=εF≠εP，其中εP=εNP。

对于含有多个二聚体微马达的系统，不同二聚体的单体之间也通过排斥LJ势相互作用，其中能量参数为εD，距离参数为σD=σλJ+σλK，J=1,…,ND和K=1,…,ND表示系统中的二聚体微马达。

采用MD-MPC方法对整个系统进行时间演化，MPC-MD方法由分子动力学(Molecular Dynamics，MD)和多粒子碰撞动力学(Multi-Particle Collision Dynamics，MPC)构成[8-9]。采用MD模拟溶液粒子运动的过程，通过牛顿运动方程演化溶液粒子与马达的相互作用、马达与马达之间的相互作用，其步长为τMD。采用MPC模拟溶液粒子之间的相互作用[10]。经过n步MD模拟后，采用MPC进行1次溶液粒子随机碰撞模拟，其步长τMPC=nτMD。用MPC模拟溶液粒子随机碰撞过程中，本文把系统分为尺寸a=1的众多立方体小格子，在每个单独的方格内，溶液粒子发生碰撞。标记为ξ号方格内的i号粒子的速度为：

(3)

数值模拟时，系统中所有的物理量均为无量纲LJ单位，分别基于能量εα，质量m0以及距离σd进行设置：r/σ→r，t(ε/(mσ2)1/2)→t和KBT/ε→T。具体参数如下：粒子F，S和P的质量均为m0=1。马达的C球、N球的直径分别为dC=2σC=2×2=4，dN=2σN=2×4=8，为了保证马达悬浮在溶液中，其质量随各自直径变化而变化。本文模拟中，mC=335，mN=2 688，每个马达的C球、N球之间的核间距固定为R=σC+σN+ξ=2+4+1=7。系统温度设定为KBT=1/6。采用MD模拟溶液粒子运动过程中，时间步长ΔτMD=0.01。采用MPC模拟溶液粒子之间相互作用过程中，纳米马达的旋转角为π/2，ΔτMPC=50×0.01=0.5，即n=50步进行1次碰撞，同时伴随1次能量转换。系统中的粒子总数NT=2.0×106，平均密度为10，即每个小方格子内含有10个溶液粒子。能量参数为εS=εF=5，εP=0.1，2个马达之间的εD=5。

初始化系统时，确保二聚体微马达位于系统中的一侧，放置规则如下：系统含有单个马达情形下，N球球心在X和Y方向上需满足10≤XN≤90,10≤YN≤90，C球相对于N球的角度范围为[0,2π]，同时满足5≤XC≤95，5≤YC≤95。当系统含有多个马达时，还需满足如下规则：任意2个马达α、β的N球之间的距离dNαNβ≥2×21/6σN，C球之间的距离dCαCβ≥2×21/6σC，C球和N球之间的距离dCαNβ≥2×21/6(σC+σN)。溶液粒子需满足条件如下：与任意马达C球的距离dSiCα≥21/6σC，与任意马达N球之间的距离dSiNα≥21/6σN。

2 模拟结果与分析

2.1 单个马达的自驱动行为

为了研究马达在不均匀环境下的群集效应与动力学性质，分别针对环境中只有惰性溶液粒子S和燃料粒子F这2种情形进行模拟实验，且只存在单个二聚体微马达。

当环境中只含惰性溶液粒子S时，溶液粒子S与马达的C球、N球之间没有化学反应，马达周围没有产生相应的梯度差，又粒子S和C球、N球之间的能量参数设置为εSC=εSN=5，所以，此时没有驱动力促使马达做定向运动，马达在系统中只做随机的布朗运动。当环境中只含燃料粒子F时，马达的C球表面发生催化反应生成产物粒子P的梯度场，虽然C球周围的P粒子场是均匀的，但在N球周围的P粒子梯度场是不均匀的，P粒子浓度沿着C-N的轴向降低，加上能量参数εF=5和εP=0.1(εF>εP)的不同，产物粒子P对二聚体N球产生的作用力小于其他溶液粒子对N球的推动，从而在N球上产生N-C指向的净力，使马达在系统中产自驱动。本文研究的二聚体微马达基于自扩散泳的机理，文献[12]通过模拟验证了自扩散泳机理，并且制备出这类二聚体微马达。

图2 马达的均方位移ΔL2随时间t变化的曲线

2.2 马达在单个圆形燃料区中的自驱动与群集

为了探究马达在不均匀环境下的行为，对单个马达(ND=1)在单个圆形燃料场中的行为进行模拟实验。以坐标(25,25)为圆心，r=25为半径建立1个圆形燃料场，区域内的溶液粒子均为具有化学活性的粒子F。初始时刻，系统中均匀分布大量惰性溶液粒子S，由于粒子S做无规则布朗运动，为了确保圆形燃料场内均为F粒子，将运动到区域边界的S粒子受到外界催化转变为粒子F。单个马达按照第1节所述规则进行放置。

图3 二聚体微马达的单个燃料场中运动及轴向速度

从图3(b)可以看出，t在0～810和16 634～17 475这2个时间段中，马达的轴向速度峰值远大于其他时间段的轴向速度，马达受到额外驱动力(FFN>FPN)的推动，在燃料区产生自驱动。

为了进一步探究马达在受限环境中的行为，在上述单个燃料场系统中放入15个二聚体微马达，即ND=15。除了按多个马达放置规则来放置马达外，还要确保惰性S溶液区和燃料区均有马达。为了观察更多二聚体微马达的构型，将马达的运行步数增加到3.0×107。图4为15个马达在单个燃料场系统中的瞬态构型。

图4 单个燃料场内，马达在不同时刻的瞬态构型

从图4可以看出，初始时刻时，15个微纳米马达不均匀地分布在系统中，大部分马达分布在惰性溶液粒子区，少数马达位于燃料区内。惰性溶液粒子区内的马达只做布朗运动，即缓慢进行自扩散，当扩散到S溶液粒子区与燃料区的交界处时，由催化反应产生的N球指向C球的力推动马达进入燃料区。当然，圆形燃料区内的马达也会跑出燃料区进入非燃料区。经过一段时间演化后，多个二聚体微马达在燃料区呈群聚状态。当燃料区中跑入多个马达后，由于化学反应，马达附近产生P粒子的梯度场，考虑到εF>εP，其他马达探测到该梯度场后，其N球向该马达的C球靠近，从而形成团簇。C球和N球的聚集排列展现出稳定的构型，跑出燃料区的马达以及原本在溶液区的马达则不会产生群集，这是一个十分有趣的现象。

在探究二聚体微马达形成团簇的过程中，对15个马达进行标记，在随机数为5时，对标号为4的二聚体微马达进行追踪，得到其轴向速度如图5所示。

从图5可以看出，马达4的初始时刻轴向速度较大，说明位于系统的燃料区内；随着模拟实验的推进，马达4跑到惰性粒子溶液区，此时的轴向速度减小；在16 256～16 666时间段内，其轴向速度处于高值，并于t=16 666时的轴向速度达到最大，说明此时马达4已重新进入燃料区。以马达4为中心，通过计算其他马达与该马达之间的距离发现，这些马达相互靠近并形成团簇。最终，系统中多个二聚体微马达以团簇运动的速度继续运动。图6为马达4的轴向速度分布直方图，通过计算得到的数学期望μ=0.002 5,σ=0.009 6。

图5 随机数为5时，4号二聚体微马达的轴向速度

图6 马达4的轴向速度分布直方图

2.3 马达在多个圆形燃料区环境中的自驱动与群集

通过改变燃料场的数目，继续研究马达在燃料区中的动力学性质与群集行为。在大小为Lx×Ly=100×100的准二维平面系统中进行模拟实验，设置5个大小相同的圆形燃料区，圆心分别为(20，20)，(20，80)，(80，20)，(80，80)，(50，50)，半径r=15，15个纳米二聚体微马达不均匀分布在整个系统中，15个马达在系统中的构型如图7所示。

图7 5个燃料区内，马达在不同时刻的瞬态构型

从图7可以看出，初始时刻时，燃料区和非燃料区均有马达。随着时间的推移，多个燃料区发生催化反应，t=30 000时，系统呈现相对稳定的马达团簇构型。由于单个燃料区域无法容纳全部马达，所以马达在不同的燃料区内形成团簇。团簇的形成区域是随机的，与燃料区的排布和大小有关。

对于富含燃料区和非燃料区的不均匀环境，马达在燃料丰富区运动时，呈现自驱动行为，有较高的活性。需要指出的是，由于本文研究的马达是化学场自驱动，产生的非均匀梯度场驱动马达，因此，马达之间可通过各自的梯度场发生相互耦合，这种耦合效应展示了与其他不同类型微纳米马达的不同性质。由于耦合梯度场的存在，马达之间相互吸引，在富集燃料区形成团簇。与其他类型马达不同，马达在快速运动区域反而出现群集现象。当燃料区较大时，可以聚集所有的马达；当燃料区较小时，可以在不同燃料区分别聚集密集的团簇；在缺乏燃料区时，马达四处散开随机运动，形成有趣的相分离现象。

3 结束语

本文运用混合的分子动力学和多粒子碰撞动力学MPC-MD方法研究由催化球和非催化球构成的二聚体微马达在不均匀燃料环境中的单体和多体动力学行为。研究发现，在燃料富集区域，微纳米马达展现自驱动行为并形成团簇；在非燃料区，马达呈布朗运动，具有分散-聚集的内在物理机制。为面向生物医疗功能的纳米机器的设计提供借鉴。但是，本文未对马达形成聚集的条件进行分析，后续将对此展开进一步研究，通过改变燃料场相关参数来调控二聚体微纳米马达的群集区域。