Molecular Dynamics Simulation of Transport and Structural Properties

Article pubs.acs.org/jced

Molecular Dynamics Simulation of Transport and Structural Properties of CO2 Using Different Molecular Models Haimin Zhong,† Shuhui Lai,† Jinyang Wang,† Wenda Qiu,† Hans-Dietrich Lüdemann,‡ and Liuping Chen*,† †

KLGHEI of Environment and Energy Chemistry, School of Chemistry and Chemical Engineering, Sun Yat-sen University, Guangzhou 510275, People’s Republic of China ‡ Universität Regensburg, Institut für Biophysik und Physikalische Biochemie, D-93040 Regensburg, Germany S Supporting Information *

ABSTRACT: The diffusion coefficients (Ds), viscosities (η), and structural properties of carbon dioxide (CO2) have been studied using molecular dynamics (MD) simulation. Three fully flexible models (MSM-flex, EPM2flex, and TraPPE-flex) from the literature are used to model CO2. Present simulations have extended the temperature range from 223 K to 450 K and pressures up to 200 MPa for the first time. Generally, the simulation results show a good agreement with the experimental ones. The overall satisfaction of the EPM2-flex model is found to be the best, with an average absolute relative deviation of 6.83 % for Ds and of 2.87 % for η, respectively. However, the TraPPE-flex model performs best at low temperatures below 273 K. Meanwhile, the lifetime of CO2 molecules in the first solvation shell (τs) is calculated, and the qualitative correlation between τs and Ds as well as τs and η is discussed. Finally, the structures of CO2 fluid in different thermodynamic states are investigated by calculating radial distribution functions and using a clustering algorithm.

1. INTRODUCTION Carbon dioxide (CO2) satisfies several principles of green chemistry and engineering, including pollution prevention, lower toxicity, abundantly available, inexpensive, easy to separate, etc.1 At present, it has many industrial applications2−6 such as supercritical extraction, dyeing, drying, and reaction engineering, etc. Knowing the diffusion coefficients and viscosity of CO2 is essential for efficient chemical process design. Therefore, there is a growing demand for transport properties of CO2 in a very broad range of temperature and pressure. However, it is difficult and expensive to determine self-diffusion coefficients7−11 and viscosities12−14 experimentally because of the rigorous experimental conditions required, especially at some extreme pressures and temperatures. It is worth mentioning that Groß et al.11 has measured the selfdiffusion coefficients of CO2 in the broadest temperature range from 223 K to 450 K and pressure range from 10 MPa to 200 MPa. Nowadays, molecular dynamics (MD) simulations have became an important method to obtain the self-diffusion coefficient and viscosity of CO2. Moreover, such simulations can provide detailed and valuable characteristic pictures of a molecular structure at the atomic level. In recent years, there have been many studies on the transport and structural properties of CO2 using MD simulation. Merker et al.15 calculated the self-diffusion coefficient and viscosity of CO2 by using MD simulation. Nieto-Draghi et al.16 used the rigid and fully flexible EMP2 model to simulate the transport properties of CO2. Bock et al.17 © XXXX American Chemical Society

studied the transport properties of CO 2 through an intermolecular potential using a classical trajectory approach. Zhang and Duan18 conducted MD simulations to obtain the self-diffusion coefficient of CO2 and proposed an optimized model for CO2. Aimoli et al.19−21 computed the transport and thermodynamics properties of CO2, methane, and their binary mixtures by using several different models for MD simulation. It is worth mentioning that their results are comparable to ours, and it will be discussed later. To the best of our knowledge, the above-mentioned simulations are performed at temperatures above 273 K and using rigid or semiflexible models to model CO2. In this work, we add flexible bonds and angles to three common rigid CO2 models and perform MD simulations to investigate the transport and structural properties of CO2 at temperatures from 223 K to 450 K and pressures from 10 MPa to 200 MPa. The simulated results are compared with experimental values to evaluate the performance of our fully flexible models and simulation techniques. The comprehensive comparisons of the performances of these models are detailedly discussed. Besides, the structures of CO2 in different thermodynamic states are explored at the atomic level. Received: October 14, 2014 Accepted: June 15, 2015

A

DOI: 10.1021/je5009526 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

Table 1. Force Field Parameters for Each CO2 Modela nonbonded interactions εC−C/kB

σC−C

εO−O/kB

σO−O

bond stretching

angle stretching

bC−O

θ0

qC

models

K

Å

K

Å

e

Å

deg

MSM22,23 EPM224 TraPPE25

29.0 28.129 27.0

2.785 2.757 2.80

83.1 80.507 79.0

3.014 3.033 3.05

+0.5957 +0.6512 +0.70

1.16 1.149 1.16

180 180 180

