郭穎旦,丁望峰,楊建宋
(杭州師范大學(xué)理學(xué)院,浙江 杭州310036)
布朗運(yùn)動(dòng)是指懸浮于氣體或液體中的微小顆粒受到周圍分子熱運(yùn)動(dòng)產(chǎn)生的撞擊而做永不停息無規(guī)則運(yùn)動(dòng)的現(xiàn)象.1827年,植物學(xué)家羅伯特·布朗首次對(duì)布朗運(yùn)動(dòng)進(jìn)行了系統(tǒng)的描述,但直到1905年才由愛因斯坦根據(jù)統(tǒng)計(jì)理論給出定量的解釋[1]. 與此同時(shí),M. von Smoluchowski 也獨(dú)立發(fā)展出一套理論來解釋布朗運(yùn)動(dòng),并給出了相似的結(jié)果[2].1908年,Paul Langevin 采用了一種更為直觀的方法來描述布朗運(yùn)動(dòng),即在牛頓第二定律的方程中引入隨機(jī)力作用項(xiàng),即著名的“朗之萬方程”[3]:
其中Mi為布朗顆粒質(zhì)量,ri為位置矢量.等式左邊第一項(xiàng)是布朗顆粒受到的其它顆?;蛲鈭?chǎng)的作用力之和;第二項(xiàng)把周圍流體分子對(duì)布朗顆粒撞擊作用等效表示為一個(gè)隨機(jī)作用力;第三項(xiàng)是微粒在流體中運(yùn)動(dòng)受到的黏滯阻力,其中γi是顆粒i 的阻力系數(shù),對(duì)于半徑為a 的球形顆粒,按Stokes 公式計(jì)算其摩擦系數(shù)γ 為
其中η 是液體在一定溫度下的黏度.
根據(jù)隨機(jī)耗散理論,隨機(jī)作用力Ri(t)與摩擦系數(shù)滿足如下關(guān)系
其中kB是玻爾茲曼常數(shù),T 是環(huán)境溫度.
在實(shí)際模擬中,往往假定隨機(jī)作用力在不同時(shí)刻是不相關(guān)的,即式(4)可采用如下式
該等式把布朗粒子受到的隨機(jī)熱擾動(dòng)與體系溫度及摩擦系數(shù)關(guān)聯(lián)起來.
當(dāng)液體的黏滯阻力很大時(shí),布朗顆粒的質(zhì)量可以忽略,其運(yùn)動(dòng)軌跡類似于無質(zhì)量的粒子.這時(shí)式(1)等式右邊可近似等于零,即
該等式也被稱為位置朗之萬方程,是朗之萬方程的“過阻尼”形式[4].基于式(6)進(jìn)行的模擬稱為布朗動(dòng)力學(xué)(Brownian Dynamics,BD)模擬,而通過式(1)進(jìn)行的模擬稱為朗之萬動(dòng)力學(xué)(Langevin Dynamics,LD)模擬.
分子動(dòng)力學(xué)一般以原子為最小研究單元,在模擬過程中需考慮原子與原子之間的相互作用力.而在布朗運(yùn)動(dòng)的模擬過程中,分子水平的作用力被簡(jiǎn)化成一項(xiàng)隨機(jī)作用力,解決了分子動(dòng)力學(xué)直接以牛頓運(yùn)動(dòng)方程追蹤系統(tǒng)中所有布朗粒子與流體分子軌跡的困境,因而被廣泛應(yīng)用于膠體科學(xué)和生物物理領(lǐng)域[5]. 然而這種簡(jiǎn)化是基于大量流體分子撞擊前提下的,即需要一定的時(shí)間以保證足夠次數(shù)的撞擊發(fā)生,這也是布朗顆粒動(dòng)量發(fā)生明顯變化所需要的弛豫時(shí)間,定義為
最近實(shí)驗(yàn)研究顯示,在極小的時(shí)間尺度下,布朗顆粒受到流體分子的等效作用并非是完全隨機(jī)的白噪聲[6].計(jì)算機(jī)數(shù)值化模擬實(shí)際上是把連續(xù)的過程分割成時(shí)間間隔為Δt 的一系列狀態(tài)點(diǎn),通過計(jì)算從t 時(shí)刻到t+Δt 時(shí)刻狀態(tài)點(diǎn)的演化來模擬真實(shí)的物理過程.當(dāng)Δt?τp時(shí),把流體分子的作用等效成隨機(jī)作用項(xiàng)Ri(t)將不能反映真實(shí)物理過程. 而在實(shí)際模擬過程中,Δt 的取值往往由顆粒之間的勢(shì)能曲面形貌所決定,所以討論不同Δt 的取值對(duì)模擬結(jié)果的影響有著重要意義.
本文對(duì)單個(gè)布朗顆粒的運(yùn)動(dòng)軌跡分在不同時(shí)間步長Δt 下分別進(jìn)行了LD 模擬和BD 模擬,并對(duì)結(jié)果做出比較分析.根據(jù)愛因斯坦理論,通過測(cè)量布朗運(yùn)動(dòng)的軌跡可以計(jì)算出阿伏加德羅常數(shù)值,我們以此為依據(jù)對(duì)LD 和BD 模擬所得軌跡的物理真實(shí)性進(jìn)行驗(yàn)證.
這里記位矢ri的三維坐標(biāo)為(x1,x2,x3),三維速度矢量(v1,v2,v3),則在根據(jù)式(1),xα(α=1,2,3)方向上粒子i 從時(shí)間t 到t+Δt 基于Euler 算法的演化方程為[7]
相應(yīng)的位置朗之萬方程可以寫成
從上式可以看出,位置朗之萬方程中并不出現(xiàn)速度,粒子下一步的走向完全取決于粒子此刻受到的作用力(當(dāng)Δt 足夠小時(shí)可看成恒力)和熱擾動(dòng)產(chǎn)生的隨機(jī)作用力.
為了滿足式(3)和(4)要求,這里的隨機(jī)作用力Rα,i(t)可以取一個(gè)期望值為0、標(biāo)準(zhǔn)差為的高斯分布,記為
如果再記布朗擴(kuò)散系數(shù)D 為
則式(9)可寫成
通過上式不難發(fā)現(xiàn),在不考慮顆粒間相互作用和外場(chǎng)影響時(shí),位移的平均平方差滿足布朗運(yùn)動(dòng)理論的推導(dǎo)結(jié)果
相對(duì)而言,從式(7 -8)并不能直接得到位移的平方差表達(dá)式,我們將在下面具體的模擬中做進(jìn)一步分析.
在下面的結(jié)果與討論中,LD 模擬與BD 模擬將分別基于式(8 -9)和式(10)來實(shí)現(xiàn).
我們選取懸浮在室溫純水中的球形SiO2顆粒作為研究對(duì)象,模擬其在二維平面內(nèi)的布朗運(yùn)動(dòng)軌跡,具體參數(shù)見表1.
為了便于比較,每次模擬均采用一串相同的“偽隨機(jī)數(shù)”,即顆粒在t 時(shí)刻受到相同的隨機(jī)作用力. 如果模擬的參數(shù)相同,則布朗顆粒的“隨機(jī)”運(yùn)動(dòng)軌跡完全相同.
計(jì)算過程中積分的時(shí)間步長選擇一般由顆粒之間的作用勢(shì)所決定,越“陡峭”的勢(shì)能曲面需要越小的時(shí)間步長.由于這里只考慮單個(gè)布朗顆粒的運(yùn)動(dòng),不涉及顆粒之間的作用力,因而時(shí)間步長可以在一個(gè)較大的范圍內(nèi)選擇.圖1 是在保證其他參量相同的前提下,使用不同的模擬時(shí)間步長Δt 分別進(jìn)行LD 和BD 模擬得到的單顆粒布朗運(yùn)動(dòng)軌跡圖,每條軌跡代表模擬中最初250 步隨機(jī)行走.
表1 SiO2 顆粒在水中作布朗運(yùn)動(dòng)的模擬參數(shù)Tab. 1 Parameters for the Brownian simulation of SiO2 spheres in water
由圖1 可見,隨著時(shí)間步長Δt 的增大,模擬的總時(shí)間增大,因而行走的空間距離也相應(yīng)增大.不同的是,BD 軌跡除了空間尺度的變化外,軌跡形狀完全相同.這從式(10)中也可以看出,時(shí)間步長Δt 僅線性地改變了行走的步長.而對(duì)于LD 模擬而言,Δt 越小,顆粒本身的慣性質(zhì)量越突顯,模擬軌跡就越接近于宏觀物體的實(shí)際運(yùn)動(dòng)軌跡,是一條較為平滑的曲線.隨著Δt 的增大,像所有的分子動(dòng)力學(xué)模擬一樣,每步行走的步長增加,LD 的軌跡變得越來越跳躍,并最終導(dǎo)致模擬數(shù)據(jù)溢出,如圖1(A)最下面一幅圖所示. 對(duì)更多的Δt 值進(jìn)行比較,我們發(fā)現(xiàn)當(dāng)Δt 接近于顆粒的弛豫時(shí)間τp時(shí),BD 和LD 軌跡最相似.
BD 是LD 的簡(jiǎn)化,進(jìn)行BD 模擬的前提條件是物理過程必須滿足布朗運(yùn)動(dòng)的特征,只要條件滿足,模擬的時(shí)間步長可由其它因素決定.而LD 模擬的Δt 取值受實(shí)際物理?xiàng)l件制約,存在上限值.這就會(huì)出現(xiàn)這樣的矛盾:隨機(jī)作用項(xiàng)Ri(t)引入的條件是時(shí)間步長Δt 足夠大(Δt ?τp),使流體分子的作用可以簡(jiǎn)化為一個(gè)隨機(jī)作用力;而對(duì)于普通的布朗顆粒,動(dòng)力學(xué)模擬中勢(shì)能曲面要求的時(shí)間步長遠(yuǎn)小于其弛豫時(shí)間τp.所以在實(shí)際計(jì)算中,Δt 的取值往往比弛豫時(shí)間小很多.那么這樣的處理方式是否科學(xué),下面通過對(duì)輸出軌跡的定量計(jì)算進(jìn)行討論.
1905年,愛因斯坦首次發(fā)表了布朗運(yùn)動(dòng)的定量分析理論,在文章的最后,他給出了驗(yàn)證理論正確性的方法,即
其中NA是阿伏加德羅常數(shù),Δt 是運(yùn)動(dòng)軌跡上鄰近兩個(gè)記錄點(diǎn)的時(shí)間間隔,R 是氣體常數(shù),T 是實(shí)驗(yàn)溫度,η 是實(shí)驗(yàn)溫度下液體的黏度,a 是球形微粒的半徑.通過測(cè)量軌跡圖上顆粒在x 方向上位移的平均平方值就可以計(jì)算出阿伏加德羅常數(shù)的實(shí)驗(yàn)測(cè)量值.這里用此公式對(duì)模擬軌跡圖進(jìn)行定量分析,并用計(jì)算所得的NA值與標(biāo)準(zhǔn)值相比較來對(duì)模擬的物理真實(shí)性進(jìn)行評(píng)判.
這里所用的“偽隨機(jī)數(shù)列”以每次運(yùn)行時(shí)計(jì)算機(jī)時(shí)間為種子產(chǎn)生,從而保證了軌跡的隨機(jī)性.對(duì)于每個(gè)時(shí)間步長Δt,我們共計(jì)算了n = 105步隨機(jī)行走的平均位移平方<λx2 >,并得到了如圖2 所示的計(jì)算結(jié)果. 需要注意的是,當(dāng)Δt >300 ns 時(shí),LD 模擬數(shù)據(jù)已經(jīng)溢出.所以對(duì)于此時(shí)間步長范圍內(nèi)的LD 模擬,采用Δt = Δt′ × N的模式,即實(shí)際模擬時(shí)間步長為較小的Δt′,但間隔N 步輸出一次坐標(biāo),這樣輸出軌跡上點(diǎn)之間的時(shí)間步長為Δt.例如Δt =1000 ns 時(shí),實(shí)際Δt′ =10 ns,N=100.計(jì)算表明,只要Δt 的取值不變,不同Δt′和N 的取值對(duì)結(jié)果幾乎沒有影響.
從圖2 可以看到,BD 軌跡與Δt 無關(guān),由此計(jì)算得到的阿伏加德羅常數(shù)值非常接近其標(biāo)準(zhǔn)值6.022 ×1023.而LD 模擬在Δt 較小時(shí)得到的NA值明顯偏大,這是因?yàn)樵诤芏痰臅r(shí)間內(nèi)要保證流體分子對(duì)布朗顆粒作用的熱力學(xué)統(tǒng)計(jì)性,勢(shì)必需要更大量的分子,這就表現(xiàn)為NA值偏大.上面提到,布朗顆粒的弛豫時(shí)間τp是其受微觀分子撞擊引起動(dòng)量變化的特征時(shí)間,當(dāng)時(shí)間步長小于該特征時(shí)間時(shí),就需要更大量的分子撞擊使之產(chǎn)生動(dòng)量上的變化,即對(duì)應(yīng)更大的阿伏加德羅常數(shù)值;當(dāng)時(shí)間步長大于τp時(shí),足夠數(shù)量的分子撞擊符合熱力學(xué)統(tǒng)計(jì)要求,由此可得到正確的NA值.圖2 中的LD 曲線的趨勢(shì)也證實(shí)了布朗顆粒的弛豫時(shí)間(τp=129.06 ns)是NA取值是否正確的“分水嶺”.同時(shí),計(jì)算結(jié)果也證實(shí)了在LD 模擬中,即使時(shí)間步長小于弛豫時(shí)間,對(duì)于長時(shí)間的宏觀觀測(cè)而言該模擬也是符合真實(shí)物理規(guī)律的.
圖2 用不同時(shí)間步長分別進(jìn)行LD 和BD 模擬,并對(duì)其輸出軌跡進(jìn)行計(jì)算得到的阿伏加德羅常數(shù)值NAFig. 2 Values of the Avogadro constant calculated from traces of Langevin dynamics and Brownian dynamics simulations performed with different time-steps
本文從朗之萬方程出發(fā),分別用朗之萬動(dòng)力學(xué)(LD)和布朗動(dòng)力學(xué)(BD)方法對(duì)單顆粒的布朗運(yùn)動(dòng)進(jìn)行了數(shù)值化模擬,并重點(diǎn)討論了不同時(shí)間步長對(duì)模擬結(jié)果的影響.LD 是在牛頓動(dòng)力學(xué)的基礎(chǔ)上加入了隨機(jī)作用項(xiàng)以等效流體分子對(duì)布朗顆粒的碰撞,必須滿足熱力學(xué)的大數(shù)量統(tǒng)計(jì).當(dāng)時(shí)間步長太短時(shí),分子碰撞次數(shù)不足以等效成隨機(jī)作用,所以在較小的時(shí)間尺度下,這樣的等效是不符合物理事實(shí)的;然而在宏觀時(shí)間尺度下,發(fā)現(xiàn)LD 模擬軌跡并不違背物理規(guī)律.BD 是對(duì)LD 的簡(jiǎn)化,其成立條件是體系滿足布朗運(yùn)動(dòng)的所有特征,因而在具體模擬過程中與時(shí)間步長無關(guān).
[1]Einstein A. Investigations on the Theory of the Brownian Movement[M]. New York:Dover Pulieation Inc,1956.
[2]Smoluchowski M von. Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen[J]. Annalen der Physik,1906,326(14):756-780.
[3]Langevin P. Sur la théorie du mouvement brownien[J]. CR Acad Sci Paris,1908,146:530-533.
[4]Coffey W T,Waldron J T,Kalmykov Y P. The Langevin Equation[M]. Singapore:World Scientific,2004.
[5]王子瑜,曹恒光.布朗運(yùn)動(dòng)、朗之萬方程式、與布朗動(dòng)力學(xué)[J].物理雙月刊,2005,27(3):456-461.
[6]Li T C,Raizen M G. Brownian motion at short time scales[J]. Annalen der Physik,2013,525(4):281-295.
[7]Leach A R.分子模擬的原理和應(yīng)用[M].2 版.倫敦:世界圖書出版公司,2001.