A global model in five steps

Download notebook View on GitHubRun when the site was built.

A global model balances, averaged over a chamber, the density of every species and the electron energy. This tutorial builds the simplest one, an argon discharge with one ion species, checks it against the closed form of the textbook model (Lieberman and Lichtenberg, Principles of Plasma Discharges, ch. 10), and looks at every term of its balances. The cross sections are illustrative step functions, not argon data.

1. Species and reactions

Rates come from cross sections averaged over a Maxwellian, with plasmacoll.

from plasmacoll import CrossSection
from plasmacoll.constants import ELECTRON_MASS

from wellbalanced import ELECTRON, Chamber, GlobalModel, Reaction, Species

argon = Species.atomic("Ar", 39.948)
argon_ion = Species.atomic("Ar+", 39.948, 1)
ionization = Reaction.from_cross_section(
    ("e", "Ar"),
    ("e", "e", "Ar+"),
    CrossSection.constant(2.0e-20, threshold=15.76),
    ELECTRON_MASS,
    electron_energy_loss=15.76,
    name="ionization",
)
[ionization.coefficient(te, 0.025) for te in (1.0, 3.0, 10.0)]  # m^3/s

[3.2105657372649555e-20, 7.538506733822054e-16, 2.2548193997125172e-14]

2. The chamber and the model

One litre, 600 cm² of wall, ions leaving at the Bohm speed with an edge-to-centre ratio of 0.5; 100 W absorbed by the electrons; the gas held at 10²⁰ m⁻³.

chamber = Chamber(volume=1.0e-3, wall_area=0.06, h_factor=0.5)
model = GlobalModel(
    species=[ELECTRON, argon, argon_ion],
    reactions=[ionization],
    chamber=chamber,
    power=100.0,
    fixed={"Ar": 1.0e20},
)
state = model.steady_state()
state.converged, state.electron_temperature, state.electron_density

(True, 3.0473789456058626, 2.0375338762410694e+17)

3. Every term of the balances

At steady state the electron power terms (W) sum to zero, and every species’ rate of change vanishes.

state.power

{'absorbed': 100.0,
 'ionization': -41.8737678082905,
 'walls': -58.12623219170946}

sum(state.power.values()), state.particle_balance

(3.552713678800501e-14,
 {'Ar': -1.6583470923687066e+22, 'Ar+': -54525952.0, 'e': -54525952.0})

4. The closed form

The particle balance alone sets the electron temperature, $K_{iz}(T_e)\, n_g = h u_B A / V$ ; the power balance then sets the density.

import math

from plasmacoll.constants import ELEMENTARY_CHARGE as e
from scipy.optimize import brentq


def bohm(te):
    return math.sqrt(e * te / argon_ion.mass)


te = brentq(
    lambda t: ionization.coefficient(t, 0.025) * 1.0e20 - bohm(t) * 0.5 * 0.06 / 1.0e-3,
    0.5,
    50.0,
)
te, state.electron_temperature

(3.04737894560586, 3.0473789456058626)

5. A scan

Each steady state starts from the last. The electron temperature does not depend on the power, the density is proportional to it.

states = model.scan("power", [50.0, 100.0, 200.0, 400.0])
[(s.electron_temperature, s.electron_density) for s in states]

[(3.0473789456058635, 1.018766938120533e+17),
 (3.0473789456058613, 2.0375338762410694e+17),
 (3.047378945605862, 4.0750677524821466e+17),
 (3.0473789456058635, 8.150135504964251e+17)]

6. As a dataset

Every result is also a labelled xarray dataset, and a scan one dataset along the scanned parameter (wellbalanced.stack): units and names travel with the numbers, for plotting (with plasma-plots or xarray’s own .plot) and for saving (.to_netcdf).

from wellbalanced import stack

powers = [50.0, 100.0, 200.0, 400.0]
scan = stack(model.scan("power", powers), "power", powers, units="W")
scan.density.sel(species="e")

<xarray.DataArray 'density' (power: 4)> Size: 32B
array([1.01876694e+17, 2.03753388e+17, 4.07506775e+17, 8.15013550e+17])
Coordinates:
  * power    (power) float64 32B 50.0 100.0 200.0 400.0
    species  <U3 12B 'e'
Attributes:
    long_name:  density
    units:      m^-3