周曉平，周 偉，郝江濤

（鄭州大學(xué)西亞斯國際學(xué)院，河南 鄭州 451150）

0 引 言

分子動力學(xué)模擬是依靠牛頓力學(xué)來模擬分子體系運(yùn)動的一套分子模擬方法，不僅可以用于模擬物質(zhì)的宏觀性質(zhì)，還可以得出微觀結(jié)構(gòu)、粒子運(yùn)動等信息[1]。一般認(rèn)為壓強(qiáng)超過100MPa就是超高壓，物質(zhì)在高壓作用下，其物理和化學(xué)性質(zhì)會發(fā)生巨大的變化，在超高壓作用下，液體的壓縮性顯著增大[2-3]。高壓物理實(shí)驗(yàn)面臨著實(shí)驗(yàn)條件、費(fèi)用、實(shí)驗(yàn)設(shè)備的限制，而分子動力學(xué)計算利用計算機(jī)模擬，可以很好地解決這些問題。乙醇分子是常見的有機(jī)、極性分子，常用于化工生產(chǎn)中，研究其熱力學(xué)性質(zhì)和動力學(xué)性質(zhì)對工業(yè)生產(chǎn)有重要意義。且乙醇和水在很多性質(zhì)上都存在相似性，被認(rèn)為是水的理想替代物[4]，關(guān)于乙醇性質(zhì)的研究也一直是熱點(diǎn)。張陽[5]用Monte Carlo方法研究了超臨界乙醇的結(jié)構(gòu)性質(zhì)和氫鍵，壓強(qiáng)范圍從常壓到高壓70 MPa；Saiz等[6]用分子動力學(xué)方法模擬研究了不同熱力學(xué)狀態(tài)下液態(tài)乙醇的結(jié)構(gòu)、氫鍵和擴(kuò)散等動態(tài)性質(zhì)；Markus等[7]研究了溫度和壓強(qiáng)（0.1～35 MPa）對甲醇和乙醇體系中氫鍵的影響。目前對乙醇的研究壓強(qiáng)一般都在100MPa以內(nèi)，更高壓強(qiáng)下乙醇性質(zhì)的文章還少見報道。因此本文對液態(tài)乙醇的研究進(jìn)行了補(bǔ)充，采用平衡分子動力學(xué)方法，模擬溫度298～500K、壓強(qiáng)從10MPa到超高壓300 MPa條件下不同狀態(tài)點(diǎn)乙醇的熱力學(xué)性質(zhì)、結(jié)構(gòu)性質(zhì)和動力學(xué)性質(zhì)。本文的乙醇分子動力學(xué)模擬值可為其傳質(zhì)性質(zhì)應(yīng)用提供理論參考。

1 模型建立

1.1 勢能模型

勢能模型的選取對模擬結(jié)果的準(zhǔn)確與否有著決定性影響。目前文獻(xiàn)中分子模擬常用的乙醇模型主要是全原子模型和聯(lián)合原子模型。鑒于全原子模型所需的計算時間較長，因此大部分文獻(xiàn)采用聯(lián)合原子模型進(jìn)行模擬。聯(lián)合原子模型有OPLS-UA（optimized potentials for liquid simulations-united atom）模型和TraPPE-UA（transferable potentials for phase equilibria-united atom）模型兩種。OPLS-UA模型由Jorgensen[8]開發(fā)，主要針對液相有機(jī)小分子，能夠很好地模擬液態(tài)乙醇的各種動態(tài)性質(zhì)。但是，在Chalaris和Samios[9-10]用OPLS-UA模型對超臨界甲醇進(jìn)行分子動力學(xué)模擬的研究中，壓強(qiáng)模擬值的準(zhǔn)確性很差，說明OPLS-UA勢能模型在壓強(qiáng)較高時并不能真實(shí)反映分子的動態(tài)性質(zhì)。TraPPE-UA模型是由Siepmann[11]小組開發(fā)，由于其模型參數(shù)已經(jīng)用實(shí)驗(yàn)的相平衡數(shù)據(jù)優(yōu)化過，因此在很大的壓強(qiáng)范圍內(nèi)都有較好的準(zhǔn)確度。本文采用TraPPE-UA勢能模型對高壓下乙醇分子團(tuán)簇的性質(zhì)進(jìn)行研究。

