Equilibrium via equilibrium constant - Computational Thermodynamics

Consider that you have a mixture with 1 kilomole of carbon monoxide (CO) that reacts with 0.5 kmol of oxygen (O $_2$ ) to form a mixture of CO, CO $_2$ , and O $_2$ , with the equilibrium conditions of 2500 K and (a) 1 atm (b) 10 atm.

Problem: Find the equilibrium composition in terms of the mole fraction.

Assume the mixture behaves as an ideal gas.

We will compare three solution methods based on the law of mass action and the equilibrium constant:

Using a tabulated equilibrium constant
Calculating equilibrium constant
Using a reaction coordinate

First, import the necessary modules, then specify the known information.

import numpy as np
import cantera as ct
from scipy.optimize import root, root_scalar

from pint import UnitRegistry
ureg = UnitRegistry()
Q_ = ureg.Quantity

# for convenience:
def to_si(quant):
    '''Converts a Pint Quantity to magnitude at base SI units.
    '''
    return quant.to_base_units().magnitude

temperature = Q_(2500, 'K')
pressures = Q_([1, 10], 'atm')

components = ['CO', 'O2', 'CO2']
moles_initial = np.array([1.0, 0.5, 0.0])
stoich_coefficients = np.array([1.0, 0.5, -1.0])

Using a tabulated equilibrium constant¶

With the pressure and temperature known, we can find the composition of the mixture by using a reaction equilibrium constant that is tabulated with respect to temperature. The primary reaction involved is

\text{CO}_2 \leftrightarrow \text{CO} + \frac{1}{2} \text{O}_2

(1)

where the equilibrium constant is

K(T) = \frac{y_{\text{CO}}^{\nu_{\text{CO}}} y_{\text{O}_2}^{\nu_{\text{O}_2}}}{y_{\text{CO}_2}^{\nu_{\text{CO}_2}}} \left(\frac{P}{P_{\text{ref}}} \right)^{ \nu_{\text{CO}} + \nu_{\text{O}_2} - \nu_{\text{CO}_2} } = \frac{y_{\text{CO}} y_{\text{O}_2}^{1/2}}{y_{\text{CO}_2}} \left(\frac{P}{P_{\text{ref}}} \right)^{1/2}

(2)

We can apply conservation of mass to find the overall balanced chemical reaction:

1\text{CO} + \frac{1}{2} \text{O}_2 \text{CO}_2 \rightarrow z \text{CO} + \frac{z}{2} \text{O}_2 + (1-z) \text{CO}_2

(3)

where $z$ is the amount of CO in kmol at equilibrium ( $0 \leq z \leq 1$ ). Then, the total number of moles $n$ in the mixture at equilibrium is:

n = z + \frac{z}{2} + (1-z) = \frac{2+z}{2}

(4)

so the mole fractions of each component at equilibrium are:

y_{\text{CO}} = \frac{2z}{2 + z} \\ y_{\text{O}_2} = \frac{z}{2+z} \\ y_{\text{CO}_2} = \frac{2(1-z)}{2+z}

(5)

Therefore, we can express the equilibrium constant as

K(T) = \frac{z}{1-z} \left(\frac{z}{2+z}\right)^{1/2} \left(\frac{P}{P_{\text{ref}}}\right)^{1/2}

(6)

At 2500 K, we can look up the tabulated value for equilibrium constant, which is

\log_{10} K (T = 2500 \, \text{K}) = -1.440

(7)

def solve_equilibrium_constant(z, pressure, equil_constant):
    pressure_ref = Q_(1, 'atm')
    K = (
        (z / (1.0 - z)) * np.sqrt(z / (2.0 + z)) *
        np.sqrt(to_si(pressure / pressure_ref))
        )
    
    return (equil_constant - K)

