Note

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

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

—

1D vs 0D on a self-consistent NRP circuit (each closes its own loop)#

The companion plot_1d_vs_0d_nrp.py drives the 1D channel with the 0D’s prescribed field, which lets the channel ionize its surroundings and broaden without limit (the integrated current/energy come out ~7x the 0D). That is the radial physics a 0D misses, but it is exaggerated because the field is never re-loaded by the broadening channel.

Here both models close the same NRP circuit on their own plasma resistance, so the comparison is quantitatively trustworthy:

  • 0D: R_p = gap/(sigma * pi*radius^2) -> circuit -> V_p -> E = V_p/gap.

  • 1D: R_p = gap / integral(sigma * 2*pi*r dr) (the radially-integrated conductance) -> the same circuit model -> V_p -> E. As the 1D channel broadens, its R_p falls, which loads the line and pulls V_p (hence the field and the Joule heating) down – self-limiting the broadening.

Both use the native C++ per-node reactor (composition-resolved sigma); the 1D is the operator-split solver (solve_split_circuit()). With the loop closed, the two currents and deposited energies are of the same order (the residual difference is the genuine 1D effect: a somewhat wider conducting channel at a correspondingly lower R_p), instead of the order-of-magnitude gap seen under a prescribed field.

Tags: plasma Cantera CH4 two-temperature transient channel 0D 1D NRP self-consistent

Import the required libraries.#

import time

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

import rizer.kin.extensible_rate  # noqa: F401 (Register CH4 custom rates)
import rizer.misc.units as u
from rizer.electric_circuit.cable import IdealCable
from rizer.electric_circuit.generator import TrapezoidalGenerator
from rizer.electric_circuit.nrp_circuit import NRPCircuit
from rizer.misc.plt_utils import set_mpl_style
from rizer.misc.simulation.simulation import (
    get_default_collision_frequency_model,
    get_momentum_transfer_collision_frequencies_list,
)
from rizer.misc.utils import get_path_to_data
from rizer.plasma.constant_mass_reactor_cpp import ConstantMassPlasmaReactorOdeCpp
from rizer.plasma.plasma_channel_split import ChannelReactor1D

set_mpl_style()

Discharge setup.#

mechanism = str(get_path_to_data("mechanisms") / "Goutier2025" / "CH4_to_C2H2.yaml")

# -- Geometry parameters ----------------------------------------------- #
gap = 3.8e-3  # Constant gap [m]
radius = 500e-6  # Plasma radius (constant, since gamma = np.inf) [m]
r_wall = 4.0e-3  # 1D domain outer radius [m]

# -- Initial conditions ------------------------------------------------ #
P0 = ct.one_atm  # Initial pressure [Pa]
Tg_0 = 1000.0  # Initial gas temperature [K]
ne_0 = 1.0e19  # Initial electron density [m^-3]

# -- Electric parameters ----------------------------------------------- #
# Number of reflections to consider
number_of_reflections = 1  # [-]

# Generator parameters (trapezoidal pulse)
U_on = 2000.0  # Voltage when on [V]
t_rise = 5e-9  # Rising time [s]
t_on = 25e-9  # On time [s]
t_fall = 15e-9  # Falling time [s]

# -- Solver parameters --------------------------------------------------#
t_end = 80e-9  # End time [s]
dt_out = 1e-9  # Output time step [s]
n_radial_points = 15  # Number of radial points in the 1D channel [-]


def make_circuit():
    """Build the NRP circuit (trapezoidal generator + ideal cable)."""
    return NRPCircuit(
        generator=TrapezoidalGenerator(
            R_g=1.0, U_on=U_on, U_off=0.0, t_rise=t_rise, t_on=t_on, t_fall=t_fall
        ),
        cable=IdealCable(L=6.2, Z_c=75.0, c=1.9e8),
    )


# Load the default per-species collision-frequency model.
cfm = get_default_collision_frequency_model()

plasma = ct.Solution(mechanism, "plasma", transport_model=None)
names = plasma.species_names
mtcf = get_momentum_transfer_collision_frequencies_list(names, cfm)
x_e = ne_0 / (P0 / (u.k_b * Tg_0))
plasma.Te = Tg_0
plasma.TPX = Tg_0, P0, f"CH4:{1 - 2 * x_e:.6e}, e-:{x_e:.6e}, CH4+:{x_e:.6e}"
rho_0 = plasma.density
V0 = gap * np.pi * radius**2

Solve the 0D reactor, closing its own circuit.#

ode = ConstantMassPlasmaReactorOdeCpp(
    mechanism,
    "plasma",
    mtcf,
    mass=V0 * rho_0,
    gap=gap,
    electric_circuit=make_circuit(),
    nb_reflections=number_of_reflections,
    polytropic_index=np.inf,
    p_ext=P0,
    species_names=names,
)
sol = scipy.integrate.ode(ode)
sol.set_integrator(
    "vode",
    method="bdf",
    order=5,
    with_jacobian=True,
    first_step=1e-13,
    max_step=1e-10,
    atol=1e-13,
    rtol=1e-7,
    nsteps=100000,
)
sol.set_initial_value(np.hstack((Tg_0, Tg_0, V0, plasma.Y.copy())), 0.0)
T0_lists = {
    k: [v] for k, v in dict(t=0.0, Te=Tg_0, ne=ne_0, Rp=0.0, I=0.0, Vp=0.0).items()
}
print("running 0D (own circuit) ...")
t_start = time.time()
for ti in np.arange(1e-10, t_end + 1e-12, 1e-10):
    sol.integrate(ti)
    if not sol.successful():
        print(f"  0D stopped at {sol.t:.2e}s")
        break
    ne = ode._r0d.electron_density(
        sol.y[0], sol.y[1], np.ascontiguousarray(sol.y[3:]), ode.mass / sol.y[2]
    )
    Rp, Vp = ode.plasma_resistance, ode.plasma_voltage
    for k, v in dict(
        t=sol.t, Te=sol.y[1], ne=ne, Rp=Rp, Vp=Vp, I=(Vp / Rp if Rp > 0 else 0.0)
    ).items():
        T0_lists[k].append(v)
T0 = {k: np.array(v) for k, v in T0_lists.items()}
Edep0 = np.concatenate(
    (
        [0.0],
        np.cumsum(
            0.5
            * (T0["Vp"][1:] * T0["I"][1:] + T0["Vp"][:-1] * T0["I"][:-1])
            * np.diff(T0["t"])
        ),
    )
)
print(
    f"  0D: {time.time() - t_start:.1f}s, peak |I|={np.max(np.abs(T0['I'])):.2f} A, "
    f"E_dep={Edep0[-1] * 1e3:.3f} mJ"
)

Solve the 1D channel, closing its own circuit.#

r = r_wall * np.sinh(3.5 * np.linspace(0, 1, n_radial_points)) / np.sinh(3.5)
gp = np.exp(-(r**2) / (2 * (radius / np.sqrt(2 * np.log(2))) ** 2))
ne1 = ne_0 * gp + 1e12
g = ct.Solution(mechanism, "plasma")
Y0 = np.zeros((n_radial_points, plasma.n_species))
for i in range(n_radial_points):
    xe = ne1[i] * u.k_b * Tg_0 / P0
    g.TPX = Tg_0, P0, f"CH4:{1 - 2 * xe:.6e}, e-:{xe:.6e}, CH4+:{xe:.6e}"
    Y0[i] = g.Y
t_out = np.arange(0.0, t_end + 1e-12, dt_out)
ch = ChannelReactor1D(
    mechanism, "plasma", r, rho_0, mtcf, names, reacting=True, kappa=0.5, T_amb=Tg_0
)
print("running 1D (own circuit, operator-split) ...")
t_start = time.time()
h = ch.solve_split_circuit(
    t_out,
    make_circuit(),
    gap,
    number_of_reflections,
    np.full(n_radial_points, Tg_0),
    np.full(n_radial_points, Tg_0),
    Y0,
)
t1 = np.array(h["t"])
Te1 = np.array([v[0] for v in h["Te"]])
ne1a = np.array([v[0] for v in h["ne"]])
I1 = np.array(h["I"])
Rp1 = np.array(h["Rp"])
Vp1 = np.array(h["Vp"])
Edep1 = np.concatenate(
    ([0.0], np.cumsum(0.5 * (Vp1[1:] * I1[1:] + Vp1[:-1] * I1[:-1]) * np.diff(t1)))
)
print(
    f"  1D: {time.time() - t_start:.1f}s, peak |I|={np.max(np.abs(I1)):.2f} A, "
    f"E_dep={Edep1[-1] * 1e3:.3f} mJ"
)

