Molecular Dynamics Simulation of a RefrigerantâLubricant Oil System

Subscriber access provided by - Access paid by the | UCSB Libraries

Article

Solubility of Carbon Dioxide in Pentaerythritol Hexanoate: Molecular Dynamics Simulation of a Refrigerant-lubricant Oil System Taisuke Sugii, Eiji Ishii, and Florian Müller-Plathe J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.5b06459 • Publication Date (Web): 19 Aug 2015 Downloaded from http://pubs.acs.org on August 24, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Solubility of Carbon Dioxide in Pentaerythritol Hexanoate: Molecular Dynamics Simulation of a Refrigerant-lubricant Oil System

Taisuke Sugii *†, Eiji Ishii †, and Florian Müller-Plathe ‡

†

Hitachi, Ltd., Research & Development Group, Center for Technology Innovation – Mechanical Engineering, 832-2, Horiguchi, Hitachinaka, Ibaraki, 312-0034, Japan.

‡

Eduard-Zintl-Institut für Anorganische und Physikalische Chemie and Center of Smart Interfaces, Technische Universität Darmstadt, Alarich-Weiss-Straße 4, D-64287 Darmstadt, Germany.

Keywords: solubility, lubricant oil, refrigerant, molecular dynamics simulation, chemical potential

ACS Paragon Plus Environment －1－


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract

We have investigated the solubility and the solvation structure between a refrigerant (carbon dioxide, CO2) and a lubricant oil (pentaerythritol hexanoate, PEC6) by molecular dynamics simulations. First, to investigate the solubility, the vapor-liquid equilibrium pressure was calculated. The chemical potential of the liquid phase and the gas phase were calculated, and the equilibrium state was obtained from the crossing point of these chemical potentials. The equilibrium pressures agreed well with experimental data over a wide range of temperature and mole fraction of CO2. Second, the solvation structure was also investigated on a molecular scale. We found the following characteristics. First, the tails of the lubricant oil are relatively rigid inside the ester groups but flexible beyond. Second, CO2 molecules barely enter the lubricant core as delimited by the ester groups. Third, the double-bonded oxygen atoms of the ester groups are good sorption sites for CO2. Fourth, only a few CO2 molecules are attached to more than one carbonyl oxygen simultaneously. Finally, there is also significant unspecific sorption of CO2 in the alkane tail region. These results indicate that increasing the size of the rigid lubricant core would probably decrease the solubility, whereas increasing the number of polar groups would increase it.


Page 2 of 36

Page 3 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION Refrigerants are widely used in many industrial products, such as air-conditioners and refrigerators. To prevent climate change, new refrigerants with low global warming potential (GWP) are required, and new products with high performance, high reliability, and low environmental burden are strongly demanded. To produce better refrigerants, many of their properties must be understood and improved. One very important property is the solubility between refrigerants and lubricant oils1,2, which are used in compressors to lubricate the sliding surfaces of the moving elements such as pistons, rods, scrolls, screws, etc. Refrigerants, however, are soluble in lubricant oils, at least to some extent. Thus, both form a mixture in the compressor and the refrigeration cycle. If the solubility between a lubricant oil and a refrigerant is not appropriate, this greatly reduces the performance and reliability. For example, if a refrigerant concentration is too large in the lubricant oil, its viscosity and lubricating performance are reduced. It might even cause foam formation. On the other hand, trace amounts of lubricant oil generally spread as oil mist from the compressor into the refrigeration cycle. If the miscibility is too low, phase separation occurs in the cycle, leading to reduced heat transfer. Moreover, the retention of oil induced by the phase separation in the cycle may cause leaking of the lubricant oil between the sliding surfaces and the seizure of the compressor. In summary, the mutual solubility of a refrigerant and a lubricant oil must be controlled, so its mechanism must be properly understood. The solubilities between refrigerants and lubricant oils have been reported in plenty of experimental papers. For example, Wahlström and Vamling reported solubilities between several hydrofluorocaobons (HFCs) and polyol ester (POE) oils3,4. They measured the equilibrium pressure and calculated the Henry’s Law constant. Yokozeki reported the solubility data of various refrigerant/lubricant oil mixtures5 and modeled them with cubic equations of state. He investigated nine refrigerants and eleven lubricant oils, including a



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