def get_mole_fractions(z):
    mole_frac_CO = 2 * z / (2 + z)
    mole_frac_O2 = z / (2 + z)
    mole_frac_CO2 = 2 * (1 - z) / (2 + z)
    return {'CO': mole_frac_CO, 'O2': mole_frac_O2, 'CO2': mole_frac_CO2}

# tabulated value of equilibrium constant at 2500 K
log10K = -1.440
equilibrium_constant = 10.0**log10K

print(f'Tabulated equilibrium constant: {equilibrium_constant: .4f}')

Tabulated equilibrium constant:  0.0363

# First pressure, 1 atm
pressure = pressures[0]

sol = root_scalar(
    solve_equilibrium_constant, x0=0.4, x1=0.5,
    args=(pressure, equilibrium_constant)
    )
mole_fractions = get_mole_fractions(sol.root)

print(f'Mole fractions at {pressure: .1f}')
for comp in components:
    print(f'{comp:3}: {mole_fractions[comp]: .3f}')

Mole fractions at  1.0 standard_atmosphere
CO :  0.121
O2 :  0.060
CO2:  0.819

# now evaluate composition at 10 atm
pressure = pressures[1]

sol = root_scalar(
    solve_equilibrium_constant, x0=0.4, x1=0.5,
    args=(pressure, equilibrium_constant)
    )
mole_fractions = get_mole_fractions(sol.root)

print(f'Mole fractions at {pressure: .1f}')
for comp in components:
    print(f'{comp:3}: {mole_fractions[comp]: .3f}')

Mole fractions at  10.0 standard_atmosphere
CO :  0.060
O2 :  0.030
CO2:  0.910

At 1 atm, the equilibrium composition has just 82% CO $_2$ by mole, while at 10 atm the mixture is 91% CO $_2$ at equilibrium.

Calculating equilibrium constant¶

Determining the equilibrium composition using this method is limited by having the tabulated equilibrium constant. However, this can be calculated by using the law of mass action and chemical species property information:

\Delta G^{\circ} = -R_{\text{univ}} T \ln K

(8)

where $\Delta G^{\circ}$ is the standard-state Gibbs free energy change of reaction and $R_{\text{univ}}$ is the universal gas constant. We can calculate $\Delta G^{\circ}$ for the above reaction:

\Delta G^{\circ} (T) = \nu_{\text{CO}} \overline{g^{\circ}}_{\text{CO}} + \nu_{\text{O}_2} \overline{g^{\circ}}_{\text{O}_2} - \nu_{\text{CO}_2} \overline{g^{\circ}}_{\text{CO}_2}\\ = \overline{g^{\circ}}_{\text{CO}} + \frac{1}{2} \overline{g^{\circ}}_{\text{O}_2} - \overline{g^{\circ}}_{\text{CO}_2}

(9)

where $\overline{g^{\circ}}_{i}$ is the molar-specific Gibbs free energy of substance $i$ at temperature $T$ and the reference pressure (1 atm).

To get the Gibbs free energy of each substance, we can use a Cantera Solution object with the state specified by the given temperature, pressure, and composition. We’ll evaluate each component separately, evaluating the property of each as a pure substance (i.e., with no other components present).

pressure = pressures[0]

# Load Cantera model for species information
gas = ct.Solution('gri30.yaml')

gas.TPX = to_si(temperature), to_si(pressure), 'CO2:1.0'
gibbs_CO2 = Q_(gas.gibbs_mole, 'J/kmol')

gas.TPX = to_si(temperature), to_si(pressure), 'CO:1.0'
gibbs_CO = Q_(gas.gibbs_mole, 'J/kmol')

gas.TPX = to_si(temperature), to_si(pressure), 'O2:1.0'
gibbs_O2 = Q_(gas.gibbs_mole, 'J/kmol')

gibbs_change_reaction = gibbs_CO + 0.5*gibbs_O2 - gibbs_CO2

equilibrium_constant = np.exp(
    -gibbs_change_reaction / 
    (Q_(ct.gas_constant, 'J/(kmol*K)') * temperature)
    )
