From f699748d0e04ac5d69a623e688f6cd881dec5c23 Mon Sep 17 00:00:00 2001 From: Santosh Shanbhogue Date: Fri, 13 May 2016 14:42:48 -0400 Subject: [PATCH] [1D/Python] Add a Twin Premixed Counterflow Flame class and example Resolves #340 --- AUTHORS | 1 + .../onedim/premixed_counterflow_twin_flame.py | 64 ++++++++++++++ interfaces/cython/cantera/onedim.py | 84 +++++++++++++++++++ 3 files changed, 149 insertions(+) create mode 100644 interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py diff --git a/AUTHORS b/AUTHORS index e85e5409a..b34023b0b 100644 --- a/AUTHORS +++ b/AUTHORS @@ -14,4 +14,5 @@ David Fronczek John Hewson Nicholas Malaya Andreas Rücker +Santosh Shanbhogue Bryan Weber diff --git a/interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py b/interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py new file mode 100644 index 000000000..36f23118f --- /dev/null +++ b/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) diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index dd8774a7a..370146a3a 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/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])