hypothetical oil. Very recently, Sun et al. investigated the solubility between dimethyl ether (DME) and two kinds of lubricant oils, and showed the effect of molecular structures of the lubricant oils on the solubility6. The studies of the solubility as well as other thermodynamic properties have been summarized in several review papers1,2. They indicate that the solubility depends highly on the molecular structures of lubricant oils and refrigerants. Meanwhile, molecular simulations using different algorithms have become a powerful tool to investigate solubility and understand its molecular mechanism. Approximately a decade ago, Maurer et al. investigated the solubility of small gases in an ionic liquid by Gibbs ensemble Monte Carlo simulations at constant pressure and temperature (NpT-GEMC)7-9. The calculated isotherms agree well with the experimental results. Shi and Maginn also investigated the solubility of CO2 in an ionic liquid by the continuous fractional component Monte Carlo method10. The isotherms, Henry’s Law constants, and partial molar enthalpies of absorption quantitatively agree with experimental data. As for the gas solubility in polymers, Eslami and Müller-Plathe studied the solubility of several gases in polystyrene by the Grand equilibrium method, and a correct trend was obtained for the solubility coefficients of the gases with temperature11. This work has been extended to the solubility study of water in poly (ethylene terephthalate)12 and water in polyamide-6,613. Mixtures of refrigerants and lubricant oils, however, have so far not been investigated for their mutual solubility or its mechanisms by molecular simulation. In this study, we extend molecular simulation to the problem of solubility between lubricant oils and refrigerants. Of the several methods that have been proposed, we chose molecular dynamics (MD) because it is a versatile and flexible simulation framework. This is because we are interested in other thermodynamic and dynamical properties beyond solubility and want to apply MD to several refrigerant/lubricant oil mixtures in future works. In the present contribution, we first evaluate the solubility of the refrigerant carbon dioxide in


Page 4 of 36

Page 5 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the lubricant oil pentaerythrol hexanoate (PEC6) by molecular simulation. To this end, we calculate the refrigerant-oil equilibrium pressure, an important characteristic of solubility, by MD simulations. In the second step, we investigate solvation structures on a molecular scale and clarify the dominant molecular structure on the solubility.

2. THEORY AND CALCULATIONS 2.1. Calculation of the equilibrium pressure. The equilibrium pressure is the pressure at which a liquid refrigerant/lubricant oil mixture – of given temperature and composition – is in equilibrium with the gas phase above it. To calculate it, we need to find the vapor-liquid equilibrium state where the chemical potentials of the gas phase and the liquid phase coincide. Several methods have been proposed to simulate vapor-liquid equilibria by molecular simulations. One prominent method is the Gibbs ensemble Monte Carlo (GEMC) method14,15, which simulates the gas and the liquid phase separately without physical contact. Pressure and chemical potential are equated by exchanging volume and molecules between the phases, while the total volume and the total number of molecules are kept constant. Another approach is the NpT plus test particle method16-18, in which the chemical potential is calculated in NpT ensemble (number of atoms N, pressure p, and temperature T are constant). Other alternative methods are the Grand equilibrium method and the grand canonical molecular dynamics method (see for example, refs. 19 and 20). In this study, we used the test particle method, but with the NVT ensemble (number of atoms N, volume V, and temperature T are constant) rather than the NpT ensemble to avoid fluctuations of the volume. We have used this approach successfully in the past21. It is similar to the approach of Okumura and Yonezawa22, which is based on the NVT plus test particle method23, but different in that we conduct only one calculation for one temperature and one molecular ratio in the liquid phase whereas several calculations are done in the gas



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

phase. This is computationally beneficial because the calculations of the chemical potential in liquid phase usually require much computational time. The concept of this procedure is shown in Figure 1. The chemical potential of the refrigerant in the liquid mixture is approximated as a linear function of the pressure, see below. With regard to the gas phase, several states are simulated to calculate the chemical potentials as a function of the pressure, and they are fitted by a non-linear function. Where the curves cross is the equilibrium state (i.e. the chemical potentials of the liquid and gas phase are equal), and the pressure at this state point is the equilibrium pressure. Note that we assume the lubricant oil to have much lower volatility than the refrigerant, such that its presence in the gas phase can safely be ignored as is the case in the experimental studies24-26. The lower volatility of the lubricant oil is evidenced by the measurement: the vapor pressures of pentaerythritol esters between 334 and 476 K were measured to be between 5.6×10-5 and 0.94 Pa27, whereas those of CO2 are in the order of MPa.

Figure 1. Schematic of calculation of the equilibrium pressure.

The chemical potential of the refrigerant in the liquid mixture, µliq(p), is obtained from a first-order Taylor series expansion11, 16, 17, 28,


Page 6 of 36

Page 7 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


µ liq ( p ) ≈ µ liq ( p ′) +

∂µ liq ( p )

( p − p ′) ∂p = µ liq ( p ′) + vliq ( p − p ′ ) =µ

ex liq

( p′) + k BT ln (ρ liq Λ q 3

−1

)+ v

