Instantiate a model
Which models are available?
In yaeos, as in thermodynamics, are two kinds of models:
Residual Helmholtz free energy models (ArModel)
Excess Gibbs free energy models (GeModel)
yaeos provides a set of predefined models that will change in the future. You can check the list of available models in the API documentation or simply executing the following command on an interactive Python session (Jupyter Notebook, Google Colab, etc):
[1]:
import yaeos
?yaeos.models
Type: module
String form: <module 'yaeos.models' from '/home/salvador/code/yaeos/python/yaeos/models/__init__.py'>
File: ~/code/yaeos/python/yaeos/models/__init__.py
Docstring:
Models module.
Yaeos models module. This module provides the following submodules:
- excess_gibbs: Excess Gibbs energy models
- NRTL: non-random two-liquid model
- UNIFACVLE: Original UNIFAC VLE model
- UNIQUAC: UNIversal QUAsiChemical Excess Gibbs free energy model
- residual_helmholtz: Residual Helmholtz energy models
- Cubic EoS:
- PengRobinson76: Peng-Robinson model (1976)
- PengRobinson78: Peng-Robinson model (1978)
- SoaveRedlichKwong: Soave-Redlich-Kwong model
- RKPR: RKPR model
- Mixing rules: mixing rules for cubic EoS
- QMR: cuadratic mixing rule
- MHV: modified Huron-Vidal mixing rule
All this models can be imported directly from yaeos with the names in the list. For example let’s import and instantiate a PengRobinson76 model:
[2]:
from yaeos import PengRobinson76
# Critical constants
Tc = [320, 375] # critical temperatures [K]
Pc = [30, 45] # critical pressures [bar]
w = [0.0123, 0.045] # acentric factors [-]
model = PengRobinson76(Tc, Pc, w)
There we have instantiated a PengRobinson76 model with a binary mixture and stored the instance in a variable called model. You can check the signature to instantiate any model with the ? operator in an interactive Python session.
[3]:
?PengRobinson76
Init signature:
PengRobinson76(
critical_temperatures,
critical_pressures,
acentric_factors,
mixrule: yaeos.models.residual_helmholtz.cubic_eos.mixing_rules.CubicMixRule = None,
) -> None
Docstring:
Peng-Robinson 1976 cubic equation of state.
Parameters
----------
critical_temperatures : array_like
Critical temperatures vector [K]
critical_pressures : array_like
Critical pressures vector [bar]
acentric_factors : array_like
Acentric factors vector
mixrule : CubicMixRule, optional
Mixing rule object. If no provided the quadratric mixing rule (QMR)
with zero for kij and lij parameters is set, by default None
Attributes
----------
nc : int
Number of components
critical_temperatures : array_like
Critical temperatures vector [K]
critical_pressures : array_like
Critical pressures vector [bar]
acentric_factors : array_like
Acentric factors vector
id : int
EoS identifier
mixrule : CubicMixRule
Mixing rule object
Example
-------
.. code-block:: python
from yaeos import PengRobinson76
tc = [190.56, 305.32] # Critical temperatures [K]
pc = [45.99, 48.72] # Critical pressures [bar]
w = [0.0115, 0.0985] # Acentric factors
pr76 = PengRobinson76(tc, pc, w)
File: ~/code/yaeos/python/yaeos/models/residual_helmholtz/cubic_eos/cubic_eos.py
Type: ABCMeta
Subclasses:
This can be done with all the listed models of yaeos.
[4]:
from yaeos import NRTL, PengRobinson76, PengRobinson78, SoaveRedlichKwong, RKPR
?RKPR
Init signature:
RKPR(
critical_temperatures,
critical_pressures,
acentric_factors,
critical_z,
k=None,
delta_1=None,
mixrule: yaeos.models.residual_helmholtz.cubic_eos.mixing_rules.CubicMixRule = None,
) -> None
Docstring:
RKPR cubic equation of state.
Parameters
----------
critical_temperatures : array_like
Critical temperatures vector [K]
critical_pressures : array_like
Critical pressures vector [bar]
acentric_factors : array_like
Acentric factors vector
critical_z : array_like
Critical compressibility factor vector
k : array_like, optional
k parameter, by default None
delta_1 : array_like, optional
delta_1 parameter, by default None
mixrule : CubicMixRule, optional
Mixing rule object. If no provided the quadratric mixing rule (QMR)
with zero for kij and lij parameters is set, by default None
Attributes
----------
nc : int
Number of components
critical_temperatures : array_like
Critical temperatures vector [K]
critical_pressures : array_like
Critical pressures vector [bar]
acentric_factors : array_like
Acentric factors vector
zc : array_like
Critical compressibility factor vector
id : int
EoS identifier
mixrule : CubicMixRule
Mixing rule object
Example
-------
.. code-block:: python
from yaeos import RKPR
tc = [190.56, 305.32] # Critical temperatures [K]
pc = [45.99, 48.72] # Critical pressures [bar]
w = [0.0115, 0.0985] # Acentric factors
zc = [0.27, 0.28] # Critical compressibility factor
rkpr = RKPR(tc, pc, w, zc)
File: ~/code/yaeos/python/yaeos/models/residual_helmholtz/cubic_eos/cubic_eos.py
Type: ABCMeta
Subclasses:
What can we do with an instantiated model?
Residual Helmholtz free energy models (ArModel)
We will get later on what we can do with the Excess Gibbs free energy models instances.
By now we will use the PengRobinson76 instantiated in the variable model to calculate properties of the mixture. For example the fugacity coefficients of the mixture’s components at a given composition, temperature and pressure.
[5]:
import numpy as np
n = np.array([7.0, 3.0]) # moles of each component [mol]
P = 2.0 # pressure [bar]
T = 300.0 # temperature [K]
# Calculate ln(phi) with the model (PengRobinson76)
ln_phi = model.lnphi_pt(moles=n, pressure=P, temperature=T)
# Use ln(phi) to calculate the fugacity
phi = np.exp(ln_phi)
z = n / np.sum(n)
fugacity = z * phi * P
# Print results
print("ln(phi) = ", ln_phi)
print("phi = ", phi)
print("fugacity = ", fugacity, "bar")
ln(phi) = [-0.03026809 -0.03078738]
phi = [0.9701854 0.96968172]
fugacity = [1.35825956 0.58180903] bar
Well, analyzing the previous code by parts we have:
import numpy:
yaeosis compatible with numpy. The functions that receives arrays arguments are able to handle numpy arrays.We create three variables:
narray with the number of moles of each componentTwith the temperature of the system in KelvinPwith the pressure of the system in bar
call the
lnphi_ptmethod inside the previously instantiatedmodelvariablemodel.<thing i want to evaluate>(arguments)
For the users not familiar with Python and object-oriented programming, the dot
.is used to access the attributes and methods of an object, in this case, themodelobject (PengRobinson76 instance).lnphi_ptis a method of thePengRobinson76class that calculates the natural logarithm of the fugacity coefficients of the components of the mixture at a given composition, temperature and pressure. In this case the syntax is:model.lnphi_pt(n, P, T)
Use the ln(phi) returned values to calculate the fugacity coefficients
First we calculate the fugacity coefficients.
python phi = np.exp(lnphi)Which means: \(\phi_i = e^{\ln(\phi_i)}\)Then, we calculate the mole fractions of the mixture.
python z = n / np.sum(n)Which means: \(z_i = \frac{n_i}{\sum n_i}\)Then we use the phi variable to calculate fugacities:
fugacity = z * phi * P
Which means: \(f = z \; \phi \; P\)
ArModel properties listing
All the properties that can be calculated with an ArModel instance can be found on the API documentation (Core Module):
https://ipqa-research.github.io/yaeos/page/python-api/core.html#yaeos.core.ArModel
There you will find all the properties that can be calculated with an ArModel, their signatures and a brief description of what they do. Also, an example to calculate each one is provided.
Some examples that we can calculate now:
[6]:
n = np.array([7.0, 3.0]) # moles of each component [mol]
P = 15.0 # pressure [bar]
V = 20.0 # volume [L]
T = 300.0 # temperature [K]
# Evaluate pressure
pressure = model.pressure(moles=n, volume=V, temperature=T)
# Evaluate ln(phi) but function of temperature and volume
ln_phi = model.lnphi_vt(moles=n, volume=V, temperature=T)
# Evaluate volume
volume = model.volume(moles=n, pressure=P, temperature=T)
And we can print the results:
[7]:
print("Pressure: ", pressure, " bar")
print("ln_phi: ", ln_phi)
print("Volume: ", volume, " L")
Pressure: 10.30361077545868 bar
ln_phi: [-0.16224585 -0.1669435 ]
Volume: 12.068572622671574 L
If you are familiar with thermodynamics, you probably may think:
"Well, which volume root is solving when `model.volume` is evaluated"
And you will be right, we didn’t tell you yet. You can get the documentation of all these properties on the API documentation (the link provided above) or by doing:
[8]:
?model.volume
Signature:
model.volume(
moles,
pressure: float,
temperature: float,
root: str = 'stable',
) -> float
Docstring:
Calculate volume given pressure and temperature [L].
Parameters
----------
moles : array_like
Moles number vector [mol]
pressure : float
Pressure [bar]
temperature : float
Temperature [K]
root : str, optional
Volume root, use: "liquid", "vapor" or "stable", by default
"stable"
Returns
-------
float
Volume [L]
Example
-------
.. code-block:: python
import numpy as np
from yaeos import PengRobinson76
tc = np.array([320.0, 375.0]) # critical temperatures [K]
pc = np.array([45.0, 60.0]) # critical pressures [bar]
w = np.array([0.0123, 0.045]) # acentric factors
model = PengRobinson76(tc, pc, w)
# Evaluating stable root volume
# will print: 23.373902973572587
print(model.volume(np.array([5.0, 5.6]), 10.0, 300.0))
# Liquid root volume (not stable)
# will print: 0.8156388756398074
print(model.volume(np.array([5.0, 5.6]), 10.0, 300.0, "liquid"))
File: ~/code/yaeos/python/yaeos/core.py
Type: method
If you check the API documentation you have an example of calculating “liquid” volume roots, “vapor” volume roots or “stable” volume roots, but we will explain it here too.
As you can see in the volume documentation above. There is another argument that we didn’t use: root. This argument is set as “stable” by default, but we can change it to “vapor” or “liquid” if we want. Let’s see which was the stable root in the previous evaluation:
[9]:
n = np.array([7.0, 3.0]) # moles of each component [mol]
P = 15.0 # pressure [bar]
T = 300.0 # temperature [K]
# Stable root
print(
"Stable: ", model.volume(moles=n, pressure=P, temperature=T, root="stable")
)
# Liquid root
print(
"Liquid: ", model.volume(moles=n, pressure=P, temperature=T, root="liquid")
)
# Vapor root
print(
"Vapor: ", model.volume(moles=n, pressure=P, temperature=T, root="vapor")
)
Stable: 12.068572622671574
Liquid: 1.2405831819647626
Vapor: 12.068572622671574
As we can see, the stable root was the vapor root. But changing the root argument we can also check the value of the liquid root if we want.
Let back to the lnphi_pt method, there we specify temperature and pressure, but again, at which volume root is solving the model? The answer is the same as the volume method, the stable root. But we can change it to:
[10]:
n = np.array([7.0, 3.0]) # moles of each component [mol]
P = 15.0 # pressure [bar]
T = 300.0 # temperature [K]
# Stable root
print(
"Stable: ",
model.lnphi_pt(moles=n, pressure=P, temperature=T, root="stable"),
)
# Liquid root
print(
"Liquid: ",
model.lnphi_pt(moles=n, pressure=P, temperature=T, root="liquid"),
)
# Vapor root
print(
"Vapor: ", model.lnphi_pt(moles=n, pressure=P, temperature=T, root="vapor")
)
Stable: [-0.24307099 -0.2526632 ]
Liquid: [-0.04887053 -0.3742866 ]
Vapor: [-0.24307099 -0.2526632 ]
[11]:
?model.lnphi_pt
Signature:
model.lnphi_pt(
moles,
pressure: float,
temperature: float,
root: str = 'stable',
dt: bool = False,
dp: bool = False,
dn: bool = False,
) -> Union[numpy.ndarray, tuple[numpy.ndarray, dict]]
Docstring:
Calculate fugacity coefficent given pressure and temperature.
Calculate :math:`ln \phi_i(n,P,T)` and its derivatives with respect to
temperature, pressure and moles number.
Parameters
----------
moles : array_like
Moles number vector [mol]
pressure : float
Pressure [bar]
temperature : float
Temperature [K]
root : str, optional
Volume root, use: "liquid", "vapor" or "stable", by default
"stable"
dt : bool, optional
Calculate temperature derivative, by default False
dp : bool, optional
Calculate pressure derivative, by default False
dn : bool, optional
Calculate moles derivative, by default False
Returns
-------
Union[np.ndarray, tuple[np.ndarray, dict]]
:math:`ln \phi_i(n,P,T)` vector or tuple with
:math:`ln \phi_i(n,P,T)` vector and derivatives dictionary if any
derivative is asked
Example
-------
.. code-block:: python
import numpy as np
from yaeos import PengRobinson76
tc = np.array([320.0, 375.0]) # critical temperatures [K]
pc = np.array([45.0, 60.0]) # critical pressures [bar]
w = np.array([0.0123, 0.045]) # acentric factors
model = PengRobinson76(tc, pc, w)
# Evaluating ln_phi only
# will print: [-0.10288733 -0.11909807]
print(model.lnphi_pt([5.0, 5.6], 10.0, 300.0))
# Asking for derivatives
# will print:
# (
# array([-0.10288733, -0.11909807]),
# {'dt': array([0.00094892, 0.00108809]), 'dp': None, 'dn': None}
# )
print(model.lnphi_pt([5.0, 5.6], 10.0, 300.0, dt=True)
File: ~/code/yaeos/python/yaeos/core.py
Type: method
Obtaining derivatives of the properties
We can obtain some derivatives of the properties that we have been calculating. For example, the fugacity coefficients derivatives with respect to the number of moles of the components of the mixture.
[12]:
lnphi, derivatives = model.lnphi_pt(
moles=n, pressure=P, temperature=T, dn=True
)
print("ln(phi): ", lnphi)
derivatives
ln(phi): [-0.24307099 -0.2526632 ]
[12]:
{'dt': None,
'dp': None,
'dn': array([[-2.56426727e-06, 5.98329030e-06],
[ 5.98329030e-06, -1.39610107e-05]])}
There we set the dn argument to True to calculate the derivatives with respect to the number of moles of the components of the mixture. When we specify any of the derivatives arguments as True we will get another return of the method. The first will be the property evaluated as usual and the second is a Python dictionary with the calculated derivatives specified as True.
If we specify all the derivatives as True:
[13]:
lnphi, derivatives = model.lnphi_pt(
moles=n, pressure=P, temperature=T, dt=True, dp=True, dn=True
)
print("ln(phi): ", lnphi)
print(derivatives)
ln(phi): [-0.24307099 -0.2526632 ]
{'dt': array([0.00247786, 0.00257727]), 'dp': array([-0.01786595, -0.01925548]), 'dn': array([[-2.56426727e-06, 5.98329030e-06],
[ 5.98329030e-06, -1.39610107e-05]])}
There we get all the derivatives of lnphi. Of course, we can access the dictionary with the keys:
[14]:
derivatives["dt"]
[14]:
array([0.00247786, 0.00257727])
Check the signature of each property to see which derivatives you can calculate. For example pressure.
[15]:
?model.pressure
Signature:
model.pressure(
moles,
volume: float,
temperature: float,
dv: bool = False,
dt: bool = False,
dn: bool = False,
) -> Union[float, tuple[float, dict]]
Docstring:
Calculate pressure given volume and temperature [bar].
Calculate :math:`P(n,V,T)` and its derivatives with respect to
volume, temperature and moles number.
Parameters
----------
moles : array_like
Moles number vector [mol]
volume : float
Volume [L]
temperature : float
Temperature [K]
dv : bool, optional
Calculate volume derivative, by default False
dt : bool, optional
Calculate temperature derivative, by default False
dn : bool, optional
Calculate moles derivative, by default False
Returns
-------
Union[float, tuple[float, dict]]
Pressure or tuple with Presure and derivatives dictionary if any
derivative is asked [bar]
Example
-------
.. code-block:: python
import numpy as np
from yaeos import PengRobinson76
tc = np.array([320.0, 375.0]) # critical temperatures [K]
pc = np.array([45.0, 60.0]) # critical pressures [bar]
w = np.array([0.0123, 0.045]) # acentric factors
model = PengRobinson76(tc, pc, w)
# Evaluating pressure only
# will print: 16.011985733846956
print(model.pressure(np.array([5.0, 5.6]), 2.0, 300.0))
# Asking for derivatives
# will print:
# (
# 16.011985733846956,
# {'dv': None, 'dt': np.float64(0.7664672352866752), 'dn': None}
# )
print(model.pressure(np.array([5.0, 5.6]), 2.0, 300.0, dt=True))
File: ~/code/yaeos/python/yaeos/core.py
Type: method
Yaeos units
If you were paying attention to the previous examples, you may have noticed that the units of yaeos are defined according to the ideal gas constant R with units:
\(R = 0.08314462618 \frac{bar \; L}{K \; mol}\)
Because of that, pressure must be specified as bar, volume as liters, temperature as Kelvin and number of moles as moles. The returns of the properties will be in the same units as the inputs. Energetic properties as enthalpy, entropy, etc will have unit of \([bar \; L]\), which equals to centijoules \([cJ]\).