Cantera does not compute a multi-temperature ideal gas law for now

[pag1pag]

Problem description

Hello, and thanks for Cantera! I think I have found a physics bug in the plasma phase, or maybe there is something that I do not understand.

A multi-temperature phase, like the Plasma phase, should take into account electron temperature when computing pressure. I think this is not the case in Cantera, when I compare the formula used in Chemkin.

Chemkin

In Chemkin Manual, equation 2-4

$$P = \sum_k n_k k_b T_k = k_b T_g \sum_{k \neq e} n_k + n_e k_b T_e$$

Cantera

In Cantera (electronic temperature is not taken into accont to compute the pressure)

$$P = \sum_k n_k k_b T_g = k_b T_g \sum_{k \neq e} n_k + n_e k_b T_g$$

For the Ideal Gas Mixture phase, the molar volume of a species is given by the ideal gas law:

$$V_k^o(T, P)=\frac{R T}{P}$$

For the Plasma phase, there is an electron pressure defined. However, the class Plasma inherits from Ideal Gas Mixture, but the method setPressure from the Ideal Gas Mixture phase is only inherited, and not redefined to consider the electron pressure.

Steps to reproduce

import cantera as ct

# Create the plasma object using a default YAML file.
plasma = ct.Solution(
    "example_data/oxygen-plasma-itikawa.yaml",
    "isotropic-electron-energy-plasma",
    transport_model=None,
)

# Initialize the plasma object with a temperature, pressure, and composition.
plasma.TDX = 300, 0.6415, {"O2": 0.5, "E": 0.5}
# Also initialize the electron temperature.
plasma.Te = 300.0

print(
    "Initializing the plasma object:\n"
    + f"T = {plasma.T} K\n"
    + f"rho = {plasma.density} kg/m^3\n"
    + f"X = {plasma.X}\n"
    + f"Te = {plasma.Te} K\n"
)

print(f"Computed (total?) pressure = {plasma.P} Pa")
# >>> 100011.93190795329 Pa, ok
print(f"Computed (electron) pressure = {plasma.Pe} Pa")
# >>> 50005.96595397664 Pa, ok

print("==========================")

# Now, changing the electron temperature.
plasma.Te = 30_000.0

print(
    f"Changing the electron temperature to T_e={plasma.Te} K:\n"
    + f"T = {plasma.T} K\n"
    + f"rho = {plasma.density} kg/m^3\n"
    + f"X = {plasma.X}\n"
)
print(f"Computed (total?) pressure = {plasma.P} Pa")
# >>> 100011.93190795329 Pa, NOT ok?
print(f"Computed (electron) pressure = {plasma.Pe} Pa")
# >>> 5000596.595397665, ok

Behavior

In first case, gas temperature and electron temperature are set, as well as density and mole fraction (50% O2 and 50 % E). The pressure computed is around 1 atm, which is correct.

Then, in the second case, electron temperature is multiply by a factor 100, leaving gas temperature, density and mole fraction unchanged. However, the pressure is still exactly the same as before, which does not seem correct.

I also think there is confusion between total pressure and gas pressure:

• If plasma.P is the total pressure
  • In case 1, $P_{tot}=1atm$ and $P_e = 0.5 atm$ --> , $P_{gas}=0.5 atm$ which is good.
  • In case 2, $P_{tot}=1atm$ and $P_e = 5 atm$ --> $P_{gas}=P_{tot}-P_e=-4.5 atm$, not good...
• If plasma.P is the gas pressure
  • In case 1, $P_{gas}=1atm$ and $P_e = 0.5 atm$ --> $P_{tot}=1.5 atm$ --> not what I expected
  • In case 2, $P_{gas}=1atm$ and $P_e = 5 atm$ --> $P_{tot}=P_{gas}+P_e= 6atm$ , maybe?

System information

• Cantera version: 3.1.0
• OS: Windows 11
• Python versions: 3.11.10

[ischoegl]

Thanks for reporting. Tagging @BangShiuh.

[BangShiuh]

@pag1pag Thanks for documenting this. There are more properties need to be updated together (to be consistent) for PlasmaPhase. For now, it is "okay" if the electron number density is low. This is on my to-do list but it would be helpful if other people can also join the development.

[erwanp]

It is "okay" as long as the electron number density is much lower than the total density, isnt it ?

[BangShiuh]

Here is a brief overview for properties that needs to be changed and its dependencies.

Proposed definition:

[speth]

Updating PlasmaPhase to implement the equation of state such that the electron's contribution is proportional to $T_e$ makes sense to me. I think doing so provides a self-consistent approach to implementing intEnergy_mole and partialMolarIntEnergies:

