<jupyter_start><jupyter_text>Adiabatic flame temperatureConsider 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.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>After importing the necessary modules, we should specify the knowns:<jupyter_code>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
    ])<jupyter_output><empty_output><jupyter_text>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} \;,$$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.<jupyter_code>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.cti')

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}')<jupyter_output>Root-finding algorithm success:  True
Function evaluation (should be small): 
0.0000e+00, 0.0000e+00, 1.1921e-10, 1.1921e-10, 2.3842e-10, 2.2352e-08

Mole fractions at equilibrium:
H2  :  0.0136
O2  :  0.6027
H2O :  0.3837
Temperature at equilibrium: 3208.46 kelvin<jupyter_text>Compare to Cantera equilibriumWe can compare this approach to the built-in equilibrium solver in Cantera, which uses  a different (but related) [element potential method](https://cantera.org/tutorials/python-tutorial.htmlchemical-equilibrium):<jupyter_code># Get all of the Species objects defined in the GRI 3.0 mechanism
species = {S.name: S for S in ct.Species.listFromFile('gri30.cti')}

# 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}')<jupyter_output>Adiabatic flame temperature:  3208.46 K
Mole fractions at equilibrium:
H2  :  0.0136
O2  :  0.6027
H2O :  0.3837
{'fraction_non_alphanumeric': 0.1359520639147803, 'fraction_numerical': 0.027430093209054592, 'mean_word_length': 3.980769230769231, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '85071', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2715, 'n_tokens_neox': 2533, 'n_words': 730, 'path': 'book/content/chemical-equilibrium/adiabatic-flame-temperature.ipynb', 'repo_name': 'kyleniemeyer/computational-thermo', 'size': '7510'}