Add dilution options to equivalence ratio functions

interfaces/cython/cantera/examples/thermo/equivalenceRatio.py

@@ -2,52 +2,131 @@
 This example demonstrates how to set a mixture according to equivalence ratio
 and mixture fraction.
 
-Requires: cantera >= 2.5.0
+Requires: cantera >= 2.6.0
 """
 
 import cantera as ct
 
 gas = ct.Solution('gri30.yaml')
 
-# fuel and oxidizer compositions
-fuel = "CH4"
-oxidizer = "O2:0.21,N2:0.79"
+# Define the oxidizer composition, here air with 21 mol-% O2 and 79 mol-% N2
+air = "O2:0.21,N2:0.79"
 
-gas.TP = 300, ct.one_atm
+# Set the mixture composition according to the stoichiometric mixture
+# (equivalence ratio phi = 1). The fuel composition in this example
+# is set to 100 mol-% CH4 and the oxidizer to 21 mol-% O2 and 79 mol-% N2.
+# This function changes the composition of the gas object and keeps temperature
+# and pressure constant
+gas.set_equivalence_ratio(phi=1.0, fuel="CH4:1", oxidizer=air)
 
-# set the mixture composition according to the stoichiometric mixture
-# (equivalence ratio = 1)
-gas.set_equivalence_ratio(1, fuel, oxidizer)
+# If fuel or oxidizer consist of a single species, a short hand notation can be
+# used, for example fuel="CH4" is equivalent to fuel="CH4:1".
+# By default, the compositions of fuel and oxidizer are interpreted as mole
+# fractions. If the compositions are given in mass fractions, an
+# additional argument can be provided. Here, the fuel is 100 mass-% CH4
+# and the oxidizer is 23.3 mass-% O2 and 76.7 mass-% N2
+gas.set_equivalence_ratio(1.0, fuel="CH4:1", oxidizer="O2:0.233,N2:0.767", basis='mass')
 
 # This function can be used to compute the equivalence ratio for any mixture.
-# An optional argument "basis" indicates if fuel and oxidizer compositions are
-# provided in terms of mass or mole fractions. Default is mole fractions.
-# If fuel and oxidizer are given in mass fractions, use basis='mass'
-phi = gas.equivalence_ratio(fuel, oxidizer)
-print("phi = {:1.3f}".format(phi))
-
-# The equivalence ratio can also be computed from the elemental composition
-# assuming that there is no oxygen in the fuel and no C,H and S elements
-# in the oxidizer so that the composition of fuel and oxidizer can be omitted
+# The first two arguments specify the compositions of the fuel and oxidizer.
+# An optional third argument "basis" indicates if fuel and oxidizer compositions
+# are provided in terms of mass or mole fractions. Default is mole fractions.
+# Note that for all functions shown here, the compositions are normalized
+# internally so the species fractions do not have to sum to unity
+phi = gas.equivalence_ratio(fuel="CH4:1", oxidizer="O2:233,N2:767", basis='mass')
+print(f"phi = {phi:1.3f}")
+
+# If the compositions of fuel and oxidizer are unknown, the function can
+# be called without arguments. This assumes that all C, H and S atoms come from
+# the fuel and all O atoms from the oxidizer. In this example, the fuel was set
+# to be pre CH4 and the oxidizer O2:0.233,N2:0.767 so that the assumption is true
+# and the same equivalence ratio as above is computed
 phi = gas.equivalence_ratio()
+print(f"phi = {phi:1.3f}")
+
+# Instead of working with equivalence ratio, mixture fraction can be used.
+# The mixture fraction is always kg fuel / (kg fuel + kg oxidizer), independent
+# of the basis argument. For example, the mixture fraction Z can be computed as
+# follows. Again, the compositions by default are interpreted as mole fractions
+Z = gas.mixture_fraction(fuel="CH4:1", oxidizer=air)
+print(f"Z = {Z:1.3f}")
 
-# In this example, the result is the same as above
-print("phi = {:1.3f}".format(phi))
+# By default, the mixture fraction is the Bilger mixture fraction. Instead,
+# a mixture fraction based on a single element can be used. In this example,
+# the following two ways of computing Z are the same:
+Z = gas.mixture_fraction(fuel="CH4:1", oxidizer=air, element="Bilger")
+print(f"Z(Bilger mixture fraction) = {Z:1.3f}")
+Z = gas.mixture_fraction(fuel="CH4:1", oxidizer=air, element="C")
+print(f"Z(mixture fraction based on C) = {Z:1.3f}")
 
-# the mixture fraction Z can be computed as follows:
-Z = gas.mixture_fraction(fuel, oxidizer)
-print("Z = {:1.3f}".format(Z))
+# Since the fuel in this example is pure methane and the oxidizer is air,
+# the mixture fraction is the same as the mass fraction of methane in the mixture
+print(f"mass fraction of CH4 = {gas['CH4'].Y[0]:1.3f}")
 
-# The mixture fraction is kg fuel / (kg fuel + kg oxidizer). Since the fuel in
-# this example is pure methane and the oxidizer is air, the mixture fraction
-# is the same as the mass fraction of methane in the mixture
-print("mass fraction of CH4 = {:1.3f}".format(gas["CH4"].Y[0]))
+# To set a mixture according to the mixture fraction, the following function
+# can be used. In this example, the final fuel/oxidizer mixture
+# contains 5.5 mass-% CH4:
+gas.set_mixture_fraction(0.055, fuel="CH4:1", oxidizer=air)
+print(f"mass fraction of CH4 = {gas['CH4'].Y[0]:1.3f}")
 
 # Mixture fraction and equivalence ratio are invariant to the reaction progress.
 # For example, they stay constant if the mixture composition changes to the burnt
-# state
+# state or for any intermediate state. Fuel and oxidizer compositions for all functions
+# shown in this example can be given as string, dictionary or numpy array
+fuel = {"CH4":1} # provide the fuel composition as dictionary instead of string
+gas.set_equivalence_ratio(1, fuel, air)
 gas.equilibrate('HP')
-phi_burnt = gas.equivalence_ratio(fuel, oxidizer)
-Z_burnt = gas.mixture_fraction(fuel, oxidizer)
-print("phi(burnt) = {:1.3f}".format(phi_burnt))
-print("Z(burnt) = {:1.3f}".format(Z_burnt))
+phi_burnt = gas.equivalence_ratio(fuel, air)
+Z_burnt = gas.mixture_fraction(fuel, air)
+print(f"phi(burnt) = {phi_burnt:1.3f}")
+print(f"Z(burnt) = {Z_burnt:1.3f}")
+
+# If fuel and oxidizer compositions are specified consistently, then
+# equivalence_ratio and set_equivalence_ratio are consistent as well, as
+# shown in the following example with arbitrary fuel and oxidizer compositions:
+gas.set_equivalence_ratio(2.5, fuel="CH4:1,O2:0.01,CO:0.05,N2:0.1",
+                          oxidizer="O2:0.2,N2:0.8,CO2:0.05,CH4:0.01")
+
+phi = gas.equivalence_ratio(fuel="CH4:1,O2:0.01,CO:0.05,N2:0.1",
+                            oxidizer="O2:0.2,N2:0.8,CO2:0.05,CH4:0.01")
+print(f"phi = {phi:1.3f}") # prints 2.5
+
+# Without specifying the fuel and oxidizer compositions, it is assumed that
+# all C, H and S atoms come from the fuel and all O atoms from the oxidizer,
+# which is not true for this example. Therefore, the following call gives a
+# different equivalence ratio based on that assumption
+phi = gas.equivalence_ratio()
+print(f"phi = {phi:1.3f}")
+
+# After computing the mixture composition for a certain equivalence ratio given
+# a fuel and mixture composition, the mixture can optionally be diluted. The
+# following function will first create a mixture with equivalence ratio 2 from pure
+# hydrogen and oxygen and then dilute it with H2O. In this example, the final mixture
+# consists of 30 mol-% H2O and 70 mol-% of the H2/O2 mixture at phi=2
+gas.set_equivalence_ratio(2.0, "H2:1", "O2:1", diluent="H2O", fraction={"diluent":0.3})
+print(f"mole fraction of H2O = {gas['H2O'].X[0]:1.3f}") # mixture contains 30 mol-% H2O
+print(f"ratio of H2/O2: {gas['H2'].X[0] / gas['O2'].X[0]:1.3f}") # according to phi=2
+
+# Another option is to specify the fuel or oxidizer fraction in the final mixture.
+# The following example creates a mixture with equivalence ratio 2 from pure
+# hydrogen and oxygen (same as above) and then dilutes it with a mixture of 50 mass-%
+# CO2 and 50 mass-% H2O so that the mass fraction of fuel in the final mixture is 0.1
+gas.set_equivalence_ratio(2.0, "H2", "O2", diluent="CO2:0.5,H2O:0.5",
+                          fraction={"fuel":0.1}, basis="mass")
+print(f"mole fraction of H2 = {gas['H2'].Y[0]:1.3f}") # mixture contains 10 mass-% fuel
+
+# To compute the equivalence ratio given a diluted mixture, a list of
+# species names can be provided which will be considered for computing phi.
+# In this example, the diluents H2O and CO2 are ignored and only H2 and O2 are
+# considered to get the equivalence ratio
+phi = gas.equivalence_ratio(fuel="H2", oxidizer="O2", include_species=["H2", "O2"])
+print(f"phi = {phi:1.3f}") # prints 2
+
+# If instead the diluent should be included in the computation of the equivalence ratio,
+# the mixture can be set in the following way. Assume the fuel is diluted with
+# 50 mol-% H2O:
+gas.set_equivalence_ratio(2.0, fuel="H2:0.5,H2O:0.5", oxidizer=air)
+
+# This creates a mixture with the specified equivalence ratio including the diluent:
+phi = gas.equivalence_ratio(fuel="H2:0.5,H2O:0.5", oxidizer=air)
+print(f"phi = {phi:1.3f}") # prints 2

interfaces/cython/cantera/test/test_thermo.py

@@ -422,6 +422,158 @@ def test_equil_results(gas, fuel, ox, Y_Cf, Y_Of, Y_Co, Y_Oo, basis):
 
         test_equil_results(gas, fuel, ox, Y_Cf, Y_Of, Y_Co, Y_Oo, 'mass')
 
+    def test_equivalence_ratio_simple_dilution(self):
+        gas = ct.Solution("gri30.yaml")
+
+        phi = 2
+        inv_afr = 2 * phi # inverse molar AFR for H2/O2
+        X = "H2:4,O2:1,CO:3,CO2:4,N2:5,CH4:6"
+        T, P = 300, 1e5
+        gas.TPX = T, P, X
+        original_X = gas.X
+        self.assertNear(gas.equivalence_ratio(include_species=["H2", "O2"]), 2)
+        self.assertNear(gas.equivalence_ratio("H2", "O2",
+                                              include_species=["H2", "O2"]), 2)
+        self.assertNear(gas.T, T)
+        self.assertNear(gas.P, P)
+        self.assertArrayNear(gas.X, original_X)
+
+        def test_simple_dilution(fraction, basis):
+            if isinstance(fraction, str):
+                fraction_dict = {fraction[:fraction.find(":")]:
+                                 float(fraction[fraction.find(":")+1:])}
+            else:
+                fraction_dict = fraction
+
+            fraction_type  = list(fraction_dict.keys())[0]
+            fraction_value = float(list(fraction_dict.values())[0])
+
+            M_H2 = gas.molecular_weights[gas.species_index("H2")]
+            M_O2 = gas.molecular_weights[gas.species_index("O2")]
+
+            gas.TP = T, P
+            gas.set_equivalence_ratio(phi, "H2", "O2", fraction=fraction,
+                                      diluent="CO2", basis=basis)
+            if basis == "mole" and fraction_type == "diluent":
+                self.assertNear(gas["H2"].X[0], (1 - fraction_value)
+                                * inv_afr / (inv_afr + 1))
+                self.assertNear(gas["O2"].X[0], (1 - fraction_value) / (inv_afr + 1))
+                self.assertNear(gas["CO2"].X[0], fraction_value)
+            elif basis == "mass" and fraction_type == "diluent":
+                self.assertNear(gas["H2"].Y[0] / gas["O2"].Y[0], inv_afr * M_H2 / M_O2)
+                self.assertNear(gas["CO2"].Y[0], fraction_value)
+            elif basis == "mole" and fraction_type == "fuel":
+                self.assertNear(gas["H2"].X[0], fraction_value)
+                self.assertNear(gas["O2"].X[0], fraction_value / inv_afr)
+                self.assertNear(gas["CO2"].X[0], 1 - fraction_value * (1 + 1 / inv_afr))
+            elif basis == "mass" and fraction_type == "fuel":
+                self.assertNear(gas["H2"].Y[0], fraction_value)
+                self.assertNear(gas["H2"].Y[0] / gas["O2"].Y[0], inv_afr * M_H2 / M_O2)
+            elif basis == "mole" and fraction_type == "oxidizer":
+                self.assertNear(gas["H2"].X[0], fraction_value * inv_afr)
+                self.assertNear(gas["O2"].X[0], fraction_value)
+                self.assertNear(gas["CO2"].X[0], 1 - fraction_value * (1 + inv_afr))
+            elif basis == "mass" and fraction_type == "oxidizer":
+                self.assertNear(gas["O2"].Y[0], fraction_value)
+                self.assertNear(gas["H2"].Y[0] / gas["O2"].Y[0], inv_afr * M_H2 / M_O2)
+
+            Y = gas.Y
+            self.assertNear(gas.equivalence_ratio("H2", "O2",
+                                                  include_species=["H2", "O2"],
+                                                  basis=basis), phi)
+            self.assertNear(gas.equivalence_ratio(include_species=["H2", "O2"],
+                                                  basis=basis), phi)
+            self.assertArrayNear(Y, gas.Y)
+            self.assertNear(gas.T, T)
+            self.assertNear(gas.P, P)
+
+        # brute force all possible input combinations
+        test_simple_dilution("diluent:0.3", "mole")
+        test_simple_dilution({"diluent": 0.3}, "mole")
+        test_simple_dilution("diluent:0.3", "mass")
+        test_simple_dilution({"diluent": 0.3}, "mass")
+
+        test_simple_dilution("fuel:0.1", "mole")
+        test_simple_dilution({"fuel": 0.1}, "mole")
+        test_simple_dilution("fuel:0.1", "mass")
+        test_simple_dilution({"fuel": 0.1}, "mass")
+
+        test_simple_dilution("oxidizer:0.1", "mole")
+        test_simple_dilution({"oxidizer": 0.1}, "mole")
+        test_simple_dilution("oxidizer:0.1", "mass")
+        test_simple_dilution({"oxidizer": 0.1}, "mass")
+
+    def test_equivalence_ratio_arbitrary_dilution(self):
+        fuel = "CH4:1,O2:0.01,N2:0.1,CO:0.05,CO2:0.02"
+        oxidizer = "O2:0.8,N2:0.2,CO:0.01,CH4:0.005,CO2:0.03"
+        diluent = "CO2:1,CO:0.025,N2:0.3,CH4:0.07,O2:0.09"
+        gas = ct.Solution("gri30.yaml")
+
+        gas.X = fuel
+        X_fuel = gas.X
+        gas.Y = fuel
+        Y_fuel = gas.Y
+        gas.X = oxidizer
+        X_oxidizer = gas.X
+        gas.Y = oxidizer
+        Y_oxidizer = gas.Y
+        gas.X = diluent
+        X_diluent = gas.X
+        gas.Y = diluent
+        Y_diluent = gas.Y
+
+        phi = 2
+        fraction = 0.6
+        gas.set_equivalence_ratio(phi, fuel, oxidizer)
+        X_Mix = gas.X
+        gas.set_equivalence_ratio(phi, fuel, oxidizer, fraction={"diluent": fraction},
+                                  diluent=diluent)
+        X_expected = X_Mix * (1 - fraction) + fraction * X_diluent
+        self.assertArrayNear(gas.X, X_expected)
+
+        gas.set_equivalence_ratio(phi, fuel, oxidizer, basis="mass")
+        Y_Mix = gas.Y
+        gas.set_equivalence_ratio(phi, fuel, oxidizer, fraction={"diluent": fraction},
+                                  diluent=diluent, basis="mass")
+        Y_expected = Y_Mix * (1 - fraction) + fraction * Y_diluent
+        self.assertArrayNear(gas.Y, Y_expected)
+
+        phi = 0.8
+        fraction = 0.05
+        gas.set_equivalence_ratio(phi, fuel, oxidizer, basis="mass")
+        AFR = gas.stoich_air_fuel_ratio(fuel, oxidizer, basis="mass") / phi
+        gas.set_equivalence_ratio(phi, fuel, oxidizer, fraction={"fuel": fraction},
+                                  diluent=diluent, basis="mass")
+        Y_expected = fraction * (Y_fuel + AFR * Y_oxidizer) \
+                     + (1 - fraction * (1 + AFR)) * Y_diluent
+        self.assertArrayNear(gas.Y, Y_expected)
+
+        gas.set_equivalence_ratio(phi, fuel, oxidizer, fraction={"oxidizer": fraction},
+                                  diluent=diluent, basis="mass")
+        Y_expected = fraction * (Y_fuel / AFR + Y_oxidizer) \
+                     + (1 - fraction * (1 + 1 / AFR)) * Y_diluent
+        self.assertArrayNear(gas.Y, Y_expected)
+
+        gas.X = fuel
+        M_fuel = gas.mean_molecular_weight
+        gas.X = oxidizer
+        M_oxidizer = gas.mean_molecular_weight
+
+        gas.set_equivalence_ratio(phi, fuel, oxidizer)
+        AFR = M_fuel / M_oxidizer * gas.stoich_air_fuel_ratio(fuel, oxidizer) / phi
+
+        gas.set_equivalence_ratio(phi, fuel, oxidizer, fraction={"fuel": fraction},
+                                  diluent=diluent)
+        X_expected = fraction * (X_fuel + AFR * X_oxidizer) \
+                     + (1 - fraction * (1 + AFR)) * X_diluent
+        self.assertArrayNear(gas.X, X_expected)
+
+        gas.set_equivalence_ratio(phi, fuel, oxidizer, fraction={"oxidizer": fraction},
+                                  diluent=diluent)
+        X_expected = fraction * (X_fuel / AFR + X_oxidizer) \
+                     + (1 - fraction * (1 + 1 / AFR)) * X_diluent
+        self.assertArrayNear(gas.X, X_expected)
+
     def test_full_report(self):
         report = self.phase.report(threshold=0.0)
         self.assertIn(self.phase.name, report)