(1) liq

( p − p ′) .

Here, kB is the Boltzmann constant, ρliq the number density of the liquid phase, Λ the thermal de Broglie wavelength, q the internal partition function, and µliqex(p') the excess part of chemical potential. According to equation (1), µliq(p) can be calculated from the excess chemical potential µliqex(p') at an arbitrary pressure p' and partial molar volume vliq, which for incompressible liquids can be assumed to be constant to a good approximation over a sizeable pressure range. Therefore, the excess chemical potential of the refrigerant in the solution needs to be calculated only once. We used the Widom insertion method29 to calculate the excess chemical potential. We cannot use the same linear approximation to calculate the chemical potential of the refrigerant in the gas phase because it changes drastically as a function of pressure. The chemical potential of the refrigerant in the gas phase µgas(p) can be calculated by,

(

)

ex µ gas ( p ) = µ gas ( p ) + k BT ln ρ gas Λ3q −1 ,

(2)

where µgasex(p) is the excess chemical potential of the gas phase and ρgas is the number density of the gas phase. The calculation of the excess chemical potential was performed at several pressures by the Widom insertion method. In general, the equilibrium pressure peq can be calculated from the equilibrium state, i.e., the crossing point between eq. (1) and (2). Practically, we define µ’ as the chemical potential except the contribution of the internal partition function to calculate the equilibrium pressure: For the liquid phase,

(

)

(3)

(

)

(4)

ex ′ ( p ) ≈ µ liq ( p ′) + k BT ln ρ liq Λ3 + vliq ( p − p ′) , µ liq

and for the gas phase, ex ′ ( p ) = µ gas µ gas ( p ) + k BT ln ρ gas Λ3 .



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 36

Note that the de Broglie wavelength Λ is equal in both coexisting phases at same temperature. The calculated µ’gas were fitted by the following equation ′ ( p ) = A0 + A1 p + A2 / p + k BT ln( p / p0 ) . µgas

(5)

Here, A0, A1, and A2 are fitting parameters and p0 is the pressure unit. A functional form similar to that of eq. (5) has been used in some previous papers18,30, where the authors analytically obtained the functional form of the chemical potential of gas phases. We consider that the contribution of the internal partition function is same in both phases and the equilibrium pressures peq were calculated from the crossing point between eq. (3) and eq. (5).

2.2. Simulation method and condition For the lubricant oil, we used pentaerythritol hexanoate (PEC6), a typical POE lubricant oil. Its molecular structure, together with the atom labels used in all discussions below, is shown in Figure 2. PEC6 has four tails, and each tail has one ester group. PEC6 was modelled using all-atom OPLS31, a standard force field for organic molecules. Carbon dioxide (CO2) was used as a prototype refrigerant because it is an important natural refrigerant that has a simple molecular structure and is easy to model. The potential parameters for carbon dioxide were taken from the literature32; the semi-flexible version of the well-known EPM2 model was used. The geometric combining rule was used for the Lennard-Jones coefficients, because both OPLS and EPM2 use this rule in the original papers and we consider that it is valid for the mixture. Note that this choice of the combing rule for the mixture is not optimized, i.e. other combing rules are not tested in this study.


Page 9 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Molecular structure of the lubricant oil used in this study (pentaerythritol hexanoate, PEC6). The atom labels used in this study are indicated in parentheses.

To calculate the long-range part of the Coulomb potential, we used the reaction field method33. The reaction field dielectric constant is taken to be 3, which roughly corresponds to the typical value of POE lubricant oils. For the Lennard-Jones part of the potential, the long-range correction was applied33. The temperatures were 263 K, 283 K, 303 K, and 323 K. The cut-off radius of non-bonded potential was 1 nm, and the time step was 1 fs (1×10-15 s). All bond lengths were constrained using the SHAKE algorithm34. For the simulations of the liquid phase, we used five different mole fractions of CO2: x = 0.3, 0.5, 0.7, 0.8, and 0.9, covering the CO2 concentration range between 8 and 43 wt.%. Table 1 lists the number of PEC6 and CO2 molecules. To check the validity of the simulation, the vapor pressure of pure CO2 at 283 K was also calculated. In the pure CO2 liquid simulation, 500 CO2 molecules were used. The systems were initially equilibrated for 1 ns at NpT conditions, with the coupling times of the Berendsen thermostat and barostat being 0.2 ps and 5 ps, respectively. The target pressures were set to close to the experimental vapor pressure. The differences between the target pressures and the vapor pressures are within a few MPa. The samples were equilibrated for a further 1 ns under NVT conditions using these dimensions. Subsequently, 10 ns NVT production runs were performed to calculate the excess chemical potential. The results were divided into ten blocks, from which the averages and the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

