有取向顆粒的布朗運(yùn)動(dòng)模擬研究

鄒全勝

1 概述

眾所周知，布朗運(yùn)動(dòng)是小顆粒在分子的無(wú)規(guī)則熱運(yùn)動(dòng)下受到來(lái)自各個(gè)方向的撞擊而引發(fā)的無(wú)規(guī)則動(dòng)力學(xué)運(yùn)動(dòng)，是分子動(dòng)理論的著名證據(jù)之一。分子在永不停息的、無(wú)規(guī)則的熱運(yùn)動(dòng)中不斷撞擊宏觀物體的表面，給予其表面上的隨機(jī)分布力。對(duì)于一般的宏觀物體，由于其表面積和質(zhì)量都較大，故在每一時(shí)間段內(nèi)受到大量分子的撞擊，且一定大小沖量所引起的速度改變更小，由于統(tǒng)計(jì)理論中隨機(jī)變量所具有的抵償性，故隨機(jī)合力的增長(zhǎng)冪次較之宏觀物體的質(zhì)量的增長(zhǎng)冪次小，很難在宏觀上觀察到其無(wú)規(guī)則運(yùn)動(dòng)。但對(duì)于小到介觀尺度的顆粒，例如，其尺度小到微米甚至納米級(jí)別，則統(tǒng)計(jì)漲落的重要性開(kāi)始顯現(xiàn)，在每一時(shí)刻，隨機(jī)分布力都無(wú)法充分抵償而形成某一隨機(jī)方向的合力，這就是布朗運(yùn)動(dòng)的成因。[1]

隨著顯微技術(shù)的發(fā)展，人類(lèi)對(duì)微小對(duì)象的認(rèn)知越來(lái)越直觀和精細(xì)。1827年英國(guó)植物學(xué)家R.布朗通過(guò)用顯微鏡觀察水溶液中的花粉顆粒，布朗運(yùn)動(dòng)由此被人們所熟知。充分體現(xiàn)了生產(chǎn)力和技術(shù)的發(fā)展對(duì)科學(xué)進(jìn)步的推動(dòng)作用。在今天，電子計(jì)算機(jī)的廣泛使用使得計(jì)算和模擬成為研究物理問(wèn)題的強(qiáng)大而高效的手段。[1]布朗運(yùn)動(dòng)可以抽象為數(shù)學(xué)中的隨機(jī)行走模型，用計(jì)算機(jī)來(lái)模擬顆粒的布朗運(yùn)動(dòng)，可以非常方便地研究諸如布朗運(yùn)動(dòng)的位移期望與時(shí)間的關(guān)系等諸多問(wèn)題。這一類(lèi)工作已經(jīng)有很多文獻(xiàn)介紹過(guò)。[2-3]但是，如果進(jìn)行布朗運(yùn)動(dòng)的小顆粒不再視為各向同性的球體，例如，納米棒材料，那么會(huì)得出什么新的結(jié)論？

納米棒材料，例如金納米棒，具有豐富的物理化學(xué)性質(zhì)，近年來(lái)在材料科學(xué)中得到了越來(lái)越廣泛的關(guān)注。在生物醫(yī)學(xué)、催化化學(xué)、傳感器和光學(xué)器件方面都有重要應(yīng)用[4]。是否能夠找到一種基本的、模擬其運(yùn)動(dòng)規(guī)律的計(jì)算方法？這是用計(jì)算機(jī)模擬來(lái)代替實(shí)驗(yàn)操作而進(jìn)行納米棒各種應(yīng)用的探索的基礎(chǔ)性工作之一。

納米棒具有取向性，其布朗運(yùn)動(dòng)的自由度更高，角量自由度在布朗運(yùn)動(dòng)中會(huì)如何變化？由角量自由度描述的旋轉(zhuǎn)是否與平動(dòng)存在耦合？本文考慮在流體介質(zhì)中納米棒的二維布朗運(yùn)動(dòng)，以動(dòng)力學(xué)方程為依據(jù)建立可在計(jì)算機(jī)上運(yùn)行的隨機(jī)行走模型，將模擬結(jié)果繪制成圖，并得出反映有取向顆粒布朗運(yùn)動(dòng)的關(guān)鍵性質(zhì)的若干重要結(jié)論。