Radially-integrated 1D quantities, to compare with the lumped 0D.#

The 1D Te/n_e vary with radius; we show both the on-axis (r=0) value AND a radial average. For Te the meaningful average is electron-density-weighted (the bulk electron temperature), <Te> = int n_e Te 2pi r dr / int n_e 2pi r dr. For n_e the domain volume-average <n_e> = int n_e 2pi r dr / (pi R_wall^2). I, R_p, V_p, E_dep are NOT axis quantities – they are the channel’s global (radially-integrated / circuit) values, e.g. I = E int sigma 2pi r dr.

r = ch.r
two_pi_r = 2.0 * np.pi * r
Te_prof = np.array(h["Te"])
ne_prof = np.array(h["ne"])
ne_line = np.trapezoid(ne_prof * two_pi_r, r, axis=1)  # int n_e 2pi r dr [1/m]
Te_neavg = np.trapezoid(ne_prof * Te_prof * two_pi_r, r, axis=1) / np.maximum(
    ne_line, 1e-300
)
ne_volavg = ne_line / (np.pi * r[-1] ** 2)  # domain volume-average

Plot the results and compare the agreement between the two implementations.#

def two(a, x0, y0, x1, y1, ylab, title, lab1="1D channel", log=False):
    """Plot a 0D vs 1D comparison panel on the given axes."""
    a.plot(x0 * 1e9, y0, "-", lw=2, label="0D reactor")
    a.plot(x1 * 1e9, y1, "--", lw=2, label=lab1)
    a.set_xlabel("Time [ns]")
    a.set_ylabel(ylab)
    a.set_title(title)
    if log:
        a.set_yscale("log")
    a.legend(fontsize="small")


fig, ax = plt.subplots(2, 3, figsize=(15, 8))
fig.suptitle(
    "CH4 NRP: 0D vs 1D, each closing its OWN NRP circuit (self-consistent field)"
)

# Te: on-axis vs the n_e-weighted radial mean (bulk electron temperature).
a = ax[0, 0]
a.plot(T0["t"] * 1e9, T0["Te"], "-", lw=2, label="0D reactor")
a.plot(t1 * 1e9, Te1, "--", lw=2, label="1D axis (r=0)")
a.plot(t1 * 1e9, Te_neavg, ":", lw=2, label=r"1D $\langle T_e\rangle_{n_e}$ (radial)")
a.set_xlabel("Time [ns]")
a.set_ylabel("Te [K]")
a.set_title("Electron temperature")
a.legend(fontsize="small")

# n_e: on-axis vs the radial volume-average over the channel domain.
a = ax[0, 1]
a.plot(T0["t"] * 1e9, T0["ne"], "-", lw=2, label="0D reactor")
a.plot(t1 * 1e9, ne1a, "--", lw=2, label="1D axis (r=0)")
a.plot(
    t1 * 1e9,
    ne_volavg,
    ":",
    lw=2,
    label=r"1D $\langle n_e\rangle$ (radial vol. avg)",
)
a.set_yscale("log")
a.set_xlabel("Time [ns]")
a.set_ylabel("n_e [m^-3]")
a.set_title("Electron density")
a.legend(fontsize="small")

# Global (radially-integrated / circuit) quantities -- "1D channel", not "axis".
# The current and resistance legends carry the radial integral that defines them.
two(ax[0, 2], T0["t"], T0["Vp"], t1, Vp1, "V_p [V]", "Plasma voltage")
two(
    ax[1, 0],
    T0["t"],
    np.abs(T0["I"]),
    t1,
    np.abs(I1),
    "|I| [A]",
    "Total current",
    lab1=r"1D channel: $I=E\!\int\!\sigma\,2\pi r\,dr$",
)
two(
    ax[1, 1],
    T0["t"],
    T0["Rp"],
    t1,
    Rp1,
    "R_p [Ohm]",
    "Plasma resistance",
    lab1=r"1D channel: $R_p=\mathrm{gap}/\!\int\!\sigma\,2\pi r\,dr$",
    log=True,
)
two(
    ax[1, 2],
    T0["t"],
    Edep0 * 1e3,
    t1,
    Edep1 * 1e3,
    "E_dep [mJ]",
    "Deposited energy",
)
fig.tight_layout(rect=(0, 0, 1, 0.96))
plt.show()