劉心爽 張超越

（安慶師范大學物理系，安徽 安慶 246133）

0 引言

計算物理是運用計算機作為工具，通過數(shù)值計算、數(shù)據(jù)處理、模擬仿真等手段解決復(fù)雜物理問題的一門應(yīng)用科學。計算物理可以處理理論物理無解析解問題，同時也可以處理實驗物理上因?qū)嶒灄l件限制而不能處理的問題。近些年來，隨著計算機技術(shù)的快速發(fā)展和學科間的進一步滲透，計算物理得以蓬勃發(fā)展，計算物理也成為國內(nèi)高校物理學及相關(guān)專業(yè)本科生的專業(yè)課程。

目前國內(nèi)計算物理教學主要集中在運用計算軟件（通常為Matlab）指令庫進行數(shù)值微分與數(shù)值積分、解常微分方程、解偏微分方程以及一些數(shù)值模擬方法（如蒙特卡羅方法）[1，2]。計算物理教學的特點是學科間交叉性強，能夠在一定程度上幫助學生對線性代數(shù)、插值、非線性方程組和計算機編程等一系列知識的掌握，同時，能夠很好地培養(yǎng)學生解決問題的能力，可見計算物理在本科教學中的重要性。

然而，在傳統(tǒng)計算物理教學中存在著一些不足。首先，教學方式單一。傳統(tǒng)教學中采取教師到學生的“單向傳授”教學模式，并不適宜于計算物理的教學。其次，教學內(nèi)容前沿性不足。目前的計算物理教學主要集中在運用數(shù)值方法求解典型的物理問題，幫助學生掌握數(shù)值處理方法及方程解法，同時提升編程能力。但是，學生對計算物理的發(fā)展和前沿領(lǐng)域知之甚少，特別是對物理、化學、生物、計算機等多學科的交叉研究領(lǐng)域和新興發(fā)展的研究方法。為了拓寬學生的視野，為了改進諸如此類的傳統(tǒng)計算物理教學上的不足、激發(fā)學生在計算物理學習中的積極性與主動性、培養(yǎng)學生的科學思維、發(fā)揮學生的主體性并提高有效解決實際問題的能力、開拓學生的視野從而幫助學生更加深入了解計算物理這門課程，筆者將多粒子碰撞動力學方法引入計算物理的教學之中，以此在教學過程中進行補充說明，提高課堂教學效果，從而滿足課堂教學效果。

多粒子碰撞動力學（Multi-particle Collision Dynamics：MPCD）方法是一門新興的計算方法，由Malevanets和Kapral于1999年提出，又被稱為隨機旋轉(zhuǎn)動力學，廣泛用于軟物質(zhì)體系的動力學模擬研究[3-5]。在軟物質(zhì)系統(tǒng)的研究中，復(fù)雜的溶劑環(huán)境是研究過程中面臨的一個重要的問題，如膠體凝膠粒子系統(tǒng)，粒子間除了存在靜電作用力以及分子間作用力和體積排斥作用外，還存在著長程流體力學作用（Hydrodynamic Interaction，HI）[6]。這種長程的相互作用是由溶劑環(huán)境所產(chǎn)生的，并伴隨著粒子運動的始終。由此可見，流體力學效應(yīng)對系統(tǒng)動力學行為的影響是一直存在的，應(yīng)該被考慮進來。而在很多時候，為了簡化動力學過程，通常會忽略流體力學效應(yīng)，將環(huán)境的影響用噪聲項來替代。多粒子碰撞動力學方法能很好地解決這個問題。

1 粒子碰撞動力學方法的基本原理

在多粒子碰撞動力學方法中，為了考慮溶劑效應(yīng)，同時又不會因粒子數(shù)的提高在很大程度上提高計算量，溶劑分子被抽象成粒子間無體積排斥的點狀粒子，粒子以連續(xù)的速度分布在連續(xù)空間運動。MPCD模擬方法中，體系的動力學演化過程分為粒子位置演化和速度演化兩步，筆者分別稱之為粒子的流動和粒子的碰撞。

1.1 粒子流動

在t時刻，粒子t的位置為ri（t），速度為vi（t），經(jīng)過ΔtMPC時間后粒子的位置更新為：

進而遍歷系統(tǒng)中所有粒子。

1.2 粒子碰撞

