張鈞惠，張欣，蘭美華

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

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

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

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

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

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

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

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

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

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

1.1 初始條件、邊界條件

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

1.2 勢(shì)函數(shù)

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

TIP4P型勢(shì)函數(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 分子間相互作用力

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

式中，第一項(xiàng)為非鍵結(jié)力，在模擬過(guò)程中只需考慮截?cái)喟霃絻?nèi)的非鍵結(jié)力即可；第二項(xiàng)為靜電力，由于靜電力隨半徑收斂速度較慢，因此引入長(zhǎng)程作用力，采用Ewald加和法來(lái)計(jì)算分子間的長(zhǎng)程作用力。

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

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

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

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

在每個(gè)時(shí)間步，速度乘以速度標(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方法對(duì)系統(tǒng)壓力進(jìn)行調(diào)整。在每個(gè)時(shí)間步，將水分子質(zhì)心坐標(biāo)乘以一個(gè)系數(shù)調(diào)整系統(tǒng)體積，進(jìn)而達(dá)到調(diào)整系統(tǒng)壓強(qiáng)的目的。

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

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

1.5 四元素法

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

這四個(gè)變量q1、q2、q3、q4滿(mǎn)足

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

1.6 粘度計(jì)算

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

式中，ηs為粘度為系統(tǒng)的應(yīng)力分量，當(dāng)α=β時(shí)，Pαβ為正應(yīng)力；當(dāng)α≠β時(shí)，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)過(guò)離散計(jì)算得到模擬計(jì)算的粘度公式：

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

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

選取7×7×7的分子動(dòng)力學(xué)模型，系統(tǒng)中含有1 688個(gè)水分子，取截?cái)喟霃饺?.5σoo。實(shí)例分別模擬了溫度為573.16 K、473.16 K，壓強(qiáng)為10 MPa、20 MPa條件下水的粘度。模擬先經(jīng)過(guò)10萬(wàn)步的弛豫過(guò)程，使分子動(dòng)力學(xué)模擬系統(tǒng)達(dá)到平衡，再循環(huán)30萬(wàn)步統(tǒng)計(jì)平均正應(yīng)力和平均剪切應(yīng)力數(shù)據(jù)，并根據(jù)應(yīng)力與粘度的關(guān)系計(jì)算水在不同溫度和壓強(qiáng)下的粘度。在對(duì)水的應(yīng)力計(jì)算時(shí)，只統(tǒng)計(jì)距離各邊界0.5α0以?xún)?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)力在模擬的初始階段波動(dòng)較大，隨著模擬時(shí)間步的推移，波動(dòng)逐漸變小，最終，平均正應(yīng)力逐步趨向期望的壓強(qiáng)值10 MPa和20 MPa，平均剪切應(yīng)力在0 MPa附近上下波動(dòng)。

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

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

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

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

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

3 結(jié)論

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

［1］徐堯，劉躍龍，劉夠生.水分子在白云母表面吸附的分子動(dòng)力學(xué)模擬［J］.化工學(xué)報(bào)，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］李健，楊延清，羅賢，等.分子動(dòng)力學(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é)實(shí)期水分供應(yīng)對(duì)寒地水稻灌漿動(dòng)態(tài)和產(chǎn)量的影響［J］.黑龍江八一農(nóng)墾大學(xué)學(xué)報(bào)，2013，25（1）：12-18.

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

［8］孫煒，黃素逸，王存文，等.超臨界水密度和自擴(kuò)散系數(shù)預(yù)測(cè)的分子動(dòng)力學(xué)模擬［J］.華中科技大學(xué)學(xué)報(bào)：自然科學(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ù)手冊(cè)（無(wú)機(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é)與研究工作。