Note

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

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

—

Field-driven ionization then relaxation of a 1D plasma channel (n_e video)#

A transient 1D radial run of PlasmaChannel showing an initially ionized channel that ionizes further under an applied electric field and then relaxes once the field switches off, rendering a video of the electron number density n_e(r, t) along the radius.

What this exercises#

Every piece added for this scenario:

  • an explicit non-uniform radial grid supplied from Python (Cantera-flame style, grid= -> Domain1D::setupGrid), clustered near the axis;

  • an arbitrary field schedule E(t) passed from Python as a (t, E) table (here a trapezoidal pulse: ramp up, hold, ramp down) – any waveform works;

  • Joule heating with sigma = n_e e^2 / (m_e nu_m) from the local n_e;

  • the stiff BDF backend (integrator="bdf", Cantera’s CVODE over the method-of-lines RHS), which integrates the field-driven plasma chemistry that the fixed-step Newton solver cannot;

  • the solver-recorded n_e(r, t) (electron_density()).

Mechanism note#

This runs on the air mechanism air_plasma_Laux2000 (native two-temperature- plasma rates), chosen here because air ionizes strongly at the few-thousand-K temperatures used, giving a clear n_e swing. The CH4 Goutier2025 mechanism also runs in this 1D channel (its electron-impact rates are finite up to >=3 eV; use integrator="bdf" with a gentle Te=Tg start and a modest field, mirroring rizer’s 0D NRP reactor) – but CH4 only ionizes weakly below ~1 eV, so a visible field-driven n_e rise needs a higher electron temperature than is comfortable for the solver here. Air makes the ionization-then-relaxation physics easiest to see; the pipeline is identical for either mechanism.

Tags: plasma Cantera two-temperature ionization transient channel 1D CVODE

Import the required libraries.#

import cantera as ct
import numpy as np

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

set_mpl_style()

Scenario setup.#

mechanism = str(get_path_to_data("mechanisms") / "air_plasma_Laux2000.yaml")
phase = "plasma"

# -- Scenario parameters ------------------------------------------------ #
wall_radius = 4.0e-3  # Wall radius [m]
wall_temperature = 2000.0  # Wall/ambient temperature [K]
initial_channel_width = 500.0e-6  # Initial channel width (FWHM) [m]
core_temperature = 7000.0  # Initial gas/electron temperature on axis [K]

# -- Field pulse parameters ---------------------------------------------- #
field_peak = 2.0e5  # Peak applied field [V/m]
t_rise = 30.0e-9  # Field pulse rise time [s]
t_hold = 90.0e-9  # Field pulse hold time [s]
t_fall = 40.0e-9  # Field pulse fall time [s]
t_end = 250.0e-9  # Total simulated time [s]

# -- Collision frequency model -------------------------------------------- #
nu_m = 1.0e12  # e-neutral momentum-transfer frequency [1/s]

# -- Solver parameters ----------------------------------------------------- #
n_points = 31  # Number of radial grid points [-]
dt = 1.0e-9  # Output cadence / initial BDF step [s]

# -- Derived quantities ---------------------------------------------------- #
sigma_gaussian = initial_channel_width / (
    2.0 * np.sqrt(2.0 * np.log(2.0))
)  # Gaussian std dev of the initial channel [m]

Define the helper functions.#

The first builds a non-uniform radial grid, the second builds the initial hot channel in local LTE, the third defines the trapezoidal field pulse, and the last renders the n_e(r, t) video from the solved channel.

def stretched_grid(R, n, beta=3.5):
    """Non-uniform radial grid on [0, R], fine near r=0 (sinh stretching)."""
    u = np.linspace(0.0, 1.0, n)
    return R * np.sinh(beta * u) / np.sinh(beta)


def build_initial_state(r):
    """Hot Gaussian channel in local LTE (consistent initial ionization)."""
    g = np.exp(-(r**2) / (2.0 * sigma_gaussian**2))
    Tg = wall_temperature + (core_temperature - wall_temperature) * g
    gas = ct.Solution(mechanism, phase)
    gas.TPX = wall_temperature, ct.one_atm, "N2:0.79, O2:0.21"
    gas.equilibrate("TP")
    Y0 = gas.Y.copy()  # ambient/wall composition
    gas.TPX = (
        0.5 * (core_temperature + wall_temperature),
        ct.one_atm,
        "N2:0.79, O2:0.21",
    )
    rho = float(gas.density)  # one (constant-volume) density
    Y_2d = np.zeros((len(r), gas.n_species))
    for i in range(len(r)):
        gas.TPX = Tg[i], ct.one_atm, "N2:0.79, O2:0.21"
        gas.equilibrate("TP")
        Y_2d[i] = gas.Y
    return rho, Y0, Tg, Y_2d


