—

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
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 = 5000  # °C
temperatures = np.linspace(T_min + 273.15, T_max + 273.15, 100)  # K

# Include solid carbon in the computation.
with_carbon = 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.

# Nice labels for the species.
labels = {
    "H": r"$\mathrm{H}$",
    "H2": r"$\mathrm{H}_2$",
    "C": r"$\mathrm{C(g)}$",
    "C(s)": r"$\mathrm{C(s)}$",
    "CH4": r"$\mathrm{CH}_4$",
    "C2": r"$\mathrm{C}_2$",
    "C2H": r"$\mathrm{C}_2\mathrm{H}$",
    "C2H2": r"$\mathrm{C}_2\mathrm{H}_2$",
    "C2H4": r"$\mathrm{C}_2\mathrm{H}_4$",
    "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.

if with_carbon:
    mechanism = "Fincke_GRC.yaml"
    reference = "Fincke2002_equilibrium_fig1.csv"
else:
    mechanism = "Fincke_GRC_no_solid_carbon.yaml"
    reference = "Fincke2002_equilibrium_fig2.csv"

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

ref_data = np.genfromtxt(
    get_path_to_data("thermo", 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(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.
    # 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=labels[species],
        )

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

plt.show()