Add condensation to Kinetics (#16)

- Add condensation class to Kinetics
- Add vapor pressures and Nucleation

include/cantera/kinetics/Condensation.h

@@ -0,0 +1,111 @@
+#ifndef CT_CONDENSATION_H
+#define CT_CONDENSATION_H
+
+#include "Kinetics.h"
+
+namespace Cantera
+{
+
+class Condensation : public Kinetics {
+ public:
+  Condensation() = default;
+
+  ~Condensation() override {}
+
+  void setQuantityMoleFraction() {
+    if (!m_use_mole_fraction) {
+      m_ROP_ok = false;
+    }
+    m_use_mole_fraction = true;
+  }
+
+  void setQuantityConcentration() {
+    if (m_use_mole_fraction) {
+      m_ROP_ok = false;
+    }
+    m_use_mole_fraction = false;
+  }
+
+  void resizeReactions() override;
+
+  string kineticsType() const override {
+    return "nucleation";
+  }
+
+  void getActivityConcentrations(double* const pdata) override;
+
+  bool isReversible(size_t i) override {
+    return true;
+  }
+
+  void getFwdRateConstants(double* kfwd) override;
+
+  void getFwdRateConstants_ddT(double* dkfwd) override;
+
+  void resizeSpecies() override;
+
+  bool addReaction(shared_ptr<Reaction> r, bool resize=true) override;
+
+  void updateROP() override;
+
+  Eigen::SparseMatrix<double> netRatesOfProgress_ddCi() override;
+  Eigen::SparseMatrix<double> netRatesOfProgress_ddX() override;
+
+ protected:
+  void _update_rates_T(double *pdata, double *pdata_ddT);
+  void _update_rates_C(double *pdata);
+  void _update_rates_X(double *pdata);
+
+  //! Process derivatives
+  //! @param stoich  stoichiometry manager
+  //! @param in  rate expression used for the derivative calculation
+  //! @param ddX true: w.r.t mole fractions false: w.r.t species concentrations
+  //! @return a sparse matrix of derivative contributions for each reaction of
+  //! dimensions nTotalReactions by nTotalSpecies
+  Eigen::SparseMatrix<double> calculateCompositionDerivatives(
+      StoichManagerN& stoich, const vector<double>& in, bool ddX=true);
+
+  //! This variable has two interpretations.
+  //! If m_use_mole_fraction is true, then it is the vector of mole fractions.
+  //! If m_use_mole_fraction is false, then it is the vector of concentrations.
+  Eigen::VectorXd m_conc;
+  Eigen::VectorXd m_intEng;
+  Eigen::VectorXd m_cv;
+
+  //! Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for
+  //! EdgeKinetics)
+  size_t m_nDim = 2;
+
+  //! Vector of rate handlers for interface reactions
+  vector<unique_ptr<MultiRateBase>> m_interfaceRates;
+  map<string, size_t> m_interfaceTypes; //!< Rate handler mapping
+
+  //! reaction index for x <=> y
+  vector<size_t> m_jxy;
+
+  //! reaction index for x1 + x2 <=> y
+  vector<size_t> m_jxxy;
+
+  //! reaction index for freezing reaction x(l) <=> x(s)
+  vector<size_t> m_jyy;
+
+  //! rate jacobian matrix
+  Eigen::SparseMatrix<double> m_jac;
+
+  //! rate jacobian with respect to temperature
+  vector<double> m_rfn_ddT;
+
+  bool m_use_mole_fraction = false;
+  bool m_ROP_ok = false;
+
+  //! Current temperature of the data
+  double m_temp = 0.0;
+
+  //! Buffers for partial rop results with length nReactions()
+  vector<double> m_rbuf0;
+  vector<double> m_rbuf1;
+};
+
+}
+
+#endif

include/cantera/kinetics/Freezing.h

@@ -0,0 +1,55 @@
+#ifndef CT_FREEZING_H
+#define CT_FREEZING_H
+
+#include <functional>
+
+#include "cantera/base/ct_defs.h"
+#include "cantera/base/Units.h"
+#include "cantera/kinetics/Arrhenius.h"
+#include "ReactionRate.h"
+#include "MultiRate.h"
+
+namespace Cantera
+{
+
+class AnyValue;
+class AnyMap;
+
+class FreezingRate : public ReactionRate {
+ public:
+  FreezingRate() = default;
+  FreezingRate(const AnyMap& node, const UnitStack& rate_units);
+
+  unique_ptr<MultiRateBase> newMultiRate() const override {
+    return make_unique<MultiRate<FreezingRate, ArrheniusData>>();
+  }
+
+  //! Set the rate parameters for this reaction.
+  void setRateParameters(const AnyValue& equation,
+                         const AnyValue& rate,
+                         const AnyMap& node);
+
+  //! return the rate coefficient type
+  const string type() const override { return "freezing"; }
+
+  void getParameters(AnyMap& rateNode, const Units& rate_units=Units(0.)) const;
+  using ReactionRate::getParameters;
+
+  void validate(const string& equation, const Kinetics& kin) override;
+
+  double evalFromStruct(const ArrheniusData& shared_data) const;
+
+  double ddTScaledFromStruct(const ArrheniusData& shared_data) const {
+    return 0.;
+  }
+
+ protected:
+  string m_formula_str = "formula";
+  std::function<double(double)> m_meltfunc;
+
+  double m_t3 = 0.0;
+};
+
+}
+
+#endif

include/cantera/kinetics/Kinetics.h

@@ -303,6 +303,11 @@ class Kinetics
      */
     size_t kineticsSpeciesIndex(const string& nm) const;
 
+    //! shorthand for kineticsSpeciesIndex
+    size_t speciesIndex(const string& nm) const {
+      return kineticsSpeciesIndex(nm);
+    }
+
     /**
      * This function looks up the name of a species and returns a
      * reference to the ThermoPhase object of the phase where the species

include/cantera/kinetics/Nucleation.h

@@ -0,0 +1,69 @@
+#ifndef CT_NUCLEATION_H
+#define CT_NUCLEATION_H
+
+#include <functional>
+
+#include "cantera/base/ct_defs.h"
+#include "cantera/base/Units.h"
+#include "cantera/kinetics/Arrhenius.h"
+#include "ReactionRate.h"
+#include "MultiRate.h"
+
+namespace Cantera
+{
+
+class AnyValue;
+class AnyMap;
+
+class NucleationRate : public ReactionRate {
+ public:
+  NucleationRate() = default;
+  NucleationRate(const AnyMap& node, const UnitStack& rate_units);
+
+  unique_ptr<MultiRateBase> newMultiRate() const override {
+    return make_unique<MultiRate<NucleationRate, ArrheniusData>>();
+  }
+
+  //! Set the rate parameters for this reaction.
+  void setRateParameters(const AnyValue& equation,
+                         const AnyValue& rate,
+                         const AnyMap& node);
+
+  //! return the rate coefficient type
+  const string type() const override { return "nucleation"; }
+
+  void getParameters(AnyMap& rateNode, const Units& rate_units=Units(0.)) const;
+  using ReactionRate::getParameters;
+
+  void validate(const string& equation, const Kinetics& kin) override;
+
+  double evalFromStruct(const ArrheniusData& shared_data) const;
+
+  double ddTScaledFromStruct(const ArrheniusData& shared_data) const;
+
+ protected:
+  string m_formula_str = "formula";
+
+  //! returns s = svp(T)/RT or svp(T)/(RT)^2
+  std::function<double(double)> m_svp;
+
+  //! returns d(log(s))/dT
+  std::function<double(double)> m_logsvp_ddT;
+
+  size_t m_order = 1;
+
+  double m_t3 = 0.0;
+  double m_p3 = 0.0;
+  double m_beta = 0.0;
+  double m_delta = 0.0;
+
+  double m_min_temp = 0.0;
+  double m_max_temp = 1.0e30;
+};
+
+std::function<double(double)> find_svp_function(const Composition& reactants, 
+                                                const string& svp_name);
+
+}
+
+#endif

include/cantera/kinetics/svp_funcs.h

@@ -0,0 +1,8 @@
+#ifndef SVP_FUNCS_INC_H
+#define SVP_FUNCS_INC_H
+
+#include "vapors/water_vapors.hpp"
+#include "vapors/ammonia_vapors.hpp"
+#include "vapors/ammonium_hydrosulfide_vapors.hpp"
+
+#endif // SVP_FUNCS_INC_H

include/cantera/kinetics/vapors/ammonia.dat

@@ -0,0 +1,59 @@
+# Ammonia Thermodynamic Properties from Engineering Toolbox
+# https://www.engineeringtoolbox.com/ammonia-d_971.html
+#
+# T - temperature (K)
+# P - pressure (bar)
+# alpha - specific volume (m^3/kg)
+# Hf(l) - liquid specific enthalpy (kJ/kg)
+# Hf(v) - vapor specific enthalpy (kJ/kg)
+# Sf(l) - liquid specific entropy (kJ/kg)
+# Sf(v) - vapor specific entropy (kJ/kg)
+#
+#                     Saturated Properties            Superheated Properties
+#                     ----------------------------    -------------------------
+# T   P       alpha   Hf(l)    Hf(v)  Sf(l)   Sf(v)   H(50K)  S(50K)  H(100K) S(100K)
+-50 	0.4089 	2.625 	-44.4 	1373.3 	-0.194 	6.159 	1479.8 	6.592 	1585.9 	6.948
+-45 	0.5454 	2.005 	-22.3 	1381.6 	-0.096 	6.057 	1489.3 	6.486 	1596.1 	6.839
+-40 	0.7177 	1.552 	0 	1390.0 	0 	5.962 	1498.6 	6.387 	1606.3 	6.736
+-35 	0.9322 	1.216 	22.3 	1397.9 	0.095 	5.872 	1507.9 	6.293 	1616.3 	6.639
+-30 	1.196 	0.9633 	44.7 	1405.6 	0.188 	5.785 	1517.0 	6.203 	1626.3 	6.547
+-28 	1.317 	0.8809 	53.6 	1408.5 	0.224 	5.751 	1520.7 	6.169 	1630.3 	6.512
+-26 	1.447 	0.8058 	62.6 	1411.4 	0.261 	5.718 	1524.3 	6.135 	1634.2 	6.477
+-24 	1.588 	0.7389 	71.7 	1414.3 	0.297 	5.686 	1527.9 	6.103 	1638.2 	6.444
+-22 	1.740 	0.6783 	80.8 	1417.3 	0.333 	5.655 	1531.4 	6.071 	1642.2 	6.411
+-20 	1.902 	0.6237 	89.8 	1420.0 	0.368 	5.623 	1534.8 	6.039 	1646.0 	6.379
+-18 	2.077 	0.5743 	98.8 	1422.7 	0.404 	5.593 	1538.2 	6.008 	1650.0 	6.347
+-16 	2.265 	0.5296 	107.9 	1425.3 	0.440 	5.563 	1541.7 	5.978 	1653.8 	6.316
+-14 	2.465 	0.4890 	117.0 	1427.9 	0.475 	5.533 	1545.1 	5.948 	1657.7 	6.286
+-12 	2.680 	0.4521 	126.2 	1430.5 	0.510 	5.504 	1548.5 	5.919 	1661.5 	6.256
+-10 	2.908 	0.4185 	135.4 	1433.0 	0.544 	5.475 	1551.7 	5.891 	1665.3 	6.227
+-8 	3.153 	0.3879 	144.5 	1435.3 	0.579 	5.447 	1554.9 	5.863 	1669.0 	6.199
+-6 	3.413 	0.3599 	153.6 	1437.6 	0.613 	5.419 	1558.2 	5.836 	1672.8 	6.171
+-4 	3.691 	0.3344 	162.8 	1439.9 	0.647 	5.392 	1561.4 	5.808 	1676.4 	6.143
+-2 	3.983 	0.3110 	172.0 	1442.2 	0.681 	5.365 	1564.6 	5.782 	1680.1 	6.116
+0 	4.295 	0.2895 	181.2 	1444.4 	0.715 	5.340 	1567.8 	5.756 	1683.9 	6.090
+2 	4.625 	0.2699 	190.4 	1446.5 	0.749 	5.314 	1570.9 	5.731 	1687.5 	6.065
+4 	4.975 	0.2517 	199.7 	1448.5 	0.782 	5.288 	1574.0 	5.706 	1691.2 	6.040
+6 	5.346 	0.2351 	209.1 	1450.6 	0.816 	5.263 	1577.0 	5.682 	1694.9 	6.015
+8 	5.736 	0.2198 	218.5 	1452.5 	0.849 	5.238 	1580.1 	5.658 	1698.4 	5.991
+10 	6.149 	0.2056 	227.8 	1454.3 	0.881 	5.213 	1583.1 	5.634 	1702.2 	5.967
+12 	6.585 	0.1926 	237.2 	1456.1 	0.914 	5.189 	1586.0 	5.611 	1705.7 	5.943
+14 	7.045 	0.1805 	246.6 	1457.8 	0.947 	5.165 	1588.9 	5.588 	1709.1 	5.920
+16 	7.529 	0.1693 	256.0 	1459.5 	0.979 	5.141 	1591.7 	5.565 	1712.5 	5.898
+18 	8.035 	0.1590 	265.5 	1461.1 	1.012 	5.118 	1594.4 	5.543 	1715.9 	5.876
+20 	8.570 	0.1494 	275.1 	1462.6 	1.044 	5.095 	1597.2 	5.521 	1719.3 	5.854
+22 	9.134 	0.1405 	284.6 	1463.9 	1.076 	5.072 	1600.0 	5.499 	1722.8 	5.832
+24 	9.722 	0.1322 	294.1 	1465.2 	1.108 	5.049 	1602.7 	5.478 	1726.3 	5.811
+26 	10.34 	0.1245 	303.7 	1466.5 	1.140 	5.027 	1605.3 	5.458 	1729.6 	5.790
+28 	10.99 	0.1173 	313.4 	1467.8 	1.172 	5.005 	1608.0 	5.437 	1732.7 	5.770
+30 	11.67 	0.1106 	323.1 	1468.9 	1.204 	4.984 	1610.5 	5.417 	1735.9 	5.750
+32 	12.37 	0.1044 	332.8 	1469.9 	1.235 	4.962 	1613.0 	5.397 	1739.3 	5.731
+34 	13.11 	0.0986 	342.5 	1470.8 	1.267 	4.940 	1615.4 	5.378 	1742.6 	5.711
+36 	13.89 	0.0931 	352.3 	1471.8 	1.298 	4.919 	1617.8 	5.358 	1745.7 	5.692
+38 	14.70 	0.0880 	362.1 	1472.6 	1.329 	4.898 	1620.1 	5.340 	1748.7 	5.674
+40 	15.54 	0.0833 	371.9 	1473.3 	1.360 	4.877 	1622.4 	5.321 	1751.9 	5.655
+42 	16.42 	0.0788 	381.8 	1473.8 	1.391 	4.856 	1624.6 	5.302 	1755.0 	5.637
+44 	17.34 	0.0746 	391.8 	1474.2 	1.422 	4.835 	1626.8 	5.284 	1758.0 	5.619
+46 	18.30 	0.0706 	401.8 	1474.5 	1.453 	4.814 	1629.0 	5.266 	1761.0 	5.602
+48 	19.29 	0.0670 	411.9 	1474.7 	1.484 	4.793 	1631.1 	5.248 	1764.0 	5.584
+50 	20.33 	0.0635 	421.9 	1474.7 	1.515 	4.773 	1633.1 	5.230 	1766.8 	5.567

include/cantera/kinetics/vapors/ammonia_vapors.hpp

@@ -0,0 +1,78 @@
+#ifndef SRC_SNAP_THERMODYNAMICS_VAPORS_AMMONIA_VAPORS_HPP_
+#define SRC_SNAP_THERMODYNAMICS_VAPORS_AMMONIA_VAPORS_HPP_
+
+// C/C++
+#include <cmath>
+
+inline double svpnh3(double t, double p, double beta, double gamma) {
+  return p * exp((1. - 1. / t) * beta - gamma * log(t));
+}
+
+inline double svp_nh3_UMich(double T) {
+  double x = -1790.00 / T - 1.81630 * log10(T) + 14.97593;
+  return pow(10.0, x) * 1.E5 / 760.;
+}
+
+// (164, 371.5)
+inline double svp_nh3_Antoine(double T) {
+  double result;
+  if (T < 239.6)
+    result = pow(10., 3.18757 - (506.713 / (T - 80.78)));
+  else
+    result = pow(10., 4.86886 - (1113.928 / (T - 10.409)));
+  return 1.E5 * result;
+}
+
+// (130, 200)
+inline double sat_vapor_p_NH3_Hubner(double T) {
+  double A = 24.3037, B = -1766.28, C = -5.64472, D = 0.00740241;
+
+  double x = A + B / T + C * log10(T) + D * T;
+  return pow(10.0, x);
+}
+
+inline double sat_vapor_p_NH3_BriggsS(double T) {
+  double a[6], x;
+  if (T < 195) {
+    a[1] = -4122.;
+    a[2] = 41.67871;
+    a[3] = -1.81630;
+    a[4] = 0.;
+    a[5] = 0.;
+  } else {
+    a[1] = -4409.3512;
+    a[2] = 76.864252;
+    a[3] = -8.4598340;
+    a[4] = 5.51029e-03;
+    a[5] = 6.80463e-06;
+  }
+  x = a[1] / T + a[2] + a[3] * log(T) + a[4] * T + a[5] * pow(T, 2);
+  return exp(x) / 10.;
+}
+
+// (15, 195.4)
+inline double sat_vapor_p_NH3_Fray(double T) {
+  double a[7], x = 0;
+  a[0] = 1.596e+01;
+  a[1] = -3.537e+03;
+  a[2] = -3.310e+04;
+  a[3] = 1.742e+06;
+  a[4] = -2.995e+07;
+  a[5] = 0.;
+  a[6] = 0.;
+
+  for (int i = 1; i < 7; i++)
+    x = x + a[i] / pow(T, i);  // best fit in [15K; 195.41K]
+
+  return 1.E5 * exp(x + a[0]);
+}
+
+inline double sat_vapor_p_NH3_Ideal(double T) {
+  double betal = 20.08, gammal = 5.62, betas = 20.64, gammas = 1.43, tr = 195.4,
+         pr = 6060.;
+
+  return T > tr ? svpnh3(T / tr, pr, betal, gammal)
+                : svpnh3(T / tr, pr, betas, gammas);
+}
+
+#endif  // SRC_SNAP_THERMODYNAMICS_VAPORS_AMMONIA_VAPORS_HPP_