a The parameters between unlike atoms are calculated by the Lorentz−Berthelot combining rule; kb = 10739.337 kJ·mol−1·Å−2 (from ref 16) and kθ = 1236 kJ·mol−1·rad−2 (from ref 24) are added to the above rigid models to make them fully flexible, as expressed in eq 1.

2. COMPUTATIONAL METHODS 2.1. Potential Model. Great efforts have been devoted to the development of accurate potential model for CO2 molecule and many models have been proposed. Site−site interaction models which designate interaction centers on the three atoms are most popular for CO2. In this study, we add bond stretching and angle bending terms to three common site−site interaction rigid models: MSM,22,23 EPM2,24 (elementary physical model 2), and TraPPE25 (transferable potentials for phase equilibria) models. Then, they are renamed as MSM-flex, EPM2-flex, and TraPPE-flex models. All potential functions calculated in the present work can be expressed as a summation of bond stretching potential, angle bending potential, Lennard-Jones (LJ) potential, and Coulombic interaction, defined as 1 1 E = ∑ kb(b − b0)2 + ∑ kθ(θ − θ0)2 bond 2 angle 2 ⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ qiqj ⎫ σij ⎥ ⎪ ⎪ ⎢ σij ⎜ ⎟ ⎜ ⎟ ⎬ + ∑ ⎨4εij⎢⎜ ⎟ − ⎜ ⎟ ⎥ + r ⎝ rij ⎠ ⎦ 4πε0rij ⎪ nonbond ⎪ ⎭ ⎩ ⎣⎝ ij ⎠

representative empirical equations, namely the NAM equation,27 EA equation,28 MCZ equation,29,30 and TBH equation31 to calculate the density of CO2, respectively. The comparison of density calculated through different empirical equations with experimental ones is plotted in Figure 1. In this figure, one can

(1)

Figure 1. Comparison between densities calculated using different empirical equations and experimental data26 at 223 and 243 K.

where kb, b, and b0 are the bond stretching force constant, bond length, and equilibrium bond length, respectively; kθ, θ, and θ0 are the angle bending force constant, bending angle, and equilibrium bending angle, respectively; rij, εij, σij, qi, and qj are the separation, LJ well depth, LJ size, and partial charges, respectively, for the pair of atoms i and j; ε0 is the dielectric constant of vacuum. All the detailed potential parameters are collected in Table 1. Here, kθ = 1236 kJ·mol−1·rad−2 is taken from Harris et al;24 kb = 10739.337 kJ·mol−1·Å−2 is taken from Nieto-Draghi et al.16 2.2. Simulation Details. Simulations were performed with GROMACS 4.0.7 software (double precision). Periodic boundary conditions in the three directions were employed for all the simulations. Leap-frog algorithm was used for integrating Newton’s equations of motion. Temperature coupling was realized through the velocity rescaling method. Fast Particle-Mesh Ewald (PME) electrostatics method was used to compute the electrostatic interactions. The cutoff distances for Lennard-Jones interactions and Coulombic interactions were both taken to be 1.4 nm. Simulations were carried out in the canonical (NVT) ensemble, so density of CO2 was needed to determine the size of simulation box. In this simulation, we mainly used the experimental density taken from National Institute of Standards and Technology database (NIST)26 to calculate the size of cubic box. However, densities in the temperature range from 223 K to 243 K and pressures up to 200 MPa had no experimental values. Therefore, densities of CO2 in these thermodynamic states were obtained from empirical equations. To obtain the most reliable density, we adopted four

see that the TBH equation represents better the experimental values at 243 K, whereas the NAM equation represents better the data at 223 K. For brevity, the details of empirical equation for calculating density and the density data applied in this simulation are presented in Supporting Information. The self-diffusion coefficient is calculated by the Einstein relation method. It can be expressed by the following equation: 1 Ds = lim ⟨|r(t ) − r(0)|2 ⟩ (2) t →∞ 6Nt where Ds is the self-diffusion coefficient of CO2, t is the time, and r(t) is the position of a molecule at time t. Viscosity is determined by transverse-current correlation functions (TACFs) for plane waves based on momentum fluctuations.32,33 The behavior of the velocity u(r, t) of a fluid is described by the Navier−Stokes equation, as follows: ρ

∂u + ρ(u·∇)u = ρa − ∇p + η∇2 u ∂t

(3)

where a is the external force per unit of mass and volume, u is the velocity, ρ is the number density, η is the viscosity, p is the pressure. Considering an incompressible liquid with an initial velocity field consisting of a plane wave in one direction (x), the y and z-component of u will remain zero. At the xcomponent, the Navier−Stokes equation can be reduced to 2 ∂ux(z , t ) η ∂ ux(z , t ) = ∂t ρ ∂z 2

B

(4) DOI: 10.1021/je5009526 J. Chem. Eng. Data XXXX, XXX, XXX−XXX


