Bulk Properties

yaeos is a library based on Equation-of-State calculations. An Equation of State (EOS) is a mathematical model that describes the relationship between thermodynamic properties of a system, typically pressure (P), volume (V), and temperature (T). These equations are fundamental in physics, chemistry, and engineering, especially in thermodynamics and fluid mechanics.

EOS models are used to predict phase behavior, thermodynamic properties, and equilibrium conditions of gases, liquids, solids, and mixtures. They are widely applied in areas such as chemical engineering, petroleum engineering, and material science. From the ideal gas law, to cubic equations of state and SAFT, to more complex models, EOS are essential tools for the study and design of processes.

The classic way of defining an EOS is though a mathematical equation that express the pressure as a function of the volume and temperature. For example, the Van der Waals EOS is defined as:

\[P = \frac{RT}{v - b} - \frac{a}{v^2}\]

Where \(P\) is the pressure, \(v\) is the molar volume, \(T\) is the temperature, \(R\) is the gas constant, and \(a\) and \(b\) are parameters that depend on the substance. For mixtures, the \(a\) and \(b\) parameters are calculated from the pure components parameters using a mixing rule.

The modern way of defining an EOS is through a mathematical model that express the residual Helmholtz free energy as a function of the mole number vector, volume and temperature. This approach is more general and allows for a more flexible definition of the EOS in the context of implementing it in code. All thermodynamic properties can be calculated from the residual Helmholtz free energy and its derivatives. yaeos is based on this approach.

\[A^{r}(\vec{n}, V, T) = -\int_{\infty}^{V} \left(P(\vec{n},V,T) - \frac{nRT}{V} \right) dV\]

If you want to know how each bulk property is calculated from the residual Helmholtz free energy, please refer to the Fortran User Documentation

Example model

First, we are going the import the yaeos library and instantiate a Peng-Robinson (1976) model with two components. Later in this tutorial you will learn got to instantiate different models, but for now we will use this one as an example.

The important thing to understand is that all the EoS can evaluate the thermodynamic properties explained in this section.

import yaeos


# Model parameters
Tc = [320, 375]      # critical temperatures [K]
Pc = [30, 45]        # critical pressures [bar]
w = [0.0123, 0.045]  # acentric factors [-]

model = yaeos.PengRobinson76(Tc, Pc, w)

Properties function of Volume and Temperature \(\left(V, T \right)\)

We can obtain different bulk properties from a \(A^r(\vec{n}, V, T)\) model specifying the number of moles, volume and temperature. But also, derivatives could be obtained if asked for.

All not asked derivatives will be set as None. The derivatives are returned as a Python dictionary with the keys dv, dt, and dn, etc. respectively. The compositional derivatives are always returned as Numpy arrays.

The following properties can be calculated:

Pressure

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.pressure(n, V, T) # Calculate pressure [bar]

56.000067835576885

You can obtain the next derivatives:

\[\left(\frac{\partial P}{\partial V}\right)_{T, \vec{n}} \quad \left(\frac{\partial P}{\partial T}\right)_{V, \vec{n}} \quad \left(\frac{\partial P}{\partial n_i}\right)_{V, T, n_{j \neq i}}\]

Ask for derivatives as follows:

# Asking for all derivatives
pressure, derivatives = model.pressure(n, V, T, dv=True, dt=True, dn=True)

# Displaying results
print("Pressure: ", pressure)
print("Derivatives: ", derivatives)

Pressure:  56.000067835576885
Derivatives:  {'dv': -647.5628198115064, 'dt': 2.900736382506016, 'dn': array([78.64160081, 50.87096316])}

Logarithm of fugacity coefficients \(ln \, \phi(\vec{n},V,T)\)

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.lnphi_vt(n, V, T)

array([-1.12084851, -1.61086   ])

You can obtain the next derivatives:

\[\left(\frac{\partial ln \, \hat{\phi}_i}{\partial T}\right)_{P, \vec{n}} \quad \left(\frac{\partial ln \, \hat{\phi}_i}{\partial P}\right)_{T, \vec{n}} \quad \left(\frac{\partial ln \, \hat{\phi}_i}{\partial n_j}\right)_{P, T, n_{k \neq j}}\]

