# Chemical equilibrium with constraints¶

Written by Allan Leal (ETH Zurich) on Jan 6th, 2022

As we mentioned in the previous section, Reaktoro can be used for a wide range of chemical equilibrium problems. In this section, we show:

the fundamentals of specifying different equilibrium constraints when calculating chemical equilibrium; and

how the equilibrium solver can be configured to allow the system to be open to one or more substances.

We learned that a chemical equilibrium solver can be created with the class `EquilibriumSolver`

, as shown in the example below:

```
from reaktoro import *
db = NasaDatabase("nasa-cea")
gases = GaseousPhase("CH4 O2 CO2 CO H2O H2")
system = ChemicalSystem(db, gases)
solver = EquilibriumSolver(system)
```

`EquilibriumSolver`

objects created this way (i.e., as in `solver = EquilibriumSolver(system)`

), are limited to chemical equilibrium calculations where:

temperature and pressure are given;

the system is closed (no mass transfer in or out of the system).

Thus, `solver = EquilibriumSolver(system)`

creates a **Gibbs energy minimization solver**.

This is not always the desired behavior, however. If we want to impose the temperature and volume of the system, we will need to construct a specialized `EquilibriumSolver`

object that can handle these constraint specifications.

Note

Chemical equilibrium problems in which temperature and volume are prescribed belong to the class of *Helmholtz energy minimization problems* [Smith and Missen, 1982]. Reaktoro, however, does not implement a Helmholtz energy minimization solver. Instead, Reaktoro implements a single **parametric Gibbs energy minimization solver that accepts general equilibrium constraints**, which can be configured to efficiently solve all other classes of equilibrium problems (including Helmholtz energy minimization problems!).

## Specifying constrained properties¶

Let’s create an equilibrium solver that can solve problems with prescribed **temperature** and **volume**:

```
specs = EquilibriumSpecs(system)
specs.temperature()
specs.volume()
solver = EquilibriumSolver(specs)
```

That’s it! We use `EquilibriumSpecs`

to provide the constraint specifications for our `EquilibriumSolver`

.

Tip

Access the link `EquilibriumSpecs`

to find out all currently existing constraints you can impose! It is even **possible to define your own constraints**, instead of predefined ones, which we’ll cover in a dedicated tutorial given this is a more advanced use case.

It’s now time to demonstrate how we use this specialized equilibrium solver. We want to model the combustion of 1 mol of CH_{4} and 0.5 mol of O_{2} in a chamber of 10 cm³ at 1000 °C. We create below an initial chemical state with this gas composition:

```
state = ChemicalState(system)
state.set("CH4", 1.0, "mol")
state.set("O2", 0.5, "mol")
print("INITIAL STATE")
print(state)
```

```
INITIAL STATE
+-----------------+--------+------+
| Property | Value | Unit |
+-----------------+--------+------+
| Temperature | 298.15 | K |
| Pressure | 100000 | Pa |
| Charge: | 0 | mol |
| Element Amount: | | |
| :: H | 4 | mol |
| :: C | 1 | mol |
| :: O | 1 | mol |
| Species Amount: | | |
| :: CH4 | 1 | mol |
| :: O2 | 0.5 | mol |
| :: CO2 | 1e-16 | mol |
| :: CO | 1e-16 | mol |
| :: H2O | 1e-16 | mol |
| :: H2 | 1e-16 | mol |
+-----------------+--------+------+
```

This state is in disequilibrium. We use `EquilibriumConditions`

below to specify the desired conditions of temperature and volume at the equilibrium state, and then equilibrate the `state`

object:

```
conditions = EquilibriumConditions(specs)
conditions.temperature(1000.0, "celsius")
conditions.volume(10.0, "cm3")
solver.solve(state, conditions)
print("FINAL STATE")
print(state)
```

```
FINAL STATE
+-----------------+-------------+------+
| Property | Value | Unit |
+-----------------+-------------+------+
| Temperature | 1273.15 | K |
| Pressure | 1.66032e+09 | Pa |
| Charge: | 0 | mol |
| Element Amount: | | |
| :: H | 4 | mol |
| :: C | 1 | mol |
| :: O | 1 | mol |
| Species Amount: | | |
| :: CH4 | 0.71576 | mol |
| :: O2 | 1e-16 | mol |
| :: CO2 | 0.225093 | mol |
| :: CO | 0.0591477 | mol |
| :: H2O | 0.490667 | mol |
| :: H2 | 0.0778142 | mol |
+-----------------+-------------+------+
```

This printed state does not show volume. Let’s check the volume in the `ChemicalProps`

object attached by `EquilibriumSolver`

to the `ChemicalState`

object `state`

at the end of the computation. Let’s also print the computed pressure:

