# Chemical equilibrium with constraints

## Contents

# 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.1500 | K |
| Pressure | 1.0000 | bar |
| Charge: | 0.0000e+00 | mol |
| Element Amount: | | |
| :: H | 4.0000e+00 | mol |
| :: C | 1.0000e+00 | mol |
| :: O | 1.0000e+00 | mol |
| Species Amount: | | |
| :: CH4 | 1.0000e+00 | mol |
| :: O2 | 5.0000e-01 | mol |
| :: CO2 | 1.0000e-16 | mol |
| :: CO | 1.0000e-16 | mol |
| :: H2O | 1.0000e-16 | mol |
| :: H2 | 1.0000e-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.1500 | K |
| Pressure | 16603.2449 | bar |
| Charge: | 0.0000e+00 | mol |
| Element Amount: | | |
| :: H | 4.0000e+00 | mol |
| :: C | 1.0000e+00 | mol |
| :: O | 1.0000e+00 | mol |
| Species Amount: | | |
| :: CH4 | 7.1576e-01 | mol |
| :: O2 | 1.0000e-16 | mol |
| :: CO2 | 2.2509e-01 | mol |
| :: CO | 5.9148e-02 | mol |
| :: H2O | 4.9067e-01 | mol |
| :: H2 | 7.7814e-02 | 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.1500 | K |
| Pressure | 10000.0000 | bar |
| Volume | 1.0000e-05 | m3 |
| Gibbs Energy | -488420.1304 | J |
| Enthalpy | -184792.5765 | J |
| Entropy | 238.4853 | J/K |
| Internal Energy | -194792.5765 | J |
| Helmholtz Energy | -498420.1304 | J |
| Charge | 0.0000e+00 | mol |
| Element Amount: | | |
| :: H | 1.5995e+00 | mol |
| :: C | 3.9988e-01 | mol |
| :: O | 1.0000e+00 | mol |
| Species Amount: | | |
| :: CH4 | 1.2748e-01 | mol |
| :: O2 | 1.0000e-16 | mol |
| :: CO2 | 2.3311e-01 | mol |
| :: CO | 3.9293e-02 | mol |
| :: H2O | 4.9449e-01 | mol |
| :: H2 | 5.0306e-02 | mol |
| Mole Fraction: | | |
| :: CH4 | 1.3495e-01 | mol/mol |
| :: O2 | 1.0586e-16 | mol/mol |
| :: CO2 | 2.4676e-01 | mol/mol |
| :: CO | 4.1594e-02 | mol/mol |
| :: H2O | 5.2345e-01 | mol/mol |
| :: H2 | 5.3252e-02 | mol/mol |
| Activity Coefficient: | | |
| :: CH4 | 1.0000 | - |
| :: O2 | 1.0000 | - |
| :: CO2 | 1.0000 | - |
| :: CO | 1.0000 | - |
| :: H2O | 1.0000 | - |
| :: H2 | 1.0000 | - |
| Activity: | | |
| :: CH4 | 1.3495e+03 | - |
| :: O2 | 1.0586e-12 | - |
| :: CO2 | 2.4676e+03 | - |
| :: CO | 4.1594e+02 | - |
| :: H2O | 5.2345e+03 | - |
| :: H2 | 5.3252e+02 | - |
| lg(Activity): | | |
| :: CH4 | 3.1302 | - |
| :: O2 | -11.9753 | - |
| :: CO2 | 3.3923 | - |
| :: CO | 2.6190 | - |
| :: H2O | 3.7189 | - |
| :: H2 | 2.7263 | - |
| ln(Activity): | | |
| :: CH4 | 7.2075 | - |
| :: O2 | -27.5741 | - |
| :: CO2 | 7.8110 | - |
| :: CO | 6.0305 | - |
| :: H2O | 8.5630 | - |
| :: H2 | 6.2776 | - |
| Chemical Potential: | | |
| :: CH4 | -2.7843e+05 | J/mol |
| :: O2 | -5.8051e+05 | J/mol |
| :: CO2 | -6.2217e+05 | J/mol |
| :: CO | -3.2477e+05 | J/mol |
| :: H2O | -4.2294e+05 | J/mol |
| :: H2 | -1.2553e+05 | J/mol |
| Standard Volume: | | |
| :: CH4 | 0.0000e+00 | m3/mol |
| :: O2 | 0.0000e+00 | m3/mol |
| :: CO2 | 0.0000e+00 | m3/mol |
| :: CO | 0.0000e+00 | m3/mol |
| :: H2O | 0.0000e+00 | m3/mol |
| :: H2 | 0.0000e+00 | m3/mol |
| Standard Gibbs Energy (formation): | | |
| :: CH4 | -354726.7026 | J/mol |
| :: O2 | -288624.1836 | J/mol |
| :: CO2 | -704857.1814 | J/mol |
| :: CO | -388605.1648 | J/mol |
| :: H2O | -513583.4128 | J/mol |
| :: H2 | -191986.0517 | J/mol |
| Standard Enthalpy (formation): | | |
| :: CH4 | -14294.0889 | J/mol |
| :: O2 | 32390.4240 | J/mol |
| :: CO2 | -344886.6903 | J/mol |
| :: CO | -79597.1200 | J/mol |
| :: H2O | -204066.9932 | J/mol |
| :: H2 | 29072.1299 | J/mol |
| Standard Entropy (formation): | | |
| :: CH4 | 267.3940 | J/(mol*K) |
| :: O2 | 252.1420 | J/(mol*K) |
| :: CO2 | 282.7400 | J/(mol*K) |
| :: CO | 242.7114 | J/(mol*K) |
| :: H2O | 243.1107 | J/(mol*K) |
| :: H2 | 173.6309 | J/(mol*K) |
| Standard Internal Energy (formation): | | |
| :: CH4 | -14294.0889 | J/mol |
| :: O2 | 32390.4240 | J/mol |
| :: CO2 | -344886.6903 | J/mol |
| :: CO | -79597.1200 | J/mol |
| :: H2O | -204066.9932 | J/mol |
| :: H2 | 29072.1299 | J/mol |
| Standard Helmholtz Energy (formation): | | |
| :: CH4 | -354726.7026 | J/mol |
| :: O2 | -288624.1836 | J/mol |
| :: CO2 | -704857.1814 | J/mol |
| :: CO | -388605.1648 | J/mol |
| :: H2O | -513583.4128 | J/mol |
| :: H2 | -191986.0517 | J/mol |
| Standard Heat Capacity (constant P): | | |
| :: CH4 | 84.1523 | J/(mol*K) |
| :: O2 | 35.9257 | J/(mol*K) |
| :: CO2 | 56.9320 | 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.9320 | 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(f"Amount of CH4 titrated in: {float(state.equilibrium().explicitTitrantAmounts())} 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.