Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp.

Over the last few decades, researchers have developed a number of empirical and theoretical models for the correlation and prediction of the thermophysical properties of pure fluids and mixtures treated as pseudo-pure fluids. In this paper, a survey of all the state-of-the-art formulations of thermophysical properties is presented. The most-accurate thermodynamic properties are obtained from multiparameter Helmholtz-energy-explicit-type formulations. For the transport properties, a wider range of methods has been employed, including the extended corresponding states method. All of the thermophysical property correlations described here have been implemented into CoolProp, an open-source thermophysical property library. This library is written in C++, with wrappers available for the majority of programming languages and platforms of technical interest. As of publication, 110 pure and pseudo-pure fluids are included in the library, as well as properties of 40 incompressible fluids and humid air. The source code for the CoolProp library is included as an electronic annex.

Introduction

A number of thermophysical property libraries exist that implement the highest-accuracy formulations for the thermophysical properties of fluids. The most widely used library is REFPROP,[1] a product of the United States National Institutes of Standards and Technology (NIST). In addition, there are a number of other thermophysical property libraries, each with varying capabilities and goals. These thermophysical property libraries are summarized in Table 1.

Table 1

Software Packages Implementing High-Accuracy Equations of State for Pure and Pseudo-pure Fluids

library name	reference	fluids	open-source	mixtures	notes
REFPROP 9.1	(1)	127	no	yes	wrappers available for numerous languages
CoolProp 4.0	(23)	110	yes	no	wrappers available for numerous languages
EES	(24)	88	no	limited
FLUIDCAL	(25)	70	no	no
Zittau	(26)	34	no	no
FPROPS	(27)	36	yes	no
HelmholtzMedia	(28)	9	yes	no	only for use with Modelica

In addition, there are a few open-source thermophysical property libraries. Unfortunately, the state-of-the-art in open-source thermophysical property libraries is not very mature, apart from the CoolProp library presented here. The primary benefit of developing an open-source thermophysical library is that it facilitates easy collaboration because the source code can be read, modified, and improved by anyone in the world. Furthermore, by developing a free, open-source, thermophysical property library, researchers all over the world can get access to state-of-the-art formulations for the thermophysical properties of fluids. Access to these high-accuracy properties will improve the quality of the research carried out in a wide range of technical fields.

The major limitation of CoolProp, and most of the other libraries as well, is that they can not handle mixtures of fluids. The treatment of mixtures of fluids introduces a great amount of complexity and numerical challenges compared with the evaluation of the thermodynamic properties of pure fluids. A description of the methods required for mixtures can be found in the literature.[2−6]

The state-of-the-art in thermodynamic property modeling is quite mature. Reference-quality equations of state, which can reproduce all experimental measurements within their experimental uncertainties, have been fit for a few pure fluids of technical interest. Methodologies have been proposed, such as the fixed form equation of state developed by Span and Wagner for polar[7] and nonpolar[8] fluids to more readily fit equations of state for other fluids for which less experimental data are available. Span et al.[9] provide a review of the state of art in the high-accuracy equations of state as of the year 2001.

Since the review of Span et al.[9] was published, high-accuracy equations of state have been published in the literature for sulfur hexafluoride,[10] para-, ortho-, and normal hydrogen,[11] propane,[12] ethane,[13]n-butane, and isobutane,[14] pentafluoroethane (R125),[15] ethanol[16] and nitrogen.[17] Additional pure fluid equations of state for cyclopentane,[18] helium,[19] propylene,[20] refrigerant R227ea,[21] refrigerant R365mfc,[21] and Solkatherm 36[22] have been constructed by other researchers that have not yet been published as of publication.

From the standpoint of transport property modeling, the state-of-the-art is less mature. Partly, this is due to the fact that, in order to develop a high-accuracy transport property correlation, a high-accuracy formulation for the thermodynamic properties is required in order to evaluate the density for given temperature and pressure. For that reason, there tends to be at least a few years lag between the publication of the equation of state and the transport property correlations. In addition, there is a general shortage of experimental data of transport properties.

In recent years, a number of high-accuracy correlations for transport properties have been developed, and as of publication, 36 fluids have fluid-specific correlations for their transport properties. These fluids are summarized in a table in the Supporting Information.

