diff --git a/.github/workflows/linters.yml b/.github/workflows/linters.yml index c6fb4497304..230badaffa3 100644 --- a/.github/workflows/linters.yml +++ b/.github/workflows/linters.yml @@ -65,3 +65,25 @@ jobs: git config --global core.whitespace \ -cr-at-eol,tab-in-indent,blank-at-eol,blank-at-eof git diff --check ${{ github.event.pull_request.base.sha }} + example-data: + name: Check for unmerged example-data pull request + runs-on: ubuntu-latest + if: github.event_name == 'pull_request' + steps: + - uses: actions/checkout@v4 + name: Checkout the repository + with: + fetch-depth: 100 + persist-credentials: true + + - name: Find matching PR in example_data + env: + BASE_PR: ${{ github.event.number }} + GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + cd data/example_data + DATA_PR=$(gh pr list --repo Cantera/cantera-example-data --jq '.[] | select(.title | test("Cantera/cantera#'${BASE_PR}'")) | .number' --json title,number) + if [ -n "$DATA_PR" ]; then + echo ":exclamation: Merge https://github.com/Cantera/cantera-example-data/pull/${DATA_PR} and update the submodule commit before merging this PR" >> $GITHUB_STEP_SUMMARY + exit 1 + fi diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 5bf79292b11..37dc676b443 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -105,6 +105,18 @@ jobs: with: submodules: recursive persist-credentials: false + - name: Check out updated example data + if: github.event_name == 'pull_request' + env: + BASE_PR: ${{ github.event.number }} + GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + cd data/example_data + gh repo set-default Cantera/cantera-example-data + DATA_PR=$(gh pr list --repo Cantera/cantera-example-data --jq '.[] | select(.title | test("Cantera/cantera#'${BASE_PR}'")) | .number' --json title,number) + if [ -n "$DATA_PR" ]; then + gh pr checkout --repo Cantera/cantera-example-data $DATA_PR + fi - name: Setup Python uses: actions/setup-python@v5 with: @@ -419,6 +431,17 @@ jobs: with: submodules: recursive persist-credentials: false + - name: Check out updated example data + if: github.event_name == 'pull_request' + env: + BASE_PR: ${{ github.event.number }} + GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + cd data/example_data + DATA_PR=$(gh pr list --repo Cantera/cantera-example-data --jq '.[] | select(.title | test("Cantera/cantera#'${BASE_PR}'")) | .number' --json title,number) + if [ -n "$DATA_PR" ]; then + gh pr checkout --repo Cantera/cantera-example-data $DATA_PR + fi - name: Set up micromamba uses: mamba-org/setup-micromamba@b09ef9b599704322748535812ca03efb2625677b # v2.0.5 with: diff --git a/AUTHORS.md b/AUTHORS.md index d776a824692..41258fa3347 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -12,6 +12,7 @@ update, please report on Cantera's - **Halla Ali** [@hallaali](https://github.com/hallaali) - **Emil Atz** [@EmilAtz](https://github.com/EmilAtz) - **Jongyoon Bae** [@jongyoonbae](https://github.com/jongyoonbae) - Brown University +- **Nicolas Barléon** [@NicolasBarleon](https://github.com/NicolasBarleon) - CERFACS - **Philip Berndt** - **Guus Bertens** [@guusbertens](https://github.com/guusbertens) - Eindhoven University of Technology - **Wolfgang Bessler** [@wbessler](https://github.com/wbessler) - Offenburg University of Applied Science @@ -49,6 +50,7 @@ update, please report on Cantera's - **Kyle Linevitch, Jr.** [@KyleLinevitchJr](https://github.com/KyleLinevitchJr) - **Christopher Leuth** - **Nicholas Malaya** [@nicholasmalaya](https://github.com/nicholasmalaya) - University of Texas at Austin +- **Quentin Malé** [@QuentinMale](https://github.com/QuentinMale) - Harvard University - **Thanasis Mattas** [@ThanasisMattas](https://github.com/ThanasisMattas) - Aristotle University of Thessaloniki - **Manik Mayur** [@manikmayur](https://github.com/manikmayur) - **Evan McCorkle** [@emccorkle](https://github.com/emccorkle) @@ -61,6 +63,7 @@ update, please report on Cantera's - **Paul Northrop** [@pwcnorthrop](https://github.com/pwcnorthrop) - **Chris Pilko** [@cpilko](https://github.com/cpilko) - **Sebastian Pinnau** [@spinnau](https://github.com/spinnau) +- **Matthew Quiram** [@mquiram](https://github.com/mquiram) - Massachusetts Institute of Technology - **Corey R. Randall** [@c-randall](https://github.com/c-randall) - Colorado School of Mines - **Andreas Rücker** [@cannondale1492](https://github.com/cannondale1492) - **Jeff Santner** [@jsantner](https://github.com/jsantner) diff --git a/doc/sphinx/yaml/index.md b/doc/sphinx/yaml/index.md index b70c7b97bd1..2209c1c2074 100644 --- a/doc/sphinx/yaml/index.md +++ b/doc/sphinx/yaml/index.md @@ -31,6 +31,8 @@ state](sec-yaml-species-eos), [transport property](sec-yaml-species-transport), YAML reaction definitions include specification of common elements such as the reaction equation and [efficiencies](sec-yaml-efficiencies), as well as parameters specific to the type of [rate parameterization](sec-yaml-rate-types). +See also the {ref}`electron collision data format ` +used in plasma-phase simulations. ## Mechanism Conversion diff --git a/doc/sphinx/yaml/phases.md b/doc/sphinx/yaml/phases.md index d61abb129bb..811a6d04ccd 100644 --- a/doc/sphinx/yaml/phases.md +++ b/doc/sphinx/yaml/phases.md @@ -862,6 +862,7 @@ Additional fields: - `isotropic` - `discretized` + - `Boltzmann-two-term` `shape-factor` : A constant in the isotropic distribution, which is shown as x in the detailed @@ -909,6 +910,9 @@ Examples: normalize: False ``` +See also {ref}`sec-yaml-electron-collisions`, +which is used to specify electron-impact cross sections relevant to plasma simulations. + :::{versionadded} 2.6 ::: diff --git a/doc/sphinx/yaml/reactions.md b/doc/sphinx/yaml/reactions.md index e1984a40451..963a5426974 100644 --- a/doc/sphinx/yaml/reactions.md +++ b/doc/sphinx/yaml/reactions.md @@ -18,6 +18,7 @@ The fields common to all `reaction` entries are: - [`Blowers-Masel`](sec-yaml-Blowers-Masel) - [`two-temperature-plasma`](sec-yaml-two-temperature-plasma) - [`electron-collision-plasma`](sec-yaml-electron-collision-plasma) + - [`electron-collisions`](sec-yaml-electron-collisions) - [`falloff`](sec-yaml-falloff) - [`chemically-activated`](sec-yaml-chemically-activated) - [`pressure-dependent-Arrhenius`](sec-yaml-pressure-dependent-Arrhenius) @@ -308,6 +309,48 @@ Example: :::{versionadded} 3.1 ::: +(sec-yaml-electron-collisions)= +### `electron-collisions` + +The `electron-collisions` field defines a list of cross-section datasets for +electron-impact processes that are used in plasma-phase simulations. These entries +are not formal reactions (they are not added to `Kinetics` objects), but serve +as data inputs for computing the electron energy distribution function. + +Each entry includes: + +`target` +: The name of the species that is the target of the collision + +`energy-levels` +: A list of electron energy values [eV] at which the cross-section is provided + +`cross-sections` +: Corresponding cross-section values [m²] for each energy level + +`kind` +: A string indicating the process type. Options include: + - `"effective"` – lumped or total effect of several channels + - `"excitation"` – electronic excitation + - `"ionization"` – electron-impact ionization + - `"attachment"` – electron attachment processes + +Example: + +```yaml +electron-collisions: +- target: N2 + energy-levels: [0.0, 0.015, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.7, 1.2, 1.5, 1.9, + 2.2, 2.8, 3.3, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 75.0, 150.0] + cross-sections: [1.1e-20, 2.55e-20, 3.4e-20, 4.33e-20, 5.95e-20, 7.1e-20, 7.9e-20, + 9e-20, 9.7e-20, 1e-19, 1.04e-19, 1.2e-19, 1.96e-19, 2.85e-19, 2.8e-19, 1.72e-19, + 1.26e-19, 1.09e-19, 1.01e-19, 1.04e-19, 1.1e-19, 1.02e-19, 9e-20, 6.6e-20, 4.9e-20] + kind: effective +``` + +:::{versionadded} 3.2 +::: + (sec-yaml-falloff)= ### `falloff` diff --git a/include/cantera/base/utilities.h b/include/cantera/base/utilities.h index 26d0c3ea817..5f6d6975769 100644 --- a/include/cantera/base/utilities.h +++ b/include/cantera/base/utilities.h @@ -179,7 +179,7 @@ void checkFinite(const double tmp); * @param values Array of *N* values to be checked * @param N Number of elements in *values* */ -void checkFinite(const string& name, double* values, size_t N); +void checkFinite(const string& name, const double* values, size_t N); //! Const accessor for a value in a map. /* diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index e627998e841..fb0a2156eca 100644 --- a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -7,10 +7,10 @@ #ifndef CT_ELECTRONCOLLISIONPLASMARATE_H #define CT_ELECTRONCOLLISIONPLASMARATE_H -#include "cantera/thermo/PlasmaPhase.h" #include "cantera/kinetics/ReactionData.h" #include "ReactionRate.h" #include "MultiRate.h" +#include "cantera/numerics/eigen_dense.h" namespace Cantera { @@ -32,19 +32,17 @@ struct ElectronCollisionPlasmaData : public ReactionData energyLevels.resize(0); distribution.resize(0); m_dist_number = -1; - m_level_number = -1; } vector energyLevels; //!< electron energy levels vector distribution; //!< electron energy distribution - bool levelChanged; + + //! integer that is incremented when electron energy levels change + int levelNumber = -1; protected: //! integer that is incremented when electron energy distribution changes int m_dist_number = -1; - - //! integer that is incremented when electron energy level changes - int m_level_number = -1; }; @@ -103,9 +101,16 @@ struct ElectronCollisionPlasmaData : public ReactionData class ElectronCollisionPlasmaRate : public ReactionRate { public: - ElectronCollisionPlasmaRate() = default; - + //! Constructor from YAML input for ElectronCollisionPlasmaRate. + /*! + * This constructor is used to initialize an electron collision plasma rate + * from an input YAML file. It extracts the energy levels, cross-sections, + * and reaction metadata used in the rate coefficient calculation. + * + * @param node The AnyMap node containing rate fields from YAML + * @param rate_units Units used for interpreting the rate fields + */ ElectronCollisionPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}) { @@ -133,6 +138,10 @@ class ElectronCollisionPlasmaRate : public ReactionRate */ double evalFromStruct(const ElectronCollisionPlasmaData& shared_data); + //! Calculate the reverse rate coefficient for super-elastic collisions + //! @param shared_data Data structure with energy levels and EEDF + //! @param kf Forward rate coefficient (input, unused) + //! @param kr Reverse rate coefficient (output, modified) void modifyRateConstants(const ElectronCollisionPlasmaData& shared_data, double& kf, double& kr); @@ -143,6 +152,47 @@ class ElectronCollisionPlasmaRate : public ReactionRate throw NotImplementedError("ElectronCollisionPlasmaRate::ddTScaledFromStruct"); } + //! The kind of the process which will be one of the following: + //! - `"effective"`: A generic effective collision + //! - `"excitation"`: Electronic or vibrational excitation + //! - `"ionization"`: Electron-impact ionization + //! - `"attachment"`: Electron attachment + //! @since New in Cantera 3.2. + const string& kind() const { + return m_kind; + } + + //! Get the target species of the electron collision process. + //! This is the name of the neutral or ionic species that the electron interacts with + //! @since New in Cantera 3.2. + const string& target() const { + return m_target; + } + + //! Get the product of the electron collision process. + //! This is the name of the species or excited state of + //! some species resulting from the process. + //! @note this may not necessarily be represented by + //! a distinct species in the mixture. + //! @since New in Cantera 3.2. + const string& product() const { + return m_product; + } + + //! Get the energy threshold of electron collision [eV] + //! + //! By default, the threshold is set to the first non-zero energy value + //! listed in the tabulated cross section data. + //! + //! @note This behavior may be subject to change. A more robust approach + //! may use the energy corresponding to the first non-zero cross section + //! value, rather than the first non-zero energy point in the data table. + //! + //! @since New in Cantera 3.2. + double threshold() const { + return m_threshold; + } + //! The value of #m_energyLevels [eV] const vector& energyLevels() const { return m_energyLevels; @@ -162,9 +212,28 @@ class ElectronCollisionPlasmaRate : public ReactionRate void updateInterpolatedCrossSection(const vector&); private: + //! The name of the kind of electron collision + string m_kind; + + //! The name of the target of electron collision + string m_target; + + //! The product of electron collision + string m_product; + + //! The energy threshold of electron collision + double m_threshold; + //! electron energy levels [eV] vector m_energyLevels; + //! Counter used to indicate when #m_energyLevels needs to be synced with the phase + int m_levelNumber = -3; + + //! Counter used to indicate when #m_crossSectionsOffset needs to be synced with the + //! phase + int m_levelNumberSuperelastic = -2; + //! collision cross sections [m2] at #m_energyLevels vector m_crossSections; diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index ffcf1827df0..f7f40d6329b 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -1474,9 +1474,6 @@ class Kinetics Eigen::SparseMatrix m_stoichMatrix; //! @} - //! Boolean indicating whether Kinetics object is fully configured - bool m_ready = false; - //! The number of species in all of the phases //! that participate in this kinetics mechanism. size_t m_kk = 0; diff --git a/include/cantera/thermo/EEDFTwoTermApproximation.h b/include/cantera/thermo/EEDFTwoTermApproximation.h new file mode 100644 index 00000000000..bd39137006e --- /dev/null +++ b/include/cantera/thermo/EEDFTwoTermApproximation.h @@ -0,0 +1,288 @@ +/** + * @file EEDFTwoTermApproximation.h EEDF Two-Term approximation solver. + */ + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_EEDF_TWO_TERM_APPROXIMATION_H +#define CT_EEDF_TWO_TERM_APPROXIMATION_H + +#include "cantera/base/ct_defs.h" +#include "cantera/numerics/eigen_sparse.h" + +namespace Cantera +{ + +class PlasmaPhase; + +//! Boltzmann equation solver for the electron energy distribution function based on +//! the two-term approximation. +/*! + * This class implements a solver for the electron energy distribution function + * based on a steady-state solution to the Boltzmann equation using the classical + * two-term expansion, applicable to weakly ionized plasmas. The numerical approach + * and theory are primarily derived from the work of Hagelaar and Pitchford + * @cite hagelaar2005. + * + * The two-term approximation assumes that the EEDF can be represented as: + * @f[ + * f(\epsilon, \mu) = f_0(\epsilon) + \mu f_1(\epsilon), + * @f] + * where @f$ \epsilon @f$ is the electron energy and @f$ \mu @f$ is the cosine + * of the angle between the electron velocity vector and the electric field. + * The Boltzmann equation is projected onto the zeroth moment over mu to obtain + * an equation for @f$ f_0(\epsilon) @f$ , the isotropic part of the distribution. + * The first-order anisotropic term @f$ f_1(\epsilon) @f$ is not solved directly, + * but is approximated and substituted into the drift and collision terms. This + * results in a second-order differential equation for @f$ f_0(\epsilon) @f$ alone. + * + * The governing equation for @f$ f_0(\epsilon) @f$ is discretized on an energy + * grid using a finite difference method and solved using a tridiagonal matrix + * algorithm. + * + * @since New in %Cantera 3.2. + * @warning This class is an experimental part of %Cantera and may be changed without + * notice. + */ +class EEDFTwoTermApproximation +{ +public: + EEDFTwoTermApproximation() = default; + + //! Constructor combined with the initialization function + /*! + * This constructor initializes the EEDFTwoTermApproximation object with everything + * it needs to start solving EEDF. + * + * @param s PlasmaPhase object that will be used in the solver calls. + */ + EEDFTwoTermApproximation(PlasmaPhase* s); + + virtual ~EEDFTwoTermApproximation() = default; + + // compute the EEDF given an electric field + // CQM The solver will take the species to consider and the set of cross-sections + // from the PlasmaPhase object. + // It will write the EEDF and its grid into the PlasmaPhase object. + // Successful returns are indicated by a return value of 0. + int calculateDistributionFunction(); + + void setLinearGrid(double& kTe_max, size_t& ncell); + + void setGridCache(); + + vector getGridEdge() const { + return m_gridEdge; + } + + vector getEEDFEdge() const { + return m_f0_edge; + } + + double getElectronMobility() const { + return m_electronMobility; + } + +protected: + + //! Formerly options for the EEDF solver + + //! The first step size + double m_delta0 = 1e14; + + //! Maximum number of iterations + size_t m_maxn = 200; + + //! The factor for step size change + double m_factorM = 4.0; + + //! The number of points in the EEDF grid + size_t m_points = 150; + + //! Error tolerance for convergence + double m_rtol = 1e-5; + + //! The growth model of EEDF + std::string m_growth = "temporal"; + + //! The threshold for species mole fractions + double m_moleFractionThreshold = 0.01; + + //! The first guess for the EEDF + std::string m_firstguess = "maxwell"; + + //! The initial electron temperature [eV] + double m_init_kTe = 2.0; + + //! Pointer to the PlasmaPhase object used to initialize this object. + /*! + * This PlasmaPhase object provides species, element, and cross-section + * data used by the EEDF solver. It is set during construction and is not + * modified afterwards. All subsequent calls to compute functions must + * use the same PlasmaPhase context. + */ + PlasmaPhase* m_phase; + + //! Iterate f0 (EEDF) until convergence + void converge(Eigen::VectorXd& f0); + + //! An iteration of solving electron energy distribution function + Eigen::VectorXd iterate(const Eigen::VectorXd& f0, double delta); + + //! The integral in [a, b] of \f$x u(x) \exp[g (x_0 - x)]\f$ + //! assuming that u is linear with u(a) = u0 and u(b) = u1 + double integralPQ(double a, double b, double u0, double u1, + double g, double x0); + + //! Vector g is used by matrix_P() and matrix_Q(). + /** + * \f[ + * g_i = \frac{1}{\epsilon_{i+1} - \epsilon_{i-1}} \ln(\frac{F_{0, i+1}}{F_{0, i-1}}) + * \f] + */ + vector vector_g(const Eigen::VectorXd& f0); + + //! The matrix of scattering-out. + /** + * \f[ + * P_{i,k} = \gamma \int_{\epsilon_i - 1/2}^{\epsilon_i + 1/2} + * \epsilon \sigma_k exp[(\epsilon_i - \epsilon)g_i] d \epsilon + * \f] + */ + Eigen::SparseMatrix matrix_P(const vector& g, size_t k); + + //! The matrix of scattering-in + /** + * \f[ + * Q_{i,j,k} = \gamma \int_{\epsilon_1}^{\epsilon_2} + * \epsilon \sigma_k exp[(\epsilon_j - \epsilon)g_j] d \epsilon + * \f] + */ + //! where the interval \f$[\epsilon_1, \epsilon_2]\f$ is the overlap of cell j, + //! and cell i shifted by the threshold energy: + /** + * \f[ + * \epsilon_1 = \min(\max(\epsilon_{i-1/2}+u_k, \epsilon_{j-1/2}),\epsilon_{j+1/2}), + * \f] + * \f[ + * \epsilon_2 = \min(\max(\epsilon_{i+1/2}+u_k, \epsilon_{j-1/2}),\epsilon_{j+1/2}) + * \f] + */ + Eigen::SparseMatrix matrix_Q(const vector& g, size_t k); + + //! Matrix A (Ax = b) of the equation of EEDF, which is discretized by the exponential scheme + //! of Scharfetter and Gummel, + /** + * \f[ + * \left[ \tilde{W} F_0 - \tilde{D} \frac{d F_0}{\epsilon} \right]_{i+1/2} = + * \frac{\tilde{W}_{i+1/2} F_{0,i}}{1 - \exp[-z_{i+1/2}]} + + * \frac{\tilde{W}_{i+1/2} F_{0,i+1}}{1 - \exp[z_{i+1/2}]} + * \f] + * where \f$ z_{i+1/2} = \tilde{w}_{i+1/2} / \tilde{D}_{i+1/2} \f$ (Peclet number). + */ + Eigen::SparseMatrix matrix_A(const Eigen::VectorXd& f0); + + //! Reduced net production frequency. Equation (10) of ref. [1] + //! divided by N. + //! @param f0 EEDF + double netProductionFrequency(const Eigen::VectorXd& f0); + + //! Diffusivity + double electronDiffusivity(const Eigen::VectorXd& f0); + + //! Mobility + double electronMobility(const Eigen::VectorXd& f0); + + //! Initialize species indices associated with cross-section data + void initSpeciesIndexCrossSections(); + + //! Update the total cross sections based on the current state + void updateCrossSections(); + + //! Update the vector of species mole fractions + void updateMoleFractions(); + + //! Compute the total elastic collision cross section + void calculateTotalElasticCrossSection(); + + //! Compute the total (elastic + inelastic) cross section + void calculateTotalCrossSection(); + + //! Compute the L1 norm of a function f defined over a given energy grid. + //! + //! @param f Vector representing the function values (EEDF) + //! @param grid Vector representing the energy grid corresponding to f + double norm(const Eigen::VectorXd& f, const Eigen::VectorXd& grid); + + //! Electron mobility [m²/V·s] + double m_electronMobility; + + //! Grid of electron energy (cell center) [eV] + Eigen::VectorXd m_gridCenter; + + //! Grid of electron energy (cell boundary i-1/2) [eV] + vector m_gridEdge; + + //! Location of cell j for grid cache + vector> m_j; + + //! Location of cell i for grid cache + vector> m_i; + + //! Cross section at the boundaries of the overlap of cell i and j + vector>> m_sigma; + + //! The energy boundaries of the overlap of cell i and j + vector>> m_eps; + + //! normalized electron energy distribution function + Eigen::VectorXd m_f0; + + //! EEDF at grid edges (cell boundaries) + vector m_f0_edge; + + //! Total electron cross section on the cell center of energy grid + vector m_totalCrossSectionCenter; + + //! Total electron cross section on the cell boundary (i-1/2) of + //! energy grid + vector m_totalCrossSectionEdge; + + //! vector of total elastic cross section weighted with mass ratio + vector m_sigmaElastic; + + //! list of target species indices in global Cantera numbering (1 index per cs) + vector m_kTargets; + + //! list of target species indices in local X EEDF numbering (1 index per cs) + vector m_klocTargets; + + //! Indices of species which has no cross-section data + vector m_kOthers; + + //! Local to global indices + vector m_k_lg_Targets; + + //! Mole fraction of targets + vector m_X_targets; + + //! Previous mole fraction of targets used to compute eedf + vector m_X_targets_prev; + + //! in factor. This is used for calculating the Q matrix of + //! scattering-in processes. + vector m_inFactor; + + double m_gamma; + + //! flag of having an EEDF + bool m_has_EEDF; + + //! First call to calculateDistributionFunction + bool m_first_call; +}; // end of class EEDFTwoTermApproximation + +} // end of namespace Cantera + +#endif \ No newline at end of file diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 266fbfe6466..7670bf0b3c0 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -10,6 +10,7 @@ #define CT_PLASMAPHASE_H #include "cantera/thermo/IdealGasPhase.h" +#include "cantera/thermo/EEDFTwoTermApproximation.h" #include "cantera/numerics/eigen_sparse.h" namespace Cantera @@ -108,9 +109,7 @@ class PlasmaPhase: public IdealGasPhase //! @param levels The vector of electron energy levels (eV). //! Length: #m_nPoints. //! @param length The length of the @c levels. - //! @param updateEnergyDist update electron energy distribution - void setElectronEnergyLevels(const double* levels, size_t length, - bool updateEnergyDist=true); + void setElectronEnergyLevels(const double* levels, size_t length); //! Get electron energy levels. //! @param levels The vector of electron energy levels (eV). Length: #m_nPoints @@ -220,11 +219,24 @@ class PlasmaPhase: public IdealGasPhase return m_nPoints; } - //! Number of collisions + //! Number of electron collision cross sections size_t nCollisions() const { return m_collisions.size(); } + //! Get the Reaction object associated with electron collision *i*. + //! @since New in %Cantera 3.2. + const shared_ptr collision(size_t i) const { + return m_collisions[i]; + } + + //! Get the ElectronCollisionPlasmaRate object associated with electron collision + //! *i*. + //! @since New in %Cantera 3.2. + const shared_ptr collisionRate(size_t i) const { + return m_collisionRates[i]; + } + //! Electron Species Index size_t electronSpeciesIndex() const { return m_electronSpeciesIndex; @@ -283,6 +295,9 @@ class PlasmaPhase: public IdealGasPhase void setParameters(const AnyMap& phaseNode, const AnyMap& rootNode=AnyMap()) override; + //! Update the electron energy distribution. + void updateElectronEnergyDistribution(); + //! Electron species name string electronSpeciesName() const { return speciesName(m_electronSpeciesIndex); @@ -298,6 +313,57 @@ class PlasmaPhase: public IdealGasPhase return m_levelNum; } + //! Get the indicies for inelastic electron collisions + //! @since New in %Cantera 3.2. + const vector& kInelastic() const { + return m_kInelastic; + } + + //! Get the indices for elastic electron collisions + //! @since New in %Cantera 3.2. + const vector& kElastic() const { + return m_kElastic; + } + + //! target of a specific process + //! @since New in %Cantera 3.2. + size_t targetIndex(size_t i) const { + return m_targetSpeciesIndices[i]; + } + + //! Get the frequency of the applied electric field [Hz] + //! @since New in %Cantera 3.2. + double electricFieldFrequency() const { + return m_electricFieldFrequency; + } + + //! Get the applied electric field strength [V/m] + double electricField() const { + return m_electricField; + } + + //! Set the absolute electric field strength [V/m] + void setElectricField(double E) { + m_electricField = E; + } + + //! Calculate the degree of ionization + //double ionDegree() const { + // double ne = concentration(m_electronSpeciesIndex); // [kmol/m³] + // double n_total = molarDensity(); // [kmol/m³] + // return ne / n_total; + //} + + //! Get the reduced electric field strength [V·m²] + double reducedElectricField() const { + return m_electricField / (molarDensity() * Avogadro); + } + + //! Set reduced electric field given in [V·m²] + void setReducedElectricField(double EN) { + m_electricField = EN * molarDensity() * Avogadro; // [V/m] + } + virtual void setSolution(std::weak_ptr soln) override; /** @@ -344,9 +410,6 @@ class PlasmaPhase: public IdealGasPhase */ void checkElectronEnergyDistribution() const; - //! Update electron energy distribution. - void updateElectronEnergyDistribution(); - //! Set isotropic electron energy distribution void setIsotropicElectronEnergyDistribution(); @@ -391,8 +454,28 @@ class PlasmaPhase: public IdealGasPhase //! Flag of normalizing electron energy distribution bool m_do_normalizeElectronEnergyDist = true; - //! Data for initiate reaction - AnyMap m_root; + //! Indices of inelastic collisions in m_crossSections + vector m_kInelastic; + + //! Indices of elastic collisions in m_crossSections + vector m_kElastic; + + //! electric field [V/m] + double m_electricField = 0.0; + + //! electric field freq [Hz] + double m_electricFieldFrequency = 0.0; + + //! Cross section data. m_crossSections[i][j], where i is the specific process, + //! j is the index of vector. Unit: [m^2] + vector> m_crossSections; + + //! Electron energy levels corresponding to the cross section data. m_energyLevels[i][j], + //! where i is the specific process, j is the index of vector. Unit: [eV] + vector> m_energyLevels; + + //! ionization degree for the electron-electron collisions (tmp is the previous one) + //double m_ionDegree = 0.0; //! Electron energy distribution Difference dF/dε (V^-5/2) Eigen::ArrayXd m_electronEnergyDistDiff; @@ -425,6 +508,10 @@ class PlasmaPhase: public IdealGasPhase void updateElasticElectronEnergyLossCoefficients(); private: + + //! Solver used to calculate the EEDF based on electron collision rates + unique_ptr m_eedfSolver = nullptr; + //! Electron energy distribution change variable. Whenever //! #m_electronEnergyDist changes, this int is incremented. int m_distNum = -1; @@ -442,10 +529,6 @@ class PlasmaPhase: public IdealGasPhase //! The collision-target species indices of #m_collisions vector m_targetSpeciesIndices; - //! Interpolated cross sections. This is used for storing - //! interpolated cross sections temporarily. - vector m_interp_cs; - //! The list of whether the interpolated cross sections is ready vector m_interp_cs_ready; @@ -454,8 +537,7 @@ class PlasmaPhase: public IdealGasPhase void setCollisions(); //! Add a collision and record the target species - void addCollision(std::shared_ptr collision); - + void addCollision(shared_ptr collision); }; } diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 2f145fd72fc..18ba16dba3a 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -232,7 +232,10 @@ class ReactorBase } //! Returns the current internal energy (J/kg) of the reactor's contents. + //! @deprecated To be removed after %Cantera 3.2. double intEnergy_mass() const { + warn_deprecated("ReactorBase::intEnergy_mass", + "To be removed after Cantera 3.2."); return m_intEnergy; } @@ -304,7 +307,10 @@ class ReactorBase double m_vol = 0.0; //!< Current volume of the reactor [m^3] double m_mass = 0.0; //!< Current mass of the reactor [kg] double m_enthalpy = 0.0; //!< Current specific enthalpy of the reactor [J/kg] - double m_intEnergy = 0.0; //!< Current internal energy of the reactor [J/kg] + + //! Current internal energy of the reactor [J/kg] + //! @deprecated To be removed after %Cantera 3.2 + double m_intEnergy = 0.0; double m_pressure = 0.0; //!< Current pressure in the reactor [Pa] vector m_state; vector m_inlet, m_outlet; diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index a22beca57f0..80fdaa4e322 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -193,6 +193,8 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": cdef cppclass CxxPlasmaPhase "Cantera::PlasmaPhase" (CxxThermoPhase): CxxPlasmaPhase() void setElectronTemperature(double) except +translate_exception + void setReducedElectricField(double) except +translate_exception + void setElectricField(double) except +translate_exception void setElectronEnergyLevels(double*, size_t) except +translate_exception void getElectronEnergyLevels(double*) void setDiscretizedElectronEnergyDist(double*, double*, size_t) except +translate_exception @@ -210,6 +212,9 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": size_t nElectronEnergyLevels() double electronPressure() string electronSpeciesName() + double reducedElectricField() + double electricField() + void updateElectronEnergyDistribution() except +translate_exception double elasticPowerLoss() except +translate_exception diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index b091f708829..cdc215b884f 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1783,6 +1783,41 @@ cdef class ThermoPhase(_SolutionBase): raise ThermoModelMethodError(self.thermo_model) return self.plasma.electronPressure() + property reduced_electric_field: + """ + Get/Set the reduced electric field (E/N) [V·m²]. + + .. versionadded:: 3.2 + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.reducedElectricField() + + def __set__(self, value): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setReducedElectricField(value) + + property electric_field: + """ + Get/Set the electric field (E) [V/m]. + + This is the absolute electric field strength. It is related to the reduced + electric field (E/N) through the number density of neutrals. + + .. versionadded:: 3.2 + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.electricField() + + def __set__(self, value): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setElectricField(value) + def set_discretized_electron_energy_distribution(self, levels, distribution): """ Set electron energy distribution. When this method is used, electron @@ -1809,6 +1844,18 @@ cdef class ThermoPhase(_SolutionBase): &data_dist[0], len(levels)) + def update_electron_energy_distribution(self): + """ + Update the electron energy distribution function to account for changes in + composition, temperature, pressure, or electric field strength. + + .. versionadded:: 3.2 + """ + if not self._enable_plasma: + raise TypeError('This method is invalid for ' + f'thermo model: {self.thermo_model}.') + self.plasma.updateElectronEnergyDistribution() + property n_electron_energy_levels: """ Number of electron energy levels """ def __get__(self): diff --git a/samples/python/reactors/nanosecond-pulse-discharge.py b/samples/python/reactors/nanosecond-pulse-discharge.py new file mode 100644 index 00000000000..0296e95591d --- /dev/null +++ b/samples/python/reactors/nanosecond-pulse-discharge.py @@ -0,0 +1,127 @@ +""" +Nanosecond Pulse Plasma Simulation +================================== + +This example simulates a nanosecond-scale pulse discharge in a reactor. +A Gaussian-shaped electric field pulse is applied over a short timescale. +The plasma reaction mechanism used is based on Colin Pavan's mechanism +for methane-air plasmas which is described in his Ph.D. dissertation and the +corresponding AIAA SciTech conference papers: + +C. A. Pavan, "Nanosecond Pulsed Plasmas in Dynamic Combustion Environments," +Ph.D. Thesis, Massachusetts Institute of Technology, 2023. Chapter 5. + +C. A. Pavan and C. Guerra-Garcia, "Modelling the Impact of a Repetitively +Pulsed Nanosecond DBD Plasma on a Mesoscale Flame," in AIAA SCITECH 2022 +Forum, Reston, Virginia: American Institute of Aeronautics and +Astronautics, Jan. 2022, pp. 1-15. DOI: 10.2514/6.2022-0975. + +C. A. Pavan and C. Guerra-Garcia, "Modeling Flame Speed Modification by +Nanosecond Pulsed Discharges to Inform Experimental Design," in AIAA +SCITECH 2023 Forum, Reston, Virginia: American Institute of Aeronautics +and Astronautics, Jan. 2023, pp. 1-15. DOI: 10.2514/6.2023-2056. + +Requires: cantera >= 3.2, matplotlib >= 2.0 + +.. tags:: Python, plasma, reactor network +""" + +import cantera as ct +import numpy as np +import matplotlib.pyplot as plt + +# Gaussian pulse parameters +EN_peak = 190 * 1e-21 # 190 Td +pulse_center = 24e-9 # 24 ns +pulse_width = 3e-9 # standard deviation (3 ns) + +def gaussian_EN(t): + return EN_peak * np.exp(-((t - pulse_center)**2) / (2 * pulse_width**2)) + +# setup +gas = ct.Solution('example_data/methane-plasma-pavan-2023.yaml') +gas.TPX = 300., 101325., 'CH4:0.095, O2:0.19, N2:0.715, e:1E-11' +gas.reduced_electric_field = gaussian_EN(0) +gas.update_electron_energy_distribution() + +r = ct.ConstPressureReactor(gas, energy="off") + +sim = ct.ReactorNet([r]) +sim.verbose = False + +# simulation parameters +t_total = 90e-9 +dt_max = 1e-10 +dt_chunk = 1e-9 # 1 ns chunk +states = ct.SolutionArray(gas, extra=['t']) + +print('{:>10} {:>10} {:>10} {:>14}'.format('t [s]', 'T [K]', 'P [Pa]', 'h [J/kg]')) + +# simulate in 1 ns chunks +t = 0.0 +while t < t_total: + + # integrate over the next chunk + t_end = min(t + dt_chunk, t_total) + while sim.time < t_end: + sim.advance(sim.time + dt_max) #use sim.step + states.append(r.thermo.state, t=sim.time) + print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format( + sim.time, r.T, r.thermo.P, r.thermo.h)) + + EN_t = gaussian_EN(t) + gas.reduced_electric_field = EN_t + gas.update_electron_energy_distribution() + + # reinitialize integrator with new source terms + sim.reinitialize() + + t = t_end + +# Plotting +fig, ax = plt.subplots(2, layout="constrained") + +ax[0].plot(states.t, states.X[:, gas.species_index('e')], label='e') +ax[0].plot(states.t, states.X[:, gas.species_index('O2+')], label='O2+') +ax[0].plot(states.t, states.X[:, gas.species_index('N2+')], label='N2+') +ax[0].plot(states.t, states.X[:, gas.species_index('H2O+')], label='H2O+') +ax[0].plot(states.t, states.X[:, gas.species_index('CH4+')], label='CH4+') +ax[0].plot(states.t, states.X[:, gas.species_index('O')], label='O') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(A)')], label='N2(A)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(B)')], label='N2(B)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(C)')], label='N2(C)') +ax[0].plot(states.t, states.X[:, gas.species_index("N2(a')")], label="N2(a')") +ax[0].plot(states.t, states.X[:, gas.species_index('CH3')], label='CH3', linestyle='--') +ax[0].plot(states.t, states.X[:, gas.species_index('CO2')], label='CO2', linestyle='--') +ax[0].plot(states.t, states.X[:, gas.species_index('CO')], label='CO', linestyle='--') +ax[0].plot(states.t, states.X[:, gas.species_index('H2O')], label='H2O', linestyle='--') +ax[0].plot(states.t, states.X[:, gas.species_index('H')], label='H', linestyle='--') +ax[0].plot(states.t, states.X[:, gas.species_index('OH')], label='OH', linestyle='--') +# N2 vibrational states +""" ax[0].plot(states.t, states.X[:, gas.species_index('N2(v1)')], label='N2(v1)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(v2)')], label='N2(v2)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(v3)')], label='N2(v3)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(v4)')], label='N2(v4)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(v5)')], label='N2(v5)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(v6)')], label='N2(v6)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(v7)')], label='N2(v7)') +ax[0].plot(states.t, states.X[:, gas.species_index('N2(v8)')], label='N2(v8)') """ + +ax[0].set_yscale('log') +ax[0].set_ylim([1e-14, 1e-3]) + +ax[1].plot(states.t, states.T, label='T') +ax2 = ax[1].twinx() +EN_values = [gaussian_EN(t) for t in states.t] +ax2.plot(states.t, EN_values, label='E/N', color='tab:red', linestyle='--') +ax2.set_ylabel('E/N', color='tab:red') +ax2.tick_params(axis='y', labelcolor='tab:red') + +for axx in ax: + axx.legend(loc='lower right', ncol=2) + axx.set_xlabel('Time [s]') + +ax[0].set_ylabel('Mole fraction [-]') +ax[1].set_ylabel('Temperature [K]') + +plt.show() diff --git a/samples/python/reactors/pfr.py b/samples/python/reactors/pfr.py index d59ca382241..2bc07c14123 100644 --- a/samples/python/reactors/pfr.py +++ b/samples/python/reactors/pfr.py @@ -136,7 +136,6 @@ t2[n] = np.sum(t_r2) # write output data states2.append(r2.thermo.state) - ##################################################################### diff --git a/samples/python/thermo/plasma-eedf.py b/samples/python/thermo/plasma-eedf.py new file mode 100644 index 00000000000..370d83e4fd5 --- /dev/null +++ b/samples/python/thermo/plasma-eedf.py @@ -0,0 +1,55 @@ +""" +EEDF calculation +================ +Compute EEDF with two term approximation solver at constant E/N. +Compare with results from BOLOS. + +Air-plasma-Phelps.yaml mechanism file is derived from Phelps cross section data +A compilation of atomic and molecular cross-section data assembled by A. V. Phelps +in the 1970s–1980s for gases such as O₂, N₂, He, Ar, etc. The compilation itself +is unpublished; data used from it are cited as: A. V. Phelps, private communication +(compilation of electron cross-sections), retrieved [8/1/25], from Phelps collection. +(see http://www.lxcat.net/contributors/#d19). + +Requires: cantera >= 3.2, matplotlib >= 2.0 + +.. tags:: Python, plasma +""" + + +import matplotlib.pyplot as plt +import cantera as ct + +gas = ct.Solution('example_data/air-plasma-Phelps.yaml') +gas.TPX = 300., 101325., 'N2:0.79, O2:0.21, N2+:1E-10, Electron:1E-10' +gas.reduced_electric_field = 200.0 * 1e-21 # Reduced electric field [V.m^2] +gas.update_electron_energy_distribution() + +grid = gas.electron_energy_levels +eedf = gas.electron_energy_distribution + +# results from BOLOS +cgrid = [6.000e-02, 6.908e-02, 7.954e-02, 9.158e-02, 1.054e-01, 1.214e-01, 1.398e-01, + 1.609e-01, 1.853e-01, 2.133e-01, 2.456e-01, 2.828e-01, 3.256e-01, 3.749e-01, + 4.317e-01, 4.970e-01, 5.723e-01, 6.589e-01, 7.586e-01, 8.735e-01, 1.006e+00, + 1.158e+00, 1.333e+00, 1.535e+00, 1.767e+00, 2.035e+00, 2.343e+00, 2.698e+00, + 3.106e+00, 3.576e+00, 4.117e+00, 4.741e+00, 5.458e+00, 6.284e+00, 7.236e+00, + 8.331e+00, 9.592e+00, 1.104e+01, 1.272e+01, 1.464e+01, 1.686e+01, 1.941e+01, + 2.235e+01, 2.573e+01, 2.962e+01, 3.411e+01, 3.927e+01, 4.522e+01, 5.206e+01, + 5.994e+01] +cf0 = [1.445e-01, 1.445e-01, 1.445e-01, 1.445e-01, 1.445e-01, 1.445e-01, 1.445e-01, + 1.445e-01, 1.445e-01, 1.444e-01, 1.444e-01, 1.444e-01, 1.443e-01, 1.442e-01, + 1.441e-01, 1.439e-01, 1.436e-01, 1.431e-01, 1.422e-01, 1.408e-01, 1.389e-01, + 1.360e-01, 1.318e-01, 1.256e-01, 1.161e-01, 9.910e-02, 7.723e-02, 6.190e-02, + 5.368e-02, 4.878e-02, 4.461e-02, 4.041e-02, 3.588e-02, 3.094e-02, 2.564e-02, + 2.009e-02, 1.446e-02, 9.423e-03, 5.364e-03, 2.571e-03, 1.085e-03, 3.935e-04, + 1.172e-04, 2.766e-05, 4.955e-06, 6.462e-07, 5.744e-08, 3.272e-09, 1.149e-10, + 4.822e-12] + +fig, ax = plt.subplots() + +ax.loglog(grid, eedf, c='k', label='Cantera') +ax.loglog(cgrid, cf0, ls='None', mfc='None', mec='k', marker='o', label='BOLOS') +ax.set(xlim=(1e-2, 1e2), ylim=(1e-10, 1e4)) +ax.legend() +plt.show() diff --git a/src/base/checkFinite.cpp b/src/base/checkFinite.cpp index cc97e2a0982..f76f7cf6abd 100644 --- a/src/base/checkFinite.cpp +++ b/src/base/checkFinite.cpp @@ -25,7 +25,7 @@ void checkFinite(const double tmp) } } -void checkFinite(const string& name, double* values, size_t N) +void checkFinite(const string& name, const double* values, size_t N) { for (size_t i = 0; i < N; i++) { if (!std::isfinite(values[i])) { diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index 10562353ce1..ba11ce331b1 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -14,6 +14,11 @@ BulkKinetics::BulkKinetics() { bool BulkKinetics::addReaction(shared_ptr r, bool resize) { + shared_ptr rate = r->rate(); + if (rate) { + rate->setContext(*r, *this); + } + bool added = Kinetics::addReaction(r, resize); if (!added) { // undeclared species, etc. @@ -29,7 +34,6 @@ bool BulkKinetics::addReaction(shared_ptr r, bool resize) m_dn.push_back(dn); - shared_ptr rate = r->rate(); string rtype = rate->subType(); if (rtype == "") { rtype = rate->type(); @@ -44,7 +48,6 @@ bool BulkKinetics::addReaction(shared_ptr r, bool resize) // Set index of rate to number of reaction within kinetics rate->setRateIndex(nReactions() - 1); - rate->setContext(*r, *this); // Add reaction rate to evaluator size_t index = m_rateTypes[rtype]; @@ -56,7 +59,6 @@ bool BulkKinetics::addReaction(shared_ptr r, bool resize) } m_concm.push_back(NAN); - m_ready = resize; return true; } diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index 68f7253199d..696fa4f5d09 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -6,6 +6,7 @@ #include "cantera/kinetics/ElectronCollisionPlasmaRate.h" #include "cantera/kinetics/Reaction.h" #include "cantera/kinetics/Kinetics.h" +#include "cantera/thermo/PlasmaPhase.h" #include "cantera/numerics/funcs.h" namespace Cantera @@ -36,13 +37,11 @@ bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase, const Kinetic pp.getElectronEnergyDistribution(distribution.data()); // Update energy levels - levelChanged = pp.levelNumber() != m_level_number; - if (levelChanged) { - m_level_number = pp.levelNumber(); + if (pp.levelNumber() != levelNumber || energyLevels.empty()) { + levelNumber = pp.levelNumber(); energyLevels.resize(pp.nElectronEnergyLevels()); pp.getElectronEnergyLevels(energyLevels.data()); } - return true; } @@ -52,22 +51,29 @@ void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, const UnitSt if (!node.hasKey("energy-levels") && !node.hasKey("cross-sections")) { return; } - if (node.hasKey("energy-levels")) { - m_energyLevels = node["energy-levels"].asVector(); + + if (node.hasKey("kind")) { + m_kind = node["kind"].asString(); } - if (node.hasKey("cross-sections")) { - m_crossSections = node["cross-sections"].asVector(); + if (node.hasKey("target")) { + m_target = node["target"].asString(); } - if (m_energyLevels.size() != m_crossSections.size()) { - throw CanteraError("ElectronCollisionPlasmaRate::setParameters", - "Energy levels and cross section must have the same length."); + if (node.hasKey("product")) { + m_product = node["product"].asString(); } + + m_energyLevels = node["energy-levels"].asVector(); + m_crossSections = node["cross-sections"].asVector(m_energyLevels.size()); + m_threshold = node.getDouble("threshold", 0.0); } void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const { node["type"] = type(); node["energy-levels"] = m_energyLevels; node["cross-sections"] = m_crossSections; + if (!m_kind.empty()) { + node["kind"] = m_kind; + } } void ElectronCollisionPlasmaRate::updateInterpolatedCrossSection( @@ -84,11 +90,20 @@ double ElectronCollisionPlasmaRate::evalFromStruct( const ElectronCollisionPlasmaData& shared_data) { // Interpolate cross-sections data to the energy levels of - // the electron energy distribution function - if (shared_data.levelChanged) { - updateInterpolatedCrossSection(shared_data.energyLevels); + // the electron energy distribution function when the EEDF from the phase changes + if (m_levelNumber != shared_data.levelNumber) { + m_crossSectionsInterpolated.clear(); + for (double level : shared_data.energyLevels) { + m_crossSectionsInterpolated.push_back(linearInterp(level, + m_energyLevels, m_crossSections)); + } + m_levelNumber = shared_data.levelNumber; } + AssertThrowMsg(m_crossSectionsInterpolated.size() == shared_data.distribution.size(), + "ECPR:evalFromStruct", "Size mismatch: len(interp) = {}, len(distrib) = {}", + m_crossSectionsInterpolated.size(), shared_data.distribution.size()); + // Map cross sections to Eigen::ArrayXd auto cs_array = Eigen::Map( m_crossSectionsInterpolated.data(), m_crossSectionsInterpolated.size() @@ -120,7 +135,7 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( // Interpolate cross-sections data to the energy levels of // the electron energy distribution function - if (shared_data.levelChanged) { + if (m_levelNumberSuperelastic != shared_data.levelNumber) { // super elastic collision energy levels and cross-sections vector superElasticEnergyLevels{0.0}; m_crossSectionsOffset.resize(shared_data.energyLevels.size()); @@ -135,6 +150,7 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( superElasticEnergyLevels, m_crossSections); } + m_levelNumberSuperelastic = shared_data.levelNumber; } // Map energyLevels in Eigen::ArrayXd @@ -155,10 +171,11 @@ void ElectronCollisionPlasmaRate::modifyRateConstants( void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) { + const ThermoPhase& thermo = kin.thermo(); // get electron species name string electronName; - if (kin.thermo().type() == "plasma") { - electronName = dynamic_cast(kin.thermo()).electronSpeciesName(); + if (thermo.type() == "plasma") { + electronName = dynamic_cast(thermo).electronSpeciesName(); } else { throw CanteraError("ElectronCollisionPlasmaRate::setContext", "ElectronCollisionPlasmaRate requires plasma phase"); @@ -177,6 +194,37 @@ void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics "ElectronCollisionPlasmaRate requires one and only one electron"); } + // Determine the "kind" of collision if not specified explicitly + if (m_kind.empty()) { + m_kind = "excitation"; // default + if (rxn.reactants == rxn.products) { + m_kind = "effective"; + } else { + for (const auto& [p, stoich] : rxn.products) { + if (p == electronName) { + continue; + } + double q = thermo.charge(thermo.speciesIndex(p)); + if (q > 0) { + m_kind = "ionization"; + } else if (q < 0) { + m_kind = "attachment"; + } + } + } + } + + if (m_threshold == 0.0 && + (m_kind == "excitation" || m_kind == "ionization" || m_kind == "attachment")) + { + for (size_t i = 0; i < m_energyLevels.size(); i++) { + if (m_energyLevels[i] > 0.0) { // Look for first non-zero cross-section + m_threshold = m_energyLevels[i]; + break; + } + } + } + if (!rxn.reversible) { return; // end checking of forward reaction } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 8b0e069d470..3ba5a515d07 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -388,6 +388,11 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base, bool resize) } size_t i = nReactions(); + shared_ptr rate = r_base->rate(); + if (rate) { + rate->setContext(*r_base, *this); + } + bool added = Kinetics::addReaction(r_base, resize); if (!added) { return false; @@ -408,9 +413,7 @@ bool InterfaceKinetics::addReaction(shared_ptr r_base, bool resize) } // Set index of rate to number of reaction within kinetics - shared_ptr rate = r_base->rate(); rate->setRateIndex(nReactions() - 1); - rate->setContext(*r_base, *this); string rtype = rate->subType(); if (rtype == "") { diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index f7287cd62ce..9a646c9f180 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -46,8 +46,6 @@ void Kinetics::resizeReactions() m_stoichMatrix = m_productStoich.stoichCoeffs(); // reactants are destroyed for positive net rate of progress m_stoichMatrix -= m_reactantStoich.stoichCoeffs(); - - m_ready = true; } void Kinetics::checkReactionArraySize(size_t ii) const @@ -700,8 +698,6 @@ bool Kinetics::addReaction(shared_ptr r, bool resize) if (resize) { resizeReactions(); - } else { - m_ready = false; } for (const auto& [id, callback] : m_reactionAddedCallbacks) { diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index cb98c81dc99..753f0d68f52 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -28,11 +28,19 @@ Reaction::Reaction(const Composition& reactants_, const Composition& products_, shared_ptr rate_, shared_ptr tbody_) - : reactants(reactants_) - , products(products_) - , m_from_composition(true) + : m_from_composition(true) , m_third_body(tbody_) { + for (auto& [species, stoich] : reactants_) { + if (stoich != 0.0) { + reactants[species] = stoich; + } + } + for (auto& [species, stoich] : products_) { + if (stoich != 0.0) { + products[species] = stoich; + } + } if (reactants.count("M") || products.count("M")) { throw CanteraError("Reaction::Reaction", "Third body 'M' must not be included in either reactant or product maps."); @@ -56,7 +64,9 @@ Reaction::Reaction(const Composition& reactants_, if (name != "M") { m_third_body->explicit_3rd = true; } - } else if (!tbody_ && third.size() == 1) { + } else if (!tbody_ && third.size() == 1 + && m_rate->type() != "electron-collision-plasma") + { // implicit third body string name = third.begin()->first; m_third_body = make_shared(name); diff --git a/src/numerics/funcs.cpp b/src/numerics/funcs.cpp index 70f9b5f14b0..c8cf46b7c03 100644 --- a/src/numerics/funcs.cpp +++ b/src/numerics/funcs.cpp @@ -12,6 +12,10 @@ namespace Cantera double linearInterp(double x, const vector& xpts, const vector& fpts) { + AssertThrowMsg(!xpts.empty(), "linearInterp", "x data empty"); + AssertThrowMsg(!fpts.empty(), "linearInterp", "f(x) data empty"); + AssertThrowMsg(xpts.size() == fpts.size(), "linearInterp", + "len(xpts) = {}, len(fpts) = {}", xpts.size(), fpts.size()); if (x <= xpts[0]) { return fpts[0]; } diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp new file mode 100644 index 00000000000..6ec29158328 --- /dev/null +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -0,0 +1,600 @@ +/** + * @file EEDFTwoTermApproximation.cpp + * EEDF Two-Term approximation solver. Implementation file for class + * EEDFTwoTermApproximation. + */ + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/thermo/EEDFTwoTermApproximation.h" +#include "cantera/base/ctexceptions.h" +#include "cantera/numerics/eigen_dense.h" +#include "cantera/numerics/funcs.h" +#include "cantera/thermo/PlasmaPhase.h" +#include "cantera/kinetics/ElectronCollisionPlasmaRate.h" + +namespace Cantera +{ + +typedef Eigen::SparseMatrix SparseMat; + +EEDFTwoTermApproximation::EEDFTwoTermApproximation(PlasmaPhase* s) +{ + // store a pointer to s. + m_phase = s; + m_first_call = true; + m_has_EEDF = false; + m_gamma = pow(2.0 * ElectronCharge / ElectronMass, 0.5); +} + +void EEDFTwoTermApproximation::setLinearGrid(double& kTe_max, size_t& ncell) +{ + m_points = ncell; + m_gridCenter.resize(m_points); + m_gridEdge.resize(m_points + 1); + m_f0.resize(m_points); + m_f0_edge.resize(m_points + 1); + for (size_t j = 0; j < m_points; j++) { + m_gridCenter[j] = kTe_max * ( j + 0.5 ) / m_points; + m_gridEdge[j] = kTe_max * j / m_points; + } + m_gridEdge[m_points] = kTe_max; + setGridCache(); +} + +int EEDFTwoTermApproximation::calculateDistributionFunction() +{ + if (m_first_call) { + initSpeciesIndexCrossSections(); + m_first_call = false; + } + + updateMoleFractions(); + updateCrossSections(); + + if (!m_has_EEDF) { + writelog("No existing EEDF. Using first guess method: {}\n", m_firstguess); + if (m_firstguess == "maxwell") { + writelog("First guess EEDF maxwell\n"); + for (size_t j = 0; j < m_points; j++) { + m_f0(j) = 2.0 * pow(1.0 / Pi, 0.5) * pow(m_init_kTe, -3. / 2.) * + exp(-m_gridCenter[j] / m_init_kTe); + } + } else { + throw CanteraError("EEDFTwoTermApproximation::calculateDistributionFunction", + " unknown EEDF first guess"); + } + } + + converge(m_f0); + + // write the EEDF at grid edges + vector f(m_f0.data(), m_f0.data() + m_f0.rows() * m_f0.cols()); + vector x(m_gridCenter.data(), m_gridCenter.data() + m_gridCenter.rows() * m_gridCenter.cols()); + for (size_t i = 0; i < m_points + 1; i++) { + m_f0_edge[i] = linearInterp(m_gridEdge[i], x, f); + } + + m_has_EEDF = true; + + // update electron mobility + m_electronMobility = electronMobility(m_f0); + return 0; +} + +void EEDFTwoTermApproximation::converge(Eigen::VectorXd& f0) +{ + double err0 = 0.0; + double err1 = 0.0; + double delta = m_delta0; + + if (m_maxn == 0) { + throw CanteraError("EEDFTwoTermApproximation::converge", + "m_maxn is zero; no iterations will occur."); + } + if (m_points == 0) { + throw CanteraError("EEDFTwoTermApproximation::converge", + "m_points is zero; the EEDF grid is empty."); + } + if (isnan(delta) || delta == 0.0) { + throw CanteraError("EEDFTwoTermApproximation::converge", + "m_delta0 is NaN or zero; solver cannot update."); + } + + for (size_t n = 0; n < m_maxn; n++) { + if (0.0 < err1 && err1 < err0) { + delta *= log(m_factorM) / (log(err0) - log(err1)); + } + + Eigen::VectorXd f0_old = f0; + f0 = iterate(f0_old, delta); + checkFinite("EEDFTwoTermApproximation::converge: f0", f0.data(), f0.size()); + + err0 = err1; + Eigen::VectorXd Df0 = (f0_old - f0).cwiseAbs(); + err1 = norm(Df0, m_gridCenter); + if (err1 < m_rtol) { + break; + } else if (n == m_maxn - 1) { + throw CanteraError("WeaklyIonizedGas::converge", "Convergence failed"); + } + } +} + +Eigen::VectorXd EEDFTwoTermApproximation::iterate(const Eigen::VectorXd& f0, double delta) +{ + // CQM multiple call to vector_* and matrix_* + // probably extremely ineficient + // must be refactored!! + + SparseMat PQ(m_points, m_points); + vector g = vector_g(f0); + + for (size_t k : m_phase->kInelastic()) { + SparseMat Q_k = matrix_Q(g, k); + SparseMat P_k = matrix_P(g, k); + PQ += (matrix_Q(g, k) - matrix_P(g, k)) * m_X_targets[m_klocTargets[k]]; + } + + SparseMat A = matrix_A(f0); + SparseMat I(m_points, m_points); + for (size_t i = 0; i < m_points; i++) { + I.insert(i,i) = 1.0; + } + A -= PQ; + A *= delta; + A += I; + + // SparseLU : + Eigen::SparseLU solver(A); + if (solver.info() == Eigen::NumericalIssue) { + throw CanteraError("EEDFTwoTermApproximation::iterate", + "Error SparseLU solver: NumericalIssue"); + } else if (solver.info() == Eigen::InvalidInput) { + throw CanteraError("EEDFTwoTermApproximation::iterate", + "Error SparseLU solver: InvalidInput"); + } + if (solver.info() != Eigen::Success) { + throw CanteraError("EEDFTwoTermApproximation::iterate", + "Error SparseLU solver", "Decomposition failed"); + return f0; + } + + // solve f0 + Eigen::VectorXd f1 = solver.solve(f0); + if(solver.info() != Eigen::Success) { + throw CanteraError("EEDFTwoTermApproximation::iterate", "Solving failed"); + return f0; + } + + checkFinite("EEDFTwoTermApproximation::converge: f0", f1.data(), f1.size()); + f1 /= norm(f1, m_gridCenter); + return f1; +} + +double EEDFTwoTermApproximation::integralPQ(double a, double b, double u0, double u1, + double g, double x0) +{ + double A1; + double A2; + if (g != 0.0) { + double expm1a = expm1(g * (-a + x0)); + double expm1b = expm1(g * (-b + x0)); + double ag = a * g; + double ag1 = ag + 1; + double bg = b * g; + double bg1 = bg + 1; + A1 = (expm1a * ag1 + ag - expm1b * bg1 - bg) / (g*g); + A2 = (expm1a * (2 * ag1 + ag * ag) + ag * (ag + 2) - + expm1b * (2 * bg1 + bg * bg) - bg * (bg + 2)) / (g*g*g); + } else { + A1 = 0.5 * (b*b - a*a); + A2 = 1.0 / 3.0 * (b*b*b - a*a*a); + } + + // The interpolation formula of u(x) = c0 + c1 * x + double c0 = (a * u1 - b * u0) / (a - b); + double c1 = (u0 - u1) / (a - b); + + return c0 * A1 + c1 * A2; +} + +vector EEDFTwoTermApproximation::vector_g(const Eigen::VectorXd& f0) +{ + vector g(m_points, 0.0); + const double f_min = 1e-300; // Smallest safe floating-point value + + // Handle first point (i = 0) + double f1 = std::max(f0(1), f_min); + double f0_ = std::max(f0(0), f_min); + g[0] = log(f1 / f0_) / (m_gridCenter[1] - m_gridCenter[0]); + + // Handle last point (i = N) + size_t N = m_points - 1; + double fN = std::max(f0(N), f_min); + double fNm1 = std::max(f0(N - 1), f_min); + g[N] = log(fN / fNm1) / (m_gridCenter[N] - m_gridCenter[N - 1]); + + // Handle interior points + for (size_t i = 1; i < N; ++i) { + double f_up = std::max(f0(i + 1), f_min); + double f_down = std::max(f0(i - 1), f_min); + g[i] = log(f_up / f_down) / (m_gridCenter[i + 1] - m_gridCenter[i - 1]); + } + return g; +} + +SparseMat EEDFTwoTermApproximation::matrix_P(const vector& g, size_t k) +{ + SparseTriplets tripletList; + for (size_t n = 0; n < m_eps[k].size(); n++) { + double eps_a = m_eps[k][n][0]; + double eps_b = m_eps[k][n][1]; + double sigma_a = m_sigma[k][n][0]; + double sigma_b = m_sigma[k][n][1]; + size_t j = m_j[k][n]; + double r = integralPQ(eps_a, eps_b, sigma_a, sigma_b, g[j], m_gridCenter[j]); + double p = m_gamma * r; + + tripletList.emplace_back(j, j, p); + } + SparseMat P(m_points, m_points); + P.setFromTriplets(tripletList.begin(), tripletList.end()); + return P; +} + +SparseMat EEDFTwoTermApproximation::matrix_Q(const vector& g, size_t k) +{ + SparseTriplets tripletList; + for (size_t n = 0; n < m_eps[k].size(); n++) { + double eps_a = m_eps[k][n][0]; + double eps_b = m_eps[k][n][1]; + double sigma_a = m_sigma[k][n][0]; + double sigma_b = m_sigma[k][n][1]; + size_t i = m_i[k][n]; + size_t j = m_j[k][n]; + double r = integralPQ(eps_a, eps_b, sigma_a, sigma_b, g[j], m_gridCenter[j]); + double q = m_inFactor[k] * m_gamma * r; + + tripletList.emplace_back(i, j, q); + } + SparseMat Q(m_points, m_points); + Q.setFromTriplets(tripletList.begin(), tripletList.end()); + return Q; +} + +SparseMat EEDFTwoTermApproximation::matrix_A(const Eigen::VectorXd& f0) +{ + vector a0(m_points + 1); + vector a1(m_points + 1); + size_t N = m_points - 1; + // Scharfetter-Gummel scheme + double nu = netProductionFrequency(f0); + a0[0] = NAN; + a1[0] = NAN; + a0[N+1] = NAN; + a1[N+1] = NAN; + + double nDensity = m_phase->molarDensity() * Avogadro; + double alpha; + double E = m_phase->electricField(); + if (m_growth == "spatial") { + double mu = electronMobility(f0); + double D = electronDiffusivity(f0); + alpha = (mu * E - sqrt(pow(mu * E, 2) - 4 * D * nu * nDensity)) / 2.0 / D / nDensity; + } else { + alpha = 0.0; + } + + double sigma_tilde; + double omega = 2 * Pi * m_phase->electricFieldFrequency(); + for (size_t j = 1; j < m_points; j++) { + if (m_growth == "temporal") { + sigma_tilde = m_totalCrossSectionEdge[j] + nu / pow(m_gridEdge[j], 0.5) / m_gamma; + } else { + sigma_tilde = m_totalCrossSectionEdge[j]; + } + double q = omega / (nDensity * m_gamma * pow(m_gridEdge[j], 0.5)); + double W = -m_gamma * m_gridEdge[j] * m_gridEdge[j] * m_sigmaElastic[j]; + double F = sigma_tilde * sigma_tilde / (sigma_tilde * sigma_tilde + q * q); + double DA = m_gamma / 3.0 * pow(E / nDensity, 2.0) * m_gridEdge[j]; + double DB = m_gamma * m_phase->temperature() * Boltzmann / ElectronCharge * m_gridEdge[j] * m_gridEdge[j] * m_sigmaElastic[j]; + double D = DA / sigma_tilde * F + DB; + if (m_growth == "spatial") { + W -= m_gamma / 3.0 * 2 * alpha * E / nDensity * m_gridEdge[j] / sigma_tilde; + } + + double z = W * (m_gridCenter[j] - m_gridCenter[j-1]) / D; + if (!std::isfinite(z)) { + throw CanteraError("matrix_A", "Non-finite Peclet number encountered"); + } + if (std::abs(z) > 500) { + writelog("Warning: Large Peclet number z = {:.3e} at j = {}. W = {:.3e}, D = {:.3e}, E/N = {:.3e}\n", + z, j, W, D, E / nDensity); + } + a0[j] = W / (1 - std::exp(-z)); + a1[j] = W / (1 - std::exp(z)); + } + + SparseTriplets tripletList; + // center diagonal + // zero flux b.c. at energy = 0 + tripletList.emplace_back(0, 0, a0[1]); + + for (size_t j = 1; j < m_points - 1; j++) { + tripletList.emplace_back(j, j, a0[j+1] - a1[j]); + } + + // upper diagonal + for (size_t j = 0; j < m_points - 1; j++) { + tripletList.emplace_back(j, j+1, a1[j+1]); + } + + // lower diagonal + for (size_t j = 1; j < m_points; j++) { + tripletList.emplace_back(j, j-1, -a0[j]); + } + + // zero flux b.c. + tripletList.emplace_back(N, N, -a1[N]); + + SparseMat A(m_points, m_points); + A.setFromTriplets(tripletList.begin(), tripletList.end()); + + //plus G + SparseMat G(m_points, m_points); + if (m_growth == "temporal") { + for (size_t i = 0; i < m_points; i++) { + G.insert(i, i) = 2.0 / 3.0 * (pow(m_gridEdge[i+1], 1.5) - pow(m_gridEdge[i], 1.5)) * nu; + } + } else if (m_growth == "spatial") { + double nDensity = m_phase->molarDensity() * Avogadro; + for (size_t i = 0; i < m_points; i++) { + double sigma_c = 0.5 * (m_totalCrossSectionEdge[i] + m_totalCrossSectionEdge[i + 1]); + G.insert(i, i) = - alpha * m_gamma / 3 * (alpha * (pow(m_gridEdge[i + 1], 2) - pow(m_gridEdge[i], 2)) / sigma_c / 2 + - E / nDensity * (m_gridEdge[i + 1] / m_totalCrossSectionEdge[i + 1] - m_gridEdge[i] / m_totalCrossSectionEdge[i])); + } + } + return A + G; +} + +double EEDFTwoTermApproximation::netProductionFrequency(const Eigen::VectorXd& f0) +{ + double nu = 0.0; + vector g = vector_g(f0); + + for (size_t k = 0; k < m_phase->nCollisions(); k++) { + if (m_phase->collisionRate(k)->kind() == "ionization" || + m_phase->collisionRate(k)->kind() == "attachment") { + SparseMat PQ = (matrix_Q(g, k) - matrix_P(g, k)) * + m_X_targets[m_klocTargets[k]]; + Eigen::VectorXd s = PQ * f0; + checkFinite("EEDFTwoTermApproximation::netProductionFrequency: s", + s.data(), s.size()); + nu += s.sum(); + } + } + return nu; +} + +double EEDFTwoTermApproximation::electronDiffusivity(const Eigen::VectorXd& f0) +{ + vector y(m_points, 0.0); + double nu = netProductionFrequency(f0); + for (size_t i = 0; i < m_points; i++) { + if (m_gridCenter[i] != 0.0) { + y[i] = m_gridCenter[i] * f0(i) / + (m_totalCrossSectionCenter[i] + nu / m_gamma / pow(m_gridCenter[i], 0.5)); + } + } + double nDensity = m_phase->molarDensity() * Avogadro; + auto f = Eigen::Map(y.data(), y.size()); + auto x = Eigen::Map(m_gridCenter.data(), m_gridCenter.size()); + return 1./3. * m_gamma * simpson(f, x) / nDensity; +} + +double EEDFTwoTermApproximation::electronMobility(const Eigen::VectorXd& f0) +{ + double nu = netProductionFrequency(f0); + vector y(m_points + 1, 0.0); + for (size_t i = 1; i < m_points; i++) { + // calculate df0 at i-1/2 + double df0 = (f0(i) - f0(i-1)) / (m_gridCenter[i] - m_gridCenter[i-1]); + if (m_gridEdge[i] != 0.0) { + y[i] = m_gridEdge[i] * df0 / + (m_totalCrossSectionEdge[i] + nu / m_gamma / pow(m_gridEdge[i], 0.5)); + } + } + double nDensity = m_phase->molarDensity() * Avogadro; + auto f = ConstMappedVector(y.data(), y.size()); + auto x = ConstMappedVector(m_gridEdge.data(), m_gridEdge.size()); + return -1./3. * m_gamma * simpson(f, x) / nDensity; +} + +void EEDFTwoTermApproximation::initSpeciesIndexCrossSections() +{ + // set up target index + m_kTargets.resize(m_phase->nCollisions()); + m_klocTargets.resize(m_phase->nCollisions()); + m_inFactor.resize(m_phase->nCollisions()); + for (size_t k = 0; k < m_phase->nCollisions(); k++) { + m_kTargets[k] = m_phase->targetIndex(k); + // Check if it is a new target or not : + auto it = find(m_k_lg_Targets.begin(), m_k_lg_Targets.end(), m_kTargets[k]); + + if (it == m_k_lg_Targets.end()){ + m_k_lg_Targets.push_back(m_kTargets[k]); + m_klocTargets[k] = m_k_lg_Targets.size() - 1; + } else { + m_klocTargets[k] = distance(m_k_lg_Targets.begin(), it); + } + + const auto& kind = m_phase->collisionRate(k)->kind(); + + if (kind == "ionization") { + m_inFactor[k] = 2; + } else if (kind == "attachment") { + m_inFactor[k] = 0; + } else { + m_inFactor[k] = 1; + } + } + + m_X_targets.resize(m_k_lg_Targets.size()); + m_X_targets_prev.resize(m_k_lg_Targets.size()); + for (size_t k = 0; k < m_X_targets.size(); k++) { + size_t k_glob = m_k_lg_Targets[k]; + m_X_targets[k] = m_phase->moleFraction(k_glob); + m_X_targets_prev[k] = m_phase->moleFraction(k_glob); + } + + // set up indices of species which has no cross-section data + for (size_t k = 0; k < m_phase->nSpecies(); k++) { + auto it = std::find(m_kTargets.begin(), m_kTargets.end(), k); + if (it == m_kTargets.end()) { + m_kOthers.push_back(k); + } + } +} + +void EEDFTwoTermApproximation::updateCrossSections() +{ + // Compute sigma_m and sigma_\epsilon + calculateTotalCrossSection(); + calculateTotalElasticCrossSection(); +} + +// Update the species mole fractions used for EEDF computation +void EEDFTwoTermApproximation::updateMoleFractions() +{ + double tmp_sum = 0.0; + for (size_t k = 0; k < m_X_targets.size(); k++) { + m_X_targets[k] = m_phase->moleFraction(m_k_lg_Targets[k]); + tmp_sum = tmp_sum + m_phase->moleFraction(m_k_lg_Targets[k]); + } + + // Normalize the mole fractions to unity: + for (size_t k = 0; k < m_X_targets.size(); k++) { + m_X_targets[k] = m_X_targets[k] / tmp_sum; + } +} + +void EEDFTwoTermApproximation::calculateTotalCrossSection() +{ + m_totalCrossSectionCenter.assign(m_points, 0.0); + m_totalCrossSectionEdge.assign(m_points + 1, 0.0); + for (size_t k = 0; k < m_phase->nCollisions(); k++) { + auto& x = m_phase->collisionRate(k)->energyLevels(); + auto& y = m_phase->collisionRate(k)->crossSections(); + + for (size_t i = 0; i < m_points; i++) { + m_totalCrossSectionCenter[i] += m_X_targets[m_klocTargets[k]] * + linearInterp(m_gridCenter[i], x, y); + } + for (size_t i = 0; i < m_points + 1; i++) { + m_totalCrossSectionEdge[i] += m_X_targets[m_klocTargets[k]] * + linearInterp(m_gridEdge[i], x, y); + } + } +} + +void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() +{ + m_sigmaElastic.clear(); + m_sigmaElastic.resize(m_points, 0.0); + for (size_t k : m_phase->kElastic()) { + auto& x = m_phase->collisionRate(k)->energyLevels(); + auto& y = m_phase->collisionRate(k)->crossSections(); + // Note: + // moleFraction(m_kTargets[k]) <=> m_X_targets[m_klocTargets[k]] + double mass_ratio = ElectronMass / (m_phase->molecularWeight(m_kTargets[k]) / Avogadro); + for (size_t i = 0; i < m_points; i++) { + m_sigmaElastic[i] += 2.0 * mass_ratio * m_X_targets[m_klocTargets[k]] * + linearInterp(m_gridEdge[i], x, y); + } + } +} + +void EEDFTwoTermApproximation::setGridCache() +{ + m_sigma.clear(); + m_sigma.resize(m_phase->nCollisions()); + m_eps.clear(); + m_eps.resize(m_phase->nCollisions()); + m_j.clear(); + m_j.resize(m_phase->nCollisions()); + m_i.clear(); + m_i.resize(m_phase->nCollisions()); + for (size_t k = 0; k < m_phase->nCollisions(); k++) { + auto& collision = m_phase->collisionRate(k); + auto& x = collision->energyLevels(); + auto& y = collision->crossSections(); + vector eps1(m_points + 1); + int shiftFactor = (collision->kind() == "ionization") ? 2 : 1; + + for (size_t i = 0; i < m_points + 1; i++) { + eps1[i] = clip(shiftFactor * m_gridEdge[i] + collision->threshold(), + m_gridEdge[0] + 1e-9, m_gridEdge[m_points] - 1e-9); + } + vector nodes = eps1; + for (size_t i = 0; i < m_points + 1; i++) { + if (m_gridEdge[i] >= eps1[0] && m_gridEdge[i] <= eps1[m_points]) { + nodes.push_back(m_gridEdge[i]); + } + } + for (size_t i = 0; i < x.size(); i++) { + if (x[i] >= eps1[0] && x[i] <= eps1[m_points]) { + nodes.push_back(x[i]); + } + } + + std::sort(nodes.begin(), nodes.end()); + auto last = std::unique(nodes.begin(), nodes.end()); + nodes.resize(std::distance(nodes.begin(), last)); + vector sigma0(nodes.size()); + for (size_t i = 0; i < nodes.size(); i++) { + sigma0[i] = linearInterp(nodes[i], x, y); + } + + // search position of cell j + for (size_t i = 1; i < nodes.size(); i++) { + auto low = std::lower_bound(m_gridEdge.begin(), m_gridEdge.end(), nodes[i]); + m_j[k].push_back(low - m_gridEdge.begin() - 1); + } + + // search position of cell i + for (size_t i = 1; i < nodes.size(); i++) { + auto low = std::lower_bound(eps1.begin(), eps1.end(), nodes[i]); + m_i[k].push_back(low - eps1.begin() - 1); + } + + // construct sigma + for (size_t i = 0; i < nodes.size() - 1; i++) { + m_sigma[k].push_back({sigma0[i], sigma0[i+1]}); + } + + // construct eps + for (size_t i = 0; i < nodes.size() - 1; i++) { + m_eps[k].push_back({nodes[i], nodes[i+1]}); + } + + // construct sigma_offset + auto x_offset = collision->energyLevels(); + for (auto& element : x_offset) { + element -= collision->threshold(); + } + } +} + +double EEDFTwoTermApproximation::norm(const Eigen::VectorXd& f, const Eigen::VectorXd& grid) +{ + string m_quadratureMethod = "simpson"; + Eigen::VectorXd p(f.size()); + for (int i = 0; i < f.size(); i++) { + p[i] = f(i) * pow(grid[i], 0.5); + } + return numericalQuadrature(m_quadratureMethod, p, grid); +} + +} diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 7bc60e04f9d..9b67d1df8c2 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -4,12 +4,15 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/thermo/PlasmaPhase.h" +#include "cantera/thermo/EEDFTwoTermApproximation.h" #include #include "cantera/thermo/Species.h" #include "cantera/base/global.h" +#include "cantera/numerics/eigen_dense.h" #include "cantera/numerics/funcs.h" #include "cantera/kinetics/Kinetics.h" #include "cantera/kinetics/Reaction.h" +#include #include "cantera/kinetics/ElectronCollisionPlasmaRate.h" namespace Cantera { @@ -22,14 +25,18 @@ PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) { initThermoFile(inputFile, id_); - // initial grid - m_electronEnergyLevels = Eigen::ArrayXd::LinSpaced(m_nPoints, 0.0, 1.0); - // initial electron temperature m_electronTemp = temperature(); - // resize vectors - m_interp_cs.resize(m_nPoints); + // Initialize the Boltzmann Solver + m_eedfSolver = make_unique(this); + + // Set Energy Grid (Hardcoded Defaults for Now) + double kTe_max = 60; + size_t nGridCells = 301; + m_nPoints = nGridCells + 1; + m_eedfSolver->setLinearGrid(kTe_max, nGridCells); + m_electronEnergyLevels = MappedVector(m_eedfSolver->getGridEdge().data(), m_nPoints); } PlasmaPhase::~PlasmaPhase() @@ -51,6 +58,30 @@ void PlasmaPhase::updateElectronEnergyDistribution() "Invalid for discretized electron energy distribution."); } else if (m_distributionType == "isotropic") { setIsotropicElectronEnergyDistribution(); + } else if (m_distributionType == "Boltzmann-two-term") { + auto ierr = m_eedfSolver->calculateDistributionFunction(); + if (ierr == 0) { + auto y = m_eedfSolver->getEEDFEdge(); + m_electronEnergyDist = Eigen::Map(y.data(), m_nPoints); + } else { + throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", + "Call to calculateDistributionFunction failed."); + } + bool validEEDF = ( + m_electronEnergyDist.size() == m_nPoints && + m_electronEnergyDist.allFinite() && + m_electronEnergyDist.maxCoeff() > 0.0 && + m_electronEnergyDist.sum() > 0.0 + ); + + if (validEEDF) { + updateElectronTemperatureFromEnergyDist(); + } else { + writelog("Skipping Te update: EEDF is empty, non-finite, or unnormalized.\n"); + } + } else { + throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution", + "Unknown method '{}' for determining EEDF", m_distributionType); } updateElectronEnergyDistDifference(); electronEnergyDistributionChanged(); @@ -71,7 +102,8 @@ void PlasmaPhase::normalizeElectronEnergyDistribution() { void PlasmaPhase::setElectronEnergyDistributionType(const string& type) { if (type == "discretized" || - type == "isotropic") { + type == "isotropic" || + type == "Boltzmann-two-term") { m_distributionType = type; } else { throw CanteraError("PlasmaPhase::setElectronEnergyDistributionType", @@ -104,22 +136,13 @@ void PlasmaPhase::setMeanElectronEnergy(double energy) { updateElectronEnergyDistribution(); } -void PlasmaPhase::setElectronEnergyLevels(const double* levels, size_t length, - bool updateEnergyDist) +void PlasmaPhase::setElectronEnergyLevels(const double* levels, size_t length) { m_nPoints = length; m_electronEnergyLevels = Eigen::Map(levels, length); checkElectronEnergyLevels(); electronEnergyLevelChanged(); - if (updateEnergyDist) updateElectronEnergyDistribution(); - m_interp_cs.resize(m_nPoints); - // The cross sections are interpolated on the energy levels - if (nCollisions() > 0) { - for (size_t i = 0; i < m_collisions.size(); i++) { - m_interp_cs_ready[i] = false; - updateInterpolatedCrossSection(i); - } - } + updateElectronEnergyDistribution(); } void PlasmaPhase::electronEnergyDistributionChanged() @@ -130,6 +153,16 @@ void PlasmaPhase::electronEnergyDistributionChanged() void PlasmaPhase::electronEnergyLevelChanged() { m_levelNum++; + // Cross sections are interpolated on the energy levels + if (m_collisions.size() > 0) { + vector energyLevels(m_nPoints); + MappedVector(energyLevels.data(), m_nPoints) = m_electronEnergyLevels; + for (shared_ptr collision : m_collisions) { + const auto& rate = boost::polymorphic_pointer_downcast + (collision->rate()); + rate->updateInterpolatedCrossSection(energyLevels); + } + } } void PlasmaPhase::checkElectronEnergyLevels() const @@ -165,7 +198,9 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(const double* levels, size_t length) { m_distributionType = "discretized"; - setElectronEnergyLevels(levels, length, false); + m_nPoints = length; + m_electronEnergyLevels = + Eigen::Map(levels, length); m_electronEnergyDist = Eigen::Map(dist, length); checkElectronEnergyLevels(); @@ -227,30 +262,28 @@ void PlasmaPhase::getParameters(AnyMap& phaseNode) const void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) { IdealGasPhase::setParameters(phaseNode, rootNode); - m_root = rootNode; if (phaseNode.hasKey("electron-energy-distribution")) { const AnyMap eedf = phaseNode["electron-energy-distribution"].as(); m_distributionType = eedf["type"].asString(); if (m_distributionType == "isotropic") { if (eedf.hasKey("shape-factor")) { - m_isotropicShapeFactor = eedf["shape-factor"].asDouble(); + setIsotropicShapeFactor(eedf["shape-factor"].asDouble()); } else { throw CanteraError("PlasmaPhase::setParameters", "isotropic type requires shape-factor key."); } - if (eedf.hasKey("energy-levels")) { - setElectronEnergyLevels(eedf["energy-levels"].asVector().data(), - eedf["energy-levels"].asVector().size(), - false); - } if (eedf.hasKey("mean-electron-energy")) { double energy = eedf.convert("mean-electron-energy", "eV"); - // setMeanElectronEnergy() calls updateElectronEnergyDistribution() setMeanElectronEnergy(energy); } else { throw CanteraError("PlasmaPhase::setParameters", "isotropic type requires electron-temperature key."); } + if (eedf.hasKey("energy-levels")) { + auto levels = eedf["energy-levels"].asVector(); + setElectronEnergyLevels(levels.data(), levels.size()); + } + setIsotropicElectronEnergyDistribution(); } else if (m_distributionType == "discretized") { if (!eedf.hasKey("energy-levels")) { throw CanteraError("PlasmaPhase::setParameters", @@ -263,9 +296,32 @@ void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) if (eedf.hasKey("normalize")) { enableNormalizeElectronEnergyDist(eedf["normalize"].asBool()); } - setDiscretizedElectronEnergyDist(eedf["energy-levels"].asVector().data(), - eedf["distribution"].asVector().data(), - eedf["energy-levels"].asVector().size()); + auto levels = eedf["energy-levels"].asVector(); + auto distribution = eedf["distribution"].asVector(levels.size()); + setDiscretizedElectronEnergyDist(levels.data(), distribution.data(), + levels.size()); + } + } + + if (rootNode.hasKey("electron-collisions")) { + for (const auto& item : rootNode["electron-collisions"].asVector()) { + auto rate = make_shared(item); + Composition reactants, products; + reactants[item["target"].asString()] = 1; + reactants[electronSpeciesName()] = 1; + if (item.hasKey("product")) { + products[item["product"].asString()] = 1; + } else { + products[item["target"].asString()] = 1; + } + products[electronSpeciesName()] = 1; + if (rate->kind() == "ionization") { + products[electronSpeciesName()] += 1; + } else if (rate->kind() == "attachment") { + products[electronSpeciesName()] -= 1; + } + auto R = make_shared(reactants, products, rate); + addCollision(R); } } } @@ -275,9 +331,10 @@ bool PlasmaPhase::addSpecies(shared_ptr spec) bool added = IdealGasPhase::addSpecies(spec); size_t k = m_kk - 1; - if (spec->composition.find("E") != spec->composition.end() && - spec->composition.size() == 1 && - spec->composition["E"] == 1) { + if ((spec->name == "e" || spec->name == "Electron") || + (spec->composition.find("E") != spec->composition.end() && + spec->composition.size() == 1 && + spec->composition["E"] == 1)) { if (m_electronSpeciesIndex == npos) { m_electronSpeciesIndex = k; } else { @@ -292,7 +349,8 @@ bool PlasmaPhase::addSpecies(shared_ptr spec) void PlasmaPhase::initThermo() { IdealGasPhase::initThermo(); - // check electron species + + // Check electron species if (m_electronSpeciesIndex == npos) { throw CanteraError("PlasmaPhase::initThermo", "No electron species found."); @@ -312,20 +370,22 @@ void PlasmaPhase::setSolution(std::weak_ptr soln) { void PlasmaPhase::setCollisions() { - m_collisions.clear(); - m_collisionRates.clear(); - m_targetSpeciesIndices.clear(); - if (shared_ptr soln = m_soln.lock()) { shared_ptr kin = soln->kinetics(); if (!kin) { return; } - // add collision from the initial list of reactions + // add collision from the initial list of reactions. Only add reactions we + // haven't seen before + set existing; + for (auto& R : m_collisions) { + existing.insert(R.get()); + } for (size_t i = 0; i < kin->nReactions(); i++) { - std::shared_ptr R = kin->reaction(i); - if (R->rate()->type() != "electron-collision-plasma") { + shared_ptr R = kin->reaction(i); + if (R->rate()->type() != "electron-collision-plasma" + || existing.count(R.get())) { continue; } addCollision(R); @@ -342,7 +402,7 @@ void PlasmaPhase::setCollisions() } } -void PlasmaPhase::addCollision(std::shared_ptr collision) +void PlasmaPhase::addCollision(shared_ptr collision) { size_t i = nCollisions(); @@ -355,13 +415,19 @@ void PlasmaPhase::addCollision(std::shared_ptr collision) }); // Identify target species for electron-collision reactions + string target; for (const auto& [name, _] : collision->reactants) { // Reactants are expected to be electrons and the target species if (name != electronSpeciesName()) { m_targetSpeciesIndices.emplace_back(speciesIndex(name)); + target = name; break; } } + if (target.empty()) { + throw CanteraError("PlasmaPhase::addCollision", "Error identifying target for" + " collision with equation '{}'", collision->equation()); + } m_collisions.emplace_back(collision); m_collisionRates.emplace_back( @@ -370,6 +436,30 @@ void PlasmaPhase::addCollision(std::shared_ptr collision) // resize parameters m_elasticElectronEnergyLossCoefficients.resize(nCollisions()); + updateInterpolatedCrossSection(i); + + // Set up data used by Boltzmann solver + auto& rate = *m_collisionRates.back(); + string kind = m_collisionRates.back()->kind(); + + if ((kind == "effective" || kind == "elastic")) { + for (size_t k = 0; k < m_collisions.size() - 1; k++) { + if (m_collisions[k]->reactants == collision->reactants && + (m_collisionRates[k]->kind() == "elastic" || + m_collisionRates[k]->kind() == "effective") && !collision->duplicate) + { + throw CanteraError("PlasmaPhase::addCollision", "Phase already contains" + " an effective/elastic cross section for '{}'.", target); + } + } + m_kElastic.push_back(i); + } else { + m_kInelastic.push_back(i); + } + + m_energyLevels.push_back(rate.energyLevels()); + m_crossSections.push_back(rate.crossSections()); + m_eedfSolver->setGridCache(); } bool PlasmaPhase::updateInterpolatedCrossSection(size_t i) diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 3ca537f7fc2..51a8bf3c4ba 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -133,7 +133,11 @@ void Reactor::syncState() m_thermo->saveState(m_state); if (m_energy) { m_enthalpy = m_thermo->enthalpy_mass(); - m_intEnergy = m_thermo->intEnergy_mass(); + try { + m_intEnergy = m_thermo->intEnergy_mass(); + } catch (NotImplementedError&) { + m_intEnergy = NAN; + } } m_pressure = m_thermo->pressure(); m_mass = m_thermo->density() * m_vol; @@ -205,7 +209,11 @@ void Reactor::updateConnected(bool updatePressure) { if (updatePressure) { m_pressure = m_thermo->pressure(); } - m_intEnergy = m_thermo->intEnergy_mass(); + try { + m_intEnergy = m_thermo->intEnergy_mass(); + } catch (NotImplementedError&) { + m_intEnergy = NAN; + } m_thermo->saveState(m_state); // Update the mass flow rate of connected flow devices diff --git a/test/data/air-plasma.yaml b/test/data/air-plasma.yaml new file mode 100644 index 00000000000..4ff6bb09284 --- /dev/null +++ b/test/data/air-plasma.yaml @@ -0,0 +1,183 @@ +description: |- + Incomplete, simplified subset of the Phelps cross sections for O2/N2. + https://www.lxcat.net/Phelps. For testing purposes only. + +units: {length: cm, quantity: molec, activation-energy: K} + +phases: +- name: air-plasma-Phelps + thermo: plasma + kinetics: gas + species: + - nasa_gas.yaml/species: [Electron, O2, N2, O2+, N2+, O, O2-] + state: {T: 300.0, P: 1 atm, X: {O2: 0.21, N2: 0.79}} + electron-energy-distribution: + type: Boltzmann-two-term + energy-levels: [0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 22.0, + 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0] + +reactions: +- equation: O2 + Electron => Electron + Electron + O2+ + type: electron-collision-plasma + energy-levels: [12.06, 13.0, 18.0, 28.0, 38.0, 48.0, 58.0, 68.0, 78.0, 88.0, 100.0, + 150.0, 200.0, 300.0, 500.0, 700.0, 1000.0, 1500.0, 2000.0, 3000.0, 5000.0, 7000.0, + 10000.0] + cross-sections: [0.0, 2.3e-22, 2e-21, 7.4e-21, 1.32e-20, 1.8e-20, 2.1e-20, 2.33e-20, + 2.5e-20, 2.6e-20, 2.7e-20, 2.7e-20, 2.5e-20, 2.17e-20, 1.66e-20, 1.35e-20, 1.04e-20, + 7.6e-21, 6e-21, 4.2e-21, 2.7e-21, 2e-21, 1.4e-21] + threshold: 12.06 +- equation: N2 + Electron => N2+ + 2 Electron + type: electron-collision-plasma + energy-levels: [0.0, 15.6, 16.0, 16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 21.0, + 22.0, 23.0, 25.0, 30.0, 34.0, 45.0, 60.0, 75.0, 100.0, 150.0, 200.0] + cross-sections: [0.0, 0.0, 1.95e-22, 4.28e-22, 6.6e-22, 9.11e-22, 1.2e-21, 1.516e-21, + 1.841e-21, 2.13e-21, 2.502e-21, 3.181e-21, 3.869e-21, 4.557e-21, 5.924e-21, + 9.579e-21, 1.1718e-20, 1.6461e-20, 2.0181e-20, 2.2134e-20, 2.3436e-20, 2.2692e-20, + 2.1018e-20] + threshold: 15.6 +- equation: O2 + Electron => O2- + type: electron-collision-plasma + energy-levels: [0.0, 0.058, 0.073, 0.083, 0.089, 0.095, 0.103, 0.109, 0.15, 0.17, 0.2, + 0.21, 0.23, 0.32, 0.33, 0.35, 0.44, 0.45, 0.47, 0.56, 0.57, 0.59, 0.68, 0.69, 0.71, + 0.79, 0.8, 0.82, 0.9, 0.91, 0.93, 1.02, 1.03, 1.05, 1.5, 100.0] + cross-sections: [0.0, 0.0, 5.60795e-41, 1.80256e-40, 4.20596e-41, 8.41193e-41, + 1.80256e-40, 0.0, 0.0, 0.0, 0.0, 3.56506e-41, 0.0, 0.0, 2.30327e-41, 0.0, 0.0, + 1.45206e-41, 0.0, 0.0, 1.10156e-41, 0.0, 0.0, 8.01136e-42, 0.0, 0.0, 7.00994e-42, + 0.0, 0.0, 5.50781e-42, 0.0, 0.0, 4.20596e-42, 0.0, 0.0, 0.0] + threshold: 0.0 + +electron-collisions: +- target: N2 + energy-levels: [0.0, 0.015, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.7, 1.2, 1.5, 1.9, + 2.2, 2.8, 3.3, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 75.0, 150.0] + cross-sections: [1.1e-20, 2.55e-20, 3.4e-20, 4.33e-20, 5.95e-20, 7.1e-20, 7.9e-20, + 9e-20, 9.7e-20, 1e-19, 1.04e-19, 1.2e-19, 1.96e-19, 2.85e-19, 2.8e-19, 1.72e-19, + 1.26e-19, 1.09e-19, 1.01e-19, 1.04e-19, 1.1e-19, 1.02e-19, 9e-20, 6.6e-20, 4.9e-20] + kind: effective +- target: N2 + product: N2(rot) + energy-levels: [0.02, 0.03, 0.4, 0.8, 1.2, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, + 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.6, 5.0] + cross-sections: [0.0, 2.5e-22, 2.5e-22, 2.5e-22, 4.7e-22, 8.6e-22, 1.5e-21, 2.35e-21, + 1.08e-20, 1.9e-20, 2.03e-20, 2.77e-20, 2.5e-20, 2.19e-20, 2.4e-20, 2.17e-20, + 1.62e-20, 1.38e-20, 1.18e-20, 1.03e-20, 8.4e-21, 6.9e-21, 5e-21, 1.7e-21, 0.0] + threshold: 0.02 + kind: excitation +- target: N2 + product: N2(v1) + energy-levels: [0.29, 0.3, 0.33, 0.4, 0.75, 0.9, 1.0, 1.1, 1.16, 1.2, 1.22, 1.4, 1.5, + 1.6, 1.65, 3.6, 4.0, 5.0, 15.0, 18.0, 20.0, 22.0, 23.0, 25.0, 29.0, 32.0, 50.0, + 80.0] + cross-sections: [0.0, 1e-23, 1.7e-23, 2.5e-23, 3.7e-23, 5.5e-23, 6.5e-23, 9e-23, + 1.1e-22, 1.25e-22, 1.35e-22, 7e-22, 1e-21, 1.5e-21, 0.0, 0.0, 5.5e-22, 3.5e-22, + 3.5e-22, 4e-22, 6.5e-22, 8.5e-22, 8.5e-22, 6e-22, 3e-22, 1.5e-22, 1.2e-22, 0.0] + threshold: 0.29 + kind: excitation +- target: N2 + product: N2(v1res) + energy-levels: [0.0, 0.291, 1.6, 1.65, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, + 2.6, 2.7, 2.75, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 4.0, 100.0] + cross-sections: [0.0, 0.0, 0.0, 2.7e-21, 3.15e-21, 5.4e-21, 1.485e-20, 4.8e-20, + 2.565e-20, 1.2e-20, 4.5e-20, 2.76e-20, 1.59e-20, 3.15e-20, 1.545e-20, 6e-21, + 1.35e-20, 5.25e-21, 8.7e-21, 1.17e-20, 8.55e-21, 6.6e-21, 6e-21, 5.85e-21, 5.7e-21, + 0.0, 0.0] + threshold: 0.291 + kind: excitation +- target: N2 + product: N2(v2) + energy-levels: [0.0, 0.59, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, + 2.75, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 100.0] + cross-sections: [0.0, 0.0, 0.0, 1.5e-22, 6.3e-21, 1.935e-20, 3.3e-20, 1.47e-20, + 5.4e-21, 2.115e-20, 3e-20, 5.4e-21, 1.05e-20, 1.725e-20, 1.275e-20, 3.3e-21, 9e-21, + 6.45e-21, 3.75e-21, 3.45e-21, 3e-21, 2.13e-21, 0.0, 0.0] + threshold: 0.59 + kind: excitation +- target: N2 + product: N2(v3) + energy-levels: [0.0, 0.88, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.75, 2.8, + 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 100.0] + cross-sections: [0.0, 0.0, 0.0, 9.6e-21, 2.055e-20, 2.7e-20, 1.695e-20, 7.5e-22, + 9.6e-21, 1.47e-20, 4.5e-21, 9.6e-21, 5.4e-21, 8.55e-21, 4.05e-21, 2.82e-21, + 2.91e-21, 6.15e-22, 0.0, 0.0] + threshold: 0.88 + kind: excitation +- target: N2 + product: N2(C3) + energy-levels: [11.03, 11.5, 12.0, 12.5, 13.0, 13.5, 13.8, 14.0, 14.2, 14.5, 15.0, + 16.0, 17.0, 18.0, 19.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 36.0, 40.0, 50.0, 70.0, + 100.0, 150.0] + cross-sections: [0.0, 2.7e-22, 6.2e-22, 1.31e-21, 2.9e-21, 4.9e-21, 6.2e-21, 6.5e-21, + 6.4e-21, 6.3e-21, 5.5e-21, 4.3e-21, 3.5e-21, 3e-21, 2.7e-21, 2.5e-21, 2.1e-21, + 1.77e-21, 1.5e-21, 1.28e-21, 1.11e-21, 7.8e-22, 6.3e-22, 3.9e-22, 1.5e-22, 1.5e-23, + 0.0] + threshold: 11.03 + kind: excitation +- target: N2 + product: N2(B3) + energy-levels: [0.0, 7.35, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, + 18.0, 20.0, 22.0, 26.0, 30.0, 34.0, 40.0, 50.0, 70.0, 150.0] + cross-sections: [0.0, 0.0, 3.62e-22, 9.38e-22, 1.508e-21, 1.863e-21, 2.003e-21, + 1.99e-21, 1.816e-21, 1.615e-21, 1.447e-21, 1.307e-21, 1.199e-21, 1.112e-21, + 9.51e-22, 8.04e-22, 6.77e-22, 5.63e-22, 4.29e-22, 2.68e-22, 6.7e-23, 0.0] + threshold: 7.35 + kind: excitation +- target: N2 + product: N2(a1) + energy-levels: [0.0, 8.55, 9.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 24.0, 26.0, + 30.0, 40.0, 50.0, 70.0, 100.0, 150.0, 200.0] + cross-sections: [0.0, 0.0, 1.27e-22, 1.474e-21, 1.715e-21, 1.916e-21, 2.023e-21, + 1.99e-21, 1.923e-21, 1.849e-21, 1.621e-21, 1.528e-21, 1.367e-21, 1.065e-21, + 8.51e-22, 6.03e-22, 4.02e-22, 2.68e-22, 2.01e-22] + threshold: 8.55 + kind: excitation +- target: N2 + product: N2(w1) + energy-levels: [0.0, 8.89, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, + 20.0, 22.0, 30.0, 38.0, 50.0, 150.0] + cross-sections: [0.0, 0.0, 1.3e-23, 2.61e-22, 4.76e-22, 6.63e-22, 7.84e-22, 7.71e-22, + 6.7e-22, 5.43e-22, 4.42e-22, 3.75e-22, 2.88e-22, 2.41e-22, 1.54e-22, 9.4e-23, + 4.7e-23, 0.0] + threshold: 8.89 + kind: excitation +- target: O2 + energy-levels: [4.4, 4.9, 5.38, 5.86, 6.1, 6.48, 6.77, 7.05, 7.3, 7.53, 7.77, 8.0, + 8.25, 8.73, 9.2, 9.68, 10.15, 11.35, 100.0] + cross-sections: [0.0, 0.0, 2.3e-23, 7.2e-23, 1.08e-22, 1.38e-22, 1.52e-22, 1.56e-22, + 1.48e-22, 1.31e-22, 1.1e-22, 8.4e-23, 5.4e-23, 2.8e-23, 1.4e-23, 8e-24, 8e-24, + 8e-24, 0.0] + threshold: 0.0 + product: O^-+O + kind: attachment +- target: O2 + energy-levels: [0.0, 0.015, 0.03, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.7, 1.2, 1.5, 1.9, + 2.2, 2.8, 3.3, 4.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 75.0, 150.0, 200.0] + cross-sections: [3.5e-21, 8.7e-21, 1.24e-20, 1.6e-20, 2.5e-20, 3.1e-20, 3.6e-20, + 4.5e-20, 5.2e-20, 6.1e-20, 7.9e-20, 7.6e-20, 6.9e-20, 6.5e-20, 5.8e-20, 5.5e-20, + 5.5e-20, 5.6e-20, 6.6e-20, 8e-20, 8.8e-20, 8.6e-20, 8e-20, 6.8e-20, 6.7e-20, 6e-20] + kind: effective +- target: O2 + product: O2(rot) + energy-levels: [0.07, 0.08, 0.1, 0.2, 0.21, 0.22, 0.32, 0.33, 0.35, 0.44, 0.45, 0.47, + 0.56, 0.57, 0.59, 0.68, 0.69, 0.71, 0.79, 0.8, 0.81, 0.9, 0.91, 0.93, 1.02, 1.03, + 1.05, 1.13, 1.14, 1.16, 1.22, 1.23, 1.26, 1.34, 1.35, 1.37, 1.44, 1.45, 1.47, 1.54, + 1.55, 1.57, 1.64, 1.65, 1.67] + cross-sections: [0.0, 5.4e-23, 0.0, 0.0, 2.16e-22, 0.0, 0.0, 3.84e-22, 0.0, 0.0, + 5.4e-22, 0.0, 0.0, 6.72e-22, 0.0, 0.0, 8.04e-22, 0.0, 0.0, 9.36e-22, 0.0, 0.0, + 8.4e-22, 0.0, 0.0, 7.2e-22, 0.0, 0.0, 4.68e-22, 0.0, 0.0, 6e-22, 0.0, 0.0, 3.6e-22, + 0.0, 0.0, 2.4e-22, 0.0, 0.0, 1.2e-22, 0.0, 0.0, 4.8e-23, 0.0] + threshold: 0.02 + kind: excitation +- target: O2 + product: O2(a1) + energy-levels: [0.977, 1.5, 3.5, 5.62, 6.53, 7.89, 13.0, 20.5, 41.0, 100.0] + cross-sections: [0.0, 5.8e-23, 4.9e-22, 8.25e-22, 9.08e-22, 8.63e-22, 5.27e-22, + 3.24e-22, 1.37e-22, 0.0] + threshold: 0.977 + kind: excitation +- target: O2 + product: O2(b1) + energy-levels: [1.627, 3.0, 4.0, 7.34, 9.26, 13.0, 17.0, 20.7, 24.0, 35.1, 45.1, 100] + cross-sections: [0.0, 9.7e-23, 1.49e-22, 1.91e-22, 1.74e-22, 1.3e-22, 1.3e-22, + 1.25e-22, 1e-22, 6.3e-23, 5e-24, 0.0] + threshold: 1.627 + kind: excitation diff --git a/test/data/consistency-cases.yaml b/test/data/consistency-cases.yaml index 94969b53f09..752318669cb 100644 --- a/test/data/consistency-cases.yaml +++ b/test/data/consistency-cases.yaml @@ -147,6 +147,12 @@ plasma: file: oxygen-plasma.yaml phase: discretized-electron-energy-plasma known-failures: + g_eq_h_minus_Ts/.+: Test does not account for distinct electron temperature + u_eq_sum_uk_Xk/.+: Test does not account for distinct electron temperature + s_eq_sum_sk_Xk/.+: Test does not account for distinct electron temperature + cp_eq_dhdT/(0|2): Test does not account for distinct electron temperature + cv_eq_dudT/(0|2): Test does not account for distinct electron temperature + c_eq_sqrt_dP_drho_const_s/1: Test does not account for distinct electron temperature cp_eq_.+/1: Test does not account for distinct electron temperature cv_eq_.+/1: Test does not account for distinct electron temperature c._eq_dsdT_const_._times_T: Test does not account for distinct electron temperature diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 1d880285612..da199ff2fe6 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1197,7 +1197,8 @@ def collision_data(self): "type": "electron-collision-plasma", "energy-levels": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], "cross-sections": [0.0, 3.83e-20, 4.47e-20, 4.79e-20, 5.07e-20, 5.31e-20, - 5.49e-20, 5.64e-20, 5.77e-20, 5.87e-20, 5.97e-20] + 5.49e-20, 5.64e-20, 5.77e-20, 5.87e-20, 5.97e-20], + "duplicate": True, } def test_converting_electron_energy_to_temperature(self, phase): @@ -1266,6 +1267,8 @@ def test_elastic_power_loss_replace_rate(self, phase): assert phase.elastic_power_loss == approx(11765800095) def test_elastic_power_loss_add_reaction(self, phase): + phase2 = ct.Solution(thermo="plasma", kinetics="bulk", + species=phase.species(), reactions=[]) phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" phase.add_reaction(ct.Reaction.from_dict(self.collision_data, phase)) assert phase.elastic_power_loss == approx(18612132428) @@ -1297,6 +1300,42 @@ def test_elastic_power_loss_change_shape_factor(self, phase): phase.isotropic_shape_factor = 1.1 assert phase.elastic_power_loss == approx(7408711810) + def test_eedf_solver(self): + + phase = ct.Solution('air-plasma.yaml') + phase.TPX = 300., 101325., 'N2:0.79, O2:0.21, N2+:1E-10, Electron:1E-10' + phase.reduced_electric_field = 200.0 * 1e-21 # Reduced electric field [V.m^2] + phase.update_electron_energy_distribution() + + grid = phase.electron_energy_levels + eedf = phase.electron_energy_distribution + + reference_grid = np.logspace(-1, np.log10(60)) + + reference_eedf = np.array([ + 9.1027381e-02, 9.1026393e-02, 9.1025267e-02, 9.1023985e-02, 9.1022523e-02, + 9.1020858e-02, 9.1015025e-02, 9.1006713e-02, 9.0997242e-02, 9.0986450e-02, + 9.0974154e-02, 9.0954654e-02, 9.0923885e-02, 9.0888824e-02, 9.0842837e-02, + 9.0775447e-02, 9.0695937e-02, 9.0578309e-02, 9.0398980e-02, 9.0118320e-02, + 8.9293838e-02, 8.7498617e-02, 8.3767419e-02, 7.5765714e-02, 6.4856820e-02, + 5.5592157e-02, 4.9309310e-02, 4.5268611e-02, 4.2261381e-02, 3.9440745e-02, + 3.6437762e-02, 3.3181527e-02, 2.9616717e-02, 2.5795007e-02, 2.1676205e-02, + 1.7347058e-02, 1.3022044e-02, 8.9705614e-03, 5.5251937e-03, 3.1894295e-03, + 1.7301525e-03, 8.4647152e-04, 3.6030983e-04, 1.2894755e-04, 3.7416645e-05, + 8.4693678e-06, 1.4299900e-06, 1.7026957e-07, 1.3992350e-08, 1.5340110e-09 + ]) + + interp = np.interp(reference_grid, grid, eedf, left=0.0, right=0.0) + + mask = reference_eedf > 1e-8 + rel_error = np.abs(interp[mask] - reference_eedf[mask]) / reference_eedf[mask] + + assert max(rel_error) < 0.01 + + l2_norm = np.linalg.norm(interp - reference_eedf) + assert l2_norm < 1e-3 + + class TestImport: """ Tests the various ways of creating a Solution object