—

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 cantera as ct
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.lines import Line2D

from rizer.misc.ct_utils import load_goutier2025_mechanism
from rizer.misc.plt_utils import (
    get_species_color,
    get_species_in_latex,
    get_text,
    set_mpl_style,
)

set_mpl_style(nb_columns=1)

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.

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

Define the gas object.

gas = load_goutier2025_mechanism("gas")

Compute and plot the equilibrium states.

# Define the figure and axis.
fig, ax = plt.subplots()


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

    # Define the initial state as pure methane at 1 atm and 273.15 K.
    gas.TPX = 273.15, P, {"CH4": 1}
    states = ct.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]
        color = get_species_color(species)
        ax.plot(
            states.T - 273.15,
            molar_fraction * 100,
            style,
            label=get_species_in_latex(species),
            color=color,
        )

        # 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
            get_text(
                float(x_val),
                float(y_val + 1),
                get_species_in_latex(species),
                ax=ax,
                color=color,
            )

    # 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
        get_text(
            float(x_val),
            float(y_val + 1),
            r"$\mathrm{others}$",
            ax=ax,
        )

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.
/home/runner/work/rizer/rizer/examples/thermodynamics/plot_methane_equilibrium_composition_high_temperature.py:139: UserWarning: ChemEquil::equilibrate: Temperature (273.14999999999986 K) outside valid range of 298.15 K to 100000 K
  gas.equilibrate("TP")
Species C(soot) not in mechanism.
Species C(s) not in mechanism.