This notebook contains course material from CBE20255 by Jeffrey Kantor (jeff at nd.edu); the content is available on Github. The text is released under the CC-BY-NC-ND-4.0 license, and code is released under the MIT license.

# PVT Computations for Non-ideal Gases¶

## Problem Statement¶

Commerical $CO_2$ cartridges have a large number of uses including bicycle tire and life jacket inflators, soda dispensers, compressed gas cleaners for electronic devices, and gas powered guns. They can be purchased at low-cost from sports store and on-line.

A popular size of $CO_2$ cartridge holds 12g in an internal volume of 17.6 cm$^3$. Estimate the pressure inside the cartridge under hot (40$^\circ$C = 104$^\circ$F = 313.15 K) conditions.

In [22]:
# problem data
n = 12/44                   # gram-moles
V = 0.0176/n                # specific volume in liters/gmol
T = 40 + 273.15             # Kelvin

# gas constant
R = 0.08206                 # gas constant liter-atm per kelvin-gmol

# critical properties of CO2
Tcritical = 304.25          # Kelvin
Pcritical = 72.9            # atm
acentric_factor  = 0.225   # Pitzer acentric factor

print('molar specific volume =', round(V,4), 'liters/gmol')
print('molar density =', round(1/V, 2), 'gmol/liter')

# a dictionary to store results
predictions = {}

molar specific volume = 0.0645 liters/gmol
molar density = 15.5 gmol/liter


## Estimating Pressure¶

### Ideal Gas Law¶

Our first attempt at computing the cartridge pressure will use the ideal gas law. Written in terms of molar specific volume

$$\hat{V} = \frac{V}{n}$$

the ideal gas law is given by

$$P \hat{V} = R T$$

Here we solve for $P$ using the problem data. Before going further, decide if the results of this calculation seem realistic to you.

In [3]:
P_ideal_gas = R*T/V

predictions['Ideal Gas'] = P_ideal_gas

print('CO2 Pressure (Ideal Gas) =', round(P_ideal_gas, 1), 'atm')
print('CO2 Pressure (Ideal Gas) =', round(14.696*P_ideal_gas, 1), 'psia')

CO2 Pressure (Ideal Gas) = 398.2 atm
CO2 Pressure (Ideal Gas) = 5851.9 psia


### Compressiblity Charts¶

Compressibility is defined as

$$z = \frac{P\hat{V}}{RT}$$

which, by definition, has a value $z=1$ for an ideal gas. Real gases display a significant deviation from 1, generally increasing with increasing pressure as the gas molecules get closer together, and decreasing with temperature.

The principle of corresponding states, developed by Johannes Diderik van der Waals in 1873, correlates the compressibility of real gases to the critical temperature $T_c$ and critical pressure $P_c$. For this purpose, the reduced temperature $T_r$ and reduced pressure $P_r$ are defined as

\begin{align*} T_r & = \frac{T}{T_c}\\ P_r & = \frac{P}{P_c} \end{align*}

A compressibility chart presents the averaged compressibility measured for a number of species. Given values of the reduced temperature and pressure, $T_r$ and $P_r$, one locates the value of compressibility factor $z$, then use the above relationship to compute molar volume $\hat{V}$.

For the case where $\hat{V}$ is one of the known variables it is convenient to define an ideal reduced volume. The first step is to define an ideal critical volume

$$\hat{V}^{ideal}_c = \frac{RT_c}{P_c}$$

The reason for the superscript 'ideal' is because this is a ficti

Then the ideal reduced volume $\hat{V}^{ideal}_r$ is

$$\hat{V}^{ideal}_r = \frac{\hat{V}}{\hat{V}_c} = \frac{P_c\hat{V}}{RT_c}$$

The compressibility chart is augmented with contours of constant $\hat{V}^{ideal}_r$. Given values for $\hat{V}^{ideal}_r$ and either $P_r$ or $T_r$, a corresponding value of $z$ is located. The compressibility definition is then used to solve for the unknown variable.

In [23]:
Tr = T/Tcritical
Vr = Pcritical*V/(R*Tcritical)

print('reduced temperature =', round(Tr,2))
print('ideal reduced volume =', round(Vr,2))

reduced temperature = 1.03
ideal reduced volume = 0.19

In [5]:
# look up values (https://pubs.acs.org/doi/pdf/10.1021/ie50467a036)
z = 0.28

# compute pressure
P_compressibility = z*R*T/V

predictions['Compressibility'] = P_compressibility

print('CO2 Pressure =', round(P_compressibility, 2), 'atm')
print('CO2 Pressure =', round(14.696*P_compressibility, 2), 'psia')

CO2 Pressure = 111.5 atm
CO2 Pressure = 1638.54 psia


While we still can't be sure this estimate is accurate, we can definitely see the estimate based on the ideal gas law is completley off-base. The reason is that the cartridge is filled to a point where the contents are close to the critical point on the phase diagram.

### Virial Model¶

The virial expansion was first proposed by a Kamerlingh Onnes, a physicist who won the Nobel prize in 1911 for his work on superconductivity and liquid helium. (He also coined the work 'enthalpy'). The key idea is to create a infinite series approximation for compressibility in the form

$$\frac{P\hat{V}}{RT} = A + \frac{B}{\hat{V}} + \frac{C}{\hat{V}^2} + \frac{D}{\hat{V}^3}\cdots$$

where $A$, $B$, $C$ are temperature dependent and known as the first, second, and third virial coefficients, respectively. The case $A=1$ and $B=C=D=\cdots = 0$ corresponds to an ideal gas.

