[1D/Python] Add a Twin Premixed Counterflow Flame class and example

Resolves #340

AUTHORS

@@ -14,4 +14,5 @@ David Fronczek
 John Hewson
 Nicholas Malaya
 Andreas Rücker
+Santosh Shanbhogue
 Bryan Weber

interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py

@@ -0,0 +1,64 @@
+# coding: utf-8
+
+"""
+Simulate two counter-flow jets of reactants shooting into each other. This
+simulation differs from the similar premixed_counterflow_flame.py example as the
+latter simulates a jet of reactants shooting into products.
+"""
+
+import cantera as ct
+import numpy as np
+
+# This function is called to run the solver
+def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1,
+                      ratio=3, slope=0.15, curve=0.25, prune=0.05):
+    """
+    Execute this function to run the Oppposed Flow Simulation This function
+    takes a CounterFlowTwinPremixedFlame object as the first argument
+    """
+
+    oppFlame.reactants.mdot = massFlux
+    oppFlame.set_refine_criteria(ratio=ratio, slope=slope, curve=curve, prune=prune)
+
+    oppFlame.show_solution()
+    oppFlame.solve(loglevel, auto=True)
+
+    # Compute the strain rate, just before the flame. It also turns out to the
+    # maximum. This is the strain rate that computations comprare against, like
+    # when plotting Su vs. K
+    peakStrain = np.max(np.gradient(oppFlame.u, np.gradient(oppFlame.grid)))
+    return np.max(oppFlame.T), peakStrain
+
+
+# Use the standard GRI3.0 Mechanism for CH4
+gas = ct.Solution('gri30.cti')
+
+# Create a CH4/Air premixed mixture with equivalence ratio=0.75, and at room
+# temperature and pressure.
+gas.set_equivalence_ratio(0.75, 'CH4', {'O2':1.0, 'N2':3.76})
+gas.TP = 300, ct.one_atm
+
+# Set the velocity
+axial_velocity = 0.25 # in m/s
+
+# Domain half-width of 2.5 cm, meaning the whole domain is 5 cm wide
+width = 0.025
+
+# Done with initial conditions
+
+# Compute the mass flux, as this is what the Flame object requires
+massFlux = gas.density * axial_velocity # units kg/m2/s
+# Create the flame object
+oppFlame = ct.CounterflowTwinPremixedFlame(gas, width=width)
+
+# Now run the solver
+
+# The solver returns the peak temperature and strain rate. You can plot/see all
+# state space variables by calling oppFlame.foo where foo is T, Y[i] or whatever
+# The spatial variable (distance in meters) is in oppFlame.grid Thus to plot
+# temperature vs distance, use oppFlame.grid and oppFlame.T
+(T, K) = solveOpposedFlame(oppFlame, massFlux)
+
+print("Peak temperature: {0}".format(T))
+print("Strain Rate: {0}".format(K))
+oppFlame.write_csv("premixed_twin_flame.csv", quiet=False)

interfaces/cython/cantera/onedim.py

@@ -952,3 +952,87 @@ class CounterflowPremixedFlame(FlameBase):
 
         self.set_profile('u', [0.0, 1.0], [uu, -ub])
         self.set_profile('V', [0.0, x0/dz, 1.0], [0.0, a, 0.0])
+
+
+class CounterflowTwinPremixedFlame(FlameBase):
+    """
+    A twin premixed counterflow flame. Two opposed jets of the same composition
+    shooting into each other.
+    """
+    __slots__ = ('reactants', 'flame', 'products')
+
+    def __init__(self, gas, grid=None, width=None):
+        """
+        :param gas:
+            `Solution` (using the IdealGas thermodynamic model) used to
+            evaluate all gas properties and reaction rates.
+        :param grid:
+            Array of initial grid points. Not recommended unless solving only on
+            a fixed grid; Use the `width` parameter instead.
+        :param width:
+            Defines a grid on the interval [0, width] with internal points
+            determined automatically by the solver.
+
+        A domain of class `AxisymmetricStagnationFlow` named ``flame`` will
+        be created to represent the flame. The three domains comprising the
+        stack are stored as ``self.reactants``, ``self.flame``, and
+        ``self.products``.
+        """
+        self.reactants = Inlet1D(name='reactants', phase=gas)
+        self.reactants.T = gas.T
+
+        self.flame = AxisymmetricStagnationFlow(gas, name='flame')
+
+        #The right boundary is a symmetry plane
+        self.products = SymmetryPlane1D(name='products', phase=gas)
+
+        if width is not None:
+            # Create grid points aligned with initial guess profile
+            grid = np.array([0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 1.0]) * width
+
+        super(CounterflowTwinPremixedFlame, self).__init__(
+                (self.reactants, self.flame, self.products), gas, grid)
+
+        # Setting X needs to be deferred until linked to the flow domain
+        self.reactants.X = gas.X
+
+    def set_initial_guess(self):
+        """
+        Set the initial guess for the solution.
+        """
+        super(CounterflowTwinPremixedFlame, self).set_initial_guess()
+
+        Yu = self.reactants.Y
+        Tu = self.reactants.T
+        self.gas.TPY = Tu, self.flame.P, Yu
+        rhou = self.gas.density
+        uu = self.reactants.mdot / rhou
+
+        self.gas.equilibrate('HP')
+        Teq = self.gas.T
+        Yeq = self.gas.Y
+
+        Tb = Teq
+        Yb = Yeq
+        self.products.T = Tb
+
+        self.gas.TPY = Tb, self.flame.P, Yb
+        rhob = self.gas.density
+        ub = self.products.mdot / rhob
+
+        locs = np.array([0.0, 0.4, 0.6, 1.0])
+        self.set_profile('T', locs, [Tu, Tu, Teq, Tb])
+        for k in range(self.gas.n_species):
+            self.set_profile(self.gas.species_name(k), locs,
+                             [Yu[k], Yu[k], Yeq[k], Yb[k]])
+
+        # estimate strain rate
+        self.gas.TPY = Teq, self.flame.P, Yeq
+        zz = self.flame.grid
+        dz = zz[-1] - zz[0]
+        a = (uu + ub)/dz
+        # estimate stagnation point
+        x0 = rhou*uu * dz / (rhou*uu + rhob*ub)
+
+        self.set_profile('u', [0.0, 1.0], [uu, -ub])
+        self.set_profile('V', [0.0, x0/dz, 1.0], [0.0, a, 0.0])