Alternative methods for chemical equilibrium#

The methods previously examined for determining the equilibrium composition rely on knowing the chemical reaction(s) occurring, and can involve highly nonlinear equations.

Fortunately, we have methods that do not require knowing what reaction(s) are occurring. We will compare two such solution methods:

Direct minimization of Gibbs free energy
Lagrange’s method of undetermined multipliers

This modules introduces these methods using the same example as Equilibrium via equilibrium constant: consider 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. Find the equilibrium composition in terms of the mole fraction. Assume the mixture behaves as an ideal gas.

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

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

Direct minimization of Gibbs free energy#

One method to finding the equilibrium composition is to directly minimize the Gibbs free energy of the mixture.

The total Gibbs free energy of the mixture is

\[ G = \sum_{i=1}^C n_i \mu_i \;, \]

where \(C\) is the number of components (i.e., chemical species), \(n_i\) is the number of moles of component \(i\), and \(\mu_i\) is the chemical potential of component \(i\). For an ideal gas in a mixture, the chemical potential can be calculated using

\[ \mu_i = \mu_i^{\circ} + R_{\text{univ}} T \ln \left( \frac{y_i P}{P^{\circ}} \right) \;, \]

where \(R_{\text{univ}}\) is the universal gas constant, \(P\) is the mixture pressure, \(P^{\circ}\) is the (standard-state) reference pressure (usually 1 atm or 100 kPa), and \(\mu_i^{\circ}\) is the chemical potential of pure substance \(i\) at temperature \(T\) and reference pressure \(P^{\circ}\), which is the same as the standard-state molar specific Gibbs free energy \(\overline{g}_i^{\circ}\):

\[ \mu_i^{\circ} = \overline{g}_i^{\circ} = \overline{h}_i^{\circ} - T \overline{s}_i^{\circ} \;. \]

This method works by treating this as an optimization problem, where the objective is to minimize \(G\), which is a function of the composition \(n_i\).

Constraints: However, this problem is constrained because the amount of each element must be balanced:

\[ E_j = E_{0, j} \]

where \(E_j = \sum_{i=1}^C n_i e_{i,j}\) is the number of moles of each element \(j\) (\(E\) is the total number of elements), \(E_{0, j} = \sum_{i=1}^C n_{0,i} e_{i,j}\) is the initial number of moles of each element, \(n_{0,i}\) is the initial number of moles of each component \(i\), and \(e_{i,j}\) is the number of moles of element \(j\) in component \(i\) (defined by the chemical formula).

In addition, the number of moles of each component must remain non-negative:

\[ n_i \geq 0 \]

This is thus a constrained optimization problem—we can solve these for simpler problems, but they can become computationally expensive for a larger number of unknowns. For now, we can use the SLSQP optimization method provided by the SciPy minimize function.

The formal statement of our problem is:

\[\begin{split} \min_{n_0, n_1, n_2} \left( n_0 \mu_0 (n_0, n_1, n_2) + n_1 \mu_1 (n_0, n_1, n_2) + n_2 \mu_2 (n_0, n_1, n_2) \right) \\ \text{subject to:} \quad \sum_{i} n_i e_{i,0} - \sum_{i} n_{0,i} e_{i,0} = 0\\ \phantom{subject to:} \quad \sum_{i} n_i e_{i,1} - \sum_{i} n_{1,i} e_{i,0} = 0\\ \phantom{subject to:} \quad n_0 \geq 0 \\ \phantom{subject to:} \quad n_1 \geq 0 \\ \phantom{subject to:} \quad n_2 \geq 0 \end{split}\]

We will need to define three functions:

Evaluate the Gibbs free energy of the mixture,
Evaluate the equality constraints on elemental balances
Evaluate the inequality constraints on numbers of moles

First, let’s input the known information:

# Known information

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

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

# elemental composition of species
elemental_comp = np.array([
    [1, 0, 1], # carbon
    [1, 2, 2], # oxygen
    ])

# initial molar amounts of each element
initial_elements = np.dot(elemental_comp, moles_initial)

