In this work, new thermodynamic models for refrigerant mixtures are provided for the binary pairs R-1234yf/134a, R-1234yf/1234ze(E), and R-134a/1234ze(E) based on new reference measurements of speed of sound, density, and bubble-point pressures. Fitting the very accurate liquid-phase speed of sound and density data reproduces the bubble-point pressures to within close to their uncertainty, yielding deviations in density less than 0.1% and speed of sound deviations less than 1% (and less than 0.1% for R-1234yf/134a). Models are also presented for the binary pairs R-125/1234yf, R-1234ze(E)/227ea, and R-1234yf/152a based solely on bubble-point measurements.

## 1. Introduction

The fourth generation of refrigerants largely comprises mixtures containing halogenated olefins,^{1} for instance, to find non-flammable blends to replace the workhorse refrigerant R-134a.^{2} In order to reliably assess novel mixtures, new mixture models in the form of equations of state (EOS) are needed for each of the binary pairs included in the candidate mixtures. Of the six binary mixtures we selected for study, current models are all based upon estimation schemes or unpublished models in REFPROP 10.0.^{3,4} New reference experimental measurements are needed to improve the mixture models. As outlined below, these measurements were obtained as part of the larger project, forming the foundation for the mixture models. This work begins the process of filling in the holes where important mixture models need to be updated or developed.^{3}

The components considered in this work are all relatively similar and some information about them is presented in Table 1. Two [R-1234yf and R-1234ze(E)] are stereoisomers, and the three main components of the comprehensive measurements have a hydrocarbon backbone with the replacement of four hydrogens by fluorines (as indicated by 4 as the right-most numerical digit). Therefore, we might expect (as we see here) that mixtures of these compounds would behave in a relatively simple way, with nearly ideal mixing behaviors. This makes refrigerant mixtures such as these ideal candidates for the development of highly accurate mixture models.

Name . | Hash^{a}
. | T_{NBP} (K)
. | T_{crit} (K)
. | p_{crit} (MPa)
. | M (kg mol^{−1})
. |
---|---|---|---|---|---|

R-125^{5} | 25c5a3a0 | 225.060 | 339.173 | 3.617 70 | 0.120 021 |

R-1234yf^{6} | 40377b40 | 243.692 | 367.850 | 3.384 40 | 0.114 042 |

R-134a^{7} | ff1c0560 | 247.076 | 374.210 | 4.059 28 | 0.102 032 |

R-152a^{8} | 63f364b0 | 249.127 | 386.411 | 4.516 75 | 0.066 051 |

R-1234ze(E)^{9} | 9905ef70 | 254.177 | 382.513 | 3.634 90 | 0.114 042 |

R-227ea^{10} | 40091ee0 | 256.810 | 374.900 | 2.925 00 | 0.170 029 |

Name . | Hash^{a}
. | T_{NBP} (K)
. | T_{crit} (K)
. | p_{crit} (MPa)
. | M (kg mol^{−1})
. |
---|---|---|---|---|---|

R-125^{5} | 25c5a3a0 | 225.060 | 339.173 | 3.617 70 | 0.120 021 |

R-1234yf^{6} | 40377b40 | 243.692 | 367.850 | 3.384 40 | 0.114 042 |

R-134a^{7} | ff1c0560 | 247.076 | 374.210 | 4.059 28 | 0.102 032 |

R-152a^{8} | 63f364b0 | 249.127 | 386.411 | 4.516 75 | 0.066 051 |

R-1234ze(E)^{9} | 9905ef70 | 254.177 | 382.513 | 3.634 90 | 0.114 042 |

R-227ea^{10} | 40091ee0 | 256.810 | 374.900 | 2.925 00 | 0.170 029 |

^{a}

The hash used in REFPROP to define the mixture models is generated from the SHA256 hash of the standard InChI key from hashlib.sha256(StdInChIkey).hexdigest()[2:9] + “0” with the hashlib module from the Python standard library.

## 2. Data

### 2.1. New measurement data

As part of the larger scope of work within this project, measurements of speed of sound (SOS) (liquid phase), density (gas and liquid phase), and vapor–liquid equilibria (VLE) were carried out for the binary mixtures R-1234yf/134a, R-1234yf/1234ze(E), and R-134a/1234ze(E). These SOS and bubble-point results are presented in a set of papers in the literature.^{11,12} The datasets are presented in Table 2 and shown graphically in Fig. 1. Most of our measurements extend up to 10 or 20 MPa, with the exception of the SOS data for R-134a/1234ze(E), which extend up to 50 MPa. The interim report from the project includes the full set of measured data,^{13} though additional analysis and screening of the data was carried out following report preparation. For three additional mixtures [R-125/1234yf, R-1234ze(E)/227ea, and R-1234yf/152a], measurements of bubble points were carried out.^{11} These datasets form the core of data used in this study because their experimental uncertainties are small and carefully assessed. Other data from the literature are compared with the models developed in this work.

Pair (1/2) . | Kind . | T (K)
. | z_{1} (mole frac.)
. | N
. | $U(\chi )\u0304$ (%) . |
---|---|---|---|---|---|

R-1234yf/1234ze(E) | VLE | 270–360 | 0.324, 0.638 | 48 | 0.219 |

R-1234yf/134a | VLE | 270–360 | 0.319 9, 0.646 7 | 39 | 0.169 |

R-134a/1234ze(E) | VLE | 270–360 | 0.334 1, 0.663 1 | 24 | 0.17 |

R-1234yf/1234ze(E) | PVT | 230–400 | 0.335 84, 0.666 6 | 225 | 0.0357 |

R-1234yf/134a | PVT | 230–400 | 0.336 34, 0.667 59 | 226 | 0.0375 |

R-134a/1234ze(E) | PVT | 230–400 | 0.332 5, 0.663 56 | 175 | 0.042 |

