[1D] Prevent two-point control from finding solutions with reverse flow

speth · BangShiuh · commit 4e7d6e29fd1e · 2024-12-11T21:10:18.000-05:00
Fixes #1787
diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h
@@ -89,6 +89,7 @@ struct ElectronCollisionPlasmaData : public ReactionData
 class ElectronCollisionPlasmaRate : public ReactionRate
 {
 public:
+
     ElectronCollisionPlasmaRate() = default;
 
     ElectronCollisionPlasmaRate(const AnyMap& node,
@@ -135,6 +136,16 @@ class ElectronCollisionPlasmaRate : public ReactionRate
         return m_crossSections;
     }
 
+    //! The value of #m_crossSectionsInterpolated [m2]
+    const vector<double>& crossSectionInterpolated() const {
+        return m_crossSectionsInterpolated;
+    }
+
+    //! Set the value of #m_crossSectionsInterpolated [m2]
+    void setCrossSectionInterpolated(vector<double>& cs) {
+        m_crossSectionsInterpolated = cs;
+    }
+
 private:
     //! electron energy levels [eV]
     vector<double> m_energyLevels;
diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h
@@ -18,6 +18,7 @@ namespace Cantera
 
 class Solution;
 class Species;
+class Kinetics;
 
 //! Class Phase is the base class for phases of matter, managing the species and
 //! elements in a phase, as well as the independent variables of temperature,
diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h
@@ -15,6 +15,9 @@
 namespace Cantera
 {
 
+class Reaction;
+class ElectronCollisionPlasmaRate;
+
 /**
  * Base class for a phase with plasma properties. This class manages the
  * plasma properties such as electron energy distribution function (EEDF).
@@ -99,13 +102,6 @@ class PlasmaPhase: public IdealGasPhase
                                           const double* distrb,
                                           size_t length);
 
-    //! Set discretized electron energy distribution.
-    //! @param  distrb The vector of electron energy distribution.
-    //!                Length: #m_nPoints.
-    //! @param  length The length of the vectors, which equals #m_nPoints.
-    void setDiscretizedElectronEnergyDist(const double* distrb,
-                                          size_t length);
-
     //! Get electron energy distribution.
     //! @param  distrb The vector of electron energy distribution.
     //!                Length: #m_nPoints.
@@ -198,6 +194,11 @@ class PlasmaPhase: public IdealGasPhase
         return m_nPoints;
     }
 
+    //! Number of collisions
+    size_t nCollisions() const {
+        return m_collisions.size();
+    }
+
     //! Electron Species Index
     size_t electronSpeciesIndex() const {
         return m_electronSpeciesIndex;
@@ -271,6 +272,25 @@ class PlasmaPhase: public IdealGasPhase
         return m_levelNum;
     }
 
+    //! Return the interpolated cross section Number #m_csNum
+    int interpCrossSectionNum() const {
+        return m_csNum;
+    }
+
+    //! add Solution object
+    virtual void addSolution(std::weak_ptr<Solution> soln) override;
+
+    /**
+     * The elastic power loss(J/s/m3)
+     *   @f[
+     *     P_k = N_A N_A C_e e \sum_k C_k K_k,
+     *   @f]
+     * where @f$ C_k @f$ and @f$ C_e @f$ are the concentration (kmol/m3) of the
+     * target species and electrons, respectively. @f$ K_k @f$ is the elastic
+     * electron energy loss coefficient (eV-m3/s).
+     */
+    double elasticPowerLoss();
+
 protected:
     void updateThermo() const override;
 
@@ -284,6 +304,10 @@ class PlasmaPhase: public IdealGasPhase
     //! the new level.
     void electronEnergyLevelChanged();
 
+    //! When the interpolated cross sections changed, plasma properties such as
+    //! elastic power loss need to be re-evaluated.
+    void interpolatedCrossSectionsChanged();
+
     //! Check the electron energy levels
     /*!
      *  The values of electron energy levels need to be positive and
@@ -317,6 +341,12 @@ class PlasmaPhase: public IdealGasPhase
     //! Electron energy distribution norm
     void normalizeElectronEnergyDistribution();
 
+    //! Update interpolated cross sections
+    void updateInterpolatedCrossSections();
+
+    //! Update electron energy distribution difference
+    void updateElectronEnergyDistDifference();
+
     // Electron energy order in the exponential term
     double m_isotropicShapeFactor = 2.0;
 
@@ -345,6 +375,28 @@ class PlasmaPhase: public IdealGasPhase
     //! Flag of normalizing electron energy distribution
     bool m_do_normalizeElectronEnergyDist = true;
 
+    //! Data for initiate reaction
+    AnyMap m_root;
+
+    //! Electron energy distribution Difference dF/dε (V^-5/2)
+    Eigen::ArrayXd m_electronEnergyDistDiff;
+
+    //! Elastic electron energy loss coefficients (eV m3/s)
+    /*! The elastic electron energy loss coefficient for species k is,
+     *   @f[
+     *     K_k = \frac{2 m_e}{m_k} \sqrt{\frac{2 e}{m_e}} \int_0^{\infty} \sigma_k
+     *           \epsilon^2 \left( F_0 + \frac{k_B T}{e}
+     *           \frac{\partial F_0}{\partial \epsilon} \right) d \epsilon,
+     *   @f]
+     * where @f$ m_e @f$ [kg] is the electron mass, @f$ \epsilon @f$ [V] is the
+     * electron energy, @f$ \sigma_k @f$ [m2] is the reaction collision cross section,
+     * @f$ F_0 @f$ [V^(-3/2)] is the normalized electron energy distribution function.
+     */
+    vector<double> m_elasticElectronEnergyLossCoefficients;
+
+    //! update elastic electron energy loss coefficients
+    void updateElasticElectronEnergyLossCoefficients();
+
 private:
     //! Electron energy distribution change variable. Whenever
     //! #m_electronEnergyDist changes, this int is incremented.
@@ -353,6 +405,28 @@ class PlasmaPhase: public IdealGasPhase
     //! Electron energy level change variable. Whenever
     //! #m_electronEnergyLevels changes, this int is incremented.
     int m_levelNum = -1;
+
+    //! Interpolated collision cross section change variable. Whenever
+    //! #ElectronCollisionPlasmaRate::m_crossSectionsInterpolated changes,
+    //! this int is incremented.
+    int m_csNum = -1;
+
+    //! The list of shared pointers of plasma collision reactions
+    vector<shared_ptr<Reaction>> m_collisions;
+
+    //! The list of shared pointers of ElectronCollisionPlasmaRate
+    vector<shared_ptr<ElectronCollisionPlasmaRate>> m_collisionRates;
+
+    //! The collision-target species indices of #m_collisions
+    vector<size_t> m_targetSpeciesIndices;
+
+    //! Interpolated cross sections. This is used for storing
+    //! interpolated cross sections temporarily.
+    vector<double> m_interp_cs;
+
+    //! Set collisions
+    void setCollisions();
+
 };
 
 }
diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h
@@ -15,6 +15,7 @@
 #include "MultiSpeciesThermo.h"
 #include "cantera/base/Units.h"
 #include "cantera/base/AnyMap.h"
+#include "cantera/base/Solution.h"
 
 namespace Cantera
 {
@@ -1957,6 +1958,11 @@ class ThermoPhase : public Phase
 
     //! @}
 
+    //! add Solution object
+    virtual void addSolution(std::weak_ptr<Solution> soln) {
+        m_soln = soln;
+    }
+
 protected:
     //! Store the parameters of a ThermoPhase object such that an identical
     //! one could be reconstructed using the newThermo(AnyMap&) function. This
@@ -1993,6 +1999,9 @@ class ThermoPhase : public Phase
 
     //! last value of the temperature processed by reference state
     mutable double m_tlast = 0.0;
+
+    //! reference to Solution
+    std::weak_ptr<Solution> m_soln;
 };
 
 }
diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd
@@ -210,6 +210,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h":
         size_t nElectronEnergyLevels()
         double electronPressure()
         string electronSpeciesName()
+        double elasticPowerLoss() except +translate_exception
 
 
 cdef extern from "cantera/cython/thermo_utils.h":
diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx
@@ -1904,6 +1904,12 @@ cdef class ThermoPhase(_SolutionBase):
                 raise ThermoModelMethodError(self.thermo_model)
             return pystr(self.plasma.electronSpeciesName())
 
+    property elastic_power_loss:
+        """ Elastic power loss (J/s/m3) """
+        def __get__(self):
+            if not self._enable_plasma:
+                raise ThermoModelMethodError(self.thermo_model)
+            return self.plasma.elasticPowerLoss()
 
 cdef class InterfacePhase(ThermoPhase):
     """ A class representing a surface, edge phase """
diff --git a/src/base/Solution.cpp b/src/base/Solution.cpp
@@ -42,6 +42,7 @@ void Solution::setName(const string& name) {
 
 void Solution::setThermo(shared_ptr<ThermoPhase> thermo) {
     m_thermo = thermo;
+    m_thermo->addSolution(weak_from_this());
     for (const auto& [id, callback] : m_changeCallbacks) {
         callback();
     }
diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp
@@ -76,10 +76,19 @@ double ElectronCollisionPlasmaRate::evalFromStruct(
     // Interpolate cross-sections data to the energy levels of
     // the electron energy distribution function
     if (shared_data.levelChanged) {
-        m_crossSectionsInterpolated.resize(0);
+        m_crossSectionsInterpolated.clear();
+        m_crossSectionsInterpolated.reserve(shared_data.energyLevels.size());
         for (double level : shared_data.energyLevels) {
-            m_crossSectionsInterpolated.push_back(linearInterp(level,
-                                                  m_energyLevels, m_crossSections));
+            m_crossSectionsInterpolated.emplace_back(linearInterp(level,
+                                                m_energyLevels, m_crossSections));
+        }
+    } else {
+        // The interpolated cross section is set by the PlasmaPhase
+        // check the size of the vector
+        if (m_crossSectionsInterpolated.size() != shared_data.energyLevels.size()) {
+            throw CanteraError("ElectronCollisionPlasmaRate::evalFromStruct",
+                               "Energy levels and interpolated cross section ",
+                               "must have the same length.");
         }
     }
 
diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp
@@ -1258,6 +1258,9 @@ void Flow1D::enableTwoPointControl(bool twoPointControl)
 {
     if (isStrained()) {
         m_twoPointControl = twoPointControl;
+        // Prevent finding spurious solutions with negative velocity (outflow) at either
+        // inlet.
+        setBounds(c_offset_V, -1e-5, 1e20);
     } else {
         throw CanteraError("Flow1D::enableTwoPointControl",
             "Invalid operation: two-point control can only be used"
diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp

Original file line number	Diff line number	Diff line change
`@@ -42,6 +42,7 @@ void Solution::setName(const string& name) {`
`42`	`42`
`43`	`43`	`void Solution::setThermo(shared_ptr<ThermoPhase> thermo) {`
`44`	`44`	`m_thermo = thermo;`
	`45`	`+ m_thermo->addSolution(weak_from_this());`
`45`	`46`	`for (const auto& [id, callback] : m_changeCallbacks) {`
`46`	`47`	`callback();`
`47`	`48`	`}`