在TraPPE-UA模型中，乙醇分子由CH3、CH2基團(tuán)和O原子、H原子合計4個作用位點(diǎn)構(gòu)成。乙醇分子間的勢能計算式為

式中：qi、qj——第i、j個原子或基團(tuán)上所帶電荷；

rij——原子或離子i與j間的距離；

σij、εij——原子間L-J作用參數(shù)。

式（1）右側(cè)第 1 項(xiàng)表示短程 Lennard-Jones（L-J）勢，第2項(xiàng)表示長程庫侖勢。TraPPE-UA模型和OPLS-UA模型的相關(guān)參數(shù)見表1。

表1 TraPPE-UA模型和OPLS-UA模型參數(shù)列表

1.2 實(shí)驗(yàn)方法

本研究采用Materials Explorer軟件，在溫度為298，400，500 K，壓強(qiáng)分別為 10，100，300 MPa 條件下對乙醇分子體系進(jìn)行了分子動力學(xué)模擬。采用Nose-Hoover熱浴法控制溫度，完成不同壓強(qiáng)條件下的等溫模擬[13]。分子動力學(xué)模擬需進(jìn)行以下假設(shè)：1）兩體有效勢能近似；2）最小鏡像準(zhǔn)則[14]。分子動力學(xué)模擬的起始構(gòu)型為面心立方晶格，將乙醇分子的起始取向定為隨機(jī)分布，采用立方周期性邊界條件。模擬過程中所用的乙醇分子總數(shù)為100個。采用球形截去法進(jìn)行位能截斷[14]，采用Ewald加和法求解長程靜電相互作用力[15]，采用Shake算法固定乙醇分子的幾何構(gòu)型[16]。分子動力學(xué)模擬過程中，每個分子都可視為是剛性體，運(yùn)用五階Gear預(yù)測-校正法求解體系的運(yùn)動方程[17]。分子動力學(xué)模擬過程中，時間步長取10-15s，單次模擬總時間取3×10-10s，采用NVT系綜進(jìn)行10-10s的初始化，使得體系達(dá)到相對平衡，隨后在NPT系綜中完成2×10-10s的模擬，完成各種參數(shù)的計算，模擬過程中每隔10-14s就輸出一次構(gòu)型，用于結(jié)果分析。

2 結(jié)果與分析

2.1 熱力學(xué)性質(zhì)

焓是物質(zhì)重要的熱力學(xué)狀態(tài)量，體系焓計算式為

式中：E——體系的內(nèi)能，kJ；

P——體系的壓強(qiáng)，MPa；

V——體系的體積，m3。

表2是本文模擬溫度為298，400，500K，壓強(qiáng)為10，100，300 MPa 時乙醇的焓值，并與文獻(xiàn)[18]的模擬值和文獻(xiàn)[19]用狀態(tài)方程推導(dǎo)的焓值進(jìn)行比較，結(jié)果較為接近，可見本文模擬結(jié)果可靠。

表2 乙醇在不同狀態(tài)點(diǎn)的焓值

從表中可以看出，當(dāng)壓強(qiáng)相同時，隨著溫度的升高乙醇體系的焓值逐漸增大。這是由于當(dāng)溫度升高，乙醇分子間的熱運(yùn)動加劇，體系內(nèi)能增大，焓也跟著增大。相同溫度下，本文模擬的焓值隨壓強(qiáng)的增大而增大，但文獻(xiàn)[18-19]中焓隨壓強(qiáng)的變化關(guān)系不明確，由于現(xiàn)有文獻(xiàn)中尚未見到高壓下乙醇焓值實(shí)驗(yàn)數(shù)據(jù)的相關(guān)報道，此計算結(jié)果有待進(jìn)一步驗(yàn)證。

2.2 結(jié)構(gòu)性質(zhì)

徑向分布函數(shù)（RDF)是反映物質(zhì)微觀結(jié)構(gòu)的重要物理量，主要用于描述距離中心粒子為r處出現(xiàn)另一粒子的概率密度與隨機(jī)分布概率密度的比值[13]，其計算式為