Thermodynamic Properties

The thermodynamic properties of all the fluids that are implemented in CoolProp are based on Helmholtz-energy-explicit equations of state. This formulation is currently employed for all the high-accuracy equations of state that are available in the literature. Span[29] provides further information on this formulation. Furthermore, equations of state based on the Bender[30] or modified Benedict–Webb–Rubin (mBWR) forms can be converted to Helmholtz-energy-explicit forms using the methods presented in Span.[29]

In the Helmholtz-energy-explicit formulation, the total nondimensionalized Helmholtz energy α can be given as the sum of two contributions: the residual (αr) and ideal-gas (α0) parts. Thus, the nondimensionalized Helmholtz energy can be given by

The elegance of this formulation is that all other thermodynamic properties can be obtained through analytic derivatives of the terms α0 and αr. For instance, the pressure can be obtained fromwhere Z is the compressibility factor, p is the pressure in kPa, ρ is the density in kg·m–3, R is the mass specific gas constant in kJ·kg–1·K–1, T is the temperature in Kelvin, the reduced density δ is given by δ = ρ/ρred, and the reciprocal reduced temperature is given by τ = Tred/T.

The reducing density ρred is generally the critical density ρc and the reducing temperature Tred is generally the critical temperature Tc. For the pseudo-pure fluids (Air, R404A, R410A, R407C, R507A, and SES36), selected siloxanes (MM, MD4M, D4, and D5), refrigerant R134a, and methanol, the reducing parameters ρred and Tred are determined as part of the fitting process.

The other fundamental thermodynamic properties can be obtained directly using the fundamental equation of state. The enthalpy is obtained fromwhere h is the enthalpy in kJ·kg–1, and the entropy is obtained fromwhere s is the entropy in kJ·kg–1·K–1.

Additionally, other thermodynamic parameters (speed of sound, specific heats, derivatives, etc.) can be obtained analytically. Lemmon et al.[5] and Span[29] provide thorough coverage of these derivatives and thermodynamic properties. Furthermore, other combinations of partial derivatives, as well as analytic derivatives along the saturation curves and in the two-phase region can be found in the work of Thorade and Sadat.[31]

The Helmholtz-energy-explicit equations of state use temperature and density as the independent variables. If other state variables are given, it is necessary to employ numerical solvers to obtain temperature and density given the other set of inputs. Span[32] provides a description of how to handle the input state variables of temperature/pressure, pressure/density, pressure/enthalpy, and pressure/entropy. Additionally, a solver for enthalpy/entropy inputs has been implemented in CoolProp.

Helmholtz Energy Components

Residual Component

In general, the form of the residual Helmholtz energy is fluid dependent and is obtained by an optimization routine that selects terms from a large library of candidate terms. This process is described in some depth in the literature.[5,15,29,33] For the residual Helmholtz energy term, there are generally six families of terms that have been employed throughout the equations of state. The residual Helmholtz energy is given by a summationwhere each term αkr is differentiable analytically with respect to δ and τ.

The types of terms that have been used in the literature in equations of state arewhere

Analytic partial derivatives of each family with respect to τ and δ up to the second order derivatives can be found in the referenced paper for each family. Furthermore, the values for the coefficients n, t, d, etc. are presented in each equation of state. All the permutations of third order partial derivatives have also been implemented in CoolProp. These higher order analytic derivatives are required in order to implement analytic derivatives for the Tabular Taylor Series Expansion (TTSE) method[35] or bicubic interpolation as described below.

Ideal Gas Component

Like the residual Helmholtz energy, the form of the ideal-gas part of the Helmholtz energy is also fluid dependent. The ideal-gas Helmholtz energy is obtained from the relationshipand thus, the ideal-gas part of the Helmholtz energy can be obtained if the reference state parameters ρ0, T0, h00, and s00 are specified and the ideal-gas isobaric specific heat cp0(T)/R relationship is known. The reference state parameters ρ0, T0, h00, and s00 are selected in order to yield the desired values for enthalpy and entropy at the reference state. The integration in eq 15 must be carried out in order to use the ideal-gas contribution in the equation of state.

Over the years numerous forms for the ideal-gas specific heat have been implemented, including Plank–Einstein terms,[33] Aly–Lee terms,[36] and polynomial terms.

