[Python] Add ThermoPhase.set_equivalence_ratio

interfaces/cython/cantera/examples/multiphase/adiabatic.py

@@ -37,15 +37,9 @@ mix = ct.Mixture(mix_phases)
 tad = np.zeros(npoints)
 xeq = np.zeros((mix.n_species,npoints))
 
-if gas.n_atoms(fuel_species,'O') > 0 or gas.n_atoms(fuel_species,'N') > 0:
-    raise "Error: only hydrocarbon fuels are supported."
-
-stoich_O2 = gas.n_atoms(fuel_species,'C') + 0.25*gas.n_atoms(fuel_species,'H')
-
 for i in range(npoints):
-    X = {fuel_species: phi[i] / stoich_O2, 'O2': 1.0, 'N2': 3.76}
     # set the gas state
-    gas.TPX = T, P, X
+    gas.set_equivalence_ratio(phi[i], fuel_species, 'O2:1.0, N2:3.76')
 
     # create a mixture of 1 mole of gas, and 0 moles of solid carbon.
     mix = ct.Mixture(mix_phases)

interfaces/cython/cantera/examples/reactors/reactor2.py

@@ -37,7 +37,8 @@ env = ct.Reservoir(ct.Solution('air.xml'))
 
 # use GRI-Mech 3.0 for the methane/air mixture, and set its initial state
 gri3 = ct.Solution('gri30.xml')
-gri3.TPX = 500.0, 0.2 * ct.one_atm, 'CH4:1.1, O2:2, N2:7.52'
+gri3.TP = 500.0, 0.2 * ct.one_atm
+gri3.set_equivalence_ratio(1.1, 'CH4:1.0', 'O2:2, N2:7.52')
 
 # create a reactor for the methane/air side
 r2 = ct.IdealGasReactor(gri3)

interfaces/cython/cantera/test/test_equilibrium.py

@@ -165,8 +165,6 @@ class TestEquil_GasCarbon(utilities.CanteraTest):
         cls.carbon = ct.Solution('graphite.xml')
         cls.fuel = 'CH4'
         cls.mix_phases = [(cls.gas, 1.0), (cls.carbon, 0.0)]
-        cls.stoich = (cls.gas.n_atoms(cls.fuel,'C') +
-                       0.25*cls.gas.n_atoms(cls.fuel,'H'))
         cls.n_species = cls.gas.n_species + cls.carbon.n_species
 
     def solve(self, solver, **kwargs):
@@ -176,9 +174,8 @@ class TestEquil_GasCarbon(utilities.CanteraTest):
         data = np.zeros((n_points, 2+self.n_species))
         phi = np.linspace(0.3, 3.5, n_points)
         for i in range(n_points):
-            X = {self.fuel: phi[i] / self.stoich, 'O2': 1.0, 'N2': 3.76}
-            self.gas.TPX = T, P, X
-
+            self.gas.set_equivalence_ratio(phi[i], self.fuel,
+                                           {'O2': 1.0, 'N2': 3.76})
             mix = ct.Mixture(self.mix_phases)
             mix.T = T
             mix.P = P

interfaces/cython/cantera/test/test_thermo.py

@@ -199,6 +199,25 @@ class TestThermoPhase(utilities.CanteraTest):
         with self.assertRaises(ValueError):
             self.phase['H2','O2'].Y = [0.1, 0.2, 0.3]
 
+    def test_set_equivalence_ratio_stoichiometric(self):
+        gas = ct.Solution('gri30.xml')
+        for fuel in ('C2H6', 'H2:0.7, CO:0.3', 'NH3:0.4, CH3OH:0.6'):
+            for oxidizer in ('O2:1.0, N2:3.76', 'H2O2:1.0'):
+                gas.set_equivalence_ratio(1.0, fuel, oxidizer)
+                gas.equilibrate('TP')
+                # Almost everything should end up as CO2, H2O and N2
+                self.assertGreater(sum(gas['H2O','CO2','N2'].X), 0.999999)
+
+    def test_set_equivalence_ratio_lean(self):
+        gas = ct.Solution('gri30.xml')
+        excess = 0
+        for phi in np.linspace(0.9, 0, 5):
+            gas.set_equivalence_ratio(phi, 'CH4:0.8, CH3OH:0.2', 'O2:1.0, N2:3.76')
+            gas.equilibrate('TP')
+            self.assertGreater(gas['O2'].X[0], excess)
+            excess = gas['O2'].X[0]
+        self.assertNear(sum(gas['O2','N2'].X), 1.0)
+
     def test_full_report(self):
         report = self.phase.report(threshold=0.0)
         self.assertIn(self.phase.name, report)

interfaces/cython/cantera/thermo.pyx

@@ -530,6 +530,65 @@ cdef class ThermoPhase(_SolutionBase):
         def __set__(self, C):
             self._setArray1(thermo_setConcentrations, C)
 
+    def set_equivalence_ratio(self, phi, fuel, oxidizer):
+        """
+        Set the composition to a mixture of *fuel* and *oxidizer* at the
+        specified equivalence ratio *phi*, holding temperature and pressure
+        constant. Considers the oxidation of C and H to CO2 and H2O. Other
+        elements are assumed not to participate in oxidation (i.e. N ends up as
+        N2)::
+
+            >>> gas.set_equivalence_ratio(0.5, 'CH4', 'O2:1.0, N2:3.76')
+            >>> gas.mole_fraction_dict()
+            {'CH4': 0.049900199, 'N2': 0.750499001, 'O2': 0.199600798}
+
+            >>> gas.set_equivalence_ratio(1.2, {'NH3;:0.8, 'CO':0.2}, 'O2:1.0')
+            >>> gas.mole_fraction_dict()
+            {'CO': 0.1263157894, 'NH3': 0.505263157, 'O2': 0.36842105}
+
+        :param phi: Equivalence ratio
+        :param fuel:
+            Fuel species name or molar composition as string, array, or dict.
+        :param oxidizer:
+            Oxidizer species name or molar composition as a string, array, or
+            dict.
+        """
+        if (isinstance(fuel, str) and ':' not in fuel
+            and fuel in self.species_names):
+            fuel += ':1.0'
+
+        if (isinstance(oxidizer, str) and ':' not in oxidizer
+            and oxidizer in self.species_names):
+            oxidizer += ':1.0'
+
+        self.TPX = None, None, fuel
+        Xf = self.X
+        self.TPX = None, None, oxidizer
+        Xo = self.X
+
+        nO = np.array([self.n_atoms(k, 'O') for k in range(self.n_species)])
+
+        if 'C' in self.element_names:
+            nC = np.array([self.n_atoms(k, 'C') for k in range(self.n_species)])
+        else:
+            nC = np.zeros(self.n_species)
+
+        if 'H' in self.element_names:
+            nH = np.array([self.n_atoms(k, 'H') for k in range(self.n_species)])
+        else:
+            nH = np.zeros(self.n_species)
+
+        Cf = nC.dot(Xf)
+        Co = nC.dot(Xo)
+        Of = nO.dot(Xf)
+        Oo = nO.dot(Xo)
+        Hf = nH.dot(Xf)
+        Ho = nH.dot(Xo)
+
+        stoichAirFuelRatio = - (Of - 2*Cf - Hf/2.0) / (Oo - 2*Co - Ho/2.0)
+        Xr = phi * Xf + stoichAirFuelRatio * Xo
+        self.TPX = None, None, Xr
+
     def elemental_mass_fraction(self, m):
         r"""
         Get the elemental mass fraction :math:`Z_{\mathrm{mass},m}` of element