Adiabatic flame temperature#

Consider a mixture of hydrogen and oxygen initially at 1000 K and 10 bar, which is ignited by a spark. This mixture reacts according to

\[ \text{H}_2 + 0.5 \text{O}_2 \leftrightarrow \text{H}_2 \text{O} \]

and proceeds to equilibrium at a constant pressure, adiabatic process. Initially the mixture has twice as much oxygen than hydrogen (by mole) and no water. Assume the mixture follows the ideal gas law.

Problem: Find the equilibrium composition and temperature, using the Lagrange multiplier method.

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

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

After importing the necessary modules, we should specify the knowns:

temperature_initial = Q_(1000, 'K')
pressure = Q_(10, 'bar')

components = ['H2', 'O2', 'H2O']
moles_initial = np.array([1.0, 2.0, 0.0]) * Q_('kmol')

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

The system of equations we solve will include the element balances and equations involving the multipliers, and also a constraint of constant enthalpy (\(H_1 = H_2\)):

\[\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 \;, \\ \sum_{i=1}^C n_{0,i} \overline{h}_{i, T_{0}} = \sum_{i=1}^C n_{i} \overline{h}_{i, T_f} \;, \end{split}\]

where the unknowns are the numbers of moles for each compound \(n_i\), the multipliers for each element \(\lambda_j\), and the final temperature \(T_f\). 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\), \(\overline{h}_{i, T}\) is the molar specific enthalpy of component \(i\) evaluated at temperature \(T\), \(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) \;, \]

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} \;. \]

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

def lagrange_system(x, pressure, components, gas, elemental_comp, 
                    temperature_initial, moles_initial):
    '''System of equations for reaction coordinate and equilibrium composition.
    '''
    moles = np.array([x[0], x[1], x[2]]) * Q_('kmol')
    multipliers = np.array([x[3], x[4]]) * Q_('J/kmol')
    temperature = Q_(x[5], 'K')
    
    mole_fractions = to_si(moles / np.sum(moles))
    
    # get standard-state Gibbs free energy and enthalpy of each component
    gibbs = np.zeros(len(components))
    enthalpies_final = np.zeros(len(components))
    enthalpies_initial = 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
        enthalpies_final[idx] = gas.enthalpy_mole
        
        gas.TPX = (
            to_si(temperature_initial), to_si(Q_(1, 'atm')),
            f'{comp}:1.0'
            )
        enthalpies_initial[idx] = gas.enthalpy_mole
        
    gibbs *= Q_('J/kmol')
    enthalpies_final *= Q_('J/kmol')
    enthalpies_initial *= Q_('J/kmol')
    
    # Calculate the chemical potentials at current pressure and temperature
    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
    # base SI units are in mol, not kmol, after conversion
    initial_moles_elements = Q_(
        np.dot(elemental_comp, to_si(moles_initial)), 'mol'
        )
    moles_elements = Q_(
        np.dot(elemental_comp, to_si(moles)), 'mol'
        )
    
    enthalpy_initial = np.sum(moles_initial * enthalpies_initial)
    enthalpy_final = np.sum(moles * enthalpies_final)
    
    return [
        to_si(moles_elements[0] - initial_moles_elements[0]),
        to_si(moles_elements[1] - initial_moles_elements[1]),
        to_si(chemical_potentials[0] + np.sum(multipliers * elemental_comp[:,0])),
        to_si(chemical_potentials[1] + np.sum(multipliers * elemental_comp[:,1])),
        to_si(chemical_potentials[2] + np.sum(multipliers * elemental_comp[:,2])),
        to_si(enthalpy_final - enthalpy_initial)
        ]

gas = ct.Solution('gri30.yaml')

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

print('Root-finding algorithm success: ', sol.success)
print('Function evaluation (should be small): \n' +
      ', '.join([f'{val:.4e}' for val in sol.fun])
      )
print()

moles = sol.x[:3]
mole_fractions = moles / np.sum(moles)
print(f'Mole fractions at equilibrium:')
for idx, comp in enumerate(components):
    print(f'{comp:4}: {mole_fractions[idx]: .4f}')
    
temperature_final = Q_(sol.x[-1], 'K')
print(f'Temperature at equilibrium: {temperature_final: .2f}')

/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)

Root-finding algorithm success:  True
Function evaluation (should be small): 
0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -2.3842e-10, -4.4703e-08

Mole fractions at equilibrium:
H2  :  0.0136
O2  :  0.6027
H2O :  0.3837
Temperature at equilibrium: 3208.46 kelvin

Compare to Cantera equilibrium#

We can compare this approach to the built-in equilibrium solver in Cantera, which uses a different (but related) element potential method:

# 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 ('H2', 'O2', 'H2O')]
gas = ct.Solution(thermo='IdealGas', species=complete_species)

gas.TPX = to_si(temperature_initial), to_si(pressure), 'O2:2.0, H2:1.0'
gas.equilibrate('HP')

print(f'Adiabatic flame temperature: {gas.T: .2f} K')

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

Adiabatic flame temperature:  3208.46 K
Mole fractions at equilibrium:
H2  :  0.0136
O2  :  0.6027
H2O :  0.3837

Both methods produce exactly the same values! 🔥

Computational Thermodynamics

Adiabatic flame temperature

Contents

Adiabatic flame temperature#

Compare to Cantera equilibrium#