Vapor–Liquid Equilibrium

In the vapor–liquid two-phase region, as well as along the saturation curves, it is necessary to evaluate the phase equilibrium between the saturated liquid and the saturated vapor.

For a pure fluid at equilibrium, the temperatures, pressures and Gibbs free energy in each phase are the same. Thus, for a given saturation temperature Ts, the system of equations to be solved iswhere the unknowns are the saturated liquid density ρ′ and the saturated vapor density ρ″.

The method proposed by Akasaka[37] is employed, which is a two-dimensional Newton–Raphson solver for the nonlinear system of equations from eqs 16 and 17. For a given temperature Ts, this method yields the solutions for the saturation pressure ps, the saturated liquid density ρ′ and the saturated vapor density ρ″. Figure 1 shows the saturation curves for all the fluids included in CoolProp in reduced coordinates.

Figure 1

Saturation curves for all fluids included in CoolProp.

This solver begins with initial guess values for ρ′(T) and ρ″(T) provided by the ancillary equations. For fluids without published ancillary curves, ancillary curves for ρ′(T), ρ″(T), and p(T) have been fit using routines provided in the CoolProp package. In general, the combination of highly accurate ancillary equations and the Newton–Raphson method yields proper convergence for temperatures where Tt < T < (Tc–0.01 K) where Tt is the triple point temperature. When the Newton–Raphson method fails with the normal method, a relaxation parameter can be introduced to yield better convergence behavior in the near-critical region.

In the near vicinity of the critical point, the behavior of the saturation solvers becomes significantly less robust, even with good guess values for the saturation densities from the ancillary equations. As a result, it is necessary to employ other methods to extend the saturation curves all the way up to the critical temperature. The solvers presented above are used to get as close to the critical temperature as possible. Beyond that point, a spline curve is used for the saturation curve, where the value and derivative constraints can be obtained from the last point that the Newton–Raphson method succeeded at temperature Tend. The constraints on the spline for the saturated liquid density arewhere the right-hand side of each constraint is evaluated analytically from the equation of state. A similar spline is constructed for the saturated vapor density as a function of the temperature. This yields a smooth (C1 continuous) transition from the EOS to the critical region spline. Furthermore, the critical spline is imposed to yield the correct value for the density at the critical point. For each fluid, the value of Tend and the saturation derivatives at Tend are precalculated and cached in order to maximize computational efficiency.

Figure 2 shows the range of the saturation curve that is treated using a spline curve as a function of the ratio of the critical temperature to the triple point temperature. For fluids with well constructed equations of state and good ancillary equations, the numerical VLE solver succeeds at temperatures within 1 × 10–9 K of the critical temperature, but for refrigerants R11 and R14, the saturation solvers fail at a distance greater than 0.1 K from the critical point.

Figure 2

Range of critical spline versus the ratio of critical to triple point temperatures for fluids with Tc – Tend > 1 × 10–7 K.

It is a common need to obtain the saturation temperature for a given saturation pressure. The saturation pressure curves as a function of temperature are continuous from the triple point temperature to the critical point temperature. Some fluids have equations of state where the minimum temperature is above the triple point temperature. Therefore, it is straightforward to obtain the saturation temperature for the given saturation pressure.

There are several means of implementing this solution procedure. The most robust is the use of Brent’s method,[38] which is a bounded one-dimensional solver with quadratic updates and guaranteed convergence. Brent’s method[38] is used to drive the residuumto zero. The saturation temperature Ts is the independent variable, which is known to lie within the closed range between the triple point temperature and the critical point temperature. The solution is found when the saturation pressure ps(Ts) (evaluated from the vapor–liquid equilibrium solver routine) is equal to the target pressure ptarget.

In the case of pseudo-pure fluids (Air, refrigerant R404A, refrigerant R410A, etc.), it is not possible to determine the vapor–liquid equilibrium with the use of the phase equilibria from eqs 16 and 17. For these mixtures, at equilibrium, the mole fractions of each component are not the same in the vapor and liquid phases and the pseudo-pure fluid equation of state can only calculate properties for the pseudo-pure fluid composition. The saturated liquid and vapor ancillary pressure equations are thus no longer optional but required to calculate the saturation pressures. The pressures calculated from the ancillary equations are then used to evaluate the saturation densities using the equation of state.

