金英俊，董康杰，溫正城

(杭州電子科技大學(xué)理學(xué)院，浙江 杭州 310018)

0 引 言

自驅(qū)動(dòng)微納米馬達(dá)是一種尺寸在微納米量級(jí)的微小器件，可以將其他形式的能量轉(zhuǎn)化為動(dòng)能，實(shí)現(xiàn)自主運(yùn)動(dòng)[1]。目前，自驅(qū)動(dòng)微納米馬達(dá)的探索尚處于早期階段，譬如Paxton等[2]研究以H2O2為燃料的Au/Pt雙元金屬棒馬達(dá)在H2O2中的運(yùn)動(dòng)；Gao等[3]先后制備了以H2O2為驅(qū)動(dòng)力的Al/Ga和Mg微米馬達(dá)。馬達(dá)動(dòng)力學(xué)性質(zhì)的研究往往在較簡單的均勻環(huán)境中進(jìn)行，主要關(guān)注的是馬達(dá)的運(yùn)動(dòng)速度、驅(qū)動(dòng)機(jī)理、功能應(yīng)用等方面，燃料時(shí)空動(dòng)力學(xué)特性引起馬達(dá)運(yùn)動(dòng)性質(zhì)變化的相關(guān)研究甚少。

眾所周知，生物系統(tǒng)中的環(huán)境絕大部分是不均勻耗散體系，展示復(fù)雜的時(shí)空動(dòng)力學(xué)現(xiàn)象。比如，在不均勻的環(huán)境中，細(xì)菌往往趨向有營養(yǎng)物質(zhì)區(qū)域，遠(yuǎn)離有毒性物質(zhì)區(qū)域，展示趨利避害的重要性質(zhì)。隨著相關(guān)研究的深入，不均勻環(huán)境對(duì)馬達(dá)動(dòng)力學(xué)的影響引起研究者關(guān)注[4-5]。本文通過構(gòu)造燃料分布不均勻環(huán)境，運(yùn)用混合的分子動(dòng)力學(xué)和多粒子碰撞動(dòng)力學(xué)(Molecular Dynamics-Multi-Particle Collision Dynamics，MD-MPC)方法，探究自驅(qū)動(dòng)二聚體微納米馬達(dá)的動(dòng)力學(xué)性質(zhì)，以及馬達(dá)間相互作用對(duì)其形成團(tuán)簇的影響。

1 系統(tǒng)模型與仿真

本文主要針對(duì)自驅(qū)動(dòng)二聚體微納米馬達(dá)在不均勻環(huán)境中的動(dòng)力學(xué)性質(zhì)和群集效應(yīng)展開研究。建立尺寸為V=Lx×Ly×Lz=100×100×20的三維系統(tǒng)進(jìn)行仿真模擬。系統(tǒng)中，隨機(jī)放置大量的燃料粒子F以及不具有化學(xué)活性的惰性溶液粒子S，粒子的總數(shù)目NT=NF+NS，其中NF為燃料粒子F的數(shù)目，NS為惰性粒子S的數(shù)目。同時(shí)，系統(tǒng)中放入ND個(gè)二聚體微納米馬達(dá)，每個(gè)二聚體微納米馬達(dá)由催化球C和非催化球N組成。當(dāng)隨機(jī)運(yùn)動(dòng)的溶劑粒子碰到馬達(dá)的催化球C時(shí)，C表面發(fā)生不可逆的催化反應(yīng)，F(xiàn)+C→P+C，將燃料粒子F轉(zhuǎn)化為產(chǎn)物粒子P；而在非催化球N處，沒有催化反應(yīng)發(fā)生。二聚體微納米馬達(dá)的構(gòu)型如圖1所示，主要包括催化球C和非催化球N。

圖1 二聚體馬達(dá)構(gòu)型

將溶液中的F粒子、S粒子以及反應(yīng)生成的P粒子作為質(zhì)量均為m的質(zhì)點(diǎn)。為了限制馬達(dá)在Z方向上的運(yùn)動(dòng)，采用9-3蘭納-瓊斯(Lennard-Jones，LJ)排斥勢進(jìn)行描述，具體表達(dá)式如下：

(1)

