Abstract
Two different methods for deriving the density and isobaric heat capacity of liquids in the subcritical pressure range, from the speed of sound, are recommended. In each method, corresponding set of differential equations relating these properties is solved as the initial boundary value problem (IBVP). The initial values are specified at the lowest pressure of the range and the boundary values along the saturation line. In the first method, numerical integration is performed along the paths connecting the Chebyshev points of the second kind between the minimum and maximum temperature at each pressure. In the second method, numerical integration is performed along the isotherms distributed in the same way, with the temperature range being extended to the saturation line after each integration step. The methods are tested with the following substances: Ar, N2, CO2, and CH4. The results obtained for the density and isobaric heat capacity have the average absolute deviation from the reference data of 0.0005% and 0.0219%, respectively. These results served as the initial values for deriving the same properties in the transcritical pressure range up to the pressure approximately twice as large as the critical pressure. The results obtained in this pressure range have respective deviations of 0.0019% and 0.1303%.
1. Introduction
The relations between the thermodynamic speed of sound (i.e., the mechanical disturbance of a small amplitude and a low frequency) and other thermodynamic properties (e.g., the density and heat capacity) comprise the set of nonlinear partial differential equations of the second-order (Trusler [1]). A general solution to this set of equations is available only for dilute imperfect gases (i.e., gases at low pressures) where the pressure effect on density may be satisfactorily described by the virial expansion truncated after the second virial coefficient, which may be obtained from the speed of sound through a model of the intermolecular potential energy. While the virial expansion does not work when applied to liquids (Allen and Tildesley [2]), it is still possible to obtain a particular solution if appropriate initial/boundary values are available from corresponding direct measurements (e.g., volumetric and caloric). In an open literature one can find several different approaches for deriving the thermodynamic properties of liquids (e.g., the density, heat capacity, isothermal compressibility, isentropic compressibility, and isobaric thermal expansivity) from the speed of sound. The majority of them are based on numerical integration of corresponding differential equations connecting these properties (Muringer et al. [3] and Sun et al. [4–6]). Other approaches include an iterative method of calculation (Petitet et al. [7]), a grid algorithm (Khasanshin et al. [8]), and a heuristic approach (Scalabrin et al. [9]).
Deriving the density and heat capacity of liquids from the speed of sound is usually performed by solving the set of differential equations, which relates these properties, as the initial value problem (IVP) for the set of ordinary differential equations. The initial values are usually specified at the lowest pressure of the range considered (Benedetto et al. [10]). This approach is applicable only to rectangular domain, that is, the one where the temperature range is constant across the entire pressure range (e.g., the supercritical pressure range). Since the IVP is highly sensitive to the initial values and since the experimental uncertainty of direct measurement of these properties is decreasing with increasing the pressure, it is preferable to specify them at pressure as close as possible to the atmospheric one. However, by decreasing the lowest pressure of the range the temperature range between a minimum temperature and that at the saturation also decreases. This may be partly overcome by imposing the initial values at the highest pressure of the range and performing the integration down to the lowest pressure along lines having similar shape to that of the saturation line. While this approach may give satisfactory results (Bijedić and Neimarlija [11]), it may not be suitable for substances with relatively high critical pressures.
In this paper an attempt is made to reconcile these two opposites, that is, to specify the initial values at the lowest pressure and, at the same time, to cover the maximum temperature range possible (i.e., to the saturation line). This lower pressure limit is chosen so as to enable sufficiently wide temperature range to accommodate a reasonable number of integration paths. Unfortunately, as it turned out, there is a compromise which has to be done. Namely, in order for the solution to be stable across the whole pressure range, the boundary values must be specified along the saturation line. Still, the majority of these values are specified at pressures which are considerably below the upper limit of the pressure range. Two different methods based on the initial boundary value problem (IBVP) are recommended. The main difference between them is in the paths of integration. In one method, these paths change their shape progressively from that of an isotherm to that of the saturation line, while in the other one they follow isotherms which are modified in each integration step to suit consecutively broader temperature ranges. To ensure that the results obtained are reliable enough they are not only compared to respective reference data but also used as the initial values for deriving the same properties in the transcritical domain up to the pressure twice of that in the critical point.
2. Theory
When the density and heat capacity are derived from the speed of sound in rectangular domain, the following set of equations may be used (Benedetto et al. [10]):where is the pressure, is the temperature, is the speed of sound, is the density, is the specific heat capacity at constant pressure, and is the thermal expansion coefficient:However, if the domain is not of rectangular shape (e.g., in the subcritical pressure range), the following set of equations is suitable (Bijedić and Neimarlija [11]):where represents a path connecting the points with temperaturewhere is the lowest temperature of the range and is the highest temperature at observed pressure (e.g., the saturation temperature). Using instead of as a path of integration enables one to shape the domain arbitrarily.
The set of (4) and (5) may be solved simultaneously for and if their initial values are specified at the lowest pressure of the range. However, for the solution to be stable across the whole pressure range, the boundary values of and need to be specified along connecting the points with the highest temperatures at each pressure (i.e., along the saturation line). If all the temperature derivatives are estimated, say, from an interpolation polynomial, the set of equations may be solved as the initial boundary value problem (IBVP) for the set of ordinary differential equations. In order to ensure that the temperature derivatives are estimated as accurate as possible, the lines of constant should be distributed so as to avoid Runge’s phenomenon (e.g., at the Chebyshev knots).
While the set of (1) and (2) is suitable for the supercritical pressure range, it may be used for the subcritical one as well, if the temperature range is being extended to the saturation line after each integration step. This set of equations may also be solved as the IBVP with the same set of initial/boundary values. However, the procedure of solution is more robust because the results from the current pressure have to be interpolated to a new set of temperatures in a wider temperature range at the next pressure.
3. Results and Discussion
The methods described are tested with several different substances. Their list, along with the ranges covered, are given in Tables 1 and 2 for the subcritical and transcritical domain, respectively. These domains are also presented graphically in coordinates at Figure 1. While this graph is given for CO2, respective graphs for other substances are qualitatively the same.
The lowest pressure of the subcritical domain is chosen so as to enable a reasonable temperature range, and the highest one corresponds to the saturation pressure for the temperature approximately 10 K below the critical point. The lowest temperature is constant across both the domains, while the highest one in the subcritical domain corresponds to the saturation temperature at each pressure.
The lowest pressure of the transcritical domain coincides with the highest pressure of the subcritical domain, and the highest one is approximately twice as large as the critical pressure. The highest temperature is constant across the domain and corresponds to the saturation temperature at the lowest pressure of the domain or to the highest temperature of the subcritical domain.
Four sets of calculations are performed for each substance. Two of them cover the subcritical domain and the other two the transcritical domain. The calculations are performed first in the subcritical domain, and the results obtained along the highest pressure of the domain served as the initial values for the transcritical domain.
For both methods to be implemented in both domains, the reference values of the density and heat capacity, as well as the speed of sound, are specified along several mainly equidistant isobars, at temperatures distributed along each isobar according to the Chebyshev points of the second kind (Berrut and Trefethen [12]):Increasing the number of these temperatures, the number of differential equations to be solved increases as well (two equations per temperature), but the temperature derivatives are estimated more accurately, and vice versa. It has been proved (by trial and error on a number of different substances and conditions) that the optimum number of these temperatures is 15. In all the cases numerical integration is performed by the Runge-Kutta-Verner (RKV) fifth-order and six-order method with adaptive step-size (Hull et al. [13]).
In the first method (IBVP-1) the initial values of and , as well as , are specified along the lowest pressure (see Table 3) at temperatures calculated from (7). The temperature derivatives of , , and and the pressure derivatives of are estimated from the Lagrangian interpolation polynomial (see Table 4). After that, there is enough information to calculate the pressure derivatives of and from (4) and (5), respectively. Then, new values of and are calculated by the RKV method at the pressure , where is the integration step chosen so as to control the norm of the local error such that the global error is proportional to the tolerance for error control (see Table 6). The procedure is repeated until the highest pressure of the domain is reached. The values of are obtained in each integration step by interpolation along s with respect to the pressure (see Table 4).
In the second method (IBVP-2) the initial values of and , as well as , are also specified along the lowest pressure (see Table 3) at temperatures calculated from (7). The temperature derivatives of and are estimated from the Lagrangian interpolation polynomial (see Table 4). After that, there is enough information to calculate the pressure derivatives of and from (1) and (2), respectively. Then, new values of and are calculated by the RKV method at the pressure , where is the integration step chosen so as to satisfy both the integration and interpolation/derivation (see Table 6). The temperature range is extended to the saturation line and new temperatures are calculated from (7). Next, and are interpolated to these new temperatures (see Table 4), and the procedure is repeated until the highest pressure of the domain is reached. The values of are obtained in each integration step also by interpolation along s with respect to the pressure (see Table 4).
In the transcritical domain, the set of (1) and (2) is solved as the initial value problem (IVP). The initial values of and are obtained from the IBVP methods at the highest pressure of the subcritical domain. The integration is performed along the same set of temperatures at which the IBVP procedures terminated (see Tables 5 and 7).
The results obtained are assessed by the average absolute deviation (AAD) from corresponding reference data (Tegeler et al. [14], Span et al. [15], Span and Wagner [16], and Setzmann and Wagner [17]). In the subcritical domain, both methods predict the density with almost the same AADs, while the results for the heat capacity obtained by the IBVP-2 method have AAD three times lower than that of the IBVP-1 method (see Table 8). However, comparing the AADs with the uncertainty of corresponding reference data (see Table 10) it is obvious that these formers are two orders of magnitude lower. Therefore, the uncertainty of the results obtained in the subcritical domain is practically the same as that of the reference data. The relative deviations of the results at the highest pressure of the subcritical domain are presented graphically as a function of temperature at Figures 2–9.
In order to test accuracy of the results additionally, their values along the highest pressure (the most unfavorable conditions) served as the initial values for deriving the same properties in the transcritical domain. Here, the results for the heat capacity obtained with the initial values from both methods have almost the same AADs, while the results for the density obtained with the initial values from the IBVP-1 method have AAD two times lower than that of the IBVP-2 method (see Table 9). However, comparing the AADs with the uncertainty of corresponding reference data (see Table 10) it is obvious that these formers are an order of magnitude lower. Therefore, the uncertainty of the results obtained in the transcritical domain is practically the same as that of the reference data (Tegeler et al. [14], Span et al. [15], Span and Wagner [16], and Setzmann and Wagner [17]).
Influence of uncertainties of the initial and boundary values and speed of sound values on the results is investigated as well. Corresponding AADs of the results obtained after changing these values in the limits of their uncertainties (see Table 10) are given in Tables 11–13. One can see that these AADs are in the limits of corresponding uncertainties from Table 10.
4. Conclusions
The density and heat capacity of a liquid may be derived from the speed of sound in the subcritical pressure range for the temperatures up to the saturation line. This may be accomplished by the use of a standard numerical procedure for solving an initial boundary value problem (IBVP) for the set of ordinary differential equations and an interpolation polynomial. The initial values are specified along the lowest pressure of the range and the boundary values along the saturation line. For the results to have uncertainty not higher than that of corresponding direct measurements, numerical derivatives must be estimated as accurately as possible prior to each integration step. This precondition may be difficult to fulfill for two reasons. The first one emerges from the fact that the set of equations is of the second-order with respect to the temperature, and the second one is that the density and heat capacity (whose temperature derivatives are estimated) change with temperature abruptly in the vicinity of the saturation line. For this to be overcome the number of these temperatures should be optimal and they have to be distributed so as to avoid Runge’s phenomenon (e.g., according to the Chebyshev knots).
Appendix
See Figures 1, 2, 3, 4, 5, 6, 7, 8, and 9 and Tables 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.