High Molecular Weight Petroleum PseudocomponentsΒΆ

Thermo is a general phase equilibrium engine, and if the user provides enough properties for the components, there is no issue adding your own components. In this basic example below, a made-up extended gas analysis is used to specify a gas consisting of the standard real components and three heavier fractions, C10+, C12+ and C15+.

A bare minimum of basic properties are estimated using the Kesler-Lee method (1976), and the estimated fraction molecular weights are turned into atomic compositions. The heat capacities of each pseudocomponent is found with the similarity variable concept of Lastovka and Shaw (2013) based on atomic composition.

This example ends with calculating a flash at 270 Kelvin and 1 bar.

[1]:
from math import log, exp
import numpy as np
from scipy.constants import psi
from thermo import *
from chemicals import *

def Tc_Kesler_Lee_SG_Tb(SG, Tb):
    r'''Estimates critical temperature of a hydrocarbon compound or petroleum
    fraction using only its specific gravity and boiling point, from
    [1]_ as presented in [2]_.

    .. math::
        T_c = 341.7 + 811.1SG + [0.4244 + 0.1174SG]T_b
        + \frac{[0.4669 - 3.26238SG]10^5}{T_b}

    Parameters
    ----------
    SG : float
        Specific gravity of the fluid at 60 degrees Farenheight [-]
    Tb : float
        Boiling point the fluid [K]

    Returns
    -------
    Tc : float
        Estimated critical temperature [K]

    Notes
    -----
    Model shows predictions for Tc, Pc, MW, and omega.
    Original units in degrees Rankine.

    Examples
    --------
    Example 2.2 from [2]_, but with K instead of R.

    >>> Tc_Kesler_Lee_SG_Tb(0.7365, 365.555)
    545.0124354151242

    References
    ----------
    .. [1] Kesler, M. G., and B. I. Lee. "Improve Prediction of Enthalpy of
       Fractions." Hydrocarbon Processing (March 1976): 153-158.
    .. [2] Ahmed, Tarek H. Equations of State and PVT Analysis: Applications
       for Improved Reservoir Modeling. Gulf Pub., 2007.
    '''
    Tb = 9/5.*Tb # K to R
    Tc = 341.7 + 811.1*SG + (0.4244 + 0.1174*SG)*Tb + ((0.4669 - 3.26238*SG)*1E5)/Tb
    Tc = 5/9.*Tc # R to K
    return Tc

def Pc_Kesler_Lee_SG_Tb(SG, Tb):
    r'''Estimates critical pressure of a hydrocarbon compound or petroleum
    fraction using only its specific gravity and boiling point, from
    [1]_ as presented in [2]_.

    .. math::
        \ln(P_c) = 8.3634 - \frac{0.0566}{SG} - \left[0.24244 + \frac{2.2898}
        {SG} + \frac{0.11857}{SG^2}\right]10^{-3}T_b
        + \left[1.4685 + \frac{3.648}{SG} + \frac{0.47227}{SG^2}\right]
        10^{-7}T_b^2-\left[0.42019 + \frac{1.6977}{SG^2}\right]10^{-10}T_b^3

    Parameters
    ----------
    SG : float
        Specific gravity of the fluid at 60 degrees Farenheight [-]
    Tb : float
        Boiling point the fluid [K]

    Returns
    -------
    Pc : float
        Estimated critical pressure [Pa]

    Notes
    -----
    Model shows predictions for Tc, Pc, MW, and omega.
    Original units in degrees Rankine and psi.

    Examples
    --------
    Example 2.2 from [2]_, but with K instead of R and Pa instead of psi.

    >>> Pc_Kesler_Lee_SG_Tb(0.7365, 365.555)
    3238323.346840464

    References
    ----------
    .. [1] Kesler, M. G., and B. I. Lee. "Improve Prediction of Enthalpy of
       Fractions." Hydrocarbon Processing (March 1976): 153-158.
    .. [2] Ahmed, Tarek H. Equations of State and PVT Analysis: Applications
       for Improved Reservoir Modeling. Gulf Pub., 2007.
    '''
    Tb = 9/5.*Tb # K to R
    Pc = exp(8.3634 - 0.0566/SG - (0.24244 + 2.2898/SG + 0.11857/SG**2)*1E-3*Tb
    + (1.4685 + 3.648/SG + 0.47227/SG**2)*1E-7*Tb**2
    -(0.42019 + 1.6977/SG**2)*1E-10*Tb**3)
    Pc = Pc*psi
    return Pc

