Note

Go to the end to download the full example code.

Isobaric chemical relaxation: non-equilibrium n_e collapsing to LTE.#

A CH4 plasma column starts already at thermal equilibrium (Te = Tg = T0, fixed P0) but with an electron density well above its chemical equilibrium value at (T0, P0) – a pure ionization/recombination relaxation, with no electric field, no circuit, and no other external driver. Isobaric + adiabatic finite-rate chemistry (IsobaricPlasma0D1T) integrates the recombination; the released energy is exactly the specific enthalpy that must be conserved. Once the composition stops changing (a Damkohler-number guard, damkohler()) the ladder hands off, smoothly, to IsobaricLTE – a stage that just holds the equilibrium state, demonstrating the switch is loss-less: the two models agree at the hand-off instant.

This exercises the same two-tier economy already validated for the isochoric deposition tier (IsochoricPlasma0D2T -> IsochoricLTE), for the isobaric cooling tier instead, built directly from the adaptive-model primitives (AdaptiveCompositeReactor) rather than the circuit-coupled production ladder (PulsedPlasmaReactor).

Tags: plasma Cantera CH4 adaptive isobaric LTE recombination

Imports.#

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

import rizer.kinetics.extensible_rate  # noqa: F401 (Register CH4 custom rates)
from rizer.adaptive_models.composite import AdaptiveCompositeReactor
from rizer.adaptive_models.models_list import IsobaricLTE, IsobaricPlasma0D1T
from rizer.adaptive_models.selector import SwitchRule, damkohler, electron_density
from rizer.adaptive_models.state import PlasmaState
from rizer.adaptive_models.transitions import IsobaricToThermochemicalEquilibrium
from rizer.misc.plt_utils import set_mpl_style
from rizer.misc.utils import get_path_to_data
from rizer.transport.loaders import (
    get_default_collision_frequency_model,
    get_momentum_transfer_collision_frequencies_list,
)
from rizer.transport.mixture_law import MixtureCollisionFrequencies

set_mpl_style()

Discharge setup.#

T0/P0/bump are chosen so thermal ionization is small-but-visible at (T0, P0) and the bumped n_e is dramatic on a log plot; not re-tuned against a reference measurement.

mechanism = str(get_path_to_data("mechanisms", "Goutier2025", "CH4_to_C2H2.yaml"))
T0 = 3000.0  # Initial (and ambient) temperature [K]
P0 = ct.one_atm  # [Pa]
BUMP = 1.0e4  # n_e seeded at BUMP x its equilibrium value at (T0, P0)
SEED_ION = "CH4+"

gap = 3.8e-3  # [m]
radius = 500.0e-6  # [m]

plasma = ct.Solution(mechanism, "plasma", transport_model=None)

# True equilibrium composition at (T0, P0): the baseline n_e this example
# perturbs away from. A fresh Solution's default composition is pure "e-"
# (no C/H at all) -- seed a real CH4 mixture first, since `equilibrate`
# conserves elemental composition exactly.
# `TPX` before `Te` here (not the usual `Te`-first order elsewhere in this
# codebase, needed whenever the mixture carries a real electron fraction: the
# two-temperature EOS's density solve, done at `TPX`-set time, depends on
# whichever `Te` is current then) is harmless in both cases below -- checked
# by direct comparison against the `Te`-first order, zero difference in the
# resulting P/density: `X_e = 0` here (pure "CH4:1.0", no ions seeded yet),
# so the EOS doesn't consult `Te` at all.
plasma.TPX = T0, P0, "CH4:1.0"
plasma.Te = T0
plasma.equilibrate("TP")  # Leaves `Te` at T0 -- equilibrate("TP") doesn't touch it.
i_e = plasma.species_index("e-")
i_ion = plasma.species_index(SEED_ION)
x_e_eq = float(plasma.X[i_e])

X = np.array(plasma.X, dtype=float)
X[i_e] = BUMP * x_e_eq
X[i_ion] = BUMP * x_e_eq
X /= X.sum()
# `X_e` is real now (`BUMP` x its equilibrium value), but `Te` is already T0
# from above (unchanged by `equilibrate`), so it isn't "stale" here either.
plasma.TPX = T0, P0, X
plasma.Te = T0

Y0 = plasma.Y.copy()
rho0 = float(plasma.density)
mass = rho0 * gap * np.pi * radius**2

s0 = PlasmaState(
    t=0.0,
    Y=Y0,
    mechanism=mechanism,
    Tg=T0,
    Te=T0,
    Tv=T0,
    rho=rho0,
    V=mass / rho0,
    P=P0,
    R=radius,
    mass=mass,
    gap=gap,
)

mtcf = get_default_collision_frequency_model()
mtcf_list = get_momentum_transfer_collision_frequencies_list(plasma.species_names, mtcf)
collision_freq = MixtureCollisionFrequencies(plasma, mtcf_list)

Build the two-stage ladder and its single switch rule.#

No circuit, no heat loss (T_wall=None – adiabatic): the only physics is isobaric finite-rate recombination releasing its stored energy as heat.

solver_input = {
    "max_order": 5,
    "max_step": 0.0,
    "atol": 1.0e-16,
    "rtol": 1.0e-10,
}
tol = 1.0e-3

stage_kin = IsobaricPlasma0D1T(plasma, mtcf_list, solver_input, T_wall=None)
stage_lte = IsobaricLTE(plasma)

TAU = 1.0e-6  # [s] Damkohler probe timescale -- tuned against the plotted n_e(t) knee
THRESHOLD = 0.05  # composition considered "stopped changing" below this Da


