楊華明，張敏華，馬靜

（1 天津大學(xué)石油化工技術(shù)開發(fā)中心，綠色合成與轉(zhuǎn)化教育部重點(diǎn)實(shí)驗(yàn)室，天津 300072； 2 天津化學(xué)化工協(xié)同創(chuàng)新中心，天津 300072）

引 言

丁烯醛又名巴豆醛，是一種無色、易燃，具有催淚性，在水中溶解度適中，在有機(jī)溶劑中易溶的醛類物質(zhì)[1]。其工業(yè)應(yīng)用廣泛，常用于食品防腐劑、酵母和霉菌抑制劑山梨酸的制備，也可用作酒精抑制劑和四乙基鉛的穩(wěn)定劑等[2]。我國現(xiàn)有的乙烯氧化法生產(chǎn)乙醛[3]裝置主要采用的德國弗里德里?！の榈掠邢薰驹O(shè)備[4]，在生產(chǎn)乙醛的同時副產(chǎn)丁烯醛。因此，有效地把副產(chǎn)品丁烯醛從乙醛中進(jìn)行回收利用是提高工廠經(jīng)濟(jì)效益的有效途徑，同時可以獲得高質(zhì)量要求的乙醛產(chǎn)品。

流體的汽液相平衡數(shù)據(jù)是分離技術(shù)以及分離設(shè)備開發(fā)和設(shè)計(jì)的理論基礎(chǔ)。因不受實(shí)驗(yàn)條件、體系穩(wěn)定性等限制，分子模擬法[5-7]成為實(shí)驗(yàn)測定方法外又一重要的相平衡數(shù)據(jù)獲取途徑。同時，其還具有計(jì)算周期短、計(jì)算成本低等優(yōu)勢。Gibbs 系綜Monte Carlo 法[8-10]（GEMC）因不考慮相邊界條件，無須計(jì)算化學(xué)勢，被廣泛地應(yīng)用于流體相平衡的計(jì)算中。Gibbs 系綜Monte Carlo 法計(jì)算汽液相平衡數(shù)據(jù)的準(zhǔn)確性直接與力場相關(guān)。常見的力場有TraPPE力場[11-12]、OPLS 力場[13]、AUA 力場[14]等。

其中，美國明尼蘇達(dá)大學(xué)Siepmann 課題組開發(fā)的TraPPE 力場參數(shù)是通過擬合液相及臨界點(diǎn)實(shí)驗(yàn)數(shù)據(jù)所得，因此在汽液相平衡計(jì)算方面較為準(zhǔn)確。同時，聯(lián)合原子力場（UA）的計(jì)算速率要高于全原子力場（AA）的計(jì)算速率。Siepmann 課題組已經(jīng)分別對正構(gòu)烷烴[8]，直鏈和支鏈烯烴[11]，醚、酮、醛類[12]等體系進(jìn)行了系統(tǒng)性研究，并補(bǔ)充了相關(guān)基團(tuán)的非鍵作用參數(shù)及鍵參數(shù)，將該力場應(yīng)用于相關(guān)物質(zhì)純組分及混合組分的汽液相平衡計(jì)算中，取得了相當(dāng)高的準(zhǔn)確度。然而，TraPPE 力場中關(guān)于烯醛類的參數(shù)還不完善，如缺少烯醛類物質(zhì)的鍵長和二面角扭轉(zhuǎn)項(xiàng)參數(shù)。本文以TraPPE-UA 力場為基礎(chǔ)，采用基于第一性原理的量子化學(xué)方法對缺少的鍵長參數(shù)和二面角扭轉(zhuǎn)項(xiàng)參數(shù)進(jìn)行了計(jì)算，補(bǔ)充了丁烯醛中缺失的力場參數(shù)。采用構(gòu)建的力場參數(shù)計(jì)算了丁烯醛純組分、丁烯醛與乙醛二元體系的汽液相平衡數(shù)據(jù)，并將計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了比較，以期為化工分離過程提供基礎(chǔ)數(shù)據(jù)。

