分子動力學(xué)模擬高溫高壓下水的粘度

張鈞惠，張欣，蘭美華

（1.黑龍江八一農(nóng)墾大學(xué)工程學(xué)院，大慶163319；2.大慶鉆探工程公司鉆井工程技術(shù)研究院；3.大慶鉆探工程公司地質(zhì)錄井一公司）

將旋轉(zhuǎn)矩陣和相關(guān)物理量分別由代換變量q1、q2、q3、q4表示。再結(jié)合剛體轉(zhuǎn)動的歐拉方程，經(jīng)過數(shù)學(xué)推導(dǎo)，可以得到加速度與四參數(shù)系統(tǒng)、角速度之間的相互關(guān)系如下：

分子動力學(xué)模擬高溫高壓下水的粘度

張鈞惠1，張欣2，蘭美華3

（1.黑龍江八一農(nóng)墾大學(xué)工程學(xué)院，大慶163319；2.大慶鉆探工程公司鉆井工程技術(shù)研究院；3.大慶鉆探工程公司地質(zhì)錄井一公司）

由于高溫高壓的極端條件實驗室很難實現(xiàn)，因此應(yīng)用分子動力學(xué)方法模擬計算水在高溫高壓下的粘度成為一種重要的手段。模擬實驗基于分子動力學(xué)方法，針對水分子的特殊結(jié)構(gòu)，引入Ewald法求長程作用力和四元素法，采用等溫等壓條件，最終模擬得到不同溫度、壓強(qiáng)下水的粘度并與實驗值進(jìn)行比較分析。

分子動力學(xué)；長程作用力；等溫等壓；四元素法；剪切應(yīng)力；粘度

分子動力學(xué)模擬（molecular dynamics simulation，MD）以經(jīng)典牛頓力學(xué)為理論基礎(chǔ)，系統(tǒng)中各粒子的運動服從拉格朗日方程或哈密頓方程，通過粒子運動與軌跡的聯(lián)系模擬系統(tǒng)中粒子的相關(guān)熱力學(xué)性質(zhì)。由于分子動力學(xué)模擬可以有效彌補(bǔ)理論分析和實驗觀察的不足，能夠獲得分子運動的各種微觀特性，揭示各種現(xiàn)象的本質(zhì)。因此，分子動力學(xué)模擬逐漸成為國內(nèi)外學(xué)術(shù)界的研究熱點，在材料科學(xué)、物理學(xué)、化學(xué)、生物學(xué)等領(lǐng)域得到廣泛應(yīng)用［1-5］，已逐步成為第三種科學(xué)研究手段。

水在農(nóng)業(yè)、工業(yè)等各個領(lǐng)域都起著不可替代的作用［6］，水在各種極端條件下性質(zhì)的研究迫在眉睫。然而超高溫高壓、超臨界狀態(tài)等極端條件在實驗室很難甚至不能達(dá)到，因此水的分子動力學(xué)模擬便成為研究極端條件下水的相關(guān)性質(zhì)的唯一方法。目前，針對極端條件下的分子動力學(xué)模擬研究日益增多，馬靜等［7］模擬了超（近）臨界水溶劑特性，給出了水在常溫、近臨界和超臨界狀態(tài)的內(nèi)聚能、溶解度參數(shù)與溫度、壓力的關(guān)系。孫煒等［8］模擬了溫度為673.15～873.15 K，壓力為22.1～131.3 MPa條件下水的密度和自擴(kuò)散系數(shù)。然而，對高溫高壓下水的粘度的研究仍未進(jìn)行，針對以上情況，推導(dǎo)了分子動力學(xué)模擬時長程作用力的計算公式和粘度計算公式，并分別模擬了溫度為473.16 K、573.16 K，壓強(qiáng)為10 MPa、20 MPa的高溫高壓條件下水的粘度。

1 模型建立及理論基礎(chǔ)

1.1 初始條件、邊界條件

采用各端自由的三維模型，系綜為NPT系綜。系統(tǒng)的初始位置采用面心立方排布，初始速度滿足Maxwell-Boltzmann分布，周期性邊界條件。

