2528df0f75
These changes make it unnecessary to copy header files around during the build process, which tends to confuse IDEs and debuggers. The headers which comprise Cantera's external C++ interface are now in the 'include' directory. All of the samples and demos are now in the 'samples' subdirectory.
241 lines
8.5 KiB
Python
241 lines
8.5 KiB
Python
"""A counterflow flame."""
|
|
|
|
from onedim import *
|
|
from Cantera.num import zeros
|
|
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):
|
|
"""The domains are [
|
|
self.fuel_inlet -- class Inlet,
|
|
self.flame -- class AxisymmetricFlow,
|
|
self.oxidizer_inlet -- class Inlet
|
|
]
|
|
"""
|
|
|
|
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 fuel species
|
|
must be specified, and the oxidizer may be
|
|
|
|
>>> f.init(fuel = 'CH4')
|
|
|
|
The initial guess is generated by assuming infinitely-fast
|
|
chemistry."""
|
|
self.getInitialSoln()
|
|
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 = zeros(nsp, 'd')
|
|
yin_o = zeros(nsp, 'd')
|
|
yst = 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 = zeros([nz,nsp],'d')
|
|
t = 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)
|
|
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
|
|
|
|
t[0] = t0f
|
|
t[-1] = t0o
|
|
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):
|
|
"""Solve the flame.
|
|
loglevel -- integer flag controlling the amount of
|
|
diagnostic output. Zero suppresses all output, and
|
|
5 produces very verbose output. Default: 1
|
|
refine_grid -- if non-zero, enable grid refinement."""
|
|
|
|
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):
|
|
"""Set the criteria used to refine the flame.
|
|
ratio -- additional points will be added if the ratio of the spacing
|
|
on either side of a grid point exceeds this value
|
|
slope -- maximum difference in value between two adjacent points,
|
|
scaled by the maximum difference in the profile
|
|
(0.0 < slope < 1.0). Adds points in regions of high slope.
|
|
curve -- maximum difference in slope between two adjacent intervals,
|
|
scaled by the maximum difference in the profile
|
|
(0.0 < curve < 1.0). Adds points in regions of high
|
|
curvature.
|
|
prune -- if the slope or curve criteria are satisfied to the level of
|
|
'prune', the grid point is assumed not to be needed and is
|
|
removed. Set prune significantly smaller than
|
|
'slope' and 'curve'. Set to zero to disable pruning
|
|
the grid.
|
|
|
|
>>> f.setRefineCriteria(ratio = 5.0, slope = 0.2, curve = 0.3,
|
|
... prune = 0.03)
|
|
"""
|
|
Stack.setRefineCriteria(self, domain = self.flame,
|
|
ratio = ratio, slope = slope, curve = curve,
|
|
prune = prune)
|
|
|
|
def setProfile(self, component, locs, vals):
|
|
"""Set a profile in the flame"""
|
|
self._initialized = 1
|
|
Stack.setProfile(self, self.flame, component, locs, vals)
|
|
|
|
def set(self, tol = None, energy = '', tol_time = None):
|
|
"""Set parameters.
|
|
tol -- (rtol, atol) for steady-state
|
|
tol_time -- (rtol, atol) for time stepping
|
|
energy -- 'on' or 'off' to enable or disable the energy equation
|
|
"""
|
|
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):
|
|
"""Set the state of the object representing the gas to the
|
|
current solution at grid point j."""
|
|
nsp = self.gas.nSpecies()
|
|
y = 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)
|
|
|
|
|
|
|
|
|
|
|