diff --git a/doc/sphinx/cython/thermo.rst b/doc/sphinx/cython/thermo.rst index 98c51eb70..00b8bb01e 100644 --- a/doc/sphinx/cython/thermo.rst +++ b/doc/sphinx/cython/thermo.rst @@ -9,3 +9,14 @@ These classes are used to describe the thermodynamic state of a system. .. autoclass:: InterfacePhase(infile='', phaseid='') .. autoclass:: PureFluid(infile='', phaseid='') .. autoclass:: Mixture + +Species Thermodynamic Properties +================================ + +These classes are used to describe the reference-state thermodynamic properties +of a pure species. + +.. autoclass:: SpeciesThermo +.. autoclass:: ConstantCp(T_low, T_high, P_ref, coeffs) +.. autoclass:: NasaPoly2(T_low, T_high, P_ref, coeffs) +.. autoclass:: ShomatePoly2(T_low, T_high, P_ref, coeffs) diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 73a5f585a..2f1f7f987 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -43,17 +43,27 @@ cdef extern from "cantera/base/smart_ptr.h": T* get() void reset(T*) +cdef extern from "cantera/thermo/SpeciesThermoInterpType.h": + cdef cppclass CxxSpeciesThermo "Cantera::SpeciesThermoInterpType": + CxxSpeciesThermo() + int reportType() + void updatePropertiesTemp(double, double*, double*, double*) except + + +cdef extern from "cantera/thermo/SpeciesThermoFactory.h": + cdef CxxSpeciesThermo* CxxNewSpeciesThermo "Cantera::newSpeciesThermoInterpType"\ + (int, double, double, double, double*) except + + cdef extern from "cantera/thermo/Species.h": cdef cppclass CxxSpecies "Cantera::Species": CxxSpecies() CxxSpecies(string, stdmap[string,double]) + shared_ptr[CxxSpeciesThermo] thermo string name stdmap[string,double] composition double charge double size - cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera": cdef cppclass CxxThermoPhase "Cantera::ThermoPhase": CxxThermoPhase() @@ -611,6 +621,11 @@ cdef class Species: cdef _assign(self, shared_ptr[CxxSpecies] other) +cdef class SpeciesThermo: + cdef shared_ptr[CxxSpeciesThermo] _spthermo + cdef CxxSpeciesThermo* spthermo + cdef _assign(self, shared_ptr[CxxSpeciesThermo] other) + cdef class _SolutionBase: cdef CxxThermoPhase* thermo cdef CxxKinetics* kinetics @@ -771,3 +786,4 @@ cdef np.ndarray get_reaction_array(Kinetics kin, kineticsMethod1d method) cdef np.ndarray get_transport_1d(Transport tran, transportMethod1d method) cdef np.ndarray get_transport_2d(Transport tran, transportMethod2d method) cdef CxxIdealGasPhase* getIdealGasPhase(ThermoPhase phase) except * +cdef wrapSpeciesThermo(shared_ptr[CxxSpeciesThermo] spthermo) diff --git a/interfaces/cython/cantera/_cantera.pyx b/interfaces/cython/cantera/_cantera.pyx index 109f0c46d..416f082d3 100644 --- a/interfaces/cython/cantera/_cantera.pyx +++ b/interfaces/cython/cantera/_cantera.pyx @@ -14,6 +14,7 @@ include "constants.pyx" include "func1.pyx" include "base.pyx" +include "speciesthermo.pyx" include "thermo.pyx" include "kinetics.pyx" include "transport.pyx" diff --git a/interfaces/cython/cantera/speciesthermo.pyx b/interfaces/cython/cantera/speciesthermo.pyx new file mode 100644 index 000000000..45ea23b2b --- /dev/null +++ b/interfaces/cython/cantera/speciesthermo.pyx @@ -0,0 +1,96 @@ +cdef extern from "cantera/thermo/speciesThermoTypes.h" namespace "Cantera": + cdef int SPECIES_THERMO_CONSTANT_CP "CONSTANT_CP" + cdef int SPECIES_THERMO_NASA2 "NASA2" + cdef int SPECIES_THERMO_SHOMATE2 "SHOMATE2" + + +cdef class SpeciesThermo: + """ + Base class for representing the reference-state thermodynamic properties of + a pure species. These properties are a function of temperature. Derived + classes implement a parameterization of this temperature dependence. This is + a wrapper for the C++ class :ct:`SpeciesThermoInterpType`. + """ + def __cinit__(self, T_low=None, T_high=None, P_ref=None, coeffs=None, *args, + init=True, **kwargs): + if not init: + return + + if len(coeffs) != self.n_coeffs: + raise ValueError("Coefficient array has incorrect length") + cdef np.ndarray[np.double_t, ndim=1] data = np.ascontiguousarray( + coeffs, dtype=np.double) + self._spthermo.reset(CxxNewSpeciesThermo(self.derived_type, T_low, + T_high, P_ref, &data[0])) + self.spthermo = self._spthermo.get() + + cdef _assign(self, shared_ptr[CxxSpeciesThermo] other): + self._spthermo = other + self.spthermo = self._spthermo.get() + + def cp(self, T): + """ Molar heat capacity at constant pressure [J/kmol/K] """ + cdef double cp_r, h_rt, s_r + self.spthermo.updatePropertiesTemp(T, &cp_r, &h_rt, &s_r) + return cp_r * gas_constant + + def h(self, T): + """ Molar enthalpy [J/kmol] """ + cdef double cp_r, h_rt, s_r + self.spthermo.updatePropertiesTemp(T, &cp_r, &h_rt, &s_r) + return h_rt * gas_constant * T + + def s(self, T): + """ Molar entropy [J/kmol/K] """ + cdef double cp_r, h_rt, s_r + self.spthermo.updatePropertiesTemp(T, &cp_r, &h_rt, &s_r) + return s_r * gas_constant + + +cdef class ConstantCp(SpeciesThermo): + """ + Thermodynamic properties for a species that has a constant specific heat + capacity. This is a wrapper for the C++ class :ct:`ConstCpPoly`. + """ + derived_type = SPECIES_THERMO_CONSTANT_CP + n_coeffs = 4 + + +cdef class NasaPoly2(SpeciesThermo): + """ + Thermodynamic properties for a species which is parameterized using the + 7-coefficient NASA polynomial form in two temperature ranges. This is a + wrapper for the C++ class :ct:`NasaPoly2`. + """ + derived_type = SPECIES_THERMO_NASA2 + n_coeffs = 15 + + +cdef class ShomatePoly2(SpeciesThermo): + """ + Thermodynamic properties for a species which is parameterized using the + Shomate equation in two temperature ranges. This is a wrapper for the C++ + class :ct:`ShomatePoly2`. + """ + derived_type = SPECIES_THERMO_SHOMATE2 + n_coeffs = 15 + + +cdef wrapSpeciesThermo(shared_ptr[CxxSpeciesThermo] spthermo): + """ + Wrap a C++ SpeciesThermoInterpType object with a Python object of the + correct derived type. + """ + cdef int thermo_type = spthermo.get().reportType() + + if thermo_type == SPECIES_THERMO_NASA2: + st = NasaPoly2(init=False) + elif thermo_type == SPECIES_THERMO_CONSTANT_CP: + st = ConstantCp(init=False) + elif thermo_type == SPECIES_THERMO_SHOMATE2: + st = ShomatePoly2(init=False) + else: + st = SpeciesThermo() + + st._assign(spthermo) + return st diff --git a/interfaces/cython/cantera/test/test_thermo.py b/interfaces/cython/cantera/test/test_thermo.py index fea93948f..bab21d297 100644 --- a/interfaces/cython/cantera/test/test_thermo.py +++ b/interfaces/cython/cantera/test/test_thermo.py @@ -657,3 +657,39 @@ class TestSpecies(utilities.CanteraTest): for name in gas.species_names: s = gas.species(name) self.assertEqual(s.name, name) + + +class TestSpeciesThermo(utilities.CanteraTest): + def setUp(self): + self.gas = ct.Solution('h2o2.xml') + self.gas.X = 'H2O:1.0' + + def test_create(self): + st = ct.NasaPoly2(300, 3500, 101325, + [1000.0, 3.03399249E+00, 2.17691804E-03, -1.64072518E-07, + -9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+00, + 4.19864056E+00, -2.03643410E-03, 6.52040211E-06, -5.48797062E-09, + 1.77197817E-12, -3.02937267E+04, -8.49032208E-01]) + + for T in [300, 500, 900, 1200, 2000]: + self.gas.TP = T, 101325 + self.assertAlmostEqual(st.cp(T), self.gas.cp_mole) + self.assertAlmostEqual(st.h(T), self.gas.enthalpy_mole) + self.assertAlmostEqual(st.s(T), self.gas.entropy_mole) + + def test_invalid(self): + with self.assertRaises(ValueError): + # not enough coefficients + st = ct.NasaPoly2(300, 3500, 101325, + [1000.0, 3.03399249E+00, 2.17691804E-03]) + + def test_wrap(self): + st = self.gas.species('H2O').thermo + + self.assertTrue(isinstance(st, ct.NasaPoly2)) + + for T in [300, 500, 900, 1200, 2000]: + self.gas.TP = T, 101325 + self.assertAlmostEqual(st.cp(T), self.gas.cp_mole) + self.assertAlmostEqual(st.h(T), self.gas.enthalpy_mole) + self.assertAlmostEqual(st.s(T), self.gas.entropy_mole) diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index 57df8bb0a..c89e390d6 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -54,6 +54,13 @@ cdef class Species: def __get__(self): return self.species.size + property thermo: + def __get__(self): + return wrapSpeciesThermo(self.species.thermo) + + def __set__(self, SpeciesThermo spthermo): + self.species.thermo = spthermo._spthermo + cdef class ThermoPhase(_SolutionBase): """