2 動(dòng)力學(xué)模型

2.1 朗之萬(wàn)方程

描述各向同性顆粒(如小球)的布朗運(yùn)動(dòng)的動(dòng)力學(xué)方程為郎之萬(wàn)方程[5]：

其中，F(xiàn)代表隨機(jī)出現(xiàn)的布朗力。用隨機(jī)行走模擬此布朗運(yùn)動(dòng)時(shí)，普遍采用的算法是設(shè)定一個(gè)步長(zhǎng)，每一步都以等概率向任何一個(gè)方向邁出此步長(zhǎng)[2-3]。

對(duì)于納米棒，一方面它的隨機(jī)行走并不能輕易斷言為各個(gè)方向等概率等距離的——垂直與取向和沿著取向的情況，從直觀上看并不會(huì)一樣。另一方面，在隨機(jī)行走的過(guò)程中，所需要關(guān)注的并不僅限于質(zhì)心的坐標(biāo)，還有取向的坐標(biāo)，或者說(shuō)取向與坐標(biāo)軸的夾角——也就是納米棒的自轉(zhuǎn)。

2.2 隨機(jī)力的簡(jiǎn)化模型

接下來(lái)將建立模擬二維平面上納米棒的布朗運(yùn)動(dòng)的模型，簡(jiǎn)明起見(jiàn)，納米棒簡(jiǎn)化為長(zhǎng)a寬b的矩形。

對(duì)于被大量分子包圍的矩形，其邊界上每一處都會(huì)承受大量分子的撞擊而有垂直于邊界的壓強(qiáng)。因此矩形邊界上每一點(diǎn)都受到一個(gè)法向隨機(jī)力，這些隨機(jī)力可以理解為成對(duì)出現(xiàn)的：矩形上每一點(diǎn)所受到的隨機(jī)力都對(duì)應(yīng)于對(duì)邊上同一點(diǎn)所受到的分布一致、反向、共線的隨機(jī)力。這一對(duì)力的合力的期望為零，呈正態(tài)分布，分別在長(zhǎng)方體的兩條軸線上鋪開(kāi)。

各向同性顆粒的隨機(jī)行走模擬中，每一幀何方向邁出一步，可以理解為每一幀隨機(jī)受到某一點(diǎn)的一個(gè)力，與此類(lèi)似，對(duì)隨機(jī)力作以下設(shè)定：在每一幀，矩形上的某一點(diǎn)受到單個(gè)的力，為了更加簡(jiǎn)化，可以直接設(shè)定每次力的大小都一樣，由此代替呈正態(tài)分布的隨機(jī)力，這在統(tǒng)計(jì)意義上并不會(huì)引起差異。

2.3 平均力的簡(jiǎn)化模型

除了漲落引起的隨機(jī)力外，流體的阻力等平均力是不可忽略的。對(duì)流場(chǎng)中物體受力的計(jì)算是一個(gè)十分復(fù)雜的問(wèn)題，流體力學(xué)中的基本方程是納維-斯托克斯方程：

▽p+pF+μΔv

要用此方程來(lái)解析求解平面平行運(yùn)動(dòng)的矩形所受的阻力和阻力矩，將會(huì)十分困難。目前主要的計(jì)算方法除了直接數(shù)值求解動(dòng)力學(xué)方程外還有格子玻爾茲曼方法等[6]。本文采用一種易于計(jì)算的簡(jiǎn)單假設(shè)，來(lái)考慮流體對(duì)納米棒的整體作用：當(dāng)固體邊界各處存在一定法向速度分布時(shí)，統(tǒng)計(jì)上看前進(jìn)方向前端的邊界將會(huì)受到更多流體分子的碰撞。假設(shè)流體分子的速度分布服從氣體的麥克斯韋分布律，與邊界法向發(fā)生一維彈性碰撞而改變動(dòng)量，只需考慮這一方向上速度的分布律：