standard deviations of the excess chemical potential were calculated. The configurations were output every 10 ps, and the insertion positions of the test particle (CO2) were determined randomly 1.6 million times per configuration. For each position, ten random orientations were tested, which means that totally 16 million insertion were done for one configuration. In our preliminary calculations, we checked the effects of the number of insertion on the accuracy of the calculation and we determined this number so that we can obtain the results with good accuracy and reasonable computational time. Note that the partial molar volume of CO2 in eq. (3) is a function of temperature, pressure, and mole fraction of the mixture and its calculation is not a trivial task. We calculated the partial molar volumes from the differences of the volumes of several NpT runs with different mole fractions over wide range of temperature and pressure. The calculated partial molar volumes are roughly in the range of 4-5×10-5 m3/mol. Hence, we used the value of 4.2×10-5 m3/mol, which is calculated under a typical condition, 303 K, 1 atm, and x ~ 0.5, throughout this study for the sake of simplicity. Note that the errors which stem from this simplification are much smaller than the typical error of the calculation of the excess chemical potential discussed below. To calculate the partial molar volume more accurately, one can use, for example, the Widom insertion method as discussed in the literatures11, 28. For the gas phase calculation, the number of CO2 molecules and the simulation cell size are listed in Table 2. The systems were equilibrated at NVT for 11 ns. The coupling time of the thermostat was chosen as 0.2 ps. A subsequent 10 ns NVT run was used to obtain the excess chemical potential. The configurations were output every 10 ps, and 100,000 insertions of the test particle were performed per configuration with ten random orientations: totally 1 million insertions were conducted. All MD simulations were carried out with the YASP package35, 36.


Page 10 of 36

Page 11 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Simulation systems for the liquid phase calculations.

Mole fraction of CO2 [x]

Number of PEC6 molecules

Number of CO2 molecules

0.3

50

21

0.5

50

50

0.7

50

117

0.8

50

200

0.9

50

450

Table 2. Simulation systems for the gas phase calculations.

Number of CO2 molecules

Simulation cell size [nm]

700

15 × 15 × 15

1400

15 × 15 × 15

830

10 × 10 × 10

150

5×5×5

200

5×5×5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3. RESULTS AND DISCUSSION 3.1. Solubility The calculated results of eqs. (3) and (4) at 283 K are shown in Figure 3. The straight lines show µ’ of CO2 in the liquid mixture calculated by eq. (3). To validate the calculation and the derivation of eq. (1), we conducted the test calculations at two different pressures at x = 0.5. The results, as plotted by the filled triangles in Figure 3, show that the line by eq. (3) is within the error of the calculations at two different pressures, indicating the validity of the derivation of eq. (1). The filled circles show µ’ in the gas phase and these points were fitted by eq. (5). The points of intersection of the lines and the fitted curve indicate the equilibrium states. The calculated equilibrium pressures as a function of the mole fraction of CO2 are collected in Figure 4, together with experimental data24. Figure 4 also shows the calculated vapor pressure of pure CO2 is 4.8 MPa (x = 1.0, T = 283 K), which agrees well with the experimental value (4.7 MPa at 285 K37). The simulation results of the mixtures (x < 1.0) also capture the characteristic behavior of the experiments: the CO2 vapor pressure becomes higher as temperature and CO2 mole fraction x increase. Moreover, the equilibrium pressures in Figure 4 agree well with the experimental data, with the deviations being slightly larger at high CO2 content. The error bars are calculated from the standard deviation of the excess chemical potential of the ten blocks described above, which can be considered as the largest part of the error. The relative deviations of the calculated equilibrium pressures from experiment are, however, +27% (T = 283 K, x = 0.7) at most and 14% on average. Considering that off-the-shelf force fields have been used without further optimization, this is encouraging for the method to be used as a predictive tool for other refrigerant-lubricant oil combinations in future research. Note that neither the experimental nor the simulated data in Figure 4 suggest a liquid-liquid phase separation at the state points studied here. A phase separation would be


Page 12 of 36

Page 13 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


evidenced by the vapor pressure being independent of CO2 concentration, i.e. by a horizontal line. Moreover, we have visually inspected snapshots of different systems in equilibrium and have never observed a sign of phase separation.

Figure 3. Calculated results of µ’ as a function of pressure. The temperature is 283 K. Filled triangles show the results of the test calculations at two different pressures.

