# Density and Heat Capacity of Liquids from Speed of Sound

Authors: Muhamed Bijedić; Sabina Begić
Journal: Journal of Thermodynamics
 (2016)

---

## 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%.


---

## Body

## 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 p
-
T 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 rectangularp
-
T domain, the following set of equations may be used (Benedetto et al. [10]):(1)
∂
ρ
∂
p
T
=
T
α
p
2
c
p
+
1
u
2
(2)
∂
c
p
∂
p
T
=
-
T
ρ
α
p
2
+
∂
α
p
∂
T
p
,where p is the pressure, T is the temperature, u is the speed of sound, ρ is the density, c
p is the specific heat capacity at constant pressure, and α
p is the thermal expansion coefficient:(3)
α
p
=
-
1
ρ
∂
ρ
∂
T
p
.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]):(4)
∂
ρ
∂
p
ϕ
=
T
α
p
2
c
p
+
1
u
2
+
∂
ρ
∂
T
p
∂
T
∂
p
ϕ
(5)
∂
c
p
∂
p
ϕ
=
-
T
ρ
α
p
2
+
∂
α
p
∂
T
p
+
∂
c
p
∂
T
p
∂
T
∂
p
ϕ
,where ϕ represents a path connecting the points with temperature(6)
T
j
=
T
m
i
n
+
j
n
T
m
a
x
-
T
m
i
n
j
=
0,1
,
2
,
…
,
n
,where T
m
i
n is the lowest temperature of the range and T
m
a
x is the highest temperature at observed pressure (e.g., the saturation temperature). Using ϕ instead of T as a path of integration enables one to shape the domain arbitrarily.The set of (4) and (5) may be solved simultaneously for ρ and c
p 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 c
p 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 thep
-
T ranges covered, are given in Tables 1 and 2 for the subcritical and transcritical domain, respectively. These domains are also presented graphically in p
-
T coordinates at Figure 1. While this graph is given for CO2, respective graphs for other substances are qualitatively the same.Table 1
p-T ranges covered in the subcritical domain.
Substance
p, MPa
T, K
min
max
min
max(
1
)
Ar
0.70
3.40
100.00
110.78/141.69
N2
0.20
2.20
70.000
83.626/117.41
CO2
1.00
6.00
220.00
233.03/295.13
CH4
0.10
3.40
100.00
111.51/181.05
(
1
)At the lowest/highest pressure of the domain.Table 2
p-T ranges covered in the transcritical domain.
Substance
p, MPa
T, K
min
max
min
max
Ar
3.40
10.0
100.00
141.69
N2
2.20
7.00
70.000
117.41
CO2
6.00
15.0
220.00
295.13
CH4
3.40
10.0
100.00
181.05Figure 1
The domains of integration for CO2: (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]):(7)
T
j
=
T
m
i
n
+
T
m
a
x
2
+
T
m
a
x
-
T
m
i
n
2
cos
⁡
j
π
14
j
=
14,13,12
,
…
,
0
.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 c
p, as well as u, are specified along the lowest pressure p
0 (see Table 3) at temperatures calculated from (7). The temperature derivatives of ρ, α
p, and c
p and the pressure derivatives of T are estimated from the Lagrangian interpolation polynomial (see Table 4). After that, there is enough information to calculate the pressure derivatives of ρ and c
p from (4) and (5), respectively. Then, new values of ρ and c
p are calculated by the RKV method at the pressure p
1
=
p
0
+
Δ
p, where Δ
p 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 u 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 c
p and of u.
Substance
Initial/boundary values
Speed of sound
At lowestp
At saturation
Subcritical
Transcritical
Ar
30
18
150
105
N2
30
20
165
120
CO2
30
20
165
135
CH4
30
28
225
105Table 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(
1
)
9
2(
2
), 11(
3
)
N2
10
14, 15
10
2, 11
CO2
10
14, 15
10
2, 11
CH4
14
14, 15
14
2, 11
(
1
)Polynomial of 15th degree is used only when the temperature range is extended.
(
2
)Polynomial of 2nd degree is used for 3 lowest temperatures.
(
3
)Polynomial of 11th degree is used for 12 highest temperatures.In the second method (IBVP-2) the initial values ofρ and c
p, as well as u, are also specified along the lowest pressure p
0 (see Table 3) at temperatures calculated from (7). The temperature derivatives of ρ and α
p are estimated from the Lagrangian interpolation polynomial (see Table 4). After that, there is enough information to calculate the pressure derivatives of ρ and c
p from (1) and (2), respectively. Then, new values of ρ and c
p are calculated by the RKV method at the pressure p
1
=
p
0
+
Δ
p, where Δ
p 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 c
p 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 u 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 c
p 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
p
T
p
T
Ar
7
14
—
2, 11
N2
8
14
—
2, 11
CO2
9
14
—
2, 11
CH4
7
14
—
2, 11Table 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(
1
)
IBVP-2(
2
)
IBVP-1
IBVP-2
Ar
0.073
0.01
37
270
N2
0.056
0.01
36
200
CO2
0.114
0.01
44
500
CH4
0.065
0.001
51
3300
(
1
)Calculated as (p
max
⁡
–
p
min
⁡)/number of steps taken.
(
2
)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(
1
)
IBVP-2(
2
)
IBVP-1
IBVP-2
Ar
0.471
0.440
14
15
N2
0.369
0.343
13
14
CO2
0.529
0.500
17
18
CH4
0.471
0.413
14
16
(
1
)Initial values obtained from the IBVP-1 method.
(
2
)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, %
ρ
c
p
IBVP-1
IBVP-2
IBVP-1
IBVP-2
Ar
0.0001
0.0001
0.0084
0.0042
N2
0.0003
0.0003
0.0066
0.0064
CO2
0.0004
0.0002
0.0903
0.0175
CH4
0.0010
0.0010
0.0260
0.0154
Average
0.0005
0.0004
0.0328
0.0109Figure 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 calculatedc
p of Ar versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.Figure 4
Relative deviation of calculatedρ of N2 versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.Figure 5
Relative deviation of calculatedc
p of N2 versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1, □  □  □: IBVP-2.Figure 6
Relative deviation of calculatedρ of CO2 versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.Figure 7
Relative deviation of calculatedc
p of CO2 versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.Figure 8
Relative deviation of calculatedρ of CH4 versus temperature at the highest pressure of the subcritical domain: ○  ○  ○: IBVP-1; □  □  □: IBVP-2.Figure 9
Relative deviation of calculatedc
p of CH4 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 Table9). 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, %
ρ
c
p
IBVP-1
IBVP-2
IBVP-1
IBVP-2
Ar
0.0003
0.0003
0.0418
0.0367
N2
0.0005
0.0016
0.0346
0.1078
CO2
0.0029
0.0032
0.2815
0.2264
CH4
0.0016
0.0046
0.0955
0.2178
Average
0.0013
0.0024
0.1134
0.1472Table 10
The uncertainty of the reference data.
Substance
Uncertainty, %
ρ
c
p
u
Ar
0.02
2.0
1.0
N2
0.02
0.8
1.0
CO2
0.04
1.5
0.5
CH4
0.03
1.0
0.3
Average
0.03
1.3
0.7Influence 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 Table10) 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, %
ρ
c
p
IBVP-1
IBVP-2
IBVP-1
IBVP-2
Ar
0.0199
0.0199
0.0085
0.0043
N2
0.0200
0.0200
0.0066
0.0064
CO2
0.0398
0.0160
0.0904
0.0036
CH4
0.0299
0.0299
0.0260
0.0154
Average
0.0274
0.0215
0.0329
0.0074Table 12
Influence of the uncertainty of the initial and boundary values ofc
p on AAD.
Substance
AAD, %
ρ
c
p
IBVP-1
IBVP-2
IBVP-1
IBVP-2
Ar
0.0076
0.0076
2.0425
2.0427
N2
0.0018
0.0018
0.8087
0.8087
CO2
0.0061
0.0007
1.5357
0.6034
CH4
0.0028
0.0028
1.0133
1.0134
Average
0.0046
0.0032
1.3501
1.1171Table 13
Influence of the uncertainty of the speed of sound on AAD.
Substance
AAD, %
ρ
c
p
IBVP-1
IBVP-2
IBVP-1
IBVP-2
Ar
0.0042
0.0042
0.0090
0.0058
N2
0.0038
0.0038
0.0067
0.0069
CO2
0.0030
0.0030
0.0905
0.0183
CH4
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).



---

*Source: 2035704-2016-01-20.xml*