Article

Figure 2. (a) Particle number dependence and (b) statistical time dependence of the simulated self-diffusion coefficient and viscosity for CO2 at 298 K.

Because inversion is asymmetry of the system, η(k) must be an even function of k. This implies that

where z is the coordinate value in z-component, ux(z,t) is the velocity in x-component at t moment and z coordinate value. The solution of eq 4 is ux(z , t ) = u0e−t / τr cos(kz) τr =

η(k) = η(0) + ak 2

(5)

ρ ηk 2

all the η(k) and k values from the simulations should be fitted to eq 12, which can be extrapolated to k = 0 in order to obtain the viscosity. To make objective opinions of the quality of a given model, the average absolute relative deviation (AARD) between the calculated value and simulated one is calculated as follows:

(6)

We see that plane waves decay exponentially with a time constant which is inversely proportional to the viscosity and proportional to the wavelength squared. However, the correlation function is not a pure exponential at short times. To account for this, a phenomenological correction is applied by incorporating a relaxation-time. This changes eq 4 to ∂ux(z , t ) η = ∂t ρ

∫o

t

ϕ(t − s)

∂ 2ux(z , s) ∂z 2

ds

AARD % =

1 −t / τm e τm

(8)

where τm is the relaxation-time. The solution of this equation with initial condition is

RD % =

−β ⎛ ⎜

⎞ 1 ux(z , t ) = u0e cosh(Ωβ) + sinh(Ωβ )⎟ cos(kz) ⎝ ⎠ Ω (9)

g (r ) =

1 − 4τ

i=1

Aicalc − Aiexpt ·100 Aiexpt

(13)

Aicalc − Aiexpt · 100 Aiexpt

(14)

dN ρ 4πr 2 dr

(15)

(10)

where g(r) is the RDF, N is the number of CO2 molecules, and r is the separation distance between mass centers for CO2 molecules. Various systems consisting of 300, 500, 700, and 1000 particles were used to analyze finite-size effects. The results are shown in Figure 2a. It can be seen that 500, 700, and 1000 molecules give almost the same accuracy. Therefore, 500 CO2 molecules are enough for our simulation. As seen in Figure 2b, we also discussed the statistical time effects in calculating sefdiffusion coeffient and viscosity. It indicates that statistical time of 3 ns is sufficient for present simulation. Hence, the simulation system contained 500 CO2 molecules. Initially, runs of 3 ns were taken to re-equilibrate the system, and then runs of 3 ns were used to calculate the transport and structural properties for CO2. The time step is 1 fs.

and Ω=

N

∑

The radial distribution function (RDF) between mass centers for different CO2 molecules is obtained as

where

t β= 2τ

1 N

where N is the number of points, Aexpt (diffusion coefficient or viscosity) is the experimental value, and Acalc is the value obtained by MD simulation. Similarly, the relative deviation (RD) between calculated value and simulated one is counted by the expression

(7)

The simplest model for the memory ϕ is ϕ(t ) =

(12)

η (k ) 2 k ρ

(11)

where u0 is the initial velocity, τ is decay constant of transverse momentum fields, k is the k-vector, η(k) is the viscosity at k value. For more details, the reader is referred to ref 33. In the present simulation, the TACFs were calculated using the k-vectors (1, 0, 0) and (2, 0, 0) each also in the y- and zdirection, (1, 1, 0) and (1, −1, 0) each also in the two other plains (these vectors are not independent), and (1, 1, 1) and the other box diagonals (also not independent). Each of these TACFs was fitted to eq 9. These provided a total of 16 different values of k and η(k) for each thermodynamic state of fluid. C



Article

3. RESULTS AND DISCUSSION 3.1. Transport Properties. First, we checked the reliability and superiority of our simulations by comparing our results with corresponding experimental values and simulated data from literature. The comparisons are presented in Table 2. Our Table 2. Comparison of Self-Diffusion Coefficient (Ds) and Viscosity (η) of CO2 Predicted by Different Models with Experimental (Expt) Data at 373 K and 30 MPa models

Ds/10−9 m2·s−1

η/10−5 Pa·s

ref

33.7 34.0 31.2 31.6 31.9 32.3 30.8 31.3

5.41 5.32

11, 26 this work 18 this work this work this work 18 15

Expt EPM2-flex EPM2-rigid TraPPE-flex MSM-flex Zhang and Duan Merker et al. a

5.22 5.31 5.19 6.25a

Figure 3. Comparison of self-diffusion coefficients for CO2 calculated using MD simulations or empirical equations with the experimental values.11

Calculated by Green−Kobu formalism.

simulation results using fully flexible models show a good agreement with experimental data, and seem to be better than the simulated values from literature. In addition, as η taken from literature was calculated by Green−Kubo formalism, the better simulation results of ours may attribute to the technique progress of TACFs method. In Table 3, we also compare the

