.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/kinetics/explosions/plot_01_h2_o2_explosion_peninsula.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_kinetics_explosions_plot_01_h2_o2_explosion_peninsula.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_kinetics_explosions_plot_01_h2_o2_explosion_peninsula.py:


The hydrogen-oxygen explosion peninsula
=======================================

In the H2/O2 chain the branching step :math:`H + O_2 \to OH + O` (net two
extra H atoms per event) competes with H loss at the vessel wall and with
the three-body :math:`H + O_2 + M \to HO_2 + M`. Semenov's net branching
factor
(:meth:`~chemistrykit.kinetics.ChainBranchingExplosion.branching_factor`),

.. math::

    \varphi = 2k_bx_{O_2}c - k_w - k_tx_{O_2}c^2,

is positive, and the mixture explodes, only between two pressures: the
first limit, where branching beats wall loss, and the second, where
three-body termination catches up. Because :math:`k_b` has a large
activation energy and :math:`k_t` has almost none, the window widens
with temperature and closes at a tip. That traces the famous explosion
peninsula. The Arrhenius parameters are of the size measured for these
reactions; the wall rate depends on the vessel. The third, thermal limit
at higher pressure comes from self-heating and is not part of this
isothermal model.

.. GENERATED FROM PYTHON SOURCE LINES 27-67

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from chemistrykit.constants import R
    from chemistrykit.kinetics import ChainBranchingExplosion, arrhenius_rate_constant

    model = ChainBranchingExplosion(
        k_branch=lambda T: arrhenius_rate_constant(A=1.9e8, Ea=70.3e3, T=T),  # m^3 mol^-1 s^-1
        k_wall=100.0,  # s^-1
        k_termination=1.0e4,  # m^6 mol^-2 s^-1
        x_O2=1.0 / 3.0,
    )

    T = np.linspace(600.0, 900.0, 301)
    limits = [model.explosion_limits(t) for t in T]
    mask = np.array([lim is not None for lim in limits])
    P1 = np.array([lim[0] if lim else np.nan for lim in limits])
    P2 = np.array([lim[1] if lim else np.nan for lim in limits])
    T_tip = T[mask][0]

    fig, axes = plt.subplots(1, 2, figsize=(12, 4.8))
    axes[0].fill_between(T[mask] - 273.15, P1[mask] / 133.322, P2[mask] / 133.322, color="darkorange", alpha=0.35, label="explosion")
    axes[0].semilogy(T - 273.15, P1 / 133.322, color="darkorange", label="first limit")
    axes[0].semilogy(T - 273.15, P2 / 133.322, color="firebrick", label="second limit")
    axes[0].set_xlabel("temperature (degC)")
    axes[0].set_ylabel("pressure (Torr)")
    axes[0].set_title("Explosion peninsula of 2H2 + O2")
    axes[0].legend(loc="lower right")

    T0 = 800.0
    p1, p2 = model.explosion_limits(T0)
    P = np.logspace(np.log10(p1) - 1, np.log10(p2) + 0.5, 400)
    axes[1].semilogx(P / 133.322, model.branching_factor(P, T0), color="black")
    axes[1].axhline(0.0, color="gray", linewidth=0.8)
    axes[1].axvspan(p1 / 133.322, p2 / 133.322, color="darkorange", alpha=0.25)
    axes[1].set_xlabel("pressure (Torr)")
    axes[1].set_ylabel(r"net branching factor $\varphi$ (s$^{-1}$)")
    axes[1].set_title(f"Semenov's criterion at {T0 - 273.15:.0f} degC")
    fig.tight_layout()




.. image-sg:: /api/gallery/kinetics/explosions/images/sphx_glr_plot_01_h2_o2_explosion_peninsula_001.png
   :alt: Explosion peninsula of 2H2 + O2, Semenov's criterion at 527 degC
   :srcset: /api/gallery/kinetics/explosions/images/sphx_glr_plot_01_h2_o2_explosion_peninsula_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 68-70

Either side of the second limit, the H-atom population either levels
off at :math:`w_0/|\varphi|` or grows exponentially:

.. GENERATED FROM PYTHON SOURCE LINES 70-78

.. code-block:: Python


    for P_test in (0.5 * p2, 1.5 * p2):
        n = model.carrier_concentration([0.0, 0.005, 0.01], P_test, T0, w0=1e-6)
        print(f"P = {P_test / 133.322:6.1f} Torr, phi = {float(model.branching_factor(P_test, T0)):8.1f} s^-1, n(t) = {np.array2string(n, precision=3)}")
    print(f"peninsula tip near {T_tip - 273.15:.0f} degC; second limit at {T0 - 273.15:.0f} degC = {p2 / 133.322:.0f} Torr")
    print(f"second-limit concentration {p2 / (R * T0):.3f} mol/m^3; wall-free limit 2 k_b / k_t = {2 * model.k_branch(T0) / 1.0e4:.3f}")

    plt.show()




.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    P =   23.6 Torr, phi =    694.1 s^-1, n(t) = [0.000e+00 4.488e-08 1.488e-06]
    P =   70.7 Torr, phi =  -2182.3 s^-1, n(t) = [0.000e+00 4.582e-10 4.582e-10]
    peninsula tip near 456 degC; second limit at 527 degC = 47 Torr
    second-limit concentration 0.945 mol/m^3; wall-free limit 2 k_b / k_t = 0.977





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.189 seconds)


.. _sphx_glr_download_api_gallery_kinetics_explosions_plot_01_h2_o2_explosion_peninsula.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_01_h2_o2_explosion_peninsula.ipynb <plot_01_h2_o2_explosion_peninsula.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_01_h2_o2_explosion_peninsula.py <plot_01_h2_o2_explosion_peninsula.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_01_h2_o2_explosion_peninsula.zip <plot_01_h2_o2_explosion_peninsula.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_