[Matlab] Make counterflow diffusion flame simulation more general
The CounterFlowDiffusionFlame (CFDF) code is able to perform more general cases of npflame_init for multiple species fuel and oxidizer streams. The stoichiometric mixture fraction in the CFDF code uses the Bilger definition of mixture fraction, using the conservation of elements C, H, and O. This method is used in the python module, but not the MATLAB npflame_init function. Also, the CFDF code uses the fuel stream density to calculate the fuel stream velocity and the oxidizer stream density to calculate the oxidizer stream velocity, where as the npflame_init code uses the fuel density for both velocity calculations. The elementMassFraction code is a MATLAB version of the python function: elemental_mass_fraction, which is needed to run the CFDF code. Update the diffflame.m example to use the more general CFDF function since the input parameters are different than the npflame_init function. This example is the same as the diffusion_flame.py sample in the Python module.
This commit is contained in:
KyleLinevitchJr
Ray Speth
parent
94a4439bd4
commit
43bd96b5ba
4 changed files with 372 additions and 90 deletions
|
|
@ -207,7 +207,8 @@ if localenv['sphinx_docs']:
|
|||
'1D/@Domain1D': ['1D/AxiStagnFlow.m', '1D/AxisymmetricFlow.m',
|
||||
'1D/Inlet.m', '1D/Outlet.m', '1D/OutletRes.m',
|
||||
'1D/Surface.m', '1D/SymmPlane.m'],
|
||||
'1D/@Stack': ['1D/FreeFlame.m', '1D/npflame_init.m'],
|
||||
'1D/@Stack': ['1D/FreeFlame.m', '1D/npflame_init.m',
|
||||
'1D/CounterFlowDiffusionFlame.m'],
|
||||
'@Interface': ['importEdge.m', 'importInterface.m'],
|
||||
'@Data': ['Air.m', 'gasconstant.m', 'GRI30.m',
|
||||
'Hydrogen.m', 'Methane.m', 'Nitrogen.m', 'oneatm.m',
|
||||
|
|
|
|||
220
interfaces/matlab/toolbox/1D/CounterFlowDiffusionFlame.m
Normal file
220
interfaces/matlab/toolbox/1D/CounterFlowDiffusionFlame.m
Normal file
|
|
@ -0,0 +1,220 @@
|
|||
function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidizer)
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% COUNTERFLOWDIFFUSIONFLAME Create a counter flow diffusion flame stack.
|
||||
% flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidizer)
|
||||
% :param left:
|
||||
% Object representing the left inlet, which must be
|
||||
% created using function :mat:func:`Inlet`.
|
||||
% :param flow:
|
||||
% Object representing the flow, created with
|
||||
% function :mat:func:`AxisymmetricFlow`.
|
||||
% :param right:
|
||||
% Object representing the right inlet, which must be
|
||||
% created using function :mat:func:`Inlet`.
|
||||
% :param tp_f:
|
||||
% Object representing the fuel inlet gas, instance of class
|
||||
% :mat:func:`Solution`, and an ideal gas.
|
||||
% :param tp_o:
|
||||
% Object representing the oxidizer inlet gas, instance of class
|
||||
% :mat:func:`Solution`, and an ideal gas.
|
||||
% :param oxidizer:
|
||||
% String representing the oxidizer species. Most commonly O2.
|
||||
% :return:
|
||||
% Instance of :mat:func:`Stack` object representing the left
|
||||
% inlet, flow, and right inlet.
|
||||
%
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Check input parameters
|
||||
%
|
||||
|
||||
if nargin ~= 6
|
||||
error('CounterFlowDiffusionFlame expects six input arguments.');
|
||||
end
|
||||
if ~isIdealGas(tp_f)
|
||||
error('Fuel gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
if ~isIdealGas(tp_o)
|
||||
error('Oxidizer gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
if ~isInlet(left)
|
||||
error('Left inlet object of wrong type.');
|
||||
end
|
||||
if ~isFlow(flow)
|
||||
error('Flow object of wrong type.');
|
||||
end
|
||||
if ~isInlet(right)
|
||||
error('Right inlet object of wrong type.');
|
||||
end
|
||||
if ~ischar(oxidizer)
|
||||
error('Oxidizer name must be of format character.');
|
||||
end
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Get the density of both fuel and oxidizer streams. To be used in
|
||||
% determining velocity of each stream. Also get the temperature of both
|
||||
% inlet streams.
|
||||
%
|
||||
|
||||
rhof = density(tp_f);
|
||||
rho0 = density(tp_o);
|
||||
tf = temperature(left);
|
||||
tox = temperature(right);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Find the species index of the oxidizer. To be used in determining initial
|
||||
% strain rate.
|
||||
%
|
||||
|
||||
ioxidizer = speciesIndex(tp_o, oxidizer);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Calculate the stoichiometric mixture fraction. Needed for determining
|
||||
% location of flame edges and composition. elMoles function used to
|
||||
% calculate the number of moles of C, H, and O atoms in the fuel and
|
||||
% oxidizer streams: elMoles = elementalMassFraction/element atomic weight.
|
||||
% From this, the stoichiometric Air/Fuel ratio can be determined.
|
||||
% 1 Mole of O needs 2 Moles of C and 0.5 Moles of H for stoichiometric
|
||||
% conditions. The stoichiometric mixture fraction, Zst, is then calculated.
|
||||
%
|
||||
|
||||
sFuel = elMoles(tp_f,'O')- 2*elMoles(tp_f,'C')- 0.5*elMoles(tp_f,'H');
|
||||
sOx = elMoles(tp_o,'O')- 2*elMoles(tp_o,'C')- 0.5*elMoles(tp_o,'H');
|
||||
phi = sFuel/sOx;
|
||||
zst = 1.0/(1.0 - phi);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Compute the stoichiometric mass fractions of each species. Use this to
|
||||
% set the fuel gas object and calculate adiabatic flame temperature and
|
||||
% equilibrium composition.
|
||||
%
|
||||
|
||||
spec = speciesNames(tp_f); % Get all of the species names in gas object.
|
||||
nsp = nSpecies(tp_f); % Get total number of species in gas object.
|
||||
% Get the current mass fractions of both fuel and inlet streams.
|
||||
yox = massFractions(tp_o);
|
||||
yf = massFractions(tp_f);
|
||||
ystoich_double = zeros(1, nsp); % Create empty vector for stoich mass frac.
|
||||
|
||||
for n = 1:nsp
|
||||
% Calculate stoichiometric mass fractions.
|
||||
ystoich_double(n) = zst*yf(n) + (1.0 - zst)*yox(n);
|
||||
% Convert mass fraction vector to string vector.
|
||||
ystoich_str{n} = num2str(ystoich_double(n));
|
||||
% Convert string vector to cell with SPECIES:MASS FRACTION format.
|
||||
y_stoich{n} = [spec{n}, ':', ystoich_str{n}];
|
||||
end
|
||||
% Initialize stoichiometric mass fraction cell with first SP:Y value.
|
||||
ystoich = [y_stoich{1}];
|
||||
for i = 2:nsp
|
||||
% Update cell to have format similar to N2:Yst,O2:Yst,...
|
||||
ystoich = [ystoich ',', y_stoich{i}];
|
||||
end
|
||||
% Set the fuel gas object as stoichiometric values and use equilibrate
|
||||
% function to determine stoichiometric equilibrium temperature and mass
|
||||
% fractions.
|
||||
set(tp_f, 'T', tf, 'P', pressure(tp_f), 'Y', ystoich);
|
||||
equilibrate(tp_f, 'HP');
|
||||
teq = temperature(tp_f);
|
||||
yeq = massFractions(tp_f);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Estimate the strain rate based on the inlet stream velocities and
|
||||
% determine initial "guess" for mixture fraction based on mass flux ratio.
|
||||
%
|
||||
|
||||
zz = gridPoints(flow);
|
||||
dz = zz(end) - zz(1);
|
||||
uleft = massFlux(left)/rhof;
|
||||
uright = massFlux(right)/rho0;
|
||||
a = (abs(uleft) + abs(uright))/dz;
|
||||
diff = mixDiffCoeffs(tp_f);
|
||||
f = sqrt(a/(2.0*diff(ioxidizer)));
|
||||
x0num = sqrt(uleft*massFlux(left))*dz;
|
||||
x0den = sqrt(uleft*massFlux(left)) + sqrt(uright*massFlux(right));
|
||||
x0 = x0num/x0den;
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Calculate initial values of temperature and mass fraction of species in
|
||||
% flame at each gridpoint. These values to be used for energy equation
|
||||
% solution. Method is based on the Burke-Schumann model.
|
||||
%
|
||||
|
||||
nz = nPoints(flow);
|
||||
zm = zeros(1, nz);
|
||||
u = zeros(1, nz);
|
||||
v = zeros(1, nz);
|
||||
y = zeros(nz, nsp);
|
||||
t = zeros(1, nz);
|
||||
for j = 1:nz
|
||||
x = zz(j);
|
||||
zeta = f*(x - x0);
|
||||
zmix = 0.5*(1.0 - erf(zeta)); % Mixture fraction in flame.
|
||||
zm(j) = zmix;
|
||||
u(j) = a*(x0 - zz(j)); % Axial velocity.
|
||||
v(j) = a; % Radial velocity.
|
||||
if zmix > zst
|
||||
for n = 1:nsp
|
||||
y(j,n) = yeq(n) + (zmix - zst)*(yf(n) - yeq(n))/(1.0 - zst);
|
||||
end
|
||||
t(j) = teq + (tf - teq)*(zmix - zst)/(1.0 - zst);
|
||||
else
|
||||
for n = 1:nsp
|
||||
y(j,n) = yox(n) + zmix*(yeq(n) - yox(n))/zst;
|
||||
end
|
||||
t(j) = tox + zmix*(teq - tox)/zst;
|
||||
end
|
||||
end
|
||||
zrel = zz/dz;
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Create the flame stack with the fuel inlet, flow object, and oxidizer
|
||||
% inlet. Set the profile of the flame with the estimated axial velocities,
|
||||
% radial velocities, temperature, and mass fractions calculated above.
|
||||
%
|
||||
|
||||
flame = Stack([left flow right]);
|
||||
setProfile(flame, 2, {'u', 'V'}, [zrel; u; v]);
|
||||
setProfile(flame, 2, 'T', [zrel; t] );
|
||||
for n = 1:nsp
|
||||
nm = speciesName(tp_f, n);
|
||||
setProfile(flame, 2, nm, [zrel; transpose(y(:,n))])
|
||||
end
|
||||
end
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Define elMoles function
|
||||
%
|
||||
|
||||
function moles = elMoles(tp, element)
|
||||
% ELMOLES Determine the elemental moles in a gas object per unit mass.
|
||||
% moles = Moles(tp, element)
|
||||
% :param tp:
|
||||
% Object representing the gas, instance of class :mat:func:`Solution`,
|
||||
% and an ideal gas. The state of this object should be set to an
|
||||
% estimate of the gas state before calling Moles.
|
||||
% :param element:
|
||||
% String representing the element name.
|
||||
% :return:
|
||||
% Elemental moles within a gas object per unit mass of mixture.
|
||||
% Units: kmol/kg
|
||||
%
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Check input parameters
|
||||
%
|
||||
|
||||
if nargin ~= 2
|
||||
error('elMoles expects two input arguments.');
|
||||
end
|
||||
if ~isIdealGas(tp)
|
||||
error('Gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
if ~ischar(element)
|
||||
error('Element name must be of format character.');
|
||||
end
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Calculate the moles per mass of mixture of an element within a gas
|
||||
% object. The equation used is: elmoles = elMassFrac/Mel where elMassFrac
|
||||
% is the elemental mass fraction within the gas object using the
|
||||
% elementalMassFraction function; Mel is the atomic mass of the element.
|
||||
%
|
||||
|
||||
elMassFrac = elementalMassFraction(tp, element);
|
||||
eli = elementIndex(tp, element);
|
||||
M = atomicMasses(tp);
|
||||
Mel = M(eli);
|
||||
moles = elMassFrac/Mel;
|
||||
end
|
||||
|
|
@ -0,0 +1,51 @@
|
|||
function elMassFrac = elementalMassFraction(tp, element)
|
||||
% ELEMENTALMASSFRACTION Determine the elemental mass fraction in gas object.
|
||||
% elMassFrac = elementalMassFraction(tp, element)
|
||||
% :param tp:
|
||||
% Object representing the gas, instance of class :mat:func:`Solution`,
|
||||
% and an ideal gas. The state of this object should be set to an
|
||||
% estimate of the gas state before calling elementalMassFraction.
|
||||
% :param element:
|
||||
% String representing the element name.
|
||||
% :return:
|
||||
% Elemental mass fraction within a gas object.
|
||||
%
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Check input parameters
|
||||
%
|
||||
|
||||
if nargin ~= 2
|
||||
error('elementalMassFraction expects two input arguments.');
|
||||
end
|
||||
if ~isIdealGas(tp)
|
||||
error('Gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
if ~ischar(element)
|
||||
error('Element name must be of format character.');
|
||||
end
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Calculate the elemental mass fraction in a gas object. Equation used is
|
||||
% elMassFrac = sum of nAtoms(k,m)*Mel(m)*Y(k)/mw(k) where nAtoms(k,m) is
|
||||
% the number of atoms of element, m, in species, k; Mel(m) is the atomic
|
||||
% weight of the element, m; Y(k) is the mass fraction of species,k, in the
|
||||
% gas object; and mw(k) is the molecular weight of species, k.
|
||||
%
|
||||
|
||||
n = nSpecies(tp);
|
||||
massFrac = massFractions(tp);
|
||||
spec = speciesNames(tp);
|
||||
eli = elementIndex(tp, element);
|
||||
M = atomicMasses(tp);
|
||||
Mel = M(eli);
|
||||
MW = molecularWeights(tp);
|
||||
% Initialize the element mass fraction as zero.
|
||||
elMassFrac = 0.0;
|
||||
% Use loop to perform summation of elemental mass fraction over all species.
|
||||
for i = 1:n
|
||||
natoms(i) = nAtoms(tp,spec{i},element);
|
||||
mw(i) = MW(i);
|
||||
Y(i) = massFraction(tp,spec{i});
|
||||
elMassFrac = elMassFrac + (natoms(i)*Mel*Y(i))/mw(i);
|
||||
end
|
||||
end
|
||||
|
|
@ -1,119 +1,129 @@
|
|||
% DIFFFLAME - A non-premixed opposed-jet flame.
|
||||
%
|
||||
% DIFFFLAME - An opposed-flow diffusion flame.
|
||||
% This example uses the CounterFlowDiffusionFlame function to solve an
|
||||
% opposed-flow diffusion flame for Ethane in Air. This example is the same
|
||||
% as the diffusion_flame.py example without radiation.
|
||||
%
|
||||
|
||||
help diffflame
|
||||
disp('press any key to begin the simulation');
|
||||
pause
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%
|
||||
|
||||
t0 = cputime; % record the starting time
|
||||
runtime = cputime; % Record the starting time
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Parameter values of inlet streams
|
||||
%
|
||||
|
||||
% parameter values
|
||||
p = oneatm; % pressure
|
||||
tin = 300.0; % inlet temperature
|
||||
mdot_o = 0.72; % air, kg/m^2/s
|
||||
mdot_f = 0.24; % fuel, kg/m^2/s
|
||||
p = oneatm; % Pressure
|
||||
tin = 300.0; % Inlet temperature
|
||||
mdot_o = 0.72; % Air mass flux, kg/m^2/s
|
||||
mdot_f = 0.24; % Fuel mass flux, kg/m^2/s
|
||||
rxnmech = 'gri30.xml'; % Reaction mechanism file
|
||||
transport = 'Mix'; % Transport model
|
||||
% NOTE: Transport model needed if mechanism file does not have transport
|
||||
% properties.
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Set-up initial grid, loglevel, tolerances. Enable/Disable grid
|
||||
% refinement.
|
||||
%
|
||||
|
||||
rxnmech = 'gri30.xml'; % reaction mechanism file
|
||||
transport = 'Mix'; % transport model
|
||||
comp1 = 'O2:0.21, N2:0.78, AR:0.01'; % air composition
|
||||
comp2 = 'C2H6:1'; % fuel composition
|
||||
|
||||
initial_grid = 0.02*[0.0 0.2 0.4 0.6 0.8 1.0]; % m
|
||||
|
||||
tol_ss = [1.0e-5 1.0e-9]; % [rtol atol] for steady-state
|
||||
% problem
|
||||
initial_grid = 0.02*[0.0 0.2 0.4 0.6 0.8 1.0]; % Units: m
|
||||
tol_ss = [1.0e-5 1.0e-9]; % [rtol atol] for steady-state problem
|
||||
tol_ts = [1.0e-3 1.0e-9]; % [rtol atol] for time stepping
|
||||
|
||||
loglevel = 1; % amount of diagnostic output (0
|
||||
% to 5)
|
||||
|
||||
refine_grid = 1; % 1 to enable refinement, 0 to
|
||||
% disable
|
||||
|
||||
%%%%%%%%%%%%%%%% create the gas object %%%%%%%%%%%%%%%%%%%%%%%%
|
||||
loglevel = 1; % Amount of diagnostic output (0 to 5)
|
||||
refine_grid = 1; % 1 to enable refinement, 0 to disable
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Create the gas objects for the fuel and oxidizer streams. These objects
|
||||
% will be used to evaluate all thermodynamic, kinetic, and transport
|
||||
% properties.
|
||||
%
|
||||
% This object will be used to evaluate all thermodynamic, kinetic,
|
||||
% and transport properties
|
||||
|
||||
fuel = GRI30('Mix');
|
||||
ox = GRI30('Mix');
|
||||
oxcomp = 'O2:0.21, N2:0.78'; % Air composition
|
||||
fuelcomp = 'C2H6:1'; % Fuel composition
|
||||
% Set each gas mixture state with the corresponding composition.
|
||||
set(fuel,'T', tin, 'P', p, 'X', fuelcomp);
|
||||
set(ox,'T',tin,'P',p,'X', oxcomp);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Set-up the flow object. For this problem, the AxisymmetricFlow model is
|
||||
% needed. Set the state of the flow as the fuel gas object. This is
|
||||
% arbitrary and is only used to initialize the flow object. Set the grid to
|
||||
% the initial grid defined prior, same for the tolerances.
|
||||
%
|
||||
gas = GRI30('Mix')
|
||||
|
||||
% set its state to that of the fuel (arbitrary)
|
||||
set(gas,'T', tin, 'P', p, 'X', comp2);
|
||||
|
||||
%%%%%%%%%%%%%%%% create the flow object %%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
f = AxisymmetricFlow(gas,'flow');
|
||||
|
||||
f = AxisymmetricFlow(fuel,'flow');
|
||||
set(f, 'P', p, 'grid', initial_grid);
|
||||
set(f, 'tol', tol_ss, 'tol-time', tol_ts);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Create the fuel and oxidizer inlet steams. Specify the temperature, mass
|
||||
% flux, and composition correspondingly.
|
||||
%
|
||||
|
||||
%%%%%%%%%%%%%%% create the air inlet %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%
|
||||
% The temperature, mass flux, and composition (relative molar) may be
|
||||
% specified.
|
||||
%
|
||||
% Set the oxidizer inlet.
|
||||
inlet_o = Inlet('air_inlet');
|
||||
set(inlet_o, 'T', tin, 'MassFlux', mdot_o, 'X', comp1);
|
||||
|
||||
%%%%%%%%%%%%%% create the fuel inlet %%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%
|
||||
set(inlet_o, 'T', tin, 'MassFlux', mdot_o, 'X', oxcomp);
|
||||
%
|
||||
% Set the fuel inlet.
|
||||
inlet_f = Inlet('fuel_inlet');
|
||||
set(inlet_f, 'T', tin, 'MassFlux', mdot_f, 'X', comp2);
|
||||
|
||||
%%%%%%%%%%%%% create the flame object %%%%%%%%%%%%
|
||||
set(inlet_f, 'T', tin, 'MassFlux', mdot_f, 'X', fuelcomp);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Once the inlets have been created, they can be assembled
|
||||
% to create the flame object. Function CounterFlorDiffusionFlame
|
||||
% (in Cantera/1D) sets up the initial guess for the solution using a
|
||||
% Burke-Schumann flame. The input parameters are: fuel inlet object, flow
|
||||
% object, oxidizer inlet object, fuel gas object, oxidizer gas object, and
|
||||
% the name of the oxidizer species as in character format.
|
||||
%
|
||||
% Once the component parts have been created, they can be assembled
|
||||
% to create the flame object. Function npflame_init (in Cantera/1D)
|
||||
% sets up the initial guess for the solution using a Burke-Schumann
|
||||
% flame.
|
||||
|
||||
fl = CounterFlowDiffusionFlame(inlet_f, f, inlet_o, fuel, ox, 'O2');
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Solve with fixed temperature profile first. Grid refinement is turned off
|
||||
% for this process in this example. To turn grid refinement on, change 0 to
|
||||
% 1 for last input is solve function.
|
||||
%
|
||||
fl = npflame_init(gas, inlet_f, f, inlet_o, 'C2H6', 'O2', 3.5);
|
||||
|
||||
% if the starting solution is to be read from a previously-saved
|
||||
% solution, uncomment this line and edit the file name and solution id.
|
||||
%restore(fl,'h2flame2.xml', 'energy')
|
||||
|
||||
% solve with fixed temperature profile first
|
||||
solve(fl, loglevel, 0); %refine_grid);
|
||||
|
||||
%%%%%%%%%%%% enable the energy equation %%%%%%%%%%%%%%%%%%%%%
|
||||
%
|
||||
% The energy equation will now be solved to compute the
|
||||
% temperature profile. We also tighten the grid refinement
|
||||
% criteria to get an accurate final solution.
|
||||
solve(fl, loglevel, 0);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Enable the energy equation. The energy equation will now be solved to
|
||||
% compute the temperature profile. We also tighten the grid refinement
|
||||
% criteria to get an accurate final solution. The explanation of the
|
||||
% setRefineCriteria function is located on cantera.org in the Matlab User's
|
||||
% Guide and can be accessed by help setRefineCriteria
|
||||
%
|
||||
|
||||
enableEnergy(f);
|
||||
setRefineCriteria(fl, 2, 200.0, 0.1, 0.2);
|
||||
solve(fl, loglevel, refine_grid);
|
||||
saveSoln(fl,'c2h6.xml','energy',['solution with energy equation']);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Show statistics of solution and elapsed time.
|
||||
%
|
||||
|
||||
%%%%%%%%%% show statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
writeStats(fl);
|
||||
elapsed = cputime - t0;
|
||||
elapsed = cputime - runtime;
|
||||
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
|
||||
disp(e);
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Make a single plot showing temperature and mass fraction of select
|
||||
% species along axial distance from fuel inlet to air inlet.
|
||||
%
|
||||
|
||||
%%%%%%%%%% make plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
z = grid(fl, 'flow'); % Get grid points of flow
|
||||
spec = speciesNames(fuel); % Get species names in gas
|
||||
T = solution(fl, 'flow', 'T'); % Get temperature solution
|
||||
for i = 1:length(spec)
|
||||
% Get mass fraction of all species from solution
|
||||
y(i,:) = solution(fl, 'flow', spec{i});
|
||||
end
|
||||
j = speciesIndex(fuel, 'O2'); % Get index of O2 in gas object
|
||||
k = speciesIndex(fuel, 'H2O'); % Get index of H2O in gas object
|
||||
l = speciesIndex(fuel, 'C2H6'); % Get index of C2H6 in gas object
|
||||
m = speciesIndex(fuel, 'CO2'); % Get index of CO2 in gas object
|
||||
clf;
|
||||
subplot(2,3,1);
|
||||
plotSolution(fl, 'flow', 'T');
|
||||
title('Temperature [K]');
|
||||
subplot(2,3,2);
|
||||
plotSolution(fl, 'flow', 'C2H6');
|
||||
title('C2H6 Mass Fraction');
|
||||
subplot(2,3,3);
|
||||
plotSolution(fl, 'flow', 'O2');
|
||||
title('O2 Mass Fraction');
|
||||
subplot(2,3,4);
|
||||
plotSolution(fl, 'flow', 'CH');
|
||||
title('CH Mass Fraction');
|
||||
subplot(2,3,5);
|
||||
plotSolution(fl, 'flow', 'V');
|
||||
title('Radial Velocity / Radius [s^-1]');
|
||||
subplot(2,3,6);
|
||||
plotSolution(fl, 'flow', 'u');
|
||||
title('Axial Velocity [m/s]');
|
||||
yyaxis left
|
||||
plot(z,T)
|
||||
xlabel('z (m)');
|
||||
ylabel('Temperature (K)');
|
||||
yyaxis right
|
||||
plot(z,y(j,:),'r',z,y(k,:),'g',z,y(l,:),'m',z,y(m,:),'b');
|
||||
ylabel('Mass Fraction');
|
||||
legend('T','O2','H2O','C2H6','CO2');
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue