mechanism.h

Go to the documentation of this file.

 
#ifndef CHEMISTRY_MECHANISM_H
#define CHEMISTRY_MECHANISM_H
 
#include <array>
#include <cstdint>
#include <string>
#include <utility>
#include <vector>
 
namespace chemistry {
 
// Physical constants (CODATA 2018 exact values, as in Cantera); quantities in kmol
constexpr double AVOGADRO = 6.02214076e26;           // 1/kmol
constexpr double BOLTZMANN = 1.380649e-23;           // J/K
constexpr double GAS_CONSTANT = AVOGADRO * BOLTZMANN; // J/(kmol K)
constexpr double ELEMENTARY_CHARGE = 1.602176634e-19; // C
constexpr double ONE_ATM = 101325.0;                 // Pa
constexpr double LIGHT_SPEED = 299792458.0;           // m/s
constexpr double EPSILON_0 = 8.854187812773345e-12;  // F/m, as in Cantera
 
double atomic_weight(const std::string & symbol);
 
enum class ThermoModel {
    NASA7,
    NASA9,
    CONSTANT_CP,
};
 
struct SpeciesThermo {
    ThermoModel model = ThermoModel::NASA7;
    std::vector<double> T_bounds;               // n_ranges + 1, increasing
    std::vector<std::array<double, 9>> coeffs;  // per range
    std::array<double, 4> constant_cp = {};     // CONSTANT_CP: T0 [K], h0 [J/kmol], s0, cp0 [J/(kmol K)]
 
    size_t range(double T) const;
};
 
enum class MoleculeGeometry {
    ATOM,
    LINEAR,
    NONLINEAR,
};
 
struct SpeciesTransport {
    MoleculeGeometry geometry = MoleculeGeometry::ATOM;
    double diameter = 0.0;               // Lennard-Jones collision diameter [m]
    double well_depth = 0.0;             // Lennard-Jones well depth [J]
    double dipole = 0.0;                 // [C m]
    double polarizability = 0.0;         // [m^3]
    double rotational_relaxation = 0.0;  // collision number at 298 K
};
 
struct Species {
    std::string name;
    std::vector<std::pair<std::string, double>> composition;  // element, atoms
    double molecular_weight = 0.0;                            // kg/kmol
    SpeciesThermo thermo;
    bool has_transport = false;
    SpeciesTransport transport;
};
 
struct Arrhenius {
    double A = 0.0;
    double b = 0.0;
    double Ea_R = 0.0;  // activation energy over the gas constant [K]
};
 
enum class ReactionType {
    ELEMENTARY,
    THREE_BODY,
    FALLOFF,
    PLOG,       // pressure-dependent Arrhenius: ln k interpolated linearly in ln p
    CHEBYSHEV,  // log10 k as a Chebyshev series in 1/T and log10 p
};
 
enum class FalloffType {
    LINDEMANN,
    TROE,
    SRI,
};
 
struct Reaction {
    std::string equation;
    ReactionType type = ReactionType::ELEMENTARY;
    std::vector<std::pair<uint32_t, double>> reactants;  // species, stoichiometric coefficient
    std::vector<std::pair<uint32_t, double>> products;
    std::vector<std::pair<uint32_t, double>> orders;     // forward reaction orders of the reactants
    bool reversible = true;
    bool duplicate = false;
    Arrhenius rate;  // falloff: the high-pressure limit
    Arrhenius low;   // falloff: the low-pressure limit
    FalloffType falloff = FalloffType::LINDEMANN;
    std::array<double, 5> falloff_params = {};  // Troe A, T3, T1, T2; SRI a, b, c, d, e
    std::vector<std::pair<uint32_t, double>> efficiencies;  // third-body efficiencies other than the default
    double default_efficiency = 1.0;
    std::vector<std::pair<double, Arrhenius>> plog;  // PLOG: (pressure [Pa], rate), by increasing pressure
    std::array<double, 4> chebyshev_range = {};     // Chebyshev: T_min, T_max [K], p_min, p_max [Pa]
    uint32_t chebyshev_n_T = 0, chebyshev_n_p = 0;
    std::vector<double> chebyshev;                  // (n_T x n_p) coefficients of log10 k, k in SI with kmol
};
 
struct Mechanism {
    std::string file;
    std::string phase;
    std::vector<std::string> elements;
    std::vector<Species> species;
    std::vector<Reaction> reactions;
 
    uint32_t n_species() const { return static_cast<uint32_t>(species.size()); }
 
    int32_t species_index(const std::string & name) const;
 
    std::vector<std::string> species_names() const;
};
 
Mechanism read_mechanism(const std::string & file, const std::string & phase = "");
 
} // namespace chemistry
 
#endif // CHEMISTRY_MECHANISM_H