def field_pulse():
    """Trapezoidal E(t): 0 ->(rise) field_peak ->(hold) field_peak ->(fall) 0 -> 0."""
    t1, t2, t3 = t_rise, t_rise + t_hold, t_rise + t_hold + t_fall
    return (
        np.array([0.0, t1, t2, t3, t_end]),
        np.array([0.0, field_peak, field_peak, 0.0, 0.0]),
    )


def make_video(ch, path="channel_ne.mp4"):
    """Animate n_e(r, t); save mp4 (ffmpeg) -> gif (pillow) -> skip. Returns path or None."""
    import matplotlib.pyplot as plt
    from matplotlib.animation import FFMpegWriter, FuncAnimation, PillowWriter

    r_mm = ch.r * 1e3
    ne = ch.electron_density() * 1e-6  # cm^-3
    t = ch.t
    t_off = t_rise + t_hold + t_fall

    fig, ax = plt.subplots(figsize=(7, 4.5))
    (line,) = ax.semilogy([], [], "o-", color="crimson", ms=3)
    ax.set_xlim(0, r_mm.max())
    ax.set_ylim(ne.min() * 0.7, ne.max() * 1.5)
    ax.set_xlabel("Radius [mm]")
    ax.set_ylabel(r"$n_e$ [cm$^{-3}$]")
    title = ax.set_title("")

    def update(f):
        line.set_data(r_mm, ne[f])
        phase = "field ON" if t[f] <= t_off else "field OFF (relaxing)"
        title.set_text(f"Air channel  t = {t[f] * 1e9:5.0f} ns   ({phase})")
        return line, title

    anim = FuncAnimation(fig, update, frames=len(t), interval=120, blit=False)
    saved = None
    for writer, p in (
        (FFMpegWriter(fps=8), path),
        (PillowWriter(fps=8), path.replace(".mp4", ".gif")),
    ):
        try:
            anim.save(p, writer=writer)
            saved = p
            break
        except Exception as exc:  # noqa: BLE001
            print(f"  ({type(writer).__name__} unavailable: {exc})")
    plt.close(fig)
    return saved

Run the transient 1D channel solve.#

print(
    f"Air channel: E = {field_peak / 1e3:.0f} kV/m pulse "
    f"({t_rise * 1e9:.0f}/{t_hold * 1e9:.0f}/{t_fall * 1e9:.0f} ns), total {t_end * 1e9:.0f} ns"
)

r = stretched_grid(wall_radius, n_points)
rho, Y0, Tg, Y_2d = build_initial_state(r)
nu_E = (2.0 * 9.1093837015e-31 / (29.0e-3 / ct.avogadro)) * nu_m
n_steps = int(round(t_end / dt))
channel = PlasmaChannel(
    mechanism,
    phase,
    wall_radius,
    rho,
    Y0,
    grid=r,
    Tg_profile=(r, Tg),
    Te_profile=(r, Tg),
    Y0_profile=(r, Y_2d),
    T_amb=wall_temperature,
    nu_m=nu_m,
    nu_E=nu_E,
    electric_field=field_pulse(),  # arbitrary E(t) table from Python
    kappa=2.0,
    kappa_e=0.1,
    reacting=True,
    integrator="bdf",  # stiff CVODE backend
    dt=dt,
    n_steps=n_steps,
    record_every=max(n_steps // 40, 1),
)

Render the electron density video.#

saved = make_video(channel)
if saved:
    print(f"saved n_e animation -> {saved}")

Plot the static summary figure.#

import matplotlib.pyplot as plt

r_mm = channel.r * 1e3
idx = np.unique(np.linspace(0, len(channel.t) - 1, 6).astype(int))
fig, (axn, axt) = plt.subplots(1, 2, figsize=(11, 4.2))
for f in idx:
    axn.semilogy(
        r_mm, channel.electron_density()[f] * 1e-6, label=f"{channel.t[f] * 1e9:.0f} ns"
    )
axn.set_title("Electron density n_e(r) (snapshots)")
axn.set_xlabel("Radius [mm]")
axn.set_ylabel(r"$n_e$ [cm$^{-3}$]")
axn.legend(fontsize="small")
axt.plot(channel.t * 1e9, ne[:, 0] * 1e-6, "o-", color="crimson", label=r"$n_e$ axis")
axt.set_xlabel("Time [ns]")
axt.set_ylabel(r"$n_e$ axis [cm$^{-3}$]")
axt.axvspan(
    0,
    (t_rise + t_hold + t_fall) * 1e9,
    color="gold",
    alpha=0.25,
    label="field on",
)
axt.set_title("Axis electron density vs time")
axt.legend(fontsize="small")
fig.tight_layout()
plt.show()