式中：ρ——系統(tǒng)密度；

dN——距離中心分子（r，r+dr）區(qū)間內(nèi)的分子數(shù)目。

圖1為T=298 K時乙醇分子的氧氫徑向分布函數(shù)隨壓強(qiáng)的變化。圖2為P=10MPa時氧氫徑向分布函數(shù)隨溫度的變化。

從圖1可看出，壓強(qiáng)變化對OH徑向分布函數(shù)第1峰影響較大，隨著壓強(qiáng)的增大，OH徑向分布函數(shù)第1峰分別出現(xiàn)在1.69，1.67，1.65?，峰位逐漸左移，峰值依次減小，第2峰也下降，但不如第1峰明顯，在Petravic[20]和Dellis[21]等的研究中，也發(fā)現(xiàn)了OH徑向分布函數(shù)的第1峰隨壓強(qiáng)增大而減小的現(xiàn)象。計算結(jié)果表明，隨著壓強(qiáng)的增加，第一配位圈內(nèi)乙醇分子間的氫原子與氧原子間距逐漸減小，氫鍵作用逐漸增強(qiáng)。

圖1 298K時乙醇的氧氫徑向分布函數(shù)

圖2 10MPa時乙醇的氧氫徑向分布函數(shù)

從圖2可看出，當(dāng)壓強(qiáng)相同時，隨著溫度的升高，OH徑向分布函數(shù)第1峰和第2峰逐漸下降，峰谷抬高，峰位右移。計算結(jié)果表明隨著溫度的升高，乙醇分子間的氫原子與氧原子間距逐漸增大，氫鍵作用逐漸減弱。表3為各狀態(tài)點(diǎn)乙醇的OH、OO、HH徑向分布函數(shù)第1特征峰值與峰位。隨著溫度的升高和壓強(qiáng)的增大，OH、OO、HH徑向分布函數(shù)第1峰峰值均下降，但OH徑向分布函數(shù)的變化較OO、HH徑向分布函數(shù)更為明顯。

表3 不同狀態(tài)下乙醇的OH、OO、HH徑向分布函數(shù)第1特征峰值與峰位

圖3、圖4為100MPa時OO、HH徑向分布函數(shù)隨溫度的變化。從圖中可看出，相同壓強(qiáng)下，隨著溫度的升高，OO、HH徑向分布函數(shù)第1峰下降，峰位右移，峰谷抬高，第2峰也下降，且在500K時，第2峰消失。這說明隨著溫度的升高，乙醇的長程有序程度逐漸下降。這與高壓下水的結(jié)構(gòu)隨溫度和壓強(qiáng)變化的趨勢相同[22]。

圖3 100MPa時乙醇的氧氧徑向分布函數(shù)

圖4 100MPa時乙醇的氫氫徑向分布函數(shù)

圖5～圖8為乙醇的MeH、MeO、MeMe徑向分布函數(shù)。隨著溫度的升高和壓強(qiáng)的增大，各徑向分布函數(shù)第1峰變化最明顯，峰高下降，峰谷抬高。500 K時，MeH徑向分布函數(shù)的第2、3、4峰峰高已超過第1峰。MeMe徑向分布函數(shù)第2峰逐漸消失。

通過對比不同原子和基團(tuán)間徑向分布函數(shù)的第一波谷位置rmin1可發(fā)現(xiàn)，rOH＜rHH＜rOO＜rMeH＜rMeO＜rMeMe。 計算結(jié)果表明：相同溫度和壓強(qiáng)下，第一配位圈內(nèi)親水的O原子與H原子距離最小，而疏水的烷基之間的距離最遠(yuǎn)。

2.3 動力學(xué)性質(zhì)

圖5 298K時乙醇的甲基氫徑向分布函數(shù)

圖6 100MPa時乙醇的甲基氧徑向分布函數(shù)

圖7 400K時乙醇的甲基甲基徑向分布函數(shù)

