超臨界二氧化碳擴散行為的分子動力學模擬

，， ，，

(1.大連理工大學 化工學院，遼寧 大連 116024；2.大連理工大學 化機學院，遼寧 大連 116024)

?開發與研究?

超臨界二氧化碳擴散行為的分子動力學模擬

王寶和1，楊劉瑋1，李群1，夏良志2，王剛1

(1.大連理工大學 化工學院，遼寧 大連 116024；2.大連理工大學 化機學院，遼寧 大連 116024)

二氧化碳分子選用EPM2模型，乙醇分子選用TraPPE-UA模型，采用分子動力學模擬技術，探討了超臨界二氧化碳的自擴散系數及乙醇—超臨界二氧化碳體系的無限稀互擴散系數。模擬結果表明，超臨界二氧化碳的自擴散系數和乙醇—超臨界二氧化碳體系的無限稀互擴散系數均隨著溫度的升高或壓力的降低而增大。

分子動力學；超臨界二氧化碳；擴散系數；模擬

超臨界二氧化碳的擴散系數與其許多物性參數相互關聯，是科學研究、工程開發和過程設計不可缺少的基礎數據。鑒于超臨界二氧化碳分子間相互作用的復雜性，完全從理論上進行準確預測仍存在較大的難度；苛刻的測試條件使得超臨界二氧化碳擴散系數的實驗測定十分困難。隨著計算機和分子動力學(MD)模擬技術的發展，從分子水平研究超臨界二氧化碳的擴散規律，已經引起國內外許多學者的極大關注[1-3]。俞聯夢等[4]采用MD模擬方法，探討了CO2-水體系的PVT性質、擴散系數、微觀結構、配位數及氫鍵作用。Higashi等[5]利用MD模擬了純CO2流體和CO2-萘混合體系的CO2擴散系數。Colina等[6]采用Monte Carlo模擬技術，計算了超臨界CO2的密度、焦耳—湯普森系數、等壓熱容、等溫壓縮系數等物性參數。Zhou等[7]利用MD模擬方法，探究了超臨界二氧化碳中有機物的擴散行為，并提出一個新的計算二氧化碳無限稀互擴散系數方程。本文擬采用MD模擬技術，研究超臨界二氧化碳的自擴散系數及乙醇—超臨界二氧化碳體系的無限稀互擴散系數。

1 模擬方法

1.1二氧化碳及乙醇的勢能模型

二氧化碳分子模型有很多，如MSM、EPM2、TraPPE等。本研究的二氧化碳選用EPM2模型，其分子結構和勢能參數如圖1和表1所示。

圖1 二氧化碳分子的結構示意圖

其中σC-C為C原子之間L-J勢能的尺度參數，σO-O為O原子之間L-J勢能的尺度參數，εC為C原子之間L-J勢能的能量參數，εO為O原子之間L-J勢能的能量參數，kB為Boltzmann常數(kB=1.380 6×10-23J/K )，rC-O為C原子與O原子之間的鍵長，θ為鍵角，qO和qC分別為二氧化碳分子中O原子和C原子所帶電荷，e為基本電荷(1e=1.6×10-19C)。

表1 二氧化碳分子的結構參數[8]

乙醇分子模型有全原子模型和聯合原子模型等。由于全原子模型的計算耗時較長，一般采用聯合原子模型。乙醇分子的聯合原子模型有OPLS-UA和TraPPE-UA模型，本研究的乙醇分子選用TraPPE-UA模型，其分子結構及其參數見圖2和表2。該模型有四個作用位點，分別為甲基CH3(1號C原子)、亞甲基CH2(2號C原子)、O原子和H原子，各個作用位點的相關勢能參數見表3。

表2 乙醇分子的結構參數[9]

表3 TraPPE-UA模型參數

圖2 乙醇分子的結構示意圖

其中rC1-C2為1號C原子與2號C原子之間的鍵長；rC2-O為2號C原子與O原子之間的鍵長；rO-H為O原子與H原子之間的鍵長；θC1-C2-O為1號C原子、2號C原子及O原子組成的鍵角；θC2-O-H為2號C原子、O原子及H原子組成的鍵角；θC1-C2-O-H為1號C原子、2號C原子、O原子及H原子之間形成的二面角。