Figure 4. Equilibrium pressure of CO2 as a function of its mole fraction. The calculated vapor pressure of pure CO2 is shown by a black filled square (x = 1.0, T = 283 K).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.2. Solvation structure To analyze the solvation structure on a molecular scale, we conducted independent 10 ns NpT runs for x = 0.5, 0.8, and 0.9. The target temperature is 303 K. The target pressures are set to close to the vapor pressure calculated in the previous section: 2.0, 5.7, and 7.6 MPa for x = 0.5, 0.8, and 0.9, respectively. Our preliminary calculation suggested that the effect of temperature and pressure is small, so we assume that we can obtain the essence of the solvation structure under this condition. The double-bonded OO atoms (for the definitions of atom labels, see Figure 2) in the carbonyl groups of PEC6 carry a negative partial charge (-0.43 e), whereas C1 atoms are positive (0.51 e). Because CO2 has a strong quadrupole moment, its interaction with the dipole of the carbonyl groups is expected to play an important role in the solvation of CO2 in this lubricant oil. Figure 5 compares the residence sites of CO2 in PEC6 with the internal structure of the lubricant oil. We show the distance distributions between the central carbon CC and other atoms in the same molecule together with the radial distribution function (RDF) between CC atoms and the carbon atom of CO2. The mole fraction of CO2 is 0.5. The distance distributions CC-OS, CC-C1, CC-OO, and CC-C2 have sharp peaks, which indicate that they form a rigid core of PEC6. The rigid core ends roughly at the ester groups: the distributions of the outer tail atoms CC-C3, CC-C4, CC-C5, and CC-C6 are broad and multi-peaked. This means that the tails beyond the ester group are flexible. Also, the first peak of the CC-CO2 RDF is clearly located just outside the CC-C1 and CC-OO distributions. This result, together with the RDF results as discussed below, implies that the CO2 molecules are likely attached to the carbonyl group, especially at the outside of them.


Page 14 of 36

Page 15 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Distance distribution functions and the radial distribution function (RDF). The distance distribution functions indicate the distance from the central carbon atom of PEC6 (CC) to the other atoms in the same molecule. The black line with markers shows RDF between CC and the carbon atom of CO2. The mole fraction of CO2 is 0.5.

Figure 6 shows RDFs in the same conditions between different lubricant atoms and carbon atoms of CO2. The RDF between OO atom and CO2 has the closest first peak among the RDFs, meaning that the CO2 molecules are strongly attached to OO atom. As shown by the vertical dashed line in Figure 6, we define the radius of the first neighbor shell from the point where the OO-CO2 RDF becomes 1 for the first time after the maximum. In this definition, we obtain r = 0.364 nm. CO2 molecules are attached to the OO atoms inside this radius. In contrast to the OO-CO2 RDF, the first peak of the CC-CO2 RDF is farthest, which indicates that CO2 is excluded from the lubricant’s core due to steric crowding, as is also evident in Figure 5. Comparing the C4-CO2 and C6-CO2 RDFs, we can see that the first peak of the C6-CO2 RDF is closer and larger than that of the C4-CO2 RDF. There are two possible reasons for this. First, C6 atoms are located at the end of the tails, so the effect of steric inhibition is relatively smaller than that of C4 atoms. Second, because of the flexibility of the tails, they are able to move into the vicinity of CO2 attached to the carbonyl group of another



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

neighboring oil molecule. For these reasons, C6 atoms have a relatively larger and closer peak than C4 atoms.

Figure 6. Radial distribution functions (RDFs) between selected PEC6 atoms and the carbon atoms of CO2. The mole fraction of CO2 is 0.5.

Figures 7(a) and 7(b) show the number of CO2 molecules within the first neighbor radius, r, from OO atoms and the number of OO atoms within r from the carbon atoms of CO2 molecules. The first neighbor radii, r, are 0.364, 0.356, and 0.346 nm for x = 0.5, 0.8, and 0.9, respectively. The molecular ratios of PEC6 and CO2 are 50:50 (x = 0.5), 50:200 (x = 0.8), and 50:450 (x = 0.9). Because each PEC6 has four carbonyl groups, the corresponding ratios of carbonyl groups to CO2 y are 200:50 = 4, 200:200 = 1, and 200:450 = 0.44, respectively. The ratio y larger than 1 corresponds to an excess of carbonyl groups but not an excess of CO2. Figure 7(a) shows that the number of CO2 molecules in the neighborhood of a carbonyl oxygen is mostly zero at x = 0.5 (y = 4), which indicates that most carbonyl groups are “free” under the condition in which the carbonyl groups are in excess. Approximately 15 % of carbonyl groups have one CO2 neighbor. To better understand this pairing, an extracted snapshot of one PEC6 molecule and its neighbor CO2 molecule within r from OO atoms is


Page 16 of 36

Page 17 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


