mirror of
https://github.com/Cantera/cantera.git
synced 2025-02-25 18:55:29 -06:00
Added and fixed Examples for the new toolbox
Also fix Examples/flame1.m
This commit is contained in:
Su Sun
Ray Speth
parent
8043b28315
commit
c9a4004ab5
239
samples/matlab_experimental/catcomb.m
Normal file
239
samples/matlab_experimental/catcomb.m
Normal file
@ -0,0 +1,239 @@
|
||||
% Catalytic combustion of a stagnation flow on a platinum surface
|
||||
%
|
||||
% This script solves a catalytic combustion problem. A stagnation flow
|
||||
% is set up, with a gas inlet 10 cm from a platinum surface at 900
|
||||
% K. The lean, premixed methane/air mixture enters at ~ 6 cm/s (0.06
|
||||
% kg/m2/s), and burns catalytically on the platinum surface. Gas-phase
|
||||
% chemistry is included too, and has some effect very near the
|
||||
% surface.
|
||||
%
|
||||
% The catalytic combustion mechanism is from Deutschman et al., 26th
|
||||
% Symp. (Intl.) on Combustion,1996 pp. 1747-1754
|
||||
|
||||
%% Initialization
|
||||
|
||||
help catcomb;
|
||||
|
||||
clear all
|
||||
close all
|
||||
cleanup
|
||||
clc
|
||||
|
||||
t0 = cputime; % record the starting time
|
||||
|
||||
%% Set parameter values
|
||||
|
||||
p = oneatm; % pressure
|
||||
tinlet = 300.0; % inlet temperature
|
||||
tsurf = 900.0; % surface temperature
|
||||
mdot = 0.06; % kg/m^2/s
|
||||
transport = 'mixture-averaged'; % transport model
|
||||
|
||||
% Solve first for a hydrogen/air case for use as the initial estimate for
|
||||
% the methane/air case.
|
||||
|
||||
% composition of the inlet premixed gas for the hydrogen/air case
|
||||
comp1 = 'H2:0.05, O2:0.21, N2:0.78, AR:0.01';
|
||||
|
||||
% composition of the inlet premixed gas for the methane/air case
|
||||
comp2 = 'CH4:0.095, O2:0.21, N2:0.78, AR:0.01';
|
||||
|
||||
% the initial grid, in meters. The inlet/surface separation is 10 cm.
|
||||
initial_grid = [0.0, 0.02, 0.04, 0.06, 0.08, 0.1]; % m
|
||||
|
||||
% numerical parameters
|
||||
tol_ss = {1.0e-8 1.0e-14}; % {rtol atol} for steady-state problem
|
||||
tol_ts = {1.0e-4 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
|
||||
%
|
||||
% This object will be used to evaluate all thermodynamic, kinetic,
|
||||
% and transport properties
|
||||
%
|
||||
% The gas phase will be taken from the definition of phase 'gas' in
|
||||
% input file 'ptcombust.yaml,' which is a stripped-down version of
|
||||
% GRI-Mech 3.0.
|
||||
|
||||
gas = Solution('ptcombust.yaml', 'gas', transport);
|
||||
gas.TPX = {tinlet, p, comp1};
|
||||
|
||||
%% create the interface object
|
||||
%
|
||||
% This object will be used to evaluate all surface chemical production
|
||||
% rates. It will be created from the interface definition 'Pt_surf'
|
||||
% in input file 'ptcombust.yaml,' which implements the reaction
|
||||
% mechanism of Deutschmann et al., 1995 for catalytic combustion on
|
||||
% platinum.
|
||||
|
||||
surf_phase = importInterface('ptcombust.yaml', 'Pt_surf', gas);
|
||||
surf_phase.T = tsurf;
|
||||
|
||||
% integrate the coverage equations in time for 1 s, holding the gas
|
||||
% composition fixed to generate a good starting estimate for the
|
||||
% coverages.
|
||||
|
||||
surf_phase.advanceCoverages(1.0);
|
||||
|
||||
% The two objects we just created are independent of the problem
|
||||
% type -- they are useful in zero-D simulations, 1-D simulations,
|
||||
% etc. Now we turn to creating the objects that are specifically
|
||||
% for 1-D simulations. These will be 'stacked' together to create
|
||||
% the complete simulation.
|
||||
|
||||
%% create the flow object
|
||||
%
|
||||
% The flow object is responsible for evaluating the 1D governing
|
||||
% equations for the flow. We will initialize it with the gas
|
||||
% object, and assign it the name 'flow'.
|
||||
|
||||
flow = AxisymmetricFlow(gas, 'flow');
|
||||
|
||||
% set some parameters for the flow
|
||||
flow.setPressure(p);
|
||||
flow.setupGrid(initial_grid);
|
||||
flow.setSteadyTolerances('default', tol_ss{:});
|
||||
flow.setTransientTolerances('default', tol_ts{:});
|
||||
|
||||
%% create the inlet
|
||||
%
|
||||
% The temperature, mass flux, and composition (relative molar) may be
|
||||
% specified. This object provides the inlet boundary conditions for
|
||||
% the flow equations.
|
||||
|
||||
inlt = Inlet('inlet');
|
||||
|
||||
% set the inlet parameters. Start with comp1 (hydrogen/air)
|
||||
inlt.T = tinlet;
|
||||
inlt.setMdot(mdot);
|
||||
inlt.setMoleFractions(comp1);
|
||||
|
||||
%% create the surface
|
||||
%
|
||||
% This object provides the surface boundary conditions for the flow
|
||||
% equations. By supplying object surface_phase as an argument, the
|
||||
% coverage equations for its surface species will be added to the
|
||||
% equation set, and used to compute the surface production rates of
|
||||
% the gas-phase species.
|
||||
|
||||
surf = Surface('surface', surf_phase);
|
||||
surf.T = tsurf;
|
||||
|
||||
%% create the stack
|
||||
%
|
||||
% Once the component parts have been created, they can be assembled
|
||||
% to create the 1D simulation.
|
||||
|
||||
sim1D = Stack([inlt, flow, surf]);
|
||||
|
||||
% set the initial profiles.
|
||||
sim1D.setProfile(2, {'u', 'V', 'T'}, [0.0, 1.0 % z/zmax
|
||||
0.06, 0.0 % u
|
||||
0.0, 0.0 % V
|
||||
tinlet, tsurf]); % T
|
||||
names = gas.speciesNames;
|
||||
|
||||
for k = 1:gas.nSpecies
|
||||
y = inlt.massFraction(k);
|
||||
sim1D.setProfile(2, names{k}, [0, 1; y, y]);
|
||||
end
|
||||
|
||||
sim1D
|
||||
|
||||
%setTimeStep(fl, 1.0e-5, [1, 3, 6, 12]);
|
||||
%setMaxJacAge(fl, 4, 5);
|
||||
|
||||
%% Solution
|
||||
|
||||
% start with the energy equation on
|
||||
flow.enableEnergy;
|
||||
|
||||
% disable the surface coverage equations, and turn off all gas and
|
||||
% surface chemistry
|
||||
|
||||
surf.setCoverageEqs('off');
|
||||
surf_phase.setMultiplier(0.0);
|
||||
gas.setMultiplier(0.0);
|
||||
|
||||
% solve the problem, refining the grid if needed
|
||||
sim1D.solve(1, refine_grid);
|
||||
|
||||
% now turn on the surface coverage equations, and turn the
|
||||
% chemistry on slowly
|
||||
|
||||
surf.setCoverageEqs('on');
|
||||
|
||||
for iter=1:6
|
||||
mult = 10.0^(iter - 6);
|
||||
surf_phase.setMultiplier(mult);
|
||||
gas.setMultiplier(mult);
|
||||
sim1D.solve(1, refine_grid);
|
||||
end
|
||||
|
||||
% At this point, we should have the solution for the hydrogen/air
|
||||
% problem. Now switch the inlet to the methane/air composition.
|
||||
inlt.setMoleFractions(comp2);
|
||||
|
||||
% set more stringent grid refinement criteria
|
||||
sim1D.setRefineCriteria(2, 100.0, 0.15, 0.2);
|
||||
|
||||
% solve the problem for the final time
|
||||
sim1D.solve(loglevel, refine_grid);
|
||||
|
||||
% show the solution
|
||||
sim1D
|
||||
|
||||
% save the solution
|
||||
sim1D.saveSoln('catcomb.xml', 'energy', ['solution with energy equation']);
|
||||
|
||||
%% Show statistics
|
||||
|
||||
sim1D.writeStats;
|
||||
elapsed = cputime - t0;
|
||||
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
|
||||
disp(e);
|
||||
|
||||
%% Make plots
|
||||
|
||||
clf;
|
||||
|
||||
subplot(3, 3, 1);
|
||||
sim1D.plotSolution('flow', 'T');
|
||||
title('Temperature [K]');
|
||||
|
||||
subplot(3, 3, 2);
|
||||
sim1D.plotSolution('flow', 'u');
|
||||
title('Axial Velocity [m/s]');
|
||||
|
||||
subplot(3, 3, 3);
|
||||
sim1D.plotSolution('flow', 'V');
|
||||
title('Radial Velocity / Radius [1/s]');
|
||||
|
||||
subplot(3, 3, 4);
|
||||
sim1D.plotSolution('flow', 'CH4');
|
||||
title('CH4 Mass Fraction');
|
||||
|
||||
subplot(3, 3, 5);
|
||||
sim1D.plotSolution('flow', 'O2');
|
||||
title('O2 Mass Fraction');
|
||||
|
||||
subplot(3, 3, 6);
|
||||
sim1D.plotSolution('flow', 'CO');
|
||||
title('CO Mass Fraction');
|
||||
|
||||
subplot(3, 3, 7);
|
||||
sim1D.plotSolution('flow', 'CO2');
|
||||
title('CO2 Mass Fraction');
|
||||
|
||||
subplot(3, 3, 8);
|
||||
sim1D.plotSolution('flow', 'H2O');
|
||||
title('H2O Mass Fraction');
|
||||
|
||||
subplot(3, 3, 9);
|
||||
sim1D.plotSolution('flow', 'H2');
|
||||
title('H2 Mass Fraction');
|
||||
27
samples/matlab_experimental/conhp.m
Normal file
27
samples/matlab_experimental/conhp.m
Normal file
@ -0,0 +1,27 @@
|
||||
function dydt = conhp(t, y, gas, mw) %#ok<INUSL>
|
||||
% CONHP - ODE system for a constant-pressure, adiabatic reactor.
|
||||
%
|
||||
% Function CONHP evaluates the system of ordinary differential
|
||||
% equations for an adiabatic, constant-pressure,
|
||||
% zero-dimensional reactor. It assumes that the 'gas' object
|
||||
% represents a reacting ideal gas mixture.
|
||||
|
||||
% Set the state of the gas, based on the current solution vector.
|
||||
gas.Y = y(2: end);
|
||||
gas.TP = {y(1), gas.P};
|
||||
nsp = gas.nSpecies;
|
||||
|
||||
% energy equation
|
||||
wdot = gas.netProdRates;
|
||||
gas.basis = 'mass';
|
||||
tdot = - gas.T * gasconstant * gas.enthalpies_RT .* wdot / (gas.D * gas.cp);
|
||||
|
||||
% set up column vector for dydt
|
||||
dydt = [tdot'
|
||||
zeros(nsp, 1)];
|
||||
|
||||
% species equations
|
||||
rrho = 1.0/gas.D;
|
||||
for i = 1:nsp
|
||||
dydt(i+1) = rrho * mw(i) * wdot(i);
|
||||
end
|
||||
31
samples/matlab_experimental/conuv.m
Normal file
31
samples/matlab_experimental/conuv.m
Normal file
@ -0,0 +1,31 @@
|
||||
function dydt = conuv(t, y, gas, mw) %#ok<INUSL>
|
||||
% CONUV ODE system for a constant-volume, adiabatic reactor.
|
||||
%
|
||||
% Function CONUV evaluates the system of ordinary differential
|
||||
% equations for an adiabatic, constant-volume,
|
||||
% zero-dimensional reactor. It assumes that the 'gas' object
|
||||
% represents a reacting ideal gas mixture.
|
||||
|
||||
|
||||
% Set the state of the gas, based on the current solution vector.
|
||||
gas.Y = y(2:end);
|
||||
gas.TD = {y(1), gas.D};
|
||||
nsp = gas.nSpecies;
|
||||
|
||||
% energy equation
|
||||
wdot = gas.netProdRates;
|
||||
gas.basis = 'mass';
|
||||
tdot = - gas.T * gasconstant * (gas.enthalpies_RT - ones(1, nsp)) ...
|
||||
.* wdot / (gas.D * gas.cv);
|
||||
|
||||
% set up column vector for dydt
|
||||
dydt = [tdot'
|
||||
zeros(nsp, 1)];
|
||||
|
||||
% species equations
|
||||
rrho = 1.0/gas.D;
|
||||
for i = 1:nsp
|
||||
dydt(i+1) = rrho * mw(i) * wdot(i);
|
||||
end
|
||||
|
||||
end
|
||||
141
samples/matlab_experimental/diff_flame.m
Normal file
141
samples/matlab_experimental/diff_flame.m
Normal file
@ -0,0 +1,141 @@
|
||||
% DIFF_FLAME - 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.
|
||||
|
||||
%% Initialization
|
||||
|
||||
help diff_flame
|
||||
|
||||
clear all
|
||||
close all
|
||||
cleanup
|
||||
clc
|
||||
|
||||
runtime = cputime; % Record the starting time
|
||||
|
||||
%% Parameter values of inlet streams
|
||||
|
||||
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
|
||||
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
|
||||
|
||||
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 objects for the fuel and oxidizer streams
|
||||
%
|
||||
% These objects will be used to evaluate all thermodynamic, kinetic, and
|
||||
% transport properties.
|
||||
|
||||
fuel = GRI30(transport);
|
||||
ox = GRI30(transport);
|
||||
oxcomp = 'O2:0.21, N2:0.78'; % Air composition
|
||||
fuelcomp = 'C2H6:1'; % Fuel composition
|
||||
% Set each gas mixture state with the corresponding composition.
|
||||
fuel.TPX = {tin, p, fuelcomp};
|
||||
ox.TPX = {tin, p, 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.
|
||||
|
||||
f = AxisymmetricFlow(fuel, 'flow');
|
||||
f.setPressure(p);
|
||||
f.setupGrid(initial_grid);
|
||||
f.setSteadyTolerances('default', tol_ss{:});
|
||||
f.setTransientTolerances('default', tol_ts{:});
|
||||
|
||||
%% Create the fuel and oxidizer inlet steams
|
||||
%
|
||||
% Specify the temperature, mass flux, and composition correspondingly.
|
||||
|
||||
% Set the oxidizer inlet.
|
||||
inlet_o = Inlet('air_inlet');
|
||||
inlet_o.T = tin;
|
||||
inlet_o.setMdot(mdot_o);
|
||||
inlet_o.setMoleFractions(oxcomp);
|
||||
|
||||
% Set the fuel inlet.
|
||||
inlet_f = Inlet('fuel_inlet');
|
||||
inlet_f.T = tin;
|
||||
inlet_f.setMdot(mdot_f);
|
||||
inlet_f.setMoleFractions(fuelcomp);
|
||||
|
||||
%% Create the flame object.
|
||||
%
|
||||
% 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.
|
||||
|
||||
fl = CounterFlowDiffusionFlame(inlet_f, f, inlet_o, fuel, ox, 'O2');
|
||||
|
||||
%% Solve with fixed temperature profile
|
||||
%
|
||||
% 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.solve(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
|
||||
|
||||
f.enableEnergy;
|
||||
fl.setRefineCriteria(2, 200.0, 0.1, 0.2);
|
||||
fl.solve(loglevel, refine_grid);
|
||||
fl.saveSoln('c2h6.xml', 'energy', ['solution with energy equation']);
|
||||
|
||||
%% Show statistics of solution and elapsed time.
|
||||
|
||||
fl.writeStats;
|
||||
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.
|
||||
%
|
||||
|
||||
z = fl.grid('flow'); % Get grid points of flow
|
||||
spec = fuel.speciesNames; % Get species names in gas
|
||||
T = fl.solution('flow', 'T'); % Get temperature solution
|
||||
|
||||
for i = 1:length(spec)
|
||||
% Get mass fraction of all species from solution
|
||||
y(i, :) = fl.solution('flow', spec{i});
|
||||
end
|
||||
|
||||
j = fuel.speciesIndex('O2'); % Get index of O2 in gas object
|
||||
k = fuel.speciesIndex('H2O'); % Get index of H2O in gas object
|
||||
l = fuel.speciesIndex('C2H6'); % Get index of C2H6 in gas object
|
||||
m = fuel.speciesIndex('CO2'); % Get index of CO2 in gas object
|
||||
|
||||
clf;
|
||||
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');
|
||||
@ -6,20 +6,20 @@ function f = flame(gas, left, flow, right)
|
||||
error('wrong number of input arguments.');
|
||||
end
|
||||
|
||||
if ~gas.thermo.isIdealGas
|
||||
if ~gas.isIdealGas
|
||||
error('gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
if ~isInlet(left)
|
||||
if ~left.isInlet
|
||||
error('burner object of wrong type.');
|
||||
end
|
||||
if ~isFlow(flow)
|
||||
if ~flow.isFlow
|
||||
error('flow object of wrong type.');
|
||||
end
|
||||
|
||||
flametype = 0;
|
||||
if isSurface(right)
|
||||
if right.isSurface
|
||||
flametype = 1;
|
||||
elseif isInlet(right)
|
||||
elseif right.isInlet
|
||||
flametype = 3;
|
||||
end
|
||||
|
||||
|
||||
122
samples/matlab_experimental/flame1.m
Normal file
122
samples/matlab_experimental/flame1.m
Normal file
@ -0,0 +1,122 @@
|
||||
% FLAME1 - A burner-stabilized flat flame
|
||||
%
|
||||
% This script simulates a burner-stablized lean hydrogen-oxygen flame
|
||||
% at low pressure.
|
||||
|
||||
%% Initialization
|
||||
|
||||
help flame1
|
||||
|
||||
clear all
|
||||
close all
|
||||
cleanup
|
||||
clc
|
||||
|
||||
t0 = cputime; % record the starting time
|
||||
|
||||
%% Set parameter values
|
||||
|
||||
p = 0.05*oneatm; % pressure
|
||||
tburner = 373.0; % burner temperature
|
||||
mdot = 0.06; % kg/m^2/s
|
||||
|
||||
rxnmech = 'h2o2.yaml'; % reaction mechanism file
|
||||
comp = 'H2:1.8, O2:1, AR:7'; % premixed gas composition
|
||||
|
||||
initial_grid = [0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.15, 0.2, 0.4,...
|
||||
0.49, 0.5]; % m
|
||||
|
||||
tol_ss = {1.0e-5, 1.0e-13}; % {rtol atol} for steady-state
|
||||
% problem
|
||||
tol_ts = {1.0e-4, 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
|
||||
max_jacobian_age = [5, 10];
|
||||
|
||||
%% Create the gas object
|
||||
%
|
||||
% This object will be used to evaluate all thermodynamic, kinetic,
|
||||
% and transport properties
|
||||
|
||||
gas = Solution(rxnmech, 'ohmech', 'Mix');
|
||||
|
||||
% set its state to that of the unburned gas at the burner
|
||||
gas.TPX = {tburner, p, comp};
|
||||
|
||||
%% Create the flow object
|
||||
|
||||
f = AxisymmetricFlow(gas, 'flow');
|
||||
f.setPressure(p);
|
||||
f.setupGrid(initial_grid);
|
||||
f.setSteadyTolerances('default', tol_ss{:});
|
||||
f.setTransientTolerances('default', tol_ts{:});
|
||||
|
||||
%% Create the burner
|
||||
%
|
||||
% The burner is an Inlet object. The temperature, mass flux,
|
||||
% and composition (relative molar) may be specified.
|
||||
burner = Inlet('burner');
|
||||
burner.T = tburner;
|
||||
burner.setMdot(mdot);
|
||||
burner.setMoleFractions(comp);
|
||||
|
||||
%% Create the outlet
|
||||
%
|
||||
% The type of flame is determined by the object that terminates
|
||||
% the domain. An Outlet object imposes zero gradient boundary
|
||||
% conditions for the temperature and mass fractions, and zero
|
||||
% radial velocity and radial pressure gradient.
|
||||
|
||||
s = Outlet('out');
|
||||
|
||||
%% Create the flame object
|
||||
%
|
||||
% Once the component parts have been created, they can be assembled
|
||||
% to create the flame object.
|
||||
%
|
||||
fl = flame(gas, burner, f, s);
|
||||
fl.setMaxJacAge(max_jacobian_age(1), max_jacobian_age(2));
|
||||
|
||||
% 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')
|
||||
|
||||
fl.solve(loglevel, 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.
|
||||
|
||||
f.enableEnergy;
|
||||
fl.setRefineCriteria(2, 200.0, 0.05, 0.1);
|
||||
fl.solve(1, 1);
|
||||
fl.saveSoln('h2fl.xml', 'energy', ['solution with energy equation']);
|
||||
|
||||
%% Show statistics
|
||||
|
||||
fl.writeStats;
|
||||
elapsed = cputime - t0;
|
||||
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
|
||||
disp(e);
|
||||
|
||||
%% Make plots
|
||||
|
||||
clf;
|
||||
subplot(2, 2, 1);
|
||||
plotSolution(fl, 'flow', 'T');
|
||||
title('Temperature [K]');
|
||||
subplot(2, 2, 2);
|
||||
plotSolution(fl, 'flow', 'u');
|
||||
title('Axial Velocity [m/s]');
|
||||
subplot(2, 2, 3);
|
||||
plotSolution(fl, 'flow', 'H2O');
|
||||
title('H2O Mass Fraction');
|
||||
subplot(2, 2, 4);
|
||||
plotSolution(fl, 'flow', 'O2');
|
||||
title('O2 Mass Fraction');
|
||||
126
samples/matlab_experimental/flame2.m
Normal file
126
samples/matlab_experimental/flame2.m
Normal file
@ -0,0 +1,126 @@
|
||||
% FLAME2 - An axisymmetric stagnation-point non-premixed flame
|
||||
%
|
||||
% This script simulates a stagnation-point ethane-air flame.
|
||||
|
||||
%% Initialization
|
||||
|
||||
help flame2
|
||||
|
||||
clear all
|
||||
close all
|
||||
cleanup
|
||||
clc
|
||||
|
||||
t0 = cputime; % record the starting time
|
||||
|
||||
%% Set 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
|
||||
|
||||
rxnmech = 'gri30.yaml'; % reaction mechanism file
|
||||
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-13}; % {rtol atol} for steady-state
|
||||
% problem
|
||||
tol_ts = {1.0e-4, 1.0e-13}; % {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
|
||||
%
|
||||
% This object will be used to evaluate all thermodynamic, kinetic,
|
||||
% and transport properties
|
||||
|
||||
gas = Solution(rxnmech, 'gri30', 'Mix');
|
||||
|
||||
% set its state to that of the fuel (arbitrary)
|
||||
gas.TPX = {tin, p, comp2};
|
||||
|
||||
%% Create the flow object
|
||||
|
||||
f = AxisymmetricFlow(gas,'flow');
|
||||
f.setPressure(p);
|
||||
f.setupGrid(initial_grid);
|
||||
f.setSteadyTolerances('default', tol_ss{:});
|
||||
f.setTransientTolerances('default', tol_ts{:});
|
||||
|
||||
%% Create the air inlet
|
||||
%
|
||||
% The temperature, mass flux, and composition (relative molar) may be
|
||||
% specified.
|
||||
|
||||
inlet_o = Inlet('air_inlet');
|
||||
inlet_o.T = tin;
|
||||
inlet_o.setMdot(mdot_o);
|
||||
inlet_o.setMoleFractions(comp1);
|
||||
|
||||
%% Create the fuel inlet
|
||||
|
||||
inlet_f = Inlet('fuel_inlet');
|
||||
inlet_f.T = tin;
|
||||
inlet_f.setMdot(mdot_f);
|
||||
inlet_f.setMoleFractions(comp2);
|
||||
|
||||
%% Create the flame object
|
||||
%
|
||||
% Once the component parts have been created, they can be assembled
|
||||
% to create the flame object.
|
||||
|
||||
fl = flame(gas, inlet_o, f, inlet_f);
|
||||
|
||||
% 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
|
||||
fl.solve(loglevel, 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.
|
||||
|
||||
f.enableEnergy;
|
||||
fl.setRefineCriteria(2, 200.0, 0.1, 0.1);
|
||||
fl.solve(loglevel, refine_grid);
|
||||
fl.saveSoln('c2h6.xml', 'energy', ['solution with energy equation']);
|
||||
|
||||
%% Show statistics
|
||||
|
||||
fl.writeStats;
|
||||
elapsed = cputime - t0;
|
||||
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
|
||||
disp(e);
|
||||
|
||||
%% Make plots
|
||||
|
||||
figure(1);
|
||||
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]');
|
||||
@ -16,20 +16,21 @@ function plotdata = ignite(g)
|
||||
|
||||
% set the initial conditions
|
||||
|
||||
gas.TPX = {1001.0, oneatm, 'X','H2:2,O2:1,N2:4'};
|
||||
gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'};
|
||||
gas.basis = 'mass';
|
||||
y0 = [gas.U
|
||||
1.0/gas.D
|
||||
gas.Y'];
|
||||
|
||||
time_interval = [0 0.001];
|
||||
options = odeset('RelTol',1.e-5,'AbsTol',1.e-12,'Stats','on');
|
||||
options = odeset('RelTol', 1.e-5, 'AbsTol', 1.e-12, 'Stats', 'on');
|
||||
|
||||
t0 = cputime;
|
||||
out = ode15s(@reactor_ode,time_interval,y0,options,gas,@vdot,@area,@heatflux);
|
||||
out = ode15s(@reactor_ode, time_interval, y0, options, gas, ...
|
||||
@vdot, @area, @heatflux);
|
||||
disp(['CPU time = ' num2str(cputime - t0)]);
|
||||
|
||||
plotdata = output(out,gas);
|
||||
plotdata = output(out, gas);
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% the functions below may be defined arbitrarily to set the reactor
|
||||
@ -71,7 +72,7 @@ function plotdata = ignite(g)
|
||||
gas.TP = {1001.0, oneatm};
|
||||
|
||||
for j = 1:n
|
||||
ss = soln(:,j);
|
||||
ss = soln(:, j);
|
||||
y = ss(3:end);
|
||||
mass = sum(y);
|
||||
u_mass = ss(1)/mass;
|
||||
@ -79,11 +80,11 @@ function plotdata = ignite(g)
|
||||
gas.Y = y;
|
||||
gas.UV = {u_mass, v_mass};
|
||||
|
||||
pv(1,j) = times(j);
|
||||
pv(2,j) = gas.T;
|
||||
pv(3,j) = gas.D;
|
||||
pv(4,j) = gas.P;
|
||||
pv(5:end,j) = y;
|
||||
pv(1, j) = times(j);
|
||||
pv(2, j) = gas.T;
|
||||
pv(3, j) = gas.D;
|
||||
pv(4, j) = gas.P;
|
||||
pv(5:end, j) = y;
|
||||
end
|
||||
|
||||
% plot the temperature and OH mass fractions.
|
||||
@ -95,7 +96,7 @@ function plotdata = ignite(g)
|
||||
|
||||
figure(2);
|
||||
ioh = gas.speciesIndex('OH');
|
||||
plot(pv(1,:),pv(4+ioh,:));
|
||||
plot(pv(1, :), pv(4+ioh, :));
|
||||
xlabel('time');
|
||||
ylabel('Mass Fraction');
|
||||
title('OH Mass Fraction');
|
||||
|
||||
38
samples/matlab_experimental/ignite_hp.m
Normal file
38
samples/matlab_experimental/ignite_hp.m
Normal file
@ -0,0 +1,38 @@
|
||||
function ignite_hp(gas)
|
||||
% IGNITE_HP Solves the same ignition problem as 'ignite', but uses
|
||||
% function conhp instead of reactor.
|
||||
|
||||
help ignite_hp
|
||||
|
||||
if nargin == 0
|
||||
gas = Solution('gri30.yaml');
|
||||
end
|
||||
|
||||
mw = gas.MolecularWeights;
|
||||
gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'};
|
||||
|
||||
y0 = [gas.T
|
||||
gas.X'];
|
||||
tel = [0, 0.001];
|
||||
options = odeset('RelTol', 1.e-5, 'AbsTol', 1.e-12, 'Stats', 'on');
|
||||
t0 = cputime;
|
||||
out = ode15s(@conhp, tel, y0, options, gas, mw);
|
||||
disp(['CPU time = ' num2str(cputime - t0)]);
|
||||
|
||||
if nargout == 0
|
||||
% plot the temperature and OH mole fractions.
|
||||
figure(1);
|
||||
plot(out.x, out.y(1,:));
|
||||
xlabel('time');
|
||||
ylabel('Temperature');
|
||||
title(['Final T = ' num2str(out.y(1, end)), ' K']);
|
||||
|
||||
figure(2);
|
||||
ioh = gas.speciesIndex('OH');
|
||||
plot(out.x, out.y(1+ioh, :));
|
||||
xlabel('time');
|
||||
ylabel('Mass Fraction');
|
||||
title('OH Mass Fraction');
|
||||
end
|
||||
|
||||
end
|
||||
38
samples/matlab_experimental/ignite_uv.m
Normal file
38
samples/matlab_experimental/ignite_uv.m
Normal file
@ -0,0 +1,38 @@
|
||||
function ignite_uv(gas)
|
||||
% IGNITE_UV Solves the same ignition problem as 'ignite2', except
|
||||
% function conuv is used instead of reactor.
|
||||
%
|
||||
help ignite_uv
|
||||
|
||||
if nargin == 0
|
||||
gas = Solution('gri30.yaml');
|
||||
end
|
||||
|
||||
mw = gas.MolecularWeights;
|
||||
gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'};
|
||||
|
||||
y0 = [gas.T
|
||||
gas.X'];
|
||||
tel = [0, 0.001];
|
||||
options = odeset('RelTol', 1.e-5, 'AbsTol', 1.e-12, 'Stats', 'on');
|
||||
t0 = cputime;
|
||||
out = ode15s(@conuv, tel, y0, options, gas, mw);
|
||||
disp(['CPU time = ' num2str(cputime - t0)]);
|
||||
|
||||
if nargout == 0
|
||||
% plot the temperature and OH mole fractions.
|
||||
figure(1);
|
||||
plot(out.x, out.y(1, :));
|
||||
xlabel('time');
|
||||
ylabel('Temperature');
|
||||
title(['Final T = ' num2str(out.y(1, end)), ' K']);
|
||||
|
||||
figure(2);
|
||||
ioh = gas.speciesIndex('OH');
|
||||
plot(out.x, out.y(1+ioh, :));
|
||||
xlabel('time');
|
||||
ylabel('Mass Fraction');
|
||||
title('OH Mass Fraction');
|
||||
end
|
||||
|
||||
end
|
||||
@ -7,6 +7,11 @@ function prandtl1(g)
|
||||
|
||||
help prandtl1
|
||||
|
||||
clear all
|
||||
close all
|
||||
cleanup
|
||||
clc
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
else
|
||||
|
||||
@ -6,6 +6,11 @@ function prandtl2(g)
|
||||
%
|
||||
help prandtl2
|
||||
|
||||
clear all
|
||||
close all
|
||||
cleanup
|
||||
clc
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
else
|
||||
|
||||
@ -21,10 +21,10 @@ gas.TPX = {t, oneatm, 'CH4:0.01, O2:0.21, N2:0.78'};
|
||||
% methane on platinum, and is from Deutschman et al., 26th
|
||||
% Symp. (Intl.) on Combustion,1996, pp. 1747-1754
|
||||
surf = importInterface('ptcombust.yaml','Pt_surf', gas);
|
||||
surf.th.T = t;
|
||||
surf.T = t;
|
||||
|
||||
nsp = gas.nSpecies;
|
||||
nSurfSp = surf.th.nSpecies;
|
||||
nSurfSp = surf.nSpecies;
|
||||
|
||||
% create a reactor, and insert the gas
|
||||
r = IdealGasReactor(gas);
|
||||
@ -60,9 +60,9 @@ nSteps = 100;
|
||||
p0 = r.P;
|
||||
names = {'CH4','CO','CO2','H2O'};
|
||||
x = zeros([nSteps 4]);
|
||||
tim = zeros(nSteps);
|
||||
temp = zeros(nSteps);
|
||||
pres = zeros(nSteps);
|
||||
tim = zeros(nSteps, 1);
|
||||
temp = zeros(nSteps, 1);
|
||||
pres = zeros(nSteps, 1);
|
||||
cov = zeros([nSteps nSurfSp]);
|
||||
t = 0;
|
||||
dt = 0.1;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user