1.2 勢函數(shù)

采用TIP4P型勢函數(shù)［9］：

TIP4P型勢函數(shù)參數(shù)為：εoo=1.077 2×10-21J，σoo= 0.315 4 nm，rOH=0.095 72 nm，rOH=0.015 nm，∠HOH= 104.52°，ε0=8.854×10-12C2·J-1·m-1，qH=0.52 e，qM= -1.04 e，e=1.602 19×10-19C。

1.3 分子間相互作用力

對勢函數(shù)求導(dǎo)可得相鄰水分子間的非鍵結(jié)力和靜電力。由疊加原理可得水分子i在x方向所受合力為

式中，第一項為非鍵結(jié)力，在模擬過程中只需考慮截斷半徑內(nèi)的非鍵結(jié)力即可；第二項為靜電力，由于靜電力隨半徑收斂速度較慢，因此引入長程作用力，采用Ewald加和法來計算分子間的長程作用力。

考慮邊長為L，含有N個水分子的立方體元胞盒子，則系統(tǒng)盒內(nèi)的水分子以Ewald方法計算的系統(tǒng)電荷作用勢的表達(dá)式為：

式中的勢能分為三項為Ur、Uk、Us，其中，Ur為系統(tǒng)中水分子間的相互作用，n=L（nx，ny，nz），nx，ny，nz=0， ±1，±2…，erfc為余誤差函數(shù)，α為控制收斂速度的經(jīng)驗參數(shù)，取α=5/L；Uk為系統(tǒng)中水分子與其他水分子的鏡像水分子間的相互作用，其中k=2π（nx，ny，nz）/L，為保證模擬結(jié)果的準(zhǔn)確性和提高模擬效率，取nx，ny，nz=0，±1，±2，±3。

分別對Ur、Uk、Us求導(dǎo)，可得水分子i在x方向所受的長程作用力Fix的表達(dá)式：

1.4 等溫等壓調(diào)節(jié)

在每個時間步，速度乘以速度標(biāo)定因子，逐步調(diào)整系統(tǒng)溫度達(dá)到恒定。速度標(biāo)定因子的表達(dá)式為：

式中，β為速度標(biāo)定因子；波爾茲曼常數(shù)k=1.38× 10-23J·K-1；T為溫度；水分子質(zhì)量m=2.987×10-23g。

采用Anderson方法對系統(tǒng)壓力進(jìn)行調(diào)整。在每個時間步，將水分子質(zhì)心坐標(biāo)乘以一個系數(shù)調(diào)整系統(tǒng)體積，進(jìn)而達(dá)到調(diào)整系統(tǒng)壓強(qiáng)的目的。

引入變量V后，系統(tǒng)的Lagrange方程調(diào)整為：

采用Velocity-Verlet算法［10］分別對S，V，s˙，V˙進(jìn)行離散，可得等壓調(diào)整的模擬計算公式：

1.5 四元素法

考慮水分子的平移運動和自身的旋轉(zhuǎn)運動建立全局坐標(biāo)系oxyz和原點在水分子質(zhì)心上的局部坐標(biāo)系o′x′y′z′。由理論力學(xué)可知，旋轉(zhuǎn)變換矩陣中含有大量歐拉角的三角函數(shù)，而計算機(jī)計算三角函數(shù)非常耗機(jī)時，因此引入四個與歐拉角θ、φ、φ相關(guān)的變量q1、q2、q3、q4，

這四個變量q1、q2、q3、q4滿足

將旋轉(zhuǎn)矩陣和相關(guān)物理量分別由代換變量q1、q2、q3、q4表示。再結(jié)合剛體轉(zhuǎn)動的歐拉方程，經(jīng)過數(shù)學(xué)推導(dǎo)，可以得到加速度與四參數(shù)系統(tǒng)、角速度之間的相互關(guān)系如下：

1.6 粘度計算

在模擬過程中，計算并記錄系統(tǒng)的應(yīng)力數(shù)據(jù)后，通過對時間取平均的方法求得物質(zhì)的粘度，以下是粘度的計算公式：