def MW_Kesler_Lee_SG_Tb(SG, Tb):
    r'''Estimates molecular weight of a hydrocarbon compound or petroleum
    fraction using only its specific gravity and boiling point, from
    [1]_ as presented in [2]_.

    .. math::
        MW = -12272.6 + 9486.4SG + [4.6523 - 3.3287SG]T_b + [1-0.77084SG
        - 0.02058SG^2]\left[1.3437 - \frac{720.79}{T_b}\right]\frac{10^7}{T_b}
        + [1-0.80882SG + 0.02226SG^2][1.8828 - \frac{181.98}{T_b}]
        \frac{10^{12}}{T_b^3}

    Parameters
    ----------
    SG : float
        Specific gravity of the fluid at 60 degrees Farenheight [-]
    Tb : float
        Boiling point the fluid [K]

    Returns
    -------
    MW : float
        Estimated molecular weight [g/mol]

    Notes
    -----
    Model shows predictions for Tc, Pc, MW, and omega.
    Original units in degrees Rankine.

    Examples
    --------
    Example 2.2 from [2]_, but with K instead of R and Pa instead of psi.

    >>> MW_Kesler_Lee_SG_Tb(0.7365, 365.555)
    98.70887589833501

    References
    ----------
    .. [1] Kesler, M. G., and B. I. Lee. "Improve Prediction of Enthalpy of
       Fractions." Hydrocarbon Processing (March 1976): 153-158.
    .. [2] Ahmed, Tarek H. Equations of State and PVT Analysis: Applications
       for Improved Reservoir Modeling. Gulf Pub., 2007.
    '''
    Tb = 9/5.*Tb # K to R
    MW = (-12272.6 + 9486.4*SG + (4.6523 - 3.3287*SG)*Tb + (1.-0.77084*SG - 0.02058*SG**2)*
    (1.3437 - 720.79/Tb)*1E7/Tb + (1.-0.80882*SG + 0.02226*SG**2)*
    (1.8828 - 181.98/Tb)*1E12/Tb**3)
    return MW

def omega_Kesler_Lee_SG_Tb_Tc_Pc(SG, Tb, Tc=None, Pc=None):
    r'''Estimates accentric factor of a hydrocarbon compound or petroleum
    fraction using only its specific gravity and boiling point, from
    [1]_ as presented in [2]_. If Tc and Pc are provided, the Kesler-Lee
    routines for estimating them are not used.

    For Tbr > 0.8:
    .. math::
        \omega = -7.904 + 0.1352K - 0.007465K^2 + 8.359T_{br}
        + ([1.408-0.01063K]/T_{br})

    Otherwise:
    .. math::
        \omega = \frac{-\ln\frac{P_c}{14.7} - 5.92714 + \frac{6.09648}{T_{br}}
        + 1.28862\ln T_{br} - 0.169347T_{br}^6}{15.2518 - \frac{15.6875}{T_{br}}
         - 13.4721\ln T_{br} + 0.43577T_{br}^6}

        K = \frac{T_b^{1/3}}{SG}

        T_{br} = \frac{T_b}{T_c}

    Parameters
    ----------
    SG : float
        Specific gravity of the fluid at 60 degrees Farenheight [-]
    Tb : float
        Boiling point the fluid [K]
    Tc : float, optional
        Estimated critical temperature [K]
    Pc : float, optional
        Estimated critical pressure [Pa]

    Returns
    -------
    omega : float
        Acentric factor [-]

    Notes
    -----
    Model shows predictions for Tc, Pc, MW, and omega.
    Original units in degrees Rankine and psi.

    Examples
    --------
    Example 2.2 from [2]_, but with K instead of R and Pa instead of psi.

    >>> omega_Kesler_Lee_SG_Tb_Tc_Pc(0.7365, 365.555, 545.012, 3238323.)
    0.306392118159797

    References
    ----------
    .. [1] Kesler, M. G., and B. I. Lee. "Improve Prediction of Enthalpy of
       Fractions." Hydrocarbon Processing (March 1976): 153-158.
    .. [2] Ahmed, Tarek H. Equations of State and PVT Analysis: Applications
       for Improved Reservoir Modeling. Gulf Pub., 2007.
    '''
    if Tc is None:
        Tc = Tc_Kesler_Lee_SG_Tb(SG, Tb)
    if Pc is None:
        Pc = Pc_Kesler_Lee_SG_Tb(SG, Tb)
    Tb = 9/5.*Tb # K to R
    Tc = 9/5.*Tc # K to R
    K = Tb**(1/3.)/SG
    Tbr = Tb/Tc
    if Tbr > 0.8:
        omega = -7.904 + 0.1352*K - 0.007465*K**2 + 8.359*Tbr + ((1.408-0.01063*K)/Tbr)
    else:
        omega = ((-log(Pc/101325.) - 5.92714 + 6.09648/Tbr + 1.28862*log(Tbr)
        - 0.169347*Tbr**6) / (15.2518 - 15.6875/Tbr - 13.4721*log(Tbr) +0.43577*Tbr**6))
    return omega

