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}(\textbf{n}, V, T) = -\int_{\infty}^{V} \left(P(\textbf{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 (Equations of State section)

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(\textbf{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, \textbf{n}} \quad \left(\frac{\partial P}{\partial T}\right)_{V, \textbf{n}} \quad \left(\frac{\partial P}{\partial n_i}\right)_{V, T}\]

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(\textbf{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, \textbf{n}} \quad \left(\frac{\partial ln \, \hat{\phi}_i}{\partial P}\right)_{T, \textbf{n}} \quad \left(\frac{\partial ln \, \hat{\phi}_i}{\partial n_j}\right)_{P, T}\]

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(\textbf{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}_{(\textbf{n}, V, T)}}{\partial V}\right)_{T, \textbf{n}} \quad \left(\frac{\partial A^{r}_{(\textbf{n}, V, T)}}{\partial T}\right)_{V, \textbf{n}} \quad \left(\frac{\partial A^{r}_{(\textbf{n}, V, T)}}{\partial n_i}\right)_{V, T} \quad\]

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

\[\left(\frac{\partial^2 A^{r}_{(\textbf{n}, V, T)}}{\partial n_i \partial n_j}\right)_{V, T} \quad \left(\frac{\partial^2 A^{r}_{(\textbf{n}, V, T)}}{\partial V \partial n_i}\right)_{T} \quad \left(\frac{\partial^2 A^{r}_{(\textbf{n}, V, T)}}{\partial T \partial n_i}\right)_{V} \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(\textbf{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}_{(\textbf{n}, V, T)}}{\partial T}\right)_{V, \textbf{n}} \quad \left(\frac{\partial G^{r}_{(\textbf{n}, V, T)}}{\partial V}\right)_{T, \textbf{n}} \quad \left(\frac{\partial G^{r}_{(\textbf{n}, V, T)}}{\partial n_j}\right)_{V, T}\]

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(\textbf{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}_{(\textbf{n}, V, T)}}{\partial T}\right)_{V, \textbf{n}} \quad \left(\frac{\partial H^{r}_{(\textbf{n}, V, T)}}{\partial V}\right)_{T, \textbf{n}} \quad \left(\frac{\partial H^{r}_{(\textbf{n}, V, T)}}{\partial n_j}\right)_{V, T}\]

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(\textbf{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}_{(\textbf{n}, V, T)}}{\partial V}\right)_{T, \textbf{n}} \quad \left(\frac{\partial U^{r}_{(\textbf{n}, V, T)}}{\partial T}\right)_{V, \textbf{n}} \quad \left(\frac{\partial U^{r}_{(\textbf{n}, V, T)}}{\partial n_i}\right)_{V, T}\]

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(\textbf{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}_{(\textbf{n}, V, T)}}{\partial T}\right)_{V, \textbf{n}} \quad \left(\frac{\partial S^{r}_{(\textbf{n}, V, T)}}{\partial V}\right)_{T, \textbf{n}} \quad \left(\frac{\partial S^{r}_{(\textbf{n}, V, T)}}{\partial n_j}\right)_{V, T}\]

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(\textbf{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(\textbf{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(\textbf{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, \textbf{n}} \quad \left(\frac{\partial ln \, \hat{\phi}_i}{\partial P}\right)_{T, \textbf{n}} \quad \left(\frac{\partial ln \, \hat{\phi}_i}{\partial n_j}\right)_{P, T}\]

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]])}

Residual Helmholtz free energy \(A^r(\textbf{n},P,T)\)

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

ar_l = model.helmholtz_residual_pt(n, P, T, root="liquid")
ar_v = model.helmholtz_residual_pt(n, P, T, root="vapor")
ar_s = model.helmholtz_residual_pt(n, P, T, root="stable")

# Displaying results
print("Residual Helmholtz PT (liquid): ", ar_l)
print("Residual Helmholtz PT (vapor): ", ar_v)
print("Residual Helmholtz PT (stable phase): ", ar_s)

Residual Helmholtz PT (liquid):  853.3194041432837
Residual Helmholtz PT (vapor):  0.017516543801525675
Residual Helmholtz PT (stable phase):  0.017516543801525675

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

\[\left(\frac{\partial A^r(\textbf{n}, P, T)}{\partial T}\right)_{P, \textbf{n}} \quad \left(\frac{\partial A^r(\textbf{n}, P, T)}{\partial P}\right)_{T, \textbf{n}} \quad \left(\frac{\partial A^r(\textbf{n}, P, T)}{\partial n_j}\right)_{P, T}\]