1 力場模型

TraPPE-UA 為TraPPE 聯(lián)合原子力場，即將CHx（0≤x≤4）基團(tuán)作為一擬原子，作用中心在C 原子上。鍵相互作用參數(shù)[11-12]中，鍵長采用固定鍵長形式（表1）。

表1 TraPPE-UA 力場中丁烯醛及乙醛鍵參數(shù)Table 1 Bond parameters of crotonaldehyde and ethanal in TraPPE-UA force field

被兩原子鍵分開的三原子間通過諧波彎曲勢能相互作用

式中，kθ、θ、θeq分別為力常數(shù)、鍵角以及平衡鍵角（表2）。

被三原子鍵分開的四原子間通過余弦形式的扭轉(zhuǎn)勢能交互作用

式中，φ和ci分別為扭轉(zhuǎn)角和系數(shù)項(xiàng)。而中間鍵為雙鍵的反式結(jié)構(gòu)，采用固定二面角，d為固定二面角扭轉(zhuǎn)項(xiàng)參數(shù)（表3）。

表2 TraPPE-UA 力場中丁烯醛及乙醛彎曲項(xiàng)參數(shù)Table 2 Bend parameters of crotonaldehyde and ethanal in TraPPE-UA force field

表3 TraPPE-UA 力場中丁烯醛及乙醛扭轉(zhuǎn)項(xiàng)參數(shù)Table 3 Torsion parameter of crotonaldehyde and ethanal in TraPPE-UA force field

非鍵交互作用是由Lennard-Jones（LJ）作用和庫侖勢能組成

式中，rij、εij、σij、qi和qj分別為原子間距、LJ阱深參數(shù)、LJ 半徑，以及原子i和j的局部電荷（表4）。本計(jì)算過程中，LJ 參數(shù)通過遷移TraPPE 力場參數(shù)[11-12]獲得，電荷進(jìn)行自動分配。

表4 TraPPE-UA 力場中丁烯醛及乙醛非鍵相互作用參數(shù)Table 4 TraPPE-UA force field parameters for nonbonded interactions of crotonaldehyde and ethanal

分子間相互作用通過Lorentz-Berthelot（LB）混合規(guī)則進(jìn)行計(jì)算

2 計(jì)算細(xì)節(jié)

2.1 鍵長及扭轉(zhuǎn)項(xiàng)參數(shù)計(jì)算

2.2 GEMC 計(jì)算

計(jì)算使用MCCCS Towhee[17]軟件，采用補(bǔ)充后的TraPPE-UA 力場參數(shù)對丁烯醛純物質(zhì)、丁烯醛與乙醛二元體系進(jìn)行了GEMC 計(jì)算。丁烯醛純組分計(jì)算使用500 個分子，丁烯醛與乙醛混合物計(jì)算共使用500 個分子。計(jì)算首先進(jìn)行2000000 步MC 運(yùn)動進(jìn)行預(yù)平衡，之后進(jìn)行3000000 步MC 運(yùn)動進(jìn)行平衡。通過比較預(yù)平衡中的塊平均密度，2000000 步足夠使體系趨于平衡狀態(tài)。模擬過程中，體積變化概率、粒子交換概率、分子再生長概率、平動和轉(zhuǎn)動概率分別為0.2%、15.8%、28%、28%和28%。對于二元組分的計(jì)算，初始體積及分子數(shù)的設(shè)定值需要同最終平衡值進(jìn)行比較，并保證它們之間的差值不大[18]。

3 結(jié)果與討論

3.1 鍵參數(shù)

分別采用B3LYP 和MP2 方法對丁烯醛分子進(jìn)行結(jié)構(gòu)優(yōu)化（圖1），其鍵長值見表5。兩種方法所得鍵長參數(shù)相差0.0001 nm，與實(shí)驗(yàn)值[1]相比相差0.0002 nm。本文選用0.1469 nm 作為TraPPE-UA 力場中丁烯醛缺少的CH CH 鍵參數(shù)，進(jìn)行GEMC計(jì)算。

