朗之萬方程在布朗運(yùn)動(dòng)數(shù)值模擬中的時(shí)間尺度分析

郭穎旦，丁望峰，楊建宋

(杭州師范大學(xué)理學(xué)院，浙江 杭州310036)

0 引 言

布朗運(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)證.

1 朗之萬方程的數(shù)值化

這里記位矢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).

2 結(jié)果與討論

我們選取懸浮在室溫純水中的球形SiO2顆粒作為研究對(duì)象，模擬其在二維平面內(nèi)的布朗運(yùn)動(dòng)軌跡，具體參數(shù)見表1.

2.1 不同模擬時(shí)間步長對(duì)比

為了便于比較，每次模擬均采用一串相同的“偽隨機(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)行討論.

2.2 布朗運(yùn)動(dòng)軌跡的定量分析

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

3 結(jié) 論

本文從朗之萬方程出發(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.