rizer.models.nrp.plasma_channel_2T1D#

Transient 1D two-temperature reacting plasma channel (Cantera-based).

PlasmaChannel is the 1-D radial, transient counterpart of rizer’s 0-D rizer.models.nrp.isomass_2T_volume_reactor.Isomass2TVolumeReactor: it solves finite-rate two-temperature plasma chemistry together with radial heat conduction, and exposes the discharge-channel diameter expansion D(t) for comparison with rizer’s engineering prediction (rizer.models.nrp.engineering_model.diffusion) and experiment.

The heavy solve lives in the compiled rizer.cantera_ext._plasma1d; this wrapper builds the initial state, calls it, and post-processes the radial temperature history into a channel radius.

Notes

The native-rate air_plasma_Laux2000 mechanism loads in C++; the Goutier CH4 mechanism has Python-only ExtensibleRate reactions and needs C++ ports for a fully-kinetic CH4 run (see ARCHITECTURE.md). For a thermal (LTE-conduction) CH4 expansion, run with reacting=False and an injected CH4 kappa_H(T).

Classes#

PlasmaChannel

Transient radial two-temperature reacting plasma channel.

Module Contents#

class rizer.models.nrp.plasma_channel_2T1D.PlasmaChannel(mech: str, phase: str, R_max: float, rho: float, Y0: numpy.ndarray | list[float], Y_amb: numpy.ndarray | list[float], *, grid: numpy.ndarray | list[float] | None = None, Tg_profile: tuple[numpy.ndarray, numpy.ndarray] | None = None, Te_profile: tuple[numpy.ndarray, numpy.ndarray] | None = None, Y0_profile: tuple[numpy.ndarray, numpy.ndarray] | None = None, T_amb: float = 300.0, Tg0: float | None = None, Te0: float | None = None, kappa_H: tuple[numpy.ndarray, numpy.ndarray] | float | None = None, kappa_e: tuple[numpy.ndarray, numpy.ndarray] | float | None = None, D_species: float = 0.0, hydrodynamics: bool = False, rho_profile: tuple[numpy.ndarray, numpy.ndarray] | None = None, u_r_profile: tuple[numpy.ndarray, numpy.ndarray] | None = None, outer_bc_mode: str = 'ambient_ghost', p_amb: float = -1.0, cfl: float = 0.9, sigma: tuple[numpy.ndarray, numpy.ndarray] | float = 0.0, nu_E: float = 0.0, nu_m: float | tuple[numpy.ndarray, numpy.ndarray] = 0.0, mtcf: list[rizer.transport.collision_frequency.MomentumTransferCollisionFrequencyModel] | None = None, Te_min: float = 200.0, Te_max: float = 100000.0, Te_n: int = 400, electric_field: float | tuple[numpy.ndarray, numpy.ndarray] = 0.0, reacting: bool = True, dt: float = 1e-08, n_steps: int = 100, record_every: int = 1, n_points: int = 81, loglevel: int = 0, gap: float | None = None, rtol: float = 1e-07, atol_T: float = 0.001, atol_Y: float = 1e-09, independent_nodes: bool = False, n_threads: int = 1, plasma_voltage: collections.abc.Callable[[float, float], float] | None = None)#

Transient radial two-temperature reacting plasma channel.

Assumptions and domain of validity

This model is a 1-D radial, axisymmetric generalisation of rizer’s 0-D Isomass2TVolumeReactor. Read these before interpreting results:

Density and flow depend on hydrodynamics. Off (the default), there is no continuity or momentum equation: rho is a fixed parameter, the state per node is (Tg, Te, Y_s) and there is no flow. Pressure is then not held constant: a hot core at fixed rho sits at an elevated local pressure, and the recorded n_e is the model’s own value at that (T, rho) state. That mode is valid when the channel evolves by conduction/diffusion and chemistry, and not valid for pressure-driven (shock / Braginskii) gas-dynamic expansion.