experimental data in different temperature ranges are listed in Table 4. Obviously, the empirical equation performs worse than the MD simulation. Moreover, the overall satisfaction of the EPM2-flex model is the best, that of the MSM-flex model is the middling, and that of the TraPPE-flex model is the worst. However, the TraPPE-flex model performs very well at low temperatures (223 K and 243 K). In fact, the EPM2 model is modified from the MSM model.24 On the other hand, we suggest that the kb and kθ for the TraPPE-flex model may need to be further optimized. As mentioned above, Aimoli et al.19 used seven different CO2 models (including single-site, rigid, and fully flexible models) to calculate the Ds of CO2 in a very wide temperature range from 273.15 K to 573.15 K. In their studies, the Cygan model (a flexible extension of the EPM2 model) was found to show the best accuracy with an AARD of 4.85 %, which is superior to our EPM2-flex model with an AARD of 6.83 %. However, our simulations are conducted at temperatures from 223 K to 450 K, which was extended to a lower temperature range. As we pointed out in this work, temperature plays an important role on the performance of a model. To further explore the temperature and pressure dependences of these models, we have calculated the relative deviations between simulated values and experimental ones. The temperature and pressure dependences of relative deviation are plotted in Figure 4, in which for a given temperature or pressure, the points are the results of different simulations for the different states of CO2 fluid. Figure 4a clearly indicates that the EPM2-flex and MSM-flex models perform badly at 223 K and 243 K, whereas the TraPPE-flex model performs well at those low temperatures. That is to say a certain model is only suitable in a limited temperature range. This may happen because all the force parameters of these models are obtained by fitting the macro properties of CO2 at a confined temperature range. Therefore, the adjustment of those parameters to make the model applicative in a broader range of temperature is still a task. In Figure 4b, it can be noticed that the distribution of relative deviation for different models varies randomly with pressure; in another word, pressure does not affect much of the behaviors of different models. 3.1.2. Comparison of Simulated η and Experimental Data. Similar to the self-diffusion coefficient, the viscosity of CO2 at

Table 3. Self-Diffusion Coefficient (Ds) and Viscosity (η) Simulated Using Fully Flexible or Rigid EPM2 Model, And Their Comparisons with the Experimental Values at 50 MPa Ds/10−9 m2·s−1

a

η/10−5 Pa·s

T/K

expta

flexible

rigid

exptb

flexible

rigid

424 373 333 298 273 243

36.6 24.3 16.0 11.0 8.42 5.19

35.0 22.1 15.3 11.0 8.39 5.56

33.8 21.0 13.0 11.4 9.00 5.94

5.88 7.61 10.0 13.5 17.3 24.3

5.87 7.81 10.1 13.8 17.1 22.4

5.74 7.33 10.4 14.1 15.6 22.0

Experimental Ds from ref 11. bExperimental η from NIST.26

simulation results between rigid and fully flexible EPM2 model. As seen, the introduction of molecular flexibility to rigid EPM2 model seems to slightly improve the accuracy of simulated Ds and η calculated by TACFs method. The effect of molecular flexibility on the Ds calculation accords with the results reported by Aimoli et al.19 3.1.1. Comparison of Simulated Ds and Experimental Data. We have obtained the simulated self-diffusion coefficients of CO2 in the temperature range from 223 K to 450 K and pressure range from 10 MPa to 200 MPa. The simulated Ds are then compared to the experimental values measured by Groß et al.11 to evaluate the performances of our fully flexible models. It is worth mentioning that we also use three representative empirical equations, namely the LJ-2 equation, 34 LJ-4 equation,35 and EWJ equation36,37 to calculate the self-diffusion coefficients of CO2. The comparison between the self-diffusion coefficient calculated by MD simulation or empirical equation and the corresponding experimental value is plotted in Figure 3. The average absolute relative deviations (AARDs) between Ds calculated using the MD simulation or empirical equation and D



Article

Table 4. Average Absolute Relative Deviations (AARDs) of Ds and η for CO2 Calculated Using MD Simulation or Empirical Equation When Compared to Experimental Data at Pressures from 10 MPa to 200 MPa AARD/% MD simulation MSM-flex

empirical equation

EPM2-flex

TraPPE-flex

LJ-2

LJ-4

EWJ

T/K

Ds

η

Ds

η

Ds

η

Ds

Ds

Ds

223 to 243 273 to 450

23.47 6.14

10.33 2.89

13.46 5.03

7.78 2.49

4.53 11.38

2.46 6.60

8.42 21.63

30.82 13.72

28.58 11.50

Figure 4. Relative deviation between simulated self-diffusion coefficients by using different CO2 models and experimental ones: (a) temperature dependence; (b) pressure dependence.