式中，r為2個(gè)粒子之間的距離，εW表示馬達(dá)與上下墻壁之間相互作用的能量參數(shù)，σW表示馬達(dá)與上下墻壁之間相互作用的距離參數(shù)。令σW=LZ/2，其中LZ表示系統(tǒng)在Z方向上的尺寸，由于這個(gè)勢的存在，馬達(dá)只能在XY平面上做準(zhǔn)二維運(yùn)動(dòng)[6]。因?yàn)楸诿娴拇嬖?，系統(tǒng)中的溶液點(diǎn)粒子不斷通過回彈碰撞與坐標(biāo)Z=0和Z=LZ處的壁相互作用，與壁面碰撞后，小球的速度不變，并沿碰撞前的鏡面對(duì)稱方向彈回。對(duì)于含有多達(dá)1023量級(jí)的微觀粒子的宏觀系統(tǒng)，通過對(duì)整個(gè)微觀系統(tǒng)的演化來獲得宏觀物理量并不容易。為了獲得較為準(zhǔn)確的物理量，在系統(tǒng)的X方向和Y方向上，采用周期性邊界條件來正確描述二聚體馬達(dá)的動(dòng)力學(xué)性質(zhì)。

對(duì)于溶液中的粒子與二聚體微馬達(dá)之間的相互作用，采用12-6形式的LJ排斥勢進(jìn)行描述，具體表達(dá)式如下：

(2)

式中，λ為C或N，表示馬達(dá)的C球或N球，μ為F，S或P，表示溶液中的粒子，rc=21/6σ為馬達(dá)的截?cái)嗑嚯x，θ(x)為Heaviside函數(shù)。當(dāng)r≤rc時(shí)，可以通過式(2)計(jì)算得出粒子與馬達(dá)之間的作用勢；當(dāng)r>rc時(shí)，該作用勢被截?cái)郲7]。除了產(chǎn)物粒子P對(duì)N球的能量參數(shù)不同外，其他所有溶液粒子與C球和N球的相互作用能量參數(shù)相同，即εCF=εCS=εNF=εNS=εF≠εP，其中εP=εNP。

對(duì)于含有多個(gè)二聚體微馬達(dá)的系統(tǒng)，不同二聚體的單體之間也通過排斥LJ勢相互作用，其中能量參數(shù)為εD，距離參數(shù)為σD=σλJ+σλK，J=1,…,ND和K=1,…,ND表示系統(tǒng)中的二聚體微馬達(dá)。

采用MD-MPC方法對(duì)整個(gè)系統(tǒng)進(jìn)行時(shí)間演化，MPC-MD方法由分子動(dòng)力學(xué)(Molecular Dynamics，MD)和多粒子碰撞動(dòng)力學(xué)(Multi-Particle Collision Dynamics，MPC)構(gòu)成[8-9]。采用MD模擬溶液粒子運(yùn)動(dòng)的過程，通過牛頓運(yùn)動(dòng)方程演化溶液粒子與馬達(dá)的相互作用、馬達(dá)與馬達(dá)之間的相互作用，其步長為τMD。采用MPC模擬溶液粒子之間的相互作用[10]。經(jīng)過n步MD模擬后，采用MPC進(jìn)行1次溶液粒子隨機(jī)碰撞模擬，其步長τMPC=nτMD。用MPC模擬溶液粒子隨機(jī)碰撞過程中，本文把系統(tǒng)分為尺寸a=1的眾多立方體小格子，在每個(gè)單獨(dú)的方格內(nèi)，溶液粒子發(fā)生碰撞。標(biāo)記為ξ號(hào)方格內(nèi)的i號(hào)粒子的速度為：

(3)

數(shù)值模擬時(shí)，系統(tǒng)中所有的物理量均為無量綱LJ單位，分別基于能量εα，質(zhì)量m0以及距離σd進(jìn)行設(shè)置：r/σ→r，t(ε/(mσ2)1/2)→t和KBT/ε→T。具體參數(shù)如下：粒子F，S和P的質(zhì)量均為m0=1。馬達(dá)的C球、N球的直徑分別為dC=2σC=2×2=4，dN=2σN=2×4=8，為了保證馬達(dá)懸浮在溶液中，其質(zhì)量隨各自直徑變化而變化。本文模擬中，mC=335，mN=2 688，每個(gè)馬達(dá)的C球、N球之間的核間距固定為R=σC+σN+ξ=2+4+1=7。系統(tǒng)溫度設(shè)定為KBT=1/6。采用MD模擬溶液粒子運(yùn)動(dòng)過程中，時(shí)間步長ΔτMD=0.01。采用MPC模擬溶液粒子之間相互作用過程中，納米馬達(dá)的旋轉(zhuǎn)角為π/2，ΔτMPC=50×0.01=0.5，即n=50步進(jìn)行1次碰撞，同時(shí)伴隨1次能量轉(zhuǎn)換。系統(tǒng)中的粒子總數(shù)NT=2.0×106，平均密度為10，即每個(gè)小方格子內(nèi)含有10個(gè)溶液粒子。能量參數(shù)為εS=εF=5，εP=0.1，2個(gè)馬達(dá)之間的εD=5。

