MD simulation estimation of saturation pressure and vapor-liquid equilibrium for binary blends of HFO-1234yf and HFO-1123

Saturation pressure and vapor-liquid equilibrium data of low GWP hydrofluoroolefin compound HFO-1234yf (2,3,3,3-tetrafluoropropene) and its binary blends with HFO-1123 (1,1,2-trifluoroethene) were estimated using molecular dynamics (MD) simulation with the aid of COMPASS force field from 214.15 K to 330.15 K. The maximum absolute deviations of present saturated vapor pressure data from published experimental results were 0.03% for HFO-1123 system and 0.017% for HFO-1234yf system. These models were used to estimate the saturated vapor pressure for six different mixtures of HFO-1234yf and HFO-1123. The saturated liquid and vapor densities, the vapor-liquid coexistence curve (VLCC), and the critical points of the blended refrigerants with various molecular mass fractions were predicted using MD simulations and Gibbs ensemble Monte Carlo (GEMC) simulations. The critical points of HFO-1234yf and HFO-1123 matched well with published data.


Introduction
The progress of science and engineering has upgraded the humans' thermal comfort of the lifestyle. Air-conditioners and refrigerators, which cool and heat residential spaces and preserve food, are usually manufactured with vapor compression refrigeration systems, wherein refrigerants such as chlorofluorocarbons (CFCs), hydrochlorofluorocarbons (HCFCs), and hydrofluorocarbons (HFCs) are working fluids. The use of conventional chemical refrigerants was significantly reduced after the discovery of the contribution of these chemical refrigerants to global warming and ozone layer depletion. Although HFC refrigerants do not have the potential to destroy the ozone layer, the EU recommended phasing out of them because they have high values of global warming potential (GWP).
Researchers are seeking next generation low GWP environment-friendly refrigerants. The carbon-carbon dual bonds in the molecular structure of Hydrofluoroolefin (HFO) were developed to make these environment-friendly refrigerants, which demonstrate small atmospheric lifetime and ensure low ODP and GWP values. Furthermore, mixtures of HFOs have been preferred for the development of refrigerants having lesser flammability and greater volumetric refrigerant capacities in comparison with those values of pure HFOs. The phase-out of high GWP refrigerants has facilitated the development of HFO blends. One low GWP HFO refrigerant, HFO-1234yf (2,3,3,3-tetrafluoropropene), was developed to satisfy EU regulations and is applicable in present vehicle air-conditioning systems. HFO-1234yf (GWP = 4) has a very short atmospheric lifetime of 11 days, much better than the lifetime of HFC-134a of approximately 13 years. HFO-1234yf is an excellent potential replacement for HFC-134a, requiring only a very minor modification to existing system designs. Richter et al. [1] experimentally measured vapor pressure and p − ρ − T values of HFO-1234yf (2,3,3,3-Tetrafluoroprop-1-ene). Lee and Jung [2] conducted a drop-in test of HFO-1234yf using a bench tester and demonstrated that it was possible to use this new material as a substitute for HFC-134a in mobile air conditioning systems. Zilio et al. [3] conducted experiments on HFO-1234yf with some modifications of conventional European automotive airconditioning systems that were originally designed to use HFC-134a. HFO-1234yf is a stable substance with no toxicity and only mild flammability under normal conditions. It has been classified as safety class A2L according to ASHRAE Standard 34. Akasaka et al. [4] performed thermodynamic property modeling for HFO1234yf. The saturation pressure of HFO-1234yf was measured experimentally by Nicola et al. [5].
Another HFO refrigerant, HFO-1123 (1,1,2-trifluoroethene CF 2 =CHF) has favorable thermodynamic properties, low GWP, and large latent heat; however, it is also flammable and has low volumetric capacity. Thus, mixtures of HFO-1234yf and HFO-1123 were recommended as substitutes for HFC-410A in air-conditioning systems because they can enhance cooling capacity and suppress flammability. This blended refrigerant can be used for domestic and commercial air conditioning applications. The vapor-liquid equilibrium properties of HFO-1123/ HFC-32 mixtures were evaluated by Miyamoto et al. [6]. Hashimoto et al. [7] experimentally investigated the thermophysical properties of HFO-1123/HFO-1234yf mixtures and found azeotropic-like behavior. Miyamoto et al. [8] measured experimentally the vapor-liquid equilibrium properties of binary blends of HFO-1123 and HFO-1234yf within a limited temperature range of 300.10 K to 329.95 K.
To study the performance of new refrigerants in various potential applications, with an eye toward commercialization, a wide range of property data such as the saturation pressure and vapor liquid equilibrium are necessary. The thermophysical properties of HFO-1234yf (2,3,3,3-tetrafluoropropene) [9] and vapor liquid equilibrium data for HFO-1234yf + HFC-32 and HFO-1234yf + HFC-161 blends [10] were obtained experimentally. The saturation pressures of HFO-1123 and HFO-1243zf were experimentally measured in [11]. To supplement experimental data, molecular dynamics simulations enable suitable prediction of saturation pressure and vapor liquid equilibrium properties of pure HFO compounds and their possible mixtures. The thermophysical properties of HFO compounds were evaluated for the first time by Raabe [12] using molecular simulations. Alam and Jeong [13] conducted MD simulations to compute the thermodynamic properties and critical parameters of HFO-1123 and of binary mixtures of HFO-1123 with HFC-32 and HFC-134a. Using MD simulations, they also evaluated the thermodynamic properties of HFO/HFC-448A and HFO/HFC-449A at temperatures from 233.15 K to 343.15 K [14].
Although blends of HFO-1123 and HFO-1234yf have been proposed as potential working fluids in commercial refrigeration systems, no macroscopic blend composition of HFO-1123 and HFO-1234yf has yet been suggested. Thus, thermophysical properties of blends having various blending ratio should be investigated. In this paper, using MD simulations, the saturation pressures of pure refrigerant HFO-1234yf and mixtures of HFO-1234yf + HFO-1123 in various compositions were evaluated at temperatures from 214.15 K to 330.15 K. The saturated densities, vapor liquid co-existence curves, and critical parameters of HFO-1234yf and HFO-1234yf + HFO-1123 mixtures were also predicted.