shown in Figure 8 (a). Figure 8 (a) also presents a terminal atom of a hydrocarbon tail (C6) that is bending back toward the core. The images in Figure 8 were rendered using the software VMD38. In the case of x = 0.9 (y = 0.44), some carbonyl groups have two CO2 molecules in their first neighbor shell, but still approximately 45% carbonyl groups carry no CO2, even though there is an excess of CO2. Figure 8 (b) shows a snapshot of a PEC6 molecule that has two neighboring CO2 molecules in one tail. Figure 7 (b) shows that the number of OO atoms around CO2 decreases as CO2 concentration increases. The data at x = 0.5 (y = 4) show that over 40% of CO2 molecules are not attached to carbonyl even under carbonyl excess. This result, together with the results of x = 0.9 (y = 0.44) in Figure 7 (a), indicates that the carbonyl group and CO2 molecules are do not fully match each other, meaning that, at all concentrations, there is unspecific CO2 sorption along with CO2 residing at carbonyl sites. From the data at x = 0.8 (y = 1) in Figures 7 (a) and (b), where the numbers of the carbonyl groups and CO2 molecules are equal, we can see that the histogram shows almost the same distribution and that the number of neighbors is mostly zero or one. This implies that one OO atom and one CO2 form a pair in the first neighbor shell in the case that they are located in the vicinity.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7. (a) Number of CO2 molecules within the first neighbor radius, r, from OO atoms. (b) Number of OO atoms within r from the carbon atoms of CO2 molecules.


Page 18 of 36

Page 19 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8. Extracted snapshots of the simulations. (a) and (b) show lubricant oil molecules with one neighboring CO2 molecule and two neighboring CO2 molecules in one hydrocarbon tail, respectively. The color code is as follows: C – gray; H – white; O – red.

From these results above, we can summarize the molecular solvation structure of the lubricant oil (PEC6) and CO2 as follows. The solvation structure is schematically shown in Figure 9. (1) The tails of the lubricant oil are rigid inside the ester groups but flexible beyond. The rigid region can be defined as a core. (2) CO2 molecules barely enter the lubricant core, meaning that the core does not contribute significantly to the solubility. (3) The double-bonded oxygens of the ester groups are good sorption sites for CO2.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(4) Few CO2 molecules are simultaneously attached to two or more carbonyl oxygens. Points (3) and (4) imply that increasing the number of polar groups would increase the solubility. (5) The largest fraction of CO2 is unspecifically absorbed in the alkane tail region. This fraction increases as the CO2 concentration increases.

Figure 9. Schematic image of the solvation structure of PEC6 and CO2 on molecular scale.

4. SUMMARY We have analyzed the solubility and the solvation structure on the molecular scale between the refrigerant carbon dioxide (CO2) and a prototype lubricant oil pentaerythritol hexanoate (PEC6) by molecular dynamics. The vapor-liquid equilibrium pressures agreed well with experimental data over a wide range of temperature and CO2 mole fraction. The solvation structure has the following characteristics: First, there is a rigid core of the lubricant oil molecule that ends approximately at the ester groups. Beyond, the alkane tails are flexible. Second, the core is practically inaccessible to CO2 molecules. Third, the CO2 are preferably attached to the double-bonded oxygen atoms of the ester groups. Fourth, only a few CO2 molecules are attached to more than one carbonyl oxygen simultaneously. Finally, there is also significant unspecific sorption of CO2 in the alkane tail region, which increases as CO2


Page 20 of 36

Page 21 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


concentration increases. These results indicate that increasing the size of the rigid lubricant core would probably decrease CO2 solubility, whereas increasing the number of polar carbonyl groups would increase it.

Corresponding author *Taisuke Sugii, Phone: +81-29-353-3449, Email: [email protected]

Notes The authors declare no competing financial interests.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

REFERENCES 1. Youbi-Idrissi, M.; Bonjour, J. The Effect of Oil in Refrigeration: Current Research Issues and Critical Review of Thermodynamic Aspects. Int. J. Refrig. 2008, 31, 165-179.

