From 3132e39e1f55dca427d4dd6ebd4022489ee22a6e Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 9 Sep 2003 16:15:11 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/python/Cantera/Kinetics.py | 9 +- .../Cantera/OneD/{flame.py => BurnerFlame.py} | 8 +- Cantera/python/Cantera/OneD/CounterFlame.py | 196 ++++++++++++++++++ Cantera/python/Cantera/OneD/StagnationFlow.py | 104 ++++++++++ Cantera/python/Cantera/OneD/__init__.py | 4 + Cantera/python/Cantera/OneD/onedim.py | 24 ++- 6 files changed, 335 insertions(+), 10 deletions(-) rename Cantera/python/Cantera/OneD/{flame.py => BurnerFlame.py} (90%) create mode 100644 Cantera/python/Cantera/OneD/CounterFlame.py create mode 100644 Cantera/python/Cantera/OneD/StagnationFlow.py diff --git a/Cantera/python/Cantera/Kinetics.py b/Cantera/python/Cantera/Kinetics.py index d0d027367..87832d132 100755 --- a/Cantera/python/Cantera/Kinetics.py +++ b/Cantera/python/Cantera/Kinetics.py @@ -237,8 +237,13 @@ class Kinetics: def multiplier(self,i): return _cantera.kin_multiplier(self.ckin,i) - def setMultiplier(self,i,v): - return _cantera.kin_setMultiplier(self.ckin,i,v) + def setMultiplier(self, value = 0.0, reaction = -1): + if reaction < 0: + nr = self.nReactions() + for i in range(nr): + _cantera.kin_setMultiplier(self.ckin,i,value) + else: + _cantera.kin_setMultiplier(self.ckin,reaction,value) def advanceCoverages(self,dt): return _cantera.kin_advanceCoverages(self.ckin,dt) diff --git a/Cantera/python/Cantera/OneD/flame.py b/Cantera/python/Cantera/OneD/BurnerFlame.py similarity index 90% rename from Cantera/python/Cantera/OneD/flame.py rename to Cantera/python/Cantera/OneD/BurnerFlame.py index d98fb807f..5dc2d474c 100644 --- a/Cantera/python/Cantera/OneD/flame.py +++ b/Cantera/python/Cantera/OneD/BurnerFlame.py @@ -54,8 +54,10 @@ class BurnerFlame(Stack): Stack.solve(self, loglevel = loglevel, refine_grid = refine_grid) - def setRefineCriteria(self, ratio = 10.0, slope = 0.8, curve = 0.8, prune = 0.0): - Stack.setRefineCriteria(self, domain = self.flame, ratio = ratio, slope = slope, curve = curve, + def setRefineCriteria(self, ratio = 10.0, slope = 0.8, + curve = 0.8, prune = 0.0): + Stack.setRefineCriteria(self, domain = self.flame, + ratio = ratio, slope = slope, curve = curve, prune = prune) def setProfile(self, component, locs, vals): @@ -89,7 +91,7 @@ class BurnerFlame(Stack): for n in range(nsp): nm = self.gas.speciesName(n) y[n] = self.solution(nm, j) - self.gas.setState_TPY(self.T(j), self.flame.pressure(), y) + self.gas.setState_TPY(self.T(j), self.pressure, y) diff --git a/Cantera/python/Cantera/OneD/CounterFlame.py b/Cantera/python/Cantera/OneD/CounterFlame.py new file mode 100644 index 000000000..9e9301638 --- /dev/null +++ b/Cantera/python/Cantera/OneD/CounterFlame.py @@ -0,0 +1,196 @@ +"""A counterflow flame.""" + +from onedim import * +import Numeric +import math + +def erfc(x): + """The complementary error function.""" + exp = math.exp + + p = 0.3275911 + a1 = 0.254829592 + a2 = -0.284496736 + a3 = 1.421413741 + a4 = -1.453152027 + a5 = 1.061405429 + + t = 1.0 / (1.0 + p*x) + erfcx = ( (a1 + (a2 + (a3 + + (a4 + a5*t)*t)*t)*t)*t ) * exp(-x*x) + return erfcx + +def erf(x): + """The error function.""" + if x < 0: + return -(1.0 - erfc(-x)) + else: + return 1.0 - erfc(x) + + +class CounterFlame(Stack): + """A non-premixed counterflow flame.""" + + def __init__(self, gas = None, grid = None): + self.fuel_inlet = Inlet('fuel inlet') + self.oxidizer_inlet = Inlet('oxidizer inlet') + self.gas = gas + self.fuel_inlet.set(temperature = gas.temperature()) + self.oxidizer_inlet.set(temperature = gas.temperature()) + self.pressure = gas.pressure() + self.flame = AxisymmetricFlow('flame',gas = gas) + self.flame.setupGrid(grid) + Stack.__init__(self, [self.fuel_inlet, self.flame, + self.oxidizer_inlet]) + self.setRefineCriteria() + self._initialized = 0 + + + def init(self, fuel = '', oxidizer = 'O2', stoich = -1.0): + """Set the initial guess for the solution. + + The initial guess is generated by assuming infinitely-fast + chemistry.""" + + gas = self.gas + nsp = gas.nSpecies() + wt = gas.molecularWeights() + + # find the fuel and oxidizer species + iox = gas.speciesIndex(oxidizer) + ifuel = gas.speciesIndex(fuel) + + # if no stoichiometric ratio was input, compute it + if stoich < 0.0: + if oxidizer == 'O2': + nh = gas.nAtoms(fuel, 'H') + nc = gas.nAtoms(fuel, 'C') + stoich = 1.0*nc + 0.25*nh + else: + raise CanteraError('oxidizer/fuel stoichiometric ratio must'+ + ' be specified, since the oxidizer is not O2') + + s = stoich*wt[iox]/wt[ifuel] + y0f = self.fuel_inlet.massFraction(ifuel) + y0ox = self.oxidizer_inlet.massFraction(iox) + phi = s*y0f/y0ox + zst = 1.0/(1.0 + phi) + + yin_f = Numeric.zeros(nsp, 'd') + yin_o = Numeric.zeros(nsp, 'd') + yst = Numeric.zeros(nsp, 'd') + for k in range(nsp): + yin_f[k] = self.fuel_inlet.massFraction(k) + yin_o[k] = self.oxidizer_inlet.massFraction(k) + yst[k] = zst*yin_f[k] + (1.0 - zst)*yin_o[k] + + gas.setState_TPY(self.fuel_inlet.temperature(), self.pressure, yin_f) + mdotf = self.fuel_inlet.mdot() + u0f = mdotf/gas.density() + t0f = self.fuel_inlet.temperature() + + gas.setState_TPY(self.oxidizer_inlet.temperature(), + self.pressure, yin_o) + mdoto = self.oxidizer_inlet.mdot() + u0o = mdoto/gas.density() + t0o = self.oxidizer_inlet.temperature() + + # get adiabatic flame temperature and composition + tbar = 0.5*(t0o + t0f) + gas.setState_TPY(tbar, self.pressure, yst) + gas.equilibrate('HP') + teq = gas.temperature() + yeq = gas.massFractions() + + # estimate strain rate + zz = self.flame.grid() + dz = zz[-1] - zz[0] + a = (u0o + u0f)/dz + diff = gas.mixDiffCoeffs() + f = math.sqrt(a/(2.0*diff[iox])) + + x0 = mdotf*dz/(mdotf + mdoto) + nz = len(zz) + + y = Numeric.zeros([nz,nsp],'d') + t = Numeric.zeros(nz,'d') + for j in range(nz): + x = zz[j] + zeta = f*(x - x0) + zmix = 0.5*(1.0 - erf(zeta)) + if zmix > zst: + for k in range(nsp): + y[j,k] = yeq[k] + (zmix - zst)*(yin_f[k] + - yeq[k])/(1.0 - zst) + t[j] = teq + (t0f - teq)*(zmix - zst)/(1.0 - zst) + print teq, t[j], t0f, zmix, zst + else: + for k in range(nsp): + y[j,k] = yin_o[k] + zmix*(yeq[k] - yin_o[k])/zst + t[j] = t0o + (teq - t0o)*zmix/zst + + zrel = zz/dz + self.setProfile('u', [0.0, 1.0], [u0f, -u0o]) + self.setProfile('V', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) + self.setProfile('T', zrel, t) + for k in range(nsp): + self.setProfile(gas.speciesName(k), zrel, y[:,k]) + + self._initialized = 1 + + + def solve(self, loglevel = 1, refine_grid = 1): + if not self._initialized: self.init() + Stack.solve(self, loglevel = loglevel, refine_grid = refine_grid) + + + def setRefineCriteria(self, ratio = 10.0, slope = 0.8, curve = 0.8, + prune = 0.0): + Stack.setRefineCriteria(self, domain = self.flame, + ratio = ratio, slope = slope, curve = curve, + prune = prune) + + def setProfile(self, component, locs, vals): + self._initialized = 1 + Stack.setProfile(self, self.flame, component, locs, vals) + + def set(self, tol = None, energy = '', tol_time = None): + if tol: + self.flame.setTolerances(default = tol) + if tol_time: + self.flame.setTolerances(default = tol_time, time = 1) + if energy: + self.flame.set(energy = energy) + + def T(self, point = -1): + """The temperature [K]""" + return self.solution('T', point) + + def u(self, point = -1): + """The axial velocity [m/s]""" + return self.solution('u', point) + + def V(self, point = -1): + """The radial velocity divided by radius [s^-1]""" + return self.solution('V', point) + + def solution(self, component = '', point = -1): + """The solution for one specified component. If a point number + is given, return the value of component 'component' at this + point. Otherwise, return the entire profile for this + component.""" + if point >= 0: return self.value(self.flame, component, point) + else: return self.profile(self.flame, component) + + def setGasState(self, j): + nsp = self.gas.nSpecies() + y = Numeric.zeros(nsp, 'd') + for n in range(nsp): + nm = self.gas.speciesName(n) + y[n] = self.solution(nm, j) + self.gas.setState_TPY(self.T(j), self.pressure, y) + + + + + diff --git a/Cantera/python/Cantera/OneD/StagnationFlow.py b/Cantera/python/Cantera/OneD/StagnationFlow.py new file mode 100644 index 000000000..25941a47d --- /dev/null +++ b/Cantera/python/Cantera/OneD/StagnationFlow.py @@ -0,0 +1,104 @@ +from onedim import * +import Numeric + +class StagnationFlow(Stack): + """An axisymmetric flow impinging on a surface at normal incidence.""" + + def __init__(self, gas = None, surfchem = None, grid = None): + self.inlet = Inlet('inlet') + self.gas = gas + self.surfchem = surfchem + self.inlet.set(temperature = gas.temperature()) + self.surface = Surface(id = 'surface', surface_mech = surfchem) + self.pressure = gas.pressure() + self.flow = AxisymmetricFlow('flow',gas = gas) + self.flow.setupGrid(grid) + Stack.__init__(self, [self.inlet, self.flow, self.surface]) + self.setRefineCriteria() + self._initialized = 0 + + + def init(self): + """Set the initial guess for the solution.""" + self.getInitialSoln() + gas = self.gas + nsp = gas.nSpecies() + yin = Numeric.zeros(nsp, 'd') + for k in range(nsp): + yin[k] = self.inlet.massFraction(k) + gas.setState_TPY(self.inlet.temperature(), self.pressure, yin) + u0 = self.inlet.mdot()/gas.density() + t0 = self.inlet.temperature() + V0 = 0.0 + + tsurf = self.surface.temperature() + + zz = self.flow.grid() + dz = zz[-1] - zz[0] + + + locs = Numeric.array([0.0, 1.0],'d') + self.setProfile('u', locs, [u0, 0.0]) + self.setProfile('V', locs, [V0, V0]) + self.setProfile('T', locs, [t0, tsurf]) + for n in range(nsp): + self.setProfile(gas.speciesName(n), locs, [yin[n], yin[n]]) + self._initialized = 1 + + + def solve(self, loglevel = 1, refine_grid = 1): + if not self._initialized: self.init() + Stack.solve(self, loglevel = loglevel, refine_grid = refine_grid) + + + def setRefineCriteria(self, ratio = 10.0, slope = 0.8, + curve = 0.8, prune = 0.0): + Stack.setRefineCriteria(self, domain = self.flow, + ratio = ratio, slope = slope, curve = curve, + prune = prune) + + def setProfile(self, component, locs, vals): + self._initialized = 1 + Stack.setProfile(self, self.flow, component, locs, vals) + + def set(self, tol = None, energy = '', tol_time = None): + if tol: + self.flow.setTolerances(default = tol) + if tol_time: + self.flow.setTolerances(default = tol_time, time = 1) + if energy: + self.flow.set(energy = energy) + + def T(self, point = -1): + return self.solution('T', point) + + def u(self, point = -1): + return self.solution('u', point) + + def V(self, point = -1): + return self.solution('V', point) + + def solution(self, component = '', point = -1): + if point >= 0: return self.value(self.flow, component, point) + else: return self.profile(self.flow, component) + + def coverages(self): + nsurf = self.surfchem.nSpecies() + cov = Numeric.zeros(nsurf,'d') + for n in range(nsurf): + nm = self.surfchem.speciesName(n) + cov[n] = self.value(self.surface, nm, 0) + return cov + + def setGasState(self, j): + nsp = self.gas.nSpecies() + y = Numeric.zeros(nsp, 'd') + for n in range(nsp): + nm = self.gas.speciesName(n) + y[n] = self.solution(nm, j) + self.gas.setState_TPY(self.T(j), self.pressure, y) + + + + + diff --git a/Cantera/python/Cantera/OneD/__init__.py b/Cantera/python/Cantera/OneD/__init__.py index 9f4e4b1e9..6f59867e4 100644 --- a/Cantera/python/Cantera/OneD/__init__.py +++ b/Cantera/python/Cantera/OneD/__init__.py @@ -2,4 +2,8 @@ The classes in this package implement one-dimensional reacting flow problems. """ from onedim import * +from BurnerFlame import BurnerFlame +from CounterFlame import CounterFlame +from StagnationFlow import StagnationFlow + diff --git a/Cantera/python/Cantera/OneD/onedim.py b/Cantera/python/Cantera/OneD/onedim.py index a69f26e23..c51b92861 100644 --- a/Cantera/python/Cantera/OneD/onedim.py +++ b/Cantera/python/Cantera/OneD/onedim.py @@ -22,7 +22,8 @@ class Domain1D: return _cantera.domain_type(self._hndl) def index(self): - """Index of this domain in a stack. Returns -1 if this domain is not part of a stack.""" + """Index of this domain in a stack. Returns -1 if this domain + is not part of a stack.""" return _cantera.domain_index(self._hndl) def nComponents(self): @@ -36,6 +37,12 @@ class Domain1D: def componentName(self, n): """Name of the nth component.""" return _cantera.domain_componentName(self._hndl, n) + + def componentNames(self): + names = [] + for n in range(self.nComponents()): + names.append(self.componentName(n)) + return names def componentIndex(self, name): """Index of the component with name 'name'""" @@ -248,6 +255,7 @@ class Surface(Bdry1D): Bdry1D.__init__(self) if surface_mech: self._hndl = _cantera.reactingsurf_new() + self.setKineticsMgr(surface_mech) else: self._hndl = _cantera.surf_new() if id: self.setID(id) @@ -256,9 +264,12 @@ class Surface(Bdry1D): def setKineticsMgr(self, kin): _cantera.reactingsurf_setkineticsmgr(self._hndl, kin.kinetics_hndl()) - def enableCoverageEqs(self, onoff=1): - _cantera.reactingsurf_enableCoverageEqs(self._hndl, onoff) - + + def setCoverageEqs(self, onoff='on'): + if onoff == 'on': + _cantera.reactingsurf_enableCoverageEqs(self._hndl, 1) + else: + _cantera.reactingsurf_enableCoverageEqs(self._hndl, 0) class AxisymmetricFlow(Domain1D): @@ -382,6 +393,9 @@ class Stack: _cantera.sim1D_setTimeStep(self._hndl, stepsize, Numeric.asarray(nsteps)) + def getInitialSoln(self): + _cantera.sim1D_getInitialSoln(self._hndl) + def solve(self, loglevel=1, refine_grid=1): return _cantera.sim1D_solve(self._hndl, loglevel, refine_grid) @@ -416,7 +430,7 @@ class Stack: for n in range(np): x[n] = self.value(domain, component, n) return x - + def workValue(self, dom, icomp, localPoint): idom = dom.index() return _cantera.sim1D_workValue(self._hndl, idom, icomp, localPoint)