初始化系統(tǒng)時(shí)，確保二聚體微馬達(dá)位于系統(tǒng)中的一側(cè)，放置規(guī)則如下：系統(tǒng)含有單個(gè)馬達(dá)情形下，N球球心在X和Y方向上需滿足10≤XN≤90,10≤YN≤90，C球相對(duì)于N球的角度范圍為[0,2π]，同時(shí)滿足5≤XC≤95，5≤YC≤95。當(dāng)系統(tǒng)含有多個(gè)馬達(dá)時(shí)，還需滿足如下規(guī)則：任意2個(gè)馬達(dá)α、β的N球之間的距離dNαNβ≥2×21/6σN，C球之間的距離dCαCβ≥2×21/6σC，C球和N球之間的距離dCαNβ≥2×21/6(σC+σN)。溶液粒子需滿足條件如下：與任意馬達(dá)C球的距離dSiCα≥21/6σC，與任意馬達(dá)N球之間的距離dSiNα≥21/6σN。

2 模擬結(jié)果與分析

2.1 單個(gè)馬達(dá)的自驅(qū)動(dòng)行為

為了研究馬達(dá)在不均勻環(huán)境下的群集效應(yīng)與動(dòng)力學(xué)性質(zhì)，分別針對(duì)環(huán)境中只有惰性溶液粒子S和燃料粒子F這2種情形進(jìn)行模擬實(shí)驗(yàn)，且只存在單個(gè)二聚體微馬達(dá)。

當(dāng)環(huán)境中只含惰性溶液粒子S時(shí)，溶液粒子S與馬達(dá)的C球、N球之間沒有化學(xué)反應(yīng)，馬達(dá)周圍沒有產(chǎn)生相應(yīng)的梯度差，又粒子S和C球、N球之間的能量參數(shù)設(shè)置為εSC=εSN=5，所以，此時(shí)沒有驅(qū)動(dòng)力促使馬達(dá)做定向運(yùn)動(dòng)，馬達(dá)在系統(tǒng)中只做隨機(jī)的布朗運(yùn)動(dòng)。當(dāng)環(huán)境中只含燃料粒子F時(shí)，馬達(dá)的C球表面發(fā)生催化反應(yīng)生成產(chǎn)物粒子P的梯度場，雖然C球周圍的P粒子場是均勻的，但在N球周圍的P粒子梯度場是不均勻的，P粒子濃度沿著C-N的軸向降低，加上能量參數(shù)εF=5和εP=0.1(εF>εP)的不同，產(chǎn)物粒子P對(duì)二聚體N球產(chǎn)生的作用力小于其他溶液粒子對(duì)N球的推動(dòng)，從而在N球上產(chǎn)生N-C指向的凈力，使馬達(dá)在系統(tǒng)中產(chǎn)自驅(qū)動(dòng)。本文研究的二聚體微馬達(dá)基于自擴(kuò)散泳的機(jī)理，文獻(xiàn)[12]通過模擬驗(yàn)證了自擴(kuò)散泳機(jī)理，并且制備出這類二聚體微馬達(dá)。

圖2 馬達(dá)的均方位移ΔL2隨時(shí)間t變化的曲線

2.2 馬達(dá)在單個(gè)圓形燃料區(qū)中的自驅(qū)動(dòng)與群集

為了探究馬達(dá)在不均勻環(huán)境下的行為，對(duì)單個(gè)馬達(dá)(ND=1)在單個(gè)圓形燃料場中的行為進(jìn)行模擬實(shí)驗(yàn)。以坐標(biāo)(25,25)為圓心，r=25為半徑建立1個(gè)圓形燃料場，區(qū)域內(nèi)的溶液粒子均為具有化學(xué)活性的粒子F。初始時(shí)刻，系統(tǒng)中均勻分布大量惰性溶液粒子S，由于粒子S做無規(guī)則布朗運(yùn)動(dòng)，為了確保圓形燃料場內(nèi)均為F粒子，將運(yùn)動(dòng)到區(qū)域邊界的S粒子受到外界催化轉(zhuǎn)變?yōu)榱Ｗ覨。單個(gè)馬達(dá)按照第1節(jié)所述規(guī)則進(jìn)行放置。