f(v)=Ae-λv2

在單位時(shí)間Δt內(nèi)，對(duì)長(zhǎng)ΔL的一小段固體壁，如果以向內(nèi)的法向速度u前進(jìn)，則對(duì)于法向速度為v的粒子群，只有面積為ΔS=ΔL(v-u)Δt的流體分子會(huì)與其發(fā)生彈性碰撞，由于流體分子質(zhì)量遠(yuǎn)小于納米棒，故單個(gè)分子的動(dòng)量改變?yōu)棣=2m(v-u)。如果流體分子的面密度為σ(在二維情況下)，則速度分布在v附近dv區(qū)間的粒子所造成的總的動(dòng)量改變dΔp=2σΔS(v-u)f(v)dv。

將不同速度對(duì)應(yīng)的面積代入，則所有速度大于u的分子的總動(dòng)量變化為：

這個(gè)積分本質(zhì)上是不完全伽馬函數(shù)，但是由于布朗運(yùn)動(dòng)的速度在數(shù)量級(jí)上遠(yuǎn)遠(yuǎn)小于分子熱運(yùn)動(dòng)速度，故可以用泰勒展開(kāi)進(jìn)行近似：

其中，第一項(xiàng)是邊界靜止時(shí)的平均壓強(qiáng)，它顯然是一個(gè)定值。對(duì)任意簡(jiǎn)單閉曲線，它本身和它的一次矩的積分都為0，這意味著它不提供合力和合力矩。第二項(xiàng)則是邊界有法向速度時(shí)對(duì)壓強(qiáng)的線性影響，運(yùn)動(dòng)方向前端的邊界所受壓強(qiáng)更大，后端的邊界所受壓強(qiáng)更小，呈現(xiàn)阻力的特征，這與預(yù)期相符。為了便于計(jì)算，本文中忽略切向的粘滯力，假定流體對(duì)顆粒的整體阻力為正比于邊界法向速度的壓強(qiáng)。

3 步長(zhǎng)的設(shè)置

3.1 平動(dòng)步長(zhǎng)

可見(jiàn)，對(duì)于不同尺寸的納米棒，步長(zhǎng)也是不同的?？梢钥紤]以1*1的納米棒(正方形納米顆粒)為基準(zhǔn)，設(shè)定一個(gè)步長(zhǎng)S。在程序中我們?cè)O(shè)定S=0.05。

通過(guò)計(jì)算機(jī)程序，可以實(shí)現(xiàn)在矩形的兩條軸線上以等概率隨機(jī)取點(diǎn)。

3.2 轉(zhuǎn)動(dòng)步長(zhǎng)

平動(dòng)在之前已經(jīng)論述過(guò)，有

F=f=Lαv

即平動(dòng)步長(zhǎng)的平方和邊長(zhǎng)成反比，S在程序中設(shè)定為0.05。例如，對(duì)于16*5的納米棒，如果力作用在長(zhǎng)邊上，則此幀的移動(dòng)距離為0.0125。

現(xiàn)在計(jì)算轉(zhuǎn)動(dòng)角速度為w時(shí)粘滯阻力所產(chǎn)生的力矩。

當(dāng)納米棒繞中心轉(zhuǎn)動(dòng)時(shí)，邊界承擔(dān)阻力，對(duì)于其中一條承擔(dān)阻力的邊界：

v(x)=ωx

df=αv(x)dx=αωxdx

dM=xdf=αωx2dx

則四條變邊界段所受合力矩：

為了達(dá)到力矩平衡，有

以上分別導(dǎo)出了對(duì)于任何給定的隨機(jī)力，這一陣幀行走時(shí)應(yīng)走的步長(zhǎng)和角度。至此，用逐幀隨機(jī)行走模擬納米棒的隨機(jī)行走的算法已經(jīng)完成。

4 模擬結(jié)果

在計(jì)算機(jī)上編寫(xiě)相應(yīng)程序，就可以用蒙特卡洛方法模擬上述的過(guò)程。以下展示一些典型的模擬結(jié)果。

4.1 納米方塊( 1*1)的一次隨機(jī)行走