print(f'Calculated equilibrium constant: {to_si(equilibrium_constant): .4f}')

Calculated equilibrium constant:  0.0368

This is very close to the value shown above obtained from tabulated data. Let’s now use this value to determine the equilibrium composition:

for pressure in pressures:
    sol = root_scalar(
        solve_equilibrium_constant, x0=0.4, x1=0.5,
        args=(pressure, to_si(equilibrium_constant))
        )
    mole_fractions = get_mole_fractions(sol.root)

    print(f'Mole fractions at {pressure: .1f}:')
    for comp in components:
        print(f'{comp}: {mole_fractions[comp]: .3f}')

    print()

Mole fractions at  1.0 standard_atmosphere:
CO:  0.122
O2:  0.061
CO2:  0.817

Mole fractions at  10.0 standard_atmosphere:
CO:  0.061
O2:  0.030
CO2:  0.909

Using a reaction coordinate¶

The methods we have used so far required reducing the three unknowns (numbers of moles) into a single unknown variable, based on conservation of mass applied to a single equation. This won’t work if multiple reactions are occuring, and we need a more-general approach.

We can use the concept of the reaction coordinate (or degree/extend of reaction), which is a proportionality constant that connects how the amount of each component changes as the reaction proceeds towards equilibrium. For the reaction and species we are considering, we can write

\frac{dn_{\text{CO}}}{\nu_{\text{CO}}} = \frac{dn_{\text{O}_2}}{\nu_{\text{O}_2}} = \frac{dn_{\text{CO}_2}}{\nu_{\text{CO}_2}} = d \epsilon \;,

(10)

where $\epsilon$ is the reaction coordinate. If we integrate that equation for each substance, where $\epsilon = 0$ at the start of reaction, we can obtain

n_i = n_{0,i} + \nu_i \epsilon \quad \text{for } i = 1, \ldots, C

(11)

where $n_i$ is the number of moles of component $i$ , $n_{0,i}$ is the initial number of moles of component $i$ , and $C$ is the number of components in the system.

By introducing one new unknown ( $\epsilon$ ), we get one additional equation for each component, which we can add to the law of mass action to set up a system of equations to solve for the unknowns: the number of moles of each component and the reaction coordinate.

Multiple reactions: This concept can be extended to multiple reactions, where we would need one reaction coordinate for each reaction:

n_i = n_{0,i} + \sum_{j=1}^R \nu_{i,j} \epsilon_j \quad \text{for } i = 1, \ldots, C

(12)

where $R$ is the number of reactions. The law of mass action then applies to each reaction:

\Delta G_j^{\circ} = -R_{\text{univ}} T \ln (K_j) \quad \text{for } j = 1, \ldots, R

(13)

def find_equilibrium_root(
    x, temperature, pressure, components, 
    moles_initial, stoich_coefficients, gas
    ):
    '''System of equations for reaction coordinate and equilibrium composition.
    '''
    epsilon = x[0]
    moles = np.array(x[1:])
    
    total_moles = np.sum(moles)
    mole_fractions = moles / total_moles
    
    # get standard-state Gibbs free energy of each component
    gibbs = np.zeros(len(components))
    for idx, comp in enumerate(components):
        gas.TPX = (
            to_si(temperature), to_si(Q_(1, 'atm')),
            f'{comp}:1.0'
            )
        gibbs[idx] = gas.gibbs_mole
        
    gibbs *= Q_('J/kmol')
        
    equil_constant = (
        np.prod([y**nu for y, nu in 
                 zip(mole_fractions, stoich_coefficients)
                ]) * 
        (pressure / Q_(1, 'atm'))**(np.sum(stoich_coefficients))
        )
    
    return [
        to_si(-np.sum(stoich_coefficients * gibbs) / 
            (Q_(ct.gas_constant, 'J/(kmol*K)') * temperature) - 
            np.log(equil_constant)
            ),
        moles[0] - (moles_initial[0] + stoich_coefficients[0] * epsilon),
        moles[1] - (moles_initial[1] + stoich_coefficients[1] * epsilon),
        moles[2] - (moles_initial[2] + stoich_coefficients[2] * epsilon),
        ]

