Cantera  3.1.0a1
PlasmaPhase.cpp
Go to the documentation of this file.
1 //! @file PlasmaPhase.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
7 #include <boost/math/special_functions/gamma.hpp>
9 #include "cantera/base/global.h"
10 #include "cantera/numerics/funcs.h"
11 
12 namespace Cantera {
13 
14 PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_)
15 {
16  initThermoFile(inputFile, id_);
17 
18  // initial grid
19  m_electronEnergyLevels = Eigen::ArrayXd::LinSpaced(m_nPoints, 0.0, 1.0);
20 
21  // initial electron temperature
23 }
24 
26 {
27  if (m_distributionType == "discretized") {
28  throw CanteraError("PlasmaPhase::updateElectronEnergyDistribution",
29  "Invalid for discretized electron energy distribution.");
30  } else if (m_distributionType == "isotropic") {
32  }
33 }
34 
36  Eigen::ArrayXd eps32 = m_electronEnergyLevels.pow(3./2.);
37  double norm = 2./3. * numericalQuadrature(m_quadratureMethod,
38  m_electronEnergyDist, eps32);
39  if (norm < 0.0) {
40  throw CanteraError("PlasmaPhase::normalizeElectronEnergyDistribution",
41  "The norm is negative. This might be caused by bad "
42  "electron energy distribution");
43  }
44  m_electronEnergyDist /= norm;
45 }
46 
48 {
49  if (type == "discretized" ||
50  type == "isotropic") {
52  } else {
53  throw CanteraError("PlasmaPhase::setElectronEnergyDistributionType",
54  "Unknown type for electron energy distribution.");
55  }
56 }
57 
59 {
61  double x = m_isotropicShapeFactor;
62  double gamma1 = boost::math::tgamma(3.0 / 2.0 * x);
63  double gamma2 = boost::math::tgamma(5.0 / 2.0 * x);
64  double c1 = x * std::pow(gamma2, 1.5) / std::pow(gamma1, 2.5);
65  double c2 = x * std::pow(gamma2 / gamma1, x);
67  c1 * m_electronEnergyLevels.sqrt() /
68  std::pow(meanElectronEnergy(), 1.5) *
69  (-c2 * (m_electronEnergyLevels /
70  meanElectronEnergy()).pow(x)).exp();
72 }
73 
74 void PlasmaPhase::setElectronTemperature(const double Te) {
75  m_electronTemp = Te;
77 }
78 
80  m_electronTemp = 2.0 / 3.0 * energy * ElectronCharge / Boltzmann;
82 }
83 
84 void PlasmaPhase::setElectronEnergyLevels(const double* levels, size_t length)
85 {
86  m_nPoints = length;
87  m_electronEnergyLevels = Eigen::Map<const Eigen::ArrayXd>(levels, length);
90 }
91 
93 {
94  Eigen::ArrayXd h = m_electronEnergyLevels.tail(m_nPoints - 1) -
96  if (m_electronEnergyLevels[0] < 0.0 || (h <= 0.0).any()) {
97  throw CanteraError("PlasmaPhase::checkElectronEnergyLevels",
98  "Values of electron energy levels need to be positive and "
99  "monotonically increasing.");
100  }
101 }
102 
104 {
105  Eigen::ArrayXd h = m_electronEnergyLevels.tail(m_nPoints - 1) -
107  if ((m_electronEnergyDist < 0.0).any()) {
108  throw CanteraError("PlasmaPhase::checkElectronEnergyDistribution",
109  "Values of electron energy distribution cannot be negative.");
110  }
111  if (m_electronEnergyDist[m_nPoints - 1] > 0.01) {
112  warn_user("PlasmaPhase::checkElectronEnergyDistribution",
113  "The value of the last element of electron energy distribution exceed 0.01. "
114  "This indicates that the value of electron energy level is not high enough "
115  "to contain the isotropic distribution at mean electron energy of "
116  "{} eV", meanElectronEnergy());
117  }
118 }
119 
121  const double* dist,
122  size_t length)
123 {
124  m_distributionType = "discretized";
125  m_nPoints = length;
127  Eigen::Map<const Eigen::ArrayXd>(levels, length);
129  Eigen::Map<const Eigen::ArrayXd>(dist, length);
133  }
136 }
137 
139 {
140  // calculate mean electron energy and electron temperature
141  Eigen::ArrayXd eps52 = m_electronEnergyLevels.pow(5./2.);
142  double epsilon_m = 2.0 / 5.0 * numericalQuadrature(m_quadratureMethod,
143  m_electronEnergyDist, eps52);
144  m_electronTemp = 2.0 / 3.0 * epsilon_m * ElectronCharge / Boltzmann;
145 }
146 
148  m_isotropicShapeFactor = x;
150 }
151 
152 void PlasmaPhase::getParameters(AnyMap& phaseNode) const
153 {
154  IdealGasPhase::getParameters(phaseNode);
155  AnyMap eedf;
156  eedf["type"] = m_distributionType;
157  vector<double> levels(m_nPoints);
158  Eigen::Map<Eigen::ArrayXd>(levels.data(), m_nPoints) = m_electronEnergyLevels;
159  eedf["energy-levels"] = levels;
160  if (m_distributionType == "isotropic") {
161  eedf["shape-factor"] = m_isotropicShapeFactor;
162  eedf["mean-electron-energy"].setQuantity(meanElectronEnergy(), "eV");
163  } else if (m_distributionType == "discretized") {
164  vector<double> dist(m_nPoints);
165  Eigen::Map<Eigen::ArrayXd>(dist.data(), m_nPoints) = m_electronEnergyDist;
166  eedf["distribution"] = dist;
167  eedf["normalize"] = m_do_normalizeElectronEnergyDist;
168  }
169  phaseNode["electron-energy-distribution"] = std::move(eedf);
170 }
171 
172 void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode)
173 {
174  IdealGasPhase::setParameters(phaseNode, rootNode);
175  if (phaseNode.hasKey("electron-energy-distribution")) {
176  const AnyMap eedf = phaseNode["electron-energy-distribution"].as<AnyMap>();
177  m_distributionType = eedf["type"].asString();
178  if (m_distributionType == "isotropic") {
179  if (eedf.hasKey("shape-factor")) {
180  setIsotropicShapeFactor(eedf["shape-factor"].asDouble());
181  } else {
182  throw CanteraError("PlasmaPhase::setParameters",
183  "isotropic type requires shape-factor key.");
184  }
185  if (eedf.hasKey("mean-electron-energy")) {
186  double energy = eedf.convert("mean-electron-energy", "eV");
187  setMeanElectronEnergy(energy);
188  } else {
189  throw CanteraError("PlasmaPhase::setParameters",
190  "isotropic type requires electron-temperature key.");
191  }
192  if (eedf.hasKey("energy-levels")) {
193  setElectronEnergyLevels(eedf["energy-levels"].asVector<double>().data(),
194  eedf["energy-levels"].asVector<double>().size());
195  }
197  } else if (m_distributionType == "discretized") {
198  if (!eedf.hasKey("energy-levels")) {
199  throw CanteraError("PlasmaPhase::setParameters",
200  "Cannot find key energy-levels.");
201  }
202  if (!eedf.hasKey("distribution")) {
203  throw CanteraError("PlasmaPhase::setParameters",
204  "Cannot find key distribution.");
205  }
206  if (eedf.hasKey("normalize")) {
207  enableNormalizeElectronEnergyDist(eedf["normalize"].asBool());
208  }
209  setDiscretizedElectronEnergyDist(eedf["energy-levels"].asVector<double>().data(),
210  eedf["distribution"].asVector<double>().data(),
211  eedf["energy-levels"].asVector<double>().size());
212  }
213  }
214 }
215 
216 bool PlasmaPhase::addSpecies(shared_ptr<Species> spec)
217 {
218  bool added = IdealGasPhase::addSpecies(spec);
219  size_t k = m_kk - 1;
220 
221  if (spec->composition.find("E") != spec->composition.end() &&
222  spec->composition.size() == 1 &&
223  spec->composition["E"] == 1) {
224  if (m_electronSpeciesIndex == npos) {
226  } else {
227  throw CanteraError("PlasmaPhase::addSpecies",
228  "Cannot add species, {}. "
229  "Only one electron species is allowed.", spec->name);
230  }
231  }
232  return added;
233 }
234 
236 {
238  // check electron species
239  if (m_electronSpeciesIndex == npos) {
240  throw CanteraError("PlasmaPhase::initThermo",
241  "No electron species found.");
242  }
243 }
244 
246 {
248  static const int cacheId = m_cache.getId();
249  CachedScalar cached = m_cache.getScalar(cacheId);
250  double tempNow = temperature();
251  double electronTempNow = electronTemperature();
252  size_t k = m_electronSpeciesIndex;
253  // If the electron temperature has changed since the last time these
254  // properties were computed, recompute them.
255  if (cached.state1 != tempNow || cached.state2 != electronTempNow) {
257  &m_cp0_R[k], &m_h0_RT[k], &m_s0_R[k]);
258  cached.state1 = tempNow;
259  cached.state2 = electronTempNow;
260  }
261  // update the species Gibbs functions
262  m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
263 }
264 
266  double value = IdealGasPhase::enthalpy_mole();
267  value += GasConstant * (electronTemperature() - temperature()) *
270  return value;
271 }
272 
273 void PlasmaPhase::getGibbs_ref(double* g) const
274 {
277 }
278 
279 void PlasmaPhase::getStandardVolumes_ref(double* vol) const
280 {
283 }
284 
286 {
289 }
290 
292 {
294  double logp = log(pressure());
295  double logpe = log(electronPressure());
296  sbar[m_electronSpeciesIndex] += GasConstant * (logp - logpe);
297 }
298 
300 {
301  const vector<double>& _h = enthalpy_RT_ref();
302  for (size_t k = 0; k < m_kk; k++) {
303  ubar[k] = RT() * (_h[k] - 1.0);
304  }
305  size_t k = m_electronSpeciesIndex;
306  ubar[k] = RTe() * (_h[k] - 1.0);
307 }
308 
309 void PlasmaPhase::getChemPotentials(double* mu) const
310 {
312  size_t k = m_electronSpeciesIndex;
313  double xx = std::max(SmallNumber, moleFraction(k));
314  mu[k] += (RTe() - RT()) * log(xx);
315 }
316 
317 void PlasmaPhase::getStandardChemPotentials(double* muStar) const
318 {
320  size_t k = m_electronSpeciesIndex;
321  muStar[k] -= log(pressure() / refPressure()) * RT();
322  muStar[k] += log(electronPressure() / refPressure()) * RTe();
323 }
324 
325 void PlasmaPhase::getEntropy_R(double* sr) const
326 {
327  const vector<double>& _s = entropy_R_ref();
328  copy(_s.begin(), _s.end(), sr);
329  double tmp = log(pressure() / refPressure());
330  for (size_t k = 0; k < m_kk; k++) {
331  if (k != m_electronSpeciesIndex) {
332  sr[k] -= tmp;
333  } else {
334  sr[k] -= log(electronPressure() / refPressure());
335  }
336  }
337 }
338 
339 void PlasmaPhase::getGibbs_RT(double* grt) const
340 {
341  const vector<double>& gibbsrt = gibbs_RT_ref();
342  copy(gibbsrt.begin(), gibbsrt.end(), grt);
343  double tmp = log(pressure() / refPressure());
344  for (size_t k = 0; k < m_kk; k++) {
345  if (k != m_electronSpeciesIndex) {
346  grt[k] += tmp;
347  } else {
348  grt[k] += log(electronPressure() / refPressure());
349  }
350  }
351 }
352 
353 }
Header file for class PlasmaPhase.
Declaration for class Cantera::Species.
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:427
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition: AnyMap.cpp:1423
double convert(const string &key, const string &units) const
Convert the item stored by the given key to the units specified in units.
Definition: AnyMap.cpp:1535
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
double enthalpy_mole() const override
Return the Molar enthalpy. Units: J/kmol.
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
double pressure() const override
Pressure.
vector< double > m_g0_RT
Temporary storage for dimensionless reference state Gibbs energies.
const vector< double > & gibbs_RT_ref() const
Returns a reference to the dimensionless reference state Gibbs free energy vector.
vector< double > m_h0_RT
Temporary storage for dimensionless reference state enthalpies.
void getGibbs_ref(double *g) const override
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
void getStandardChemPotentials(double *mu) const override
Get the array of chemical potentials at unit activity for the species at their standard states at the...
void getStandardVolumes_ref(double *vol) const override
Get the molar volumes of the species reference states at the current T and P_ref of the solution.
virtual void updateThermo() const
Update the species reference state thermodynamic functions.
vector< double > m_s0_R
Temporary storage for dimensionless reference state entropies.
const vector< double > & enthalpy_RT_ref() const
Returns a reference to the dimensionless reference state enthalpy vector.
vector< double > m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
const vector< double > & entropy_R_ref() const
Returns a reference to the dimensionless reference state Entropy vector.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
virtual void update_single(size_t k, double T, double *cp_R, double *h_RT, double *s_R) const
Get reference-state properties for a single species.
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
Definition: Phase.h:822
size_t m_kk
Number of species in the phase.
Definition: Phase.h:842
double temperature() const
Temperature (K).
Definition: Phase.h:562
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:439
void checkElectronEnergyDistribution() const
Check the electron energy distribution.
double meanElectronEnergy() const
Mean electron energy [eV].
Definition: PlasmaPhase.h:148
double enthalpy_mole() const override
Return the Molar enthalpy. Units: J/kmol.
size_t m_nPoints
Number of points of electron energy levels.
Definition: PlasmaPhase.h:292
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
void normalizeElectronEnergyDistribution()
Electron energy distribution norm.
Definition: PlasmaPhase.cpp:35
void getStandardChemPotentials(double *muStar) const override
Get the array of chemical potentials at unit activity for the species at their standard states at the...
void updateThermo() const override
Update the species reference state thermodynamic functions.
void setElectronTemperature(double Te) override
Set the internally stored electron temperature of the phase (K).
Definition: PlasmaPhase.cpp:74
void getEntropy_R(double *sr) const override
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
void getGibbs_ref(double *g) const override
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
Eigen::ArrayXd m_electronEnergyDist
Normalized electron energy distribution vector [-] Length: m_nPoints.
Definition: PlasmaPhase.h:299
Eigen::ArrayXd m_electronEnergyLevels
electron energy levels [ev]. Length: m_nPoints
Definition: PlasmaPhase.h:295
void setDiscretizedElectronEnergyDist(const double *levels, const double *distrb, size_t length)
Set discretized electron energy distribution.
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
string type() const override
String indicating the thermodynamic model implemented.
Definition: PlasmaPhase.h:74
void checkElectronEnergyLevels() const
Check the electron energy levels.
Definition: PlasmaPhase.cpp:92
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
void updateElectronTemperatureFromEnergyDist()
Update electron temperature (K) From energy distribution.
string m_distributionType
Electron energy distribution type.
Definition: PlasmaPhase.h:308
void getStandardVolumes_ref(double *vol) const override
Get the molar volumes of the species reference states at the current T and P_ref of the solution.
void updateElectronEnergyDistribution()
Update electron energy distribution.
Definition: PlasmaPhase.cpp:25
string m_quadratureMethod
Numerical quadrature method for electron energy distribution.
Definition: PlasmaPhase.h:311
PlasmaPhase(const string &inputFile="", const string &id="")
Construct and initialize a PlasmaPhase object directly from an input file.
Definition: PlasmaPhase.cpp:14
double m_electronTemp
Electron temperature [K].
Definition: PlasmaPhase.h:305
double RTe() const
Return the Gas Constant multiplied by the current electron temperature.
Definition: PlasmaPhase.h:176
void setElectronEnergyLevels(const double *levels, size_t length)
Set electron energy levels.
Definition: PlasmaPhase.cpp:84
void getGibbs_RT(double *grt) const override
Get the nondimensional Gibbs functions for the species in their standard states at the current T and ...
bool m_do_normalizeElectronEnergyDist
Flag of normalizing electron energy distribution.
Definition: PlasmaPhase.h:314
void getPartialMolarIntEnergies(double *ubar) const override
Return an array of partial molar internal energies for the species in the mixture.
void setIsotropicElectronEnergyDistribution()
Set isotropic electron energy distribution.
Definition: PlasmaPhase.cpp:58
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
void setParameters(const AnyMap &phaseNode, const AnyMap &rootNode=AnyMap()) override
Set equation of state parameters from an AnyMap phase description.
void setMeanElectronEnergy(double energy)
Set mean electron energy [eV].
Definition: PlasmaPhase.cpp:79
size_t m_electronSpeciesIndex
Index of electron species.
Definition: PlasmaPhase.h:302
void setElectronEnergyDistributionType(const string &type)
Set electron energy distribution type.
Definition: PlasmaPhase.cpp:47
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
virtual double electronPressure() const
Electron pressure.
Definition: PlasmaPhase.h:184
double electronTemperature() const override
Electron Temperature (K)
Definition: PlasmaPhase.h:168
void setIsotropicShapeFactor(double x)
Set the shape factor of isotropic electron energy distribution.
void enableNormalizeElectronEnergyDist(bool enable)
Set flag of automatically normalize electron energy distribution Flag: m_do_normalizeElectronEnergyDi...
Definition: PlasmaPhase.h:154
virtual void setParameters(const AnyMap &phaseNode, const AnyMap &rootNode=AnyMap())
Set equation of state parameters from an AnyMap phase description.
virtual void getParameters(AnyMap &phaseNode) const
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
double RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:1062
virtual void initThermo()
Initialize the ThermoPhase object after all species have been set up.
void initThermoFile(const string &inputFile, const string &id)
Initialize a ThermoPhase object using an input file.
MultiSpeciesThermo m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
Definition: ThermoPhase.h:1962
virtual double refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:436
CachedScalar getScalar(int id)
Get a reference to a CachedValue object representing a scalar (double) with the given id.
Definition: ValueCache.h:161
int getId()
Get a unique id for a cached value.
Definition: ValueCache.cpp:21
Header for a file containing miscellaneous numerical functions.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
double numericalQuadrature(const string &method, const Eigen::ArrayXd &f, const Eigen::ArrayXd &x)
Numerical integration of a function.
Definition: funcs.cpp:112
const double Boltzmann
Boltzmann constant [J/K].
Definition: ct_defs.h:84
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:120
const double ElectronCharge
Elementary charge [C].
Definition: ct_defs.h:90
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition: global.h:267
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:158
A cached property value and the state at which it was evaluated.
Definition: ValueCache.h:33
double state2
Value of the second state variable for the state at which value was evaluated, for example density or...
Definition: ValueCache.h:106
double state1
Value of the first state variable for the state at which value was evaluated, for example temperature...
Definition: ValueCache.h:102