Journal of Thermodynamics
Volume 2016 (2016), Article ID 2035704, 8 pages
http://dx.doi.org/10.1155/2016/2035704

Research Article

Density and Heat Capacity of Liquids from Speed of Sound

Muhamed Bijedić and Sabina Begić

Faculty of Technology, University of Tuzla, 75000 Tuzla, Bosnia and Herzegovina

Received 5 October 2015; Revised 25 December 2015; Accepted 30 December 2015

Academic Editor: Felix  Sharipov

Copyright © 2016 Muhamed Bijedić and Sabina Begić. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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, N₂, CO₂, and CH₄. 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 temperature

where 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 CO₂, respective graphs for other substances are qualitatively the same.

Table 1: p-T ranges covered in the subcritical domain.

Substance
, MPa
, K

min
max
min
max

Ar
0.70
3.40
100.00
110.78/141.69

N₂
0.20
2.20
70.000
83.626/117.41

CO₂
1.00
6.00
220.00
233.03/295.13

CH₄
0.10
3.40
100.00
111.51/181.05

At the lowest/highest pressure of the domain.

Table 2: p-T ranges covered in the transcritical domain.

Substance
, MPa
, K

min
max
min
max

Ar
3.40
10.0
100.00
141.69

N₂
2.20
7.00
70.000
117.41

CO₂
6.00
15.0
220.00
295.13

CH₄
3.40
10.0
100.00
181.05

Figure 1: The domains of integration for CO₂: (red line): the boiling line, (blue line): the melting line, C.P.: the critical point, and T.P.: the triple point.

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).

Table 3: Number of initial/boundary data-points of and and of .

Substance
Initial/boundary values
Speed of sound

At lowest
At saturation
Subcritical
Transcritical

Ar
30
18
150
105

N₂
30
20
165
120

CO₂
30
20
165
135

CH₄
30
28
225
105

Table 4: The Lagrange polynomial orders used in the subcritical domain.

Substance
Interpolation with respect to
Derivation with respect to

p
T
p
T

Ar
9
14, 15
9
2, 11

N₂
10
14, 15
10
2, 11

CO₂
10
14, 15
10
2, 11

CH₄
14
14, 15
14
2, 11

Polynomial of 15th degree is used only when the temperature range is extended.

Polynomial of 2nd degree is used for 3 lowest temperatures.

Polynomial of 11th degree is used for 12 highest temperatures.

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).

Table 5: The Lagrange polynomial orders used in the transcritical domain.

Substance
Interpolation with respect to
Derivation with respect to

Ar
7
14
—
2, 11

N₂
8
14
—
2, 11

CO₂
9
14
—
2, 11

CH₄
7
14
—
2, 11

Table 6: The average step-size and number of steps taken in the subcritical domain.

Substance
Average step-size, MPa
Number of steps taken

IBVP-1
IBVP-2
IBVP-1
IBVP-2

Ar
0.073
0.01
37
270

N₂
0.056
0.01
36
200

CO₂
0.114
0.01
44
500

CH₄
0.065
0.001
51
3300

Calculated as ()/number of steps taken.

Constant step-size is used across the domain.

Table 7: The average step-size and number of steps taken in the transcritical domain.

Substance
Average step-size, MPa
Number of steps taken

IBVP-1
IBVP-2
IBVP-1
IBVP-2

Ar
0.471
0.440
14
15

N₂
0.369
0.343
13
14

CO₂
0.529
0.500
17
18

CH₄
0.471
0.413
14
16

Initial values obtained from the IBVP-1 method.

Initial values obtained from the IBVP-2 method.

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.

Table 8: The average absolute deviation (AAD) in the subcritical domain.

Substance
AAD, %

IBVP-1
IBVP-2
IBVP-1
IBVP-2

Ar
0.0001
0.0001
0.0084
0.0042

N₂
0.0003
0.0003
0.0066
0.0064

CO₂
0.0004
0.0002
0.0903
0.0175

CH₄
0.0010
0.0010
0.0260
0.0154

Average
0.0005
0.0004
0.0328
0.0109

Figure 2: Relative deviation of calculated of Ar versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.

Figure 3: Relative deviation of calculated of Ar versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.

Figure 4: Relative deviation of calculated of N₂ versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.

Figure 5: Relative deviation of calculated of N₂ versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1, □  □  □: IBVP-2.

Figure 6: Relative deviation of calculated of CO₂ versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.

Figure 7: Relative deviation of calculated of CO₂ versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.

Figure 8: Relative deviation of calculated of CH₄ versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.

Figure 9: Relative deviation of calculated of CH₄ versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.

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]).

Table 9: The average absolute deviation (AAD) in the transcritical domain.

Substance
AAD, %

IBVP-1
IBVP-2
IBVP-1
IBVP-2

Ar
0.0003
0.0003
0.0418
0.0367

N₂
0.0005
0.0016
0.0346
0.1078

CO₂
0.0029
0.0032
0.2815
0.2264

CH₄
0.0016
0.0046
0.0955
0.2178

Average
0.0013
0.0024
0.1134
0.1472

Table 10: The uncertainty of the reference data.

Substance
Uncertainty, %

Ar
0.02
2.0
1.0

N₂
0.02
0.8
1.0

CO₂
0.04
1.5
0.5

CH₄
0.03
1.0
0.3

Average
0.03
1.3
0.7

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.

Table 11: Influence of the uncertainty of the initial and boundary values of on AAD.