結(jié)合國內(nèi)外研究現(xiàn)狀，自擴(kuò)散系數(shù)D的計算方法主要包括Einstein法和Green-Kubo法。Eintein法通過對均方位移（MSD）求斜率得到擴(kuò)散系數(shù)。鑒于分子動力系統(tǒng)中的原子在永不停息地做無規(guī)則運(yùn)動，各原子在不同時刻的位置都不相同。將r→（t）、r→（0）作為粒子在t時刻和初始時刻的位置，則粒子均方位移計算式[13]為

圖8 300MPa時乙醇的甲基甲基徑向分布函數(shù)

式中＜＞表示平均值。

根據(jù)統(tǒng)計學(xué)原理，只要分子數(shù)目足夠多，計算時間足夠長，系統(tǒng)的任一瞬間都可以當(dāng)作時間的零點(diǎn)，所計算的平均值應(yīng)該相同。根據(jù)愛因斯坦的擴(kuò)散定律，擴(kuò)散系數(shù)（diffusion constant）與均方位移的關(guān)系是：

式中D即為粒子的擴(kuò)散系數(shù)。

當(dāng)模擬的時間很長時，均方位移對時間曲線的斜率的1/6就是擴(kuò)散系數(shù)。表4是本文用均方位移計算的乙醇的自擴(kuò)散系數(shù)。

表4 乙醇在不同狀態(tài)下的自擴(kuò)散系數(shù)

由表可知，在相同壓強(qiáng)下，隨溫度的升高，乙醇的自擴(kuò)散系數(shù)增大。這與文獻(xiàn)[23]測定的常壓下乙醇自擴(kuò)散系數(shù)隨溫度的變化趨勢相同，且隨著壓強(qiáng)的增大，這種變化趨勢逐漸減緩。而相同溫度下，隨著壓強(qiáng)的增加，乙醇自擴(kuò)散系數(shù)逐漸減小。這是因?yàn)閴簭?qiáng)增大使乙醇分子之間的距離變小，分子間排列更為緊密，相互碰撞現(xiàn)象明顯增加，從而阻礙了分子的擴(kuò)散運(yùn)動；當(dāng)溫度不斷上升時，乙醇分子間的距離逐漸增加，相互作用力減弱，此時將有利于分子的擴(kuò)散運(yùn)動。本文模擬值與文獻(xiàn)[18]的模擬結(jié)果基本一致，但在500K，10MPa時相差較大，由于實(shí)驗(yàn)條件的限制，高壓下乙醇分子團(tuán)簇的自擴(kuò)散系數(shù)實(shí)驗(yàn)數(shù)據(jù)匱乏，有待進(jìn)一步驗(yàn)證。

3 結(jié)束語

本文模擬了溫度298～500K，壓強(qiáng)從高壓10MPa到超高壓300MPa條件下乙醇分子體系的熱力學(xué)性質(zhì)、結(jié)構(gòu)性質(zhì)和動力學(xué)性質(zhì)，并得出以下結(jié)論：

1）相同壓強(qiáng)下，乙醇體系的焓值隨溫度的升高而增大，本文模擬值與文獻(xiàn)模擬值接近，本文所用TraPPE-UA模型可靠。

2）通過對徑向分布函數(shù)的研究，發(fā)現(xiàn)在壓強(qiáng)很大的情況下，乙醇體系的結(jié)構(gòu)依然很有規(guī)律。隨著壓強(qiáng)和溫度的升高，各徑向分布函數(shù)的第1峰高度明顯下降，乙醇的長程有序程度下降。

3）當(dāng)壓強(qiáng)相同時，乙醇的自擴(kuò)散系數(shù)隨溫度的升高而增大，但壓強(qiáng)增大時，這種變化趨勢逐漸減緩。乙醇的自擴(kuò)散系數(shù)隨溫度的升高而增大，而壓強(qiáng)對擴(kuò)散系數(shù)的影響與溫度正好相反。

[1]RAHMAN A.Correlations in the motion of atoms in liquid argon[J].Phys Rev，1964（136）：405-411.

[2]曾勇平，朱曉敏，楊正華.水、甲醇和乙醇液體微結(jié)構(gòu)性質(zhì)的Car-Parrinello分子動力學(xué)模擬[J].物理化學(xué)學(xué)報，2011，27（12）：2779-2785.