Ask for derivatives as follows

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

# Displaying results
print("Residual Helmholtz PT (stable phase): ", ar_s)
print("Derivatives: ", derivatives)

Residual Helmholtz PT (stable phase):  0.017516543801525675
Derivatives:  {'dt': -0.0003279930026914013, 'dp': 0.03560764705662445, 'dn': array([0.00153596, 0.00196735])}

Residual Gibbs free energy \(G^r(\textbf{n}, P, T)\)

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

gr_l = model.gibbs_residual_pt(n, P, T, root="liquid")
gr_v = model.gibbs_residual_pt(n, P, T, root="vapor")
gr_s = model.gibbs_residual_pt(n, P, T, root="stable")

# Displaying results
print("Residual Gibbs PT (liquid): ", gr_l)
print("Residual Gibbs PT (vapor): ", gr_v)
print("Residual Gibbs PT (stable phase): ", gr_s)

Residual Gibbs PT (liquid):  605.0401264789099
Residual Gibbs PT (vapor):  -3.78952905165867
Residual Gibbs PT (stable phase):  -3.78952905165867

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

\[\left(\frac{\partial G^r(\textbf{n}, P, T)}{\partial T}\right)_{P, \textbf{n}} \quad \left(\frac{\partial G^r(\textbf{n}, P, T)}{\partial P}\right)_{T, \textbf{n}} \quad \left(\frac{\partial G^r(\textbf{n}, P, T)}{\partial n_j}\right)_{P, T}\]

Ask for derivatives as follows

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

# Displaying results
print("Residual Gibbs PT (stable phase): ", gr_s)
print("Derivatives: ", derivatives)

Residual Gibbs PT (stable phase):  -3.78952905165867
Derivatives:  {'dt': 0.020464395699106415, 'dp': -3.807045595410901, 'dn': array([-0.3759291 , -0.38197671])}

Residual enthalpy \(H^r(\textbf{n}, P, T)\)

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

hr_l = model.enthalpy_residual_pt(n, P, T, root="liquid")
hr_v = model.enthalpy_residual_pt(n, P, T, root="vapor")
hr_s = model.enthalpy_residual_pt(n, P, T, root="stable")

# Displaying results
print("Residual Enthalpy PT (liquid): ", hr_l)
print("Residual Enthalpy PT (vapor): ", hr_v)
print("Residual Enthalpy PT (stable phase): ", hr_s)

Residual Enthalpy PT (liquid):  -1180.945594725174
Residual Enthalpy PT (vapor):  -9.928847761390575
Residual Enthalpy PT (stable phase):  -9.928847761390575

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

\[\left(\frac{\partial H^r(\textbf{n}, P, T)}{\partial T}\right)_{P, \textbf{n}} \quad \left(\frac{\partial H^r(\textbf{n}, P, T)}{\partial P}\right)_{T, \textbf{n}} \quad \left(\frac{\partial H^r(\textbf{n}, P, T)}{\partial n_j}\right)_{P, T}\]

Ask for derivatives as follows

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

# Displaying results
print("Residual Enthalpy PT (stable phase): ", hr_s)
print("Derivatives: ", derivatives)

Residual Enthalpy PT (stable phase):  -9.928847761390575
Derivatives:  {'dt': 0.04583315148296449, 'dp': -10.044762206001554, 'dn': array([-0.99541312, -0.99035643])}

Residual internal energy \(U^r(\textbf{n}, P, T)\)

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

ur_l = model.internal_energy_residual_pt(n, P, T, root="liquid")
ur_v = model.internal_energy_residual_pt(n, P, T, root="vapor")
ur_s = model.internal_energy_residual_pt(n, P, T, root="stable")

# Displaying results
print("Residual Internal Energy PT (liquid): ", ur_l)
print("Residual Internal Energy PT (vapor): ", ur_v)
print("Residual Internal Energy PT (stable phase): ", ur_s)

Residual Internal Energy PT (liquid):  -932.6663170608002
Residual Internal Energy PT (vapor):  -6.121802165930379
Residual Internal Energy PT (stable phase):  -6.121802165930379

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

\[\left(\frac{\partial U^r(\textbf{n}, P, T)}{\partial T}\right)_{P, \textbf{n}} \quad \left(\frac{\partial U^r(\textbf{n}, P, T)}{\partial P}\right)_{T, \textbf{n}} \quad \left(\frac{\partial U^r(\textbf{n}, P, T)}{\partial n_j}\right)_{P, T}\]

