rizer.plasma.reactor#

Custom Cantera reactor solving energy equation.

Here, Cantera is used for evaluating thermodynamic properties and kinetic rates while an external ODE solver is used to integrate the resulting equations. In this case, the SciPy wrapper for VODE is used, which uses the same variable-order BDF methods as the Sundials CVODES solver used by Cantera.

Initial discussions with Bang-Shiuh Chen, 03/05/2023

> > BOLOS is in python. It is slower than ZDPlasKin(bolsig). > You will need to write a customized solver in python (use scipy ode solver). > BOLOS will output plasma properties such as electron temperature, electron mobility, > and collision reaction rates. You will then use function such as gas.multiplier to set the reaction rates. > > > https://cantera.org/examples/python/reactors/custom.py.html > > This is an example of using python to solve a customized problem. > I have a similar script for NRP plasma. It was built with my professor > (Sally Bane) at Purdue University. We need to collaborate officially if you > want to use it. I would say zero D is not too difficult to build. You should > give it a shot.

Attributes#

`logger`
`u`

Classes#

`PlasmaExtension`	Extend a Cantera plasma object.
`PlasmaReactorOde`
`ConstantPressurePlasmaReactorOde`	Parameters of the ODE system and auxiliary data are stored in the ReactorOde object.

Module Contents#

rizer.plasma.reactor.logger#

rizer.plasma.reactor.u#

class rizer.plasma.reactor.PlasmaExtension(plasma: cantera.Solution, name: str = 'gas', heavy_species_mechanism: str | pathlib.Path = 'gri30.yaml', discard_species: list[str] = ['e-'], missing_mole_fraction_warning_threshold: float = 0.01, bypass_temperature_clipping: bool = False, threshold_main_species_by_mole: float = 0.9)#

Extend a Cantera plasma object.

Extend a Cantera plasma object to include properties not currently implemented in Cantera, such as cp_mole.

Parameters:

plasma (ct.Solution) – Cantera plasma object
heavy_species_mechanism (str | Path, optional) – Mechanism of heavy species, by default “gri30.yaml”
name (str, optional) –
Name of the heavy species object, by default “gas”.

From Cantera.Solution documentation:
”If an input file defines multiple phases, the corresponding key in the phases map can be used to specify the desired phase via the name keyword argument of the constructor”
discard_species (list[str], optional) – List of species to discard from the Plasma object when computing properties not implemented in PlasmaPhase (ex: discard electrons e when computing cp_mass), by default [“e-“]
missing_mole_fraction_warning_threshold (float, optional) – Warn if the sum of mole fractions of species not taken into account is above this threshold, by default 1e-2.
bypass_temperature_clipping (bool, optional) –
Bypass the temperature clipping when updating the heavies object, by default False. Species has valid thermo data between a certain temperature range. If the gas temperature is above this range, the temperature is clipped to the maximum temperature. Note that this temperature range depends on the species (see threshold_main_species_by_mole).
- If True, the temperature of the heavies object is not clipped.
- If False, it is clipped to min(max temperature of main species, gas temperature).
threshold_main_species_by_mole (float, optional) –
Threshold for the main species by mole fraction, by default 0.9. Only used if bypass_temperature_clipping is False, to compute the maximum temperature of the heavies object. For instance, thermo data for CH4 are valid up to 6000 K ; and up 20 000 K for H.
- If there are more than 90%mol of CH4, and the gas temperature is 10 000 K, the temperature of the heavies object will be clipped to 6000 K.
- If there are more than 90%mol of H, and the gas temperature is 10 000 K, the temperature of the heavies object will not be clipped.
- If there are 50%mol of CH4 and 50%mol of H, and the gas temperature is 10 000 K, the temperature of the heavies object will be clipped to 6000 K.

Raises:

ValueError – The mechanism of the heavies object should not be a Plasma phase.

plasma#

heavies#

missing_mole_fraction_warning_threshold = 0.01#

bypass_temperature_clipping = False#

threshold_main_species_by_mole = 0.9#