A commonly used version of this expansion is to assume $A=1$, $B(T)$, and $C=D=\cdots = 0$, with the further approximation $\hat{V} = \frac{RT}{P}$. Then

$$\frac{P\hat{V}}{RT} = 1 + \frac{BP}{RT}$$

This can be simplified to read

$$P= \frac{RT}{\hat{V}-B}$$

The temperature dependent value of $B$ is estimated by

\begin{align*} B_0 & = 0.083 - \frac{0.422}{T_r^{1.6}} \\ B_1 & = 0.139 - \frac{0.172}{T_r^{4.2}} \\ B & = \frac{RT_c}{P_c}\left(B_0 + \omega B_1\right) \end{align*}

$\omega$ is the Pitzer acentric factor, values of which are tabulated in standard sources for chemical data.

In [18]:
Tr = T/Tcritical
B0 = 0.083 - 0.422/Tr**1.6
B1 = 0.139 - 0.172/Tr**4.2
B = (R*Tcritical/Pcritical)*(B0 + acentric_factor*B1)
print('B =', B)

P_virial = R*T/(V - B)

predictions['Virial'] = P_virial

print('CO2 Pressure =', round(P_virial, 2), 'atm')
print('CO2 Pressure =', round(14.696*P_virial, 2), 'psia')

B = -0.11061596372197276
CO2 Pressure = 146.72 atm
CO2 Pressure = 2156.13 psia


### van der Waals Equation of State¶

The van der Waals equation of state has the form

$$P = \frac{RT}{\hat{V} - b} - \frac{a}{\hat{V}^2}$$

where values for the coefficients $a$ and $b$ are determined by the critical point temperature and pressure

$$a = \frac{27 R^2 T_c^2}{64 P_c}\qquad b = \frac{RT_c}{8P_c}$$

The parameter $a$ accounts for long-range attractive forces acting between molecules, and $b$ for short-range repulsive forces.

In [19]:
a = 27*R**2*Tcritical**2/(64*Pcritical)
b = R*Tcritical/(8*Pcritical)
print('a =', a, '  b =', b)

P_vdw = R*T/(V-b) - a/V**2

predictions['van der Waals'] = P_vdw

print('CO2 Pressure =', round(P_vdw, 2), 'atm')
print('CO2 Pressure =', round(14.696*P_vdw, 2), 'psia')

a = 3.6072850418404214   b = 0.042809936556927296
CO2 Pressure = 316.73 atm
CO2 Pressure = 4654.72 psia


### Soave-Redlich-Kwong Equation of State¶

The Soave-Redlick-Kwong equation of state is one of most widely used equations of states, and proven to be applicable to a wide variety of systems. The general expression

$$P = \frac{RT}{\hat{V}-b} - \frac{\alpha a}{\hat{V}(\hat{V}+b)}$$

where the parameters are given by

\begin{align*} a & = 0.42747 \frac{(RT_c)^2}{P_c} \\ b & = 0.08664 \frac{RT_c}{P_c} \\ m & = 0.48508 + 1.55171\omega - 0.1561\omega^2\\ T_r & = \frac{T}{T_c}\\ \alpha & = \left[1 + m(1 - \sqrt{T_r})\right]^2 \end{align*}
In [8]:
from math import sqrt

a = 0.42747*(R*Tcritical)**2/Pcritical
b = 0.08664*R*Tcritical/Pcritical
m = 0.48508 + 1.5517*acentric_factor - 0.1561*acentric_factor**2
Tr = T/Tcritical
alpha = (1 + m*(1-sqrt(Tr)))**2

P_srk = R*T/(V-b) - alpha*a/V/(V+b)

predictions['SRK'] = P_srk

print('CO2 Pressure =', round(P_srk, 2), 'atm')
print('CO2 Pressure =', round(14.696*P_srk, 2), 'psia')

CO2 Pressure = 150.24 atm
CO2 Pressure = 2207.97 psia


### Reference Data from NIST Webbook¶

The National Institute of Standards and Technology (NIST) maintains a web site devoted to the distribution of standard reference data. The NIST Chemistry WebBook is an excellent source of carefully curated data on over 7000 organic and small inorganic compounds.

In [9]:
print('density =', 1/V, 'mol/liter')

density = 15.495867768595039 mol/liter

In [20]:
P_nist = 108.13  # atm

predictions['NIST Webbook'] = P_nist

print('CO2 Pressure =', round(P_nist, 2), 'atm')
print('CO2 Pressure =', round(14.696*P_nist, 2), 'atm')

CO2 Pressure = 108.13 atm
CO2 Pressure = 1589.08 atm


### Comparison of Estimates¶

In [21]:
print('{0:15s}  {1:9s}   {2:5s}'.format('EOS', 'Pressure', 'Error'))

for key,val in predictions.items():
err = 100*(val-P_nist)/P_nist
print('{0:15s}  {1:5.1f} atm   {2:5.1f}%'.format(key, val, err))


EOS              Pressure    Error
Ideal Gas        398.2 atm   268.3%
Compressibility  111.5 atm     3.1%
Virial           146.7 atm    35.7%
van der Waals    316.7 atm   192.9%
SRK              150.2 atm    38.9%
NIST             108.1 atm     0.0%
NIST Webbook     108.1 atm     0.0%


## Non-Ideal Behavior¶

How does the $CO_2$ pressure depend on temperature in the range from 30C to 50C?

## Discussion Points¶

1. Why is the ideal gas law so wrong?
2. Why do the compressibility charts work so well in this case? Do you expect them to work well in others?
3. How would the calculations change if we were given pressure and temperature, and asked to compute molar volume?