Adiabatic flame temperature - Computational Thermodynamics

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}

(1)

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$ ):

\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} \;,

(2)

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) \;,

(3)

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

(4)

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}')

Root-finding algorithm success:  True
Function evaluation (should be small): 
2.2737e-13, 0.0000e+00, -1.1921e-10, 1.1921e-10, 0.0000e+00, -1.4901e-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='ideal-gas', 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 (to two decimal points)! 🔥