[3]曾勇平，時榮，楊正華.Be～（2+）在水、甲醇和乙醇中結(jié)構(gòu)性質(zhì)的從頭算分子動力學(xué)模擬[J].物理化學(xué)學(xué)報，2013，2（10）：2180-2186.

[4]羅輝，王睿，范維玉，等.分子動力學(xué)模擬計算超臨界CO的溶解度參數(shù)[J].石油學(xué)報（石油加工），2015，31（1）：78-83.

[5]張陽.分子模擬在純超臨界流體及其二元混合物體系中的應(yīng)用[D].北京：清華大學(xué)，2005.

[6]SAIZ L，PADRO J A，GUARDIA E.Structure and dyna mics of liquid ethanol[J].Phys Chem B，1997（101）：78-86.

[7]MARKUS M，HOFFMANN，MARK S C.Are there hydrogen bonds in supercritical methanol and ethanol[J].Phys Chem B，1998（102）：263-271.

[8]JORGENSEN W L.Optimized intermolecular potential functions for liquid alcohols[J].Phys Chem，1986（90）：1276-1284.

[9]CHALARIS M，SAMIOS J.Hydrogen bonding in supercritical methanol A molecular dynamics investigation[J].Phys Chem B，1999（103）：1161-1166.

[10]CHALARIS M，SAMIOS J.Translational and rotational dynamics in supercritical methanol from molecular dynamics simulation[J].Pure and Applied Chemistry，2004（76）：203-213.

[11]CHEN B，POTOFF J J，SIEPMANN J I.Monte carlo calculations for alcohols and their mixtures with alkanes.transferable potentials for phase equilibria.5.United-Atom description of primary secondary and tertiary alcohols[J].Phys Chem B，2001（105）：3093-3104.

[12]AND I Y G，PIOTROVSKAYA E M.Properties of coexisting phases for the ethanol-ethane binary system by computer simulation[J].Journal of Physical Chemistry B，1999，103（36）：7681-7686.

[13]陳正隆，徐為人，湯立達(dá).分子模擬的理論與實(shí)踐[M].北京：化學(xué)工業(yè)出版社，2007：6.

[14]ALLEN M P，TILDESLEY D J.Computer simulation of liquids[M].Oxford： Clareendon Press，1987：327-341.

[15]DE-LEEUW S W，PERRAM J W，SMITH E R.Simulation ofelectrostatic systemsin periodic boundary conditions I Lattice sums and dielectric constant[J].Proc R Soc Lond A，1980（373）：27-56.

[16]RYCHAERT J P，CICCOTTI G，BERENDSEN H J C.Numercial intergration of the cartesian equations of motion of a system with constraints:molecular dynamics of n-alkanes[J].J Chem Phys，1977（23）：327-341.

[17]GEAR C W.Numerical integration of ordinary differential equations[J].Mathemetics of Computation，1967，21（98）：146-156.

[18]李勇，劉錦超，蘆鵬飛，等.從常溫常壓到超臨界乙醇的分子動力學(xué)模擬[J].物理學(xué)報，2010，59（7）：4880-4887.

[19]DILLON H E，PENONCELLO S G.A fundamental equation for calculation of the thermodynamic properties ofethanol[J].International Journal of Thermophysics，2004，25（2）：321-335.

[20]PETRAVIC J，DELHOMMELLE J.Influence of temperature，pressure and internal degrees of freedom on hydrogen bonding and diffusion in liquid ethanol[J].Chemical Physics，2003，286（2/3）：303-314.

[21]DELLIS D，MICHALIS C A，SAMIOS J.Pressure and temperature dependence of the hydrogen bonding in supercriticalethanol： A computersimulation study[J].Journal of Physical Chemistry B，2005，109（39）：18575-18590.

[22]周曉平，楊向東，劉錦超.超高壓水的分子動力學(xué)模擬[J].高壓物理學(xué)報，2009，23（4）：1-5.

[23]KARGER N，VARDAG T.Temperature dependence of self iffusion in compressed monohydric alcohols[J].Journal of Chemical Physics，1990，93（93）：3437-3444.