這是正方形的納米棒，并且作為校對(duì)步長(zhǎng)的基準(zhǔn)。以下是2萬(wàn)步質(zhì)心坐標(biāo)路徑圖和角度變化趨勢(shì)(圖1，圖2)。

圖1

圖2

由于步長(zhǎng)設(shè)置得很短，因此2萬(wàn)步并沒(méi)有走太遠(yuǎn)，覆蓋面積和1*1的納米方塊尺寸比起來(lái)只有大約幾十倍。角度的變化范圍大約在2-3個(gè)圓周內(nèi)。

將取向作為z軸，畫(huà)三維圖(圖3)。

圖3

4.2 小納米棒(4*1)的一次隨機(jī)行走

作出4*1小納米棒質(zhì)心坐標(biāo)路徑(圖4)，角度變化趨勢(shì)(圖5)及三維圖(圖6)，可以看出同樣是兩萬(wàn)步的結(jié)果，因?yàn)槌叽缱兇?，所以角度的變化?*1納米方塊緩和。

圖4

圖5

圖6

4.3 大納米棒(40*10)的一次隨機(jī)行走

將4*1納米棒換為形狀相同、100倍大的，觀察隨機(jī)行走，作出相應(yīng)模擬圖(圖7，8，9)。

圖7

圖8

圖9

可以明顯看到移動(dòng)范圍小得多，最遠(yuǎn)距離原點(diǎn)不到1。角度變化和之前相比也非常小，2萬(wàn)步最大轉(zhuǎn)動(dòng)不到十分之一弧度。

4.4 隨機(jī)行走與顆粒尺寸的關(guān)系

對(duì)比4*1納米棒和40*10納米棒，可以很明顯地看到：同樣長(zhǎng)寬比的納米棒，當(dāng)尺度變大時(shí)，布朗運(yùn)動(dòng)會(huì)變得緩和。這很容易理解，因?yàn)樵谇拔牡姆治鲋幸呀?jīng)指出步長(zhǎng)的平方與邊長(zhǎng)成反比。從物理角度上說(shuō)，考慮到隨機(jī)事件的抵償性，當(dāng)研究對(duì)象的尺寸變大時(shí)，分子撞擊的統(tǒng)計(jì)漲落將會(huì)越來(lái)越大的樣本容量所覆蓋。很顯然，我們只能在花粉等小顆粒上觀察到明顯的布朗運(yùn)動(dòng)，而較大的物體則幾乎體現(xiàn)不出布朗運(yùn)動(dòng)。

5 重要性質(zhì)

在繪制了納米棒布朗運(yùn)動(dòng)的圖像后，進(jìn)一步對(duì)數(shù)據(jù)進(jìn)行定量分析，可以得出一些重要結(jié)論。

5.1 平均行走距離

驗(yàn)證平均行走距離隨時(shí)間的演化。選取不同尺寸的納米棒分別進(jìn)行1000次8萬(wàn)步的模擬，并分別對(duì)這1000次中每一幀的質(zhì)心到原點(diǎn)距離的平方。

圖10

如圖所示，R2的期望與t大體上成正比。這正是數(shù)學(xué)中隨機(jī)行走理論的一個(gè)重要結(jié)論。[7]與勻速運(yùn)動(dòng)中R與t的正比關(guān)系不同，隨機(jī)行走的移動(dòng)效率將越來(lái)越低。在有取向的顆粒的布朗運(yùn)動(dòng)中，平方正比的關(guān)系并沒(méi)有被破壞，說(shuō)明轉(zhuǎn)動(dòng)并不影響平動(dòng)的整體效應(yīng)。

另一方面，對(duì)比4*1、8*2、40*10三個(gè)納米棒的斜率滿足10:5:1的關(guān)系，而且8*2和6*4幾乎重合，說(shuō)明在本文所使用的模型中，R2的期望與納米棒的周長(zhǎng)成反比，與形狀無(wú)關(guān)。

5.2 取向的自關(guān)聯(lián)系數(shù)