[2]:
# Basic composition and names. All pure component properties are obtained from Chemicals and Thermo.
pure_constants = ChemicalConstantsPackage.constants_from_IDs(
    ['water', 'hydrogen', 'helium', 'nitrogen', 'carbon dioxide', 'hydrogen sulfide', 'methane',
    'ethane', 'propane', 'isobutane', 'n-butane', 'isopentane', 'n-pentane', 'hexane',
    'heptane', 'octane', 'nonane'])
pure_fractions = [.02, .00005, .00018, .009, .02, .002, .82, .08, .031,
                  .009, .0035, .0033, .0003, .0007, .0004, .00005, .00002]
[3]:
pseudo_names = ['C10-C11', 'C12-C14', 'C15+']
pseudo_carbon_numbers = [10.35, 12.5, 16.9]
pseudo_SGs = [.73, .76, .775] # Specific gravity values are based of the alkane series
pseudo_Tbs = [447, 526, 589]

# Using the estimation methods defined earlier, we obtain some critical properties
pseudo_Tcs = [Tc_Kesler_Lee_SG_Tb(SG, Tb) for SG, Tb in zip(pseudo_SGs, pseudo_Tbs)]
pseudo_Pcs = [Pc_Kesler_Lee_SG_Tb(SG, Tb) for SG, Tb in zip(pseudo_SGs, pseudo_Tbs)]
pseudo_MWs = [MW_Kesler_Lee_SG_Tb(SG, Tb) for SG, Tb in zip(pseudo_SGs, pseudo_Tbs)]
pseudo_omegas = [omega_Kesler_Lee_SG_Tb_Tc_Pc(SG, Tb) for SG, Tb in zip(pseudo_SGs, pseudo_Tbs)]
# Estimate the hydroen counts
hydrogen_counts = [(MW - C*periodic_table.C.MW)/periodic_table.H.MW
                  for C, MW in zip(pseudo_carbon_numbers, pseudo_MWs)]
# Get the atomic compositions
pseudo_atoms = [{'C': C, 'H': H} for C, H in zip(pseudo_carbon_numbers, hydrogen_counts)]
# Calculate the similarity variable of each species
similarity_variables = [similarity_variable(atoms=atoms) for atoms in pseudo_atoms]

pseudo_fractions = [.0003, .00015, .00005]
[4]:
pseudos = ChemicalConstantsPackage(names=pseudo_names, MWs=pseudo_MWs, Tbs=pseudo_Tbs,
                                   atomss=pseudo_atoms,
                                  Tcs=pseudo_Tcs, Pcs=pseudo_Pcs, omegas=pseudo_omegas,
                                  similarity_variables=similarity_variables)
# Add the pure components and the pseudocomponents to create a new package of constant values
# which will be used by the phase and flash objects
constants = pure_constants + pseudos
# Obtain the temperature and pressure dependent objects
properties = PropertyCorrelationsPackage(constants=constants)
# This is the feed composition
zs = normalize(pure_fractions + pseudo_fractions)
T = 270 # K
P = 1e5 # bar
[5]:
kijs = np.zeros((constants.N, constants.N)).tolist() # kijs left as zero in this example
eos_kwargs = dict(Tcs=constants.Tcs, Pcs=constants.Pcs, omegas=constants.omegas, kijs=kijs)

# The API SRK equation of state is used, but other cubic equations of state can be uesd instead
gas = CEOSGas(APISRKMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)
liq = CEOSLiquid(APISRKMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)
liq2 = CEOSLiquid(APISRKMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)
phase_list = [gas, liq, liq]

# Set up the three phase flash engine
flashN = FlashVLN(constants, properties, liquids=[liq, liq2], gas=gas)
[6]:
# Do the flash, and get some properties
res = flashN.flash(T=T, P=P, zs=zs)
res.phase_count, res.gas_beta, res.liquids_betas
[6]:
(3, 0.9827041561275567, [0.016838840039984482, 0.00045700383245886604])
[7]:
res.H(), res.Cp_mass(), res.MW(), res.gas.mu(), res.gas.k()
[7]:
(-1961.5089633224945,
 1989.3915447041693,
 19.675910651652533,
 1.0015480015807005e-05,
 0.027074133060773903)
[8]:
res.heaviest_liquid.rho_mass(), res.lightest_liquid.rho_mass()
[8]:
(769.2525386053419, 599.2086838769081)