圖3 二聚體微馬達(dá)的單個(gè)燃料場中運(yùn)動(dòng)及軸向速度

從圖3(b)可以看出，t在0～810和16 634～17 475這2個(gè)時(shí)間段中，馬達(dá)的軸向速度峰值遠(yuǎn)大于其他時(shí)間段的軸向速度，馬達(dá)受到額外驅(qū)動(dòng)力(FFN>FPN)的推動(dòng)，在燃料區(qū)產(chǎn)生自驅(qū)動(dòng)。

為了進(jìn)一步探究馬達(dá)在受限環(huán)境中的行為，在上述單個(gè)燃料場系統(tǒng)中放入15個(gè)二聚體微馬達(dá)，即ND=15。除了按多個(gè)馬達(dá)放置規(guī)則來放置馬達(dá)外，還要確保惰性S溶液區(qū)和燃料區(qū)均有馬達(dá)。為了觀察更多二聚體微馬達(dá)的構(gòu)型，將馬達(dá)的運(yùn)行步數(shù)增加到3.0×107。圖4為15個(gè)馬達(dá)在單個(gè)燃料場系統(tǒng)中的瞬態(tài)構(gòu)型。

圖4 單個(gè)燃料場內(nèi)，馬達(dá)在不同時(shí)刻的瞬態(tài)構(gòu)型

從圖4可以看出，初始時(shí)刻時(shí)，15個(gè)微納米馬達(dá)不均勻地分布在系統(tǒng)中，大部分馬達(dá)分布在惰性溶液粒子區(qū)，少數(shù)馬達(dá)位于燃料區(qū)內(nèi)。惰性溶液粒子區(qū)內(nèi)的馬達(dá)只做布朗運(yùn)動(dòng)，即緩慢進(jìn)行自擴(kuò)散，當(dāng)擴(kuò)散到S溶液粒子區(qū)與燃料區(qū)的交界處時(shí)，由催化反應(yīng)產(chǎn)生的N球指向C球的力推動(dòng)馬達(dá)進(jìn)入燃料區(qū)。當(dāng)然，圓形燃料區(qū)內(nèi)的馬達(dá)也會(huì)跑出燃料區(qū)進(jìn)入非燃料區(qū)。經(jīng)過一段時(shí)間演化后，多個(gè)二聚體微馬達(dá)在燃料區(qū)呈群聚狀態(tài)。當(dāng)燃料區(qū)中跑入多個(gè)馬達(dá)后，由于化學(xué)反應(yīng)，馬達(dá)附近產(chǎn)生P粒子的梯度場，考慮到εF>εP，其他馬達(dá)探測到該梯度場后，其N球向該馬達(dá)的C球靠近，從而形成團(tuán)簇。C球和N球的聚集排列展現(xiàn)出穩(wěn)定的構(gòu)型，跑出燃料區(qū)的馬達(dá)以及原本在溶液區(qū)的馬達(dá)則不會(huì)產(chǎn)生群集，這是一個(gè)十分有趣的現(xiàn)象。

在探究二聚體微馬達(dá)形成團(tuán)簇的過程中，對(duì)15個(gè)馬達(dá)進(jìn)行標(biāo)記，在隨機(jī)數(shù)為5時(shí)，對(duì)標(biāo)號(hào)為4的二聚體微馬達(dá)進(jìn)行追蹤，得到其軸向速度如圖5所示。

從圖5可以看出，馬達(dá)4的初始時(shí)刻軸向速度較大，說明位于系統(tǒng)的燃料區(qū)內(nèi)；隨著模擬實(shí)驗(yàn)的推進(jìn)，馬達(dá)4跑到惰性粒子溶液區(qū)，此時(shí)的軸向速度減小；在16 256～16 666時(shí)間段內(nèi)，其軸向速度處于高值，并于t=16 666時(shí)的軸向速度達(dá)到最大，說明此時(shí)馬達(dá)4已重新進(jìn)入燃料區(qū)。以馬達(dá)4為中心，通過計(jì)算其他馬達(dá)與該馬達(dá)之間的距離發(fā)現(xiàn)，這些馬達(dá)相互靠近并形成團(tuán)簇。最終，系統(tǒng)中多個(gè)二聚體微馬達(dá)以團(tuán)簇運(yùn)動(dòng)的速度繼續(xù)運(yùn)動(dòng)。圖6為馬達(dá)4的軸向速度分布直方圖，通過計(jì)算得到的數(shù)學(xué)期望μ=0.002 5,σ=0.009 6。