Ask for derivatives as follows

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

# Displaying results
print("Residual Internal Energy PT (stable phase): ", ur_s)
print("Derivatives: ", derivatives)

Residual Internal Energy PT (stable phase):  -6.121802165930379
Derivatives:  {'dt': 0.025040762781166678, 'dp': -6.202108963534028, 'dn': array([-0.61794807, -0.60641237])}

Residual entropy \(S^r(\textbf{n}, P, T)\)

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

sr_l = model.entropy_residual_pt(n, P, T, root="liquid")
sr_v = model.entropy_residual_pt(n, P, T, root="vapor")
sr_s = model.entropy_residual_pt(n, P, T, root="stable")

# Displaying results
print("Residual Entropy PT (liquid): ", sr_l)
print("Residual Entropy PT (vapor): ", sr_v)
print("Residual Entropy PT (stable phase): ", sr_s)

Residual Entropy PT (liquid):  -5.953285737346947
Residual Entropy PT (vapor):  -0.020464395699106352
Residual Entropy PT (stable phase):  -0.020464395699106352

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

\[\left(\frac{\partial S^r(\textbf{n}, P, T)}{\partial T}\right)_{P, \textbf{n}} \quad \left(\frac{\partial S^r(\textbf{n}, P, T)}{\partial P}\right)_{T, \textbf{n}} \quad \left(\frac{\partial S^r(\textbf{n}, P, T)}{\partial n_j}\right)_{P, T}\]

Ask for derivatives as follows

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

# Displaying results
print("Residual Entropy PT (stable phase): ", sr_s)
print("Derivatives: ", derivatives)

Residual Entropy PT (stable phase):  -0.020464395699106352
Derivatives:  {'dt': 0.00015277717160988143, 'dp': -0.02079238870196884, 'dn': array([-0.00206495, -0.00202793])}

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

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

cv_l = model.cv_residual_pt(n, P, T, root="liquid")
cv_v = model.cv_residual_pt(n, P, T, root="vapor")
cv_s = model.cv_residual_pt(n, P, T, root="stable")

# Displaying results
print("Residual Heat Capacity CV (liquid): ", cv_l)
print("Residual Heat Capacity CV (vapor): ", cv_v)
print("Residual Heat Capacity CV (stable phase): ", cv_s)

Residual Heat Capacity CV (liquid):  0.5870249332719473
Residual Heat Capacity CV (vapor):  0.0038530934828700093
Residual Heat Capacity CV (stable phase):  0.0038530934828700093

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

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

cp_l = model.cp_residual_pt(n, P, T, root="liquid")
cp_v = model.cp_residual_pt(n, P, T, root="vapor")
cp_s = model.cp_residual_pt(n, P, T, root="stable")

# Displaying results
print("Residual Heat Capacity CP (liquid): ", cp_l)
print("Residual Heat Capacity CP (vapor): ", cp_v)
print("Residual Heat Capacity CP (stable phase): ", cp_s)

Residual Heat Capacity CP (liquid):  7.758917618189358
Residual Heat Capacity CP (vapor):  0.04583315148296441
Residual Heat Capacity CP (stable phase):  0.04583315148296441

Excess properties (properties of mixing)

Excess volume \(V^E(\textbf{n}, P, T)\)

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

ve_l = model.volume_excess(n, P, T, root="liquid")
ve_v = model.volume_excess(n, P, T, root="vapor")
ve_s = model.volume_excess(n, P, T, root="stable")

# Displaying results
print("Excess Volume (liquid): ", ve_l)
print("Excess Volume (vapor): ", ve_v)
print("Excess Volume (stable phase): ", ve_s)

Excess Volume (liquid):  -244.47214785405885
Excess Volume (vapor):  8.42149039748108e-05
Excess Volume (stable phase):  8.42149039748108e-05

Logarithm of activity coefficients \(ln \, \gamma(\textbf{n},P,T)\)

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

lngamma_l = model.ln_gamma(n, P, T, root="liquid")
lngamma_v = model.ln_gamma(n, P, T, root="vapor")
lngamma_s = model.ln_gamma(n, P, T, root="stable")

# Displaying results
print("Logarithm of Activity Coefficients (liquid): ", lngamma_l)
print("Logarithm of Activity Coefficients (vapor): ", lngamma_v)
print("Logarithm of Activity Coefficients (stable phase): ", lngamma_s)