Figure 5. Comparison of simulated viscosity using different models with the experimental data taken from NIST26 at (a) 450 K and 424 K, (b) 373 K and 333 K, (c) 298 K and 273 K, (d) 243 K and 223 K, in the pressure range from 10 MPa to 200 MPa.

without experimental values. The comparisons between simulated η and experimental values are plotted in Figure 5. Generally speaking, the simulated values agree very well with the experimental ones. The AARDs between η calculated using MD simulation and experimental data in different temperature ranges are also collected in Table 4. As seen, the EMP2-flex

temperatures from 223 K to 450 K and pressures from 10 MPa to 200 MPa is acquired in this simulation. Then the simulated viscosity is also compared with the experimental values taken from the NIST.26 It is worth mentioning that we also predict vicosities for CO2 in certain thermodynamic states, such as at temperatures of 223 K and 243 K and pressures up to 200 MPa E



Article

Figure 6. Relative deviation between simulated viscosities by using different CO2 models and experimental ones: (a) temperature dependence; (b) pressure dependence.

model achieves the best overall satisfaction. Compared with the prior work of Aimoli et al.,19 in which rigid EPM2 was reported to show the best estimation for viscosity of an AARD of 4.69 %, the AARD of EPM2-flex model in the present work is only 2.87 %. However, as discussed before, this may be because the temperature range and force parameters in the two simulations are quite different. Additionally, the TACFs method in the present work is also different from the Green−Kubo method used by Aimoli et al.19 Meanwhile, the distributions of relative deviations between simulated viscosities and experimental values for different models over temperature and pressure are also presented in Figure 6. As seen in Figure 6a, EMP2-flex and MSM-flex models show a good estimation for viscosity, but they behave significantly worse at 223 K and 243 K, whereas the TraPPEflex model performs well at those temperatures. In Figure 6b, it can be seen that the relative deviations vary randomly with pressure. These findings correspond with the change rule obsverved in temperature and pressure dependences of selfdiffusion coefficients. 3.1.3. τs. To further reveal the transport properties of CO2 fluid, we calculated the lifetime of the CO2 molecule in the first solvation shell (τs). In this work τs can be defined as38 τ s = t ⟨N (t , r0) → 0⟩

Figure 7. Temperature and pressure dependences of the lifetime of CO2 molecule in the first solvation shell (τs).

practical because we can directly obtain the Ds which is very useful in industrial chemistry only by using MD simulation, instead of complicated empirical equation correlated to τs. Besides, Ds decreases sharply with τs until τs = 40 ps, after which the decrease of Ds becomes not so obvious. In Figure 8b, by contrast, η increases with increased τs. It should be emphasized that a sharp increase of η at 223 K (τs > 27 ps) can be noticed, which we believe may attributed to the CO2 having a strong tendency of forming the solid phase at such a low temperature if the pressure is slightly increased, and resulting in a large viscosity for CO2. On the whole, τs provides us an another understanding of how the movement of CO2 molecule in the first solvation shell can be qualitatively correlated to the transport properties Ds and η. 3.2. Structural Properties. 3.2.1. Radial Distribution Functions. Furthermore, the radial distribution functions (RDFs) of mass centers for different CO2 molecules are calculated to explore the structure of CO2 fluid. It is worth mentioning that RDFs of mass centers for different CO2 molecules can be experimentally determined.39−42 Prior to investigate the structure of CO2, we checked the reliability of our simulated RDFs by comparing them with the experimental values.39,40 As shown in Figure 9, the simulated RDFs accord well with measured ones. The simulated RDFs show a maximum at the location of approximately 4 Å, which corresponds to the experimental observation by Ishii et al.41 Then the RDFs of mass centers for different CO2 molecules in a wide range of temperatures and pressures are computed,

(16)

where r0 is the radius of the first solvation shell, which has been taken as 6 Å corresponding to the first minimum of g(r), and N(t, r0) is the number of contacts of CO2 molecules in the first salvation shell. τs is counted when the number of contact is zero. The calculated τs values for CO2 at different temperatures and pressures are shown in Figure 7. It can be noted that τs increases with increased pressure, but decreases with increased temperature. The difference between τs at 450 K (supercritical state) and 298 K (near critical sate) is not very large, while the τs at 450 K (supercritical state) is significantly smaller than that at 223 K (subcritical state). That is to say, the CO2 molecules move and escape from the first solvation shell more quickly at high temperatures especially at the supercritical state, which makes supercritical CO2 (scCO2) show very good mass transfer property. To see the correlation between Ds or η and τs, the experimental Ds and η are plotted as a function of τs in Figure 8. As seen in Figure 8a, Ds decreases with increased τs, and they seem to have an inverse relationship. We speculate an empirical equation between Ds and τs may be derived, but this may be not F



Article

