Computational DePriester Chart

depriester.py

from math import log, exp

tolerance = 1e-5

coefficients = {
    "methane": (-292860.0, 8.24450, -0.89510, 59.8465, 0),
    "ethene": (-600076.875, 7.90595, -0.84677, 42.94594),
    "ethane": (-687248.25, 7.90699, -0.88600, 49.02654),
    "propene": (-923484.6875, 7.71725, -0.87871, 47.67624),
    "propane": (-970688.5625, 7.15059, -0.76984, 0, 6.90224),
    "ibutane": (-1166846.0, 7.72668, -0.92213),
    "nbutane": (-1280557.0, 7.994986, -0.96455),
    "ipentane": (-1481583.0, 7.58071, -0.93159),
    "npentane": (-1524891.0, 7.33129, -0.89143),
    "nhexane": (-1778901.0, 6.96783, -0.84634),
    "nheptane": (-2013803.0, 6.52914, -0.79543)
}


def comp_derivative(h, function, x, *args):
    return (function(x + h, *args) - function(x - h, *args)) / (2 * h)


def newton_raphson(h, function, initial, *args):
    derivative = comp_derivative(h, function, initial, *args)
    new = initial - function(initial, *args) / derivative
    if abs(initial - new) <= tolerance:
        return new
    return newton_raphson(h, function, new, *args)


def correlation(t, p, at1, at6, ap1, ap2=0.0, ap3=0.0):
    """ T (R); P (psia) """
    return exp(at1 / t**2 + at6 + ap1 * log(p) + ap2 / p**2 + ap3 / p)


def celsius_to_rankine(t):
    return (t + 273.15) * 9 / 5


def rankine_to_celsius(t):
    return (t - 491.67) * 5 / 9


def kpa_to_psia(p):
    return p / 6.8947572932


def psia_to_kpa(p):
    return 6.8947572932 * p


def bubble_point_from_pressure_iter(t, p, data):
    return sum([v * correlation(t, p, *coefficients[k])
                for k, v in data.items()]) - 1


def bubble_point_from_pressure(p, data, guess=600, h=0.01):
    """ P (kPa) """
    return rankine_to_celsius(newton_raphson(
        h, bubble_point_from_pressure_iter, guess, kpa_to_psia(p), data))


def bubble_point_from_temp_iter(p, t, data):
    return sum([v * correlation(t, p, *coefficients[k])
                for k, v in data.items()]) - 1


def bubble_point_from_temp(t, data, guess=50, h=0.01):
    """ T (C) """
    return psia_to_kpa(newton_raphson(
        h, bubble_point_from_temp_iter, guess, celsius_to_rankine(t), data))


def dew_point_from_pressure_iter(t, p, data):
    return sum([v / correlation(t, p, *coefficients[k])
                for k, v in data.items()]) - 1


def dew_point_from_pressure(p, data, guess=600, h=0.01):
    """ P (kPa) """
    return rankine_to_celsius(newton_raphson(
        h, dew_point_from_pressure_iter, guess, kpa_to_psia(p), data))


def dew_point_from_temp_iter(p, t, data):
    return sum([v / correlation(t, p, *coefficients[k])
                for k, v in data.items()]) - 1


def dew_point_from_temp(t, data, guess=50, h=0.01):
    """ T (C) """
    return psia_to_kpa(newton_raphson(
        h, dew_point_from_temp_iter, guess, celsius_to_rankine(t), data))


data = {"nbutane": 0.4,
        "npentane": 0.25,
        "nhexane": 0.2,
        "nheptane": 0.15
        }
# Usage:
#   bubble_point_from_temp(t)     (C)
#   bubble_point_from_pressure(t) (kPa)
#
#   dew_point_from_temp(t)        (C)
#   dew_point_from_pressure(t)    (kPa)

# e.g:

pressure = 405.3
bubble_temp = bubble_point_from_pressure(pressure, data)
print(
    "bubble temperature: {:.2f}C; bubble pressure {:.2f} kPa".format(
        bubble_temp,
        pressure))

temperature = 65
bubble_pressure = bubble_point_from_temp(
    temperature, data)
print(
    "bubble temperature: {:.2f}C; bubble pressure {:.2f} kPa".format(
        temperature,
        bubble_pressure))

dew_temp = dew_point_from_pressure(pressure, data)
print(
    "dew temperature: {:.2f}C; dew pressure {:.2f} kPa".format(
        dew_temp,
        pressure))

dew_pressure = dew_point_from_temp(
    temperature, data)
print(
    "dew temperature: {:.2f}C; dew pressure {:.2f} kPa".format(
        temperature,
        dew_pressure))