式中，ηs為粘度為系統(tǒng)的應(yīng)力分量，當(dāng)α=β時，Pαβ為正應(yīng)力；當(dāng)α≠β時，Pαβ為剪切應(yīng)力，〈Pαβ（0）Pαβ（t）〉表示應(yīng)力的自相關(guān)函數(shù)。

把n組應(yīng)力模擬數(shù)據(jù)平分為兩組：Pαβ（1）、Pαβ（2）…Pαβ（n/2）和Pαβ（n/2+1）、Pαβ（n/2+2）…Pαβ（n）。令Qαβ（t）=〈Pαβ（0）Pαβ（t）〉，經(jīng)過離散計算得到模擬計算的粘度公式：

由于粘度是一個宏觀統(tǒng)計量，需要進(jìn)行足夠時間的循環(huán)，記錄大量的剪切應(yīng)力值，才能得到比較準(zhǔn)確的粘度值。

2 模型結(jié)果與分析

選取7×7×7的分子動力學(xué)模型，系統(tǒng)中含有1 688個水分子，取截斷半徑取為2.5σoo。實例分別模擬了溫度為573.16 K、473.16 K，壓強(qiáng)為10 MPa、20 MPa條件下水的粘度。模擬先經(jīng)過10萬步的弛豫過程，使分子動力學(xué)模擬系統(tǒng)達(dá)到平衡，再循環(huán)30萬步統(tǒng)計平均正應(yīng)力和平均剪切應(yīng)力數(shù)據(jù)，并根據(jù)應(yīng)力與粘度的關(guān)系計算水在不同溫度和壓強(qiáng)下的粘度。在對水的應(yīng)力計算時，只統(tǒng)計距離各邊界0.5α0以內(nèi)區(qū)域的受力。

圖1 平均正應(yīng)力Fig.1Average normal stress

圖2 平均剪切應(yīng)力Fig.2Average shear stress

圖1、2分別給出不同溫度、壓強(qiáng)下，系統(tǒng)平均正應(yīng)力、平均剪切應(yīng)力的變化情況。由圖可得：在四種給定的溫度、壓強(qiáng)條件下，系統(tǒng)的平均正應(yīng)力和平均剪切應(yīng)力在模擬的初始階段波動較大，隨著模擬時間步的推移，波動逐漸變小，最終，平均正應(yīng)力逐步趨向期望的壓強(qiáng)值10 MPa和20 MPa，平均剪切應(yīng)力在0 MPa附近上下波動。

圖3給出不同溫度、壓強(qiáng)下水的模擬粘度的變化情況，圖中圖像表明：模擬的初始階段水的模擬粘度值波動較大，隨著模擬時間步的推移，水的模擬粘度值上、下波動越來越緩，逐步趨向穩(wěn)定。

表1給出不同的溫度、壓力下水的粘度模擬值及其與化學(xué)化工物性數(shù)據(jù)手冊［11］中的試驗值的比較分析。

圖3 水的粘度Fig.3Viscosity of water

表1 不同溫度、壓強(qiáng)下系統(tǒng)的剪切粘度Table 1Viscosity under different temperatures and pressure

綜合觀察圖3和表1可得，當(dāng)模擬系統(tǒng)的溫度保持不變時，壓強(qiáng)增大則水的模擬粘度值增大；當(dāng)模擬系統(tǒng)的壓強(qiáng)保持不變時，溫度升高則水的模擬粘度值減小。這與化學(xué)化工物性數(shù)據(jù)手冊中的粘度的變化規(guī)律一致。

3 結(jié)論