```
V = state.props().volume() # volume in m³
P = state.pressure() # pressure in Pa
V = units.convert(V, "m3", "cm3") # convert volume from m³ to cm³
P = units.convert(P, "Pa", "GPa") # convert pressure from Pa to GPa
print("Volume at equilibrium state is", V, "cm³!")
print("Pressure at equilibrium state is", P, "GPa!")
```

```
Volume at equilibrium state is 10 cm³!
Pressure at equilibrium state is 1.66032 GPa!
```

Note that the volume at the computed equilibrium state is exactly what we enforced in the `EquilibriumConditions`

object: 10 cm³.

Attention

Don’t forget to set the expected conditions in the `EquilibriumConditions`

object. For every condition marked to be known in the `EquilibriumSpecs`

object (e.g., `specs.temperature()`

and `specs.volume()`

), make sure you give a value for each of them in the associated `EquilibriumConditions`

object (e.g., `conditions.temperature(1000.0, "celsius")`

and `conditions.volume(10.0, "cm3")`

)! Not doing this can cause unexpected behavior or a runtime error.

Did you know?

The Gibbs energy minimization computation performed with:

```
solver = EquilibriumSolver(system)
solver.solve(state)
```

is executed in the background with the following code:

```
specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
conditions = EquilibriumConditions(specs)
conditions.temperature(state.temperature())
conditions.pressure(state.pressure())
solver = EquilibriumSolver(specs)
solver.solve(state, conditions)
```

so you can type a lot less if all you need is equilibrium calculations with specified temperatures and pressures!

## Specifying open conditions¶

Let’s say we would like to find out how many mols of CH_{4} must react with 0.5 mol of O_{2} in the same chamber (with volume 10 cm³ and temperature 1000 °C) to obtain pressure 1 GPa (which, let’s assume, is the maximum pressure supported by the chamber material).

In this new equilibrium calculation, the system is open to the mass transfer of CH_{4}, and this creates a new *degree of freedom* in the problem: the amount added/removed of this substance. That’s why we can simultaneously impose a value for temperature (1000 °C), volume (10 cm³), and pressure (1 GPa), which would be thermodynamically difficult otherwise.

The code below starts by creating a chemical state in which there is only 0.5 mols of O_{2}:

```
state = ChemicalState(system)
state.set("O2", 0.5, "mol")
```

Next, we create an `EquilibriumSpecs`

object in which we specify that **temperature**, **volume**, and **pressure** are constrained properties. We also specify that the system is open to CH_{4}. This is all followed by the creation of a specialized and optimized `EquilibriumSolver`

object to handle this kind of equilibrium calculations:

```
specs = EquilibriumSpecs(system)
specs.temperature()
specs.volume()
specs.pressure()
specs.openTo("CH4")
solver = EquilibriumSolver(specs)
```

We can now create the `EquilibriumConditions`

object in which we provide the values for temperature, volume, and pressure, and then perform the calculation:

```
conditions = EquilibriumConditions(specs)
conditions.temperature(1000.0, "celsius")
conditions.volume(10.0, "cm3")
conditions.pressure(1.0, "GPa")
solver.solve(state, conditions)
print(state.props())
```

