—

Plot CH₄ chemical equilibrium vs temperature for various pressures (including solid carbon).

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 Fincke GRC mechanism, which includes solid carbon (C(s)). The script can be modified to exclude solid carbon from the computation.

Thermodynamic data comes from [GRCDatabase] and [Burcat].

Tags: cantera chemical equilibrium methane CH4

Import the required libraries.

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

from rizer.misc.plt_utils import get_species_in_latex, set_mpl_style
from rizer.misc.utils import get_path_to_data

set_mpl_style()

Define the parameters for the plot.

# Define the temperature range.
T_min = 0  # °C
T_max = 5000  # °C
temperatures = np.linspace(T_min + 273.15, T_max + 273.15, 100)  # K

# Include solid carbon in the computation.
with_carbon = False

# Plot reference data from Fincke 2002.
plot_reference_data = False

# Define the species to plot.
species_to_plot = ["CH4", "C6H6", "C2H2", "C2H", "C2", "H2", "H", "C", "C(s)"]

if not with_carbon and "C(s)" in species_to_plot:
    # Remove solid carbon if not included in the computation.
    print("Solid carbon is not included in the computation. Removing it from the plot.")
    species_to_plot.remove("C(s)")


# Define the pressures to plot.
pressures = [
    1,
    # 10,
]  # 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."
    )

Solid carbon is not included in the computation. Removing it from the plot.

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")
    for i, p in enumerate(pressures)
]

Define the gas object.

mechanism = "Fincke_GRC.yaml"
if with_carbon:
    phase = "gas_with_soot"
    reference = "papers/Fincke2002/equilibrium_fig1.csv"
else:
    phase = "gas_without_soot"
    reference = "papers/Fincke2002/equilibrium_fig2.csv"

gas = Solution(get_path_to_data("mechanisms", mechanism), name=phase)

if plot_reference_data:
    ref_data = np.genfromtxt(
        get_path_to_data(reference),
        delimiter=",",
        skip_header=3,
        encoding="utf-8",
        missing_values="nan",
    )
    ref_data_dict: dict[str, np.ndarray] = {}
    if with_carbon:
        ref_data_dict["T"] = ref_data[:, 0]  # [°C]
        ref_data_dict["C"] = ref_data[:, 1]  # moles of C(g)
        ref_data_dict["C(s)"] = ref_data[:, 2]  # moles of C(s)
        ref_data_dict["C2"] = ref_data[:, 3]  # moles of C2
        ref_data_dict["C2H"] = ref_data[:, 4]  # moles of C2H
        ref_data_dict["C2H2"] = ref_data[:, 5]  # moles of C2H2
        ref_data_dict["CH4"] = ref_data[:, 6]  # moles of CH4
        ref_data_dict["H"] = ref_data[:, 7]  # moles of H
        ref_data_dict["H2"] = ref_data[:, 8]  # moles of H2
        # Normalize the reference data to get molar fractions.
        total_moles = np.nansum(
            [ref_data_dict[s] for s in ref_data_dict if s != "T"], axis=0
        )
        for species in ref_data_dict:
            if species != "T":
                ref_data_dict[species] /= total_moles
    else:
        ref_data_dict["T"] = ref_data[:, 0]  # [°C]
        ref_data_dict["C"] = ref_data[:, 1]  # moles of C(g)
        ref_data_dict["C2"] = ref_data[:, 2]  # moles of C2
        ref_data_dict["C2H"] = ref_data[:, 3]  # moles of C2H
        ref_data_dict["C2H2"] = ref_data[:, 4]  # moles of C2H2
        ref_data_dict["C2H4"] = ref_data[:, 5]  # moles of C2H4
        ref_data_dict["C6H6"] = ref_data[:, 6]  # moles of C6H6
        ref_data_dict["CH4"] = ref_data[:, 7]  # moles of CH4
        ref_data_dict["H"] = ref_data[:, 8]  # moles of H
        ref_data_dict["H2"] = ref_data[:, 9]  # moles of H2
        # Normalize the reference data to get molar fractions.
        total_moles = np.nansum(
            [ref_data_dict[s] for s in ref_data_dict if s != "T"], axis=0
        )
        for species in ref_data_dict:
            if species != "T":
                ref_data_dict[species] /= total_moles

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 * 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.
    # Total molar fraction of the other species.
    others_x = np.ones_like(temperatures, dtype=float)

    for species in species_to_plot:
        # Molar fraction of the species.
        molar_fraction = states(species).X[:, 0]
        if species == "C(s)":
            # Include soot carbon as solid carbon in the molar fraction.
            molar_fraction += states("C(soot)").X[:, 0]
        ax.plot(
            states.T - 273.15,
            molar_fraction * 100,
            style,
            label=get_species_in_latex(species),
        )

        if plot_reference_data and species in ref_data_dict:
            # Plot reference data as circles.
            ax.plot(
                ref_data_dict["T"][::10],
                ref_data_dict[species][::10] * 100,
                "o",
                color=ax.lines[-1].get_color(),
                markersize=5,
                label=None,
            )

        # 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:
            # Take the maximum value between ref_data and computed data.
            y_val_computed: float = np.max(molar_fraction) * 100
            x_val_computed: float = states.T[np.argmax(states(species).X)] - 273.15
            y_val_ref: float | None = None
            if plot_reference_data and species in ref_data_dict:
                y_val_ref = np.nanmax(ref_data_dict[species]) * 100
                x_val_ref = ref_data_dict["T"][np.nanargmax(ref_data_dict[species])]
            if y_val_ref is not None and y_val_ref > y_val_computed:
                y_val = y_val_ref
                x_val = x_val_ref
            else:
                y_val = y_val_computed
                x_val = x_val_computed
            ax.text(
                x=x_val,
                y=y_val + 2,
                s=get_species_in_latex(species),
                color=ax.lines[-1].get_color(),
                fontsize=30,
                ha="center",
                zorder=10,  # Ensure the text is on top of the curves.
                # Add a white background to the text for better visibility.
                bbox=dict(facecolor="white", edgecolor="none", alpha=0.7, pad=0.5),
            )

    # 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",
            zorder=10,  # Ensure the text is on top of the curves.
            # Add a white background to the text for better visibility.
            bbox=dict(facecolor="white", edgecolor="none", alpha=0.7, pad=0.5),
        )

    # 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 "
title += "(with graphite)" if with_carbon else "(without graphite)"
ax.set_title(title, fontsize=20)
ax.set_xlim(T_min, T_max)

plt.show()