On, the state gains (rho, u_r) and a radial hydrodynamic sub-step advances them in a Strang split around the reaction step, which is what represents that expansion. rho then sets the initial density and rho_profile a radial one; p_amb must be given, since it has no default that is safe to assume.
Prescribed electric field E(t); no self-consistent field. Poisson is not solved and the plasma is treated quasi-neutral, so the dielectric-relaxation timescale never enters. Valid for a given/measured applied field; not valid for self-consistent sheath/streamer dynamics (use a drift-diffusion-Poisson solver for those).
Two-temperature (Te != Tg): separate electron and heavy-species energy equations coupled by elastic/inelastic exchange, assuming a Maxwellian EEDF at Te.
Heat conduction: two-temperature Fourier — radial conduction for both Tg (kappa_H(Tg)) and Te (kappa_e(Te)), each active only if its conductivity table is supplied.
Species transport: Fickian with a single constant scalar D_species. No thermal diffusion (Soret) or Dufour, not mixture-averaged / multicomponent (one D for every species), and no ambipolar diffusion for charged species (electrons and ions share D, no field coupling). Species diffusion is frozen when D_species == 0 (the default). Species diffusion also feeds back into both energy equations: a compressibility term (from the two-temperature ideal-gas EOS) and an enthalpy-diffusion (species heat transport) term, each nonzero whenever a composition gradient coincides with species of different heat capacity – not just the mass-fraction equations (see ARCHITECTURE.md).
Joule heating sigma*E^2. Two mutually exclusive closures for sigma (and, for the electron energy equation, the elastic electron-heavy exchange frequency): (a) mtcf, a composition-resolved collision model evaluated live from the evolving composition every step (the same model the 0-D Isomass2TVolumeReactor uses) – preferred, since sigma then tracks the reacting mixture instead of a frozen snapshot; or (b) the legacy closure: sigma from a sigma(Te) table, or – if nu_m > 0 – sigma = n_e e^2/(m_e nu_m) from the local n_e, plus a separately-supplied nu_E. Passing both raises ValueError.
Geometry/BC: cylindrical 1-D radial, axisymmetric; r=0 symmetry; r=R_max Dirichlet wall at T_amb (and ambient composition Y_amb for diffusion).
Numerics: a stiff adaptive-order BDF integrator (Cantera CVODE over the method-of-lines RHS), error-controlled in both order and step, which is what handles the field-driven plasma chemistry. The stiff limit is the electron-energy equation (tiny electron heat capacity -> sub-ns Te dynamics): start gently (Te=Tg) and use a modest field, as rizer’s 0D NRP reactor does. Both air and the CH4 Goutier2025 mechanism run this way.

Parameters:

mech (str) – Cantera mechanism file and phase name (a PlasmaPhase with e-).
phase (str) – Cantera mechanism file and phase name (a PlasmaPhase with e-).
R_max (float) – Outer radius of the computational domain [m] (ignored if grid is given).
rho (float) – (Fixed) mass density [kg/m^3].
Y0 (array_like) – Initial composition inside the domain where no Y0_profile is supplied (length = number of species). The electron mass fraction must be positive at every node, here and in Y0_profile: the electron temperature is undefined where there are no electrons, so an electron-free start raises. Seed a small fraction with a matching positive ion. Y_amb may hold no electrons.
Y_amb (array_like) –
Ambient mass fractions OUTSIDE the domain (length = number of species). Used for the outer-wall ghost of the diffusion terms, the wall-ghost density, and – with hydrodynamics – the composition of gas flowing inward across the wall face.

Required, and never derived from Y0. The two are independent physical quantities: a seeded Y0 left to stand in for the ambient is a plasma flowing inward at the boundary, which on the reference CH4 discharge puts the wall node at n_e 2.3e20 instead of 2.9e14 and moves the plasma resistance by up to a factor of two, since the outermost cell carries the largest area weight in integral(sigma dA). Pass Y_amb=Y0 only when the boundary genuinely holds the initial composition.
grid (array_like or None) – Explicit, increasing radial grid [m] (grid[0] == 0). Overrides n_points/R_max; installed directly (Cantera-flame style). If None, a uniform n_points grid on [0, R_max] is built.
Tg_profile (tuple(array, array) or None) – (r, T) initial radial temperature profiles [m, K]. If None, uniform T_amb.
Te_profile (tuple(array, array) or None) – (r, T) initial radial temperature profiles [m, K]. If None, uniform T_amb.
Y0_profile (tuple(array, 2d-array) or None) – (r, Y) initial radial composition, Y shaped (len(r), nsp). If None, the uniform Y0 is used everywhere.
T_amb (float) – Wall/ambient temperature [K]: the outer-wall Dirichlet ghost value for both Tg and Te, and the default for Tg0/Te0.
Tg0 (float or None) –
Initial gas and electron temperature [K], uniform where no Tg_profile/Te_profile is supplied. None (the default) starts them at T_amb.