pressure = pressures[0]
x0 = [-0.5, 0.5, 0.5, 0.5]
gas = ct.Solution('gri30.yaml')

sol = root(
    find_equilibrium_root, x0, method='lm',
    args=(temperature, pressure, components, moles_initial, stoich_coefficients, gas)
    )

print(f'Root-finding success: {sol.success}\n')

epsilon = sol.x[0]
moles = sol.x[1:]
mole_fractions = moles / np.sum(moles)


# Check constraints:
for idx, mole in enumerate(moles):
    if mole < 0:
        print(f'Error: moles of {components[idx]} below zero.')
        break
else:
    print(f'Mole fractions at {pressure: .1f}:')
    for idx, comp in enumerate(components):
        print(f'{comp:3}: {mole_fractions[idx]: .3f}')

Root-finding success: True

Mole fractions at  1.0 standard_atmosphere:
CO :  0.122
O2 :  0.061
CO2:  0.817

pressure = pressures[1]

x0 = [-0.5, 0.1, 0.1, 0.9]
gas = ct.Solution('gri30.yaml')

sol = root(
    find_equilibrium_root, x0, method='lm',
    args=(temperature, pressure, components, moles_initial, stoich_coefficients, gas)
    )

print(f'Root-finding success: {sol.success}\n')

epsilon = sol.x[0]
moles = sol.x[1:]
mole_fractions = moles / np.sum(moles)

# Check constraints:
for idx, mole in enumerate(moles):
    if mole < 0:
        print(f'Error: moles of {components[idx]} below zero.')
        break
else:
    print(f'Mole fractions at {pressure: .1f}:')
    for idx, comp in enumerate(components):
        print(f'{comp:3}: {mole_fractions[idx]: .3f}')

Root-finding success: True

Mole fractions at  10.0 standard_atmosphere:
CO :  0.061
O2 :  0.030
CO2:  0.909

Compare to Cantera equilibrium solution¶

Cantera has a handy built-in equilibrium solver, which we can compare our calculated equilibrium compositions against.

When using the equilibrate() function, you have to specify which two properties to hold constant. The options are ['TP', 'TV', 'HP', 'SP', 'SV', 'UV']; our current problem is examining an isothermal, isobaric process (TP).

We can use the gri30.yaml model, but by default it contains many more species than we want to consider (53 in total). So, we need to extract the three species we want (CO, O $_2$ , and CO $_2$ ), and construct a new Solution object with just those species.

temperature = Q_(2500, 'K')
pressures = Q_([1, 10], 'atm')

# Get all of the Species objects defined in the GRI 3.0 mechanism
species = {S.name: S for S in ct.Species.list_from_file('gri30.yaml')}

# Create an IdealGas object with species representing complete combustion
complete_species = [species[S] for S in ('CO', 'O2', 'CO2')]
gas = ct.Solution(thermo='ideal-gas', species=complete_species)

for pressure in pressures:
    gas.TPX = to_si(temperature), to_si(pressure), 'CO:1.0, O2:0.5'

    # Find equilibrium state holding temperature and pressure constant
    gas.equilibrate('TP')

    print(f'Mole fractions at {pressure: .1f}:')
    for sp, mole_fraction in zip(gas.species_names, gas.X):
        print(f'{sp:4}: {mole_fraction: .3f}')
    print()

Mole fractions at  1.0 standard_atmosphere:
CO  :  0.122
O2  :  0.061
CO2 :  0.817

Mole fractions at  10.0 standard_atmosphere:
CO  :  0.061
O2  :  0.030
CO2 :  0.909

These values match what we found using the various methods above.