在進行超臨界二氧化碳的自擴散系數及乙醇—超臨界二氧化碳體系的無限稀互擴散系數的MD模擬過程中，每個分子(包括二氧化碳分子和乙醇分子)被視為剛性結構，忽略分子內的作用勢能，只考慮分子間作用勢能，分子間的勢能函數如式(1)～(3)所示[10]。

U=Ue+ULJ

(1)

(2)

(3)

兩種不同粒子間的L-J勢能參數，用Lorentz-Berthelot規則計算，如式(4)所示。

(4)

對于短程L-J勢能函數，為了避免計算過程中由于邊界位置的不連續，導致Hamiltonian函數的不守恒，從而引起誤差，因此采用勢能截斷的方法來處理此問題，如式(5)所示[11]。

(5)

采用作用場法來處理長程靜電勢能。將靜電力分為兩部分，位于截斷半徑內的計算所有電荷間的庫倫靜電力，將截斷半徑外的庫倫作用看作介電常數為εs的均勻介質對電荷的作用。依照作用場法，長程靜電勢能可表示為式(6)[11]。

Ue=

(6)

通常，環境介電常數εs取無窮大，因此，式(6)可簡化為式(7)

Ue=

(7)

1.2模擬體系的建立

1.2.1 超臨界二氧化碳自擴散的模擬體系

將超臨界二氧化碳分子的初始結構按照面心立方(FCC)的形式均勻排布于邊長為L的立方體盒子里，L的大小由一定的溫度和壓力下超臨界二氧化碳的實驗密度來確定。且體系的質心位于盒子中心處。如圖3所示。

圖3 超臨界二氧化碳自擴散的模擬體系示意圖

1.2.2 乙醇—超臨界二氧化碳無限稀互擴散的模擬體系

對于溶質乙醇在溶劑超臨界二氧化碳中的二元混合體系，無限稀意味著溶質乙醇分子四周幾乎都充滿了溶劑超臨界二氧化碳分子。在工程實際計算中，當溶質濃度達到5%，甚至10%時，仍可按照無限稀狀態進行計算。在本文的MD模擬中，為了實現無限稀狀態，模擬體系的總分子數為256，其中乙醇分子數為1，超臨界二氧化碳分子數為255，初始結構為面心立方(FCC)結構，如圖4所示。由于整個體系為無限稀狀態，我們假設乙醇分子的引入不影響整個體系的密度及體積，因此，立方體的邊長L仍由一定溫度T和壓力p下超臨界二氧化碳的實驗密度來確定。

圖4 乙醇—二氧化碳無限稀互擴散的模擬體系示意圖

1.3模擬細節

本模擬在x、y、z三個方向上均采用周期性邊界條件；模擬時間步長為1 fs，總模擬時間為1 000 ps，前500 ps使得系統達到平衡，后500 ps統計計算并輸出系統的均方位移MSD。選取正則系綜(NVT)，并采用Woodcock控溫法維持體系溫度的恒定；依照設定的溫度，隨機分布分子的初始平動速度；為了保證系統的總動量為0，每隔1萬步矯正體系的質心，使之始終處于盒子的中心處；本文模擬數據均采用LAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)軟件計算得到。

2 結果與討論

2.1超臨界二氧化碳的自擴散系數

2.1.1 超臨界二氧化碳自擴散系數的確定方法

采用Einstein法，計算超臨界二氧化碳自擴散系數的方程，如式(8)所示[12]。

(8)

根據MD模擬的計算結果，利用方程(9)可計算出均方位移(MSD)。

(9)

根據方程(8)和(9)，可得式(10)。

(10)

由式(10)可見，當模擬時間t很長時，MSD與時間t的關系曲線應該為一條直線，根據其斜率可求得自擴散系數。

