From 829c9b38cfce0a7723f3ad626099b1cf2a4b5b63 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 10 Jul 2014 22:34:12 +0000 Subject: [PATCH] [Python] Move some classes from Cython to pure Python The classes implementing specific flame geometries don't directly interact with any of the underlying C++ implementation, so they can be moved out of the compiled Cython extension. This reduces the size of the compiled extension, and makes it easier to implement additional flame types in Python. --- interfaces/cython/cantera/__init__.py | 1 + interfaces/cython/cantera/onedim.py | 498 ++++++++++++++++++++++++++ interfaces/cython/cantera/onedim.pyx | 491 ------------------------- 3 files changed, 499 insertions(+), 491 deletions(-) create mode 100644 interfaces/cython/cantera/onedim.py diff --git a/interfaces/cython/cantera/__init__.py b/interfaces/cython/cantera/__init__.py index d1dfcb99a..491e8ece9 100644 --- a/interfaces/cython/cantera/__init__.py +++ b/interfaces/cython/cantera/__init__.py @@ -1,6 +1,7 @@ from ._cantera import * from ._cantera import __version__, _have_sundials from .liquidvapor import * +from .onedim import * from .utils import * import os as _os diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py new file mode 100644 index 000000000..bb9e50ade --- /dev/null +++ b/interfaces/cython/cantera/onedim.py @@ -0,0 +1,498 @@ +import numpy as np +from ._cantera import * + +try: + # Python 2.7 or 3.2+ + from math import erf +except ImportError: + from scipy.special import erf + + +class FlameBase(Sim1D): + """ Base class for flames with a single flow domain """ + + def __init__(self, domains, gas, grid=None): + """ + :param gas: + object to use to evaluate all gas properties and reaction rates + :param grid: + array of initial grid points + """ + if grid is None: + grid = np.linspace(0.0, 0.1, 6) + self.flame.grid = grid + super(FlameBase, self).__init__(domains) + self.gas = gas + self.flame.P = gas.P + + def set_refine_criteria(self, ratio=10.0, slope=0.8, curve=0.8, prune=0.0): + super(FlameBase, self).set_refine_criteria(self.flame, ratio, slope, + curve, prune) + + def set_profile(self, component, locations, values): + super(FlameBase, self).set_profile(self.flame, component, locations, + values) + + @property + def transport_model(self): + return self.gas.transport_model + + @transport_model.setter + def transport_model(self, model): + self.gas.transport_model = model + self.flame.set_transport(self.gas) + + @property + def energy_enabled(self): + return self.flame.energy_enabled + + @energy_enabled.setter + def energy_enabled(self, enable): + self.flame.energy_enabled = enable + + @property + def soret_enabled(self): + return self.flame.soret_enabled + + @soret_enabled.setter + def soret_enabled(self, enable): + self.flame.soret_enabled = enable + + @property + def grid(self): + """ Array of grid point positions along the flame. """ + return self.flame.grid + + @property + def P(self): + return self.flame.P + + @P.setter + def P(self, P): + self.flame.P = P + + @property + def T(self): + """ Array containing the temperature [K] at each grid point. """ + return self.profile(self.flame, 'T') + + @property + def u(self): + """ + Array containing the velocity [m/s] normal to the flame at each point. + """ + return self.profile(self.flame, 'u') + + @property + def V(self): + """ + Array containing the tangential velocity gradient [1/s] at each point. + """ + return self.profile(self.flame, 'V') + + @property + def L(self): + """ + Array containing the radial pressure gradient (1/r)(dP/dr) [N/m^4] at + each point. Note: This value is named 'lambda' in the C++ code. + """ + return self.profile(self.flame, 'lambda') + + def solution(self, component, point=None): + if point is None: + return self.profile(self.flame, component) + else: + return self.value(self.flame, component, point) + + def set_gas_state(self, point): + k0 = self.flame.component_index(self.gas.species_name(0)) + Y = [self.solution(k, point) + for k in range(k0, k0 + self.gas.n_species)] + self.gas.TPY = self.value(self.flame, 'T', point), self.P, Y + + def write_csv(self, filename, species='X', quiet=True): + """ + Write the velocity, temperature, density, and species profiles + to a CSV file. + + :param filename: + Output file name + :param species: + Attribute to use obtaining species profiles, e.g. ``X`` for + mole fractions or ``Y`` for mass fractions. + """ + + z = self.grid + T = self.T + u = self.u + V = self.V + + csvfile = open(filename, 'w') + writer = csv.writer(csvfile) + writer.writerow(['z (m)', 'u (m/s)', 'V (1/s)', + 'T (K)', 'rho (kg/m3)'] + self.gas.species_names) + for n in range(self.flame.n_points): + self.set_gas_state(n) + writer.writerow([z[n], u[n], V[n], T[n], self.gas.density] + + list(getattr(self.gas, species))) + csvfile.close() + if not quiet: + print("Solution saved to '{0}'.".format(filename)) + + +def _trim(docstring): + """Remove block indentation from a docstring.""" + if not docstring: + return '' + lines = docstring.splitlines() + # Determine minimum indentation (first line doesn't count): + indent = 999 + for line in lines[1:]: + stripped = line.lstrip() + if stripped: + indent = min(indent, len(line) - len(stripped)) + # Remove indentation (first line is special): + trimmed = [lines[0].strip()] + if indent < 999: + for line in lines[1:]: + trimmed.append(line[indent:].rstrip()) + + # Return a single string, with trailing and leading blank lines stripped + return '\n'.join(trimmed).strip('\n') + + +def _array_property(attr, size=None): + """ + Generate a property that retrieves values at each point in the flame. The + 'size' argument is the attribute name of the gas object used to set the + leading dimension of the resulting array. + """ + def getter(self): + if size is None: + # 1D array for scalar property + vals = np.empty(self.flame.n_points) + else: + # 2D array + vals = np.empty((getattr(self.gas, size), self.flame.n_points)) + + for i in range(self.flame.n_points): + self.set_gas_state(i) + vals[...,i] = getattr(self.gas, attr) + + return vals + + if size is None: + extradoc = "\nReturns an array of length `n_points`." + else: + extradoc = "\nReturns an array of size `%s` x `n_points`." % size + + doc = _trim(getattr(Solution, attr).__doc__) + extradoc + return property(getter, doc=doc) + +# Add scalar properties to FlameBase +for attr in ['density', 'density_mass', 'density_mole', 'volume_mass', + 'volume_mole', 'int_energy_mole', 'int_energy_mass', 'h', + 'enthalpy_mole', 'enthalpy_mass', 's', 'entropy_mole', + 'entropy_mass', 'g', 'gibbs_mole', 'gibbs_mass', 'cv', + 'cv_mole', 'cv_mass', 'cp', 'cp_mole', 'cp_mass', + 'isothermal_compressibility', 'thermal_expansion_coeff', + 'viscosity', 'thermal_conductivity']: + setattr(FlameBase, attr, _array_property(attr)) +FlameBase.volume = _array_property('v') # avoid confusion with velocity gradient 'V' +FlameBase.int_energy = _array_property('u') # avoid collision with velocity 'u' + +# Add properties with values for each species +for attr in ['X', 'Y', 'concentrations', 'partial_molar_enthalpies', + 'partial_molar_entropies', 'partial_molar_int_energies', + 'chemical_potentials', 'electrochemical_potentials', 'partial_molar_cp', + 'partial_molar_volumes', 'standard_enthalpies_RT', + 'standard_entropies_R', 'standard_int_energies_RT', + 'standard_gibbs_RT', 'standard_cp_R', 'creation_rates', + 'destruction_rates', 'net_production_rates', 'mix_diff_coeffs', + 'mix_diff_coeffs_mass', 'mix_diff_coeffs_mole', 'thermal_diff_coeffs']: + setattr(FlameBase, attr, _array_property(attr, 'n_species')) + +# Add properties with values for each reaction +for attr in ['forward_rates_of_progress', 'reverse_rates_of_progress', 'net_rates_of_progress', + 'equilibrium_constants', 'forward_rate_constants', 'reverse_rate_constants', + 'delta_enthalpy', 'delta_gibbs', 'delta_entropy', + 'delta_standard_enthalpy', 'delta_standard_gibbs', + 'delta_standard_entropy']: + setattr(FlameBase, attr, _array_property(attr, 'n_reactions')) + + +class FreeFlame(FlameBase): + """A freely-propagating flat flame.""" + + def __init__(self, gas, grid=None): + """ + A domain of type FreeFlow named 'flame' will be created to represent + the flame. The three domains comprising the stack are stored as + ``self.inlet``, ``self.flame``, and ``self.outlet``. + """ + self.inlet = Inlet1D(name='reactants', phase=gas) + self.outlet = Outlet1D(name='products', phase=gas) + self.flame = FreeFlow(gas, name='flame') + + super(FreeFlame, self).__init__((self.inlet, self.flame, self.outlet), + gas, grid) + + def set_initial_guess(self): + """ + Set the initial guess for the solution. The adiabatic flame + temperature and equilibrium composition are computed for the inlet gas + composition. The temperature profile rises linearly over 20% of the + domain width to Tad, then is flat. The mass fraction profiles are set + similarly. + """ + super(FreeFlame, self).set_initial_guess() + self.gas.TPY = self.inlet.T, self.P, self.inlet.Y + Y0 = self.inlet.Y + u0 = self.inlet.mdot/self.gas.density + T0 = self.inlet.T + + # get adiabatic flame temperature and composition + self.gas.equilibrate('HP') + Teq = self.gas.T + Yeq = self.gas.Y + u1 = self.inlet.mdot/self.gas.density + + locs = [0.0, 0.3, 0.5, 1.0] + self.set_profile('u', locs, [u0, u0, u1, u1]) + self.set_profile('T', locs, [T0, T0, Teq, Teq]) + self.set_fixed_temperature(0.5 * (T0 + Teq)) + for n in range(self.gas.n_species): + self.set_profile(self.gas.species_name(n), + locs, [Y0[n], Y0[n], Yeq[n], Yeq[n]]) + + +class BurnerFlame(FlameBase): + """A burner-stabilized flat flame.""" + + def __init__(self, gas, grid=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 + + A domain of class `AxisymmetricStagnationFlow` named ``flame`` will + be created to represent the flame. The three domains comprising the + stack are stored as ``self.burner``, ``self.flame``, and + ``self.outlet``. + """ + self.burner = Inlet1D(name='burner', phase=gas) + self.burner.T = gas.T + self.outlet = Outlet1D(name='outlet', phase=gas) + self.flame = AxisymmetricStagnationFlow(gas, name='flame') + + super(BurnerFlame, self).__init__((self.burner, self.flame, self.outlet), + gas, grid) + + def set_initial_guess(self): + """ + Set the initial guess for the solution. The adiabatic flame + temperature and equilibrium composition are computed for the burner + gas composition. The temperature profile rises linearly in the first + 20% of the flame to Tad, then is flat. The mass fraction profiles are + set similarly. + """ + super(BurnerFlame, self).set_initial_guess() + + self.gas.TPY = self.burner.T, self.P, self.burner.Y + Y0 = self.burner.Y + u0 = self.burner.mdot/self.gas.density + T0 = self.burner.T + + # get adiabatic flame temperature and composition + self.gas.equilibrate('HP') + Teq = self.gas.T + Yeq = self.gas.Y + u1 = self.burner.mdot/self.gas.density + + locs = [0.0, 0.2, 1.0] + self.set_profile('u', locs, [u0, u1, u1]) + self.set_profile('T', locs, [T0, Teq, Teq]) + for n in range(self.gas.n_species): + self.set_profile(self.gas.species_name(n), + locs, [Y0[n], Yeq[n], Yeq[n]]) + + +class CounterflowDiffusionFlame(FlameBase): + """ A counterflow diffusion flame """ + + def __init__(self, gas, grid=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 + + A domain of class `AxisymmetricStagnationFlow` named ``flame`` will + be created to represent the flame. The three domains comprising the + stack are stored as ``self.fuel_inlet``, ``self.flame``, and + ``self.oxidizer_inlet``. + """ + self.fuel_inlet = Inlet1D(name='fuel_inlet', phase=gas) + self.fuel_inlet.T = gas.T + + self.oxidizer_inlet = Inlet1D(name='oxidizer_inlet', phase=gas) + self.oxidizer_inlet.T = gas.T + + self.flame = AxisymmetricStagnationFlow(gas, name='flame') + + super(CounterflowDiffusionFlame, self).__init__( + (self.fuel_inlet, self.flame, self.oxidizer_inlet), gas, grid) + + def set_initial_guess(self, fuel, oxidizer='O2', stoich=None): + """ + Set the initial guess for the solution. The fuel species must be + specified: + + >>> f.set_initial_guess(fuel='CH4') + + The oxidizer and corresponding stoichiometry must be specified if it + is not 'O2'. The initial guess is generated by assuming infinitely- + fast chemistry. + """ + + super(CounterflowDiffusionFlame, self).set_initial_guess() + + if stoich is None: + if oxidizer == 'O2': + stoich = 0.0 + if 'H' in self.gas.element_names: + stoich += 0.25 * self.gas.n_atoms(fuel, 'H') + if 'C' in self.gas.element_names: + stoich += self.gas.n_atoms(fuel, 'C') + else: + raise Exception('oxidizer/fuel stoichiometric ratio must be ' + 'specified since the oxidizer is not O2') + + kFuel = self.gas.species_index(fuel) + kOx = self.gas.species_index(oxidizer) + + s = stoich * self.gas.molecular_weights[kOx] / self.gas.molecular_weights[kFuel] + phi = s * self.fuel_inlet.Y[kFuel] / self.oxidizer_inlet.Y[kOx] + zst = 1.0 / (1.0 + phi) + + Yin_f = self.fuel_inlet.Y + Yin_o = self.oxidizer_inlet.Y + Yst = zst * Yin_f + (1.0 - zst) * Yin_o + + self.gas.TPY = self.fuel_inlet.T, self.P, Yin_f + mdotf = self.fuel_inlet.mdot + u0f = mdotf / self.gas.density + T0f = self.fuel_inlet.T + + self.gas.TPY = self.oxidizer_inlet.T, self.P, Yin_o + mdoto = self.oxidizer_inlet.mdot + u0o = mdoto/self.gas.density + T0o = self.oxidizer_inlet.T + + # get adiabatic flame temperature and composition + Tbar = 0.5 * (T0f + T0o) + self.gas.TPY = Tbar, self.P, Yst + self.gas.equilibrate('HP') + Teq = self.gas.T + Yeq = self.gas.Y + + # estimate strain rate + zz = self.flame.grid + dz = zz[-1] - zz[0] + a = (u0o + u0f)/dz + f = np.sqrt(a / (2.0 * self.gas.mix_diff_coeffs[kOx])) + + x0 = mdotf * dz / (mdotf + mdoto) + nz = len(zz) + + Y = np.zeros((nz, self.gas.n_species)) + T = np.zeros(nz) + for j in range(nz): + x = zz[j] + zeta = f * (x - x0) + zmix = 0.5 * (1.0 - erf(zeta)) + if zmix > zst: + Y[j] = Yeq + (Yin_f - Yeq) * (zmix - zst) / (1.0 - zst) + T[j] = Teq + (T0f - Teq) * (zmix - zst) / (1.0 - zst) + else: + Y[j] = Yin_o + zmix * (Yeq - Yin_o) / zst + T[j] = T0o + (Teq - T0o) * zmix / zst + + T[0] = T0f + T[-1] = T0o + zrel = zz/dz + + self.set_profile('u', [0.0, 1.0], [u0f, -u0o]) + self.set_profile('V', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) + self.set_profile('T', zrel, T) + for k,spec in enumerate(self.gas.species_names): + self.set_profile(spec, zrel, Y[:,k]) + + +class ImpingingJet(FlameBase): + """An axisymmetric flow impinging on a surface at normal incidence.""" + def __init__(self, gas, grid=None, surface=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 + :param surface: + A Kinetics object used to compute any surface reactions. + + A domain of class `AxisymmetricStagnationFlow` named ``flame`` will be + created to represent the flow. The three domains comprising the stack + are stored as ``self.inlet``, ``self.flame``, and ``self.surface``. + """ + self.inlet = Inlet1D(name='inlet', phase=gas) + self.inlet.T = gas.T + self.flame = AxisymmetricStagnationFlow(gas, name='flame') + + if surface is None: + self.surface = Surface1D(name='surface', phase=gas) + self.surface.T = gas.T + else: + self.surface = ReactingSurface1D(name='surface', phase=gas) + self.surface.set_kinetics(surface) + self.surface.T = surface.T + + super(ImpingingJet, self).__init__( + (self.inlet, self.flame, self.surface), gas, grid) + + def set_initial_guess(self, products='inlet'): + """ + Set the initial guess for the solution. If products = 'equil', then + the equilibrium composition at the adiabatic flame temperature will be + used to form the initial guess. Otherwise the inlet composition will + be used. + """ + super(ImpingingJet, self).set_initial_guess() + + Y0 = self.inlet.Y + T0 = self.inlet.T + self.gas.TPY = T0, self.flame.P, Y0 + u0 = self.inlet.mdot / self.gas.density + + if products == 'equil': + self.gas.equilibrate('HP') + Teq = self.gas.T + Yeq = self.gas.Y + locs = np.array([0.0, 0.3, 0.7, 1.0]) + self.set_profile('T', locs, [T0, Teq, Teq, self.surface.T]) + for k in range(self.gas.n_species): + self.set_profile(self.gas.species_name(k), locs, + [Y0[k], Yeq[k], Yeq[k], Yeq[k]]) + else: + locs = np.array([0.0, 1.0]) + self.set_profile('T', locs, [T0, self.surface.T]) + for k in range(self.gas.n_species): + self.set_profile(self.gas.species_name(k), locs, + [Y0[k], Y0[k]]) + + locs = np.array([0.0, 1.0]) + self.set_profile('u', locs, [u0, 0.0]) + self.set_profile('V', locs, [0.0, 0.0]) diff --git a/interfaces/cython/cantera/onedim.pyx b/interfaces/cython/cantera/onedim.pyx index 0358de3b2..f35ed2f52 100644 --- a/interfaces/cython/cantera/onedim.pyx +++ b/interfaces/cython/cantera/onedim.pyx @@ -1,12 +1,6 @@ import csv import interrupts -try: - # Python 2.7 or 3.2+ - from math import erf -except ImportError: - from scipy.special import erf - cdef class Domain1D: def __cinit__(self, *args, **kwargs): self.domain = NULL @@ -791,488 +785,3 @@ cdef class Sim1D: def __dealloc__(self): del self.sim - - -class FlameBase(Sim1D): - """ Base class for flames with a single flow domain """ - - def __init__(self, domains, gas, grid=None): - """ - :param gas: - object to use to evaluate all gas properties and reaction rates - :param grid: - array of initial grid points - """ - if grid is None: - grid = np.linspace(0.0, 0.1, 6) - self.flame.grid = grid - super().__init__(domains) - self.gas = gas - self.flame.P = gas.P - - def set_refine_criteria(self, ratio=10.0, slope=0.8, curve=0.8, prune=0.0): - super().set_refine_criteria(self.flame, ratio, slope, curve, prune) - - def set_profile(self, component, locations, values): - super().set_profile(self.flame, component, locations, values) - - @property - def transport_model(self): - return self.gas.transport_model - - @transport_model.setter - def transport_model(self, model): - self.gas.transport_model = model - self.flame.set_transport(self.gas) - - @property - def energy_enabled(self): - return self.flame.energy_enabled - - @energy_enabled.setter - def energy_enabled(self, enable): - self.flame.energy_enabled = enable - - @property - def soret_enabled(self): - return self.flame.soret_enabled - - @soret_enabled.setter - def soret_enabled(self, enable): - self.flame.soret_enabled = enable - - @property - def grid(self): - """ Array of grid point positions along the flame. """ - return self.flame.grid - - @property - def P(self): - return self.flame.P - - @P.setter - def P(self, P): - self.flame.P = P - - @property - def T(self): - """ Array containing the temperature [K] at each grid point. """ - return self.profile(self.flame, 'T') - - @property - def u(self): - """ - Array containing the velocity [m/s] normal to the flame at each point. - """ - return self.profile(self.flame, 'u') - - @property - def V(self): - """ - Array containing the tangential velocity gradient [1/s] at each point. - """ - return self.profile(self.flame, 'V') - - @property - def L(self): - """ - Array containing the radial pressure gradient (1/r)(dP/dr) [N/m^4] at - each point. Note: This value is named 'lambda' in the C++ code. - """ - return self.profile(self.flame, 'lambda') - - def solution(self, component, point=None): - if point is None: - return self.profile(self.flame, component) - else: - return self.value(self.flame, component, point) - - def set_gas_state(self, point): - k0 = self.flame.component_index(self.gas.species_name(0)) - Y = [self.solution(k, point) - for k in range(k0, k0 + self.gas.n_species)] - self.gas.TPY = self.value(self.flame, 'T', point), self.P, Y - - def write_csv(self, filename, species='X', quiet=True): - """ - Write the velocity, temperature, density, and species profiles - to a CSV file. - - :param filename: - Output file name - :param species: - Attribute to use obtaining species profiles, e.g. ``X`` for - mole fractions or ``Y`` for mass fractions. - """ - - z = self.grid - T = self.T - u = self.u - V = self.V - - csvfile = open(filename, 'w') - writer = csv.writer(csvfile) - writer.writerow(['z (m)', 'u (m/s)', 'V (1/s)', - 'T (K)', 'rho (kg/m3)'] + self.gas.species_names) - for n in range(self.flame.n_points): - self.set_gas_state(n) - writer.writerow([z[n], u[n], V[n], T[n], self.gas.density] + - list(getattr(self.gas, species))) - csvfile.close() - if not quiet: - print("Solution saved to '{0}'.".format(filename)) - - -def _trim(docstring): - """Remove block indentation from a docstring.""" - if not docstring: - return '' - lines = docstring.splitlines() - # Determine minimum indentation (first line doesn't count): - indent = 999 - for line in lines[1:]: - stripped = line.lstrip() - if stripped: - indent = min(indent, len(line) - len(stripped)) - # Remove indentation (first line is special): - trimmed = [lines[0].strip()] - if indent < 999: - for line in lines[1:]: - trimmed.append(line[indent:].rstrip()) - - # Return a single string, with trailing and leading blank lines stripped - return '\n'.join(trimmed).strip('\n') - -def _array_property(attr, size=None): - """ - Generate a property that retrieves values at each point in the flame. The - 'size' argument is the attribute name of the gas object used to set the - leading dimension of the resulting array. - """ - def getter(self): - if size is None: - # 1D array for scalar property - vals = np.empty(self.flame.n_points) - else: - # 2D array - vals = np.empty((getattr(self.gas, size), self.flame.n_points)) - - for i in range(self.flame.n_points): - self.set_gas_state(i) - vals[...,i] = getattr(self.gas, attr) - - return vals - - if size is None: - extradoc = "\nReturns an array of length `n_points`." - else: - extradoc = "\nReturns an array of size `%s` x `n_points`." % size - - doc = _trim(getattr(Solution, attr).__doc__) + extradoc - return property(getter, doc=doc) - -# Add scalar properties to FlameBase -for attr in ['density', 'density_mass', 'density_mole', 'volume_mass', - 'volume_mole', 'int_energy_mole', 'int_energy_mass', 'h', - 'enthalpy_mole', 'enthalpy_mass', 's', 'entropy_mole', - 'entropy_mass', 'g', 'gibbs_mole', 'gibbs_mass', 'cv', - 'cv_mole', 'cv_mass', 'cp', 'cp_mole', 'cp_mass', - 'isothermal_compressibility', 'thermal_expansion_coeff', - 'viscosity', 'thermal_conductivity']: - setattr(FlameBase, attr, _array_property(attr)) -FlameBase.volume = _array_property('v') # avoid confusion with velocity gradient 'V' -FlameBase.int_energy = _array_property('u') # avoid collision with velocity 'u' - -# Add properties with values for each species -for attr in ['X', 'Y', 'concentrations', 'partial_molar_enthalpies', - 'partial_molar_entropies', 'partial_molar_int_energies', - 'chemical_potentials', 'electrochemical_potentials', 'partial_molar_cp', - 'partial_molar_volumes', 'standard_enthalpies_RT', - 'standard_entropies_R', 'standard_int_energies_RT', - 'standard_gibbs_RT', 'standard_cp_R', 'creation_rates', - 'destruction_rates', 'net_production_rates', 'mix_diff_coeffs', - 'mix_diff_coeffs_mass', 'mix_diff_coeffs_mole', 'thermal_diff_coeffs']: - setattr(FlameBase, attr, _array_property(attr, 'n_species')) - -# Add properties with values for each reaction -for attr in ['forward_rates_of_progress', 'reverse_rates_of_progress', 'net_rates_of_progress', - 'equilibrium_constants', 'forward_rate_constants', 'reverse_rate_constants', - 'delta_enthalpy', 'delta_gibbs', 'delta_entropy', - 'delta_standard_enthalpy', 'delta_standard_gibbs', - 'delta_standard_entropy']: - setattr(FlameBase, attr, _array_property(attr, 'n_reactions')) - - -class FreeFlame(FlameBase): - """A freely-propagating flat flame.""" - - def __init__(self, gas, grid=None): - """ - A domain of type FreeFlow named 'flame' will be created to represent - the flame. The three domains comprising the stack are stored as - ``self.inlet``, ``self.flame``, and ``self.outlet``. - """ - self.inlet = Inlet1D(name='reactants', phase=gas) - self.outlet = Outlet1D(name='products', phase=gas) - self.flame = FreeFlow(gas, name='flame') - - super().__init__((self.inlet, self.flame, self.outlet), gas, grid) - - def set_initial_guess(self): - """ - Set the initial guess for the solution. The adiabatic flame - temperature and equilibrium composition are computed for the inlet gas - composition. The temperature profile rises linearly over 20% of the - domain width to Tad, then is flat. The mass fraction profiles are set - similarly. - """ - super().set_initial_guess() - self.gas.TPY = self.inlet.T, self.P, self.inlet.Y - Y0 = self.inlet.Y - u0 = self.inlet.mdot/self.gas.density - T0 = self.inlet.T - - # get adiabatic flame temperature and composition - self.gas.equilibrate('HP') - Teq = self.gas.T - Yeq = self.gas.Y - u1 = self.inlet.mdot/self.gas.density - - locs = [0.0, 0.3, 0.5, 1.0] - self.set_profile('u', locs, [u0, u0, u1, u1]) - self.set_profile('T', locs, [T0, T0, Teq, Teq]) - self.set_fixed_temperature(0.5 * (T0 + Teq)) - for n in range(self.gas.n_species): - self.set_profile(self.gas.species_name(n), - locs, [Y0[n], Y0[n], Yeq[n], Yeq[n]]) - - -class BurnerFlame(FlameBase): - """A burner-stabilized flat flame.""" - - def __init__(self, gas, grid=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 - - A domain of class `AxisymmetricStagnationFlow` named ``flame`` will - be created to represent the flame. The three domains comprising the - stack are stored as ``self.burner``, ``self.flame``, and - ``self.outlet``. - """ - self.burner = Inlet1D(name='burner', phase=gas) - self.burner.T = gas.T - self.outlet = Outlet1D(name='outlet', phase=gas) - self.flame = AxisymmetricStagnationFlow(gas, name='flame') - - super().__init__((self.burner, self.flame, self.outlet), gas, grid) - - def set_initial_guess(self): - """ - Set the initial guess for the solution. The adiabatic flame - temperature and equilibrium composition are computed for the burner - gas composition. The temperature profile rises linearly in the first - 20% of the flame to Tad, then is flat. The mass fraction profiles are - set similarly. - """ - super().set_initial_guess() - - self.gas.TPY = self.burner.T, self.P, self.burner.Y - Y0 = self.burner.Y - u0 = self.burner.mdot/self.gas.density - T0 = self.burner.T - - # get adiabatic flame temperature and composition - self.gas.equilibrate('HP') - Teq = self.gas.T - Yeq = self.gas.Y - u1 = self.burner.mdot/self.gas.density - - locs = [0.0, 0.2, 1.0] - self.set_profile('u', locs, [u0, u1, u1]) - self.set_profile('T', locs, [T0, Teq, Teq]) - for n in range(self.gas.n_species): - self.set_profile(self.gas.species_name(n), - locs, [Y0[n], Yeq[n], Yeq[n]]) - - -class CounterflowDiffusionFlame(FlameBase): - """ A counterflow diffusion flame """ - - def __init__(self, gas, grid=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 - - A domain of class `AxisymmetricStagnationFlow` named ``flame`` will - be created to represent the flame. The three domains comprising the - stack are stored as ``self.fuel_inlet``, ``self.flame``, and - ``self.oxidizer_inlet``. - """ - self.fuel_inlet = Inlet1D(name='fuel_inlet', phase=gas) - self.fuel_inlet.T = gas.T - - self.oxidizer_inlet = Inlet1D(name='oxidizer_inlet', phase=gas) - self.oxidizer_inlet.T = gas.T - - self.flame = AxisymmetricStagnationFlow(gas, name='flame') - - super().__init__((self.fuel_inlet, self.flame, self.oxidizer_inlet), - gas, grid) - - def set_initial_guess(self, fuel, oxidizer='O2', stoich=None): - """ - Set the initial guess for the solution. The fuel species must be - specified: - - >>> f.set_initial_guess(fuel='CH4') - - The oxidizer and corresponding stoichiometry must be specified if it - is not 'O2'. The initial guess is generated by assuming infinitely- - fast chemistry. - """ - - super().set_initial_guess() - - if stoich is None: - if oxidizer == 'O2': - stoich = 0.0 - if 'H' in self.gas.element_names: - stoich += 0.25 * self.gas.n_atoms(fuel, 'H') - if 'C' in self.gas.element_names: - stoich += self.gas.n_atoms(fuel, 'C') - else: - raise Exception('oxidizer/fuel stoichiometric ratio must be ' - 'specified since the oxidizer is not O2') - - kFuel = self.gas.species_index(fuel) - kOx = self.gas.species_index(oxidizer) - - s = stoich * self.gas.molecular_weights[kOx] / self.gas.molecular_weights[kFuel] - phi = s * self.fuel_inlet.Y[kFuel] / self.oxidizer_inlet.Y[kOx] - zst = 1.0 / (1.0 + phi) - - Yin_f = self.fuel_inlet.Y - Yin_o = self.oxidizer_inlet.Y - Yst = zst * Yin_f + (1.0 - zst) * Yin_o - - self.gas.TPY = self.fuel_inlet.T, self.P, Yin_f - mdotf = self.fuel_inlet.mdot - u0f = mdotf / self.gas.density - T0f = self.fuel_inlet.T - - self.gas.TPY = self.oxidizer_inlet.T, self.P, Yin_o - mdoto = self.oxidizer_inlet.mdot - u0o = mdoto/self.gas.density - T0o = self.oxidizer_inlet.T - - # get adiabatic flame temperature and composition - Tbar = 0.5 * (T0f + T0o) - self.gas.TPY = Tbar, self.P, Yst - self.gas.equilibrate('HP') - Teq = self.gas.T - Yeq = self.gas.Y - - # estimate strain rate - zz = self.flame.grid - dz = zz[-1] - zz[0] - a = (u0o + u0f)/dz - f = np.sqrt(a / (2.0 * self.gas.mix_diff_coeffs[kOx])) - - x0 = mdotf * dz / (mdotf + mdoto) - nz = len(zz) - - Y = np.zeros((nz, self.gas.n_species)) - T = np.zeros(nz) - for j in range(nz): - x = zz[j] - zeta = f * (x - x0) - zmix = 0.5 * (1.0 - erf(zeta)) - if zmix > zst: - Y[j] = Yeq + (Yin_f - Yeq) * (zmix - zst) / (1.0 - zst) - T[j] = Teq + (T0f - Teq) * (zmix - zst) / (1.0 - zst) - else: - Y[j] = Yin_o + zmix * (Yeq - Yin_o) / zst - T[j] = T0o + (Teq - T0o) * zmix / zst - - T[0] = T0f - T[-1] = T0o - zrel = zz/dz - - self.set_profile('u', [0.0, 1.0], [u0f, -u0o]) - self.set_profile('V', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) - self.set_profile('T', zrel, T) - for k,spec in enumerate(self.gas.species_names): - self.set_profile(spec, zrel, Y[:,k]) - - -class ImpingingJet(FlameBase): - """An axisymmetric flow impinging on a surface at normal incidence.""" - def __init__(self, gas, grid=None, surface=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 - :param surface: - A Kinetics object used to compute any surface reactions. - - A domain of class `AxisymmetricStagnationFlow` named ``flame`` will be - created to represent the flow. The three domains comprising the stack - are stored as ``self.inlet``, ``self.flame``, and ``self.surface``. - """ - self.inlet = Inlet1D(name='inlet', phase=gas) - self.inlet.T = gas.T - self.flame = AxisymmetricStagnationFlow(gas, name='flame') - - if surface is None: - self.surface = Surface1D(name='surface', phase=gas) - self.surface.T = gas.T - else: - self.surface = ReactingSurface1D(name='surface', phase=gas) - self.surface.set_kinetics(surface) - self.surface.T = surface.T - - super().__init__((self.inlet, self.flame, self.surface), - gas, grid) - - def set_initial_guess(self, products='inlet'): - """ - Set the initial guess for the solution. If products = 'equil', then - the equilibrium composition at the adiabatic flame temperature will be - used to form the initial guess. Otherwise the inlet composition will - be used. - """ - super().set_initial_guess() - - Y0 = self.inlet.Y - T0 = self.inlet.T - self.gas.TPY = T0, self.flame.P, Y0 - u0 = self.inlet.mdot / self.gas.density - - if products == 'equil': - self.gas.equilibrate('HP') - Teq = self.gas.T - Yeq = self.gas.Y - locs = np.array([0.0, 0.3, 0.7, 1.0]) - self.set_profile('T', locs, [T0, Teq, Teq, self.surface.T]) - for k in range(self.gas.n_species): - self.set_profile(self.gas.species_name(k), locs, - [Y0[k], Yeq[k], Yeq[k], Yeq[k]]) - else: - locs = np.array([0.0, 1.0]) - self.set_profile('T', locs, [T0, self.surface.T]) - for k in range(self.gas.n_species): - self.set_profile(self.gas.species_name(k), locs, - [Y0[k], Y0[k]]) - - locs = np.array([0.0, 1.0]) - self.set_profile('u', locs, [u0, 0.0]) - self.set_profile('V', locs, [0.0, 0.0])