—

Plot CH₄ chemical equilibrium vs temperature for various pressures (at high temperatures).

This script computes the chemical equilibrium of methane (CH₄) for various temperatures and pressures using Cantera. It plots the molar concentration of the main species in the equilibrium state as a function of temperature.

The mechanism used is the one from ./data/mechanisms/Goutier2025. Thermodynamic data comes from [GRCDatabase] and [Burcat], and have been extended to 50 000 K, by fixing Cantera’s extrapolation at high temperature. Moreover, C+, C++, C+++ and electrons have been added to the mechanism, with thermodynamic data computed using [MinPlasCalc] (& [minplascalcPaper]).

Note that Cantera extrapolates the thermodynamic data if T is outside the range of validity of the thermodynamic data. This can lead to unphysical results, especially at high temperatures. See the following issue for more details: Cantera/cantera#270

Tags: cantera chemical equilibrium methane CH4 high temperature

Import the required libraries.

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from cantera import Solution, SolutionArray
from matplotlib.lines import Line2D

from rizer.misc.utils import get_path_to_data

sns.set_theme("talk")

Define the parameters for the plot.

# Define the temperature range.
T_min = 0  # [°C]
T_max = 50_000  # [°C]
temperatures = np.linspace(T_min + 273.15, T_max + 273.15, 1000)  # [K]


# Define the species to plot.
species_to_plot = [
    "e-",
    "H",
    "H+",
    "H2",
    # "H2+",
    "C",
    "C(soot)",
    "C(s)",
    "C+",
    "C++",
    "C+++",
    # "CH",
    # "CH+",
    # "CH2",
    # "CH2+",
    # "CH3",
    # "CH3+",
    "CH4",
    # "CH4+",
    # "C2",
    # "C2+",
    # "C2H",
    # "C2H+",
    "C2H2",
    # "C2H2+",
    # "C2H3",
    # "C2H3+",
    # "C2H4",
    # "C2H4+",
    # "C2H5",
    # "C2H5+",
    # "C2H6",
    # "C2H6+",
    # "C6H6",
]


# Define the pressures to plot.
pressures = [1, 100]  # atm

# Style linestyle for the various pressures.
linestyles = ["-", "--"]

if len(linestyles) != len(pressures):
    raise ValueError(
        "The number of linestyles must be equal to the number of pressures."
    )

Define some parameters for the plot.

# Nice labels for the species.
labels = {
    "e-": r"$\mathrm{e}^-$",
    "H": r"$\mathrm{H}$",
    "H+": r"$\mathrm{H}^+$",
    "H2": r"$\mathrm{H}_2$",
    "H2+": r"$\mathrm{H}_2^+$",
    "C": r"$\mathrm{C}$",
    "C(soot)": r"$\mathrm{C}(\mathrm{soot})$",
    "C(s)": r"$\mathrm{C}(s)$",
    "C+": r"$\mathrm{C}^+$",
    "C++": r"$\mathrm{C}^{++}$",
    "C+++": r"$\mathrm{C}^{3+}$",
    "CH": r"$\mathrm{CH}$",
    "CH+": r"$\mathrm{CH}^+$",
    "CH2": r"$\mathrm{CH}_2$",
    "CH2+": r"$\mathrm{CH}_2^+$",
    "CH3": r"$\mathrm{CH}_3$",
    "CH3+": r"$\mathrm{CH}_3^+$",
    "CH4": r"$\mathrm{CH}_4$",
    "CH4+": r"$\mathrm{CH}_4^+$",
    "C2": r"$\mathrm{C}_2$",
    "C2+": r"$\mathrm{C}_2^+$",
    "C2H": r"$\mathrm{C}_2\mathrm{H}$",
    "C2H+": r"$\mathrm{C}_2\mathrm{H}^+$",
    "C2H2": r"$\mathrm{C}_2\mathrm{H}_2$",
    "C2H2+": r"$\mathrm{C}_2\mathrm{H}_2^+$",
    "C2H3": r"$\mathrm{C}_2\mathrm{H}_3$",
    "C2H3+": r"$\mathrm{C}_2\mathrm{H}_3^+$",
    "C2H4": r"$\mathrm{C}_2\mathrm{H}_4$",
    "C2H4+": r"$\mathrm{C}_2\mathrm{H}_4^+$",
    "C2H5": r"$\mathrm{C}_2\mathrm{H}_5$",
    "C2H5+": r"$\mathrm{C}_2\mathrm{H}_5^+$",
    "C2H6": r"$\mathrm{C}_2\mathrm{H}_6$",
    "C2H6+": r"$\mathrm{C}_2\mathrm{H}_6^+$",
    "C6H6": r"$\mathrm{C}_6\mathrm{H}_6$",
}