e_index#

heavies_index = []#

plasma_species_in_heavies_mask#

momentum_transfer_collision_frequencies_list: list[rizer.plasma.collision_frequency.MomentumTransferCollisionFrequencyModel] = []#

update_heavies() → None#: Update the heavies object with the mole fractions from plasma.

get_main_species_by_mole(threshold_main_species_by_mole: float) → list#

density() → float#

Compute the density of the plasma in kg/m^3.

Cantera does not take into account the fact that electrons are at a different temperature than the heavy species. This method updates the electron concentration and computes the density of the plasma.

Returns:: Density of the plasma in kg/m^3.
Return type:: float

Notes

The density of the plasma is given (eq. 2-5 in [Chemkin]):

\[\rho = \sum_k \[ X_k \] M_k\]

with:

\(\[ X_k \]\) the molar concentration of species \(k\),
\(M_k\) the molar mass of species \(k\).

The molar concentration can be expressed using the mole fraction (eq. 2-13 in [Chemkin]):

\[\[ X_k \] = \frac{X_k P}{R \sum_i X_i T_i}\]

In our case, all heavies species are at temperature \(T_g\), but the electrons are at temperature \(T_e\).

Therefore, the density is:

\[\rho = \sum_k \frac{X_k P}{R \sum_i X_i T_i} M_k = \frac{P}{R \sum_i X_i T_i} \sum_k X_k M_k = \frac{P}{R \left( \sum_{i \neq e) X_i T_i + X_e T_e \right) } \sum_k X_k M_k\]

mean_temperature() → float#

Compute the mean temperature of the plasma in K.

This is the equation of state used in [Chemkin]. It is not yet implemented in Cantera.

Returns:: Mean temperature of the plasma in K.
Return type:: float

Notes

The mean temperature of the plasma is given by the formula:

\[\overline{T} = \frac{\sum_k X_k T_k}{\sum_k X_k} = \sum_k X_k T_k\]

with:

\(X_k\) the molar fraction of species \(k\),
\(T_k\) the temperature of species \(k\).

This can be rewritten as:

\[\overline{T} = \sum_{k \neq e} X_k T_g + X_e T_e = (1 - X_e) T_g + X_e T_e\]

This comes from the gas equation of state (equation 2-4 of [Chemkin]):

\[P = \sum_k \[ X_k \] R T_k = \sum_k n_k k_b T_k = k_b n_tot \overline{T}\]

with:

\(P\) the pressure,
\(\[ X_k \]\) the molar concentration of species \(k\),
\(n_k\) the number density of species \(k\),
\(n_tot\) the total number density of the plasma.

total_number_density() → float#

Compute the total number density of the plasma in m^-3.

Returns:: Total number density of the plasma in m^-3.
Return type:: float

Notes

The total number density of the plasma is given by the formula:

\[\rho_{\text{tot}} = \frac{P}{k_b \overline{T}}\]

with:

\(P\) the pressure, in Pa,
\(k_b\) the Boltzmann constant, in J/K,
\(\overline{T}\) the mean temperature of the plasma, in K.

References

Equation 2-4 of [Chemkin].

electron_number_density() → float#

Compute the electron number density of the plasma in m^-3.

Returns:: Electron density of the plasma in m^-3.
Return type:: float

Notes

The electron number density of the plasma is given by the formula:

\[n_e = X_{e^-} n_{\text{tot}} = X_{e^-} \frac{P}{k_b \overline{T}}\]

with:

\(X_{e^-}\) the electron mole fraction,
\(P\) the pressure, in Pa,
\(k_b\) the Boltzmann constant, in J/K,
\(\overline{T}\) the mean temperature of the plasma, in K.

References

Equation 2-4 of [Chemkin].

cp_mass_mean_heavies() → float#

Mean massic heat capacity at constant pressure of the heavies species in J/kg/K.

Returns:: Mean massic heat capacity at constant pressure of the heavies species in J/kg/K.
Return type:: float

Notes