在二氧化碳分子數N=256，溫度T=323 K，壓力p=10.99 MPa的條件下，MD模擬得到的均方位移(MSD)隨時間t的變化及線性擬合曲線，如圖5所示。可見，均方位移(MSD)與時間t具有很好的線性關系，從而得到超臨界二氧化碳的自擴散系數Dself=5.634×10-8m2/s。

圖5 均方位移與時間關系曲線

2.1.2 壓力對二氧化碳自擴散系數的影響

在二氧化碳分子數N=256，T=323 K的條件下，MD模擬得到的二氧化碳自擴散系數與壓力之間的關系曲線，并與文獻值[13]進行了比較，如圖6所示。由圖6可見，二氧化碳的自擴散系數隨著壓力的升高而減小；在低密度的亞臨界區，流體很容易被壓縮，擴散系數隨著體系壓力的升高而快速減小；在靠近臨界點附近，二氧化碳分子的擴散系數對體系壓力的變化也很敏感；當體系壓力大于15 MPa時，隨著壓力的升高，二氧化碳的自擴散系數變化緩慢，呈現出液體的特性。

圖6 不同壓力下二氧化碳自擴散系數的MD模擬結果

2.1.3 溫度對二氧化碳自擴散系數的影響

在二氧化碳分子數N=256，壓力p=8 MPa的條件下，MD模擬得到的二氧化碳自擴散系數與溫度之間的關系曲線，并與文獻值[14]進行了比較，如圖7所示。由圖7可知，隨著溫度的升高，二氧化碳自擴散系數增大，并在不同的溫度段展現出不同的增大速度；當溫度較低時，即T<300 K，模擬系統為液相，隨著溫度的升高，自擴散系數增大比較緩慢，當T>300 K時，隨著溫度的升高，自擴散系數顯著增大，并在臨界點附近有明顯的轉折點；可見，在超臨界狀態下，二氧化碳分子具有較強的擴散能力。

圖7 不同溫度下二氧化碳自擴散系數的MD模擬結果

2.2乙醇—超臨界二氧化碳體系的無限稀互擴散系數

2.2.1 乙醇—超臨界二氧化碳無限稀體系的互擴散系數確定方法

采用Einstein法，確定乙醇—超臨界二氧化碳無限稀體系的互擴散系數，如式(11)所示[12]。

(11)

與自擴散系數確定方法類似，根據MD模擬的計算結果，可計算出均方位移(MSD)，再由MSD與時間t關系曲線的斜率可求得無限稀互擴散系數。

在溫度T=323 K，壓力p=8.41 MPa，二氧化碳分子數為255，乙醇分子數為1的條件下，MD模擬得到的均方位移(MSD)隨時間t的變化及線性擬合曲線，如圖8所示，從而得到乙醇—超臨界二氧化碳無限稀體系的互擴散系數D12=6.048×10-8m2/s。

圖8 均方位移與時間的擬合曲線

2.2.2 壓力對該體系無限稀互擴散系數的影響

在T=323 K，二氧化碳分子數為255，乙醇分子數為1的條件下，MD模擬得到的乙醇—超臨界二氧化碳無限稀體系的互擴散系數與壓力之間的關系曲線，并與文獻值[14]進行了比較，如圖9所示。與二氧化碳自擴散系數的變化規律相似，隨著壓力的增大，互擴散系數不斷減小，擴散能力變弱；在低壓區，互擴散系數隨壓力的變化很大，在高壓區，壓力的影響變得很微弱，擴散系數的降低速率變得很小。

圖9 不同壓力下互擴散系數的MD模擬結果

2.2.3 溫度對該體系無限稀互擴散系數的影響

在二氧化碳分子數為255，乙醇分子數為1，壓力p=8 MPa，MD模擬得到的乙醇—超臨界二氧化碳無限稀體系的互擴散系數與溫度之間的關系曲線，如圖10所示。可見，隨著溫度的升高，互擴散系數增大，擴散能力增強；當溫度小于臨界溫度時，互擴散系數隨溫度的增大而緩慢地增加，當溫度位于臨界溫度附近時，互擴散系數的增加速率突然增大，曲線有一個明顯的轉折點。

圖10 不同溫度下互擴散系數的MD模擬結果