Interpolation Methods

When using equations of state in engineering applications, computational efficiency is of the utmost importance. In order to improve the speed of evaluation of the equation of state, interpolation methods can be used. While a comparison of interpolation methods is beyond the scope of this work, two interpolation methods that have been found to yield excellent behavior are the Tabular Taylor Series Expansion (TTSE) method and the bicubic interpolation method. These two methods share the requirement that values of state variables are tabulated on a regularly (either linearly or logarithmically) spaced grid, as well as derivatives of the state variable with respect to the two independent variables.

Using the TTSE method, with pressure and enthalpy as independent variables, the temperature can be obtained from the expansionwhere the derivatives are evaluated at the point i,j, and the differences are given by Δp = p – p and Δh = h – h. The nearest state point can be found directly due to the regular spacing of the grid of points. Pressure and enthalpy are very common state variable inputs in the simulation of thermal engineering systems.

For an improved representation of the p–v–T surface, bicubic interpolation can be used. In the bicubic interpolation method, the state variable and its derivatives are known at each grid point. This information is used to generate a bicubic representation for the property in the cell, which could be expressed aswhere a are constants based on the cell boundary values and x and y are normalized values for the enthalpy and pressure, for instance. The constants a in each cell are cached for additional computational speed.

As an example of the increase in computational speed possible through the use of these interpolation methods, the density is calculated as a function of the pressure and enthalpy for subcooled water. The IAPWS 1995 formulation for the equation of state of ordinary water[33] is one of the most involved equations of state in the literature. For subcooled water at a pressure of 10 MPa and an enthalpy of 475 kJ·kg–1 (where the reference enthalpy is 0.611872 J·kg–1 for the saturated liquid at the triple point), both the TTSE method and the bicubic interpolation method are more than 120 times faster than the evaluation of the density from the equation of state (it takes approximately 1 μs to evaluate density using the TTSE or bicubic interpolation methods). Thus, using one of these interpolation methods can yield a reduction in computational time of greater than 98%.

Practical implementation of these methods involves building tables with the fluid properties and their derivatives at each grid point. This task is only performed at the first property call and takes only a few seconds. The tables are then cached for further use.

As another example of the accuracy of these interpolation methods, the density of refrigerant R245fa is evaluated at 40000 data points covering the entire fluid surface. Figure 3 shows the results of this analysis. These data show that the accuracy of the bicubic interpolation method is generally several orders of magnitude better than that of the TTSE method, though both yield acceptable accuracy for most technical needs.

Figure 3

Comparison of the accuracy of TTSE and bicubic interpolation methods for refrigerant R245fa (interpolation grid is 200 × 200, enthalpy spaced linearly, pressure spaced logarithmically).

Transport Properties

For the transport properties (here viscosity, thermal conductivity, and surface tension), the state-of-the-art is less mature. A wider range of methodologies have been employed to correlate and/or predict these properties. For a number of fluids, high-accuracy fluid-specific correlations have been developed based on wide-ranging experimental data, but for others, little or no experimental data are available and predictive or empirical methods must be employed.

Pure Fluid Correlations

Viscosity

Correlation of the viscosity of pure fluids is typically divided into two contributions: one part provides the temperature-dependent viscosity in the zero-density limit (dilute-gas), and the second part considers the temperature- and density-dependent residual viscosity, as in

For a very restricted subset of fluids, there is sufficient information about viscosity in the critical region to consider the critical enhancement of the viscosity. In general, the critical enhancement of viscosity is not considered. Of all the pure fluid viscosity correlations developed, the only ones with a critical enhancement term for the viscosity are ordinary water[39] and carbon dioxide.[40]

It is possible to theoretically treat the zero-density viscosity using Chapman–Enskog theory, which yields the dilute gas viscosity ofwhere η(0) is the viscosity in the limit of zero density in μPa·s, M is the molar mass in kg·kmol–1, T is the temperature in Kelvin, ση is the size parameter of the Lennard-Jones model in nm, and Ω(2,2) is the empirical collision integral given by the form from Neufeld[41]where T* is the reduced temperature defined by T* = kT/εη and where the ratio εη/k of the pair potential energy to Boltzmann’s constant is in Kelvin and is fluid dependent. For fluids that are well characterized by experimental data, it is possible to fit the term Ω(2,2) to experimental data. Also, for fluids for which the terms ση and εη/k are unknown, they can be estimated based on the method from Chung et al. in eqs 51 and 52.