Simulation method
Six different blends were considered by changing molecule numbers and ratios of molecular mass. Combinations of molecule number, mole ratio, and ratio of molecular mass of these blends are presented in Table 1. The molecular structures of HFO-1123 and HFO-1234yf were constructed in three-dimensions and then optimized geometrically using the energy minimization principle, as shown in Fig. 1. The COMPASS II [15] force field was applied to all of the energy minimization and dynamics computations in this work. The molecular topologies for C, H, and F atoms were fitted with the aid of the COMPASS II force field. The potential function of the COMPASS II force field with modified LJ-9-6 potentials is defined as follows: The first and second terms on the right-side of Eq. (1) are the modified LJ-9-6 potential for the van der Waals term and the Coulombic function for electrostatic interaction, respectively. Here, r ij represents the distance between the ith and jth atoms; ε ij is the potential energy parameter; σ ij represents the fixed distance at which the non-bonded potential turn into zero; and q i , q j are the corresponding atomic partial charges. The remaining three terms in Eq. (1) denote bonded interactions (bond stretching, angle bending, and dihedral torsions). b 0 and θ 0 denote the equilibrium bond length and angle, respectively, and k b , k θ , k φ are the force factors. The off-diagonal LJ-9-6 parameters (ε, σ) were predicted according to a 6th-order combination law for unlike atom pairs: The electrostatic interaction was computed using partial charges of atoms. The bond-increment δ ij indicates the charge separation between valence-bonded atoms i and j needed to ensure that the charge parameters are transferable. The partial charge for atom i is the sum of the entire charge bond increment, so that The Amorphous cell and Forcite modules of BIOVIA Materials Studio 7.1 [16] were applied for simulation and result analysis. The atom types and force field parameters are given in Table 2. Force field parameters relating to both bonded terms and LJ-9-6 terms are transferable for blend components.
Eight individual cubic cells with periodic boundary conditions in all directions were built for pure HFO-1234yf, pure HFO-1123, and the six binary blends of HFO-1234yf and HFO-1123. The molecule numbers in these cells were varied from 284 to 297 to meet the blending ratios. Cells enclosing 285 HFO-1123 vapor molecules with density of 0.0472 g cm −3 and lattice constant of 93.70 Å are presented in Fig. 2. Eleven distinct PBC cells were built with initial vapor phase densities determined according to the eleven different temperature values from 214.15 K to 330.15 K for each system. These cells were optimized geometrically using the Smart minimizer algorithm; then, Materials Studio Forcite module was used to perform MD-NVT simulations for estimation of the saturated vapor pressures. Simulations were performed for 2.5 ns to equilibrate the cells; the following 2.5 ns were a production simulation. The time step was set to 1 fs. A Nosé-Hoover-Langevin (NHL) thermostat [17] was used with an effective relaxation time of 0.1 ps to maintain the temperature of each cell. The van der Waals interactions and electrostatic interactions were calculated via Groupbased technique and Ewald summation technique, respectively. The cut-off radius of 12.5 Å with L-J tail corrections  [17] and a Berendsen barostat [18], with corresponding effective relaxation times of 0.1 ps. and 1 ps., were applied to maintain the temperature and pressure of these cells.
In the present work, using the Materials Studio EQUI-LIBRIA package, Gibbs ensemble Monte Carlo (GEMC) [19] simulations were performed to estimate the VLCC with critical temperatures and critical densities. From previous MD-NPT simulations, the eight saturated liquid and vapor densities computed at 273.15 K, 283.15 K, 290.15 K, and 300.15 K were fitted using the EQUILIB-RIA module to predict the VLCC and the critical points of each pure and binary mixture. The critical points were predicted based on the scaling law [19] and the law of rectilinear diameters, where A and B are constants and β is the scaling exponent. The scaling exponent was set to 0.30, and there were 20 coexistence points computed.