圖1 丁烯醛優(yōu)化后結(jié)構(gòu)Fig.1 Optimizational structure of crotonaldehyde

表5 TraPPE-UA 力場中CH CH 鍵長參數(shù)Table 5 Bond parameter of CH CH in TraPPE-UA force field

對CH CH CH O 二面角進(jìn)行扭轉(zhuǎn)并分別采用B3LYP 和MP2 方法計(jì)算新構(gòu)型能量，將計(jì)算所得扭轉(zhuǎn)能代入式(2)進(jìn)行擬合獲得扭轉(zhuǎn)項(xiàng)參數(shù)。擬合結(jié)果如圖2、表6所示。

表6 TraPPE-UA 力場中CH CH CH O 二面角扭轉(zhuǎn)項(xiàng)參數(shù)Table 6 Torsional parameters of CH CH CH O in TraPPE-UA force field

圖2中采用B3LYP 和MP2 方法獲得的扭轉(zhuǎn)能曲線與擬合曲線吻合度較高，相關(guān)系數(shù)分別為99.862%和99.897%。本文選用擬合度稍高的MP2方法計(jì)算獲得的二面角參數(shù)作為TraPPE-UA 力場中丁烯醛缺少的CH CH CH O 二面角參數(shù)，進(jìn) 行GEMC 計(jì)算。

圖2 丁烯醛分子中CHCHCHO 二面角 扭轉(zhuǎn)能擬合曲線Fig.2 Torsional energies and fitting curves of CHCHCHO in crotonaldehyde

3.2 丁烯醛純組分汽液相平衡模擬

本文采用TraPPE-UA 力場，并采用新補(bǔ)充的鍵長和扭轉(zhuǎn)項(xiàng)參數(shù)進(jìn)行丁烯醛純組分汽液相平衡模擬。臨界溫度和臨界密度通過式(5)、式(6)[19-20]計(jì)算得到

式中，ρl、ρv、ρc分別為液相、氣相和臨界密度，T、Tc分別為對應(yīng)溫度和臨界溫度，A、B為可調(diào)參數(shù)，β為臨界指數(shù)，一般設(shè)定為0.32。并將結(jié)果與實(shí)驗(yàn)值進(jìn)行比較。計(jì)算結(jié)果如圖3、表7所示。

從圖3、表7中可以看出，丁烯醛液相密度的模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)[21]相比，其誤差最低為0.28%，最高為1.23%。產(chǎn)生誤差的原因有：丁烯醛體系含有順式丁烯醛、反式丁烯醛兩種構(gòu)型，其中反式丁烯醛含量高達(dá)95%以上，模擬過程中對此做近似處理；量子化學(xué)計(jì)算過程中，選用的方法本身存在微小誤差。臨界溫度模擬結(jié)果為562.77 K，與實(shí)驗(yàn)值571.15 K 相比誤差為1.46%，臨界密度模擬值為288.23 kg·m-3，與實(shí)驗(yàn)值280.36 kg·m-3相比誤差為2.81%。通過將模擬值與實(shí)驗(yàn)值比較發(fā)現(xiàn)，采用TraPPE-UA 力場計(jì)算的汽液相平衡數(shù)據(jù)具有較高的準(zhǔn)確性。

圖3 丁烯醛純組分汽液相平衡數(shù)據(jù)Fig.3 Vapor-liquid equilibria data of pure crotonaldehyde

表7 采用TraPPE-UA 力場計(jì)算的丁烯醛純組分 汽液相平衡數(shù)據(jù)Table 7 Vapor-liquid equilibria data of pure crotonaldehyde calculated by TraPPE-UA force field

3.3 丁烯醛、乙醛二元體系的汽液相平衡模擬