Figure 8. Experimental (a) Ds and (b) η plotted vs τs at various temperatures: for a certain temperature, these four points correspond to pressures of (10, 50, 100, and 200) MPa, respectively. Ds at 450 K and 10 MPa, as well as η at 223 K and (50, 100, or 200) MPa are simulated values of this work.

of the cluster is less than the cutoff distance. Here, the cutoff distance is defined as 4 Å.44 That is to say, CO2 dimer is formed if the distance between two carbon atoms is less than 4 Å. The percentages of multimers as a function of the cluster size are calculated to know the microscopic distributions of CO2 molecules at various temperatures and pressures, and the results are plotted in Figure 11. As we have seen, when the temperature is lowered or the pressure is increased, the distributions of multimers are pushed out to favor larger clusters at the cost of depleting the dimers. The forming of large clusters will result in an effective solute molecule size greater than one CO2 molecule. As a consequence, the mobility of CO2 molecules are hindered, thus molecular simulations can give us a view of how the self-diffusion coefficient of CO2 would be reduced at low temperatures or high pressures in the atomic details.

4. CONCLUSIONS In summary, extensive MD simulations are conducted to compute the transport and structural properties of CO2 in a wide range of temperatures and pressures. Carbon dioxide is modeled using three fully flexible models, the MSM-flex model, the EPM2-flex model, and the TraPPE-flex model, respectively. The simulated self-diffusion coefficient and viscosity are compared against the experimental data. Generally, the EPM2-flex model achieves the best overall satisfaction and can simultaneously reproduce both the Ds and η. However, although the TraPPE-flex model shows the worst overall estimation, it is found to perform the best below 273 K. To this point, present simulation reveals the temperature effect on the performances of different CO2 models, and is helpful to the development toward new improved models which will be applicative in a broader range of temperature to predict the transport properties for CO2. The lifetime of a CO2 molecule in the first solvation shell (τs) at high temperatures or low pressures is found to be small, in which CO2 fluid exhibits a high self-diffusion coefficient and low viscosity. Further analysis illustrates that τs can provide us another viewpoint of understanding how the movement of a CO2 molecule in the first solvation shell can be qualitatively correlated to the transport properties Ds and η of fluids. Additionally, from the calculation results of g(r), the structure of CO2 at low pressure or high temperature is noted to be very sensitive to T and p. The calculated clustering results reveal how the self-diffusion coefficient of CO2 would be reduced at low temperatures or high pressures in the atomic details.

Figure 9. Comparison of simulated RDFs using different CO2 models with experimental data:39,40 (a) 473 K, 0.8331 g·cm−3; (b) 312 K, 0.8331 g·cm−3; and (c) 240 K, 1.0889 g·cm−3.

and the results are presented in Figure 10. In Figure 10 panels a and b, one can observe two peaks located at 4 Å and 7.3 Å for g(r) at subcritical states (223 K and 10 MPa or 100 MPa), which accords with the experimental observation by van Tricht et al.43 However, when the temperature is increased, the second peak is fading and finally disappearing at 450 K. This phenomenon is more obvious at the low pressure of 10 MPa. As seen in Figure 10 panels c and d, the pressure almost does not affect the intensity of the two peaks of g(r) except when the temperature is 450 K. Therefore, from the observation of g(r) changing with the temperature and pressure, we find that the fluid structure of CO2 in low pressure or high temperature is very sensitive to T and p. The change of transport property for CO2 is accompanied by the variation of fluid structure. As a result, the Ds and η of scCO2 could be drastically altered by slightly changing the temperature or pressure. 3.2.2. Clustering Results. Also, a single linkage clustering algorithm38 is used to see the multimer in CO2 fluid, which is still not able to be observed by experimental method at present. In this simulation, the single linkage clustering algorithm means adding a structure to a cluster when its distance to any element G



Article

Figure 10. (a, b) Temperature dependence and (c, d) pressure dependence of simulated RDFs.

Figure 11. Percentage of multimer in CO2 fluid as a function of cluster size N (number of molecules) at pressures of (a) 10 MPa, (b) 100 MPa, (c) 200 MPa, and temperatures from 223 K to 450 K.

■

On the whole, the present work provides a powerful tool to obtain the self-diffusion coefficient and viscosity of CO2 of which experimentally determined values may be difficult and expensive, and it is helpful for us to have a deep insight into the transport and structural properties of CO2 at the atomic levels by using MD simulation.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +8620-84115559. Fax: +8620-84112245. Email: cesclp@ mail.sysu.edu.cn. Notes

The authors declare no competing financial interest.

■

ASSOCIATED CONTENT

S Supporting Information *

Empirical equations for calculating densities and self-diffusion coefficients, calculated density and self-diffusion coefficient data, along with simulated self-diffusion coefficient and viscosity data. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/je5009526.

ACKNOWLEDGMENTS