They are the initial state, not the boundary. A channel that starts hot behind a cold wall sets these above T_amb; one that starts at ambient leaves them None. Giving T_amb both jobs is what makes the two easy to confuse, which is why they are separate arguments here. In C++ they are separate arguments of solve_1D_unsteady_discharge rather than fields of PlasmaChannelConfig: the binding turns each into a flat two-point initTg_r/initTg_v (or initTe_*) profile, and only T_amb reaches the config as a temperature of its own.
Te0 (float or None) –
Initial gas and electron temperature [K], uniform where no Tg_profile/Te_profile is supplied. None (the default) starts them at T_amb.

They are the initial state, not the boundary. A channel that starts hot behind a cold wall sets these above T_amb; one that starts at ambient leaves them None. Giving T_amb both jobs is what makes the two easy to confuse, which is why they are separate arguments here. In C++ they are separate arguments of solve_1D_unsteady_discharge rather than fields of PlasmaChannelConfig: the binding turns each into a flat two-point initTg_r/initTg_v (or initTe_*) profile, and only T_amb reaches the config as a temperature of its own.
kappa_H (tuple(array, array) or float or None) – Heavy-pool / electron thermal conductivity, as (T, value) tables or a constant [W/m/K]. None disables that conduction.
kappa_e (tuple(array, array) or float or None) – Heavy-pool / electron thermal conductivity, as (T, value) tables or a constant [W/m/K]. None disables that conduction.
D_species (float) – Constant species mass diffusivity [m^2/s] (0 = frozen transport).
sigma (tuple(array, array) or float) – Electrical conductivity sigma(Te) [S/m] for Joule heating (legacy closure; mutually exclusive with mtcf).
nu_E (float) – Effective electron-heavy elastic energy-exchange frequency [1/s] (legacy closure; mutually exclusive with mtcf).
nu_m (float or tuple(array, array)) – Electron momentum-transfer frequency [1/s]; if > 0, sigma = n_e e^2/(m_e nu_m) from the local n_e (overrides the sigma table). Legacy closure; mutually exclusive with mtcf.
mtcf (list of MomentumTransferCollisionFrequencyModel or None) – Per-species momentum-transfer collision-frequency models (as built by rizer.transport.loaders.get_momentum_transfer_collision_frequencies_list()), one per species in mechanism order. When given, sigma and the electron-heavy elastic exchange are evaluated live from the evolving composition every step via the same composition-resolved collision model Isomass2TVolumeReactor uses, instead of the frozen sigma/nu_m/nu_E closure below. Mutually exclusive with sigma, nu_m, and nu_E: passing a non-default value of any of those together with mtcf raises ValueError.
Te_min (float, float, int) – Log-spaced electron-temperature grid [K],[K],[-] the collision model (mtcf) precomputes its per-species cross-section tables on. Only used when mtcf is given.
Te_max (float, float, int) – Log-spaced electron-temperature grid [K],[K],[-] the collision model (mtcf) precomputes its per-species cross-section tables on. Only used when mtcf is given.
Te_n (float, float, int) – Log-spaced electron-temperature grid [K],[K],[-] the collision model (mtcf) precomputes its per-species cross-section tables on. Only used when mtcf is given.
electric_field (float or tuple(array, array)) – Constant field E [V/m], or an arbitrary schedule (t, E) [s, V/m] (linearly interpolated, e.g. a pulse) supplied from Python.
reacting (bool) – Include finite-rate chemistry.
dt (float, int, int) – Output cadence / initial step [s] (dt*n_steps = total time), number of output steps, recording cadence.
n_steps (float, int, int) – Output cadence / initial step [s] (dt*n_steps = total time), number of output steps, recording cadence.
record_every (float, int, int) – Output cadence / initial step [s] (dt*n_steps = total time), number of output steps, recording cadence.
n_points (int) – Radial grid points (ignored if grid is given).
rtol (float) – Error tolerances: relative [-], a temperature floor [K], and a mass-fraction floor [-]. Loosening them buys much less than the step count suggests – bigger steps are harder steps, so the per-step cost rises as the count falls and the two nearly cancel.
atol_T (float) – Error tolerances: relative [-], a temperature floor [K], and a mass-fraction floor [-]. Loosening them buys much less than the step count suggests – bigger steps are harder steps, so the per-step cost rises as the count falls and the two nearly cancel.
atol_Y (float) – Error tolerances: relative [-], a temperature floor [K], and a mass-fraction floor [-]. Loosening them buys much less than the step count suggests – bigger steps are harder steps, so the per-step cost rises as the count falls and the two nearly cancel.
independent_nodes (bool) –
False (the default) solves the reaction sub-step as one banded system over the whole domain; True gives every node its own integrator and step size. Legal only with every transport term off, and rejected otherwise.

