cantera/samples/matlab/ignite.m
Ray Speth 2528df0f75 Reorganized source tree structure
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.
2012-02-12 02:27:14 +00:00

105 lines
2.6 KiB
Matlab
Executable file

function plotdata = ignite(g)
% IGNITE Zero-dimensional kinetics: adiabatic, constant pressure.
%
% This example solves the same problem as 'reactor1,' but does
% it using on of MATLAB's ODE integrators, rather than using the
% Cantera Reactor class.
%
help ignite
if nargin == 1 & isa(g,'solution')
gas = g;
else
gas = IdealGasMix('gri30.xml');
end
nsp = nSpecies(gas);
% set the initial conditions
set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
y0 = [intEnergy_mass(gas)
1.0/density(gas)
massFractions(gas)];
time_interval = [0 0.001];
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);
disp(['CPU time = ' num2str(cputime - t0)]);
plotdata = output(out,gas);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the functions below may be defined arbitrarily to set the reactor
% boundary conditions - the rate of change of volume, the heat
% flux, and the area.
% Rate of change of volume. Any arbirtrary function may be implemented.
% Input arguments:
% t time
% vol volume
% gas ideal gas object
function v = vdot(t, vol, gas)
%v = 0.0; %uncomment for constant volume
v = 1.e11 * (pressure(gas) - 101325.0); % holds pressure very
% close to 1 atm
% heat flux (W/m^2).
function q = heatflux(t, gas)
q = 0.0; % adiabatic
% surface area (m^2). Used only to compute heat transfer.
function a = area(t,vol)
a = 1.0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Since the solution variables used by the 'reactor' function are
% not necessarily those desired for output, this function is called
% after the integration is complete to generate the desired
% outputs.
function pv = output(s, gas)
times = s.x;
soln = s.y;
[m n] = size(times);
pv = zeros(nSpecies(gas) + 4, n);
set(gas,'T',1001.0,'P',oneatm);
for j = 1:n
ss = soln(:,j);
y = ss(3:end);
mass = sum(y);
u_mass = ss(1)/mass;
v_mass = ss(2)/mass;
setMassFractions(gas, y);
setState_UV(gas, [u_mass v_mass]);
pv(1,j) = times(j);
pv(2,j) = temperature(gas);
pv(3,j) = density(gas);
pv(4,j) = pressure(gas);
pv(5:end,j) = y;
end
% plot the temperature and OH mass fractions.
figure(1);
plot(pv(1,:),pv(2,:));
xlabel('time');
ylabel('Temperature');
title(['Final T = ' num2str(pv(2,end)) ' K']);
figure(2);
ioh = speciesIndex(gas,'OH');
plot(pv(1,:),pv(4+ioh,:));
xlabel('time');
ylabel('Mass Fraction');
title('OH Mass Fraction');