R-1234yf/1234ze(E) | SOS | 230–345 | 0.335 84, 0.666 6 | 131 | 0.0801 |

R-1234yf/134a | SOS | 230–345 | 0.336 34, 0.667 59 | 118 | 0.0793 |

R-134a/1234ze(E) | SOS | 230–345 | 0.329 16, 0.671 02 | 304 | 0.0624 |

Pair (1/2) . | Kind . | T (K)
. | z_{1} (mole frac.)
. | N
. | $U(\chi )\u0304$ (%) . |
---|---|---|---|---|---|

R-1234yf/1234ze(E) | VLE | 270–360 | 0.324, 0.638 | 48 | 0.219 |

R-1234yf/134a | VLE | 270–360 | 0.319 9, 0.646 7 | 39 | 0.169 |

R-134a/1234ze(E) | VLE | 270–360 | 0.334 1, 0.663 1 | 24 | 0.17 |

R-1234yf/1234ze(E) | PVT | 230–400 | 0.335 84, 0.666 6 | 225 | 0.0357 |

R-1234yf/134a | PVT | 230–400 | 0.336 34, 0.667 59 | 226 | 0.0375 |

R-134a/1234ze(E) | PVT | 230–400 | 0.332 5, 0.663 56 | 175 | 0.042 |

R-1234yf/1234ze(E) | SOS | 230–345 | 0.335 84, 0.666 6 | 131 | 0.0801 |

R-1234yf/134a | SOS | 230–345 | 0.336 34, 0.667 59 | 118 | 0.0793 |

R-134a/1234ze(E) | SOS | 230–345 | 0.329 16, 0.671 02 | 304 | 0.0624 |

### 2.2. Existing literature data

The collection of experimental data was initially based upon the survey of Bell *et al.*,^{3} followed by the addition of one dataset from the work of Tomassetti.^{14} All the datasets are included in the SOURCE database, accessible through NIST TDE,^{15} and are listed in Tables 7 and 8. Aside from the sources listed, there is additionally one SOS dataset,^{16} two references reporting critical loci,^{17,18} and one reporting specific heats.^{19}

### 2.3. Error metrics

In this work, the relative deviation in an arbitrary quantity *χ* is defined by

where the subscript TW indicates the value obtained from the model developed in this work and the subscript exp indicates the experimental value. The average absolute relative deviation (AAD) in a quantity *χ* is defined by

where $\Delta \u20d7%(\chi )$ is the vector of deviations and AAD is therefore on a percentage basis.

## 3. Modeling

### 3.1. Multi-fluid model

The most accurate thermodynamic mixture models available today in the reference software libraries (NIST REFPROP,^{4} CoolProp,^{20} and TREND^{21}) are based upon the Helmholtz energy. Derivatives of the Helmholtz energy can be used to obtain all other thermodynamic properties. The Helmholtz energy [divided by the molar gas constant *R* and the temperature; *α* = *a*/(*RT*)] is given as the sum of residual and ideal gas contributions,

The ideal gas portion does not enter into the fitting and is therefore not further discussed here. The molar gas constant *R*, which is now an exactly defined value according to CODATA,^{22} is instead implemented as the mole-fraction-weighted average of the molar gas constants used in developing the EOS for the pure fluids.

The mixture model for the residual portion is given by

where $\alpha CSr$ is the corresponding-state contribution given by