本文采用TraPPE-UA 力場，并采用新補(bǔ)充的鍵長和扭轉(zhuǎn)項(xiàng)參數(shù)對丁烯醛、乙醛二元體系的汽液相平衡進(jìn)行了計(jì)算，壓力條件為101325 Pa，并將模擬結(jié)果與實(shí)驗(yàn)值[3]進(jìn)行比較。實(shí)驗(yàn)數(shù)據(jù)是在101325 Pa條件下，通過Rose 釜測量所得，經(jīng)過活度系數(shù)方程驗(yàn)證，平均誤差最高為0.56%，計(jì)算結(jié)果如圖4、表8所示。

圖4 丁烯醛與乙醛組分汽液相平衡數(shù)據(jù)Fig.4 Vapor-liquid equilibria data of crotonaldehyde and ethanal mixture

表8 采用TraPPE-UA 力場計(jì)算的丁烯醛與乙醛二元體系汽液相平衡數(shù)據(jù)Table 8 Vapor-liquid equilibria data of crotonaldehyde and ethanal mixture calculated by TraPPE-UA force field

從圖4、表8中可以看出，丁烯醛與乙醛二元體系的汽液相平衡數(shù)據(jù)的模擬結(jié)果與實(shí)驗(yàn)值相比，液相誤差范圍為0.82%～7.44%，氣相誤差范圍為 0.13%～7.09%。相比于丁烯醛純組分，丁烯醛與乙醛二元體系模擬結(jié)果的準(zhǔn)確性略差，這與二元體系的相互作用相比于純組分更為復(fù)雜有關(guān)，需要更準(zhǔn)確的力場參數(shù)對其進(jìn)行描述。

4 結(jié) 論

本文采用基于第一性原理的量子化學(xué)方法，對丁烯醛在TraPPE-UA 力場缺失的鍵參數(shù)和二面角扭轉(zhuǎn)項(xiàng)參數(shù)進(jìn)行了計(jì)算，完善了丁烯醛分子的TraPPE-UA力場參數(shù)。采用GEMC 方法，對丁烯醛純組分、丁烯醛與乙醛二元體系進(jìn)行了汽液相平衡計(jì)算，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對比。模擬結(jié)果顯示了較高的準(zhǔn)確性。對于丁烯醛與乙醛二元混合體系，基于二元體系相互作用的復(fù)雜性，模擬結(jié)果的準(zhǔn)確性要低于純組分。在今后的研究中，還需要更準(zhǔn)確的力場參數(shù)來對其進(jìn)行描述，其普遍性也有待進(jìn)一步的驗(yàn)證。

符 號 說 明

A，B——參數(shù)

ci——二面角扭轉(zhuǎn)項(xiàng)公式系數(shù)，J

d——固定二面角扭轉(zhuǎn)項(xiàng)參數(shù)，K

kθ——鍵角彎曲力常數(shù)，J

r——鍵長，nm

rij——原子間距，nm

T——溫度，K

ubend——諧波彎曲勢能，K

u(rij)——非鍵相互作用能，K

utors——扭轉(zhuǎn)勢能，K

x——液相摩爾分?jǐn)?shù)

y——?dú)庀嗄柗謹(jǐn)?shù)

β——臨界指數(shù)

Δdev——誤差值，%

εij——LJ 阱深參數(shù)，J

θ——鍵角彎曲度，（°）

θeq——平衡鍵角，（°）

ρ——密度，kg·m-3

σij——LJ 半徑，nm

?——二面角角度，rad

下角標(biāo)

bend——鍵角

cal.——模擬計(jì)算值

dev——誤差

eq——平衡態(tài)

exp.——實(shí)驗(yàn)值

i，j——組分i，j

tors——二面角

[1]Jayaprakash A,Arjunan V,Jose S P,Mohan S.Vibrational and electronic investigations,thermodynamic parameters,HOMO and LUMO analysis on crotonaldehyde byab initioand DFT methods [J].Spectrochimica Acta Part A:Molecular and Biomolecular Spectroscopy,2011,83:411-419

[2]Fernandez J E,Solomons T W G.Crotonaldehyde [J].Chem.Rev.,1962,62:485-502