3 結論

采用分子動力學模擬方法，考察了溫度和壓力對超臨界二氧化碳自擴散系數和乙醇—二氧化碳無限稀體系互擴散系數的影響。模擬結果表明，超臨界二氧化碳的自擴散系數和乙醇—二氧化碳無限稀體系互擴散系數，均隨著溫度的升高或壓力的降低而增大。

[1]張新波.水和二氧化碳體系的分子模擬研究[D].廈門:廈門大學,2007.

[2]Chatzis G,Samios J.Binary mixtures of supercritical carbon dioxide with methanol. a molecular dynamics simulation study[J].Chem Phys Lett,2003,374(1-2):187-193.

[3]Stubbs J M,Siepmann J I.Binary phase behavior and aggregation of dilute methanol in supercritical carbon dioxide:a monte carlo simulation study[J].J Chem Phys,2004,121(3):1525-1534.

[4]俞聯夢,劉錦超.二氧化碳—水二元混合流體的分子動力學模擬研究[J].西南民族大學學報：自然科學版,2009,35(6):1248-1254.

[5]Higashi H,Iwai Y,Arai Y.Calculation of self-diffusion and tracer diffusion coefficients near the critical point of carbon dioxide using molecular dynamics simulation[J].Industrial & Engineering Chemistry Research,2000,39(12):4567-4570.

[6]Colina C M,Olivera-Fuentes C G,Siperstein F R,et al.Thermal properties of supercritical carbon dioxide by Monte Carlo simulations[J].Molecular Simulation,2003,29(6-7):405-412.

[7]Zhou J,Lu X,Wang Y,et al.Molecular dynamics investigation on the infinite dilute diffusion coefficients of organic compounds in supercritical carbon dioxide[J].Fluid Phase Equilibria,2000,172(2):279-291.

[8]張 陽.分子模擬在純超臨界流體及其二元混合物體系中的應用[D].北京:清華大學,2005.

[9]李 勇,劉錦超,蘆鵬飛,等.從常溫常壓到超臨界乙醇的分子動力學模擬[J].物理學報,2010(7):4880-4887.

[10]Rapaport D C.The art of molecular dynamics simulation[M].Cambridge:Cambridge University Press,1995.

[11]陳正隆,徐為人,湯立達.分子模擬的理論與實踐[M].北京:化學工業出版社,2007.

[12]Allen M P,Tildesley D J.Computer simulation of liquids[M].Oxford:Oxford University Press,1989.

[13]Coelho L A F,Oliveira J V,Tavares F W.Dense fluid self-diffusion coefficient calculations using perturbation theory and molecular dynamics[J].Brazilian Journal of Chemical Engineering,1999,16(3):319-329.

[14]何 巖,張敏華,姜浩錫.純CO2體系擴散性質的分子動力學模擬[J].化學工業與工程,2009,26(1):57-61.

MolecularDynamicsSimulationofDiffusionBehaviorforSupercriticalCarbonDioxide

WANGBao-he1,YANGLiu-wei1,LIQun1,XIALiang-zhi2,WANGGang1

(1.School of Chemical Engineering,Dalian University of Technology,Dalian 116024,China;2.School of Chemical Machinery, Dalian University of Technology,Dalian 116024,China)

Using carbon dioxide EPM2 potential energy model and ethanol TraPPE-UA potential energy model, self-diffusion coefficients of supercritical carbon dioxide and infinite dilute inter-diffusion coefficients of ethanol in supercritical carbon dioxide are investigated by molecular dynamics simulation. The simulation results indicate that the self- diffusion coefficient of supercritical carbon dioxide increases with the temperature increasing or pressure decreasing,as well as the inter-diffusion coefficient of ethanol in supercritical carbon dioxide system.

molecular dynamics;supercritical carbon dioxide;diffusion coefficient;simulation

2014-06-18

王寶和(1959-)，男，副教授，主要從事不同形貌微納結構的制備、干燥及分子動力學模擬研究，電話：(0411)84986167，E-mail：wbaohe@163.com。

O561

A

1003-3467(2014)08-0031-05