Thermodynamic properties for a 2-temperature plasma

[pag1pag]

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). See also Cantera/enhancements#258 for details about the physics behind the code.

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}$.

Changes proposed in this pull request

• Define or Redefine for a 2-temperature system:
  • Molar Thermodynamic Properties: enthalpy_mole, entropy_mole, gibbs_mole, intEnergy_mole
  • Mechanical Equation of State: meanTemperature, pressure, setPressure, electronPressure
  • Chemical Potentials and Activities: standardConcentration
  • Partial Molar Properties: getChemPotentials, getPartialMolarEnthalpies, getPartialMolarIntEnergies, getPartialMolarVolumes
  • Properties of the Standard State of the Species: getStandardChemPotentials, getStandardVolumes
  • Thermodynamic Values for the Species Reference: getGibbs_ref, getStandardVolumes_ref
• Update the setState method to account for additional electron temperature. Also introduce setState_TgTeP and setState_TgTeD to set electron temperature along with other quantities.
• Looking now at consistency tests, all of them pass when the electron temperature is equal to gas temperature. New cases were added with different electron and gas temperatures. The skipped 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$ ?). Some Python tests required to update the electron temperature after setting TPX, and doing that gives the expected value.
• Refactor and document plasma thermo files by regrouping related method in the same section (no change of code, just move some parts for better lisibility)

If applicable, fill in the issue number this pull request is fixing

Closes #1872

If applicable, provide an example illustrating new features this pull request is introducing

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, at the same temperature as the heavy particles.
plasma.Te = 300.0

print(
    "Initializing the plasma object:\n"
    + f"T = {plasma.T} K\n"  # >>> 300 K, ok
    + f"rho = {plasma.density} kg/m^3\n"  # >>> 0.6415 kg/m^3, ok
    + f"X = {plasma.X}\n"  # >>> [0.5 0.  0.5 0.  0.  0. ], ok
    + f"Te = {plasma.Te} K\n"  # >>> 300 K, ok
)

print(f"Computed 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"  # >>> 300 K, ok
    + f"rho = {plasma.density} kg/m^3\n"  # >>> 0.6415 kg/m^3, ok
    + f"X = {plasma.X}\n"  # >>> [0.5 0.  0.5 0.  0.  0. ], ok
    + f"Te = {plasma.Te} K\n"  # >>> 30000 K, ok
)
print(f"Computed (total) pressure = {plasma.P} Pa")
# >>> 5050602.561351642 Pa, ok
print(f"Computed (electron) pressure = {plasma.Pe} Pa")
# >>> 5000596.595397665, ok
# The total pressure is now dominated by the electron pressure, as expected.

AI Statement (required)

• No generative AI was used. This contribution was written entirely without AI
  assistance.

Checklist

• The pull request includes a clear description of this code change
• Commit messages have short titles and reference relevant issues
• Build passes (scons build & scons test) and unit tests address code coverage
• Style & formatting of contributed code follows contributing guidelines
• AI Statement is included
• The pull request is ready for review

[codecov]

Codecov Report

❌ Patch coverage is 69.14498% with 83 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.65%. Comparing base (44732e7) to head (c626963).
⚠️ Report is 9 commits behind head on main.

Files with missing lines	Patch %	Lines
src/thermo/PlasmaPhase.cpp	68.42%	35 Missing and 43 partials ⚠️
include/cantera/thermo/PlasmaPhase.h	73.68%	5 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2101      +/-   ##
==========================================
- Coverage   77.71%   77.65%   -0.06%     
==========================================
  Files         452      452              
  Lines       53316    53412      +96     
  Branches     8894     8914      +20     
==========================================
+ Hits        41432    41479      +47     
- Misses       8870     8899      +29     
- Partials     3014     3034      +20     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[pag1pag]

Hello !

I have reread the Contributing guide, and I just seen that I have used merge instead of rebase in this work. Is it still possible to use git rebase, or is it too late?

Also, I have opened an enhancement request, to discuss the physics behind the code.

[pag1pag]

As for the CI/CD tests,

Locally, after running scons build and scons test, I see that two tests are failing, but they do not seem related to this PR:

Tests passed: 7378
Up-to-date tests skipped: 0
Tests failed: 2

Failed tests:
    - python: test_reactor.TestMoleReactor.test_tolerances
    - cxx-openmp-ignition

As for the examples, samples\python\thermo\plasma-eedf.py seems to be failing online. It seems to be due to the UserWarning being raised during the execution of the file, but I do not understand why the code crashes on a warning.