The authors gratefully acknowledge Dr. Huajie Feng (School of Chemistry and Chemical Engineering, Hainan Normal University, Haikou 571158, P. R. China) and Dr. Gao Wei (College of Pharmacy, Guangdong Pharmaceutical University, H



Article

(24) Harris, J. G.; Yung, K. H. Carbon dioxide’s liquid−vapor coexistence curve and critical properties as predicted by a simple molecular model. J. Phys. Chem. 1995, 99, 12021−12024. (25) Potoff, J. J.; Siepmann, J. I. Vapor−liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen. AIChE J. 2001, 47, 1676−1682. (26) Lemmon, E. W.; McLinden, M. O.; Friend, D. G. Thermophysical Properties of Fluid Systems. In NIST Chemistry WebBook, NIST Standard Reference Database Number 69; National Institute of Standards and Technology: Gaithersburg, MD, 2011. (27) Nasrifar, K.; Ayatollahi, S.; Moshfeghian, M. A compressed liquid density correlation. Fluid Phase Equilib. 2000, 168, 149−163. (28) Eslami, H.; Azin, R. Corresponding-states correlation for compressed liquid densities. Fluid Phase Equilib. 2003, 209, 245−254. (29) Chang, C. H.; Zhao, X. A new generalized equation for predicting volumes of compressed liquids. Fluid Phase Equilib. 1990, 58, 231−238. (30) Aalto, M.; Keskinen, K. I.; Aittamaa, J.; Liukkonen, S. An improved correlation for compressed liquid densities of hydrocarbons. Part 1. Pure compounds. Fluid Phase Equilib. 1996, 114, 1−19. (31) Thomson, G. H.; Brobst, K. R.; Hankinson, R. W. An improved correlation for densities of compressed liquids and liquid mixtures. AIChE J. 1982, 28, 671−676. (32) Palmer, B. Transverse-current autocorrelation-function calculations of the shear viscosity for molecular liquids. J. Phys. Rev. E 1994, 49, 359−366. (33) Hess, B. Determining the shear viscosity of model liquids from molecular dynamics simulations. J. Chem. Phys. 2002, 116, 209−217. (34) Liu, H.; silva, C. M.; Macedo, E. A. Unified approach to the selfdiffusion coefficients of dense fluids over wide ranges of temperature and pressurehard-sphere, square-well, Lennard−Jones and real substances. Chem. Eng. Sci. 1998, 53, 2403−2422. (35) Silva, C. M.; Liu, H.; Macedo, E. A. Models for self-diffusion coefficients of dense fluids, including hydrogen-bonding substances. Chem. Eng. Sci. 1998, 53, 2423−2429. (36) Easteal, A. J.; Woolf, L. A.; Jolly, D. L. Self-diffusion in a dense hard-sphere fluid: A molecular dynamics simulation. Physica A 1983, 121, 286−292. (37) Easteal, A. J.; Woolf, L. A.; Jolly, D. L. On the number dependence of hard-spheres diffusion of coefficients from molecular dynamics simulations. Physica A 1984, 127, 344−346. (38) van der Spoel, D.; Lindahl, E.; Hess, B.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Feenstra, K. A.; van Drunen, R.; Berendsen, H. J. C. Gromacs User Manual, version 4.0; Groningen University: 2005; www.gromacs.org. (39) Cipriani, P.; Nardone, M.; Ricci, F. P. Neutron diffraction measurements on CO2 in both undercritical and supercritical states. Physica B 1997, 241, 940−946. (40) Chiappini, S.; Nardone, M.; Ricci, F. P.; Bellissent-Funel, M. C. Neutron diffraction measurements on high pressure supercritical CO2. Mol. Phys. 1996, 89, 975−987. (41) Ishii, R.; Okazaki, S.; Odawara, O.; Okada, I.; Misawa, M.; Fukunaga, T. Structural study of supercritical carbon dioxide by neutron diffraction. Fluid Phase Equilib. 1995, 104, 291−304. (42) Nishikawa, K.; Takematsu, M. X-ray scattering study of carbon dioxide at supercritical states. Chem. Phys. Lett. 1994, 226, 359−363. (43) van Tricht, J. B.; Fredrikze, H.; van der Laan, J. Neutron diffraction study of liquid carbon dioxide at two thermodynamic states. Mol. Phys. 1984, 52, 115−127. (44) Zhang, Y.; Yang, J.; Yu, Y. X. Dielectric constant and density dependence of the structure of supercritical carbon dioxide using a new modified empirical potential model: A Monte Carlo simulation study. J. Phys. Chem. B 2005, 109, 13375−13382.

Guangzhou 510006, P. R. China) for their kind help in MD simulation techniques.

■

REFERENCES