圖5 隨機(jī)數(shù)為5時(shí)，4號(hào)二聚體微馬達(dá)的軸向速度

圖6 馬達(dá)4的軸向速度分布直方圖

2.3 馬達(dá)在多個(gè)圓形燃料區(qū)環(huán)境中的自驅(qū)動(dòng)與群集

通過改變?nèi)剂蠄龅臄?shù)目，繼續(xù)研究馬達(dá)在燃料區(qū)中的動(dòng)力學(xué)性質(zhì)與群集行為。在大小為Lx×Ly=100×100的準(zhǔn)二維平面系統(tǒng)中進(jìn)行模擬實(shí)驗(yàn)，設(shè)置5個(gè)大小相同的圓形燃料區(qū)，圓心分別為(20，20)，(20，80)，(80，20)，(80，80)，(50，50)，半徑r=15，15個(gè)納米二聚體微馬達(dá)不均勻分布在整個(gè)系統(tǒng)中，15個(gè)馬達(dá)在系統(tǒng)中的構(gòu)型如圖7所示。

圖7 5個(gè)燃料區(qū)內(nèi)，馬達(dá)在不同時(shí)刻的瞬態(tài)構(gòu)型

從圖7可以看出，初始時(shí)刻時(shí)，燃料區(qū)和非燃料區(qū)均有馬達(dá)。隨著時(shí)間的推移，多個(gè)燃料區(qū)發(fā)生催化反應(yīng)，t=30 000時(shí)，系統(tǒng)呈現(xiàn)相對(duì)穩(wěn)定的馬達(dá)團(tuán)簇構(gòu)型。由于單個(gè)燃料區(qū)域無法容納全部馬達(dá)，所以馬達(dá)在不同的燃料區(qū)內(nèi)形成團(tuán)簇。團(tuán)簇的形成區(qū)域是隨機(jī)的，與燃料區(qū)的排布和大小有關(guān)。

對(duì)于富含燃料區(qū)和非燃料區(qū)的不均勻環(huán)境，馬達(dá)在燃料豐富區(qū)運(yùn)動(dòng)時(shí)，呈現(xiàn)自驅(qū)動(dòng)行為，有較高的活性。需要指出的是，由于本文研究的馬達(dá)是化學(xué)場自驅(qū)動(dòng)，產(chǎn)生的非均勻梯度場驅(qū)動(dòng)馬達(dá)，因此，馬達(dá)之間可通過各自的梯度場發(fā)生相互耦合，這種耦合效應(yīng)展示了與其他不同類型微納米馬達(dá)的不同性質(zhì)。由于耦合梯度場的存在，馬達(dá)之間相互吸引，在富集燃料區(qū)形成團(tuán)簇。與其他類型馬達(dá)不同，馬達(dá)在快速運(yùn)動(dòng)區(qū)域反而出現(xiàn)群集現(xiàn)象。當(dāng)燃料區(qū)較大時(shí)，可以聚集所有的馬達(dá)；當(dāng)燃料區(qū)較小時(shí)，可以在不同燃料區(qū)分別聚集密集的團(tuán)簇；在缺乏燃料區(qū)時(shí)，馬達(dá)四處散開隨機(jī)運(yùn)動(dòng)，形成有趣的相分離現(xiàn)象。

3 結(jié)束語

本文運(yùn)用混合的分子動(dòng)力學(xué)和多粒子碰撞動(dòng)力學(xué)MPC-MD方法研究由催化球和非催化球構(gòu)成的二聚體微馬達(dá)在不均勻燃料環(huán)境中的單體和多體動(dòng)力學(xué)行為。研究發(fā)現(xiàn)，在燃料富集區(qū)域，微納米馬達(dá)展現(xiàn)自驅(qū)動(dòng)行為并形成團(tuán)簇；在非燃料區(qū)，馬達(dá)呈布朗運(yùn)動(dòng)，具有分散-聚集的內(nèi)在物理機(jī)制。為面向生物醫(yī)療功能的納米機(jī)器的設(shè)計(jì)提供借鑒。但是，本文未對(duì)馬達(dá)形成聚集的條件進(jìn)行分析，后續(xù)將對(duì)此展開進(jìn)一步研究，通過改變?nèi)剂蠄鱿嚓P(guān)參數(shù)來調(diào)控二聚體微納米馬達(dá)的群集區(qū)域。