[speth]

I have reread the Contributing guide, and I just seen that I have used merge instead of rebase in this work. Is it still possible to use git rebase, or is it too late?

I think you can rebase your local branch on top of our main and it will eliminate the merge commit. Once that is done, you can do a "force push" to this PR.

One other thing: in commit 41ea41e, the files consistency.cpp and thermo.pyx got converted to Windows-style line endings, which makes them show up as completely changed files and it's impossible to see the actual code changes (if any). Could you please fix those files? The best way would be to do so as part of an interactive rebase so the erroneous change never occurs. Otherwise, we can just squash the PR later if you don't want to worry about the detailed commit history.

[pag1pag]

Hello, I have rebased, rewrote the history and force push.
I hope it is clearer now :D

[speth]

Thanks for working on this, @pag1pag. Overall I think this is a significant improvement in the consistency of how we handle the thermo properties when there are two distinct temperatures. However, there is still one discrepancy that we need to resolve, which involves the definitions of the standard concentration, the activities, and the species standard molar volumes. The implementation here defines a single value for the standard concentration, $C^o \equiv P / R \bar{T}$, and specifies that the activities are equal to the mole fractions, $a_k \equiv X_k$. However, this is not consistent with the definition of the partial molar volume $V^o_k \equiv R T_k / P$ when $T_e \ne T_g$.

Since the current kinetics implementation relies on there being only a single standard concentration for a phase, I think the best fix for this is to define the activity for the electron as $a_e \equiv X_e T_e / T$ while leaving $a_k = X_k$ for the gas phase species. This is approximate, but correct in the limit of small $X_e$, since $\bar{T} \approx T_g$.

I'm not quite sure what the value of introducing the methods setState_TgTeP and setState_TgTeD is, over just calling setElectronTemperature followed by setState_TP or setState_TD (other than the fact that the latter currently overwrite the electron temperature). If these methods left the electron temperature alone, I think we could avoid adding these specialized methods (which are difficult to access in any case).

As an editorial note, the combination of formatting updates (specifically, rearranging the methods in both the header and implementation) with substantive code changes made this quite challenging to review. It's very difficult to identify where there are actual implementation changes versus just moving functions around. I would strongly recommend not mixing such changes in a single PR in the future.

[speth]

Suggested change

	* For an ideal gas, enthalpy is indepedant of pressure.
	* For an ideal gas, internal energy is independent of pressure.

[speth]

Suggested change

	* For an ideal gas, enthalpy is indepedant of pressure.
	* For an ideal gas, enthalpy is independent of pressure.

[speth]

Add here:

.. versionadded:: 4.0

[speth]

These default to false, so presumably these lines can just be deleted now.

[speth]

I'd prefer to leave out these decorative dividers; the named method groups can stay.

[speth]

I know we aren't broadly consistent about this, but I don't think we need a different decoration for mixture molar enthalpy and partial molar enthalpy. My preference would be to keep $\hat{h}_k$ as the partial molar enthalpies, as before. Likewise, I'm not sure it's worth the churn on changing the notation for the standard-state from superscript 0 to superscript ref here. I think it would be better to stay notationally consistent with IdealGasPhase so the meaningful differences stand out more clearly.

[speth]

I'm not sure this is the best choice for how to handle setting T for either the isotropic or full EEDF-based cases. In the latter case, clearly we should not set the EEDF based on the gas temperature. But then there seems to be no way to set the temperature of the gas without triggering this warning (as now seen in various tests and the examples). On the other hand, it also seems strange to reset the electron state by setting T only in the case of an isotropic distribution.

My intuition would be that we should leave the electron state alone when setting T -- let that just be the gas temperature, and have that be true regardless of how the EEDF is specified. That also seems more consistent with introduction of the meanTemperature property.

[speth]

I wonder if we can support the full set of options allowed at the ThermoPhase level, and avoid reimplementing much of that method, by simply setting the electron temperature first and then calling the base class method.

[speth]

This guard should apply to setMeanElectronEnergy as well, which can be achieved by having that method call setElectronTemperature with the computed value.

[speth]

I think this description is a bit misleading -- the reason these methods do not need to be overridden is the same as for getEnthalpy_RT and friends: because the electron temperature correction is handled by updateThermo. More generally, I think it would be more clear to document this all in the active sense as documentation for the updateThermo method, explaining that this method handles getting properties at the correct temperature for each species, rather than explaining the things we don't have to do here.