def _chemistry_converged(state: PlasmaState, t: float) -> float:
    """Decreasing indicator: composition has stopped changing (Da -> 0)."""
    return damkohler(state, plasma, tau=TAU)


rule = SwitchRule(
    outgoing=stage_kin.name,
    incoming=stage_lte.name,
    indicator=_chemistry_converged,
    threshold=THRESHOLD,
    operators=(IsobaricToThermochemicalEquilibrium(atol_rel=tol),),
)


class RecombinationRelaxation(AdaptiveCompositeReactor):
    """Minimal composite: two stages, one rule, no circuit."""

    def __init__(self, plasma, collision_freq, s0, stages, rules):
        self._s0 = s0
        super().__init__(plasma, collision_freq, stages, rules)

    def _initial_state(self) -> PlasmaState:
        return self._s0.copy()

    def _macro_grid(self, t_end: float, stage1: object) -> np.ndarray:
        """Geometric grid sized to this example's own recombination timescale."""
        t0 = float(self._s0.t)
        return np.geomspace(t0 + 1.0e-13, t_end, 600)


T_END = 5.0e-5  # [s]
reactor = RecombinationRelaxation(
    plasma, collision_freq, s0, [stage_kin, stage_lte], [rule]
)
result = reactor.solve(t_end=T_END)

/home/runner/work/rizer/rizer/rizer/adaptive_models/composite.py:187: UnreviewedPhysicsWarning: model selector 'ModelSelector' (ModelSelector) uses physics that no human has reviewed: see rizer/adaptive_models/PHYSICS.md for its equations and known limitations, and set ModelSelector.reviewed_by once validated.
  self._selector = ModelSelector(
/home/runner/work/rizer/rizer/examples/plasma/plot_isobaric_recombination_relaxation.py:177: UnreviewedPhysicsWarning: model 'isobaric_plasma_0d1t' (IsobaricPlasma0D1T) uses physics that no human has reviewed: see rizer/adaptive_models/PHYSICS.md for its equations and known limitations, and set IsobaricPlasma0D1T.reviewed_by once validated.
  result = reactor.solve(t_end=T_END)
/home/runner/work/rizer/rizer/rizer/adaptive_models/selector.py:452: UnreviewedPhysicsWarning: transition operator 'IsobaricToThermochemicalEquilibrium' (IsobaricToThermochemicalEquilibrium) uses physics that no human has reviewed: see rizer/adaptive_models/PHYSICS.md for its equations and known limitations, and set IsobaricToThermochemicalEquilibrium.reviewed_by once validated.
  seed = op.apply(seed, self._plasma)
/home/runner/work/rizer/rizer/rizer/adaptive_models/selector.py:578: UnreviewedPhysicsWarning: model 'isobaric_lte' (IsobaricLTE) uses physics that no human has reviewed: see rizer/adaptive_models/PHYSICS.md for its equations and known limitations, and set IsobaricLTE.reviewed_by once validated.
  self._arm_if_ready(state, t)

Plot.#

MW = np.asarray(plasma.molecular_weights, dtype=float)
n_e = np.array(
    [
        electron_density(s0.copy(Y=Y_row, rho=mass / V_row), i_e, MW)
        for Y_row, V_row in zip(result.Y, result.V)
    ]
)  # [m^-3]

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(9, 9))
ax1.plot(result.t, result.Tg, label="T [K]")
ax1.set_ylabel("Temperature [K]")
ax1.legend(loc="best")
ax2.semilogy(result.t, n_e, label="n_e [m^-3]")
ax2.set_xlabel("t [s]")
ax2.set_ylabel("Electron density [m^-3]")
ax2.legend(loc="best")
for t_sw, name in zip(result.switch_times, result.active_model_log[1:]):
    for ax in (ax1, ax2):
        ax.axvline(t_sw, color="k", ls=":", lw=0.8)
    ax2.text(t_sw, n_e.max(), name, rotation=90, va="top", ha="right", fontsize=8)
plt.show()
plot isobaric recombination relaxation

Sanity checks.#

T rises only slightly (a fraction of a Kelvin): the excess electrons recombining are still a tiny fraction of the gas’s total particle density even at BUMP=1e4, so their recombination energy is a tiny fraction of the gas’s total thermal energy. The dramatic signal here is n_e collapsing several orders of magnitude, not T.

assert result.active_model_log == [stage_kin.name, stage_lte.name]
assert len(result.switch_times) == 1
np.testing.assert_allclose(result.Y.sum(axis=1), 1.0, atol=1.0e-9)
assert np.all(np.diff(n_e[result.t > result.switch_times[0]]) <= 0.0), (
    "n_e must be non-increasing after the switch (recombination, not re-ionization)"
)
assert result.Tg[-1] > T0, "recombination heat release must raise T above T0"

plasma.TP = result.Tg[-1], P0
plasma.equilibrate("TP")
n_e_eq_final = electron_density(
    PlasmaState(
        t=0.0,
        Y=np.array(plasma.Y, dtype=float),
        mechanism=mechanism,
        Tg=result.Tg[-1],
        Te=result.Tg[-1],
        rho=float(plasma.density),
        V=1.0,
        P=P0,
        R=radius,
        mass=mass,
        gap=gap,
    ),
    i_e,
    MW,
)
assert abs(n_e[-1] / n_e_eq_final - 1.0) < 0.05, (
    f"final n_e ({n_e[-1]:.3e}) should match equilibrium at the final T "
    f"({n_e_eq_final:.3e}) to within a few percent"
)
print("All checks passed.")

All checks passed.

Total running time of the script: (0 minutes 6.363 seconds)