Substance
AAD, %

IBVP-1
IBVP-2
IBVP-1
IBVP-2

Ar
0.0199
0.0199
0.0085
0.0043

N₂
0.0200
0.0200
0.0066
0.0064

CO₂
0.0398
0.0160
0.0904
0.0036

CH₄
0.0299
0.0299
0.0260
0.0154

Average
0.0274
0.0215
0.0329
0.0074

Table 12: Influence of the uncertainty of the initial and boundary values of on AAD.

Substance
AAD, %

IBVP-1
IBVP-2
IBVP-1
IBVP-2

Ar
0.0076
0.0076
2.0425
2.0427

N₂
0.0018
0.0018
0.8087
0.8087

CO₂
0.0061
0.0007
1.5357
0.6034

CH₄
0.0028
0.0028
1.0133
1.0134

Average
0.0046
0.0032
1.3501
1.1171

Table 13: Influence of the uncertainty of the speed of sound on AAD.

Substance
AAD, %

IBVP-1
IBVP-2
IBVP-1
IBVP-2

Ar
0.0042
0.0042
0.0090
0.0058

N₂
0.0038
0.0038
0.0067
0.0069

CO₂
0.0030
0.0030
0.0905
0.0183

CH₄
0.0019
0.0019
0.0261
0.0154

Average
0.0032
0.0032
0.0331
0.0116

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.

References

J. P. M. Trusler, Physical Acoustics and Metrology of Fluids, Adam Hilger, Bristol, UK, 1991.

M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, UK, 1987.

M. J. Muringer, N. J. Trappeniers, and S. N. Biswas, “The effect of pressure on the sound velocity and density of toluene and n-heptane up to 2600 bar,” Physics and Chemistry of Liquids, vol. 14, no. 4, pp. 273–296, 2007.

T. F. Sun, P. J. Kortbeek, N. J. Trappeniers, and S. N. Biswas, “Acoustic and thermodynamic properties of benzene and cyclohexane as a function of pressure and temperature,” Physics and Chemistry of Liquids, vol. 16, no. 3, pp. 163–178, 1987.

T. F. Sun, J. A. Schouten, and S. N. Biswas, “Determination of the thermodynamic properties of liquid methanol from 203 to 263 K and up to 280 MPa from speed of sound measurements,” Berichte der Bunsengesellschaft für physikalische Chemie, vol. 94, no. 4, pp. 528–534, 1990.

T. F. Sun, J. A. Schouten, and S. N. Biswas, “Determination of the thermodynamic properties of liquid ethanol from 193 to 263 K and up to 280 MPa from speed-of-sound measurements,” International Journal of Thermophysics, vol. 12, no. 2, pp. 381–395, 1991.

J. P. Petitet, R. Tufeu, and B. Le Neindre, “Determination of the thermodynamic properties of water from measurements of the speed of sound in the temperature range 251.15–293.15 K and the pressure range 0.1–350 MPa,” International Journal of Thermophysics, vol. 4, no. 1, pp. 35–50, 1983.

T. S. Khasanshin, A. P. Shchamialiou, and O. G. Poddubskij, “Thermodynamic properties of heavy n-alkanes in the liquid state: n-dodecane,” International Journal of Thermophysics, vol. 24, no. 5, pp. 1277–1289, 2003.

G. Scalabrin, M. G. Concion, P. Marchi, and S. Lago, “Development of a fundamental equation of state for the liquid region of a pure fluid from speed of sound measurements. Application to water,” Experimental Thermal and Fluid Science, vol. 31, no. 6, pp. 539–549, 2007.

G. Benedetto, R. M. Gavioso, P. A. G. Albo, S. Lago, D. M. Ripa, and R. Spagnolo, “Speed of sound in pure water at temperatures between 274 and 394 K and at pressures up to 90 MPa,” International Journal of Thermophysics, vol. 26, no. 6, pp. 1667–1680, 2005.

M. Bijedić and N. Neimarlija, “Thermodynamic properties of liquids from speed of sound measurements,” International Journal of Thermodynamics, vol. 15, no. 2, pp. 61–68, 2012.

J.-P. Berrut and L. N. Trefethen, “Barycentric Lagrange interpolation,” SIAM Review, vol. 46, no. 3, pp. 501–517, 2004.

T. E. Hull, W. H. Enright, and K. R. Jackson, “User's guide for DVERK—a subroutine for solving non-stiff ODEs,” Tech. Rep. 100, University of Toronto, Toronto, Canada, 1976.

C. Tegeler, R. Span, and W. Wagner, “A new equation of state for argon covering the fluid region for temperatures from the melting line to 700 K at pressures up to 1000 MPa,” Journal of Physical and Chemical Reference Data, vol. 28, no. 3, pp. 779–850, 1999.

R. Span, E. W. Lemmon, R. T. Jacobsen, W. Wagner, and A. Yokozeki, “A reference quality thermodynamic property formulation for nitrogen,” Journal of Physical and Chemical Reference Data, vol. 29, no. 6, pp. 1361–1433, 2000.

R. Span and W. Wagner, “A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa,” Journal of Physical and Chemical Reference Data, vol. 25, no. 6, pp. 1509–1596, 1996.

U. Setzmann and W. Wagner, “A new equation of state and tables of thermodynamic properties for methane covering the range from the melting line to 625 K at pressures up to 1000 MPa,” Journal of Physical and Chemical Reference Data, vol. 20, no. 6, pp. 1061–1155, 1991.