Note
Go to the end to download the full example code.
Does the radial expansion change the discharge? 7 kV into CH4 through a 6.2 m cable#
A nanosecond pulse heats a thin channel at nearly constant volume, so the core ends up over-pressured and pushes outward. The question this example answers is whether that motion feeds back on the electrical behaviour, or whether it is something that merely happens afterwards.
It is answered by running the same discharge through the same circuit
– 7 kV trapezoid, 1 ohm source, 6.2 m of 75 ohm cable with reflections, 3.8 mm
of CH4 at one atmosphere, 300 ns – under four models: the 1-D radial channel
with \(\rho\) and \(u_r\) solved (hydrodynamics=True), the same
channel with the density frozen, and a 0-D lumped reactor at constant volume
and on a polytropic expansion. Everything else is identical, so each difference
belongs to one modelling choice.
Within the 1-D channel, the answer has two halves, and only the first is the obvious one:
The first current pulse does not care. Peak current is 184.089 A with and without the sub-step, and the resistance minimum is 0.04937 vs 0.04942 ohm, a part in 2000. The channel breaks down in about 12 ns; the gas needs hundreds of nanoseconds to move a kernel radius. The expansion has not happened yet while the current is being drawn.
The cable echoes do care. The 6.2 m line returns a reflection every \(2L/c = 65.3\) ns, so the third arrives at 196 ns, by which time the channel has rarefied substantially on axis. The current it drives peaks a few nanoseconds later at 118.9 A into the expanded channel against 125.9 A into the frozen one, and over the 300 ns window the expanded channel takes 3850 uJ against 4090 uJ, 5.9 % less.
So the expansion is real – 1461 m/s, 20.6 atm on axis at peak, 56.5 % of the axis density gone by 300 ns – and it is invisible to the pulse that caused it. What it changes is what the line sees afterwards. Without a long cable there would be nothing for it to change.
Against the 0-D reactor, two more things separate:
Geometry dominates the resistance. The lumped cell holds one conductivity over a fixed 0.60 mm disc and bottoms out at 6.16 ohm; the 1-D channel integrates \(\int \sigma \, 2 \pi r \, \mathrm{d}r\) over a core that ionizes far harder than the mean, and reaches 0.049 ohm. That is a factor of 125, and it is the radial resolution, not the flow.
The current barely notices, because the line sets it. Peak current is 170.2 A in 0-D against 184.1 A in 1-D, 8 % apart despite those two decades of resistance: once the channel conducts, 75 ohms of cable is the larger impedance and the plasma stops being what limits the current.
A lumped expansion is not a small version of a radial one. The polytropic cell grows by 15 % in volume over the window, a 13 % density drop, against 56.5 % on the 1-D axis. A volume law driven by the cell’s mean pressure against ambient cannot represent a rarefaction running inward from the core, which is the mechanism that empties the axis.
Conduction and species diffusion are off (kappa = 0), so what is plotted is
the pulse, the chemistry and the hydrodynamics alone. The 0-D radius is a
choice rather than a measurement, and the resistance comparison moves with it;
the derivation used is stated where it is set.
Import the required libraries.#
import pathlib
import tempfile
import time
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
import rizer.kinetics.extensible_rate # noqa: F401 (registers the CH4 rates)
import rizer.misc.units as u
from rizer.electrical_model.circuit.transmission_line_circuit import (
TransmissionLineResistiveCircuit,
)
from rizer.electrical_model.components.cable import IdealCable
from rizer.electrical_model.components.generator import TrapezoidalGenerator
from rizer.electrical_model.components.source_circuit import ResistiveSourceCircuit
from rizer.misc.ct_utils import partially_ionized_mole_fractions
from rizer.misc.plt_utils import set_mpl_style
from rizer.misc.utils import get_path_to_data
from rizer.models.nrp.isomass_2T_volume_reactor_cpp import Isomass2TVolumeReactor
from rizer.models.nrp.plasma_channel_split_2T1D import ChannelReactor1D
from rizer.transport.loaders import get_default_momentum_transfer_collision_frequencies
set_mpl_style()
Discharge, circuit and grid.#
The applied field is not U_on / gap: the circuit sets the plasma voltage
against the channel’s own resistance
\(R_p = \mathrm{gap} / \int \sigma \,\mathrm{d}A\), so a more conductive
channel loads the line and lowers its own field. The cable length is what
makes this example about more than the first 12 ns: at c = 1.9e8 m/s a
6.2 m line brings each reflection back after 2 L / c.
mechanism = str(get_path_to_data("mechanisms") / "Goutier2025" / "CH4_to_C2H2.yaml")
gap = 3.8e-3 # Inter-electrode gap [m]
kernel_radius = 500e-6 # Seed ionization width (FWHM) [m]
r_wall = 4.0e-3 # Domain outer radius [m]
p_0 = ct.one_atm # Initial pressure [Pa]
Tg_0 = 1000.0 # Initial gas and electron temperature [K]
ne_0 = 1.0e19 # Peak seed electron density [m^-3]
U_on = 7000.0 # Generator voltage when on [V]
t_rise, t_on, t_fall = 3e-9, 6e-9, 3e-9 # Trapezoid timings [s]
cable_length, cable_speed = 6.2, 1.9e8 # [m], [m/s]
t_end = 300e-9 # End time [s]
dt_out = 1e-9 # Output cadence [s]
n_radial_points = 31 # Radial nodes [-]
# The 0D reactor lumps the channel into one cell and so needs a single radius,
# which a Gaussian seed does not uniquely supply. The top-hat radius carrying
# the same electron inventory at the peak density is used: for a Gaussian of
# standard deviation `sigma`, int(n_e * 2 pi r dr) = n_e0 * 2 pi sigma^2, which
# a disc of radius sigma*sqrt(2) at n_e0 matches. Any other defensible choice
# moves the 0D resistance, which is itself part of what the comparison shows.
seed_sigma = kernel_radius / np.sqrt(2.0 * np.log(2.0)) # [m]
radius_0d = seed_sigma * np.sqrt(2.0) # [m]
#: Heat-capacity ratio of the seeded CH4 at the initial state [-], used as the
#: 0D polytropic index for its isentropic case. The native reactor takes a
#: number rather than the symbolic "gamma", so it is held fixed for the run
#: while the real ratio falls to about 1.07 by 5000 K.
polytropic_index_0d = 1.1272
round_trip = 2.0 * cable_length / cable_speed # Cable echo period [s]
print(f"cable round trip 2L/c = {round_trip * 1e9:.1f} ns")
def make_circuit():
"""Build the NRP circuit (trapezoidal generator + ideal cable).
Returns
-------
TransmissionLineResistiveCircuit
The 7 kV source behind 1 ohm, feeding 6.2 m of 75 ohm line with
reflections enabled.
"""
return TransmissionLineResistiveCircuit(
source=ResistiveSourceCircuit(
TrapezoidalGenerator(
U_on=U_on, U_off=0.0, t_rise=t_rise, t_on=t_on, t_fall=t_fall
),
R_g=1.0,
),
cable=IdealCable(L=cable_length, Z_c=75.0, c=cable_speed),
include_reflections=True,
)
One seed, two runs.#
Both channels start from the identical Gaussian ionized core on the identical
stretched grid, so nothing but hydrodynamics separates them. Y_amb is
the ambient composition the outer ghost cell is held at; it is required once
the sub-step is on, because the species fluxes only telescope into the mass
flux at the wall when it sums to one.
plasma = ct.Solution(mechanism, "plasma", transport_model=None)
names = list(plasma.species_names)
mtcf = get_default_momentum_transfer_collision_frequencies(names)
r = r_wall * np.sinh(3.5 * np.linspace(0.0, 1.0, n_radial_points)) / np.sinh(3.5)
width = kernel_radius / np.sqrt(2.0 * np.log(2.0))
ne_profile = ne_0 * np.exp(-(r**2) / (2.0 * width**2)) + 1.0e12
gas = ct.Solution(mechanism, "plasma")
Y0 = np.zeros((n_radial_points, plasma.n_species))
for node in range(n_radial_points):
x_e = ne_profile[node] * u.k_b * Tg_0 / p_0
gas.TPX = Tg_0, p_0, partially_ionized_mole_fractions(x_e)
Y0[node] = gas.Y
gas.TPX = Tg_0, p_0, "CH4:1"
Y_ambient = gas.Y.copy()
x_e = ne_0 * u.k_b * Tg_0 / p_0
plasma.TPX = Tg_0, p_0, partially_ionized_mole_fractions(x_e)
rho_0 = float(plasma.density)
t_out = np.arange(0.0, t_end + 1e-12, dt_out)
def solve(hydrodynamics):
"""Run the circuit-coupled channel with the radial sub-step on or off.
Parameters
----------
hydrodynamics : bool
Whether to solve continuity, radial momentum, species advection and the
two pool energies alongside the chemistry.
Returns
-------
dict
The recorded history, as :meth:`ChannelReactor1D.solve_split_circuit`
returns it.
"""
channel = ChannelReactor1D(
mechanism,
"plasma",
r,
rho_0,
mtcf,
names,
reacting=True,
kappa=0.0,
T_amb=Tg_0,
p_ext=p_0,
hydrodynamics=hydrodynamics,
Y_amb=Y_ambient,
)
start = time.time()
history = channel.solve_split_circuit(
t_out,
make_circuit(),
gap,
np.full(n_radial_points, Tg_0),
np.full(n_radial_points, Tg_0),
Y0,
)
label = "expanding" if hydrodynamics else "frozen"
print(f" {label:<10} done in {time.time() - start:.1f} s")
return history
print("running the same pulse twice ...")
expanding = solve(hydrodynamics=True)
frozen = solve(hydrodynamics=False)
t_ns = np.asarray(expanding["t"]) * 1e9
r_mm = r * 1e3
Tg = np.asarray(expanding["Tg"])
Te = np.asarray(expanding["Te"])
electron_density = np.asarray(expanding["ne"])
pressure = np.asarray(expanding["p"])
velocity = np.asarray(expanding["u_r"])
density = np.asarray(expanding["rho"])
def deposited(history):
"""Electrical energy delivered to the plasma over the window [J]."""
return float(
np.trapezoid(np.asarray(history["I"]) * np.asarray(history["Vp"]), t_out)
)
print(
f" peak |I| expanding {np.max(np.abs(expanding['I'])):.3f} A, "
f"frozen {np.max(np.abs(frozen['I'])):.3f} A"
)
print(
f" min R_p expanding {min(expanding['Rp']):.5f} ohm, "
f"frozen {min(frozen['Rp']):.5f} ohm"
)
print(
f" energy expanding {deposited(expanding) * 1e6:.1f} uJ, "
f"frozen {deposited(frozen) * 1e6:.1f} uJ "
f"({100 * (deposited(expanding) / deposited(frozen) - 1):+.1f} %)"
)
print(
f" peak Te {Te.max() * u.k_b / u.e:.2f} eV, peak n_e "
f"{electron_density.max():.3e} m^-3, peak |u_r| {np.abs(velocity).max():.0f} m/s"
)
print(
f" axis density {density[0, 0]:.4f} -> {density[-1, 0]:.4f} kg/m^3 "
f"({100 * (density[-1, 0] / density[0, 0] - 1):+.1f} %)"
)
The same circuit seen by a 0D reactor.#
Isomass2TVolumeReactor
lumps the channel into one cell of fixed mass and closes the identical
circuit on \(R_p = \mathrm{gap} / (\sigma \pi R^2)\). It is run twice as
well: at constant volume (polytropic_index = inf) and letting the cell
expand against ambient on a polytropic law. That second run is the lumped
analogue of the radial sub-step, which is what makes it worth plotting
beside it.
def solve_0d(polytropic_index):
"""Run the 0D reactor on the same circuit, at a given polytropic index.
Parameters
----------
polytropic_index : float
``numpy.inf`` holds the volume fixed; a finite value lets the cell
expand against ``p_ext`` on :math:`p V^n = \\mathrm{const}`.
Returns
-------
dict of str to numpy.ndarray
``t`` [s], ``Tg``/``Te`` [K], ``ne`` [m^-3], ``Rp`` [ohm], ``I`` [A],
``Vp`` [V] and the cell volume ``V`` [m^3], on a 0.1 ns grid.
"""
cell = ct.Solution(mechanism, "plasma")
cell.TPX = Tg_0, p_0, partially_ionized_mole_fractions(x_e)
volume = gap * np.pi * radius_0d**2
reactor = Isomass2TVolumeReactor(
cell,
mechanism,
"plasma",
mtcf,
mass=volume * float(cell.density),
gap=gap,
electric_circuit=make_circuit(),
polytropic_index=polytropic_index,
p_ext=p_0,
)
net = ct.ReactorNet([reactor])
net.max_order = 5
net.max_time_step = 1e-10
net.atol, net.rtol, net.max_steps = 1e-13, 1e-7, 100000
columns = dict(t=0.0, Tg=Tg_0, Te=Tg_0, ne=ne_0, Rp=0.0, I=0.0, Vp=0.0, V=volume)
out = {key: [value] for key, value in columns.items()}
start = time.time()
for target in np.arange(1e-10, t_end + 1e-12, 1e-10):
net.advance(target)
state = np.asarray(reactor.get_state())
electrons = reactor._r0d.electron_density(
state[0],
state[1],
np.ascontiguousarray(state[3:]),
reactor.plasma_mass / state[2],
)
resistance, voltage = reactor.plasma_resistance, reactor.plasma_voltage
for key, value in dict(
t=net.time,
Tg=state[0],
Te=state[1],
ne=electrons,
Rp=resistance,
Vp=voltage,
I=(voltage / resistance if resistance > 0 else 0.0),
V=state[2],
).items():
out[key].append(value)
label = "constant volume" if np.isinf(polytropic_index) else "polytropic"
print(f" 0D {label:<16} done in {time.time() - start:.1f} s")
return {key: np.asarray(value) for key, value in out.items()}
zero_d_frozen = solve_0d(np.inf)
zero_d_expanding = solve_0d(polytropic_index_0d)
for name, run in (
("0D constant volume", zero_d_frozen),
("0D polytropic", zero_d_expanding),
):
energy = float(np.trapezoid(run["I"] * run["Vp"], run["t"]))
live = run["Rp"] > 0.0
print(
f" {name:<20} peak |I| {np.max(np.abs(run['I'])):.2f} A, "
f"min R_p {run['Rp'][live].min():.4f} ohm, "
f"energy {energy * 1e6:.1f} uJ, V/V0 {run['V'][-1] / run['V'][0]:.2f}"
)
The current and the resistance: identical, then not.#
The first pulse lands on top of itself. Afterwards the two channels drift
apart, and the cable echoes arriving at multiples of 2L/c fall on
measurably different plasmas. The dashed vertical lines are those echo
arrivals, drawn from the cable length rather than read off the trace.
fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))
echoes = np.arange(round_trip, t_end, round_trip) * 1e9
for ax, (label, key, scale) in zip(
axes,
(
("Current I [A]", "I", 1.0),
(r"$R_p$ [$\Omega$]", "Rp", 1.0),
("Plasma voltage [kV]", "Vp", 1e-3),
),
strict=True,
):
ax.plot(t_ns, np.asarray(expanding[key]) * scale, label="1D expanding", lw=1.4)
ax.plot(t_ns, np.asarray(frozen[key]) * scale, "--", label="1D frozen", lw=1.2)
ax.plot(
zero_d_expanding["t"] * 1e9,
zero_d_expanding[key] * scale,
":",
color="0.35",
lw=1.2,
label="0D polytropic",
)
for echo in echoes:
ax.axvline(echo, color="0.8", lw=0.6, zorder=0)
ax.set_xlabel("t [ns]")
ax.set_ylabel(label)
axes[1].set_yscale("log")
axes[0].legend(fontsize="small")
fig.suptitle("Same circuit, same seed: the expansion only shows up after the pulse")
Zoom on the breakdown: the expansion is not there yet.#
Over the first 20 ns the two traces are indistinguishable. This is the half of the answer that says the radial motion cannot be inferred from the first current pulse.
early = t_ns <= 20.0
fig, axes = plt.subplots(1, 2, figsize=(10, 3.6))
axes[0].plot(t_ns[early], np.asarray(expanding["I"])[early], lw=1.6, label="expanding")
axes[0].plot(t_ns[early], np.asarray(frozen["I"])[early], "--", lw=1.2, label="frozen")
axes[0].set_ylabel("Current I [A]")
axes[0].legend(fontsize="small")
relative = (np.asarray(expanding["I"])[early] - np.asarray(frozen["I"])[early]) / max(
np.max(np.abs(np.asarray(frozen["I"])[early])), 1e-30
)
axes[1].plot(t_ns[early], 100 * relative, color="firebrick", lw=1.4)
axes[1].set_ylabel("Current difference [% of peak]")
for ax in axes:
ax.set_xlabel("t [ns]")
fig.suptitle("First 20 ns: breakdown finishes before the gas moves")
Radial profiles through the pulse and the expansion.#
The five fields the channel actually carries, at six times spanning the pulse, the relaxation and the rarefaction. The pressure panel is where the mechanism is clearest: a core at 20 atm at 10 ns, flat by 500 ns, with the velocity panel showing where the gas went.
snapshot_times = [0.0, 10e-9, 50e-9, 150e-9, 300e-9, t_end]
snapshots = sorted(
{min(int(round(t / dt_out)), len(t_out) - 1) for t in snapshot_times}
)
fields = (
(r"$T_\mathrm{H}$ [K]", Tg, None),
(r"$T_\mathrm{e}$ [K]", Te, None),
(r"$n_\mathrm{e}$ [m$^{-3}$]", electron_density, "log"),
(r"$p$ [atm]", pressure / ct.one_atm, None),
(r"$u_r$ [m s$^{-1}$]", velocity, None),
)
fig, axes = plt.subplots(1, 5, figsize=(21, 3.8))
for ax, (label, field, y_scale) in zip(axes, fields, strict=True):
for frame in snapshots:
ax.plot(r_mm, field[frame], label=f"{t_ns[frame]:.0f} ns", lw=1.3)
if y_scale == "log":
ax.set_yscale("log")
ax.set_xlabel("r [mm]")
ax.set_ylabel(label)
axes[0].legend(fontsize="small")
fig.suptitle("Radial profiles of the expanding channel")
What the frozen run cannot show.#
Overlaying the two at the same late time separates “the model has no flow” from “the flow does nothing”. The frozen channel holds its density and keeps its heat on the axis; the expanding one has moved three quarters of the axis mass outward and cooled by doing the work.
late = len(t_out) - 1
frozen_Tg = np.asarray(frozen["Tg"])
frozen_ne = np.asarray(frozen["ne"])
frozen_rho = np.asarray(frozen["rho"])
fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))
for ax, (label, a, b, y_scale) in zip(
axes,
(
(r"$T_\mathrm{H}$ [K]", Tg[late], frozen_Tg[late], None),
(r"$n_\mathrm{e}$ [m$^{-3}$]", electron_density[late], frozen_ne[late], "log"),
(r"$\rho$ [kg m$^{-3}$]", density[late], frozen_rho[late], None),
),
strict=True,
):
ax.plot(r_mm, a, lw=1.6, label="expanding")
ax.plot(r_mm, b, "--", lw=1.3, label="frozen")
if y_scale == "log":
ax.set_yscale("log")
ax.set_xlabel("r [mm]")
ax.set_ylabel(label)
axes[0].legend(fontsize="small")
fig.suptitle(f"At t = {t_ns[late]:.0f} ns: with and without the radial sub-step")
What a lumped expansion cannot do.#
The 0D polytropic cell expands against ambient pressure, which is the right
idea in one dimension too few. Its density falls as V0/V, driven by the
mean pressure of the cell; the 1D core evacuates into a rarefaction running
inward from a face that is not there in the lumped model. Plotting the two
density ratios against each other is the sharpest statement of the
difference, and the right-hand panels show what it costs: the 0D electron
temperature and resistance follow a cell that never really rarefies.
fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))
axes[0].plot(t_ns, density[:, 0] / density[0, 0], lw=1.6, label="1D, on axis")
axes[0].plot(
zero_d_expanding["t"] * 1e9,
zero_d_expanding["V"][0] / zero_d_expanding["V"],
":",
color="0.35",
lw=1.5,
label="0D polytropic cell",
)
axes[0].plot(
zero_d_frozen["t"] * 1e9,
zero_d_frozen["V"][0] / zero_d_frozen["V"],
"-.",
color="0.6",
lw=1.2,
label="0D constant volume",
)
axes[0].axhline(1.0, color="0.8", lw=0.6)
axes[0].set_ylabel(r"$\rho / \rho_0$ [-]")
axes[0].legend(fontsize="small")
axes[1].plot(t_ns, Te[:, 0], lw=1.6, label="1D, on axis")
axes[1].plot(
zero_d_expanding["t"] * 1e9,
zero_d_expanding["Te"],
":",
color="0.35",
lw=1.5,
label="0D polytropic",
)
axes[1].set_ylabel(r"$T_\mathrm{e}$ [K]")
axes[1].legend(fontsize="small")
# `Rp_state`, not `Rp`: the 0D resistance is read after that model's own
# advance, so the comparable 1D quantity is the one belonging to the recorded
# state rather than the one the circuit was driven with over the step.
axes[2].plot(t_ns, np.asarray(expanding["Rp_state"]), lw=1.6, label="1D expanding")
axes[2].plot(t_ns, np.asarray(frozen["Rp_state"]), "--", lw=1.2, label="1D frozen")
axes[2].plot(
zero_d_expanding["t"] * 1e9,
zero_d_expanding["Rp"],
":",
color="0.35",
lw=1.5,
label="0D polytropic",
)
axes[2].plot(
zero_d_frozen["t"] * 1e9,
zero_d_frozen["Rp"],
"-.",
color="0.6",
lw=1.2,
label="0D constant volume",
)
axes[2].set_yscale("log")
axes[2].set_ylabel(r"$R_p$ [$\Omega$]")
axes[2].legend(fontsize="small")
for ax in axes:
ax.set_xlabel("t [ns]")
fig.suptitle("Radially resolved against lumped: the same pulse, two geometries")
plt.show()
An animation of the expansion.#
Written to a file rather than embedded: the gallery does not execute the
1dplasma examples, so this renders only when the script is run directly.
# sphinx_gallery_thumbnail_number = 3
def make_video(path=None):
"""Animate the radial profiles of the expanding channel over time.
Parameters
----------
path : str, optional
Output path. A ``.gif`` of the same stem is written instead when no
ffmpeg writer is available. The default writes into the system
temporary directory rather than the working directory: a bare relative
name drops a megabyte of binary beside the source whenever the example
is run from the repository, and sphinx-gallery executes example
scripts with ``__name__ == "__main__"``.
Returns
-------
str or None
The path actually written, or None when no writer succeeded.
"""
if path is None:
path = str(pathlib.Path(tempfile.gettempdir()) / "nrp_expansion.mp4")
from matplotlib.animation import FFMpegWriter, FuncAnimation, PillowWriter
panels = (
(r"$T_\mathrm{H}$ [K]", Tg, None),
(r"$T_\mathrm{e}$ [K]", Te, None),
(r"$n_\mathrm{e}$ [m$^{-3}$]", electron_density, "log"),
(r"$p$ [atm]", pressure / ct.one_atm, None),
(r"$u_r$ [m s$^{-1}$]", velocity, None),
)
fig, axes = plt.subplots(1, 5, figsize=(21, 3.8))
lines = []
for ax, (label, field, y_scale) in zip(axes, panels, strict=True):
(line,) = ax.plot(r_mm, field[0], lw=1.6, color="teal")
lines.append(line)
if y_scale == "log":
ax.set_yscale("log")
positive = field[field > 0.0]
ax.set_ylim(positive.min() * 0.5, field.max() * 2.0)
else:
margin = 0.05 * (field.max() - field.min())
ax.set_ylim(field.min() - margin, field.max() + margin)
ax.set_xlim(0.0, r_mm.max())
ax.set_xlabel("r [mm]")
ax.set_ylabel(label)
title = fig.suptitle("")
def update(frame):
"""Redraw every panel at output frame `frame`."""
for line, (_, field, _) in zip(lines, panels, strict=True):
line.set_ydata(field[frame])
title.set_text(f"7 kV CH4 channel, expanding: t = {t_ns[frame]:5.0f} ns")
return (*lines, title)
anim = FuncAnimation(fig, update, frames=len(t_out), interval=40, blit=False)
for writer, out in (
(FFMpegWriter(fps=25), path),
(PillowWriter(fps=25), path.replace(".mp4", ".gif")),
):
try:
anim.save(out, writer=writer)
plt.close(fig)
return out
except Exception as exc: # noqa: BLE001 - writer availability varies
print(f" ({type(writer).__name__} unavailable: {exc})")
plt.close(fig)
return None
if __name__ == "__main__":
written = make_video()
print(f"animation written to {written}" if written else "no animation writer")