Representation of thermodynamic states

Goal of this notebook

• Learn about the State object, how to construct, and how to use it.

from feos.si import *
from feos.pcsaft import *
from feos.eos import *

The State object

The State object is the most important object type in \(\text{FeO}_\text{s}\). It defines a thermodynamic state in the natural variables of the Helmholtz energy - the amount of substance of each component, \(\mathbf{N}\), the volume, \(V\), and the temperature, \(T\).

Once a State object is constructed, we can calculate thermodynamic properties. Internally, \(\text{FeO}_\text{s}\) transforms the state variables to generalized hyper dual numbers (see the separate tutorial on the topic of dual numbers) with which partial derivatives of the Helmholtz energy are computed.

There are several ways to construct State objects:

1. Given the natural variables, \(\mathbf{N}, V, T\).

2. Given a combination of other state variables, such as \(\mathbf{N}, p, T\) or \(\mathbf{N}, p, h\).

3. At critical conditions.

4. At phase equilibrium (this will generate multiple State objects, one for each phase).

Constructor methods need the EquationOfState or HelmholtzEnergyFunctional object as input, since, except for given \(\mathbf{N}, V, T\), the density and/or temperature has to be iteratively determined for which derivatives of the Helmholtz energy with respect to volume and temperature are utilized.

# Equation of state object.
parameters = PcSaftParameters.from_json(
    ['hexane'],
    '../parameters/pcsaft/gross2001.json'
)
pcsaft = EquationOfState.pcsaft(parameters)

The default constructor

The default constructor, State(...), takes a combination of input state variables. The first argument, however, is always the equation of state. For all cases, if we do not define the amount of substance, it is set to the inverse of Avogradro’s number, \(N_\text{AV}^{-1}\). The default constructor takes SINumber (and SIArray1 for e.g. the partial densities) as input. If you want to learn more about dimensioned quantities, please take a look at the respective tutorial.

state_nvt = State(
    pcsaft,
    temperature=300.15*KELVIN,
    density=7.5182*KILO*MOL/METER**3,
    total_moles=100.0*MOL
)
state_nvt

temperature	density
300.15000 K	7.51820 kmol/m³

Since internally, only \(\mathbf{N}, V\), and \(T\) are stored, all other properties have to be computed even if they were used to create the State.

For example, we can create a State for given temperature and pressure (using the default amount of substance for a pure component). The density is thus determined iteratively and stored as a field which can be accessed via state.density. If we are interested in the pressure of the state, we have to call the state.pressure() method which computes the pressure as partial derivative even though we used the pressure to create the State in the first place.

state_npt = State(
    pcsaft,
    temperature=300.15*KELVIN,
    pressure=1.0*BAR
)
print('density : ', state_npt.density)
print('pressure: ', state_npt.pressure())

density :  7.518194138679665 kmol/m³
pressure:  100.00000000009686 kPa

In the above case, specifying temperature and pressure may not yield the expected result. Consider thermodynamic conditions near phase equilibrium. The resulting density (for given tempreature and pressure) can be that of a meta-stable liquid or vapor phase depending on the initial density for the iteration.

To control the initial values for the density iteration, you can use the density_initialization keyword. Below, we create two State objects for the same temperature and pressure with different initial densities denoted by the vapor and liquid keywords for density_initialization. Alternatively, a starting density (as SINumber) can be provided.

s_vapor = State(
    pcsaft,
    temperature=335.0*KELVIN,
    pressure=1.0*BAR,
    density_initialization='vapor'
)
print('mass density: ', s_vapor.mass_density())

mass density:  3.2263994087922345 kg/m³

s_liquid = State(
    pcsaft,
    temperature=335.0*KELVIN,
    pressure=1.0*BAR,
    density_initialization='liquid'
)
print('mass density: ', s_liquid.mass_density())

mass density:  616.3096597655958 kg/m³

If no value for density_initialization is provided, both a low and high density is used to as starting point for the iteration and only the stable phase is returned.

s = State(
    pcsaft,
    temperature=335.0*KELVIN,
    pressure=1.0*BAR
)
print('mass density:', s.mass_density())

mass density: 616.3096597655958 kg/m³

You can run a stability analysis for each state using the is_stable() method.

print('Vapor stable? ', s_vapor.is_stable())
print('Liquid stable?', s_liquid.is_stable())

Vapor stable?  False
Liquid stable? True

---

Stored information

Once we create a State object, we can access its fields withouth further computations:

• density: molar density of the thermodynamic state for the given substance(s)

• molefracs: molar fractions for each substance

• moles: amount of substance for each substance

