夏益祺 諶莊琳 郭永坤

(蘇州大學(xué)軟凝聚態(tài)物理及交叉研究中心,蘇州大學(xué)物理科學(xué)與技術(shù)學(xué)院,蘇州 215006)

1 引 言

從微觀尺度的細(xì)胞骨架到宏觀尺度的魚(yú)群、鳥(niǎo)群[1],活性物質(zhì)(active matter)普遍存在于自然界中.人工也能設(shè)計(jì)和制造活性物質(zhì)體系,如磁場(chǎng)[2]、聲場(chǎng)[3]或電場(chǎng)[4]驅(qū)動(dòng)的活性體系.這些系統(tǒng)的共同特征是組成單元可以通過(guò)消耗能量進(jìn)行運(yùn)動(dòng).此類(lèi)系統(tǒng)展現(xiàn)出許多有趣的非平衡物理現(xiàn)象,如體系密度巨漲落(giant fluctuations)[5]、自發(fā)相分離(spontaneous phase separation)[6,7]、斑圖形成(pattern formation)[8]和自組裝行為(selfassembly behavior)[9].

除了理解此類(lèi)體系豐富現(xiàn)象背后的物理機(jī)理,科學(xué)家們開(kāi)始思索如何將活性物質(zhì)體系的能量提取出來(lái)為人類(lèi)所用,比如細(xì)菌自發(fā)運(yùn)動(dòng)的能量[10-13].細(xì)菌的運(yùn)動(dòng)是通過(guò)消耗化學(xué)能來(lái)實(shí)現(xiàn),其自驅(qū)運(yùn)動(dòng)本身已經(jīng)打破系統(tǒng)的時(shí)間對(duì)稱(chēng)性.為了實(shí)現(xiàn)定向旋轉(zhuǎn),還需要打破空間對(duì)稱(chēng)性,常見(jiàn)的方法是使用形狀不對(duì)稱(chēng)的棘輪.比如,Angelani等[14]的模擬研究表明,非對(duì)稱(chēng)棘輪浸入棒狀粒子能夠產(chǎn)生定向轉(zhuǎn)動(dòng).Sokolov等[15]在細(xì)菌作為活性物質(zhì)的實(shí)驗(yàn)中也觀察到了類(lèi)似的現(xiàn)象.近年來(lái),將柔性鏈狀結(jié)構(gòu)的軟邊界[16-23]置于活性粒子浴中的動(dòng)力學(xué)行為越來(lái)越受研究者的關(guān)注.由于活性粒子的碰撞,軟邊界會(huì)發(fā)生非熱波動(dòng)進(jìn)而出現(xiàn)反常的現(xiàn)象.最近,我們課題組研究發(fā)現(xiàn),柔性[21]和半剛性鏈[22]嫁接的圓盤(pán)在活性粒子浴中會(huì)因?yàn)殒湗?gòu)象的對(duì)稱(chēng)性被打破,自發(fā)產(chǎn)生結(jié)構(gòu)的不對(duì)稱(chēng)性,從而產(chǎn)生長(zhǎng)時(shí)間持續(xù)的定向運(yùn)動(dòng).

我們之前的研究工作表明[21],活性粒子傾向于聚集在鏈上高曲率區(qū)域,而柔性鏈的形變對(duì)粒子的積聚存在正反饋效果[23],從而進(jìn)一步加強(qiáng)柔性鏈形變,使粒子大量聚集.上述工作中我們通過(guò)多根鏈嫁接到納米盤(pán)來(lái)實(shí)現(xiàn)自發(fā)轉(zhuǎn)動(dòng),本文中我們?cè)O(shè)計(jì)更簡(jiǎn)單的具有柔性邊界的對(duì)稱(chēng)性棘輪,研究其在非平衡粒子浴內(nèi)的行為,進(jìn)一步探討柔性邊界是否能讓棘輪轉(zhuǎn)動(dòng).我們發(fā)現(xiàn)粒子的聚集導(dǎo)致柔性邊界的結(jié)構(gòu)發(fā)生改變,力分布的不對(duì)稱(chēng)性誘導(dǎo)棘輪產(chǎn)生扭矩,進(jìn)而實(shí)現(xiàn)棘輪的定向轉(zhuǎn)動(dòng).活性粒子的自驅(qū)力大小、旋轉(zhuǎn)擴(kuò)散能力和體系粒子密度都對(duì)棘輪定向轉(zhuǎn)動(dòng)產(chǎn)生影響: 棘輪的旋轉(zhuǎn)速率隨活性力和活性粒子數(shù)密度的增大而增加; 隨活性粒子的旋轉(zhuǎn)擴(kuò)散系數(shù)的增大而減小.