Worth 2.19x / 3.41x / 4.57x at 15 / 41 / 81 nodes on a discharge with a hot core and a cold edge, and SLOWER on a uniform relaxation: the gain is in how unlike the nodes are, not their number.
n_threads (int) – Threads advancing the per-node systems. Needs independent_nodes=True; each thread costs one mechanism parse before the first step. Measured 2.57x at 41 nodes on four threads, and bit-identical at every count.
plasma_voltage (typing.Callable or None) –
(t [s], R_p [Ohm]) -> V_p [V], closing an electrical circuit on the channel’s own radially-integrated conductance, R_p = gap / integral(sigma dA). Needs gap. None (the default) leaves the field prescribed by electric_field, which is a field the plasma cannot push back on.

R_p is frozen over each global step and V_p interpolated across it, so the callable sees two monotonically increasing times per step and never the integrator’s internal trial times – which is what a stateful circuit, such as a transmission line keeping a wave history, requires. Adds Rp, Rp_state, Vp, I and E to the recorded history.
hydrodynamics (bool) – Run the radial hydrodynamic sub-step in a Strang split around the reaction step. Off by default, and off means rho is a fixed parameter and there is no flow. On, p_amb must be given, record_every must divide n_steps, and the grid needs at least 2 points.
rho_profile (tuple(array, array) or None) – Initial radial density (r, rho) [m],[kg/m^3], every value > 0. None holds the scalar rho at every node.
u_r_profile (tuple(array, array) or None) – Initial radial velocity (r, u_r) [m],[m/s]. None starts at rest. Only read when hydrodynamics is set.
outer_bc_mode (str) –
"ambient_ghost" (default, and the accurate choice while R_max sits beyond the front) or "pressure_outlet", which pins p at the wall and is measurably reflecting. Only read when hydrodynamics is set, but validated either way.

"pressure_outlet" is a documented non-default, kept as the comparison point that justifies the default rather than as a treatment to reach for: validation case 21 measures it returning 33% to 62% of an incident wave against ambient_ghost’s 1.3% to 2.4%, and it rings instead of settling. Removing it would delete the only evidence in the repo that the default is the right one. No example or driver selects it.
p_amb (float) – Ambient pressure the outer boundary imposes [Pa]. No default that silently validates: a zero is a vacuum, so hydrodynamics requires it to be set.
cfl (float) – Courant number for the hydro sub-step’s own CFL sub-cycling [-], in (0, 1].

T_amb#

property t: numpy.ndarray#: Recorded times [s].

property r: numpy.ndarray#: Radial grid [m].

property species: list[str]#: Species names, in the mechanism order used by mass_fractions().

temperature_profiles() → tuple[numpy.ndarray, numpy.ndarray]#

Return the recorded gas- and electron-temperature histories.

Returns:

Tg (numpy.ndarray) – Gas (heavy-species) temperature [K], shape (nframes, npts).
Te (numpy.ndarray) – Electron temperature [K], shape (nframes, npts).

electron_density() → numpy.ndarray#

Electron number density n_e(r, t), shape (nframes, npts) [1/m^3].

Computed inside the solver as X_e * P / (k_B * Tmean) per node (the same n_e used in the electron energy balance), so it is the model’s own value rather than a post-hoc reconstruction.

mass_fractions() → numpy.ndarray#: Species mass fractions Y_s(r, t), shape (nframes, npts, nsp) [-].

density_profiles() → numpy.ndarray#

Mass density rho(r, t), shape (nframes, npts) [kg/m^3].

