多粒子碰撞動力學方法在計算物理教學中的應用

劉心爽 張超越

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

0 引言

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

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

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

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

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

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

1.1 粒子流動

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

進而遍歷系統中所有粒子。

1.2 粒子碰撞

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

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

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

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

2 多粒子碰撞動力學方法的應用

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

圖1 泊肅葉流場示意圖

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

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

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

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

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

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

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

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

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

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

3 結語

針對安慶師范大學物理學專業教學的實際情況，對計算物理教學內容進行了改進完善。在教學中引入多粒子碰撞動力學方法的介紹，并結合具體例子—泊肅葉流的產生來說明該方法的實際應用。通過模擬泊肅葉流流場分布以及流場中高分子鏈的分布來刻畫流體力學效應。“多粒子碰撞動力學”這一新興方法的引入，將對構建更加完整的計算物理教學體系發揮重要的作用。更重要的是極大地拓展了學生的知識面和認知，深化了學生從介觀體系對統計物理問題的理解，激發了學生對新事物的好奇心，培養了學生的科學思維，充分體現了教學中學生為中心的主線和持續改進的理念。