2 模型和研究方法

類(lèi)似于之前的工作[19-23],我們用二維模型來(lái)進(jìn)行研究(模擬活性物質(zhì)在界面或薄膜中的情況).體系的初始結(jié)構(gòu)為Np個(gè)自驅(qū)動(dòng)活性粒子圍繞在柔性邊界周?chē)?其中活性粒子(active particles)的力驅(qū)動(dòng)方向?yàn)?(cosθ,sinθ),該方向可以隨時(shí)間改變; 柔性邊界的初始狀態(tài)為一個(gè)圓形由Nb=1256個(gè)珠子(passive beads)用彈簧勢(shì)連接的柔性鏈組成.整個(gè)體系中所有粒子之間采用相互排斥的WCA勢(shì),其勢(shì)函數(shù)表示為

其中,彈性系數(shù)取k=5000ε/σ2,為了防止活性粒子穿過(guò)柔性鏈,平衡鍵長(zhǎng)取r0=0.25σ,由鍵連接的近鄰珠子之間的對(duì)相互作用被忽略.棘輪由柔性鏈結(jié)構(gòu)和剛性支架組成.柔性鏈的初始構(gòu)型為圓形,半徑R0≈ 50σ.鏈到圓心由剛性支架連接,結(jié)構(gòu)類(lèi)似于自行車(chē)輪.四支架體系,初始結(jié)構(gòu)示意圖如圖1.

體系中的活性粒子遵循郎之萬(wàn)動(dòng)力學(xué),其運(yùn)動(dòng)方程可以表示為:

其中,(3)和(4)式分別表示活性粒子的平移和旋轉(zhuǎn)運(yùn)動(dòng); 柔性鏈上的珠子只遵循運(yùn)動(dòng)方程(3),且F=0; 平動(dòng)摩擦系數(shù)為γ=kBT/D0,D0為平移擴(kuò)散系數(shù);Dr為旋轉(zhuǎn)擴(kuò)散系數(shù);=(cosθ,sinθ)是一個(gè)單位向量指向沿著活性力的推進(jìn)軸,由旋轉(zhuǎn)擴(kuò)散控制;ηi(t)和ξi(t)為白噪聲項(xiàng),滿(mǎn)足漲落-耗散定理:

在我們的分子動(dòng)力學(xué)模擬過(guò)程中,周期性盒子大小為 200σ×200σ.我們采用約化單位,令m=1,σ=1,kBT=1; 時(shí)間單位為此外,本文中我們令ε=10,γ=10,此時(shí)D0=0.1.我們采用修改后的LAMMPS軟件[24,25]進(jìn)行模擬,研究不同的活性力(F)、粒子數(shù)密度(φ)、旋轉(zhuǎn)擴(kuò)散系數(shù)(Dr)及支架數(shù)(Ns)對(duì)棘輪運(yùn)動(dòng)行為的影響.運(yùn)動(dòng)軌跡采用VMD軟件進(jìn)行可視化處理和分析[26].

圖1 棘輪在活性粒子浴中(綠色部分)的初始示意結(jié)構(gòu)虛線部分表示剛性支架,紅色表示柔性邊界,左上部分綠色粒子中的黑色箭頭表示活性粒子的推進(jìn)力方向,右上部分為一段邊界上的粒子排布,粒子間間距為0.25σFig.1.The initial structure of flexible boundary(red)in the active particle bath(green).The dashed lines represent four rigid trestles.The black arrow in the green particle indicates the active force direction of the active particle in the upper left.The upper right part is the particle arrangement on a boundary.The spacing between the particles is 0.25σ.

