# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# # 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$):
#
# $$
# \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.

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

# ## 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](https://cantera.org/tutorials/python-tutorial.html#chemical-equilibrium):

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

# Both methods produce exactly the same values! 🔥

StarcoderdataJupyterScriptsDedupFiltered