由于粒子被抽象成點狀粒子，粒子間無相互作用力，為了實現(xiàn)粒子間的動量交換，在MPCD方法中采用了一種粗?；呐鲎蔡幚聿僮?。具體來說，首先將系統(tǒng)分割成一系列網(wǎng)格并做好編號，系統(tǒng)中粒子劃分到具體網(wǎng)格中，同一網(wǎng)格中的粒子間相互碰撞交換動量，不同網(wǎng)格間粒子之間互不干擾不進行動量交換。在模擬過程中，通常選定常數(shù)作為網(wǎng)格的邊長，網(wǎng)格中粒子平均數(shù)目M∈[3，20]。在實際碰撞過程中，由于每次進出格子的粒子數(shù)目并不一定相等，因而網(wǎng)格中真實的粒子數(shù)目Nc并不固定，存在一定漲落，一般在平均值M附近波動，其原因是粒子在碰撞后進入鄰近格子之中。在t時刻，由碰撞網(wǎng)格ξ中所有粒子的瞬時速度vi可計算出網(wǎng)格的質(zhì)心速度vξcm=t），則在t+ΔtMPC時刻粒子的速度大小為：

式中：δvi（t）=vi（t）-vξcm（t）為粒子在碰撞過程中速度的變化量。(R為旋轉(zhuǎn)矩陣，在每個時間步長，它將每個單元格的速度圍繞隨機生成的軸旋轉(zhuǎn)一個角度α。需要注意的是，當溶劑粒子平均自由程λ=ΔtMPCkB為玻爾茲曼常數(shù)，T為溫度，m為溶劑粒子質(zhì)量）相對網(wǎng)格尺寸比較小時，粒子需要經(jīng)過反復(fù)的碰撞才會進入或者流出網(wǎng)格，這種反復(fù)的碰撞會使得粒子間的時間關(guān)聯(lián)加強，導(dǎo)致伽利略不變性被破壞。為了保證伽利略不變性，通常在每次碰撞之前對網(wǎng)格進行一次隨機移動，這樣粒子會重新劃分格子，與新的格子中的粒子進行碰撞，碰撞結(jié)束后格子再移動回來。這種移動網(wǎng)格的方式，既可以保證伽利略不變性，還可以加速動量在格子間的傳遞。當λ的值大于a0/2時，不需要移動網(wǎng)格。

在MPCD算法中，通過粒子流動和碰撞步的交替演化，粒子之間實現(xiàn)了動量的傳遞，同時每個格子中的動量和能量守恒，因而體系可以通過溶劑產(chǎn)生流體力學相互作用，進而通過此方法可以模擬體系中的流體力學。

當系統(tǒng)中除了溶劑粒子外，還存在溶質(zhì)粒子（包括膠體粒子、納米粒子、高分子鏈等）時，特別是需要研究溶劑的流體力學效應(yīng)對溶質(zhì)的動力學行為影響時，需要在溶質(zhì)和溶劑間建立相互作用關(guān)系，溶質(zhì)粒子與溶劑粒子間需要進行耦合。為了實現(xiàn)溶質(zhì)與溶劑間耦合，Kapral等提出了多粒子碰撞結(jié)合粗?；肿觿恿W（MD-MPCD）的方法。在該方法中，溶質(zhì)與溶劑的動力學演化過程是分開的，溶質(zhì)粒子采用分子動力學進行演化，而溶劑粒子采用上面介紹的流動-碰撞交替的方式進行，這種方式有效地避免了由溶劑粒子數(shù)目產(chǎn)生的分子動力學的計算，極大地節(jié)約了計算成本。需要注意的是ΔtMD＜＜ΔtMPC，也就是說，在動力學模擬過程中，當溶質(zhì)粒子進行多次分子動力學演化后，溶劑分子才進行一次演化。對于溶劑和溶質(zhì)間的耦合，通常采用碰撞耦合的方式，在該過程中，和溶劑粒子一樣，溶質(zhì)粒子也被視為點狀粒子，動量交換過程中溶質(zhì)和溶劑共同參與碰撞。具體說來就是，在碰撞步前，同時將溶質(zhì)粒子和溶劑粒子視為點狀粒子劃分到網(wǎng)格之中，當網(wǎng)格ξ中含有Nc個質(zhì)量為m的溶劑粒子和Nm個質(zhì)量為M的粒子時，網(wǎng)格的質(zhì)心速度為其中，vi和Vj分別為第i號溶劑粒子和第j號溶質(zhì)粒子的速度。碰撞時，按照公式（2）的碰撞方式，溶質(zhì)粒子和溶劑粒子隨機碰撞，從而實現(xiàn)溶質(zhì)和溶劑間動量交換，進而達到溶質(zhì)和溶劑耦合的效果。該方法相比于相互作用力耦合和熱邊界耦合簡單實用，已被廣泛應(yīng)用于軟物質(zhì)系統(tǒng)。

2 多粒子碰撞動力學方法的應(yīng)用

作為一個具體例子，首先運用多粒子碰撞動力學模擬了二維泊肅葉流，再將高分子鏈置于二維泊肅葉流中。如圖1所示，在x方向選取Lx=100a0，滿足周期邊界條件，y方向兩平板間距設(shè)定為Ly=30a0。當流體粒子在x方向受到重力場g驅(qū)動時，由于在邊界處（y=0和y=Ly），流體粒子與邊界的摩擦效應(yīng)使得沿x方向速度vx減小為零，而在兩平行板中心，vx取得最大值。兩平行板間流場形成拋物形的輪廓線，理論上滿足v（xy）=Ly-y）y，其中，η為流體粘滯系數(shù)，g為施加在流體上的力產(chǎn)生的加速度大小。

圖1 泊肅葉流場示意圖

在模擬過程中，溶劑粒子和平行板之間采用無滑邊界。模擬結(jié)果如圖2所示，在不同的驅(qū)動力作用下，在邊界處流場速度趨向于0，在平行板中間位置流體速度達到最大值，這與理論上保持一致。模擬中，逐漸增大重力場大小0.000 5、0.00 1、0.001 5，發(fā)現(xiàn)單調(diào)上升，三條曲線對應(yīng)的最大值分別為0.044 9、0.090 6、0.133 1，值的大小比例與保持一致。實際中，通過模擬流體在管道中的速度場，可以運用理論公式間接的測量流體黏滯系數(shù)，這里就不做過多陳述。

圖2 不同驅(qū)動力作用下泊肅葉流場平均速度輪廓線（自上而下g=0.000 5、0.00 1、0.001 5）

當該系統(tǒng)中含有溶質(zhì)粒子時（如高分子鏈），在通道中，高分子鏈在泊肅葉流的驅(qū)動下將沿著流體的驅(qū)動方向運動，高分子鏈采用鏈—珠模型，臨近高分子鏈單體之間通過有限張伸非線性彈性勢（FENE）連接。

式中：k是彈性耦合常數(shù)，Rmax表示鍵長，即高分子鏈單體間所能允許的最遠距離。高分子鏈單體之間的即體積排斥作用通過WCA（Weeks-Chandler-Andersen）勢表示。

式中：ε表示勢能的強度大??；σ為高分子鏈單體直徑；rij為任意兩單體之間的距離。關(guān)于高分子鏈，采用標準的分子動力學進行演化。

根據(jù)Verlet算法求解（5）式牛頓運動方程為：

式中，Ri（t）、Vi（t）、Fi（t）分別為i號溶質(zhì)粒子的位置、速度、所受力。溶質(zhì)與溶劑耦合時，溶質(zhì)和溶劑之間相互碰撞并交換動量。

筆者研究了高分子鏈在平行板管道中的密度分布，如圖3所示，橫坐標為管道寬度，縱坐標為高分子鏈的分布。觀察到，高分子鏈主要分布在管道中間，很難在管道壁中長時間停留，在邊界處的分布趨向于零。有趣的是，高分子鏈在管道中的密度分布不是呈現(xiàn)單峰分布的，在管道中間位置的分布不是最多，在管道中心的附近呈現(xiàn)雙峰分布。

圖3 高分子鏈在管道中的密度分布

將這些偏離中心的峰值歸因于相互競爭的遷移效應(yīng)。其一，二維系統(tǒng)中泊肅葉流的速度分布是拋物形的，中間快兩邊慢，這種剪切的速度梯度在很大程度上影響高分子鏈的構(gòu)型（受到高剪切速率的高分子會被拉伸的更多），使得高分子鏈活性產(chǎn)生一個運動梯度（高分子鏈擴散）。根據(jù)動力學理論，高分子鏈擴散率的梯度會導(dǎo)致鏈遷移，并且鏈遷移的方向是遠離泊肅葉流的通道中心[7]。其二，熱擴散和流體力學效應(yīng)使得高分子鏈向管道中間運動[8]。正是由于這兩種效應(yīng)的相互競爭，使得高分子鏈在管道中的分布呈現(xiàn)有趣的雙峰分布。

3 結(jié)語

針對安慶師范大學物理學專業(yè)教學的實際情況，對計算物理教學內(nèi)容進行了改進完善。在教學中引入多粒子碰撞動力學方法的介紹，并結(jié)合具體例子—泊肅葉流的產(chǎn)生來說明該方法的實際應(yīng)用。通過模擬泊肅葉流流場分布以及流場中高分子鏈的分布來刻畫流體力學效應(yīng)?！岸嗔Ｗ优鲎矂恿W”這一新興方法的引入，將對構(gòu)建更加完整的計算物理教學體系發(fā)揮重要的作用。更重要的是極大地拓展了學生的知識面和認知，深化了學生從介觀體系對統(tǒng)計物理問題的理解，激發(fā)了學生對新事物的好奇心，培養(yǎng)了學生的科學思維，充分體現(xiàn)了教學中學生為中心的主線和持續(xù)改進的理念。