Properties of Pure Fluids#

import numpy as np
import cantera as ct
from pint import UnitRegistry
ureg = UnitRegistry()
Q_ = ureg.Quantity

We can obtain critical properties of fluids using Cantera:

f = ct.Nitrogen()
print(f'Critical temperature: {f.critical_temperature: .2f} K')
print(f'Critical pressure: {f.critical_pressure: .2e} Pa')

density_crit = f.critical_density
mol_weight = f.mean_molecular_weight
z_crit = (
    f.critical_pressure * mol_weight / 
    (density_crit * ct.gas_constant * f.critical_temperature)
    )
print(f'Critical compressibility factor: {z_crit: .2f}')

Critical temperature:  126.20 K
Critical pressure:  3.40e+06 Pa
Critical compressibility factor:  0.29

van der Waals equation of state#

The van der Waals equation of state is

\[ P = \frac{RT}{v - b} - \frac{a}{v^2} \;, \]

where \(a\) and \(b\) are substance-specific constants.

We can find these constants by applying the contraints on equations of state along the critical isotherms:

\[\begin{split} \left( \frac{\partial P}{\partial v} \right)_T = 0 \\ \left( \frac{\partial^2 P}{\partial v^2} \right)_T = 0 \;, \end{split}\]

which both apply at the critical point. These will give us closed-form expressions for \(a\) and \(b\) based on the critical properties.

Let’s use SymPy to help us find these expressions, by taking these derivatives and applying the constraints at the critical point:

import sympy
sympy.init_printing(use_latex='mathjax')

R, Tc, vc, Pc, a, b = sympy.symbols('R, T_c, v_c, P_c, a, b', real=True)

equation_state = (R*Tc)/(vc - b) - (a/vc**2)

# first derivative with respect to v
first_deriv = sympy.diff(equation_state, vc)
display(sympy.Eq(first_deriv, 0))

second_deriv = sympy.diff(first_deriv, vc)
display(sympy.Eq(second_deriv, 0))

\[\displaystyle - \frac{R T_{c}}{\left(- b + v_{c}\right)^{2}} + \frac{2 a}{v_{c}^{3}} = 0\]

\[\displaystyle \frac{2 R T_{c}}{\left(- b + v_{c}\right)^{3}} - \frac{6 a}{v_{c}^{4}} = 0\]

Then, we can solve this system of two equations for our two unknowns:

sol = sympy.solve((
    first_deriv, 
    second_deriv, 
    ), (a, b), dict=True)

display(sol[0])

a_expr = sol[0][a]
b_expr = sol[0][b]

\[\displaystyle \left\{ a : \frac{9 R T_{c} v_{c}}{8}, \ b : \frac{v_{c}}{3}\right\}\]

This gives us the coefficients \(a\) and \(b\) as a function of \(T_{\text{crit}}\), \(P_{\text{crit}}\), and \(R\). However, while the critical temperature can be measured experimentally relatively easily, the critical specific volume is harder to measure accurately.

So, we can introduce the equation of state applied at the critical point,

\[ P_{\text{crit}} = \frac{R T_{\text{crit}}}{v_{\text{crit}} - b} - \frac{a}{v_{\text{crit}}^2} \;, \]

and then solve the system of three equations for \(a\), \(b\), and \(v_{\text{crit}}\) to eliminate the critical specific volume:

sol = sympy.solve((
    a - a_expr, b - b_expr, Pc - ((R*Tc / (vc-b)) - (a / vc**2))
    ), (a, b, vc), dict=True)
display(sol[0])

\[\displaystyle \left\{ a : \frac{27 R^{2} T_{c}^{2}}{64 P_{c}}, \ b : \frac{R T_{c}}{8 P_{c}}, \ v_{c} : \frac{3 R T_{c}}{8 P_{c}}\right\}\]

Thus, we can now determine the coefficients \(a\) and \(b\) based on measured values of the critical temperature and pressure.

We can now use the equation of state to get \(P = f(T,v)\), but what about determining \(v\) based on \(T\) and \(P\)? It turns out that the van der Waals equation of state gives a cubic function of \(v\), which when expressed in terms of reduced properties is

(27)#\[\begin{equation} v_r^3 - \left( \frac{8}{3} \frac{T_r}{P_r} + \frac{1}{3} \right) + \frac{3}{P_r} v_r - \frac{1}{P_r} = 0 \end{equation}\]

This means that we need to find the roots of this cubic polynomial to get specific volume. As a consequence, there will be three roots and potential solutions for \(v_r\). At minimum, we should only consider real roots; imaginary roots are not physically meaningful.

First, consider when \(T_r > 1.0\), such as for \(T_r = 2.5\) and \(P_r = 2.0\):

temp_r = 2.5
pres_r = 2.0
coeffs = [
    1.0, 
    -((8.0*temp_r / (3.0*pres_r)) + (1.0 / 3.0)),
    3.0 / pres_r,
    -1.0 / pres_r
    ]
roots = np.roots(coeffs)
print(f'Roots: {roots}')
# get only real roots
roots = roots.real[np.abs(roots.imag)<1e-5]

print(f'Number of real roots: {len(roots)}')

print(f'T_r = {temp_r: .2f}')
print(f'P_r = {pres_r: .2f}')

reduced_volumes = ','.join([f'{vol_r: .3f}' for vol_r in roots])
print(f'v_r = {reduced_volumes}')

Roots: [3.25277894+0.j         0.20694386+0.33299993j 0.20694386-0.33299993j]
Number of real roots: 1
T_r =  2.50
P_r =  2.00
v_r =  3.253

So, for temperatures above the critical temperature, there will be only one real value for reduced specific volume.

But, if \(T_r < 1\), we can get three real roots:

temp_r = 0.90
pres_r = 0.61
coeffs = [
    1.0, 
    -((8.0*temp_r / (3.0*pres_r)) + (1.0 / 3.0)),
    3.0 / pres_r,
    -1.0 / pres_r
    ]
roots = np.roots(coeffs)
roots = roots.real[np.abs(roots.imag)<1e-5]

print(f'Number of real roots: {len(roots)}')

print(f'T_r = {temp_r: .2f}')
print(f'P_r = {pres_r: .2f}')

reduced_volumes = ','.join([f'{vol_r: .3f}' for vol_r in roots])
print(f'v_r = {reduced_volumes}')

Number of real roots: 3
T_r =  0.90
P_r =  0.61
v_r =  2.640, 1.017, 0.610

In this case, we need to apply some additional knowledge to understand these roots and identify any non-physical values, such as considering the sign of \(\frac{\partial P_r}{\partial v_r}\) for these points.

Computational Thermodynamics

Properties of Pure Fluids

Contents

Properties of Pure Fluids#

van der Waals equation of state#