Constant at the configured rho unless hydrodynamics is on, in which case the hydro sub-step solves it.

velocity_profiles() → numpy.ndarray#

Radial velocity u_r(r, t), shape (nframes, npts) [m/s].

Zero throughout unless hydrodynamics is on.

pressure_profiles() → numpy.ndarray#

Total pressure p(r, t), shape (nframes, npts) [Pa].

Derived from the two-temperature equation of state rather than held as a state variable, and computed inside the solver so it comes from the same state-setting the source terms use.

field() → numpy.ndarray#

Return the applied electric field E(t) [V/m] at the recorded output times.

With plasma_voltage set, this is the circuit’s field as the solver recorded it, V_p / gap at the end of each step (0 at t = 0). Otherwise it reconstructs the prescribed field the C++ solver evaluated: a constant, or linear interpolation of the (t, E) schedule supplied at construction (see ChannelOptions::Et_t/ Et_v in PlasmaChannelSolver.h).

Returns:: Field [V/m] at each recorded time in t.
Return type:: numpy.ndarray

total_current() → numpy.ndarray#

Total current I(t) [A] through the channel cross-section.

I = E(t) \int_0^R \sigma(r) 2\pi r\,dr. \sigma follows the same fallback as the solver’s own Joule-heating term: \sigma = n_e e^2/(m_e \nu_m) (from the local n_e) when nu_m > 0 or a nu_m(Te) table is given; otherwise the sigma(Te) conductivity table supplied at construction (matching PlasmaCell::transport, so Joule heating configured purely via sigma is still diagnosable). Assumes the recorded radial profile is uniform along the inter-electrode gap.

The integral is the cell sum of _cell_areas, which is the quadrature the plasma resistance uses, so the current reported here and the R_p a circuit closes on describe the same channel.

With plasma_voltage set, this is the current the circuit drove and the solver recorded, V_p / R_p with R_p frozen at the start of each step, for either closure. Rebuilding it from each frame’s own conductivity would divide by the end-of-step resistance instead, which across a breakdown differs from the one that carried the current.

Raises:: NotImplementedError – If the channel was built with mtcf (the composition-resolved collision model) and no circuit: reconstructing sigma(r, t) after the fact would need the live collision-model evaluation at each recorded state, which isn’t exposed from Python yet – reconstructing it from a frozen table here would silently misrepresent the solve.

joule_power(gap: float | None = None) → numpy.ndarray#

Total Joule power dissipated in the channel cross-section.

P(t) = E(t) * gap * I(t) = U(t) * I(t) [W], i.e. the 1-D radial current profile from total_current() extruded uniformly over the inter-electrode gap (the same constant-along-the-gap assumption as total_current()).

Parameters:: gap (float or None, optional) – Inter-electrode gap [m]. If None (default), uses the gap passed at construction.
Returns:: Joule power [W] at each recorded time in t.
Return type:: numpy.ndarray
Raises:: ValueError – If no gap was supplied here or at construction.

deposited_energy(gap: float | None = None) → float#

Total deposited (Joule) energy over the run.

E_dep = \int_0^{t_{end}} E(t) * gap * I(t)\,dt [J], the time integral of joule_power() (trapezoidal rule over the recorded output times).

Parameters:: gap (float or None, optional) – Inter-electrode gap [m]. If None (default), uses the gap passed at construction.
Returns:: Total deposited energy [J].
Return type:: float

channel_radius(method: str = 'half_max') → numpy.ndarray#

Channel radius R(t) [m] from the gas-temperature profile.

Parameters:

method ({"half_max", "sigma_weighted"}) –

half_max: radius where Tg - T_amb drops to half its axis value.
sigma_weighted: sqrt(<r^2>) weighted by Tg - T_amb.

Returns:

Channel radius [m] at each recorded time in t.

Return type:

numpy.ndarray

channel_diameter(method: str = 'half_max') → numpy.ndarray#

Channel diameter D(t) = 2 R(t) [m].

See channel_radius() for the method options; this is the quantity directly comparable to rizer’s engineering diffusion model (rizer.models.nrp.engineering_model.diffusion) and experiment.

Parameters:: method ({"half_max", "sigma_weighted"}) – Passed through to channel_radius().
Returns:: Channel diameter [m] at each recorded time in t.
Return type:: numpy.ndarray