15. ThermoSTEAM 101#

15.1. Pure component chemical models#

Thermosteam packages chemical and mixture thermodynamic models in a flexible framework that allows users to fully customize and extend the models, as well as create new models. Central to all thermodynamic algorithms is the Chemical object, which contain the same thermodynamic models and data as thermo’s Chemical objects, but present an API more suitable for BioSTEAM’s needs.

import thermosteam as tmo
# Initialize chemical with an identifier (e.g. by name, CAS, InChI...)
Water = tmo.Chemical('Water')
Chemical: Water (phase_ref='l')
[Names]  CAS: 7732-18-5
         InChI: H2O/h1H2
         InChI_key: XLYOFNOQVPJJNP-U...
         common_name: water
         iupac_name: ('oxidane',)
         pubchemid: 962
         smiles: O
         formula: H2O
[Groups] Dortmund: <1H2O>
         UNIFAC: <1H2O>
         PSRK: <1H2O>
         NIST: <Empty>
[Data]   MW: 18.015 g/mol
         Tm: 273.15 K
         Tb: 373.12 K
         Tt: 273.16 K
         Tc: 647.1 K
         Pt: 611.65 Pa
         Pc: 2.2064e+07 Pa
         Vc: 5.5948e-05 m^3/mol
         Hf: -2.8582e+05 J/mol
         S0: 70 J/K/mol
         LHV: -44011 J/mol
         HHV: -0 J/mol
         Hfus: 6010 J/mol
         Sfus: None
         omega: 0.3443
         dipole: 1.85 Debye
         similarity_variable: 0.16653
         iscyclic_aliphatic: 0
         combustion: {'H2O': 1.0}

Access pure component data:

# CAS number
# Molecular weight (g/mol)
# Boiling point (K)

Temperature (in Kelvin) and pressure (in Pascal) dependent properties can be computed:

# Vapor pressure (Pa)
# Surface tension (N/m)
# Liquid molar volume (m^3/mol)
Water.V(phase='l', T=298.15, P=101325)
# Vapor molar volume (m^3/mol)
Water.V(phase='g', T=298.15, P=101325)
# With user-specified units of measure:
Water.get_property('rho', 'kg/m3', phase='l', T=298.15, P=101325)

Temperature dependent properties are managed by objects:

VaporPressure(CASRN="7732-18-5", Tb=373.124295848, Tc=647.096, Pc=22064000.0, omega=0.3443, extrapolation="AntoineAB|DIPPR101_ABC", method="IAPWS")

Phase dependent properties have attributes with model handles for each phase:

PhaseTPHandle(phase, T, P=None) -> V [m^3/mol]
VolumeLiquid(CASRN="7732-18-5", MW=18.01528, Tb=373.124295848, Tc=647.096, Pc=22064000.0, Vc=5.59480372671e-05, Zc=0.22943845208106295, omega=0.3443, dipole=1.85, Psat=VaporPressure(CASRN="7732-18-5", Tb=373.124295848, Tc=647.096, Pc=22064000.0, omega=0.3443, extrapolation="AntoineAB|DIPPR101_ABC", method="IAPWS"), eos=[PR(Tc=647.096, Pc=22064000.0, omega=0.3443, T=298.15, P=101325.0)], extrapolation="constant", method="HEOS_FIT", method_P="COSTALD_COMPRESSED", tabular_extrapolation_permitted=True)

A new model can be added easily using add_method, for example:

def User_antoine_model(T):
     return 10.0**(10.116 -  1687.537 / (T - 42.98))
Water.Psat.add_method(f=User_antoine_model, Tmin=273.20, Tmax=473.20)

The add_method method is a high level interface that even lets you create a constant model:

Water.Cn('l', 350)

Choose what model to use through the method attribute:

Water.Cn.l.method = 'ZABRANSKY_SPLINE_C'
Water.Cn('l', 350)

15.1.1. Managing chemical sets#

Define multiple chemicals as a Chemicals object:

chemicals = tmo.Chemicals(['Water', 'Ethanol'])
Chemicals([Water, Ethanol])

The chemicals are attributes:

(chemicals.Water, chemicals.Ethanol)
(Chemical('Water'), Chemical('Ethanol'))

Chemicals are indexable:

Water = chemicals['Water']
chemicals['Ethanol', 'Water']
[Chemical('Ethanol'), Chemical('Water')]

Chemicals are also iterable:

for chemical in chemicals:

More chemicals can also be appended:

Propanol = tmo.Chemical('Propanol')
Chemicals([Water, Ethanol, Propanol])

The main benefit of using a Chemicals object, is that they can be compiled and used as part of a thermodynamic property package, as defined through a Thermo object:

# A Thermo object is built with an iterable of Chemicals or their IDs.
# Default mixture, thermodynamic equilibrium models are selected.
thermo = tmo.Thermo(chemicals)
    chemicals=CompiledChemicals([Water, Ethanol, Propanol]),