• partial_density: molar density for each substance

• temperature: temperature

• total_moles: total amount of substance

• volume: volume

For an equation of state that implements a molar weight (i.e. stores the molar weight in the parameter set), mass specific properties are also available as methods:

• mass(): mass for each substance

• mass_density(): total mass density

• massfracs(): mass fractions for each substance

• total_mass(): total mass

• total_molar_weight(): total molar weight

state = State(
    pcsaft,
    temperature=335.0*KELVIN,
    pressure=1.0*BAR,
    total_moles=200.0*MOL
)
print('total moles: ', state.total_moles)
print('total mass : ', state.total_mass())

total moles:  200  mol
total mass :  17.235400000000002 kg

---

Computing properties

Thermodynamic properties can be computed by invoking the appropriate method. For example, we can compute the total system pressure via the pressure() method.

pressure = state.pressure()
print('pressure: ', pressure)

pressure:  99.99999999988452 kPa

For a full list of possible thermodyanmic properties, please refer to the API documentation of the State object or take a look at the very bottom of this notebook.

Some properties accept an optional Contributions object which we can use to compute specific contributions to the property. The Contributions object allows for four options:

• Contributions.IdealGas: only the ideal gas contribution is considered (which is defined by the ideal gas model for the de Broglie wavelength)

• Contributions.ResidualNvt: only the residual contributions to the Helmholtz energy with respect to an ideal gas for given \(\mathbf{N}, V, T\) are considered.

• Contributions.ResidualNpt: only the residual contributions to the Helmholtz energy with respect to an ideal gas for given \(\mathbf{N}, p, T\) are considered.

• Contributions.Total: all contributions to the Helmholtz energy (and thus the property of interest) are considered, i.e. ideal gas plus residual. This is the default for most properties if no argument is provided. Please refer to the method documentation if you are not sure about the contributions used.

print('entropy (default)    :', state.molar_entropy())
print('entropy (total)      :', state.molar_entropy(Contributions.Total))
print('entropy (ideal gas)  :', state.molar_entropy(Contributions.IdealGas))
print('entropy (residual)   :', state.molar_entropy(Contributions.ResidualNvt))
print('entropy (residual p) :', state.molar_entropy(Contributions.ResidualNpt))

entropy (default)    : -63.82223112938985  J/mol/K
entropy (total)      : -63.82223112938985  J/mol/K
entropy (ideal gas)  : -20.979592032157587  J/mol/K
entropy (residual)   : -42.84263909723226  J/mol/K
entropy (residual p) : -86.86192380683525  J/mol/K

state.molar_entropy() - (state.molar_entropy(Contributions.IdealGas) + state.molar_entropy(Contributions.ResidualNvt))

$0\,\mathrm{\frac{ J}{molK}}$

---

State at critical conditions

\(\text{FeO}_\text{s}\) provides constructors for State objects at critical conditions as well in form of static class methods.

• State.critical_point(...): critial point of the system

• State.critical_point_pure(...): critical point for each substance in the system

• State.critical_point_binary_p(...): critical point for binary system, given pressure

• State.critical_point_binary_t(...): critical point for binary system, given temperature

Optional keywords are:

• moles: amount of substance for each component. For mixtures this is mandatory.

• initial_temperature : initial value for temperature. Can be used to speed up / increase convergence.

• max_iter: number of allowed iterations. Can be increased if convergence is an issue.

• tol: tolerance for the solution

• verbosity: a Verbosity object can be used to print information of the computation. Can be used if convergence is an issue.

critical_point = State.critical_point(pcsaft)
critical_point

temperature	density
519.33427 K	2.65414 kmol/m³

print('Critical conditions for hexane')
print('temperature :', critical_point.temperature)
print('pressure    :', critical_point.pressure())
print('density     :', critical_point.mass_density())

Critical conditions for hexane
temperature : 519.3342707319016 K
pressure    : 3.5427176263083453 MPa
density     : 228.72574604854387 kg/m³

---

PhaseEquilibrium: State objects at phase equilibrium

Another common use case for equations of state is the computation of phase equilibria. In \(\text{FeO}_\text{s}\), we can generate multiple State objects that are in equilibrium using a PhaseEquilibrium object. We will not discuss PhaseEquilibrium objects in detail in this notebook but merely consider it as another way to generate State objects. Please refer to the tutorial about phase diagrams if you want to learn more.

For pure substances, we can generate two states in equilibrium using the PhaseEquilibrium.pure(...) static method. The resulting object contains two (or more) State objects which we in this case can access via the liquid and vapor fields. As before, we can now compute properties using these objects.

