# Properties of Pure Fluids



```python
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:


```python
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:

$$
\left( \frac{\partial P}{\partial v} \right)_T = 0 \\
\left( \frac{\partial^2 P}{\partial v^2} \right)_T = 0 \;,
$$

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:


```python
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:


```python
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:


```python
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
\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$:


```python
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:


```python
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.

872a495d02bd3ab29bda336d801c99739c2c9d14
ipynb
Jupyter Notebook
book/content/properties-pure/properties-pure.ipynb
kyleniemeyer/computational-thermo
3f0d1d4a6d4247ac3bf3b74867411f2090c70cbd
CC-BY-4.0 BSD-3-Clause
2020-04-01T05:52:06.000Z
2022-03-27T20:25:59.000Z
book/content/properties-pure/properties-pure.ipynb
kyleniemeyer/computational-thermo
3f0d1d4a6d4247ac3bf3b74867411f2090c70cbd
CC-BY-4.0 BSD-3-Clause
2020-04-28T04:02:05.000Z
2020-04-29T17:49:52.000Z
book/content/properties-pure/properties-pure.ipynb
kyleniemeyer/computational-thermo
3f0d1d4a6d4247ac3bf3b74867411f2090c70cbd
CC-BY-4.0 BSD-3-Clause
2020-04-03T14:52:24.000Z
2022-03-29T02:29:43.000Z
Qwen/Qwen-72B
1. YES
2. YES