—

Plot the energy cost of producing H₂ from CH₄ as a function of temperature.

The energy cost of producing hydrogen (H₂) from methane (CH₄) is computed as the ratio of the enthalpy change to the amount of hydrogen produced. It is plotted in kWh/kgH₂ as a function of temperature for different cases.

There are several ways to compute the enthalpy change, depending on the final state of the system.

The initial state is defined as pure methane at 1 atm and 298.15 K. The final states are:

• case 1a: equilibrium mixture at T [K] and 1 atm.

• case 1b: equilibrium mixture at T [K] and 1 atm, cooled down to 298.15 K.

• case 2a: equilibrium mixture at T [K] and 1 atm, without solid carbon in the equilibrium computation.

• case 2b: equilibrium mixture at T [K] and 1 atm, without solid carbon in the equilibrium computation, cooled down to 298.15 K.

Cases 1b and 2b include the heat recovery from cooling the products down to 298.15 K, which reduces the energy cost.

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

Tags: cantera chemical equilibrium methane CH4

Import the required libraries.

from dataclasses import dataclass

import cantera as ct
import matplotlib.pyplot as plt
import numpy as np

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

set_mpl_style()

Define the parameters for the plot.

# Define the reference temperature and pressure.
T_ref = 298.15  # K
pressure_ref = ct.one_atm  # Pa

# Initial composition.
initial_composition = {"CH4": 1.0}

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


@dataclass
class Case:
    with_carbon: bool
    heat_recovery: bool
    energy_costs: list[float]
    label: str
    xy: tuple[float, float]


# Cases to plot.
cases: dict[str, Case] = {
    "case 1a": Case(
        with_carbon=True,
        heat_recovery=False,
        energy_costs=[],
        label="C(s)",
        xy=(500, 10),
    ),
    "case 1b": Case(
        with_carbon=True,
        heat_recovery=True,
        energy_costs=[],
        label="C(s) + heat recovery",
        xy=(1500, 3),
    ),
    "case 2a": Case(
        with_carbon=False,
        heat_recovery=False,
        energy_costs=[],
        label="without C(s)",
        xy=(1500, 30),
    ),
    "case 2b": Case(
        with_carbon=False,
        heat_recovery=True,
        energy_costs=[],
        label="without C(s) + heat recovery",
        xy=(1800, 19),
    ),
}

# Define the gas object.
path_to_mechanism = get_path_to_data("mechanisms", "Fincke_GRC.yaml")
gas_with_carbon = ct.Solution(path_to_mechanism, name="gas_with_soot")
gas_without_carbon = ct.Solution(path_to_mechanism, name="gas_without_soot")

Compute the energy cost for each case.

for case_name, case_params in cases.items():
    print(f"Computing energy cost for {case_name}...")
    with_carbon = case_params.with_carbon
    heat_recovery = case_params.heat_recovery

    gas = gas_with_carbon if with_carbon else gas_without_carbon

    energy_costs = []

    # Set the initial state.
    gas.TPX = T_ref, pressure_ref, initial_composition

    # Get the initial enthalpy of the system.
    h_initial = gas.enthalpy_mass  # [J/kg]

    for T in temperatures:
        # Set the temperature and pressure for the equilibrium computation.
        gas.TP = T, pressure_ref

        # Compute the equilibrium state at the final temperature and pressure.
        gas.equilibrate("TP")

        # Get the enthalpy at the equilibrium state at T.
        h_T_eq = gas.enthalpy_mass  # [J/kg]

        # Compute the enthalpy change of the reaction (eq. @298.15 K to eq. @T).
        Δh = gas.enthalpy_mass - h_initial  # [J/kg]

        # Get the mass fraction of hydrogen produced.
        y_H2 = gas["H2"].Y[0]  # mass fraction of H2

        # Compute the energy cost.
        energy_cost = Δh / y_H2 if y_H2 > 0 else np.inf  # [J/kgH2]

        # Convert the energy cost to kWh/kgH2.
        energy_cost_kWh_per_kgH2 = energy_cost / 3.6e6  # [kWh/kgH2]

        if heat_recovery:
            # Set the temperature and pressure for the reference state (298.15 K and 1 atm).
            gas.TP = T_ref, pressure_ref

            # Compute the enthalpy change of heat_recovery down the products from T to 298.15 K.
            h_cooling = h_T_eq - gas.enthalpy_mass  # [J/kg]

            # Subtract the heat_recovery enthalpy from the energy cost.
            energy_cost_kWh_per_kgH2 -= h_cooling / (y_H2 * 3.6e6)  # [kWh/kgH2]

        energy_costs.append(energy_cost_kWh_per_kgH2)

    cases[case_name].energy_costs = energy_costs

Computing energy cost for case 1a...
Computing energy cost for case 1b...
Computing energy cost for case 2a...
Computing energy cost for case 2b...

Plot the energy cost as a function of temperature.

fig, ax = plt.subplots()
for case_name, case_params in cases.items():
    with_carbon = case_params.with_carbon
    heat_recovery = case_params.heat_recovery
    energy_costs = case_params.energy_costs
    label = case_params.label
    xy = case_params.xy
    x, y = xy

    ax.plot(temperatures - 273.15, energy_costs)
    ax.text(
        x,
        y,
        label,
        fontsize=30,
        color=ax.get_lines()[-1].get_color(),
        ha="center",
        va="center",
        bbox=dict(facecolor="white", alpha=0.7, edgecolor="none"),
    )

ax.set_xlabel("Temperature (°C)")
ax.set_ylabel("Energy cost (kWh/kgH₂)")
ax.set_title("Energy cost of producing H₂ from CH₄", fontsize=32)
# Set x and y axis limits for better visualization.
ax.set_ylim(0, 40)
ax.set_xlim(0, 2500)
plt.show()