# Create legends for the various pressures.
legend_elements = [
    Line2D([0], [0], color="k", linestyle=linestyles[i], lw=2, label=f"{p} atm")
    for i, p in enumerate(pressures)
]

Define the gas object.

gas = Solution(
    get_path_to_data(
        "mechanisms",
        "Goutier2025",
        "CH4_to_C2H2.yaml",
    ),
    name="gas",
)

Compute and plot the equilibrium states.

# Define the figure and axis.
fig, ax = plt.subplots(figsize=(12, 8), layout="constrained")


for i, pressure_atm in enumerate(pressures):
    # Convert the pressure to Pa.
    P = pressure_atm * 101325  # Pa

    # Define the initial state as pure methane at 1 atm and 273.15 K.
    gas.TPX = 273.15, P, {"CH4": 1}
    states = SolutionArray(gas)

    # Calculate the equilibrium states.
    for T in temperatures:
        gas.TP = T, P
        gas.equilibrate("TP")
        states.append(gas.state)

    # Plot the results.
    style = linestyles[i]  # Linestyle for the current pressure.
    others_x = np.ones_like(temperatures)  # Total molar fraction of the other species.

    for species in species_to_plot:
        # Molar fraction of the species.
        if species not in gas.species_names:
            print(f"Species {species} not in mechanism.")
            continue

        molar_fraction = states(species).X[:, 0]
        ax.plot(
            states.T - 273.15,
            molar_fraction * 100,
            style,
            label=labels[species],
        )

        # Update the total molar fraction of the other species.
        others_x -= molar_fraction

        # Add labels on curves, at the maximum of the curves,
        # with the same color as the curve, only for the first pressure.
        if i == 0:
            y_val = np.max(molar_fraction) * 100
            x_val = states.T[np.argmax(states(species).X)] - 273.15
            ax.text(
                x=x_val,
                y=y_val + 1,
                s=labels[species],
                color=ax.lines[-1].get_color(),
                fontsize=30,
                ha="center",
            )

    # Plot the total molar fraction of the other species.
    ax.plot(
        states.T - 273.15,
        others_x * 100,
        style,
        color="k",
        label=r"$\mathrm{others}$",
    )
    # Add label on the curve, at the maximum of the curve,
    # with the same color as the curve, only for the first pressure.
    if i == 0:
        y_val = float(np.max(others_x) * 100)
        x_val = states.T[np.argmax(others_x)] - 273.15
        ax.text(
            x=x_val,
            y=y_val + 1,
            s=r"$\mathrm{others}$",
            color=ax.lines[-1].get_color(),
            fontsize=30,
            ha="center",
        )

    # Reset the color cycle.
    ax.set_prop_cycle(None)  # type: ignore[call-overload]

ax.legend(handles=legend_elements, loc="upper right")
ax.set_xlabel("T [°C]")
ax.set_ylabel("x [%]")

title = "Equilibrium molar concentration [%] vs temperature [°C] of methane "
ax.set_title(title, fontsize=20)
ax.set_xlim(T_min, T_max)

plt.show()

Species C(soot) not in mechanism.
Species C(s) not in mechanism.
Species C(soot) not in mechanism.
Species C(s) not in mechanism.