vle = PhaseEquilibrium.pure(
    pcsaft,
    temperature_or_pressure=1.0*BAR
)
vle

	temperature	density
phase 1	341.53511 K	36.63788 mol/m³
phase 2	341.53511 K	7.07977 kmol/m³

vle.liquid

temperature	density
341.53511 K	7.07977 kmol/m³

vle.vapor

temperature	density
341.53511 K	36.63788 mol/m³

enthalpy_of_vaporization = vle.vapor.molar_enthalpy() - vle.liquid.molar_enthalpy()
print(f'enthalpy of vaporization (T = {vle.vapor.temperature}): {enthalpy_of_vaporization}')

enthalpy of vaporization (T = 341.53510965735256 K): 29.123040330216202 kJ/mol

---

Caching partial derivatives of the Helmholtz energy

A State object caches partial derivatives of the Helmholtz energy. If efficiency is a concern, you might want to consider the order in which you compute properties for a given state. If a method is called multiple times, only the first call will invoke a computation while additional calls will pull prior results from the cache.

---

Contributions to the Helmholtz energy

If you are interested in developing an equation of state, you might find the ..._contributions() methods of a State object useful. These methods return the contributions to a property which can be insightful (or a useful debugging tool).

state = State(
    pcsaft,
    temperature=300.0*KELVIN,
    pressure=1.0*BAR,
    total_moles=25.0*MOL
)

state.helmholtz_energy_contributions()

[('Ideal gas (QSPR)', 264.00707466668604 kJ),
 ('Hard Sphere', 549.608863830199 kJ),
 ('Hard Chain', -159.16321691198146 kJ),
 ('Dispersion', -750.1133884410586 kJ)]

state.pressure_contributions()

[('Ideal gas (QSPR)', 18.75674450579379 MPa),
 ('Hard Sphere', 304.57586929403857 MPa),
 ('Hard Chain', -63.0137117835539 MPa),
 ('Dispersion', -260.2189020162783 MPa)]

state.chemical_potential_contributions(0)

[('Ideal gas (QSPR)', 13.054621772113414 kJ/mol),
 ('Hard Sphere', 62.48794000517027 kJ/mol),
 ('Hard Chain', -14.746316825114631 kJ/mol),
 ('Dispersion', -64.60937326973789 kJ/mol)]

---

Dynamic properties via entropy scaling

If an equation of state implements correlation functions for entropy scaling, it can be used to compute dynamic properties via entropy scaling. For more information, see the respective tutorial for entropy scaling.

---

List of State methods and fields

Constructors

• critical_point

• critical_point_binary_p

• critical_point_binary_t

• critical_point_pure

• tp_flash

Fields

• density,

• molefracs,

• moles,

• partial_density,

• temperature,

• total_moles,

• volume

Stability analysis

• is_stable,

• stability_analysis

Thermodynamic properties

• c_p,

• c_v,

• chemical_potential,

• chemical_potential_contributions,

• compressibility,

• d2p_drho2,

• d2p_dv2,

• dc_v_dt,

• dln_phi_dnj,

• dln_phi_dp,

• dln_phi_dt,

• dmu_dni,

• dmu_dt,

• dp_dni,

• dp_drho,

• dp_dt,

• dp_dv,

• ds_dt,

• enthalpy,

• entropy,

• gibbs_energy,

• helmholtz_energy,

• helmholtz_energy_contributions,

• internal_energy,

• isentropic_compressibility,

• isothermal_compressibility,

• joule_thomson,

• ln_phi,

• ln_phi_pure_liquid,

• ln_symmetric_activity_coefficient,

• molar_enthalpy,

• molar_entropy,

• molar_gibbs_energy,

• molar_helmholtz_energy,

• molar_internal_energy,

• partial_molar_volume,

• partial_molar_enthalpy,

• partial_molar_entropy,

• pressure,

• pressure_contributions

• structure_factor,

• thermodynamic_factor

Mass related and mass specific properties

• specific_enthalpy,

• specific_entropy,

• specific_gibbs_energy,

• specific_helmholtz_energy,

• specific_internal_energy,

• speed_of_sound,

• thermodynamic_factor,

• total_mass,

• total_molar_weight

Dynamic properties (entropy scaling)

• diffusion,

• diffusion_reference,

• ln_diffusion_reduced,

• ln_thermal_conductivity_reduced,

• ln_viscosity_reduced,

• thermal_conductivity,

• thermal_conductivity_reference,

• viscosity,

• viscosity_reference

Concluding remkars

Hopefully you found this example helpful. If you have comments, critique or feedback, please let us know and consider opening an issue on github.