```
+----------------------------------------+-------------+-----------+
| Property | Value | Unit |
+----------------------------------------+-------------+-----------+
| Temperature | 1273.15 | K |
| Pressure | 1e+09 | Pa |
| Volume | 1e-05 | m3 |
| Gibbs Energy | -488420 | J |
| Enthalpy | -184793 | J |
| Entropy | 238.485 | J/K |
| Internal Energy | -194793 | J |
| Helmholtz Energy | -498420 | J |
| Charge | 0 | mol |
| Element Amount: | | |
| :: H | 1.59954 | mol |
| :: C | 0.399884 | mol |
| :: O | 1 | mol |
| Species Amount: | | |
| :: CH4 | 0.127484 | mol |
| :: O2 | 1e-16 | mol |
| :: CO2 | 0.233107 | mol |
| :: CO | 0.0392932 | mol |
| :: H2O | 0.494494 | mol |
| :: H2 | 0.050306 | mol |
| Mole Fraction: | | |
| :: CH4 | 0.134949 | mol/mol |
| :: O2 | 1.05856e-16 | mol/mol |
| :: CO2 | 0.246756 | mol/mol |
| :: CO | 0.041594 | mol/mol |
| :: H2O | 0.523449 | mol/mol |
| :: H2 | 0.0532517 | mol/mol |
| Activity Coefficient: | | |
| :: CH4 | 1 | - |
| :: O2 | 1 | - |
| :: CO2 | 1 | - |
| :: CO | 1 | - |
| :: H2O | 1 | - |
| :: H2 | 1 | - |
| Activity: | | |
| :: CH4 | 1349.49 | - |
| :: O2 | 1.05856e-12 | - |
| :: CO2 | 2467.56 | - |
| :: CO | 415.94 | - |
| :: H2O | 5234.49 | - |
| :: H2 | 532.517 | - |
| lg(Activity): | | |
| :: CH4 | 3.13017 | - |
| :: O2 | -11.9753 | - |
| :: CO2 | 3.39227 | - |
| :: CO | 2.61903 | - |
| :: H2O | 3.71887 | - |
| :: H2 | 2.72633 | - |
| ln(Activity): | | |
| :: CH4 | 7.20748 | - |
| :: O2 | -27.5741 | - |
| :: CO2 | 7.81099 | - |
| :: CO | 6.03054 | - |
| :: H2O | 8.56302 | - |
| :: H2 | 6.27761 | - |
| Chemical Potential: | | |
| :: CH4 | -278431 | J/mol |
| :: O2 | -580512 | J/mol |
| :: CO2 | -622174 | J/mol |
| :: CO | -324769 | J/mol |
| :: H2O | -422939 | J/mol |
| :: H2 | -125534 | J/mol |
| Standard Volume: | | |
| :: CH4 | 0 | m3/mol |
| :: O2 | 0 | m3/mol |
| :: CO2 | 0 | m3/mol |
| :: CO | 0 | m3/mol |
| :: H2O | 0 | m3/mol |
| :: H2 | 0 | m3/mol |
| Standard Gibbs Energy (formation): | | |
| :: CH4 | -354727 | J/mol |
| :: O2 | -288624 | J/mol |
| :: CO2 | -704857 | J/mol |
| :: CO | -388605 | J/mol |
| :: H2O | -513583 | J/mol |
| :: H2 | -191986 | J/mol |
| Standard Enthalpy (formation): | | |
| :: CH4 | -14294.1 | J/mol |
| :: O2 | 32390.4 | J/mol |
| :: CO2 | -344887 | J/mol |
| :: CO | -79597.1 | J/mol |
| :: H2O | -204067 | J/mol |
| :: H2 | 29072.1 | J/mol |
| Standard Entropy (formation): | | |
| :: CH4 | 267.394 | J/(mol*K) |
| :: O2 | 252.142 | J/(mol*K) |
| :: CO2 | 282.74 | J/(mol*K) |
| :: CO | 242.711 | J/(mol*K) |
| :: H2O | 243.111 | J/(mol*K) |
| :: H2 | 173.631 | J/(mol*K) |
| Standard Internal Energy (formation): | | |
| :: CH4 | -14294.1 | J/mol |
| :: O2 | 32390.4 | J/mol |
| :: CO2 | -344887 | J/mol |
| :: CO | -79597.1 | J/mol |
| :: H2O | -204067 | J/mol |
| :: H2 | 29072.1 | J/mol |
| Standard Helmholtz Energy (formation): | | |
| :: CH4 | -354727 | J/mol |
| :: O2 | -288624 | J/mol |
| :: CO2 | -704857 | J/mol |
| :: CO | -388605 | J/mol |
| :: H2O | -513583 | J/mol |
| :: H2 | -191986 | J/mol |
| Standard Heat Capacity (constant P): | | |
| :: CH4 | 84.1523 | J/(mol*K) |
| :: O2 | 35.9257 | J/(mol*K) |
| :: CO2 | 56.932 | J/(mol*K) |
| :: CO | 34.4671 | J/(mol*K) |
| :: H2O | 44.7425 | J/(mol*K) |
| :: H2 | 31.3023 | J/(mol*K) |
| Standard Heat Capacity (constant V): | | |
| :: CH4 | 84.1523 | J/(mol*K) |
| :: O2 | 35.9257 | J/(mol*K) |
| :: CO2 | 56.932 | J/(mol*K) |
| :: CO | 34.4671 | J/(mol*K) |
| :: H2O | 44.7425 | J/(mol*K) |
| :: H2 | 31.3023 | J/(mol*K) |
+----------------------------------------+-------------+-----------+
```

The above table shows that we obtained an equilibrium state with desired values for temperature, pressure, and volume (in SI units). We can also see below the amount of CH_{4} that entered the system so that all these prescribed conditions could be attained:

```
print("Amount of CH4 titrated in:", state.equilibrium().explicitTitrantAmounts()[0], "mol")
```

```
Amount of CH4 titrated in: 0.3998837695340861 mol
```

## Recap¶

Reaktoro can solve chemical equilibrium problems with a variety of specifications. This guide demonstrated how this can be done for some relatively simple cases using classes:

We will later learn how these classes can be used to model more complicated equilibrium problems.