[3]Shanghai Wusong Chemical Plant (上海吳淞化工廠),Department of Chemical Engineering,Zhejiang University (浙江大學(xué)化工系).Vapor-liquid equilibria in ethanol-crotonaldehyde-water ternary system [J].Chemical Engineering(China)（化學(xué)工程），1976 (6):1-31

[4]Fan Wenge (范文革),Zhang Deguang (張德光),Hu Jiwei (胡繼偉).Discussion on way out for China’s existing acetaldehyde production plants using ethylene oxidation process [J].Refining and Chemical Industry(煉油與化工),2007,18 (1):8-9

[5]Li Xiaofeng,Zhao Lifeng,Cheng Tao,Liu Lianchi,Sun Huai.One force field for predicting multiple thermodynamic properties of liquid and vapor ethylene oxide [J].Fluid Phase Equilb.,2008,274 (1/2):36-43

[6]Xu W H,Yang J C.Computer simulations on aggregation of acetic acid in the gas phase,liquid phase,and supercritical carbon dioxide [J].J.Phys.Chem.A,2010,114:5377-5388

[7]Shi Xiaoming (施曉明).Vapor-liquid equilibria simulation by GEMC [D].Chongqing:Chongqing University,2006

[8]Martin M G,Siepmann J I.Transferable potentials for phase equilibria (Ⅰ):United-atom description ofn-alkanes [J].Journal of Physical Chemistry B,1998,102 (14):2569-2577

[9]Panagiotopoulos A Z.Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble [J].Molecular Physics,1987,61 (4):813-826

[10]Panagiotopoulos A Z,Quirke N,Stapleton M,Tildesley D J.Phase equilibria by simulation in the Gibbs ensemble:alternative derivation,generalization and application to mixture and membrane equilibria [J].Molecular Physics,1988,63 (4):527-545

[11]Wick C D,Martin M G,Siepmann J I.Transferable potentials for phase equilibria (Ⅳ):United-atom description of linear and branched alkenes and alkylbenzenes [J].J.Phys.Chem.B,1998,102 (38):2569-2577

[12]Stubbs J M,Potoff J J,Siepmann J I.Transferable potentials for phase equilibria (Ⅵ):United-atom description for ethers,glycols,ketones and aldehydes [J].J.Phys.Chem.B,2004,108:17596-17605

[13]Jorgensen W L,Maxwell D S,Tirado-Rives J.Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids [J].J.Am.Chem.Soc.,1996,118 (45):11225-11236

[14]Ungerer P,Nieto-Draghi C,Rousseau B.Molecule simulation of the thermophysical properties of fluids:from understanding toward quantitative predictions [J].J.Mol.Liq.,2007,134:71-89

[15]Becke A D.Density functional thermochemistry (Ⅲ):The role of exact exchange [J].J.Chem.Phys.,1993,98:5648-5652

[16]Stephens P J,Devlin F J,Chabalowski C F,Frisch M J.Ab initiocalculation of vibration absorption and circular dichroism spectra using density functional force fields [J].J.Phys.Chem.,1994,98 (45):11623-11627

[17]Martin M G.MCCCS Towhee:a tool for Monte Carlo molecular simulation [J].Molecular Simulation,2013,39 (14/15):1212-1222

[18]McKnight T J,Vlugt T J H,Ramjugernath D,Starzak M,Ahlstr?m P,Bolton K.Simulation of 1-alkene and alkane binary vapour-liquid equilibrium using different united-atom transferable force fields [J].Fluid Phase Equilibria,2005,232 (1/2):136-148

[19]Rowildson J S,Widom B.Molecular Theory of Capillarity [M].New York:Oxford University Press,1989:261

[20]Rowildson J S,Widom B.Liquids and Liquid Mixtures [M].3rd ed.London:Butterworth,1982:70-75

[21]Yaws C L.Chemical Properties Handbook [M].New York:McGraw-Hill,1999