Note

You can download below the full example code. and run it online in Codespaces.

https://github.com/codespaces/badge.svg

—

Computing and comparing the properties of an \(CH_4\) plasma.#

Compute the properties of an \(CH_4\) plasma using minplascalc. Data are compared with reference data from the literature.

Properties calculated include:

  • Density

  • Enthalpy

  • Heat capacity

  • Viscosity

  • Thermal conductivity

  • Electrical conductivity

Species data is loaded from ./data/minplascalc/species folder. These data have been generated using the generate_minplascalc_data.py script.

Reference data are [Wu2016] and [Niu2016].

Import the required libraries.#

import cantera as ct
import matplotlib.pyplot as plt
import minplascalc as mpc
import numpy as np
from adjustText import adjust_text

import rizer.misc.units as u
from rizer.io.thermo_transport_data_reader import ThermoTransportDataReader
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,
)
from rizer.misc.utils import get_path_to_data

set_mpl_style(nb_columns=1)

Create mixture object for the species we’re interested in.#

Next, we create a minplascalc LTE mixture object. Here we use a helper function in minplascalc which creates the object directly from a list of the species names.

species_and_concentrations = {
    "CH4": 1,
    "CH3": 0,
    "CH2": 0,
    "CH": 0,
    "C": 0,
    "C+": 0,
    "C2H6": 0,
    "C2H4": 0,
    "C2H2": 0,
    "C2H": 0,
    "C2": 0,
    "H2": 0,
    "H": 0,
    "H+": 0,
}

SPECIES_PATH = get_path_to_data("minplascalc", "species")
species_list = [
    mpc.species.from_file(SPECIES_PATH / f"{species_name}.json")
    for species_name in species_and_concentrations
]

methane_mixture = mpc.mixture.LTE(
    species=species_list,
    x0=species_and_concentrations.values(),
    T=300,
    P=101325,
    gfe_initial_particles=1e23,
    gfe_rtol=1e-10,
    gfe_max_iter=1000,
)

Set a range of temperatures to calculate the equilibrium compositions at.#

Next, set a range of temperatures to calculate the equilibrium compositions at - in this case we’re going from 300 to 6000 K in 100 K steps. Also initialise a list to store the property values at each temperature.

temperatures = np.arange(300, 6000, 100)
species_names = [sp.name for sp in methane_mixture.species]
ni_list: list[np.ndarray] = []
density: list[float] = []
h: list[float] = []
cp: list[float] = []
viscosity: list[float] = []
thermal_conductivity: list[float] = []
electrical_conductivity: list[float] = []

Perform the composition calculations.#

Now we can perform the property calculations.

Note that execution of this calculation is fairly compute intensive and the following code snippet may take several minute to complete.

for T in temperatures:
    print(f"Calculating properties at T = {T} K...")
    methane_mixture.T = T
    # give back the composition of the species in the mixture?
    ni_list.append(methane_mixture.calculate_composition())

    density.append(methane_mixture.calculate_density())
    h.append(methane_mixture.calculate_enthalpy())
    cp.append(methane_mixture.calculate_heat_capacity())
    viscosity.append(methane_mixture.calculate_viscosity())
    thermal_conductivity.append(methane_mixture.calculate_thermal_conductivity())
    electrical_conductivity.append(methane_mixture.calculate_electrical_conductivity())

ni = np.array(ni_list).transpose()

print("Done!")

Load reference data.#

# Load reference data from Wu2016 and Niu2016 for comparison.
data_ref = ThermoTransportDataReader(
    gas_name="CH4", pressure_atm=1, source="Wu2016_Niu2016", ignore_missing_values=True
)
temperatures_ref = data_ref.temperature
enthalpy_ref = data_ref.enthalpy
thermal_conductivity_ref = data_ref.thermal_conductivity

Plot the results.#

Now we can visualise the properties by plotting them against temperature, to see how they vary.

# Plot enthalpy vs. temperature.
fig, ax = plt.subplots()
ax.plot(temperatures, h, "k", label="minplascalc")
ax.plot(temperatures_ref, enthalpy_ref, "k--", label="Wu2016")
ax.legend()
ax.set_title(r"$\mathregular{CH_4}$ plasma enthalpy at 1 atm")
ax.set_xlabel("T [K]")
ax.set_ylabel(r"$\mathregular{H [J/kg]}$")
ax.set_xlim(300, 6000)
ax.set_ylim(-0.1e8, 2e8)
plt.show()

# Plot thermal conductivity vs. temperature.
fig, ax = plt.subplots()
ax.plot(temperatures, thermal_conductivity, "k", label="minplascalc")
ax.plot(temperatures_ref, thermal_conductivity_ref, "k--", label="Wu2016")
ax.legend()
ax.set_title(r"$\mathregular{CH_4}$ plasma thermal conductivity at 1 atm")
ax.set_xlabel("T [K]")
ax.set_ylabel(r"$\mathregular{kappa [W/(m.K)]}$")
ax.set_xlim(300, 6000)

plt.show()

Compute the equilibrium composition with Cantera.#

gas = load_goutier2025_mechanism("gas")
states = ct.SolutionArray(gas, shape=(len(temperatures),))
states.TPX = temperatures, ct.one_atm, {"CH4": 1}
states.equilibrate("TP")
xi_cantera = states.X.transpose()
species_names_cantera = states.species_names

# ni_cantera = xi_cantera * n_total
# P = n_total * k_b * T
n_total = ct.one_atm / (u.k_b * temperatures)
ni_cantera = xi_cantera * n_total

Plot the species composition.#

fig, ax = plt.subplots()
ax.set_title(r"$\mathregular{CH_4}$ LTE plasma composition with temperature at 1 atm")
ax.set_xlabel("T [K]")
ax.set_ylabel(r"$\mathregular{n_i [m^{-3}]}$")
# ax.set_ylim(1e15, 5e25)
ax.set_ylim(1e-6, 100)
ax.set_xlim(300, 6000)


texts = []
for species_density, species_name in zip(ni, species_names):
    if species_name not in species_names_cantera:
        if species_name == "e":
            species_name = "e-"
        else:
            continue

    n_tot = ct.one_atm / (u.k_b * temperatures) / 100

    ax.semilogy(
        temperatures, species_density / n_tot, color=get_species_color(species_name)
    )
    texts.append(
        get_text(
            x=temperatures[np.argmax(species_density)],
            y=np.max(species_density / n_tot),
            text=get_species_in_latex(species_name),
            ax=ax,
            color=get_species_color(species_name),
        )
    )
    idx_ct = species_names_cantera.index(species_name)
    ax.scatter(
        temperatures[::5],
        ni_cantera[idx_ct, ::5] / n_tot[::5],
        color=get_species_color(species_name),
        marker="x",
        s=100,
    )

adjust_text(texts, avoid_self=False)
plt.show()