Ask for derivatives as follows:

# Asking for all the derivatives
lnphi, derivatives = model.lnphi_vt(n, V, T, dt=True, dp=True, dn=True)

# Displaying results
print("Logarithm of Fugacity Coefficients: ", lnphi)
print("Derivatives: ", derivatives)

Logarithm of Fugacity Coefficients:  [-1.12084851 -1.61086   ]
Derivatives:  {'dt': array([0.01380036, 0.01926589]), 'dp': array([-0.0129884 , -0.01470769]), 'dn': array([[-0.00904545,  0.00904545],
       [ 0.00904545, -0.00904545]])}

You can evaluate fugacity by simply:

import numpy as np


p = model.pressure(n, V, T)  # Calculate pressure [bar]

z = n / np.sum(n)            # Calculate mole fractions

fugacities = np.exp(lnphi) * z * p

print("Fugacities: ", fugacities)

Fugacities:  [9.12809679 5.59204872]

Residual Helmholtz free energy \(A^r(\vec{n},V,T)\)

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.helmholtz_residual_vt(n, V, T) # Calculate residual Helmholtz energy [bar L]

-519.8710586930237

You can obtain the next derivatives:

\[\left(\frac{\partial A^{r}_{(\vec{n}, V, T)}}{\partial V}\right)_{T, \vec{n}} \quad \left(\frac{\partial A^{r}_{(\vec{n}, V, T)}}{\partial T}\right)_{V, \vec{n}} \quad \left(\frac{\partial A^{r}_{(\vec{n}, V, T)}}{\partial n_i}\right)_{V, T, n_{j \neq i}} \quad\]

\[\left(\frac{\partial^2 A^{r}_{(\vec{n}, V, T)}}{\partial V \partial T}\right)_{ \vec{n}} \quad \left(\frac{\partial^2 A^{r}_{(\vec{n}, V, T)}}{\partial V^2}\right)_{T, \vec{n}} \quad \left(\frac{\partial^2 A^{r}_{(\vec{n}, V, T)}}{\partial T^2}\right)_{V, \vec{n}} \quad\]

\[\left(\frac{\partial^2 A^{r}_{(\vec{n}, V, T)}}{\partial n_i \partial n_j}\right)_{V, T, n_{k \neq i,j}} \quad \left(\frac{\partial^2 A^{r}_{(\vec{n}, V, T)}}{\partial V \partial n_i}\right)_{T, n_{j \neq i}} \quad \left(\frac{\partial^2 A^{r}_{(\vec{n}, V, T)}}{\partial T \partial n_i}\right)_{V, n_{j \neq i}} \quad\]

Ask for derivatives as follows:

# Asking for all the derivatives
Ar , derivatives = model.helmholtz_residual_vt(
    n, V, T, dv=True, dt=True, dn=True, dtv=True, dvn=True, dtn=True, dn2=True,
    dv2=True, dt2=True
)

# Displaying results
print("Residual Helmholtz free energy: ", Ar)
print("Derivatives: ", derivatives)

Residual Helmholtz free energy:  -519.8710586930237
Derivatives:  {'dt': 1.746237812013294, 'dv': 193.43381070442314, 'dn': array([-65.21921351, -77.44176037]), 'dtv': -2.069290120706016, 'dv2': 398.1289412715064, 'dt2': -0.0021897892913883494, 'dvn': array([-53.69821295, -25.9275753 ]), 'dtn': array([0.39595811, 0.36714748]), 'dn2': array([[6.8304639 , 3.90917869],
       [3.90917869, 1.27633637]])}

Residual Gibbs free energy \(G^r(\vec{n}, V, T)\)

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.gibbs_residual_vt(n, V, T) # Calculate residual Gibbs energy [bar L]

-713.3048693974467

You can obtain the next derivatives:

\[\left(\frac{\partial G^{r}_{(\vec{n}, V, T)}}{\partial T}\right)_{V, \vec{n}} \quad \left(\frac{\partial G^{r}_{(\vec{n}, V, T)}}{\partial V}\right)_{T, \vec{n}} \quad \left(\frac{\partial G^{r}_{(\vec{n}, V, T)}}{\partial n_j}\right)_{V, T, n_{i \neq j}}\]