3 模擬結(jié)果與討論

3.1 壓力不對(duì)稱(chēng)分布誘導(dǎo)轉(zhuǎn)動(dòng)

首先考察了支架數(shù)為Ns=4,活性粒子的數(shù)密度為φ=0.025、旋轉(zhuǎn)擴(kuò)散系數(shù)取Dr=0.0001、驅(qū)動(dòng)力F=40時(shí)棘輪的轉(zhuǎn)動(dòng)行為(見(jiàn)附屬材料中的視頻F=40.avi).如圖2(a)所示,當(dāng)棘輪處于活性粒子浴中會(huì)發(fā)生定向旋轉(zhuǎn)的現(xiàn)象.

經(jīng)過(guò)觀察棘輪及活性粒子的運(yùn)動(dòng)軌跡,我們發(fā)現(xiàn)開(kāi)始時(shí)刻棘輪不轉(zhuǎn)動(dòng),活性粒子在驅(qū)動(dòng)力作用下,快速聚集到柔性邊界周?chē)?并誘導(dǎo)其發(fā)生形變,如圖2(a)的左上圖.當(dāng)粒子聚集到一定程度,集聚在棘輪上曲率較大處的活性粒子出現(xiàn)崩塌的情況:在棘輪上接近支架頂端處的大部分活性粒子會(huì)被棘輪拋出,也就是甩尾現(xiàn)象(圖2(a)中紅色圓圈處).這種運(yùn)動(dòng)導(dǎo)致粒子在邊界上的分布不均勻,進(jìn)而對(duì)棘輪施加一定的扭矩,導(dǎo)致棘輪的旋轉(zhuǎn).圖2(a)的右下圖為旋轉(zhuǎn)角度θ隨著時(shí)間的變化曲線,里面的插圖為對(duì)應(yīng)的角速度ω隨時(shí)間的變化曲線.旋轉(zhuǎn)角速度的計(jì)算依據(jù)以下公式:

其中,旋轉(zhuǎn)角速度ω為正值時(shí)表示棘輪沿逆時(shí)針?lè)较蛐D(zhuǎn),負(fù)值時(shí)表示棘輪沿順時(shí)針?lè)较蛐D(zhuǎn).

圖2 (a)四支架棘輪隨時(shí)間轉(zhuǎn)動(dòng)的示意圖,圖示中剛性支架被忽略,左上圖黑點(diǎn)代表棘輪質(zhì)心,其中黑色箭頭表示從質(zhì)心到一支點(diǎn)的方向以此表明棘輪發(fā)生轉(zhuǎn)動(dòng),右下圖是該情況下的旋轉(zhuǎn)角度與角速度隨著時(shí)間的變化曲線圖;(b)上圖為棘輪柔性邊界的曲率分布圖,下圖為壓力分布圖,其中棘輪支架數(shù)為Ns=4,活性粒子的數(shù)密度為φ=0.025,旋轉(zhuǎn)擴(kuò)散系數(shù)取Dr=0.0001,驅(qū)動(dòng)力F=40Fig.2.(a)Schematic diagram of the ratchet rotation with time for Ns=4,φ=0.025,Dr=0.0001,F=40.The black dot on the top left figure represents the center of mass of the ratchet; the black arrows denote the orientation from the center of mass of ratchet to the fixed point of boundary to indicate the rotation of the ratchet.The right-down figure shows the time evolution of angle and angular velocity of ratchet.(b)The curvature distribution(up)of boundary and the pressure distribution(down)around it.

為了更好地探究棘輪旋轉(zhuǎn)的機(jī)制,計(jì)算了構(gòu)成棘輪的柔性鏈的曲率和所受壓力的分布情況,如圖2(b)所示.該圖的上圖為棘輪的柔性鏈的曲率分布圖,其中散點(diǎn)顏色由深藍(lán)到紅色代表曲率由小到大變化; 下圖為組成棘輪的柔性鏈所受壓力的分布散點(diǎn)圖,其中深藍(lán)色部分代表組成棘輪的柔性鏈所受壓力較小,由藍(lán)色到紅色部分,表示組成棘輪的柔性鏈所受壓力逐漸增大.經(jīng)過(guò)觀察運(yùn)動(dòng)軌跡,發(fā)現(xiàn)棘輪沿順時(shí)針?lè)较蛐D(zhuǎn),與曲率、壓力分布圖所表示的效果一致.壓力分布在高曲率位置,且具有不對(duì)稱(chēng)性,進(jìn)而產(chǎn)生扭轉(zhuǎn)力,驅(qū)使其轉(zhuǎn)動(dòng).