If you define the partial pressures as
$P_k = R T_k X_k \rho_m$
where $\rho_m$ is the mixture molar density, then $P = \sum P_k$. Further,
$u_k = h_k - \frac{P_k}{X_k \rho_m} = h_k - R T_k$
is also consistent with
$u = h - P/\rho_m$

I think defining the normal ThermoPhase methods to exclude the electron and introducing others specifically for the electron makes it hard to achieve consistency among the properties, and would also make the interface much more complicated.

[pag1pag]

Hello all!

I have tried to implement this equation of state (see here, should I open a PR now?). The main changes are the following:

Note that the following assumes that a temperature can be defined for all species. While this may be true for all heavy species, electrons may not be characterized by an electron temperature (for instance, if the electron energy distribution function (EEDF) is far away from a Maxwellian). It is therefore assumed that the implementation below holds if, for instance, the EEDF is Maxwellian (which may be seen as a strong approximation in plasma).

The partial pressure is defined as: $P_k = R T_k X_k \rho_m$ (with $\rho_m$ the mixture molar density).
The total pressure is the sum of all partial pressure: $P = \sum_k P_k = R \rho_m \sum_k X_k T_k$
Introducing a mean temperature, defined by $\overline{T}=\sum_k X_k T_k$, the total pressure is rewritten $P = R \rho_m \overline{T}$.

The standard concentration (standardConcentration) uses this mean temperature: $c^°=\frac{P}{R \overline{T}}$.
The standard volume (getStandardVolumes and getPartialMolarVolumes) is $v_k = \frac{R T_k}{P}$.

In the implementation before, some properties (like getPartialMolarEntropies, which depend on total pressure) were redefined using the electron partial pressure (for electron species). It may not be correct, since this would suggest to use the partial pressure of the other species (which is not the case). Therefore, the use of electron pressure was removed.

As was already implemented, $c_{p, k}$, $h_k$, $s_k$, $g_k$ are evaluated at the respective species temperature $T_k$.
With the definition of the molar volume ($v_k = \frac{R T_k}{P}$), the internal energy is defined as $u_k = h_k - P v_k = h_k - R T_k$.

---

Setting the state

Moreover, a new setState function is required, to directly define the electron temperature. It is only done if the distribution is isotropic (for discretized or Boltzmann-term, an electron temperature can be obtained from the distribution function, but not the other way around).
Regarding the method TPX, for a two temperature plasma, what is T? Is it the gas temperature or the electron temperature? I assumed both are equal, and introduced a new TgTeP method to set separetely gas and electron temperature.

---

Test

Looking now at consistency tests, all of them pass when the electron temperature is equal to gas temperature.

A new case was added with different electron and gas temperatures. The failing tests are due to temperature differential, which are hard to define when there are multiple temperature (When $T_e \neq T_g$, $c_p=\frac{dh}{dT}$ may not be correct: is $dT$ equal to $dT_h$ or $dT_e$ ?)

Six tests are failing for now (test/python/test_thermo/TestPlasmaPhase). Maybe it is due to the change in molar concentration, introduced by a different electron temperature? I have not yet fix this issue.

---

New methods and file refactoring

I also introduced some new functions to get the heat capacity of electrons (only) and heavy species (only).
And I also refactor plasma thermo files by regrouping related method in the same section (no change of code, just move some part for lisibility)

---

A multi-temperature class?

While doing this, an idea came to my mind. Maybe a full multi-temperature class should be implemented, from which a two-temperature plasma could be derived? IdealGasPhase could also derived from this multi-temperature class, and set all temperature equal to the gas temperature?

[speth]

Thanks for digging into this. There's some overlap with what @mquiram has done in #2068. That PR is more focused on getting the reactor and equilibrium calculations to work and give reasonable results for low electron concentrations, but I think you're right that there are a few additional considerations beyond just the definitions for the partial molar species properties.

I agree with you that the usage of the "electron pressure" is currently a bit problematic.

Updating setState to accept an electron temperature makes sense. I'm less keen on the introduction of the TgTeP methods for setting the full state -- a full implementation of this would quickly propagate to variants for also setting the mass or mole fractions from strings, maps, or arrays. In the case of an isotropic EEDF, I think just letting the TP method set both the gas and electron temperature is fine, with the setElectronTemperature method being called afterward to modify the state. Likewise, I think this would let you leave the rest of the state setter functions intact, with the assumption being that they set the electron temperature equal to the gas temperature for an isotropic EEDF, or leave the electron temperature fixed otherwise.

I'd like to avoid adding the functions like cp_mass_e and friends if possible, unless it turns out that they're really necessary for some purpose.