Ask for derivatives as follows:

# Asking for all the derivatives
Gr, derivatives = model.gibbs_residual_vt(n, V, T, dt=True, dv=True, dn=True)

# Displaying results
print("Residual Gibbs Energy: ", Gr)
print("Derivatives: ", derivatives)

Residual Gibbs Energy:  -713.3048693974467
Derivatives:  {'dt': 3.81552793271931, 'dv': -398.1289412715064, 'dn': array([-11.52100055, -51.51418507])}

Residual enthalpy \(H^r(\vec{n}, V, T)\)

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.enthalpy_residual_vt(n, V, T) # Calculate residual enthalpy [bar L]

-1237.176213001435

You can obtain the next derivatives:

\[\left(\frac{\partial H^{r}_{(\vec{n}, V, T)}}{\partial T}\right)_{V, \vec{n}} \quad \left(\frac{\partial H^{r}_{(\vec{n}, V, T)}}{\partial V}\right)_{T, \vec{n}} \quad \left(\frac{\partial H^{r}_{(\vec{n}, V, T)}}{\partial n_j}\right)_{V, T, n_{i \neq j}}\]

Ask for derivatives as follows:

# Asking for all the derivatives
Hr, derivatives = model.enthalpy_residual_vt(n, V, T, dt=True, dv=True, dn=True)

# Displaying results
print("Residual Enthalpy: ", Hr)
print("Derivatives: ", derivatives)

Residual Enthalpy:  -1237.176213001435
Derivatives:  {'dt': 2.7262269081225208, 'dv': 222.65809494029838, 'dn': array([-130.30843347, -161.65842812])}

Residual internal energy \(U^r(\vec{n}, V, T)\)

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.internal_energy_residual_vt(n, V, T) # Calculate residual internal energy [bar L]

-1043.7424022970117

You can obtain the next derivatives:

\[\left(\frac{\partial U^{r}_{(\vec{n}, V, T)}}{\partial V}\right)_{T, \vec{n}} \quad \left(\frac{\partial U^{r}_{(\vec{n}, V, T)}}{\partial T}\right)_{V, \vec{n}} \quad \left(\frac{\partial U^{r}_{(\vec{n}, V, T)}}{\partial n_i}\right)_{V, T, n_{j \neq i}}\]

Ask for derivatives as follows:

# Asking for all the derivatives
Ur, derivatives = model.internal_energy_residual_vt(n, V, T, dt=True, dv=True, dn=True)

# Displaying results
print("Residual Internal Energy: ", Ur)
print("Derivatives: ", derivatives)

Residual Internal Energy:  -1043.7424022970117
Derivatives:  {'dt': 0.6569367874165049, 'dv': 814.2208469162279, 'dn': array([-184.00664642, -187.58600342])}

Residual entropy \(S^r(\vec{n}, V, T)\)

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.entropy_residual_vt(n, V, T) # Calculate residual entropy [bar L / K]

-1.746237812013294

You can obtain the next derivatives:

\[\left(\frac{\partial S^{r}_{(\vec{n}, V, T)}}{\partial T}\right)_{V, \vec{n}} \quad \left(\frac{\partial S^{r}_{(\vec{n}, V, T)}}{\partial V}\right)_{T, \vec{n}} \quad \left(\frac{\partial S^{r}_{(\vec{n}, V, T)}}{\partial n_j}\right)_{V, T, n_{i \neq j}}\]

Ask for derivatives as follows:

# Asking for all the derivatives
Sr, derivatives = model.entropy_residual_vt(n, V, T, dt=True, dv=True, dn=True)

# Displaying results
print("Residual Entropy: ", Sr)
print("Derivatives: ", derivatives)

Residual Entropy:  -1.746237812013294
Derivatives:  {'dt': 0.0021897892913883494, 'dv': 2.069290120706016, 'dn': array([-0.39595811, -0.36714748])}

