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):

import yaeos

?yaeos.models

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:

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.

?PengRobinson76

This can be done with all the listed models of yaeos.

from yaeos import NRTL, PengRobinson76, PengRobinson78, SoaveRedlichKwong, RKPR


?RKPR

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.

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: yaeos is compatible with numpy. The functions that receives arrays arguments are able to handle numpy arrays.

• We create three variables:

  • n array with the number of moles of each component

  • T with the temperature of the system in Kelvin

  • P with the pressure of the system in bar

• call the lnphi_pt method inside the previously instantiated model variable

  model.<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, the model object (PengRobinson76 instance). lnphi_pt is a method of the PengRobinson76 class 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:

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:

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:

?model.volume

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:

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:

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 ]

?model.lnphi_pt

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.

lnphi, derivatives = model.lnphi_pt(
    moles=n, pressure=P, temperature=T, dn=True
)

print("ln(phi): ", lnphi)

derivatives

ln(phi):  [-0.24307099 -0.2526632 ]

{'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:

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:

derivatives["dt"]

array([0.00247786, 0.00257727])

Check the signature of each property to see which derivatives you can calculate. For example pressure.

?model.pressure

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]\).