超臨界二氧化碳擴(kuò)散行為的分子動(dòng)力學(xué)模擬

，， ，，

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

?開(kāi)發(fā)與研究?

超臨界二氧化碳擴(kuò)散行為的分子動(dòng)力學(xué)模擬

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

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

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

分子動(dòng)力學(xué)；超臨界二氧化碳；擴(kuò)散系數(shù)；模擬

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

1 模擬方法

1.1二氧化碳及乙醇的勢(shì)能模型

二氧化碳分子模型有很多，如MSM、EPM2、TraPPE等。本研究的二氧化碳選用EPM2模型，其分子結(jié)構(gòu)和勢(shì)能參數(shù)如圖1和表1所示。

圖1 二氧化碳分子的結(jié)構(gòu)示意圖

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

表1 二氧化碳分子的結(jié)構(gòu)參數(shù)[8]

乙醇分子模型有全原子模型和聯(lián)合原子模型等。由于全原子模型的計(jì)算耗時(shí)較長(zhǎng)，一般采用聯(lián)合原子模型。乙醇分子的聯(lián)合原子模型有OPLS-UA和TraPPE-UA模型，本研究的乙醇分子選用TraPPE-UA模型，其分子結(jié)構(gòu)及其參數(shù)見(jiàn)圖2和表2。該模型有四個(gè)作用位點(diǎn)，分別為甲基CH3(1號(hào)C原子)、亞甲基CH2(2號(hào)C原子)、O原子和H原子，各個(gè)作用位點(diǎn)的相關(guān)勢(shì)能參數(shù)見(jiàn)表3。

表2 乙醇分子的結(jié)構(gòu)參數(shù)[9]

表3 TraPPE-UA模型參數(shù)

圖2 乙醇分子的結(jié)構(gòu)示意圖

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

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

U=Ue+ULJ

(1)

(2)

(3)

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

(4)

對(duì)于短程L-J勢(shì)能函數(shù)，為了避免計(jì)算過(guò)程中由于邊界位置的不連續(xù)，導(dǎo)致Hamiltonian函數(shù)的不守恒，從而引起誤差，因此采用勢(shì)能截?cái)嗟姆椒▉?lái)處理此問(wèn)題，如式(5)所示[11]。

(5)

采用作用場(chǎng)法來(lái)處理長(zhǎng)程靜電勢(shì)能。將靜電力分為兩部分，位于截?cái)喟霃絻?nèi)的計(jì)算所有電荷間的庫(kù)倫靜電力，將截?cái)喟霃酵獾膸?kù)倫作用看作介電常數(shù)為εs的均勻介質(zhì)對(duì)電荷的作用。依照作用場(chǎng)法，長(zhǎng)程靜電勢(shì)能可表示為式(6)[11]。

Ue=

(6)

通常，環(huán)境介電常數(shù)εs取無(wú)窮大，因此，式(6)可簡(jiǎn)化為式(7)

Ue=

(7)

1.2模擬體系的建立

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

將超臨界二氧化碳分子的初始結(jié)構(gòu)按照面心立方(FCC)的形式均勻排布于邊長(zhǎng)為L(zhǎng)的立方體盒子里，L的大小由一定的溫度和壓力下超臨界二氧化碳的實(shí)驗(yàn)密度來(lái)確定。且體系的質(zhì)心位于盒子中心處。如圖3所示。

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

1.2.2 乙醇—超臨界二氧化碳無(wú)限稀互擴(kuò)散的模擬體系

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

圖4 乙醇—二氧化碳無(wú)限稀互擴(kuò)散的模擬體系示意圖

1.3模擬細(xì)節(jié)

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

2 結(jié)果與討論

2.1超臨界二氧化碳的自擴(kuò)散系數(shù)

2.1.1 超臨界二氧化碳自擴(kuò)散系數(shù)的確定方法

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

(8)

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

(9)

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

(10)

由式(10)可見(jiàn)，當(dāng)模擬時(shí)間t很長(zhǎng)時(shí)，MSD與時(shí)間t的關(guān)系曲線應(yīng)該為一條直線，根據(jù)其斜率可求得自擴(kuò)散系數(shù)。

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

圖5 均方位移與時(shí)間關(guān)系曲線

2.1.2 壓力對(duì)二氧化碳自擴(kuò)散系數(shù)的影響