3.2 自驅(qū)力效應(yīng)

接下來(lái)研究轉(zhuǎn)動(dòng)速度和驅(qū)動(dòng)力之間的關(guān)系.選取了支架數(shù)為Ns=4的棘輪,體系中活性粒子的數(shù)密度為φ=0.025,其中φ=σ2Np/4R02.如圖3(a)所示,隨著活性力的增加,棘輪的平均旋轉(zhuǎn)角速度單調(diào)增加.圖3(b)是體系進(jìn)入穩(wěn)態(tài)后在不同活性力的作用下,棘輪的瞬時(shí)旋轉(zhuǎn)角速度.從圖3(b)中可以看出,當(dāng)活性力較小時(shí)(F=20)棘輪的旋轉(zhuǎn)角速度較小且瞬時(shí)角速度有正有負(fù),棘輪的旋轉(zhuǎn)方向不固定.隨著活性力的進(jìn)一步增加(從F=20增大到F=80),棘輪上聚集的活性粒子增多且扭矩也會(huì)變得越來(lái)越大,從而實(shí)現(xiàn)定向轉(zhuǎn)動(dòng).通過(guò)圖3也可以發(fā)現(xiàn),棘輪的瞬時(shí)角速度漲落隨驅(qū)動(dòng)力的增大而在一定程度上增大,這是因?yàn)轵?qū)動(dòng)力增大,粒子運(yùn)動(dòng)速度較快,粒子碰撞棘輪時(shí),對(duì)棘輪的擾動(dòng)性增強(qiáng).

3.3 密度效應(yīng)

在我們之前的工作中[20]發(fā)現(xiàn)囊泡內(nèi)的活性粒子的非均勻分布會(huì)對(duì)囊泡的形狀產(chǎn)生一定的影響,且囊泡的形狀與活性粒子的密度成正相關(guān)關(guān)系[27].那么對(duì)于本文所模擬的體系,活性粒子的密度也勢(shì)必會(huì)對(duì)柔性邊界有一定影響,且可能會(huì)對(duì)棘輪的旋轉(zhuǎn)有促進(jìn)或者抑制作用.為了更好地驗(yàn)證我們的設(shè)想,在模擬過(guò)程中我們保持活性粒子的旋轉(zhuǎn)擴(kuò)散系數(shù)為Dr=0.0001,活性粒子施加的驅(qū)動(dòng)力選取F=80,研究不同粒子數(shù)密度φ下棘輪的旋轉(zhuǎn)角速度的變化.從圖4中可以看出,隨著活性粒子數(shù)密度的增加,棘輪的角速度也增加,由此可以理解為棘輪的轉(zhuǎn)動(dòng)速率隨著活性粒子數(shù)密度的提高而提高.從圖中可以直觀地看出,當(dāng)活性粒子的數(shù)密度較小時(shí),活性粒子浴中活性粒子的集聚效果沒(méi)有粒子數(shù)密度大時(shí)明顯(圖4插圖).大量粒子的聚集有利于棘輪扭矩的增加,進(jìn)而增加角速度.因邊界周?chē)枷輩^(qū)域聚集的數(shù)目有限,從圖中可以看出,隨著粒子數(shù)密度的增加,角速度增幅降低; 大量的增加粒子數(shù)密度,可能出現(xiàn)飽和現(xiàn)象.

圖3 支架數(shù)為4、不同活性力的作用下,棘輪的旋轉(zhuǎn)角速度(a)平均角速度,(b)瞬時(shí)角速度,誤差棒表示平均角速度的標(biāo)準(zhǔn)方差Fig.3.(a)Average angular velocity and(b)instantaneous angular velocity for the four rigid trestles with various active forces.The error bars represent the standard the deviations of angular velocity of rigid trestles.