Residual heat capacity at constant volume \(C_V^r(\vec{n}, V, T)\)

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.cv_residual_vt(n, V, T) # Calculate residual heat capacity [bar L / K]

0.6569367874165049

Residual heat capacity at constant pressure \(C_P^r(\vec{n}, V, T)\)

n = [5.0, 5.0]  # number of moles [mol]
V = 1.0         # volume [L]
T = 300.0       # temperature [K]

model.cp_residual_vt(n, V, T) # Calculate residual heat capacity [bar L / K]

3.72361653132665

Properties function of Pressure and Temperature \(\left(P, T \right)\)

We can obtain different bulk properties from a \(A^r\) model specifying the number of moles, pressure and temperature. But also, derivatives could be obtained if asked for.

All not asked derivatives will be set as None. The derivatives are returned as a Python dictionary with the keys dp, dt, and dn, etc. respectively. The compositional derivatives are always returned as Numpy arrays.

When we are calculating properties with pressure as specified variable, we have to ask for a specific volume root. The possible roots are:

• “liquid”

• “vapor”

• “stable” (default if not specified)

The “stable” root refers to the stable phase of the system. If this options is set, both liquid and vapor roots will be calculated and the phase with the lower residual Gibbs free energy \((G^{r}_{(\vec{n},P,T)})\) will be returned as the stable phase.

The following properties can be calculated:

Volume

n = [5.0, 5.0]  # number of moles [mol]
P = 1.0         # pressure [bar]
T = 300.0       # temperature [K]

volume_liquid = model.volume(n, P, T, root="liquid")
volume_vapor = model.volume(n, P, T, root="vapor")
volume_stable = model.volume(n, P, T, root="stable")

# Displaying results
print("Volume (liquid) in liters: ", volume_liquid)
print("Volume (vapor) in liters: ", volume_vapor)
print("Volume (stable phase) in liters: ", volume_stable)

Volume (liquid) in liters:  1.1546008756262833
Volume (vapor) in liters:  245.62683294458913
Volume (stable phase) in liters:  245.62683294458913

Logarithm of fugacity coefficients \(ln \, \phi(\vec{n},P,T)\)

n = [5.0, 5.0]  # number of moles [mol]
P = 1.0         # pressure [bar]
T = 300.0       # temperature [K]

ln_phi_l = model.lnphi_pt(n, P, T, root="liquid")
ln_phi_v = model.lnphi_pt(n, P, T, root="vapor")
ln_phi_s = model.lnphi_pt(n, P, T, root="stable")

# Displaying results
print("Logarithm of fugacity coefficients (liquid): ", ln_phi_l)
print("Logarithm of fugacity coefficients (vapor): ", ln_phi_v)
print("Logarithm of fugacity coefficients (stable phase): ", ln_phi_s)

Logarithm of fugacity coefficients (liquid):  [2.6074475  2.24385923]
Logarithm of fugacity coefficients (vapor):  [-0.01507129 -0.01531375]
Logarithm of fugacity coefficients (stable phase):  [-0.01507129 -0.01531375]

You can obtain the next derivatives (on the specified volume root):

\[\left(\frac{\partial ln \, \hat{\phi}_i}{\partial T}\right)_{P, \vec{n}} \quad \left(\frac{\partial ln \, \hat{\phi}_i}{\partial P}\right)_{T, \vec{n}} \quad \left(\frac{\partial ln \, \hat{\phi}_i}{\partial n_j}\right)_{P, T, n_{k \neq j}}\]

Ask for derivatives as follows:

# Asking for all the derivatives
ln_phi_s, derivatives = model.lnphi_pt(n, P, T, root="stable", dt=True, dp=True, dn=True)

# Displaying results
print("Logarithm of fugacity coefficients (stable phase): ", ln_phi_s)
print("Derivatives: ", derivatives)

Logarithm of fugacity coefficients (stable phase):  [-0.01507129 -0.01531375]
Derivatives:  {'dt': array([0.00013302, 0.00013235]), 'dp': array([-0.01513287, -0.01539262]), 'dn': array([[-8.47067747e-08,  8.47067748e-08],
       [ 8.47067748e-08, -8.47067747e-08]])}