在二氧化碳分子數(shù)N=256，T=323 K的條件下，MD模擬得到的二氧化碳自擴(kuò)散系數(shù)與壓力之間的關(guān)系曲線，并與文獻(xiàn)值[13]進(jìn)行了比較，如圖6所示。由圖6可見(jiàn)，二氧化碳的自擴(kuò)散系數(shù)隨著壓力的升高而減?。辉诘兔芏鹊膩喤R界區(qū)，流體很容易被壓縮，擴(kuò)散系數(shù)隨著體系壓力的升高而快速減?。辉诳拷R界點(diǎn)附近，二氧化碳分子的擴(kuò)散系數(shù)對(duì)體系壓力的變化也很敏感；當(dāng)體系壓力大于15 MPa時(shí)，隨著壓力的升高，二氧化碳的自擴(kuò)散系數(shù)變化緩慢，呈現(xiàn)出液體的特性。

圖6 不同壓力下二氧化碳自擴(kuò)散系數(shù)的MD模擬結(jié)果

2.1.3 溫度對(duì)二氧化碳自擴(kuò)散系數(shù)的影響

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

圖7 不同溫度下二氧化碳自擴(kuò)散系數(shù)的MD模擬結(jié)果

2.2乙醇—超臨界二氧化碳體系的無(wú)限稀互擴(kuò)散系數(shù)

2.2.1 乙醇—超臨界二氧化碳無(wú)限稀體系的互擴(kuò)散系數(shù)確定方法

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

(11)

與自擴(kuò)散系數(shù)確定方法類似，根據(jù)MD模擬的計(jì)算結(jié)果，可計(jì)算出均方位移(MSD)，再由MSD與時(shí)間t關(guān)系曲線的斜率可求得無(wú)限稀互擴(kuò)散系數(shù)。

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

圖8 均方位移與時(shí)間的擬合曲線

2.2.2 壓力對(duì)該體系無(wú)限稀互擴(kuò)散系數(shù)的影響

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

圖9 不同壓力下互擴(kuò)散系數(shù)的MD模擬結(jié)果

2.2.3 溫度對(duì)該體系無(wú)限稀互擴(kuò)散系數(shù)的影響

在二氧化碳分子數(shù)為255，乙醇分子數(shù)為1，壓力p=8 MPa，MD模擬得到的乙醇—超臨界二氧化碳無(wú)限稀體系的互擴(kuò)散系數(shù)與溫度之間的關(guān)系曲線，如圖10所示?？梢?jiàn)，隨著溫度的升高，互擴(kuò)散系數(shù)增大，擴(kuò)散能力增強(qiáng)；當(dāng)溫度小于臨界溫度時(shí)，互擴(kuò)散系數(shù)隨溫度的增大而緩慢地增加，當(dāng)溫度位于臨界溫度附近時(shí)，互擴(kuò)散系數(shù)的增加速率突然增大，曲線有一個(gè)明顯的轉(zhuǎn)折點(diǎn)。

圖10 不同溫度下互擴(kuò)散系數(shù)的MD模擬結(jié)果

3 結(jié)論

采用分子動(dòng)力學(xué)模擬方法，考察了溫度和壓力對(duì)超臨界二氧化碳自擴(kuò)散系數(shù)和乙醇—二氧化碳無(wú)限稀體系互擴(kuò)散系數(shù)的影響。模擬結(jié)果表明，超臨界二氧化碳的自擴(kuò)散系數(shù)和乙醇—二氧化碳無(wú)限稀體系互擴(kuò)散系數(shù)，均隨著溫度的升高或壓力的降低而增大。

[1]張新波.水和二氧化碳體系的分子模擬研究[D].廈門:廈門大學(xué),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]俞聯(lián)夢(mèng),劉錦超.二氧化碳—水二元混合流體的分子動(dòng)力學(xué)模擬研究[J].西南民族大學(xué)學(xué)報(bào)：自然科學(xué)版,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]張 陽(yáng).分子模擬在純超臨界流體及其二元混合物體系中的應(yīng)用[D].北京:清華大學(xué),2005.

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

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

[11]陳正隆,徐為人,湯立達(dá).分子模擬的理論與實(shí)踐[M].北京:化學(xué)工業(yè)出版社,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體系擴(kuò)散性質(zhì)的分子動(dòng)力學(xué)模擬[J].化學(xué)工業(yè)與工程,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-)，男，副教授，主要從事不同形貌微納結(jié)構(gòu)的制備、干燥及分子動(dòng)力學(xué)模擬研究，電話：(0411)84986167，E-mail：wbaohe@163.com。

O561

A

1003-3467(2014)08-0031-05