圖4 棘輪的平均角速度與粒子數(shù)密度之間的關(guān)系,內(nèi)部?jī)蓮埐鍒D分別對(duì)應(yīng)活性粒子密度φ=0.025和φ=0.1,誤差棒表示平均角速度的標(biāo)準(zhǔn)方差Fig.4.Average angular velocity vs the number density of active particles.The insets are typical snapshots for φ=0.025 and φ=0.1.The error bars represent the standard the deviations of angular velocity.

3.4 旋轉(zhuǎn)擴(kuò)散系數(shù)效應(yīng)

研究過(guò)程中,我們把旋轉(zhuǎn)擴(kuò)散系數(shù)Dr作為可調(diào)參數(shù)來(lái)探究其對(duì)棘輪轉(zhuǎn)動(dòng)速率的影響[23].在模擬時(shí),為了便于觀察,給活性粒子施加的驅(qū)動(dòng)力取較大的一組F=80.得到不同Dr下棘輪的旋轉(zhuǎn)角速度與時(shí)間的關(guān)系圖,如圖5(a)所示.

從旋轉(zhuǎn)角速度的變化曲線可以看出,在旋轉(zhuǎn)擴(kuò)散系數(shù)Dr較小時(shí),棘輪的轉(zhuǎn)動(dòng)速率較高,且隨著旋轉(zhuǎn)擴(kuò)散系數(shù)Dr的增大,棘輪的轉(zhuǎn)動(dòng)速率會(huì)降低,當(dāng)增大到一定程度時(shí),棘輪的轉(zhuǎn)動(dòng)速率基本為零且棘輪出現(xiàn)不定向旋轉(zhuǎn)的情況.圖5(b)為棘輪的旋轉(zhuǎn)角速度平均值與旋轉(zhuǎn)擴(kuò)散系數(shù)Dr的關(guān)系圖,由此可以得出,棘輪的轉(zhuǎn)動(dòng)速率隨著Dr的變大而降低,并且當(dāng)其達(dá)到一定值時(shí),定向旋轉(zhuǎn)消失.因粒子脫離柔性邊界的能力隨著Dr的增加而增加[28],換句話說(shuō),Dr的增加導(dǎo)致邊界聚集的粒子數(shù)目減少(圖5(b)插圖),從而導(dǎo)致棘輪扭矩降低,定向轉(zhuǎn)動(dòng)速度下降.當(dāng)Dr增加到一定程度,活性粒子在邊界滯留時(shí)間大大降低,粒子對(duì)邊界的碰撞和熱碰撞類(lèi)似,定向轉(zhuǎn)動(dòng)消失.

3.5 支架數(shù)效應(yīng)

最后考慮支架數(shù)對(duì)棘輪旋轉(zhuǎn)速率的影響.因?yàn)榛钚粤^小時(shí)棘輪的旋轉(zhuǎn)方向不固定,為了避免棘輪的旋轉(zhuǎn)方向?qū)ρ芯拷Y(jié)果的影響,選取活性力為F=80的情況.圖6為棘輪旋轉(zhuǎn)角速度的平均值與棘輪支架數(shù)的影響,從圖像中可以直觀地看出,在活性力一定時(shí),棘輪的旋轉(zhuǎn)角速度隨著支架數(shù)的增大而減小.當(dāng)支架數(shù)增加時(shí),凹陷區(qū)增加,大量粒子或聚集到這些區(qū)域,減弱了粒子分布的不均勻性(圖6插圖),進(jìn)而使扭矩降低,從而引起轉(zhuǎn)動(dòng)角速度下降.

圖5 (a)不同旋轉(zhuǎn)擴(kuò)散系數(shù)Dr下棘輪的旋轉(zhuǎn)角速度與時(shí)間的關(guān)系;(b)平均旋轉(zhuǎn)角速度與Dr的關(guān)系,內(nèi)部?jī)刹鍒D分別對(duì)應(yīng)旋轉(zhuǎn)擴(kuò)散系數(shù)Dr=0.0001和Dr=0.01,誤差棒表示平均角速度的標(biāo)準(zhǔn)方差Fig.5.(a)The angular velocity as a function of time for various rotational diffusion coefficients,Drs.(b)The average angular velocity as a function of Dr.The insets are typical snapshots for Dr=0.0001和and Dr=0.1.The error bars represent the standard the deviations of angular velocity.