2. Marsh, K. N.; Kandil, M. E. Review of Thermodynamic Properties of Refrigerants + Lubricant Oils. Fluid Phase Equil. 2002, 199, 319-334. 3. Wahlström, Å.; Vamling, L. Solubility of HFC32, HFC125, HFC134a, HFC143a, and HFC152a in a Pentaerythritol Tetrapentanoate Ester. J. Chem. Eng. Data 1999, 44, 823-828. 4. Wahlström, Å.; Vamling, L. Solubility of HFCs in Pentaerythritol Tetraalkyl Esters. J. Chem. Eng. Data 2000, 45, 97-103. 5. Yokozeki, A. Solubility of Refrigerants in Various Lubricants. Int. J. Thermophys. 2001, 22, 1057-1071. 6. Sun, Y.; Wang, X.; Gong, N.; Liu, Z. Solubility of Dimethyl Ether in Pentaerythritol Tetrahexanoate (PEC6) and in Pentaerythritol Tetraoctanoate (PEC8) Between (283.15 and 353.15) K. J. Chem. Eng. Data 2014, 59, 3791-3797. 7. Kumelan, J.; Pérez-Salado Kamps, Á.; Urukova, I.; Tuma, D.; Maurer, G. Solubility of Oxygen in the Ionic Liquid [bmim][PF6]: Experimental and Molecular Simulation Results. J. Chem. Thermodynamics 2005, 37, 595-602. 8. Kumelan, J.; Pérez-Salado Kamps, Á.; Urukova, I.; Tuma, D.; Maurer, G. Corrigendum to: “Solubility of Oxygen in the Ionic Liquid [bmim][PF6]: Experimental and Molecular Simulation Results”. J. Chem. Thermodynamics 2007, 39, 335. 9. Urukova, I.; Vorholz, J.; Maurer, G. Solubility of CO2, CO, and H2 in the Ionic Liquid


Page 22 of 36

Page 23 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[bmim][PF6] from Monte Carlo Simulations. J. Phys. Chem. B 2005, 109, 12154-12159. 10. Shi, W.; Maginn, E. J. Atomistic Simulation of the Absorption of Carbon Dioxide and

Water

in

the

Ionic

Liquid

1-n-Hexyl-3-methylimidazolium