Creating a thermo property package, may be a little challenging if some chemicals cannot be found in the database, in which case they can be built from scratch. A complete example on how this can be done is available in another tutorial.

15.2. Material and energy balance#

A Stream object is the main interface for estimating thermodynamic properties, vapor-liquid equilibrium, and material and energy balances. First set the thermo property package and we can start creating streams:

# This also works: tmo.settings.set_thermo(['Water', 'Ethanol', 'Propanol'])
s1 = tmo.Stream('s1', Water=20, Ethanol=20, units='kg/hr')
Stream: s1
phase: 'l', T: 298.15 K, P: 101325 Pa
flow (kg/hr): Water    20
              Ethanol  20

Create another stream at a higher temperature:

s2 = tmo.Stream('s2', Water=10, units='kg/hr', T=350, P=101325)
Stream: s2
phase: 'l', T: 350 K, P: 101325 Pa
flow (kg/hr): Water  10

Mix both streams into a new one:

s_mix = tmo.Stream('s_mix')
s_mix.mix_from([s1, s2])
Stream: s_mix
phase: 'l', T: 310.54 K, P: 101325 Pa
flow (kg/hr): Water    30
              Ethanol  20

Check the energy balance through enthalpy:

s_mix.H - (s1.H + s2.H)

Note that the balance is not perfect as the solver stops within a small temperature tolerance. However, the approximation is less than 0.01% off:

error = s_mix.H - (s1.H + s2.H)
percent_error = 100 * error / (s1.H + s2.H)

Split the mixture to two streams by defining the component splits:

# First define an array of component splits
component_splits = s_mix.chemicals.array(['Water', 'Ethanol'], [0, 1])
s_mix.split_to(s1, s2, component_splits)
s1.T = s2.T = s_mix.T # Take care of energy balance
Stream: s1
phase: 'l', T: 310.54 K, P: 101325 Pa
flow (kg/hr): Ethanol  20
Stream: s2
phase: 'l', T: 310.54 K, P: 101325 Pa
flow (kg/hr): Water  30

Separate out stream from mixture:

s_mix.show(flow='kg/hr') # Only enthanol will remain
Stream: s_mix
phase: 'l', T: 310.54 K, P: 101325 Pa
flow (kg/hr): Ethanol  20

Note that the energy balance still holds:

error = s_mix.H - s1.H
percent_error = 100 * error / s2.H

15.3. Flow rates#

The most convenient way to get and set flow rates is through the get_flow and set_flow methods:

# Set and get flow of a single chemical
# in gallons per minute
s1.set_flow(1, 'gpm', 'Water')
s1.get_flow('gpm', 'Water')
# Set and get flows of many chemicals
# in kilograms per hour
s1.set_flow([10, 20], 'kg/hr', ('Ethanol', 'Water'))
s1.get_flow('kg/hr', ('Ethanol', 'Water'))
array([10., 20.])

It is also possible to index flow rate data using chemical IDs through the imol, imass, and ivol indexers:

ChemicalMolarFlowIndexer (kmol/hr):
(l) Water    1.11
    Ethanol  0.217
s1.imol['Ethanol', 'Water']
array([0.217, 1.11 ])

All flow rates are stored as a sparse array in the mol attribute. These arrays work just like numpy arrays, but are more scalable (saving memory and increasing speed) for sparse chemical data:

s1.mol # Molar flow rates [kmol/hr]
sparse([1.11 , 0.217, 0.   ])

Mass and volumetric flow rates are also available for convenience:

sparse([20., 10.,  0.])
sparse([0.02 , 0.013, 0.   ])

The data of these arrays are linked to the molar flows:

# Mass flows are always up to date with molar flows
s1.mol[0] = 1
# Changing mass flows changes molar flows
s1.mass[0] *= 2
# New arrays are not linked to molar flows
s1.mass + 2 # A new array is created
sparse([38.031, 12.   ,  2.   ])
# Array methods are also the same

15.4. Thermal condition#

Temperature and pressure can be get and set through the T and P attributes:

s1.T = 400.
s1.P = 2 * 101325.
Stream: s1
phase: 'l', T: 400 K, P: 202650 Pa
flow (kmol/hr): Water    2
                Ethanol  0.217

The phase may also be changed (‘s’ for solid, ‘l’ for liquid, and ‘g’ for gas):

s1.phase = 'g'

Notice that VLE is not enforced, but it is possible to perform. For now, just check that the dew point is lower than the actual temperature to assert it must be gas:

dp = s1.dew_point_at_P() # Dew point at constant pressure
DewPointValues(T=390.81, P=202650, IDs=('Water', 'Ethanol'), z=[0.902 0.098], x=[0.991 0.009])
dp.T < s1.T