The residual viscosity η( can be treated in a variety of different ways. In older viscosity correlations, it was common practice to develop an empirical correlation for η( directly. In the last 15 years, the preponderance of pure fluid viscosity correlations[42−44] have been based on the division of the residual viscosity into a theoretically derived initial-density term from Rainwater–Friend theory[45,46] and a higher-order correction term. Thus, the residual viscosity is given bywhere Bη is the second viscosity virial coefficient in L·mol–1, ρ̅ is the molar density of the fluid in mol·L–1, and Δηh is the higher order correction term in μPa·s.

The second viscosity virial coefficient is given bywhere ση is the molecular size in nm, and T* = T/(εη/k) and with

The coefficients b are from Vogel et al.[42] and are duplicated in Table 2 for completeness.

Table 2

Coefficients for the Second Viscosity Virial Coefficient in Equation 30

b0	–19.572881
b1	219.73999
b2	–1015.3226
b3	2471.01251
b4	–3375.1717
b5	2491.6597
b6	–787.26086
b7	14.085455
b8	–0.34664158

The higher-order term is often of a form similar to the free-volume term proposed by Batschinski[47] and Hildebrand.[48] A general form of the higher-order term is given bywhere Δηh is in μPa·s, Tr = T/Tred, and the coefficients e and f1 are fit for each fluid. Furthermore, the term δ0(Tr) is given by the formwhere the coefficients g are also fit for the given fluid.

It should also be mentioned that the generalized friction theory model has been successfully applied to the prediction of the viscosity of some fluids, notably the n-alkanols,[49] hydrogen sulfide,[50] and sulfur hexafluoride.[51] Currently, the generalized friction theory method remains less used in the reference literature than the viscosity correlation method outlined here.

Thermal Conductivity

The correlations for thermal conductivity are decomposed into three terms, yielding the following form:where each term is in mW·m–1·K–1.

Unlike viscosity (where the critical enhancement term is very small except in the immediate vicinity of the critical point), the critical enhancement term for thermal conductivity λ( is non-negligible well away from the critical point.

The dilute gas term in the limit of zero-density λ(0) is typically correlated with the temperature using a body of low-density thermal conductivity measurements, which usually results in a short polynomial form similar to

The residual term is often given by a form similar to

Finally the critical enhancement term needs to be considered. The most commonly used critical enhancement term used is the simplified critical enhancement term of Olchowy and Sengers:[52]where λ( is in mW·m–1·K–1, ζ is in m, c and c are in kJ·kg–1·K–1, p and pc are in kPa, ρ and ρc are in kg·m–3, η is the viscosity in μPa·s, and the remaining parameters are defined in Table 3. The factor 1012 is a unit conversion parameter that yields a thermal conductivity in mW·m–1·K–1.

Table 3

Coefficients for Use in the Simplified Olchowy-Sengers Critical Term in Equations 36 to 39

Universal Constants
Boltzmann constant	k	1.3806488 × 10–23 J·K–1
universal amplitude	RD	1.03
critical exponent	ν	0.63
critical exponent	γ	1.239
reference temp.	TR	1.5Tc

Surface Tension

Mulero et al.[54] have recently refit correlations for the surface tension of nearly all the fluids in REFPROP 9.0.[55] These correlations are each of the formwhere σ is the surface tension in mN·m–1, Tc is the critical temperature in Kelvin, and a and n are correlation constants. This formulation ensures that the surface tension goes to zero at the critical point. The mean absolute percentage difference of each of these correlations is less than 6%, and most are below 3%.

For fluids that are not included in the database of Mulero, the following general form from Miqueu et al.[56] is employedwhere σ is the surface tension in N·m–1, k is the Boltzmann constant (k = 1.3806488 × 10–23 J·K–1), Tc is the critical temperature in Kelvin, NA is Avogadro’s number (NA = 6.02214129 × 1023 mol–1), Vc is the critical molar volume in m3·mol–1, ω is the accentric factor, and t = 1–T/Tc. This equation predicts the surface tension of the fluids used to develop the correlation within an average absolute difference (AAD) of 3.5%.

Extended Corresponding States

For many fluids, high-accuracy viscosity and thermal conductivity correlations are not available because these fluids have not been experimentally studied in great enough depth.

For these less-studied fluids, it is still necessary to be able to provide reasonable predictions of the viscosity and thermal conductivity over the whole fluid surface, and one method that has been used successfully is the method of extended corresponding states. In this method, the transport properties for the fluid of interest are obtained from the transport properties for a well-characterized reference fluid. The reference fluid selected should have high-accuracy transport property measurements as well as have a p–v–T surface that is similar in shape to the fluid of interest.

The analysis in this section follows the method proposed by Huber et al.,[53] which has been implemented in REFPROP.[1] The primary contribution of this section on the extended corresponding states is the presentation of a small set of example data that can be used to validate the implementation of the extended corresponding states method. No validation data for extended corresponding states has been published before. These example data are provided in Table 4 to allow for proper validation of the implementation of the extended corresponding states method.

Table 4

Data to Check ECS Implementationa

fluid of interest (EOS[57])	R124
reference fluid (EOS,[12] λ,[58] η[42])	propane
state	saturated liquid
T◇ [K]	350.000
ρ◇ [kg·m–3]	1143.37994
ρ̅◇ [mol·L–1]	8.378
Conformal State
T⊥ [K]	321.054
ρ⊥ [kg·m–3]	453.03224
ρ̅⊥ [mol·L–1]	10.274
Viscosity
ψη [−]	1.0454
η◇(0) [μPa·s]	13.617
Fη [−]	1.60328
η⊥(r) (T⊥,ρ̅⊥ψη) [μPa·s]	77.61535
η [μPa·s]	138.056
η [μPa·s] (REFPROP 9.1)	138.056
Conductivity
ψλ [−]	1.0583
fint [−]	0.0014
λ◇int [mW·m–1·K–1]	12.411
λ◇* [mW·m–1·K–1]	3.111
Fλ [−]	0.51802
λ⊥(r) (T⊥,ρ̅⊥ψλ) [mW·m–1·K–1]	70.24348
λ◇c [mW·m–1·K–1]	0.884
λ [mW·m–1·K–1]	52.794
λ [mW m–1·K–1] (REFPROP 9.1)	52.794
Correlations[59]
ψλ = 1.0898 – 1.54229 × 10–2δ
fint = 1.17690 × 10–3 + 6.78397 × 10–7T
ψη = 1.04253 + 1.38528 × 10–3δ

Note: Both CoolProp and REFPROP implement the EOS for propane from Lemmon et al.,[12] which causes errors in viscosity prediction of propane of up to 2%.

In the analysis that follows in this section, the subscript ⊥ refers to the reference fluid, and ◊ refers to the fluid of interest. Molar specific quantities are given with an overbar, and mass-specific quantities do not have an overbar.

Conformal State

The conformal state is the thermodynamic state point for the reference fluid that is used to evaluate the reference-fluid contribution to the extended corresponding states method. This conformal state point is determined based on the equivalent substance reducing ratios f and h ofor alternatively expressed in terms of shape factors θ and ϕ

The corresponding states method is most accurately applied to monatomic gases, and the shape factor can be thought of as a term that accounts for deviation from spherical molecular geometry. The shape factor can be approximated based on one of several empirical forms that have been proposed, such as those from Erickson and Ely[60] (for the reference fluid propane) ofwith a1 = 0.5202976, a2 = −0.7498189, a3 = 0.1435971, and a4 = −0.2821562 or by the more general solution (of a similar form) from Estella–Uribe and Trusler,[61] which provides higher fidelity predictions in the critical region.

For the highest accuracy and generality (shape factors independent of the reference fluid selected), it is preferable to use the “exact” shape factors, which are obtained through the use of the equations of state of the fluid of interest and the reference fluid.

The “exact” shape factors are defined based on the conformal state T⊥, ρ⊥ of the reference fluid. The conformal state is defined by equating the compressibility factor and the residual component of the Helmholtz energy of the reference fluid and the fluid of interest,[53,62,63]and

The right-hand side of each equation is known for the fluid of interest. Thus, it simply remains to obtain the conformal state point T⊥, ρ⊥ from a simultaneous solution of the two equations. The most straightforward solver to be used is a conventional two-dimensional Newton nonlinear system of equations solver.

The Newton method for the conformal state solver can be given by x = x + v where x is the vector ⟨τ⊥,, δ⊥,⟩ and where v is obtained by solving the system of equations Jv = −r and the Jacobian matrix for the solver can be given analytically bywhere each of the partial derivatives are evaluated at the state point T⊥, ρ⊥. The residual vector r is given by

This solver generally yields good convergence behavior when started at the initial guess value defined by θ = 1 and ϕ = 1.

At very low densities, the conformal solver may fail, which can be avoided by only evaluating the conformal state for densities of the fluid of interest above 1.0 kg·m–3. Below this density, the conformal state is determined by assuming θ = 1 and ϕ = 1. This treatment introduces a small discontinuity in the conformal state at the density of 1.0 kg·m–3, but in the dilute gas domain, both the viscosity and thermal conductivity are dominated by the dilute gas contribution.

In the extended corresponding states method, the viscosity is divided into two contributions, one for the dilute gas contribution in the limit of zero density, and another for the residual contribution. This division is analogous to the separation employed for fluid-specific viscosity correlations(see eq 25). Thus the viscosity can be given bywhere η◊ is the viscosity of the fluid of interest in μPa·s, η◊(0) is the dilute-gas viscosity contribution of the fluid of interest in μPa·s, and ηECS( is the residual contribution from the extended corresponding states method, in μPa·s.

The dilute gas contribution η◊(0) can be treated theoretically and is obtained from the eqs 26 and 27. If the Lennard-Jones parameters εη/k and ση are unknown for the fluid of interest, they can be obtained from the method from Chung et al.[64]where ση is in nm, ρ̅c is in mol·L–1 and Tc and εη/k are in Kelvin. If εη/k and ση are known for a reference fluid but not the fluid of interest, these values for the fluid of interest can be obtained from the method proposed in Huber et al.[53] ofwhich is simply the application of the method of Chung et al.[64] to both the reference fluid and the fluid of interest. It should be emphasized that molar densities must be used in eq 54.

The Lennard-Jones parameters εη/k and ση for a number of fluids can be found in the works of Chichester and Huber[65] and Poling et al.[66]

The residual contribution to the viscosity is obtained using the residual viscosity of the reference fluid. To begin with, the conformal temperature T⊥ and conformal molar density ρ̅⊥ are obtained using the methods presented in the conformal state section. For some fluids, there are sufficient experimental data in order to fit a simple polynomial correction in the reduced density of the fluid of interest of the form

This correction term shifts the density of the reference fluid used in the viscosity correlation away from the conformal density. If no experimental information is available to obtain ψη, ψη is assumed to be equal to 1.0. Huber et al.,[53] McLinden et al.,[62] and Klein et al.[67] provide some of the only published values for these correction polynomials. Significant work has been carried out by the authors of REFPROP to develop correction polynomials, but the density correction polynomials in REFPROP are not in the public domain.

Thus, the extended corresponding states contribution to the viscosity is obtained fromwhere η⊥((T⊥, ρ̅⊥ψη) is the contribution of the residual viscosity from the reference fluid evaluated at the temperature T⊥ and the molar density ρ̅⊥ψη. The residual viscosity of the reference fluid includes all the density-dependent terms of the viscosity correlation, which, based on the formulation in the prior section, would be the contribution from eq 28. Fη is a factor that arrives from the fact that the corresponding states theory states that the viscosity of two fluids at the same reduced state are equivalent.[67]Fη can be given bywhere h and f are the equivalent substance reducing ratios obtained from the conformal state solver, and M◊ and M⊥ are the molar masses of the fluid of interest and the reference fluid, respectively, each in kg·kmol–1.

A similar protocol is used to calculate the thermal conductivity using extended corresponding states. Again, the division of terms for thermal conductivity is similar to that of fluid-specific correlations(see eq 33). The thermal conductivity is divided into four terms, as inwhere λ◊int is the internal thermal conductivity contribution of the fluid of interest due to internal motion of the molecules, λ◊* is the dilute gas contribution from the fluid of interest, λECS( is the contribution from extended corresponding states, and λ◊( is the critical enhancement term for the fluid of interest. Each term is in mW·m–1·K–1.

The internal thermal conductivity is given bywhere λ◊int is in mW·m–1·K–1, η◊(0) is the dilute-gas viscosity in μPa·s evaluated from eqs 26 and 27, cp,◊(0) is the ideal-gas specific heat in kJ·kg–1·K–1, and R◊ is the mass-specific gas constant in kJ·kg–1·K–1. The factor fint is taken to be equal to 1.32 × 10–3, or if sufficient experimental data are available, fint is fit as a linear function of the temperature as in Huber et al.[53]

The dilute gas component is evaluated fromwhere λ◊ is in mW·m–1·K–1, R◊ is the mass-specific gas constant in kJ·kg–1·K–1, and η◊(0) is the dilute-gas viscosity in μPa·s evaluated from eqs 26 and 27.

Thus, the residual component of the thermal conductivity is obtained fromwhere λ⊥((T⊥, ρ⊥ψλ) is the contribution of the residual thermal conductivity from the reference fluid evaluated at the temperature T⊥ and the density ρ⊥ψλ. As with viscosity, for some fluids there is sufficient experimental data to fit a simple polynomial correction in the reduced density. If no experimental information is available to obtain ψλ, ψλ is assumed to be equal to 1.0.

Fλ can be given bywhere h and f are the equivalent substance reducing ratios obtained from the conformal state solver, and M◊ and M⊥ are the molar masses of the fluid of interest and the reference fluid, respectively, each in kg·kmol–1. It should be noted that Fλ differs from Fη from eq 57 in that the molar masses of each fluid are inverted.

Finally, the last component in the thermal conductivity is the critical enhancement λ◊( evaluated for the fluid of interest given by eqs 36 to 39.

CoolProp

The CoolProp library currently provides thermophysical data for 110 pure and pseudo-pure working fluids. The literature sources for the thermodynamic and transport properties of each fluid are summarized in a table in the Supporting Information available online.

The code of CoolProp is written in C++ to utilize modern C++ language features and the functionalities inherent in object oriented programming. In addition, as the code of CoolProp has been written in C++, Simplified Wrapper and Interface Generator (SWIG) can be used to readily generate an interface to any programming language that SWIG supports. As a result, fully featured high-level interfaces have been developed for most programming languages of technical interest, including Microsoft Excel, Labview, MATLAB, Python, C#, Engineering Equation Solver and many others. In addition the C++ code is cross-platform and has been successfully compiled and tested on Windows, Linux, and Mac OSX.

In addition to the inclusion of the most accurate equations of state of pure and pseudo-pure fluids, CoolProp provides the properties of eight secondary working fluids and thirteen aqueous solutions from Melinder[68] and a selection of fourteen other secondary working fluids and five brines, as well as the most accurate thermodynamic properties of humid air from Herrmann et al.[69]

The interface to the library is very straightforward. For example, from most programming languages, the code to obtain the density (’D’) of nitrogen at standard temperature (’T’) and pressure (’P’) (298.15 K and 101.325 kPa) is given by a variation on

Conclusions

In this paper, the state-of-the-art of the thermophysical properties of pseudo-pure and pure fluids has been summarized. The state-of-the-art in thermodynamic property evaluation is quite mature, with more than 100 fluids with Helmholtz-energy-explicit formulations for their equations of state. The transport properties of these fluids have been less studied, and for that reason, fluid-specific correlations for their viscosity and thermal conductivity are only available for 36 fluids. The extended corresponding states method can be used for fluids that do not have fluid-specific correlations for the transport properties.

Furthermore, all the methodologies presented above have been implemented into an open-source thermophysical library CoolProp. The current version of CoolProp as of publication is included as an electronic annex. This library is free to use and is finding increasingly wide application in a range of technical fields.

The primary limitation of this library is that it does not include mixture thermophysical properties. Mixtures of fluids are of great technical interest, and further work is ongoing to add mixture properties to this library.