某一時(shí)刻的自關(guān)聯(lián)系數(shù)定義為此時(shí)取向方向單位向量與初始取向方向單位向量的內(nèi)積，表征對(duì)初始取向的“記憶”隨時(shí)間的衰減。由于初始條件都設(shè)為納米棒的取向指向x軸正方向，故計(jì)算自關(guān)聯(lián)系數(shù)即是計(jì)算取向的余弦的平均值。

圖11

左圖分別畫(huà)出了不同尺寸的納米棒在隨機(jī)行走中取向的自關(guān)聯(lián)系數(shù)的變化規(guī)律。定性來(lái)說(shuō)，前半段是上凹函數(shù)，從1衰減到0附近，之后開(kāi)始震蕩。尺寸越小的納米棒，衰減的速度越快，衰減完畢后震蕩得越劇烈。例如1*1納米方塊在幾千步長(zhǎng)之內(nèi)迅速衰減到0，而8*2/6*4納米棒則一直到走完8萬(wàn)步都只是衰減到0.6-0.7左右，40*10在8萬(wàn)步之內(nèi)幾乎沒(méi)有衰減。

用半對(duì)數(shù)坐標(biāo)重繪左圖。取縱坐標(biāo)為自關(guān)聯(lián)系數(shù)的負(fù)對(duì)數(shù)，得到右圖?？梢钥闯觯陂_(kāi)始震蕩之前，自關(guān)聯(lián)系數(shù)的負(fù)對(duì)數(shù)與t大體呈線性的線性關(guān)系。這表明在本文的模型中，自關(guān)聯(lián)函數(shù)隨時(shí)間呈指數(shù)衰減。

[r(0)·r(t)]=exp(-kt)

k是比例系數(shù)。

衰減到0附近后，開(kāi)始無(wú)規(guī)則震蕩。

容易理解，越小的納米棒，則衰減得越快，k越大。正如前文所述，越小的納米棒布朗運(yùn)動(dòng)越劇烈，則原來(lái)的取向在布朗運(yùn)動(dòng)中會(huì)快速丟失。而大的納米棒，例如本實(shí)驗(yàn)中40*10納米棒，布朗運(yùn)動(dòng)很小，取向的隨機(jī)改變也很慢，能夠保存很長(zhǎng)時(shí)間。至于宏觀尺度上的物體，布朗運(yùn)動(dòng)不可能導(dǎo)致其取向改變。

在原先取向的“記憶”衰減殆盡后，納米棒的取向就與原始取向不再有關(guān)聯(lián)性。不過(guò)自相關(guān)函數(shù)之后的噪聲又提出了一個(gè)新的問(wèn)題，這可能是計(jì)算機(jī)性能有限、樣本容量不夠大所致的統(tǒng)計(jì)漲落，不應(yīng)存在于解析結(jié)果中。

[參 考 文 獻(xiàn)]

[1] 郝柏林.布朗運(yùn)動(dòng)理論一百年[J].物理,2011,40(1):1-7.

[2] 謝征微,李玲.布朗運(yùn)動(dòng)的計(jì)算機(jī)模擬[J].四川師范大學(xué)學(xué)報(bào)(自然科學(xué)版),1997(1):97-100.

[3] 王萊.布朗運(yùn)動(dòng)和隨機(jī)行走的計(jì)算機(jī)仿真[J].工科物理,1998(6):47-49.

[4] 楊玉東,劉公召,徐菁華,楊林梅,李冬至.金納米棒:合成、修飾、自組裝、SERS及生物醫(yī)學(xué)應(yīng)用[J].中國(guó)科學(xué):化學(xué),2015,45(6):581-596.

[5] 王竹溪.統(tǒng)計(jì)物理學(xué)導(dǎo)論[M].北京：高等教育出版社,1965.

[6] 聶德明,林建忠.模擬顆粒布朗運(yùn)動(dòng)的格子Boltzmann模型[J].浙江大學(xué)學(xué)報(bào)(工學(xué)版),2009,43(8):1438-1442.

[7] 張波,張景肖.應(yīng)用隨機(jī)過(guò)程[M].北京：清華大學(xué)出版社,2004.