# Chemical equilibrium with fixed volume and internal energy¶

Written by Allan Leal (ETH Zurich) on Jan 7th, 2022 (last revised on Jan 21th, 2022)

Let’s consider the combustion of CH_{4} again. However, this time, the combustion process happens in a rigid and adiabatic chamber of 10 cm^{3}. Therefore, we will solve a chemical equilibrium problem in which we specify **volume** and **internal energy** as known properties at the equilibrium state (since these properties must be preserved within the chamber!).

In this problem, both temperature and pressure are unknown and will be calculated together with the amounts of species that bring the system to a state of chemical equilibrium.

Let’s create our `ChemicalSystem`

object (one with a single gaseous phase):

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

Next, we create an initial chemical state for this system in which CH_{4} and O_{2} exist in a 10 cm^{3} volume with mole fractions 0.75 and 0.25 respectively. You’ll notice that we accomplish this by setting the amounts of CH_{4} and O_{2} to 0.75 and 0.25 moles, respectively, followed by **scaling the volume of the chemical state** to our desired 10 cm^{3}.

```
state0 = ChemicalState(system)
state0.temperature(25.0, "celsius")
state0.pressure(1.0, "bar")
state0.set("CH4", 0.75, "mol")
state0.set("O2", 0.25, "mol")
state0.scaleVolume(10.0, "cm3")
print("INITIAL STATE")
print(state0)
```

```
INITIAL STATE
+-----------------+-------------+------+
| Property | Value | Unit |
+-----------------+-------------+------+
| Temperature | 298.15 | K |
| Pressure | 100000 | Pa |
| Charge: | 0 | mol |
| Element Amount: | | |
| :: H | 0.00121019 | mol |
| :: C | 0.000302547 | mol |
| :: O | 0.000201698 | mol |
| Species Amount: | | |
| :: CH4 | 0.000302547 | mol |
| :: O2 | 0.000100849 | mol |
| :: CO2 | 4.03395e-20 | mol |
| :: CO | 4.03395e-20 | mol |
| :: H2O | 4.03395e-20 | mol |
| :: H2 | 4.03395e-20 | mol |
+-----------------+-------------+------+
```

Let’s check the volume and internal energy in this initial chemical state. The next code block computes the chemical properties of the system in this initial state (using `ChemicalProps`

). It also stores the value of volume and internal energy in the variables `V0`

and `U0`

respectively, which we will use later.

```
props0 = ChemicalProps(state0)
V0 = props0.volume() # the initial volume of the gases
U0 = props0.internalEnergy() # the initial internal energy of the gases
print("Initial volume of the gases is", V0, "m3")
print("Initial internal energy of the gases is", U0, "J")
```

```
Initial volume of the gases is 1e-05 m3
Initial internal energy of the gases is -23.5698 J
```

Our chemical equilibrium problem (to model this combustion process) needs to impose the following properties at equilibrium:

volume; and

internal energy.

The code below creates an `EquilibriumSpecs`

object to reflect these specifications for the equilibrium constraints, which is then used to create our `EquilibriumSolver`

object.

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

The next step is to create an `EquilibriumConditions`

object in which we specify the values of volume and internal energy at equilibrium (which should be those values at the initial state, and that’s why we created variables `V0`

and `U0`

before!). We do this in the following code block, which also sets lower and upper bounds for temperature and pressure (to avoid negative or extremely large values for these properties during the equilibrium calculation!):

```
conditions = EquilibriumConditions(specs)
conditions.volume(V0)
conditions.internalEnergy(U0)
conditions.setLowerBoundTemperature(25.0, "celsius")
conditions.setUpperBoundTemperature(2000.0, "celsius")
conditions.setLowerBoundPressure(1.0, "bar")
conditions.setUpperBoundPressure(1000.0, "bar")
```

Tip

It’s good to set lower and upper bounds for temperature and pressure when solving chemical equilibrium problems in which these properties are unknown. In the course of the algorithm execution, these properties could become negative or extremely large, causing errors in the algorithm. At the moment, no default lower and upper bounds are set for temperature and pressure because sensible values do not exist for every possible application. For example, consider we had opted for 0.1 K and 0.1 Pa as default lower bound values for temperature and pressure. What would happen if you relied on a thermodynamic model (e.g., an equation of state) that only works for temperatures and pressures above 10 °C and 1 atm, respectively?!

We have everything we need now to perform the equilibrium calculation, which is done next:

```
state = ChemicalState(state0) # let's create a copy of the initial state
solver.solve(state, conditions)
print("FINAL STATE")
print(state)
```

```
FINAL STATE
+-----------------+-------------+------+
| Property | Value | Unit |
+-----------------+-------------+------+
| Temperature | 1080.39 | K |
| Pressure | 583621 | Pa |
| Charge: | 0 | mol |
| Element Amount: | | |
| :: H | 0.00121019 | mol |
| :: C | 0.000302547 | mol |
| :: O | 0.000201698 | mol |
| Species Amount: | | |
| :: CH4 | 0.000128967 | mol |
| :: O2 | 1e-16 | mol |
| :: CO2 | 9.68738e-06 | mol |
| :: CO | 0.000163893 | mol |
| :: H2O | 1.84304e-05 | mol |
| :: H2 | 0.000328729 | mol |
+-----------------+-------------+------+
```

Note that O_{2} was entirely consumed in the process. Its amount (`1e-16`

moles) is not zero because Reaktoro’s equilibrium solver has a default lower bound of `1e-16`

moles for the species amounts. Thus, the equilibrium calculation converged with O_{2} *attached to its lower bound*.

Tip

You can use `EquilibriumOptions`

to change this default lower bound for species amounts:

```
options = EquilibriumOptions()
options.epsilon = 1e-30
solver = EquilibriumSolver(specs)
solver.setOptions(options)
```

The previously printed state does not show volume and internal energy. To verify our equilibrium state stored in `state`

has volume and internal energy equal to `V0`

and `U0`

respectively, we do the following:

```
V = state.props().volume()
U = state.props().internalEnergy()
print("Volume at initial state:", V0, "m3")
print("Volume at final state:", V, "m3")
print("Internal energy at initial state:", U0, "J")
print("Internal energy at final state:", U, "J")
```

```
Volume at initial state: 1e-05 m3
Volume at final state: 1e-05 m3
Internal energy at initial state: -23.5698 J
Internal energy at final state: -23.5698 J
```

The verification above demonstrates we found an equilibrium state in which volume and internal energy were preserved.

Let’s now check what is the final temperature and pressure when burning CH_{4} in that rigid and adiabatic chamber with our given initial conditions:

```
T = units.convert(state.temperature(), "K", "degC") # convert from K to °C
P = units.convert(state.pressure(), "Pa", "bar") # convert from Pa to bar
print("Temperature at final state:", T, "°C")
print("Pressure at final state:", P, "bar")
```

```
Temperature at final state: 807.237 °C
Pressure at final state: 5.83621 bar
```

This concludes this tutorial, which demonstrated how Reaktoro can be used to perform chemical equilibrium calculations with given volume and internal energy.

Keep on reading to learn how to set your own constraints from scratch!