Results and discussion
Molecular dynamics simulations were performed with the aid of the COMPASS II force field model to compute the vapor pressures and saturated densities of HFO-1123 and HFO-1234yf refrigerants. Table 3 shows these simulation results. These values were compared with previous MD simulation values of Raabe [12] and experimentally measured data of Higashi et al. [11] for temperatures from 214.15 K to 330 K. On the other hand, the vapor pressures of HFO-1234yf obtained by the present MD simulation were compared with the experimental data of Richter et al. [1]. All these data are plotted in Fig. 3a. The comparison shows excellent agreement. The present MD simulation model yielded consistent estimations of the vapor pressures for HFO-1123 and HFO-1234yf; these values are believed to be useful for computation of the thermodynamic properties of the refrigerant blends. The COMPASS II force field parameters with regard to both the intramolecular terms and the LJ parameters are also transferable for mixture components. Moreover, our previous MD simulations [13,14,[20][21][22] using these force field potentials predicted the saturation vapor pressure and vapor-liquid equilibrium properties well for pure and blend components. Figure 3a shows that the saturation pressures of both HFO-1123 and HFO-1234yf increased with temperature, and the saturation pressure of HFO-1123 was significantly greater than that of HFO-1234yf. Pressure deviations of HFO-1123 and HFO-1234yf in the present data, experimental data [11], EOS [1] data, and Raabe [12] simulated data are plotted in Fig. 3b to clearly show the differences.
The COMPASS II force field model was extensively used to run MD simulations of six binary blends A−F of HFO-1123 and HFO-1234yf with various compositions of molecule mass and mole fractions. The computed vapor pressures of HFO-1123, HFO-1234yf, and their six binary blends A−F at temperatures from 214.15 K to 330.15 K are presented in Table 4 and plotted in Fig. 4. It can be seen in Fig. 4 that the vapor pressure curves varied consistently among HFO-1123, the binary blends, and  Table 5. The densities of HFO-1123 estimated by Raabe [12] and the densities of HFO-1234yf from REFPROP [23] are also included in this table for comparison with the present work. Table 5 shows that the saturated density of HFO-1123 predicted by the present MD simulation was slightly larger than that in the estimation by Raabe [12]. The density of HFO-1234yf estimated by the present MD simulation was also higher than that in the REFPROP data. Table 5 shows that saturated liquid density of pure HFO-1234yf was slightly larger than that of pure HFO-1123, though saturated vapor density of pure HFO-1234yf was smaller than that of pure HFO-1123. The density of the saturated vapor decreased significantly as the content of HFO-1234yf increased, while the density of the saturated liquid changed only slightly. It should be noted, however, that the density of the saturated liquid increased when the molecular mass fraction of HFO-1234yf increased from 0 (pure HFO-1123) to 10% (blends A ), and decreased to approach the value of pure HFO-1234yf as the content of HFO-1234yf increased. This trend is in line with the trend in vapor pressure that can be observed in Fig. 4, in which the vapor pressure significantly decreased as the molecular mass fraction of HFO-1234yf increased from 0 (pure HFO-1123) to 10% (blend A ), and then decreased slightly as the content of HFO-1234yf increased further. This trend is thought to be due to interactions between molecules resulting from the mixing of two pure substances. When two pure components are mixed at various composition ratios, inter-and intra-molecular interactions occur between pure components and a new mixed compound with new molecular configuration is formed. The individual thermodynamic properties of the two pure compounds mutate in these new mixed compounds and, finally, the blended compounds exhibit modified azeotropic-like properties. Vapor-liquid coexistence curves (VLCC) of pure HFO-1234yf and of blends A−F were obtained by conducting GEMC simulations with twenty data points; the curves are presented in Figs. 5, 6, 7, 8, 9, 10 and 11. The saturated vapor and liquid densities, as computed from the MD simulation data shown in Table 5, are also depicted in Figs. 5, 6, 7, 8, 9, 10 and 11 for comparison between the present MD and GEMC simulations. In Fig. 5, the VLCC of pure HFO-1234yf is presented and compared with the phase diagram obtained by REFPROP. Figure 5 demonstrates a good agreement of the vapor densities between values computed from the present GEMC simulations and the REFPROP data, while the computed liquid densities are slightly larger than those of the REFPROP data. The saturated liquid Table 3 Saturated pressures (P sat ) of HFO-1123 and HFO-1234yf from present MD simulations in comparison with MD simulated data from Raabe [12], and experimental data from Higashi et al. [11] and Richter et al. [  and the REFPROP version 9.1 values for HFO-1234yf, as shown in Fig. 5. The temperatures and densities at the critical points for HFO-1123 and HFO-1234yf were computed using the present GEMC simulation and are presented in Table 6. These critical values are also presented in Fig. 5 for HFO-1234yf. The critical temperature and density of HFO-1123 obtained by the present GEMC simulation agreed well with the experimental data of Higashi et al. [11] and the simulated values of Raabe [12]. In this regard, it was supposed that the present force field model can be applied to calculate the VLCC of blends A−F . Table 6 displays the critical points of HFO-1123, HFO-1234yf, and blends A−F as estimated by the present GEMC simulations. The EQUILIBRIA module was applied to estimate the critical points shown in Table 6. Figure 6 shows values of VLCC computed by GEMC simulation and the saturated densities obtained from MD simulations of blend A . The critical temperature value and the critical density from the present GEMC simulations of blend A are presented in Fig. 6 and Table 6. Figures 7, 8, 9, 10 and 11 show VLCC values predicted by GEMC simulations  and the saturated densities of the MD simulations for blend B , blend C , blend D , blend E , and blend F , respectively. The critical temperature values and the critical densities from the present GEMC simulations of blend B , blend C , blend D , blend E , and blend F are also presented in Figs. 7, 8, 9, 10, 11 and Table 6. Figure 12 shows the combined VLCC curves of pure HFO-1123, blends A−F , and pure HFO-1234yf obtained by the present GEMC simulations.
This figure indicates that the critical points and both vapor and liquid densities of blends A−F were predicted appropriately, with appropriate deviation between HFO-1234yf and HFO-1123. Moreover, it should be noticed that the critical temperature increased as the molecule mass fraction of HFO-1234yf increased from blend A to blend F ; conversely, the critical densities decreased as the molecular mass fraction of HFO-1234yf increased. The   critical point of pure HFO-1123 matched well with the experimental data and the critical point of pure HFO-1234yf agreed with the REFPROP data, with only minor deviations, as shown in Table 6 and Fig. 12. A uniform and consistent variation among the critical points of blends A−F can be noticed in Fig. 12. It should be noted that the saturated liquid densities for the HFO-1234yf and HFO-1123 mixture were higher than those for pure HFO-1123 and pure HFO-1234yf in a certain range of temperature. This observation is not common but particular for this blend. The VLCCs of blends A−F that we obtained in this research can be used to understand the nature of multi-compound systems and to investigate whether additional factors specific to multi-compound systems need to be determined. The internal energy and enthalpy of HFO-1123, HFO-1234yf, and blends A−F , calculated from the present MD-NVT simulations, are given in Table 7 in the Appendix. These MD simulation data are expected to supplement the experimental data, as there are no existing data on the thermodynamic properties of blends of HFO-1123 and HFO-1234yf.  The saturated vapor-liquid densities of the six binary blends were computed at various temperatures. GEMC simulations were conducted to estimate the vapor-liquid coexistence curves and critical parameters of pure refrigerants, as well as of their blends. The critical temperature increased as the content of HFO-1234yf increased, while the critical density decreased. The vapor-liquid coexistence curves (VLCC) of the pure refrigerants and their mixtures were obtained by conducting GEMC simulations. VLCC of the binary blends obtained in this research can be used to understand the nature of multicompound systems.
In addition, internal energy, enthalpy, and density of various blends of HFO-1123 and HFO-1234yf are presented in tabular forms. These data are expected to supplement the experimental data because there are currently no published data. The saturation pressure and vapor liquid equilibrium properties of the binary mixtures of HFO-1123 and HFO-1234yf will be valuable for investigating the performance of refrigeration systems working with the new mixture refrigerants.