基于分子動力學(xué)模擬理論，針對水分子特殊的極性分子結(jié)構(gòu)，引入Ewald加和法計算長程作用力，為避免水分子平移和旋轉(zhuǎn)運動相關(guān)的大量三角函數(shù)的計算，引入四元素法。模擬通過溫度壓力調(diào)整達(dá)到等溫等壓，最終根據(jù)粘度的模擬計算公式，由剪切應(yīng)力的統(tǒng)計數(shù)據(jù)求得在高溫高壓下水的粘度值。取7× 7×7的模擬系統(tǒng)，系統(tǒng)溫度、壓強(qiáng)分別為：473.16、573.16 K，10、20 MPa時，水的模擬粘度值與實驗粘度數(shù)據(jù)的相對誤差分別在1.88%～6.97%之間，且水的模擬粘度隨著溫度的升高而減小，隨著壓強(qiáng)的增大而增大，與實驗所得粘度的變化規(guī)律一致。

［1］徐堯，劉躍龍，劉夠生.水分子在白云母表面吸附的分子動力學(xué)模擬［J］.化工學(xué)報，2014（12）：4814-4822.

［2］Song H L.The Rheology of Liquid Pentane Isomers by Non-Equilibrium Molecular Dynamics Simulations［J］. Molecular Simulation，2000，23（4）：243-256.

［3］李健，楊延清，羅賢，等.分子動力學(xué)模擬在復(fù)合材料界面研究中的進(jìn)展［J］.稀有金屬材料與工程，2013，42（3）：644-648.

［4］Wei L，Yang C L，Zhu Y T，et al.Interactions Between Single-Walled Carbon Nanotubes and Polyethylene/Polypropylene/Polystyrene/Poly（phenylacetylene）/Poly（p-phenylenevinylene）Considering Repeat UnitArrangements and Conformations：A Molecular Dynamics Simulation Study［J］.The Journal of Chemical Physics，2008，112：1803-1811.

［5］Susumu Kotake.Molecular Simulation Engineering［J］. JSME International Journal，1995（1）：38-48.

［6］黃富才，呂艷東，郭曉紅，等.結(jié)實期水分供應(yīng)對寒地水稻灌漿動態(tài)和產(chǎn)量的影響［J］.黑龍江八一農(nóng)墾大學(xué)學(xué)報，2013，25（1）：12-18.

［7］馬靜，董秀芹，張敏華.超（近）臨界水溶劑特性的分子動力學(xué)模擬［J］.計算機(jī)與應(yīng)用化學(xué)，2007，24（7）：872-874.

［8］孫煒，黃素逸，王存文，等.超臨界水密度和自擴(kuò)散系數(shù)預(yù)測的分子動力學(xué)模擬［J］.華中科技大學(xué)學(xué)報：自然科學(xué)版，2008，36（5）：103-105.［9］Chen C K，Lin D T.TIP4P Potential For Lid-Driven Cavity Flow［J］.Acta Mechanica，2004（178）：223-237.

［10］Swope W C，Anderson H C.Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules：Application to Small Water Clusters［J］.The Journal of Chemical Physics，1982（76）：637-649.

［11］劉光啟.化學(xué)化工物性數(shù)據(jù)手冊（無機(jī)卷）［M］.北京：化學(xué)工業(yè)出版社，2002.

Water Viscosity of Molecular Dynamics Calculations under Different Temperatures and Pressure

Zhang Junhui1，Zhang Xin2，Lan Meihua3
（1.College of Engineering，Heilongjiang Bayi Agricultural University，Daqing 163319；2.Daqing Oilfield Drilling Engineering and Technology Institute；3.No.1 Logging Company，Daqing Drilling and Exploration Engineering Corporation）

It is difficult to achieve in the lab under the extreme conditions of high temperature and high pressure，therefore，it is an important method to calculate the water viscosity in high temperature and high pressure by using molecular dynamics simulation. Based on molecular dynamics simulation and special structure of water molecules，this paper introduced Ewald method for long-range force and four-element method.Finally，the viscosity under different temperature and pressure was simulated in the isothermal isobaric system，and the results were compared with the experimental value.

molecular dynamics；long-range force；isothermal pressure control；four-element method；shear stress；viscosity

O313

A

1002-2090（2016）04-0089-04

10.3969/j.issn.1002-2090.2016.04.020

2015-06-06

張鈞惠（1981-），女，助教，燕山大學(xué)畢業(yè)，現(xiàn)主要從事力學(xué)方向的教學(xué)與研究工作。