3.6 討 論

本文初始的柔性邊界是圓形的,在這種對(duì)稱(chēng)結(jié)構(gòu)下,活性粒子開(kāi)始時(shí)刻的隨機(jī)堆積也是比較對(duì)稱(chēng)的,粒子分布的不對(duì)稱(chēng)性來(lái)源于粒子分布的漲落被柔性邊界加強(qiáng),因?yàn)榱Ｗ觾A向于堆積在高曲率地方[29].這種事件也是隨機(jī)發(fā)生的,導(dǎo)致棘輪發(fā)生逆時(shí)針和順時(shí)針?lè)较蚨ㄏ蜣D(zhuǎn)動(dòng)的概率在統(tǒng)計(jì)下應(yīng)該是相同的.此外,我們采用的布朗動(dòng)力學(xué)方法,忽略了流體力學(xué)的相互作用,流體力學(xué)效應(yīng)需要進(jìn)一步的研究.棘輪的轉(zhuǎn)動(dòng)依賴(lài)于壓力的不對(duì)稱(chēng)分布,壓力的不對(duì)稱(chēng)分布不但和粒子驅(qū)動(dòng)力有關(guān),而且和柔性邊界的設(shè)計(jì)有關(guān).本文中我們?cè)O(shè)計(jì)的柔性邊界大于棘輪固定點(diǎn)的直線鏈接長(zhǎng)度,利于邊界高曲率區(qū)域的形成和活性粒子聚集,較短的邊界也可能導(dǎo)致定向旋轉(zhuǎn)消失.此外,我們也發(fā)現(xiàn)在小自驅(qū)動(dòng)力下,棘輪有頻繁反轉(zhuǎn)現(xiàn)象,定向轉(zhuǎn)動(dòng)時(shí)間較短,我們一般不定義其為定向轉(zhuǎn)動(dòng).我們定義的定向轉(zhuǎn)動(dòng)是在模擬時(shí)間尺度下沒(méi)有出現(xiàn)頻繁反轉(zhuǎn)現(xiàn)象的情況.

圖6 棘輪的平均角速度與支架數(shù)之間的關(guān)系,內(nèi)部?jī)刹鍒D分別對(duì)應(yīng)支架數(shù)為4和6,誤差棒表示平均角速度的標(biāo)準(zhǔn)方差Fig.6.The average angular velocity of ratchet vs the number of rigid trestle.The insets are typical snapshots for the number of ratchet wheel.The error bars represent the standard the deviations of angular velocity.

4 結(jié) 論

本文中用計(jì)算機(jī)模擬方法設(shè)計(jì)了一種柔性棘輪,并研究了其在活性粒子浴內(nèi)的動(dòng)力學(xué)行為.發(fā)現(xiàn)棘輪結(jié)構(gòu)在活性粒子浴中會(huì)自發(fā)產(chǎn)生結(jié)構(gòu)不對(duì)稱(chēng)性實(shí)現(xiàn)定向旋轉(zhuǎn)運(yùn)動(dòng).研究了棘輪的旋轉(zhuǎn)角速度與棘輪支架數(shù)、自驅(qū)力、數(shù)密度以及旋轉(zhuǎn)擴(kuò)散系數(shù)的關(guān)系,發(fā)現(xiàn): 1)棘輪的轉(zhuǎn)動(dòng)速度隨著活性力的增強(qiáng)而變大; 2)隨著活性粒子密度的提高,棘輪旋轉(zhuǎn)速度提高; 3)當(dāng)活性粒子的旋轉(zhuǎn)擴(kuò)散系數(shù)提高時(shí),旋轉(zhuǎn)速度下降,并且當(dāng)其超過(guò)一定值時(shí),定向旋轉(zhuǎn)消失;4)棘輪的轉(zhuǎn)動(dòng)速率隨著棘輪支架數(shù)的增加而降低.