def calc_total_gibbs(moles, temperature, pressure, components, gas):
    '''Evaluate Gibbs free energy of mixture, based on component numbers of moles.
    '''
    moles = Q_(moles, 'kmol')
    mole_fractions = moles / np.sum(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')
    
    gas_constant = Q_(ct.gas_constant, 'J/(kmol*K)')
    chemical_potentials = (
        gibbs + gas_constant * temperature * np.log(
            mole_fractions * pressure / Q_(1.0, 'atm')
            )
        )
    
    # scale this result down
    return to_si(np.sum(moles * chemical_potentials)) / 1e6


# We need to define functions for the constraints:

def inequality_cons(x):
    '''Inequality constraint: all numbers of moles must be ≥ 0.
    '''
    return x
    

def equality_cons(x):
    '''Equality constraint: Number of moles of each element remain constant.
    '''
    return np.dot(elemental_comp, x) - initial_elements

# Solve for first pressure

pressure = pressures[0]
gas = ct.Solution('gri30.yaml')

x0 = np.array([0.5, 0.5, 0.5])
sol = minimize(
    calc_total_gibbs, x0, method='SLSQP',
    args=(temperature, pressure, components, gas),
    constraints=[
        {'type': 'eq','fun': equality_cons},
        {'type': 'ineq','fun': inequality_cons}
        ],
    options={'maxiter': 1000}
    )

moles = sol.x
mole_fractions = moles / np.sum(moles)

print('Successful convergence: ', sol.success)

# check constraints
print('All moles non-negative: ', all(moles > 0))
print('All elements balanced: ', all(equality_cons(moles) == 0))

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

Successful convergence:  True
All moles non-negative:  True
All elements balanced:  False

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

/usr/share/miniconda/envs/thermo-env/lib/python3.9/site-packages/pint/numpy_func.py:303: RuntimeWarning: invalid value encountered in log
  result_magnitude = func(*stripped_args, **stripped_kwargs)
/usr/share/miniconda/envs/thermo-env/lib/python3.9/site-packages/pint/numpy_func.py:303: RuntimeWarning: divide by zero encountered in log
  result_magnitude = func(*stripped_args, **stripped_kwargs)
/usr/share/miniconda/envs/thermo-env/lib/python3.9/site-packages/pint/quantity.py:1313: RuntimeWarning: invalid value encountered in multiply
  magnitude = magnitude_op(new_self._magnitude, other._magnitude)

# Now try next pressure

pressure = pressures[1]
gas = ct.Solution('gri30.yaml')

x0 = np.array([0.5, 0.5, 0.5])
sol = minimize(
    calc_total_gibbs, x0, method='SLSQP',
    args=(temperature, pressure, components, gas),
    constraints=[
        {'type': 'eq','fun': equality_cons},
        {'type': 'ineq','fun': inequality_cons}
        ],
    options={'maxiter': 1000}
    )

moles = sol.x
mole_fractions = moles / np.sum(moles)

print('Successful convergence: ', sol.success)

# check constraints
print('All moles non-negative: ', all(moles > 0))
print('All elements balanced: ', all(equality_cons(moles) == 0))

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

Successful convergence:  True
All moles non-negative:  True
All elements balanced:  True

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

These results match the values we found previously—whew! 😅

Lagrange’s method of undetermined multipliers#

This method converts the problem into a system of algebraic equations, where the number of equations equal the number of unknowns. It does this by introducing a set of unknown multipliers, \(\lambda_j\), with one for each element in the system.

Then, the system of equations we need to solve includes the element balances and equations involving the multipliers:

\[\begin{split} \sum_{i=1}^C n_i e_{i,j} - \sum_{i=1}^C n_{0,i} e_{i,j} = 0 \quad \text{for } j=1, \ldots, E \;, \\ \mu_i + \sum_{j=1}^E \lambda_j e_{i,j} = 0 \quad \text{for } i=1, \ldots, C \;, \\ \end{split}\]

where the unknowns are the numbers of moles for each compound \(n_i\) where \(i = 1, \ldots, C\) and the multipliers for each element \(\lambda_j\) where \(j = 1, \ldots, E\). In this system, \(e_{i,j}\) is the number of moles of element \(j\) in component \(i\), \(n_{0,i}\) is the initial number of moles of component \(i\), \(\mu_i\) is the chemical potential of component \(i\), \(E\) is the number of elements, and \(C\) is the number of components (chemical species).

The chemical potentials can be calculated for each component of an ideal gas:

\[ \mu_i = \mu_i^{\circ} + R_{\text{univ}} T \ln \left( \frac{y_i P}{P^{\circ}} \right) \;, \]

\[ \mu_i^{\circ} = \overline{g}_i^{\circ} = \overline{h}_i^{\circ} - T \overline{s}_i^{\circ} \;. \]

We can evaluate \(\overline{g}_i^{\circ} (T)\) using a Cantera Solution object and specifying the temperature, pressure (using the 1 atm reference), and composition of each component as a pure substance.

# Known information

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

# Elemental makeup of components
elemental_comp = np.array([
    [1, 0, 1], # carbon
    [1, 2, 2], # oxygen
    ])

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

def lagrange_system(x, temperature, pressure, components, 
                    gas, elemental_comp, moles_initial):
    '''System of equations for Lagrange multiplier approach.
    '''
    moles = np.array([x[0], x[1], x[2]])
    multipliers = np.array([x[3], x[4]])
    
    mole_fractions = moles / np.sum(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')
    
    gas_constant = Q_(ct.gas_constant, 'J/(kmol*K)')
    chemical_potentials = (
        gibbs + gas_constant * temperature * np.log(
            mole_fractions * pressure / Q_(1.0, 'atm')
            )
        )

    # initial molar amounts of each element
    initial_moles_elements = np.dot(elemental_comp, moles_initial)
    moles_elements = np.dot(elemental_comp, moles)
    
    # We can take advantage of element-wise operations with these arrays,
    # and concisely evaluate all the equations
    element_equations = moles_elements - initial_moles_elements
    multiplier_equations = to_si(
        chemical_potentials + 
        np.dot(multipliers, elemental_comp) * Q_('J/kmol')
        )
    
    # Return the set of equations joined together
    return np.concatenate((element_equations, multiplier_equations))

After setting up the function to evaluate the system of equations, we can solve for the equilibrium composition at the first pressure using the root function, with the lm (Levenberg-Marquardt) method.

We do need to specify some initial guess values for each of the unknowns; while guess values for the numbers of moles of each component may be straightforward (e.g., typically around one), the Lagrange multipliers are more abstract and may take some trial and error.

# Solve at first pressure

pressure = pressures[0]
gas = ct.Solution('gri30.yaml')

# initial guesses
x0 = [1.0, 1.0, 1.0, 1e6, 1e6]

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

print('Root-finding algorithm success: ', sol.success)
print(f'Function evaluation (should be small): {sol.fun}')
print('Number of function evaluations: ', sol.nfev)
print()

moles = sol.x[0:3]
mole_fractions = moles / np.sum(moles)
print(f'Mole fractions at {pressure: .1f}:')
for idx, comp in enumerate(components):
    print(f'{comp}: {mole_fractions[idx]: .3f}')

Root-finding algorithm success:  True
Function evaluation (should be small): [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 2.38418579e-10]
Number of function evaluations:  84

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

pressure = pressures[1]
gas = ct.Solution('gri30.yaml')

x0 = [1.0, 1.0, 1.0, 1e6, 1e6]
sol = root(
    lagrange_system, x0, method='lm',
    args=(temperature, pressure, components, gas, elemental_comp, moles_initial)
    )

print('Root-finding algorithm success: ', sol.success)
print(f'Function evaluation (should be near zero): {sol.fun}')
print('Number of function evaluations: ', sol.nfev)
print()

moles = sol.x[0:3]
mole_fractions = moles / np.sum(moles)
print(f'Mole fractions at {pressure: .1f}:')
for idx, comp in enumerate(components):
    print(f'{comp}: {mole_fractions[idx]: .3f}')

Root-finding algorithm success:  True
Function evaluation (should be near zero): [-1.11022302e-16  0.00000000e+00  1.19209290e-10  0.00000000e+00
  1.19209290e-10]
Number of function evaluations:  69

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

As expected, this approach also produces the same equilibrium composition! 🎉

Computational Thermodynamics

Alternative methods for chemical equilibrium

Contents

Alternative methods for chemical equilibrium#

Direct minimization of Gibbs free energy#

Lagrange’s method of undetermined multipliers#