It is also possible to get and set in other units of measure:

s1.set_property('P', 1, 'atm')
s1.get_property('P', 'atm')
s1.set_property('T', 125, 'degC')
s1.get_property('T', 'degF')

Enthalpy and entropy can also be set. In both cases, an energy balance is made to solve for temperature at isobaric conditions:

s1.H = 1.05 * s1.H
s1.get_property('T', 'degC') # Temperature should go up
s1.S = 0.95 * s1.S
s1.get_property('T', 'degC') # Temperature should go down

15.5. Thermal properties#

Thermodynamic properties are pressure, temperature and phase dependent. In the following examples, let’s just use water as it is easier to check properties:

s_water = tmo.Stream('s_water', Water=1, units='kg/hr')
s_water.rho # Density [kg/m^3]
s_water.T = 350
s_water.rho # Density changes

Get properties in different units:

s_water.get_property('sigma', 'N/m') # Surface tension
s_water.get_property('V', 'm3/kmol') # Molar volume

15.6. Flow properties#

Several flow properties are available, such as net material and energy flow rates:

# Net molar flow rate [kmol/hr]
# Net mass flow rate [kg/hr]
# Net volumetric flow rate [m3/hr]
# Enthalpy flow rate [kJ/hr]
# Entropy flow rate [kJ/hr]
# Capacity flow rate [J/K]

15.7. Phase equilibrium#

When streams are created, by default, they do not perform phase equilibrium. To perform phase equilibrium assuming 2 liquid phases and 1 gas phase is possible, pass vlle=True:

s = tmo.Stream('s', Water=1, Ethanol=1, T=400, P=101325,
               phase='l', # Guess phase
s.show() # Notice how stream is a gas (not a liquid)
Stream: s
phase: 'g', T: 400 K, P: 101325 Pa
flow (kmol/hr): Water    1
                Ethanol  1

Performing VLLE is quite slow and it is not recommended to perform VLLE calculations unless 3 phases can exist. When possible, perform phase equilibrium assuming 2 phases instead. Before moving into vapor-liquid and liquid-liquid equilibrium calculations, it may be useful to have a look at the phase envelopes to understand chemical interactions and ultimately how they partition between phases.

Plot the binary phase evelope of two chemicals in vapor-liquid equilibrium at constant pressure:

import matplotlib.pyplot as plt
eq = tmo.equilibrium # Thermosteam's equilibrium module
eq.plot_vle_binary_phase_envelope(['Ethanol', 'Water'], P=101325)

Plot the ternary phase diagram of three chemicals in liquid-liquid equilibrium at constant pressure:

# You'll need to "pip install python-ternary" to run this line
eq.plot_lle_ternary_diagram('Water', 'Ethanol', 'EthylAcetate', T=298.15)

15.8. Vapor-liquid equilibrium#

Vapor-liquid equilibrium can be performed by setting 2 degrees of freedom from the following list: T (Temperature; in K), P (Pressure; in Pa), V (Vapor fraction), H (Enthalpy; in kJ/hr), and S (Entropy; in kJ/K/hr).

For example, set vapor fraction and pressure:

s_eq = tmo.Stream('s_eq', Water=10, Ethanol=10)
s_eq.vle(V=0.5, P=101325)
MultiStream: s_eq
phases: ('g', 'l'), T: 353.94 K, P: 101325 Pa
composition (%): (g) Water    38.7
                     Ethanol  61.3
                     -------  10 kmol/hr
                 (l) Water    61.3
                     Ethanol  38.7
                     -------  10 kmol/hr

Note that the stream is now a MultiStream to manage multiple phases. Each phase can be accessed separately too:

phase: 'l', T: 353.94 K, P: 101325 Pa
flow (kmol/hr): Water    6.13
                Ethanol  3.87
phase: 'g', T: 353.94 K, P: 101325 Pa
flow (kmol/hr): Water    3.87
                Ethanol  6.13

Note that the phase of these substreams cannot be changed:

s_eq['g'].phase = 'l'
AttributeError                            Traceback (most recent call last)
Cell In[68], line 1
----> 1 s_eq['g'].phase = 'l'

File ~\code\biosteam\thermosteam\thermosteam\_stream.py:973, in Stream.phase(self, phase)
    971 @phase.setter
    972 def phase(self, phase):
--> 973     self._imol._phase.phase = phase

File ~\code\biosteam\thermosteam\thermosteam\_phase.py:189, in LockedPhase.__setattr__(self, name, value)
    187 def __setattr__(self, name, value):
    188     if value != self.phase:
--> 189         raise AttributeError('phase is locked')

AttributeError: phase is locked

Again, the most convenient way to get and set flow rates is through the get_flow and set_flow methods:

# Set flow of liquid water
s_eq.set_flow(1, 'gpm', ('l', 'Water'))
s_eq.get_flow('gpm', ('l', 'Water'))
# Set multiple liquid flows
key = ('l', ('Ethanol', 'Water'))
s_eq.set_flow([10, 20], 'kg/hr', key)
s_eq.get_flow('kg/hr', key)
array([10., 20.])

Chemical flows across all phases can be retrieved if no phase is given:

# Get water and ethanol flows summed across all phases
s_eq.get_flow('kg/hr', ('Water', 'Ethanol'))
array([ 89.791, 292.216])

However, setting chemical data of MultiStream objects requires the phase to be specified:

s_eq.set_flow([10, 20], 'kg/hr', ('Water', 'Ethanol'))
IndexError                                Traceback (most recent call last)
Cell In[72], line 1
----> 1 s_eq.set_flow([10, 20], 'kg/hr', ('Water', 'Ethanol'))

File ~\code\biosteam\thermosteam\thermosteam\_multi_stream.py:458, in MultiStream.set_flow(self, data, units, key)
    456 name, factor = self._get_flow_name_and_factor(units)
    457 indexer = getattr(self, 'i' + name)
--> 458 indexer[key] = np.asarray(data, dtype=float) / factor

File ~\code\biosteam\thermosteam\thermosteam\indexer.py:1036, in MaterialIndexer.__setitem__(self, key, data)
   1034 index, kind, sum_across_phases = self._get_index_data(key)
   1035 if sum_across_phases:
-> 1036     raise IndexError("multiple phases present; must include phase key "
   1037                      "to set chemical data")
   1038 if kind is None:
   1039     if index is None:

IndexError: multiple phases present; must include phase key to set chemical data

Similar to Stream objects, all flow rates can be accessed through the imol, imass, and ivol attributes:

s_eq.imol # Molar flow rates
MolarFlowIndexer (kmol/hr):
(g) Water     3.87
    Ethanol   6.13
(l) Water     1.11
    Ethanol   0.217
# Index a single chemical in the liquid phase
s_eq.imol['l', 'Water']
# Index multiple chemicals in the liquid phase
s_eq.imol['l', ('Ethanol', 'Water')]
array([0.217, 1.11 ])
# Index the vapor phase
sparse([3.874, 6.126, 0.   ])
# Index flow of chemicals summed across all phases
s_eq.imol['Ethanol', 'Water']
array([6.343, 4.984])

Because multiple phases are present, overall chemical flows in MultiStream objects cannot be set like in Stream objects:

s_eq.imol['Ethanol', 'Water'] = [1, 0]
IndexError                                Traceback (most recent call last)
Cell In[78], line 1
----> 1 s_eq.imol['Ethanol', 'Water'] = [1, 0]

File ~\code\biosteam\thermosteam\thermosteam\indexer.py:1036, in MaterialIndexer.__setitem__(self, key, data)
   1034 index, kind, sum_across_phases = self._get_index_data(key)
   1035 if sum_across_phases:
-> 1036     raise IndexError("multiple phases present; must include phase key "
   1037                      "to set chemical data")
   1038 if kind is None:
   1039     if index is None:

IndexError: multiple phases present; must include phase key to set chemical data

Chemical flows must be set by phase:

s_eq.imol['l', ('Ethanol', 'Water')] = [1, 0]

One main difference between a MultiStream object and a Stream object is that the mol attribute no longer stores any data, it simply returns the total flow rate of each chemical. Setting an element of the array raises an error to prevent the wrong assumption that the data is linked:

sparse([3.874, 7.126, 0.   ])
s_eq.mol[0] = 1
ValueError                                Traceback (most recent call last)
Cell In[81], line 1
----> 1 s_eq.mol[0] = 1

File ~\code\biosteam\thermosteam\thermosteam\base\sparse.py:1646, in SparseVector.__setitem__(self, index, value)
   1645 def __setitem__(self, index, value):
-> 1646     if self.read_only: raise ValueError('assignment destination is read-only')
   1647     dct = self.dct
   1648     if index.__class__ is tuple:

ValueError: assignment destination is read-only

Note that for both Stream and MultiStream objects, get_flow, imol, and mol return chemical flows across all phases when given only chemical IDs.

15.9. Liquid-liquid equilibrium#

Liquid-liquid equilibrium (LLE) only requires the temperature. Pressure is not a significant variable as liquid fungacity coefficients are not a strong function of pressure.

tmo.settings.set_thermo(['Water', 'Butanol', 'Octane'])
liquid_mixture = tmo.Stream('liquid_mixture', Water=100, Octane=100, Butanol=5)
MultiStream: liquid_mixture
phases: ('L', 'l'), T: 300 K, P: 101325 Pa
flow (kmol/hr): (L) Water    98.5
                    Butanol  1.21
                    Octane   0.00198
                (l) Water    1.46
                    Butanol  3.79
                    Octane   100