Bis(trifluoromethylsulfonyl)imide ([hmim][Tf2N]. J. Phys. Chem. B 2008, 112, 2045-2055. 11. Eslami, H.; Müller-Plathe, F. Molecular Dynamics Simulation of Sorption of Gases in Polystyrene. Macromolecules 2007, 40, 6413-6421. 12. Eslami, H.; Müller-Plathe, F. Water Permeability of Poly(ethylene terephthalate): A Grand Canonical Ensemble Molecular Dynamics Simulation Study. J. Chem. Phys. 2009, 131, 234904. 13. Eslami, H.; Mehdipour, N. Grand Canonical Ensemble Molecular Dynamics Simulation of Water Solubility in Polyamide-6,6. Phys. Chem. Chem. Phys. 2011, 13, 669-673. 14. Panagiotopoulos, A. Z. Direct Determination of Phase Coexistence Properties of Fluids by Monte Carlo Simulation in a New Ensemble. Mol. Phys. 1987, 61, 813-826. 15. 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. Mol. Phys. 1988, 63, 527-545. 16. Möller, D.; Fischer, J. Vapour Liquid Equilibrium of a Pure Fluid from Test Particle Method in Combination with NpT Molecular Dynamics Simulations. Mol. Phys. 1990, 69, 463-473. Erratum: Möller, D.; Fischer, J. Mol. Phys. 1992, 75, 1461. 17. Lotfi, A.; Vrabec, J.; Fischer, J. Vapour Liquid Equilibria of the Lennard-Jones Fluid



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

from the NpT plus Test Particle Method. Mol. Phys. 1992, 76, 1319-1333. 18. Vrabec, J.; Fischer, J. Vapour Liquid Equilibria of Mixtures from the NpT plus Test Particle Method. Mol. Phys. 1995, 85, 781-792. 19. Guevara-Carrion, G.; Hasse, H.; Vrabec, J. In Multiscale Molecular Methods in Applied Chemistry; Kirchner, B., Vrabec, J., Eds.; Topics in Current Chemistry; Springer: Berlin, 2012; 307, pp 201-250. 20. Frenkel, D.; Smit, B. Understanding molecular simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, CA, 2002. 21. Müller-Plathe, F. Calculation of the Free Energy for Gas Absorption in Amorphous Polypropylene. Macromolecules 1991, 24, 6475-6479. 22. Okumura, H.; Yonezawa, F. Reliable Determination of the Liquid-Vapor Critical Point by the NVT plus Test Particle Method. J. Phys. Soc. Jpn. 2001, 70, 1990-1994. 23. Szalai, I.; Liszi, J.; Boda, D. The NVT plus Test Particle Method for the Determination of the Vapour-liquid Equilibria of Pure Fluids. Chem. Phys. Lett. 1995, 246, 214-220. 24. Bobbo, S.; Pernechele, F.; Fedele, L.; Stryjek, R. Solubility Measurements and Data Correlation of Carbon Dioxide in Pentaerythritol Tetrahexanoate (PEC6). J. Chem. Eng. Data 2008, 53, 2581-2585. 25. Fedele, L.; Pernechele, F.; Bobbo, S.; Scattolini, M.; Stryjek, R. Solubility of Carbon Dioxide in Pentaerythritol Tetraoctanoate. Fluid Phase Equil. 2009, 277, 55-60. 26. Sun, Y.; Wang, X.; Gong, N.; Liu, Z. Solubility of Dimethyl Ether in Pentaerythritol Tetrahexanoate (PEC6) and in Pentaerythritol Tetraoctanoate (PEC8) Between (283.15


Page 24 of 36

Page 25 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and 353.15) K. J. Chem. Eng. Data 2014, 59, 3791-3797. 27. Razzouk, A.; Mokbel, I.; García, J.; Fernandez, J.; Msakni, N.; Jose, J. Vapor Pressure Measurements in the Range 10−5 Pa to 1 Pa of Four Pentaerythritol Esters: Density and Vapor–liquid Equilibria Modeling of Ester Lubricants. Fluid Phase Equil. 2007, 260, 248-261. 28. Vrabec, J.; Hasse, H. Grand Equilibrium: Vapour-Liquid Equilibria by a New Molecular Simulation Method. Mol. Phys. 2002, 100, 3375-3383. 29. Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808-2812. 30. Vrabec, J.; Lotfi, A.; Fischer, J. Vapour Liquid Equilibria of Lennard-Jones Model Mixtures from the NpT plus Test Particle Method. Fluid Phase Equil. 1995, 112, 173-197. 31. 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. Am. Chem. Soc. 1996, 118, 11225-11236. 32. 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. 33. Barker, J. A.; Watts, R.O. Monte Carlo Studies of the Dielectric Properties of Water-like Models. Mol. Phys. 1973, 26, 789-792. 34. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

n-alkanes. J. Comput. Phys. 1977, 23, 327-341. 35. Müller-Plathe, F. YASP: A Molecular Simulation Package. Comput. Phys. Comm. 1993, 78, 77-94. 36. Tarmyshov, K. B.; Müller-Plathe, F. Parallelizing a Molecular Dynamics Algorithm on a Multiprocessor Workstation Using OpenMP. J. Chem. Inform. Model. 2005, 45, 1943-1952. 37. Haynes, W. M., Ed. CRC Handbook of Chemistry and Physics, 95th ed.; CRC Press: Boca Raton, FL, 2014. 38. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Molec. Graphics 1996, 14, 33-38.


Page 26 of 36

Page 27 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For Table of Contents Only



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1. Schematic of calculation of the equilibrium pressure. 44x24mm (600 x 600 DPI)

ACS Paragon Plus Environment

Page 28 of 36

Page 29 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Molecular structure of the lubricant oil used in this study (pentaerythritol hexanoate, PEC6). The atom labels used in this study are indicated in parentheses. 43x24mm (600 x 600 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. Calculated results of µ' as a function of pressure. The temperature is 283 K. Filled triangles show the results of the test calculations at two different pressures. 58x45mm (600 x 600 DPI)


Page 30 of 36

Page 31 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Equilibrium pressure of CO2 as a function of its mole fraction. The calculated vapor pressure of pure CO2 is shown by a black filled square (x = 1.0, T = 283 K). 62x51mm (600 x 600 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. Distance distribution functions and the radial distribution function (RDF). The distance distribution functions indicate the distance from the central carbon atom of PEC6 (CC) to the other atoms in the same molecule. The black line with markers shows RDF between CC and the carbon atom of CO2. The mole fraction of CO2 is 0.5. 43x23mm (600 x 600 DPI)


Page 32 of 36

Page 33 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Radial distribution functions (RDFs) between selected PEC6 atoms and the carbon atoms of CO2. The mole fraction of CO2 is 0.5. 51x32mm (600 x 600 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7. (a) Number of CO2 molecules within the first neighbor radius, r, from OO atoms. (b) Number of OO atoms within r from the carbon atoms of CO2 molecules. 105x174mm (600 x 600 DPI)


Page 34 of 36

Page 35 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8. Extracted snapshots of the simulations. (a) and (b) show lubricant oil molecules with one neighboring CO2 molecule and two neighboring CO2 molecules in one hydrocarbon tail, respectively. The color code is as follows: C – gray; H – white; O – red. 114x213mm (600 x 600 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9. Schematic image of the solvation structure of PEC6 and CO2 on molecular scale. 57x40mm (600 x 600 DPI)


Page 36 of 36

Molecular Dynamics Simulation of a RefrigerantâLubricant Oil System

Recommend Documents

Molecular Dynamics Simulation of a RefrigerantâLubricant Oil System