The mean massic heat capacity at constant pressure of the heavies species is given by the formula:

\[\overline{c_{p, mass}} = \sum_{k \neq e} Y_k c_{p, k}\]

with:

\(Y_k\) the mass fraction of species \(k\), in the plasma phase,
\(c_{p, k}\) the heat capacity of species \(k\).

It can be rewritten as:

\[\]

overline{c_{p, mass}} = (1 - Y_e) sum_{k} Y_k’ c_{p, k}

with:

\(Y_k'\) the mass fraction of species \(k\), in the heavy phase,
\(Y_e\) the electron mass fraction, in the plasma phase.
\(c_{p, k}\) the heat capacity of species \(k\).

Indeed, one the one hand, \(\sum_{k} Y_k = \sum_{\neq e} Y_k + Y_e = 1.0\). On the other hand, \(\sum_{k} Y_k' = \sum_{\neq e} Y_k' = 1.0\), since Y_e’ = 0. In the heavy phase, the mass fraction are normalized to 1.0, thus \(Y_k' = \frac{Y_k}{Y_e}\).

References

Equation 30 of [Aurora].

See also

cp_mass(), cp_e_mass()

cv_mass_mean_heavies() → float#

Mean massic heat capacity at constant volume of the heavies species in J/kg/K.

Returns:: Mean massic heat capacity at constant volume of the heavies species in J/kg/K.
Return type:: float

See also

cp_mass()

cp_mole_mean_heavies() → float#

Mean molar heat capacity at constant pressure of the heavies species in J/kmol/K.

Returns:: Mean molar heat capacity at constant pressure of the heavies species in J/kmol/K.
Return type:: float

See also

cp_mole()

cp_mass() → float#

Massic heat capacity at constant pressure of the plasma in J/kg/K.

Returns:: Massic heat capacity at constant pressure of the plasma in J/kg/K.
Return type:: float

Notes

The massic heat capacity at constant pressure of the plasma is given by the formula:

\[c_{p, mass} = (1 - Y_{e-}) \cdot c_{p, mass, heavy} + Y_{e-} \cdot c_{p, mass, e^-}\]

with:

\(Y_{e-}\) the electron mass fraction,
\(c_{p, mass, heavy}=\sum_{k \neq e^-} Y_k c_{p, k}\) the heat capacity of the heavy species,
\(c_{p, mass, e^-}\) the heat capacity of electrons.

References

Equation 30 of [Aurora].

See also

cp_e_mass()

cp_mole() → float#

Molar heat capacity at constant pressure of the plasma in J/kmol/K.

Returns:: Molar heat capacity at constant pressure of the plasma in J/kmol/K.
Return type:: float

Notes

The molar heat capacity at constant pressure of the plasma is given by the formula:

\[c_{p, mol} = (1 - X_{e-}) \cdot c_{p, mol, heavy} + X_{e-} \cdot c_{p, mol, e^-}\]

with:

\(X_{e-}\) the electron mole fraction,
\(c_{p, mol, heavy}=\sum_{k \neq e^-} X_k c_{p, k}\) the heat capacity of the heavy species,
\(c_{p, mol, e^-}\) the heat capacity of electrons.

References

Equation 30 of [Aurora].

See also

cp_e_mole()

cp_e_mass() → float#

Massic heat capacity at constant pressure of electrons in J/kg/K.

Returns:: Massic heat capacity at constant pressure of electrons in J/kg/K.
Return type:: float

Notes

The massic heat capacity at constant pressure of electrons is given by the formula:

\[c_{p, mass, e} = \frac{5}{2} \frac{R}{M_e}\]

with:

\(R\) the ideal gas constant (J/mol/K),
\(M_e\) the electron molar mass (kg/mol).

References

Equation 33 of [Aurora].

cv_e_mass() → float#

Massic heat capacity at constant volume of electrons in J/kg/K.

Returns:: Massic heat capacity at constant volume of electrons in J/kg/K.
Return type:: float

Notes

The massic heat capacity at constant volume of electrons is given by the formula:

\[c_{v, mass, e} = \frac{3}{2} \frac{R}{M_e}\]

with:

\(R\) the ideal gas constant (J/mol/K),
\(M_e\) the electron molar mass (kg/mol).

References

Equation 33 of [Aurora].

cp_e_mole() → float#

Molar heat capacity at constant pressure of electrons in J/kmol/K.

Returns:: Molar heat capacity at constant pressure of electrons in J/kmol/K.
Return type:: float

Notes

The molar heat capacity at constant pressure of electrons is given by the formula:

\[c_{p, mol, e} = \frac{5}{2} R_{\text{kmol}}\]

with:

\(R_{\text{kmol}}\) the ideal gas constant, expressed in J/kmol/K.

References

Equation 33 of [Aurora].

cv_e_mole() → float#

Molar heat capacity at constant volume of electrons in J/kmol/K.

Returns:: Molar heat capacity at constant volume of electrons in J/kmol/K.
Return type:: float

Notes

The molar heat capacity at constant volume of electrons is given by the formula:

\[c_{v, mol, e} = \frac{3}{2} R_{\text{kmol}}\]

with:

\(R_{\text{kmol}}\) the ideal gas constant, expressed in J/kmol/K.

References

Equation 33 of [Aurora].

plasma_power_E(E: float) → float#

Return the Joule heating power density in W/m3 as a function of E.

For now, it is assumed that the plasma is weakly ionized.

Parameters:: E (float) – External electric field in V/m
Returns:: Joule heating power density in W/m3
Return type:: float

Notes

The power density is given by the formula:

\[P_W = j \cdot E\]

with:

\(j=\sigma \cdot E\) the current density,
\(E\) the electric field.

The electrical conductivity \(\sigma\) is given by the formula:

\[\sigma = \frac{1}{ \frac{1}{\sigma^{weakly}} + \frac{1}{\sigma^{strongly}}}\]

with:

\(\sigma^{weakly}\) the weakly ionized electrical conductivity,
\(\sigma^{strongly}\) the strongly ionized electrical conductivity.

See also

electrical_conductivity(), weakly_ionized_electrical_conductivity(), weakly_ionized_electrical_conductivity(), strongly_ionized_electrical_conductivity(), strongly_ionized_electrical_conductivity()

plasma_power_sfactor(δ: float = 1000.0) → float#

Return the Joule heating power (elastic+inelastic) per unit volume of the plasma in W/m^3.

The inelastic power is scaled from the elastic power by a default factor of δ=1000.

For details about the formula used, see [LauxLecture].

Parameters:: δ (float, optional) – Scaling factor for inelastic cross-section, by default 1000.
Returns:: Joule heating power density in W/m3
Return type:: float

Notes

The power density is given by the formula [LauxLecture]:

\[P = n_e \cdot (\delta + 1) \cdot 2 \cdot \frac{m_e}{m_h} \cdot \frac{3}{2} \cdot k_b \cdot (T_e - T) \cdot n_h \cdot \sigma \cdot v_e\]

with:

\(n_e\) the electron density,
\(\delta\) the scaling factor for inelastic cross-section,
\(m_e\) the electron mass,
\(m_h\) the mass of heavy particles,
\(T_e\) the electron temperature,
\(T\) the heavy species temperature,
\(n_h\) the heavy species density,
\(\sigma\) the cross-section,
\(v_e\) the electron thermal velocity.

plasma_power_elastic() → float#

Compute the elastic power loss density in W/m^3.

Return the rate of electron energy loss per unit volume as a result of elastic collisions with heavy particles in W/m^3.

Returns:: Elastic power loss density in W/m^3.
Return type:: float

Notes

The power density is given by equation (VI 5.1) of [Mitchner1973].

\[P_el = \sum_h \frac{2 m_e}{m_h} \frac{3}{2} k_b\left(T_e-T_g\right) \bar{\nu}_{e h} n_e\]

with:

\(m_e\) the electron mass, in kg,
\(m_h\) the mass of heavy particles, in kg,
\(k_b\) the Boltzmann constant, in J/K,
\(T_e\) the electron temperature, in K,
\(T_g\) the heavy species temperature, in K,
\(\bar{\nu}_{e h}\) the energy-weighted average momentum transfer collision frequency between
electrons and heavy particles, in s^-1, as defined in (II 6.29) of [Mitchner1973].

References

Equation 35 of [Aurora]
End of chapter 5 of [LauxLecture]
equation (VI 5.1) of [Mitchner1973]

plasma_power_inelastic() → float#

Return the inelastic Joule heating power per unit volume of the plasma in W/m^3.

Returns:: Inelastic Joule heating power density in W/m^3.
Return type:: float

Notes

The power density is given by equation 36 of [Aurora]:

\[P_{inel} = \sum_{i}^{I_{ei}} \Delta H_i \cdot R_i\]

with:

\(I_{ei}\) the number of reactions that depend on electron-temperature,
\(\Delta H_i\) the enthalpy change of reaction \(i\),
\(R_i\) the net rate of progress of reaction \(i\).

electrical_conductivity() → float#

Return the electrical conductivity in S/m.

No assumption is made on wether the plasma is weakly or strongly ionized.

Returns:: electrical conductivity in S/m
Return type:: float

Notes

The electrical conductivity \(\sigma\) of a plasma is given by (II 13.7b) of [Mitchner1973]:

\[\sigma_e = \frac{n_e e^2}{m_e \bar{\nu}_{eH}}\]

where:

\(n_e\) is the electron number density in m^-3,
\(e\) is the elementary charge in C,
\(m_e\) is the electron mass in kg,
\(\bar{\nu}_{eH}\) is the average momentum transfer collision frequency of an electron with all heavy particle species.

This last term is defined in (II 13.3) of [Mitchner1973] as the following sum:

\[\bar{\nu}_{eH} = \bar{\nu}_{en} + \bar{\nu}_{ei}\]

where:

\(\bar{\nu}_{en}\) is the average momentum transfer collision frequency of an electron with neutral heavy particles,
\(\bar{\nu}_{ei}\) is the average momentum transfer collision frequency of an electron with ionized heavy particles.

thermal_conductivity() → float#

Return the thermal conductivity of the plasma in W/m/K.

The mixture thermal conductivity is computed using [Mathur1967] formula. It is computed at an average temperature between the gas temperature and the external temperature.

For now, only the conductivities of H2 and CH4 are available. For other species, the thermal conductivity is assumed to be the same as CH4. Therefore, if the mixture contains other gases, the thermal conductivity is not accurate.

TODO: Implement the thermal conductivity of other species. TODO: Use https://cantera.org/dev/python/transport.html

Returns:: thermal conductivity in W/m/K
Return type:: float

Notes

The thermal conductivity \(\lambda\) is given by [Mathur1967]:

\[\lambda = \frac{1}{2} \left \sum_k \frac{X_k}{\lambda_k} + \frac{1}{\sum_k X_k / \lambda_k} \right)\]

with:

\(X_k\) the mole fraction of species \(k\),
\(\lambda_k\) the thermal conductivity of species \(k\).

The thermal conductivity of a species is for now found in the following references:

https://ntrs.nasa.gov/api/citations/19950021761/downloads/19950021761.pdf

conductive_heat_loss() → float#

Return the heat loss per unit volume of the plasma in W/m^3.

Returns:: Heat loss in W/m^3.
Return type:: float

Notes

The heat loss is given by the formula (p199 of [Raizer1991], eq. 8 of [Heijkers2020]):

\[q = 8 \cdot \frac{\lambda_{\text{mix}} \cdot (T_{\text{gas}} - T_{\text{ext}})}{r^2}\]

with:

\(\lambda_{\text{mix}}\) the thermal conductivity of the mixture,
\(T_{\text{gas}}\) the gas temperature,
\(T_{\text{ext}}\) the external temperature,
\(r\) the radius of the plasma.

See also

thermal_conductivity()

radiative_heat_loss(epsilon: float = 0.1) → float#

Return the radiative heat loss per unit volume of the plasma in W/m^3.

Simple formula using the Stefan-Boltzmann law. A warning is raised to inform the user that this formula is a simple approximation.

Parameters:: epsilon (float, optional) – Emissivity of the plasma, by default 0.1
Returns:: Radiative heat loss in W/m^3.
Return type:: float

Warning

The current implementation of the radiative heat loss formula uses a black body models. This radiation model is far away real thermal plasma emitting in UV and visible range. Use this model with caution.

Notes

The radiative heat loss is given by the formula:

\[q = \epsilon \cdot \sigma^{Stefan} \cdot (T_{\text{gas}}^4 - T_{\text{ext}}^4)\]

with:

\(\epsilon\) the emissivity of the plasma,
\(\sigma^{Stefan}\) the Stefan-Boltzmann constant,
\(T_{\text{gas}}\) the gas temperature,
\(T_{\text{ext}}\) the external temperature.

References

See spark-cleantech-l3/rizer#83

heat_loss(include_radiation: bool = False) → float#

Return the total heat loss per unit volume of the plasma in W/m^3.

Parameters:: include_radiation (bool, optional) – Include radiative heat loss, by default False until a proper model is implemented.
Returns:: Total heat loss in W/m^3.
Return type:: float

Notes

The total heat loss is given by the formula:

\[q = q_{\text{conductive}} + q_{\text{radiative}}\]

with:

\(q_{\text{conductive}}\) the conductive heat loss,
\(q_{\text{radiative}}\) the radiative heat loss.

class rizer.plasma.reactor.PlasmaReactorOde(plasma: cantera.Solution, plasma_extension: PlasmaExtension, energy_enabled=True)#

plasma#

P#

plasma_ext#

energy_enabled = True#

class rizer.plasma.reactor.ConstantPressurePlasmaReactorOde(plasma: cantera.Solution, plasma_extension: PlasmaExtension, electric_field: Callable[[float], float] | None = None)#

Parameters of the ODE system and auxiliary data are stored in the ReactorOde object.

Parameters:

plasma (ct.Solution) – Cantera plasma object.
plasma_extension (PlasmaExtension) – Object that extends the plasma object to include properties not currently implemented in Cantera, such as cp_mole.
electric_field (Callable[[float], float] | None, optional) – Electric field, as a function of time, by default None. If None, the electric field is set to 0.0 V/m at all times. The function should take a time in seconds as input and return the electric field in V/m.

Notes

Equations solved implies no mass flow, and homogeneous quantities only.

The equations solved are based on [Aurora].

The conservation of species:

\[\frac{dY_k}{dt} = \frac{W_k \dot{\omega_k}}{\rho}\]

with \(Y_k\) the mass fraction of species \(k\), \(W_k\) the molecular weight of species \(k\), \(\omega_k\) the production rate of species \(k\), and \(\rho\) the density.

The conservation of energy for the electrons:

\[Y_e c_{v, e} \frac{dT_e}{dt} = \frac{R T_e}{M_e} \frac{d Y_e}{dt} + \frac{\dot{\omega} M_e c_{p, e} (T_g - T_e)}{\rho} + \frac{\vec{j} \cdot \vec{E} - \dot{Q}_{elas} - \dot{Q}_{inel}}{\rho}\]

The conservation of energy for the gas:

\[\bar{c_p} \frac{dT_g}{dt} = - \frac{R}{M_e} \frac{d (T_e Y_e) }{dt} - \sum_{k \neq e} \frac{W_k h_k \dot{\omega_k}}{\rho} - \frac{W_e c_{p, e} T_g \dot{\omega_e}}{\rho} + \frac{\dot{Q}_{elas} + \dot{Q}_{inel} - \dot{Q}_{loss}}{\rho}\]

See also

equations

Examples

plasma#

P#

plasma_ext#