Calculates point until which ideal gas calculations are accurate to within one percent for various gases. https://physikmix.lewinb.net for more!

gasapproximation.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat May  4 11:36:15 2019

Results for:

    He  - 1.015 MPa (+/- 0.99%)
    CO2 -   191.754 kPa (+/- 1%)
    N2  - 1.002 MPa (+/- 0.99%)
    H2  - 1.413 MPa (+/- 1%)
    Air -   974.579 kPa (+/- 0.99%)

(= pressure at which ideal and van-der-Waals equations differ by less than
or exactly 1%)
"""

import matplotlib.pyplot as plt
import numpy as np

import math
from collections import namedtuple

R = 8.314472

# Gas parameters

Gas = namedtuple('gas', ('a', 'b'))


# Common gases

GASES = {
        "He": Gas(3.45e-3, 23.7e-6),
        "CO2": Gas(363.7e-3, 42.7e-6),
        "N2": Gas(140.8e-3, 39.1e-6),
        "H2": Gas(24.7e-3, 26.6e-6),
        "Air": Gas(135.8e-3, 36.4e-6),
}

def ideal_pressure(V, n=1, T=273):
    """Calculates the pressure of an ideal gas."""
    return (n * R * T)/V

def vdw_pressure(gas, V, n=1, T=273):
    """Calculates the approx. pressure of an ideal gas."""
    return (n * R * T)/(V - n*gas.b) - (n*n*gas.a)/(V*V)

def inaccuracy(gas, V, **kwargs):
    """Returns difference between ideal and approximate pressure in percent."""
    return abs(100 * (1 - vdw_pressure(gas, V, **kwargs)/ideal_pressure(V, **kwargs)))

def find_limit_volume(gas, Vmax=1, Vmin=None, max_diff=0.01, **kwargs):
    """Finds limit volume at which inaccuracy (rel. difference) between
    ideal and approximate pressure is just less than max_diff, by bisecting
    the range between 0 and Vmax.
    
    Also takes keyword arguments T, n.
    """
    INACCURACY_ACCURACY = 0.01 # percent
    Vmin = Vmin if Vmin else Vmax/2
    # We expect vmax_inacc < vmin_inacc.
    vmax_inacc = inaccuracy(gas, Vmax, **kwargs)
    vmin_inacc = inaccuracy(gas, Vmin, **kwargs)

    if vmin_inacc < 1:
        # Continue search with Vmin as upper limit.
        return find_limit_volume(gas, Vmax=Vmin, **kwargs)
    elif vmax_inacc > 1:
        # Give up if the upper inaccuracy is already too bad.
        return -1
    # In this case, the limit point of ~1 percent inaccuracy is located between
    # Vmin and Vmax.
    elif vmin_inacc > 1 > vmax_inacc:
        # Check if Vmin or Vmax are already sufficiently close
        # (within +- INACCURACY_ACCURACY percent)
        if vmin_inacc - 1 < INACCURACY_ACCURACY:
            return Vmin
        if 1 - vmax_inacc < INACCURACY_ACCURACY:
            return Vmax
        # If neither is close, continue bisection by cutting the current
        # (Vmin, Vmax) interval in two halves. We always return the upper
        # result.
        return max(find_limit_volume(gas, Vmin=Vmin, Vmax=(Vmax+Vmin)/2),
                   find_limit_volume(gas, Vmin=(Vmax+Vmin)/2, Vmax=Vmax))

def find_limit_pressure(gas, **kwargs):
    limit_volume = find_limit_volume(gas, **kwargs)
    vdw = vdw_pressure(gas, limit_volume, **kwargs)
    return (vdw, inaccuracy(gas, limit_volume))

def plot_pressure_figure(gas, name='', Vmin=1e-4, Vmax=.002, T=273, n=1):
    fig = plt.figure(dpi=100)
    x = np.logspace(math.log10(Vmin), math.log10(Vmax), 100)
    ideal = np.vectorize(lambda v: ideal_pressure(v, T=T, n=n))(x)
    vdw = np.vectorize(lambda v: vdw_pressure(gas, v, T=T, n=n))(x)
    inacc = np.vectorize(lambda v: inaccuracy(gas, v, T=T, n=n))(x)
    
    axes1 = fig.add_subplot(111)
    axes1.grid(True)
    axes1.set_yscale('log')
    axes1.set_xscale('log')
    axes1.set_ylabel('p/Pa')
    axes1.set_xlabel('V/m3')
    axes1.set_title(name)
    
    axes2 = fig.add_subplot(111, sharex=axes1, frameon=False)
    axes2.set_ylabel('%')
    axes2.yaxis.set_label_position('right')
    axes2.yaxis.tick_right()

    lines = []
    lines += axes1.plot(x, ideal, 'b', label='ideal')
    lines += axes1.plot(x, vdw, 'g', label='vdw')
    lines += axes2.plot(x, inacc, 'r', label='inacc %')

    fig.legend()
    plt.show()

for name, gas in GASES.items():
    (limit_pressure, acc) = find_limit_pressure(gas)
    print('{}: limit pressure at 273K: {:.0f} Pa (within {:.2f}%)'.format(name, limit_pressure, acc))
    
    Tcrit = 8/27 * gas.a/(R*gas.b)
    Vcrit = 3*gas.b
    plot_pressure_figure(gas, name=(name + ' at critical point'), T=Tcrit, Vmin=Vcrit/2, Vmax=100*Vcrit)