Logarithm of Activity Coefficients (liquid):  [2.62251921 2.2591734 ]
Logarithm of Activity Coefficients (vapor):  [4.23355361e-07 4.23716698e-07]
Logarithm of Activity Coefficients (stable phase):  [4.23355361e-07 4.23716698e-07]

Helmholtz free energy of excess \(A^E(\textbf{n}, P, T)\)

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

ae_l = model.helmoltz_excess(n, P, T, root="liquid")
ae_v = model.helmoltz_excess(n, P, T, root="vapor")
ae_s = model.helmoltz_excess(n, P, T, root="stable")

# Displaying results
print("Excess Helmholtz Free Energy (liquid): ", ae_l)
print("Excess Helmholtz Free Energy (vapor): ", ae_v)
print("Excess Helmholtz Free Energy (stable phase): ", ae_s)

Excess Helmholtz Free Energy (liquid):  853.3019090288403
Excess Helmholtz Free Energy (vapor):  2.1429330636774715e-05
Excess Helmholtz Free Energy (stable phase):  2.1429330636774715e-05

Gibbs free energy of excess \(G^E(\textbf{n}, P, T)\)

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

ge_l = model.gibbs_excess(n, P, T, root="liquid")
ge_v = model.gibbs_excess(n, P, T, root="vapor")
ge_s = model.gibbs_excess(n, P, T, root="stable")

# Displaying results
print("Excess Gibbs Free Energy (liquid): ", ge_l)
print("Excess Gibbs Free Energy (vapor): ", ge_v)
print("Excess Gibbs Free Energy (stable phase): ", ge_s)

Excess Gibbs Free Energy (liquid):  608.8297611747815
Excess Gibbs Free Energy (vapor):  0.00010564423461158551
Excess Gibbs Free Energy (stable phase):  0.00010564423461158551

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

\[\left(\frac{\partial G^E(\textbf{n}, P, T)}{\partial T}\right)_{P, \textbf{n}} \quad \left(\frac{\partial G^E(\textbf{n}, P, T)}{\partial P}\right)_{T, \textbf{n}} \quad \left(\frac{\partial G^E(\textbf{n}, P, T)}{\partial n_j}\right)_{P, T}\]

Ask for derivatives as follows

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

# Displaying results
print("Excess Gibbs Free Energy (stable phase): ", ge_s)
print("Derivatives: ", derivatives)

Excess Gibbs Free Energy (stable phase):  0.00010564423461158551
Derivatives:  {'dt': -3.4411866979980597e-07, 'dp': 8.42149039748108e-05, 'dn': array([1.05599170e-05, 1.05689299e-05])}

Excess enthalpy \(H^E(\textbf{n}, P, T)\)

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

he_l = model.enthalpy_excess(n, P, T, root="liquid")
he_v = model.enthalpy_excess(n, P, T, root="vapor")
he_s = model.enthalpy_excess(n, P, T, root="stable")

# Displaying results
print("Excess Enthalpy (liquid): ", he_l)
print("Excess Enthalpy (vapor): ", he_v)
print("Excess Enthalpy (stable phase): ", he_s)

Excess Enthalpy (liquid):  -1171.0165380839467
Excess Enthalpy (vapor):  0.00020887983555152728
Excess Enthalpy (stable phase):  0.00020887983555152728

Excess internal energy \(U^E(\textbf{n}, P, T)\)

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

ue_l = model.internal_energy_excess(n, P, T, root="liquid")
ue_v = model.internal_energy_excess(n, P, T, root="vapor")
ue_s = model.internal_energy_excess(n, P, T, root="stable")

# Displaying results
print("Excess Internal Energy (liquid): ", ue_l)
print("Excess Internal Energy (vapor): ", ue_v)
print("Excess Internal Energy (stable phase): ", ue_s)

Excess Internal Energy (liquid):  -926.5443902298878
Excess Internal Energy (vapor):  0.00012466493157671648
Excess Internal Energy (stable phase):  0.00012466493157671648

Excess entropy \(S^E(\textbf{n}, P, T)\)

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

se_l = model.entropy_excess(n, P, T, root="liquid")
se_v = model.entropy_excess(n, P, T, root="vapor")
se_s = model.entropy_excess(n, P, T, root="stable")

# Displaying results
print("Excess Entropy (liquid): ", se_l)
print("Excess Entropy (vapor): ", se_v)
print("Excess Entropy (stable phase): ", se_s)

Excess Entropy (liquid):  -5.9328209975290935
Excess Entropy (vapor):  3.4411866979980597e-07
Excess Entropy (stable phase):  3.4411866979980597e-07