where **z** is the vector of mole fractions. The pure-fluid contributions $\alpha 0,ir$ are given at the mixture reduced states *τ* and *δ*. The reduced density *δ* = *ρ*/*ρ*_{red}(**z**) and reciprocal reduced temperature *τ* = *T*_{red}(**z**)/*T* are defined based on the reducing functions given in a common form by

where *Y* is the parameter of interest, either molar volume *v* = 1/*ρ* or temperature *T*. The necessary parameters are given by

The departure contribution in Eq. (4) is given by

where *F*_{ij} is the scaling factor, normally equal to 1.0 if a departure term has been fit and zero if not.

The departure function $\alpha ijr(\tau ,\delta )$ for the binary pair *ij* is, in principle, an arbitrary mathematical function that takes the value of zero at zero density, yields a good representation of the experimental data for the mixture, and has reasonable extrapolation behavior.

As summarized in Ref. 23, the standard thermodynamic properties may be expressed in terms of derivatives of the Helmholtz energy with a concise derivative representation given by

where $*$ is one of ig (ideal gas), r (residual), or tot (total). Thermodynamic properties can be given by combinations of derivatives of the Helmholtz energy. For instance, pressure is given by residual contributions only,

and the SOS *w* is obtained from

which contains mostly residual properties, except for $\Lambda 20(ig)$ that comes from the ideal gas. The quantity *R* is the molar gas constant, *M* is the molar mass, and all quantities on the right-hand side are non-dimensional.

### 3.2. Implementation

In the course of this project and in the development of routines for tracing critical curves,^{24} it became clear that a novel and more flexible approach was needed to obtain the thermodynamic derivatives of a mathematical model. The design constraints are that the obtained values should be fast to evaluate but also allow for rapid prototyping of new modeling approaches. Thus, the teqp library was established;^{23} teqp uses automatic differentiation to obtain the derivatives needed to calculate thermodynamic properties without any hand-written derivatives, dramatically speeding up the development process. Even though automatic differentiation is used, the obtained values are in many cases still faster to evaluate than the hand-written derivatives in REFPROP while suffering only a negligible loss in numerical precision. As of publication, teqp does not include any iterative routines, rather the focus is on forward derivatives of the residual Helmholtz energy with respect to temperature, density, and compositions, which meets the needs of this optimization campaign.

## 4. Model Optimization

The fitting approach for the binary mixtures with comprehensive measurements considered only SOS and densimetry data (no phase equilibria data). SOS and densimetry data are well-suited to optimization because they do not require a full phase equilibrium calculation, the phase equilibrium calculation being an unreliable (and slow) mixture calculation in general. The overall cost function *C*_{$} was defined based upon differences between model predictions and the experimental data and is given by

where the weighting factor *W* = 0.25 is used to balance deviations between the two kinds of data and the residua *r*_{ρ,i} and *r*_{w,i} are relative deviations.

### 4.1. Density

For density data, it would be ideal to use the residuum

but a downside of this approach is that the calculation *ρ*_{model}(*T*, *p*, **z**) requires computationally costly iteration (which also can fail if the initial guess is not accurate enough). Instead, we would like a non-iterative proxy that is a good stand-in for the iterative calculation. Starting from an isothermal series expansion of molar density around the experimental density yields

Truncation of higher-order terms yields the difference in density of

and the proxy residue for *p*-*ρ*-*T* data is therefore

in which *p*(*T*, *ρ*) and the isothermal derivative of *p* with respect to *ρ* are evaluated simultaneously (and directly) from the equation of state. This cost function contribution has a similar behavior to Eq. (14), though it is non-iterative and, thus, very fast to evaluate (on the order of a few μs/evaluation in C++).

Data points along the curve where $\u2202p/\u2202\rho T$ approaches zero (it is zero at the critical point for a pure species) result in large contributions to the cost function, but most of the data points are relatively far from this locus.

### 4.2. Speed of sound

The multi-fluid Helmholtz-explicit model has as independent variables temperature *T*, density *ρ*, and the vector of mole fractions **z**. Unfortunately, SOS is usually measured quasi-isochorically with temperature, pressure, and composition as independent variables, so internal iteration is required to solve for the mixture density, given the temperature, pressure, and composition (what was attempted to be avoided in the density deviation term). In the case of SOS data, the initial guess density was that obtained from the model in REFPROP 10.0.^{4} The deviations between the guess density and the final density were very small, small enough to serve as a reliable-enough starting value in any case.

Once the density has been obtained from one step of iteration, the SOS *w* is obtained from Eq. (12). In teqp , $\Lambda 01r$ and $\Lambda 02r$ are obtained from a single call and $\Lambda 11r$ and $\Lambda 20r$ are obtained from two further calls. The ideal gas contribution $\Lambda 20(ig)$ is not implemented in teqp as of publication, and so this quantity was taken from the implementation in the HEOS backend of CoolProp.^{20} The quantity $\Lambda 20(ig)$ does not depend on the departure function or the interaction parameters.

The SOS residue is therefore defined by

### 4.3. Optimization

As has been successfully applied by the author previously,^{25} stochastic (random) global optimization provides a reliable approach that is robust to failures and computationally efficient enough for practical application.

While experiments with new model formulations were carried out (for instance, the invariant reducing functions from the GERG-2004 model^{26} were added to teqp), the end goal was always to develop thermodynamic models compatible with REFPROP 10.0.^{4} The full set of departure functions available in REFPROP 10 is relatively narrow, with terms of the following kinds available:

$nk\tau tk\delta dk$,

$nk\tau tk\delta dk\u2061exp(\u2212\delta lk)$ (see Ref. 27),

$nk\tau tk\delta dk\u2061exp(\u2212\eta k(\delta \u2212\epsilon k)2\u2212\beta k(\delta \u2212\gamma k))$ (see Eq. 7.8 from Ref. 26), and

$nk\tau tk\delta dk\u2061exp(\u2212\eta k(\delta \u2212\epsilon k)2\u2212\beta k(\tau \u2212\gamma k)2)$ (Gaussian terms like those used for pure fluids).

The parameters *l*_{k} and *d*_{k} must be integers in order to yield finite contributions to the second (and higher) virial coefficients. Thus, optimization of the departure function is a mixed-integer optimization problem in which the values of some parameters are integers and others are floating point values. The conventional approach for this problem is to “fit” the *d*_{k} values by manual optimization while allowing the optimizer to fit the floating point values. A similar approach was taken here. The values of *d*_{k} were set to *d*_{k} = *k* + 1 for *k* starting at 0 and the same for *l*_{k}: *l*_{k} = *k* + 1.

The cost function defined above is minimized with global optimization. The variables *β*_{T,ij}, *γ*_{T,ij}, *β*_{v,ij}, and *γ*_{v,ij} were all bounded to be in [0.75, 1.25], *n*_{k} was bounded to be in [−3, 3], and *t*_{k} was bounded to be in [0, 4]. The CEGO library^{28} was specifically designed for this particular problem, but the widely available differential evolution algorithm implemented in scipy^{29,30} proved to be adequate for the optimization task and was used for the optimization.

The departure term used was very simple,

where the sgn function is the sign of the value, zero for an argument of zero and 1 for positive arguments.

In a single evaluation of the cost function, first a Python data structure is constructed that converts an array of double-precision numbers (the variables in the optimization) to the necessary values of *β*_{T,ij}, *γ*_{T,ij}, **n**, **d**, and so on in Eqs. (19) and (6). This Python data structure is dumped to a string in the JSON format, which is then unpacked in the mutant construction routines of teqp. One important implementation note is that the mutant models only hold *references* to the pure fluid EOS because construction of the pure fluid EOS is the most costly part of the mixture model construction. Thus, users should be careful to ensure that the pure fluid models do not fall out of scope and get prematurely de-allocated as this will result in a dangling reference to the pure fluid EOS. Construction of the mutant usually takes on the order of tens of microseconds.

## 5. Comprehensive Results

The first set of model results are those for the binary pairs R-1234yf/134a, R-1234yf/1234ze(E), and R-134a/1234ze(E). For each of these binary pairs, sufficient data were available to fit a reducing function in the form of Eq. (6) and a small departure function in the form of Eq. (19). In order to avoid model over-fitting, a graduated optimization approach was taken for each mixture model. The departure function defined above was given an increasing number of terms, and the statistics of fit were monitored in order to determine how many terms would be required. In all cases, the variables *β*_{T,ij}, *γ*_{T,ij}, *β*_{v,ij}, and *γ*_{v,ij} were fitted. This process was done in a deterministic fashion. For each of the three binary pairs with the full complement of data, 0 to *N*_{dep} terms were added to the departure function for the pair. For each term count in the departure term, the optimization was repeated five times to hopefully find the global minimum of the cost function (although that cannot be *guaranteed* in general).

The optimized cost function as a function of the number of terms in the departure function is shown in Fig. 2 for each of the binary pairs. The cost function values for *N*_{dep} = 0 are obtained by optimizing only *β*_{T,ij}, *γ*_{T,ij}, *β*_{v,ij}, and *γ*_{v,ij}. It is difficult to make out at the scale of the figure without zooming in, but all five replicates of the optimization result are shown, highlighting that the method quite reliably finds close to the same minimum of the cost function. There is a relatively large step decrease when adding a single term to the departure function. This result demonstrates that adjusting the four interaction parameters alone is not sufficient to obtain a good representation of the data; the departure function is needed. The reduction in cost function slows down after more than two terms are included, so the (somewhat arbitrary) decision was made to use two terms in the departure function.

The obtained values of the interaction parameters are shown in Table 3, and the departure functions are given in Tables 4–6. The optimization result with two terms in the departure function appears to be a good compromise of flexibility and model fidelity.

Pair (1/2) . | β_{T,ij}
. | γ_{T,ij}
. | β_{v,ij}
. | γ_{v,ij}
. | F_{ij}
. |
---|---|---|---|---|---|

R-1234yf/1234ze(E) | 0.998 886 | 0.993 309 | 0.999 302 | 0.998 590 | 1.0 |

R-1234yf/134a | 1.000 026 | 0.987 057 | 1.000 272 | 1.003 747 | 1.0 |

R-134a/1234ze(E) | 0.998 593 | 0.992 009 | 0.998 995 | 0.998 621 | 1.0 |

Pair (1/2) . | β_{T,ij}
. | γ_{T,ij}
. | β_{v,ij}
. | γ_{v,ij}
. | F_{ij}
. |
---|---|---|---|---|---|

R-1234yf/1234ze(E) | 0.998 886 | 0.993 309 | 0.999 302 | 0.998 590 | 1.0 |

R-1234yf/134a | 1.000 026 | 0.987 057 | 1.000 272 | 1.003 747 | 1.0 |

R-134a/1234ze(E) | 0.998 593 | 0.992 009 | 0.998 995 | 0.998 621 | 1.0 |

k
. | n
. | t
. | d
. | l
. |
---|---|---|---|---|

0 | 0.051 900 | 2.477 314 | 1 | 1 |

1 | −0.011 472 | 0.070 541 | 2 | 2 |

k
. | n
. | t
. | d
. | l
. |
---|---|---|---|---|

0 | 0.051 900 | 2.477 314 | 1 | 1 |

1 | −0.011 472 | 0.070 541 | 2 | 2 |

k
. | n
. | t
. | d
. | l
. |
---|---|---|---|---|

0 | 0.072 640 | 0.012 643 | 1 | 1 |

1 | −0.024 746 | 3.992 829 | 2 | 2 |

k
. | n
. | t
. | d
. | l
. |
---|---|---|---|---|

0 | 0.072 640 | 0.012 643 | 1 | 1 |

1 | −0.024 746 | 3.992 829 | 2 | 2 |

### 5.1. NIST results

First, we consider deviations in the SOS data in Fig. 3. For the mixture R-1234yf/134a, the SOS deviations are mostly significantly below 0.1%. To be sure, some of the quality of fit can be attributed to the symmetry of the mixture interactions, as partially evidenced by the near perfect overlaying of their vapor pressure curves. For the two binary pairs containing R-1234ze(E), the deviations are significantly larger and appear to, in general, increase as the amount of R-1234ze(E) increases. Model deviations are only meaningful when compared with the measurement uncertainty, and in this case we have experimental uncertainties that have been carefully assessed. Figure 4 shows the model deviations compared with the combined relative uncertainties of the experimental measurements at each state point. The model deviations for R-1234yf/134a are mostly within two times the combined expanded relative uncertainties, even in the critical region, which is excellent agreement, especially considering that there is still uncertainty in the pure fluid equation of state’s representation of the SOS. As before, the SOS deviations for the binary mixtures containing R-1234ze(E) are much larger than the experimental uncertainty in the binary mixture measurements.

In order to understand the deviations for R-134a/1234ze(E) and R-1234yf/1234ze(E), the SOS data were collected for pure R-1234ze(E). In total, three datasets were used to fit the pure fluid EOS of Thol and Lemmon:^{9} Lago *et al.*^{31} (liquid phase), Perkins and McLinden^{32} (vapor phase), and Kano *et al.*^{33} (vapor phase). A recent liquid-phase dataset is available from our group.^{34} Figure 5 shows the deviations for the liquid-phase data. For the data obtained in McLinden and Perkins,^{34} the path length obtained from direct measurement of the spacer length at atmospheric pressure and 293 K agreed with that obtained from SOS measurements with propane to within 0.03% at the same conditions, giving confidence in the obtained values. These two datasets show that while the EoS fits the single liquid-phase dataset available at the time (except for the highest temperature of 360 K), the new data from our group contradict the data from Lago *et al.*^{31} Indeed, the same problem is seen with the measurements taken from Ref. 31 for R-1234yf; their R-1234yf measurements deviate by more than 2.5% from the new EoS for R-1234yf,^{6} calling into question their liquid-phase SOS data for pure R-1234ze(E) that appear to deviate systematically.

Figure 6 shows the deviations in density for the three mixtures. Away from the critical region (see Sec. 5.2.3), the deviations are nearly all less than 0.1%, with AAD all significantly less than 0.1%. The pure-fluid densities are in general not better reproduced than 0.1% by the pure-fluid EOS, so this is about the best that can be achieved without overfitting. Just like for SOS, deviations relative to the uncertainty are considered, and shown in Fig. 7. Most of the density data are represented within five times the experimental uncertainty, except for in the critical region. The density deviations appear to have less sensitivity to defects in the pure-fluid EoS for R-1234ze(E) than for SOS.

Finally, Fig. 8 shows the bubble-point deviations. Almost all of the data are represented within 1%, except for two data points for R-134a/1234ze(E) for the composition more rich in R-1234ze(E). The deviations divided by the respective combined expanded uncertainties at each individual state point are shown in Fig. 9. The deviations in pressure are mostly within two times the combined expanded uncertainties, again except for the mixture R-134a/1234ze(E). These deviations are remarkable because the bubble-point pressures were not included in the optimization procedure; only the SOS and density data were included.

Overall, the set of new measurements appear to be very consistent, as evidenced by the fact that fitting the SOS and densimetry data allows the vapor–liquid-equilibrium data to be represented to within close to their experimental uncertainties. From a certain standpoint, this is not terribly surprising, as a brief thought experiment will indicate. Let us suppose that we can measure the pressure of a pure fluid along a subcritical isotherm as a function of density, including the unstable portion of the isotherm between the spinodals. This pressure curve is all that is needed to obtain the vapor pressure after applying the Maxwell conditions for phase equilibrium. Therefore, for a pure fluid, perfect knowledge of the density along an isotherm would yield the correct vapor pressure by default. The extension of this argument to mixtures is somewhat more challenging because the energetic portion of the phase equilibrium is no longer defined by the lever rule (alternatively, equality of the Gibbs energy); rather, it is necessary to equate chemical potentials of each component in each phase.

### 5.2. Other literature data

Aside from the data from our group (see above), there are a few other datasets in the literature for the properties we measured. Some significant discrepancies can be identified. In each case, our VLE data (not included in the fitting) agrees with one of the datasets, but not the other one. The datasets are listed in Tables 7 and 8, and the statistics of the fit are also included.

Pair (1/2) . | Author . | N
. | z_{1} (mole frac.)
. | T (K)
. | $AADp\sigma $ (%) . |
---|---|---|---|---|---|

R-1234yf/134a | Kamiaka et al.^{35} | 67 | 0.00–1.00 | 273–333 | 0.17 |

R-1234yf/134a | Chen et al.^{36} | 41 | 0.48–0.58 | 268–323 | 4.10 |

R-1234yf/1234ze(E) | Al Ghafri et al.^{19} | 3 | 0.66–0.73 | 274–342 | 3.58 |

R-1234yf/1234ze(E) | Ye et al.^{37} | 77 | 0.00–1.00 | 284–333 | 0.34 |

R-134a/1234ze(E) | Al Ghafri et al.^{19} | 3 | 0.49–0.57 | 274–341 | 1.79 |

R-134a/1234ze(E) | Kou et al.^{38} | 40 | 0.00–1.00 | 293–323 | 0.89 |

Pair (1/2) . | Author . | N
. | z_{1} (mole frac.)
. | T (K)
. | $AADp\sigma $ (%) . |
---|---|---|---|---|---|

R-1234yf/134a | Kamiaka et al.^{35} | 67 | 0.00–1.00 | 273–333 | 0.17 |

R-1234yf/134a | Chen et al.^{36} | 41 | 0.48–0.58 | 268–323 | 4.10 |

R-1234yf/1234ze(E) | Al Ghafri et al.^{19} | 3 | 0.66–0.73 | 274–342 | 3.58 |

R-1234yf/1234ze(E) | Ye et al.^{37} | 77 | 0.00–1.00 | 284–333 | 0.34 |

R-134a/1234ze(E) | Al Ghafri et al.^{19} | 3 | 0.49–0.57 | 274–341 | 1.79 |

R-134a/1234ze(E) | Kou et al.^{38} | 40 | 0.00–1.00 | 293–323 | 0.89 |

Pair (1/2) . | Author . | N
. | z_{1} (mole frac.)
. | T (K)
. | AAD_{ρ} (%)
. | N_{fail}
. |
---|---|---|---|---|---|---|

R-1234yf/134a | Yotsumoto et al.^{39} | 575 | 0.00–0.82 | 263–323 | 0.02 | 0 |

R-1234yf/134a | Akasaka et al.^{17} | 22 | 0.32–0.72 | 350–371 | 1.66 | 6 |

R-1234yf/134a | Chen et al.^{40} | 94 | 0.04–0.86 | 299–403 | 0.39 | 0 |

R-1234yf/1234ze(E) | Higashi^{18} | 14 | 0.50–0.50 | 355–374 | 2.53 | 5 |

R-1234yf/1234ze(E) | Higashi^{41} | 52 | 0.50–0.50 | 340–430 | 1.56 | 5 |

R-1234yf/1234ze(E) | Al Ghafri et al.^{19} | 37 | 0.50–0.50 | 252–404 | 0.18 | 0 |

R-134a/1234ze(E) | Zhang et al.^{42} | 101 | 0.36–0.57 | 270–300 | 0.24 | 0 |

R-134a/1234ze(E) | Al Ghafri et al.^{19} | 59 | 0.50–0.50 | 252–403 | 0.11 | 0 |

Pair (1/2) . | Author . | N
. | z_{1} (mole frac.)
. | T (K)
. | AAD_{ρ} (%)
. | N_{fail}
. |
---|---|---|---|---|---|---|

R-1234yf/134a | Yotsumoto et al.^{39} | 575 | 0.00–0.82 | 263–323 | 0.02 | 0 |

R-1234yf/134a | Akasaka et al.^{17} | 22 | 0.32–0.72 | 350–371 | 1.66 | 6 |

R-1234yf/134a | Chen et al.^{40} | 94 | 0.04–0.86 | 299–403 | 0.39 | 0 |

R-1234yf/1234ze(E) | Higashi^{18} | 14 | 0.50–0.50 | 355–374 | 2.53 | 5 |

R-1234yf/1234ze(E) | Higashi^{41} | 52 | 0.50–0.50 | 340–430 | 1.56 | 5 |

R-1234yf/1234ze(E) | Al Ghafri et al.^{19} | 37 | 0.50–0.50 | 252–404 | 0.18 | 0 |

R-134a/1234ze(E) | Zhang et al.^{42} | 101 | 0.36–0.57 | 270–300 | 0.24 | 0 |

R-134a/1234ze(E) | Al Ghafri et al.^{19} | 59 | 0.50–0.50 | 252–403 | 0.11 | 0 |

#### 5.2.1. VLE

For the mixture R-1234yf/134a, the VLE dataset from Shimoura *et al.*^{16} has been dropped because the deviations are greater than 30%. These large deviations suggest a systematic error in the measurements and also call into question the SOS data from the same reference, perhaps explaining the significant deviations from the measurements of Rowane and Perkins.^{12} Otherwise, the VLE dataset from Chen *et al.*^{36} shows large deviations, some greater than 10% compared with the model from this work, while that of Kamiaka *et al.*^{35} agrees well with the model developed in this work, with deviations mostly significantly less than 1% and an AAD of 0.2%.

#### 5.2.2. PVT

The comparisons against existing experimental data are more challenging for density because of the near-critical region. These calculations are complicated by the fact that a solution may not exist thermodynamically or that the initial guesses may not be adequate to converge to the correct solution.

As is discussed in the work of Bell *et al.*,^{3} mixture models may cause problems for the iterative calculations in the critical region. Figure 11 provides a demonstration. The isopleth of the phase envelope of the new mixture model (the solid curve) falls just below the points in the critical region (worst deviations on the order of 0.3 K in the temperature direction). Saturation temperatures above the maximum of the curve are not accessible by the mixture model. The model from this work but with the departure term disabled (with *F*_{ij} = 0), the dashed curve, is in better agreement with the data in the critical region but provides a worse representation of the other high-accuracy data.

Taking into consideration the challenges of assessing density data, deviation plots for the existing data are shown in Fig. 12.

For R-1234yf/134a, the datasets of Yotsumoto *et al.*^{39} and Chen *et al.*^{40} are in good agreement with the new model, although systematic composition-dependent errors can be seen for Chen *et al.*^{40} Although the deviations for Akasaka *et al.*^{17} are generally larger, this error metric exaggerates the errors^{3} as a consequence of the flatness of the saturation curve near the critical point (see Fig. 11).

For R-1234yf/1234ze(E), all the measurements are at the same bulk composition. The deviations of Al Ghafri *et al.*^{19} are small, while Higashi^{41} scatters more, though the same discussion about the exaggeration of the error metric applies.

For R-134a/1234ze(E), the two other datasets^{19,42} are in good agreement with the data. The datasets cover a relatively limited range of composition.

#### 5.2.3. Critical region

Although mixture critical locus data were not included in the model development, there are a few experimental data points available in the literature.^{17,41} These critical loci data are shown in Fig. 13, along with curves calculated from the new mixture model in this work. The curves were traced with the teqp^{23} library. No VLE data were included in the model development for either binary pair, so it is comforting to see that the temperature and pressure of the critical loci are well represented with the new mixture models according to the admittedly very limited data in the literature.

An interesting point about the critical locus for R-1234yf/134a is that it has a temperature minimum at 367.69 K, which is ∼0.1 K below the critical point of R-1234yf. This means that for temperatures between the temperature minimum of the critical curve and the critical temperature of R-1234yf, there are two critical points and the *p*–*x* curve is in two unconnected portions originating from the pure fluids. This behavior is not uncommon;^{43} the mixture carbon dioxide + ethane also shows this behavior.^{44}

### 5.3. Extrapolation

A new means of assessing the extrapolation behavior of EOS is that of the effective hardness of interaction,^{45} which we denote as *n*_{eff}. In the dilute-gas limit, this quantity is defined by second virial coefficients *B*_{2} and their temperature derivatives,

and is plotted for two models in Fig. 14 in the dilute-gas limit as a function of temperature and composition. While the EOS for R-134a has the wrong infinite temperature limit (which should approach 3/2 for all potentials that are finitely valued at all separations; derivation in appendix of Ref. 46), the mixture models appear to smoothly transition between the two pure fluids. This is as expected, but comforting to see nonetheless. The EOS for R-1234ze(E) and R-1234yf demonstrate a small bump in the vicinity of 200 K, which is suspect because neither the Stockmayer fluid (a common model potential for molecules with dipole–dipole interactions) nor molecular models for water^{47} (a molecule with strong gas-phase association) show the bumps.^{45} Otherwise, the shapes of the *n*_{eff} curves are qualitatively similar to those of the small rigid molecules for which *ab initio* data are available.^{45}

## 6. VLE-Only Binary Pairs

The second set of model results are for the binary pairs for which we included in the optimization only bubble-point pressure measurements from our group. In this case, an alternative model fitting procedure was employed. The method described in detail in Ref. 25 was used, which uses an evolutionary optimization approach in concert with iterative bubble-point pressure calculations from REFPROP 10.0.^{4} The interaction parameters *β*_{T,ij} and *γ*_{T,ij} were optimized, the values of *β*_{v,ij} and *γ*_{v,ij} were set to 1.0, and no departure function was used. The obtained parameters are in Table 9, and the deviations are shown in Fig. 15. The AAD for the bubble-point pressures from our measurements are all less than 0.2%. For the binary pair R-125/R-1234yf, the model is able to fully capture the bubble-point data by fitting *β*_{T,ij} and *γ*_{T,ij}, showing no systematic errors and an AAD less than 0.1% (which is much smaller than the experimental uncertainty). For the other two binary pairs, fitting only *β*_{T,ij} and *γ*_{T,ij} does not appear to be sufficient to remove the systematic temperature deviation, although the AAD are all still below 0.2%.

Pair (1/2) . | β_{T,ij}
. | γ_{T,ij}
. | β_{v,ij}
. | γ_{v,ij}
. | F_{ij}
. |
---|---|---|---|---|---|

R-125/1234yf | 0.999 637 | 0.999 356 | 1.0 | 1.0 | 0.0 |

R-1234yf/152a | 1.002 918 | 0.983 928 | 1.0 | 1.0 | 0.0 |

R-1234ze(E)/227ea | 1.000 895 | 0.993 523 | 1.0 | 1.0 | 0.0 |

Pair (1/2) . | β_{T,ij}
. | γ_{T,ij}
. | β_{v,ij}
. | γ_{v,ij}
. | F_{ij}
. |
---|---|---|---|---|---|

R-125/1234yf | 0.999 637 | 0.999 356 | 1.0 | 1.0 | 0.0 |

R-1234yf/152a | 1.002 918 | 0.983 928 | 1.0 | 1.0 | 0.0 |

R-1234ze(E)/227ea | 1.000 895 | 0.993 523 | 1.0 | 1.0 | 0.0 |

The other datasets from the literature are listed in Table 10 and their deviations are plotted in Fig. 15. For R-125/1234yf and R-1234yf/152a, the existing data do not completely agree with our measured data, demonstrating scatter in pressure up to 2%.

Pair (1/2) . | Kind . | Author . | N
. | T (K)
. | AAD (%) . |
---|---|---|---|---|---|

R-125/1234yf | VLE | Kamiaka et al.^{49} | 28 | 273–333 | 0.55 |

R-125/1234yf | VLE | Kamiaka et al.^{35} | 84 | 273–333 | 0.28 |

R-125/1234yf | VLE | Yang et al.^{50} | 35 | 283–323 | 0.16 |

R-125/1234yf | PVT | Dang et al.^{48} | 27 | 284–318 | 0.47 |

R-125/1234yf | PVT | Al Ghafri et al.^{19} | 40 | 252–383 | 0.26 |

R-1234yf/152a | VLE | Hu et al.^{51} | 60 | 283–323 | 0.70 |

R-1234yf/152a | VLE | Yang et al.^{52} | 25 | 283–323 | 0.60 |

R-1234yf/152a | PVT | Tomassetti^{14} | 136 | 268–373 | 0.16 |

Pair (1/2) . | Kind . | Author . | N
. | T (K)
. | AAD (%) . |
---|---|---|---|---|---|

R-125/1234yf | VLE | Kamiaka et al.^{49} | 28 | 273–333 | 0.55 |

R-125/1234yf | VLE | Kamiaka et al.^{35} | 84 | 273–333 | 0.28 |

R-125/1234yf | VLE | Yang et al.^{50} | 35 | 283–323 | 0.16 |

R-125/1234yf | PVT | Dang et al.^{48} | 27 | 284–318 | 0.47 |

R-125/1234yf | PVT | Al Ghafri et al.^{19} | 40 | 252–383 | 0.26 |

R-1234yf/152a | VLE | Hu et al.^{51} | 60 | 283–323 | 0.70 |

R-1234yf/152a | VLE | Yang et al.^{52} | 25 | 283–323 | 0.60 |

R-1234yf/152a | PVT | Tomassetti^{14} | 136 | 268–373 | 0.16 |

The existing density datasets for these three binary mixtures are those of Dang *et al.*,^{48} Al Ghafri *et al.*,^{19} and Tomassetti.^{14} The data from Dang *et al.*^{48} deviate up to 1% in density. Al Ghafri *et al.*^{19} claimed a combined expanded uncertainty of 0.2% in the liquid phase and 2% in the gas phase; their density data are reproduced within this band. Tomassetti^{14} included some two-phase points, resulting in nonsensical density deviation values. The two-phase points were dropped, and the phase was specified to be gas.

## 7. Validation and Verification

The models presented in this work were implemented in REFPROP and CoolProp. The necessary files are provided in the supplementary material. Calculated values are provided in Table 11, and it was confirmed that the upcoming version of REFPROP and CoolProp version 6.4.2 yielded the same results to 13 digits. The new equation of state for R-1234yf of Lemmon and Akasaka^{6} is also provided in REFPROP and CoolProp formats. An updated fluid file for REFPROP for R-152a is also provided, with the coefficients rounded to agree with the original publication^{8} in order to yield the precise values in the table.

Names (1/2) . | z_{1} (mole frac.)
. | T (K)
. | ρ (mol/m^{3})
. | T_{red} (K)
. | ρ_{red} (mol/m^{3})
. | α^{r}
. |
---|---|---|---|---|---|---|

R125 | 1.0 | 424 | 3823 | 339.173 000 000 000 0 | 4779.000 000 000 000 0 | −0.455 060 052 344 49 |

NEWR1234YF | 1.0 | 460 | 3344 | 367.850 000 000 000 0 | 4180.000 000 000 000 0 | −0.468 353 699 040 81 |

R134A | 1.0 | 468 | 3983 | 374.180 000 000 000 0 | 4978.830 171 000 001 4 | −0.466 824 484 145 93 |

R152A | 1.0 | 483 | 4457 | 386.411 000 000 000 0 | 5571.449 999 999 999 8 | −0.507 421 495 701 51 |

R1234ZEE | 1.0 | 478 | 3432 | 382.513 000 000 000 0 | 4290.000 000 000 000 0 | −0.463 409 784 472 30 |

R227EA | 1.0 | 469 | 2796 | 374.900 000 000 000 0 | 3495.000 000 000 000 0 | −0.442 385 761 979 82 |

NEWR1234YF+R1234ZEE | 0.4 | 469 | 3399 | 375.368 708 235 417 6 | 4248.595 802 001 349 5 | −0.460 594 641 762 52 |

NEWR1234YF+R134A | 0.4 | 462 | 3698 | 369.337 535 207 333 2 | 4622.435 740 472 200 5 | −0.465 508 591 288 31 |

R134A+R1234ZEE | 0.4 | 472 | 3639 | 377.666 767 145 286 7 | 4548.679 872 122 646 9 | −0.462 451 303 341 93 |

R125+NEWR1234YF | 0.4 | 445 | 3523 | 356.118 063 976 404 5 | 4403.741 839 301 608 7 | −0.463 908 051 082 07 |

NEWR1234YF+R152A | 0.4 | 470 | 3947 | 376.126 334 562 794 2 | 4933.399 297 446 838 3 | −0.492 524 178 122 31 |

R1234ZEE+R227 EA | 0.4 | 471 | 3025 | 376.790 938 517 181 3 | 3781.010 485 718 356 4 | −0.451 298 816 594 40 |

Names (1/2) . | z_{1} (mole frac.)
. | T (K)
. | ρ (mol/m^{3})
. | T_{red} (K)
. | ρ_{red} (mol/m^{3})
. | α^{r}
. |
---|---|---|---|---|---|---|

R125 | 1.0 | 424 | 3823 | 339.173 000 000 000 0 | 4779.000 000 000 000 0 | −0.455 060 052 344 49 |

NEWR1234YF | 1.0 | 460 | 3344 | 367.850 000 000 000 0 | 4180.000 000 000 000 0 | −0.468 353 699 040 81 |

R134A | 1.0 | 468 | 3983 | 374.180 000 000 000 0 | 4978.830 171 000 001 4 | −0.466 824 484 145 93 |

R152A | 1.0 | 483 | 4457 | 386.411 000 000 000 0 | 5571.449 999 999 999 8 | −0.507 421 495 701 51 |

R1234ZEE | 1.0 | 478 | 3432 | 382.513 000 000 000 0 | 4290.000 000 000 000 0 | −0.463 409 784 472 30 |

R227EA | 1.0 | 469 | 2796 | 374.900 000 000 000 0 | 3495.000 000 000 000 0 | −0.442 385 761 979 82 |

NEWR1234YF+R1234ZEE | 0.4 | 469 | 3399 | 375.368 708 235 417 6 | 4248.595 802 001 349 5 | −0.460 594 641 762 52 |

NEWR1234YF+R134A | 0.4 | 462 | 3698 | 369.337 535 207 333 2 | 4622.435 740 472 200 5 | −0.465 508 591 288 31 |

R134A+R1234ZEE | 0.4 | 472 | 3639 | 377.666 767 145 286 7 | 4548.679 872 122 646 9 | −0.462 451 303 341 93 |

R125+NEWR1234YF | 0.4 | 445 | 3523 | 356.118 063 976 404 5 | 4403.741 839 301 608 7 | −0.463 908 051 082 07 |

NEWR1234YF+R152A | 0.4 | 470 | 3947 | 376.126 334 562 794 2 | 4933.399 297 446 838 3 | −0.492 524 178 122 31 |

R1234ZEE+R227 EA | 0.4 | 471 | 3025 | 376.790 938 517 181 3 | 3781.010 485 718 356 4 | −0.451 298 816 594 40 |

^{a}

CoolProp version 6.4.2dev, revision be204ce8d5b877246c2954a9948c3d72302ff587 (SHA1 hash).

## 8. Conclusions

Most importantly, this study demonstrates that high-accuracy liquid-phase SOS and density data are largely sufficient to fit highly accurate thermodynamic models for mixtures of chemically similar HFO- and HFC-containing binary mixtures. With the fitting of four interaction parameters and small departure functions, the experimental data can be reproduced to nearly their experimental uncertainty. In addition, the availability of highly accurate reference mixture measurements allows the quality of the datasets in the literature to be assessed, identifying a number of discrepancies in the mixture data and recommending which datasets should be considered. The highly accurate mixture data show that the equation of state for R-1234ze(E) needs to be refit in order to better reproduce the new liquid-phase SOS data.

In a parallel project, further measurements are under way on some of the binary mixtures for which only VLE measurements were available in this work. These measurements will allow refinement of the interim models up to the level of the other mixture models obtained in this work.

## 9. 9. Supplementary Material

In order to ensure reproducibility of the results, the supplementary material includes (1) the Python code used to generate the check values as well as the fluid files needed for the new EOS for R-1234yf, (2) models in the formats needed for REFPROP and CoolProp/teqp, and (3) the fitting code and the experimental data used in the fitting process.

## Acknowledgments

This work was partially supported by the Strategic Environmental Research and Development Program (Project WP19-1385: WP-2740 Follow-On: Low-GWP Alternative Refrigerant Blends for HFC-134a). The author would like to thank Demian Riccardi (of NIST) for help accessing data of TDE.

## 10. AUTHOR DECLARATIONS

### 10.1. Conflict of Interest

The author has no conflicts to disclose

## 11. DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material. Additional information is available from the corresponding author upon reasonable request.