# Chemical equilibrium with constraints

## Contents

# Chemical equilibrium with constraints#

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

Attention

Always make sure you are using the latest version of Reaktoro. Otherwise, some new features documented on this website will not work on your machine and you may receive unintuitive errors. Follow these update instructions to get the latest version of Reaktoro!

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: {state.equilibrium().titrantAmount('CH4')} mol")
```

```
Amount of CH4 titrated in: 0.3998837695333316 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.