(1) Leitner, W. Green chemistry: Designed to dissolve. Nature 2000, 405, 129−130. (2) Halim, R.; Danquah, M. K.; Webley, P. A. Extraction of oil from microalgae for biodiesel production: A review. Biotechnol. Adv. 2012, 30, 709−732. (3) Long, J. J.; Xu, H. M.; Cui, C. L.; Wei, X. C.; Chen, F.; Cheng, A. K. A novel plant for fabric rope dyeing in supercritical carbon dioxide and its cleaner production. J. Clean. Prod. 2014, 65, 574−582. (4) Zhang, X. T.; Chang, D. W.; Liu, J. R.; Luo, Y. J. Conducting polymer aerogels from supercritical CO2 drying PEDOT-PSS hydrogels. J. Mater. Chem. 2010, 20, 5080−5085. (5) Desimone, J. M.; Guan, Z.; Elsbernd, C. S. Synthesis of fluoropolymers in supercritical carbon dioxide. Science 1992, 257, 945−947. (6) Kim, S. W.; Ahn, J.-P. Polycrystalline nanowires of gadoliniumdoped ceria via random alignment mediated by supercritical carbon dioxide. Sci. Rep. 2013, 3, 1606−1610. (7) Winn, E. B. The temperature dependence of the self-diffusion coefficients of argon, neon, nitrogen, oxygen, carbon dioxide, and methane. Phys. Rev. 1950, 80, 1024−1027. (8) O’Hern, H. A.; Martin, J. J. Diffusion in carbon dioxide at elevated pressures. Ind. Eng. Chem. 1955, 42, 2081−2085. (9) Robinson, R. C.; Stewart, W. E. Self-diffusion in liquid carbon dioxide and propane. Ind. Eng. Chem. Fundam. 1968, 7, 90−95. (10) Krynicki, K.; Meragi, A.-L.; Powles, J. G. Self-diffusion in carbon dioxide near the critical point. Ber. Bunsen. Phys. Chem. 1981, 85, 1153−1154. (11) Groß, T.; Buchhauser, J.; Lüdemann, H.-D. Self-diffusion in fluid carbon dioxide at high pressures. J. Chem. Phys. 1998, 109, 4518− 4522. (12) Iwasaki, H.; Takahashi, M. Viscosity of carbon dioxide and ethane. J. Chem. Phys. 1981, 74, 1930−1943. (13) Hunter, I. N.; Marsh, G.; Matthews, G. P.; Smith, E. B. Argon +carbon dioxide gaseous mixture viscosities and anisotropic pair potential energy functions. Int. J. Thermophys. 1993, 14, 819−833. (14) Fenghour, A.; Wakeham, W. A.; Vesovic, V. The viscosity of carbon dioxide. J. Phys. Chem. Ref. Data 1998, 27, 31−44. (15) Merker, T.; Engin, C.; Vrabec, J.; Hasse, H. Molecular model for carbon dioxide optimized to vapor−liquid equilibria. J. Chem. Phys. 2010, 132, 234512−7. (16) Nieto-Draghi, C.; de Bruin, T.; Pérez-Pellitero, J.; Avalos, J. B.; Mackie, A. D. Thermodynamic and transport properties of carbon dioxide from molecular simulation. J. Chem. Phys. 2007, 126, 064509− 8. (17) Bock, S.; Bich, E.; Vogel, E. A new intermolecular potential energy surface for carbon dioxide from ab initio calculations. J. Chem. Phys. 2002, 117, 2151−2160. (18) Zhang, Z.; Duan, Z. An optimized molecular potential for carbon dioxide. J. Chem. Phys. 2005, 122, 214507−15. (19) Aimoli, C. G.; Maginn, E. J.; Abreu, C. R. A. Transport properties of carbon dioxide and methane from molecular dynamics simulations. J. Chem. Phys. 2014, 141, 134101−13. (20) Aimoli, C. G.; Maginn, E. J.; Abreu, C. R. A. Force field comparison and thermodynamic property calculation of supercritical CO2 and CH4 using molecular dynamics simulations. Fluid Phase Equilib. 2014, 368, 80−90. (21) Aimoli, C. G.; Maginn, E. J.; Abreu, C. R. A. Thermodynamic properties of supercritical mixtures of carbon dioxide and methane: A molecular simulation study. J. Chem. Eng. Data 2014, 59, 3041−3054. (22) Murthy, C. S.; Singer, K.; McDonald, I. R. Interaction site models for carbon dioxide. Mol. Phys. 1981, 44, 135−143. (23) Murthy, C. S.; Oshea, S. F.; McDonald, I. R. Electrostatic interactions in molecular crystals: Lattice dynamics of solid nitrogen and carbon dioxide. Mol. Phys. 1983, 50, 531−541. I


Molecular Dynamics Simulation of Transport and Structural Properties

Recommend Documents