GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

[MATLAB-legacy] Remove samples

Browse Source
This commit is contained in:
Ingmar Schoegl 2024-02-20 19:30:02 -06:00 committed by Ray Speth
parent e9215e6600
commit 4b4e790f5a
32 changed files with 2 additions and 2647 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
samples/matlab/PFR_solver.m
View File

@ -1,53 +0,0 @@
function [F] = PFR_solver(x,soln_vector,gas,mdot,A_in,dAdx,k)
% This function defines the spatial derivatives for an ideal gas plug-flow
% reactor, where the cross-sectional area and pressure are allowed to vary,
% axially. The model is set up by the example file 'Plug_Flow_Reactor.m',
% which points the integrator to this function. The integrator integrates the
% derivatives spatially, to solve the density, temperature, and species mass
% fraction profiles as a function of distance x.
rho = soln_vector(1);
T = soln_vector(2);
Y = soln_vector(3:end);
if k==1
A = A_in+k*dAdx*x;
elseif k==-1
A = A_in+k*dAdx*x;
dAdx = -dAdx;
else
A = A_in+k*dAdx*x;
end
% the gas is set to the corresponding properties during each iteration of the ode loop
set(gas,'Temperature',T,'Density',rho,'MassFractions',Y);
MW_mix = meanMolecularWeight(gas);
Ru = gasconstant();
R = Ru/MW_mix;
nsp = nSpecies(gas);
vx = mdot/(rho*A);
P = rho*R*T;
MW = molecularWeights(gas);
h = enthalpies_RT(gas).*R.*T;
w = netProdRates(gas);
Cp = cp_mass(gas);
%--------------------------------------------------------------------------
%---F(1), F(2) and F(3:end) are the differential equations modelling the---
%---density, temperature and mass fractions variations along a plug flow---
%-------------------------reactor------------------------------------------
%--------------------------------------------------------------------------
F(1) = ((1-R/Cp)*((rho*vx)^2)*(1/A)*(dAdx)...
+ rho*R*sum(MW.*w.*(h-MW_mix*Cp*T./MW))/(vx*Cp) )...
/ (P*(1+vx^2/(Cp*T)) - rho*vx^2);
F(2) = (vx*vx/(rho*Cp))*F(1) + vx*vx*(1/A)*(dAdx)/Cp...
- (1/(vx*rho*Cp))*sum(h.*w.*MW);
F(3:nsp+2) = w(1:nsp).*MW(1:nsp)./(rho*vx);
F = F';
end

165
samples/matlab/Plug_Flow_Reactor.m
View File

@ -1,165 +0,0 @@
% General_Plug_Flow_Reactor (PFR) - to solve PFR equations for reactors
%
% This code snippet is to model a constant area and varying area
% (converging and diverging) nozzle as Plug Flow Reactor with given
% dimensions and an incoming gas. The pressure is not assumed to be
% constant here, as opposed to the Python Version of the
% Plug Flow Reactor model.
%
% The reactor assumes that the flow follows the Ideal Gas Law.
%
% The governing equations used in this code can be referenced at:
% *S.R Turns, An Introduction to Combustion - Concepts and Applications,
% McGraw Hill Education, India, 2012, 206-210.*
%
% The current example is written for methane combustion, but can be readily
% adapted for other chemistries.
%
% Developed by Ashwin Kumar/Dr.Joseph Meadows (mgak@vt.edu/jwm84@vt.edu) on 3-June-2020
% Research Assistant/Assistant Professor
% Advanced Propulsion and Power Laboratory
% Virginia Tech
%
% Keywords: combustion, user-defined model, compressible flow, plotting
%% Clear all variables, close all figures, clear the command line:
clear all
close all
clc
% Temperature of gas, in K
T0 = 1473;
% Pressure of gas, in Pa
P0 = 4.47*101325;
% Equivalence Ratio
Phi = 0.2899;
% Import the gas phase, read out key species indices:
gas_calc = Solution('gri30.yaml');
ich4 = speciesIndex(gas_calc,'CH4');
io2 = speciesIndex(gas_calc,'O2');
in2 = speciesIndex(gas_calc,'N2');
nsp = nSpecies(gas_calc);
x = zeros(nsp,1);
% Change the below values for different Phi values of methane Combustion
x(ich4,1) = Phi;
x(io2,1) = 2.0;
x(in2,1) = 7.52;
% Set the initial state and then equilibrate for a given enthalpy and pressure:
set(gas_calc,'T',T0,'P',P0,'MoleFractions',x);
gas_calc = equilibrate(gas_calc,'HP');
%% Calculation of properties along the reactor length
% The Dimensions and conditions of the reactor are given below
% Inlet Area, in m^2
A_in = 0.018;
% Exit Area, in m^2
A_out = 0.003;
% Length of the reactor, in m
L = 1.284*0.0254;
% The whole reactor is divided into n small reactors
n = 100;
% Mass flow rate into the reactor, in kg/s
mdot_calc = 1.125;
% Flag to indicate whether the area is converging, diverging, or constant:
% k = -1 makes the solver solve for converging area.
% k = +1 makes the solver solve for diverging area.
% k = 0 makes the solver solve for constant cross sectional area
if A_in>A_out
k = -1;
elseif A_out>A_in
k = 1;
else
k = 0;
end
dAdx = abs(A_in-A_out)/L;
% The whole length of the reactor is divided into n small lengths
dx = L/n;
x_calc = 0:dx:L;
nsp = nSpecies(gas_calc);
% Initialize arrays for T, Y, and rho at each location:
T_calc = zeros(length(x_calc),1);
Y_calc = zeros(length(x_calc),nsp);
rho_calc = zeros(length(x_calc),1);
T_calc(1) = temperature(gas_calc);
Y_calc(1,:) = massFractions(gas_calc);
rho_calc(1) = density(gas_calc);
for i = 2:length(x_calc)
%Solver location indicator
fprintf('Solving reactor %d of %d\n', i, length(x_calc))
%--------------------------------------------------------------------------
%------The values of variables at previous location are given as initial---
%------values to the current iteration and the limits of the current-------
%--------------reactor and the gas entering it are being set---------------
%--------------------------------------------------------------------------
inlet_soln(1) = rho_calc(i-1);
inlet_soln(2) = T_calc(i-1);
inlet_soln(3:nsp+2) = Y_calc(i-1,:);
limits = [x_calc(i-1),x_calc(i)];
set(gas_calc,'T',T_calc(i-1),'Density',rho_calc(i-1),'MoleFractions',Y_calc(i-1,:));
options = odeset('RelTol',1.e-10,'AbsTol',1e-10,'InitialStep',1e-8,'NonNegative',1);
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% These values are passed onto the ode15s solver
[~,y] = ode15s(@PFR_solver,limits,inlet_soln,options,gas_calc,mdot_calc,A_in,dAdx,k);
T_calc(i) = y(end,2);
rho_calc(i) = y(end,1);
Y_calc(i,:) = y(end,3:nsp+2);
end
A_calc = A_in+k.*x_calc*dAdx;
vx_calc = zeros(length(x_calc),1);
R_calc = zeros(length(x_calc),1);
M_calc = zeros(length(x_calc),1);
P_calc = zeros(length(x_calc),1);
for i=1:length(x_calc)
% The gas is set to the solved property values at each location
set(gas_calc,'Temperature',T_calc(i),'Density',rho_calc(i),'MassFractions',Y_calc(i,:));
% Velocity is calculated from Mass flow rate, Area and Density
vx_calc(i) = mdot_calc./(A_calc(i)*rho_calc(i));
% Specific Gas Constant
R_calc(i) = gasconstant()/meanMolecularWeight(gas_calc);
% Mach No. is calculated from local velocity and local speed of sound
M_calc(i) = vx_calc(i)/soundspeed(gas_calc);
% Pressure is calculated from density, temeprature and gas constant
P_calc(i) = rho_calc(i)*R_calc(i)*T_calc(i);
end
%% Plotting
plot(x_calc,M_calc)
xlabel('X-Position (m)')
ylabel('Mach No')
title('Mach No. Variation')
figure(2)
plot(x_calc,A_calc)
xlabel('X-Position (m)')
ylabel('Area (m^2)')
title('Reactor Profile')
figure(3)
plot(x_calc,T_calc)
xlabel('X-Position (m)')
ylabel('Temperature')
title('Temperature Variation')
figure(4)
plot(x_calc,rho_calc)
xlabel('X-Position (m)')
ylabel('Density (kg/m^3)')
title('Density Variation')
figure(5)
plot(x_calc,P_calc)
xlabel('X-Position (m)')
ylabel('Pressure (Pa)')
title('Pressure Variation')

219
samples/matlab/catcomb.m
View File

@ -1,219 +0,0 @@
% 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 Deutschmann et al., 26th
% Symp. (Intl.) on Combustion,1996 pp. 1747-1754
%
% Keywords: combustion, catalysis, 1D flow, surface chemistry
help catcomb;
clear all;
cleanup;
t0 = cputime; % record the starting time
% Parameter values are collected here to make it easier to modify
% them
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
% We will solve first for a hydrogen/air case to
% 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
%%%%%%%%%%%%%%% end of parameter list %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% 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);
set(gas,'T',tinlet,'P',p,'X',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);
setTemperature(surf_phase, tsurf);
% integrate the coverage equations in time for 1 s, holding the gas
% composition fixed to generate a good starting estimate for the
% coverages.
advanceCoverages(surf_phase, 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
set(flow, 'P', p, 'grid', initial_grid, 'tol', tol_ss, 'tol-time', 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)
set(inlt, 'T', tinlet, 'MassFlux', mdot, 'X', 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);
setTemperature(surf,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.
setProfile(sim1D, 2, {'velocity', 'spread_rate', 'T'}, ...
[0.0 1.0 % z/zmax
0.06 0.0 % velocity (u)
0.0 0.0 % spread rate (V)
tinlet tsurf]); % T
names = speciesNames(gas);
for k = 1:nSpecies(gas)
y = massFraction(inlt, k);
setProfile(sim1D, 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
enableEnergy(flow);
% disable the surface coverage equations, and turn off all gas and
% surface chemistry
setCoverageEqs(surf, 'off');
setMultiplier(surf_phase, 0.0);
setMultiplier(gas, 0.0);
% solve the problem, refining the grid if needed
solve(sim1D, 1, refine_grid);
% now turn on the surface coverage equations, and turn the
% chemistry on slowly
setCoverageEqs(surf, 'on');
for iter=1:6
mult = 10.0^(iter - 6);
setMultiplier(surf_phase, mult);
setMultiplier(gas, mult);
solve(sim1D, 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.
set(inlt,'X',comp2);
% set more stringent grid refinement criteria
setRefineCriteria(sim1D, 2, 100.0, 0.15, 0.2);
% solve the problem for the final time
solve(sim1D, loglevel, refine_grid);
% show the solution
sim1D
% save the solution
saveSoln(sim1D,'catcomb.yaml','energy',['solution with energy equation']);
%%%%%%%%%% show statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
writeStats(sim1D);
elapsed = cputime - t0;
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
disp(e);
%%%%%%%%%% make plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;
subplot(3,3,1);
plotSolution(sim1D, 'flow', 'T');
title('Temperature [K]');
subplot(3,3,2);
plotSolution(sim1D, 'flow', 'u');
title('Axial Velocity [m/s]');
subplot(3,3,3);
plotSolution(sim1D, 'flow', 'V');
title('Radial Velocity / Radius [1/s]');
subplot(3,3,4);
plotSolution(sim1D, 'flow', 'CH4');
title('CH4 Mass Fraction');
subplot(3,3,5);
plotSolution(sim1D, 'flow', 'O2');
title('O2 Mass Fraction');
subplot(3,3,6);
plotSolution(sim1D, 'flow', 'CO');
title('CO Mass Fraction');
subplot(3,3,7);
plotSolution(sim1D, 'flow', 'CO2');
title('CO2 Mass Fraction');
subplot(3,3,8);
plotSolution(sim1D, 'flow', 'H2O');
title('H2O Mass Fraction');
subplot(3,3,9);
plotSolution(sim1D, 'flow', 'H2');
title('H2 Mass Fraction');

28
samples/matlab/conhp.m
View File

@ -1,28 +0,0 @@
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.
setMassFractions(gas, y(2:end), 'nonorm');
set(gas, 'T', y(1), 'P', pressure(gas));
nsp = nSpecies(gas);
% energy equation
wdot = netProdRates(gas);
tdot = - temperature(gas) * gasconstant * enthalpies_RT(gas)' ...
* wdot / (density(gas)*cp_mass(gas));
% set up column vector for dydt
dydt = [ tdot
zeros(nsp, 1) ];
% species equations
rrho = 1.0/density(gas);
for i = 1:nsp
dydt(i+1) = rrho*mw(i)*wdot(i);
end

28
samples/matlab/conuv.m
View File

@ -1,28 +0,0 @@
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.
setMassFractions(gas, y(2:end), 'nonorm');
set(gas, 'T', y(1), 'Rho', density(gas));
nsp = nSpecies(gas);
% energy equation
wdot = netProdRates(gas);
tdot = - temperature(gas) * gasconstant * (enthalpies_RT(gas) - ones(nsp,1))' ...
* wdot / (density(gas)*cv_mass(gas));
% set up column vector for dydt
dydt = [ tdot
zeros(nsp, 1) ];
% species equations
rrho = 1.0/density(gas);
for i = 1:nsp
dydt(i+1) = rrho*mw(i)*wdot(i);
end

129
samples/matlab/diffflame.m
View File

@ -1,129 +0,0 @@
% 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.
%
% Keywords: combustion, 1D flow, diffusion flame, plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
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 = 'mixture-averaged'; % 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.
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.
%
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.
%
% Set the oxidizer inlet.
inlet_o = Inlet('air_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', 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.
%
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.
%
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.yaml','energy',['solution with energy equation']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Show statistics of solution and elapsed time.
%
writeStats(fl);
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 = 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;
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');

62
samples/matlab/equil.m
View File

@ -1,62 +0,0 @@
function equil(g)
% EQUIL a chemical equilibrium example.
%
% This example computes the adiabatic flame temperature and
% equilibrium composition for a methane/air mixture as a function of
% equivalence ratio.
%
% Keywords: combustion, equilibrium, plotting
help equil
if nargin == 1
gas = g;
else
gas = Solution('gri30.yaml');
end
nsp = nSpecies(gas);
% find methane, nitrogen, and oxygen indices
ich4 = speciesIndex(gas,'CH4');
io2 = speciesIndex(gas,'O2');
in2 = speciesIndex(gas,'N2');
nPhis = 50;
phi = linspace(0.2, 2.70, nPhis);
tad(nPhis) = 0;
xeq(nsp,nPhis) = 0;
for i = 1:nPhis
x = zeros(nsp,1);
x(ich4,1) = phi(i);
x(io2,1) = 2.0;
x(in2,1) = 7.52;
set(gas,'Temperature',300.0,'Pressure',101325.0,'MoleFractions',x);
equilibrate(gas,'HP');
tad(i) = temperature(gas);
xeq(:,i) = moleFractions(gas);
end
% make plots
clf;
subplot(1,2,1);
plot(phi,tad);
xlabel('Equivalence Ratio');
ylabel('Temperature (K)');
title('Adiabatic Flame Temperature');
subplot(1,2,2);
semilogy(phi,xeq);
axis([phi(1) phi(50) 1.0e-14 1]);
%legend(speciesName(gas,1:nsp),1);
j = 10;
for k = 1:nsp
text(phi(j),1.5*xeq(k,j),speciesName(gas,k))
j = j + 2;
if j > 46
j = 10;
end
end
xlabel('Equivalence Ratio');
ylabel('Mole Fraction');
title('Equilibrium Composition');

86
samples/matlab/flame.m
View File

@ -1,86 +0,0 @@
function f = flame(gas, left, flow, right)
% FLAME - create a flame object.
%
% gas -- object representing the gas. This object will be used to
% compute all required thermodynamic, kinetic, and transport
% properties. The state of this object should be set
% to an estimate of the gas state emerging from the
% burner before calling StagnationFlame.
%
% left -- object representing the burner, which must be
% created using function Inlet.
%
% flow -- object representing the flow, created with
% function AxisymmetricFlow.
%
% right -- object representing the surface.
%
% Check input parameters
if nargin ~= 4
error('wrong number of input arguments.');
end
if ~isIdealGas(gas)
error('gas object must represent an ideal gas mixture.');
end
if ~isInlet(left)
error('burner object of wrong type.');
end
if ~isFlow(flow)
error('flow object of wrong type.');
end
flametype = 0;
if isSurface(right)
flametype = 1;
elseif isInlet(right)
flametype = 3;
end
% create the container object
f = Stack([left flow right]);
% set default initial profiles.
rho0 = density(gas);
% find the adiabatic flame temperature and corresponding
% equilibrium composition
equilibrate(gas, 'HP');
teq = temperature(gas);
yeq = massFractions(gas);
z1 = 0.2;
mdot0 = massFlux(left);
mdot1 = massFlux(right);
t0 = temperature(left);
if flametype == 0
t1 = teq;
mdot1 = -mdot0;
else
t1 = temperature(right);
end
setProfile(f, 2, {'velocity', 'spread_rate'}, ...
[0.0 1.0
mdot0/rho0 -mdot1/rho0
0.0 0.0]);
setProfile(f, 2, 'T', [0.0 z1 1.0; t0 2000.0 t1]);
for n = 1:nSpecies(gas)
nm = speciesName(gas,n);
if strcmp(nm,'H') | strcmp(nm,'OH') | strcmp(nm,'O') | strcmp(nm,'HO2')
yint = 1.0*yeq(n);
else
yint = yeq(n);
end
if flametype == 3
y1 = massFraction(right, n);
else
y1 = yeq(n);
end
setProfile(f, 2, nm, [0 z1 1
massFraction(left, n) yint y1]);
end
% set minimal grid refinement criteria
setRefineCriteria(f, 2, 10.0, 0.8, 0.8);

113
samples/matlab/flame1.m
View File

@ -1,113 +0,0 @@
% FLAME1 - A burner-stabilized flat flame
%
% This script simulates a burner-stablized lean hydrogen-oxygen flame
% at low pressure.
%
% Keywords: combustion, 1D flow, burner-stabilized flame, plotting
help flame1
t0 = cputime; % record the starting time
% 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', 'mixture-averaged');
% set its state to that of the unburned gas at the burner
set(gas,'T', tburner, 'P', p, 'X', comp);
%%%%%%%%%%%%%%%% create the flow object %%%%%%%%%%%%%%%%%%%%%%%
f = AxisymmetricFlow(gas,'flow');
set(f, 'P', p, 'grid', initial_grid);
set(f, 'tol', tol_ss, 'tol-time', 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');
set(burner, 'T', tburner, 'MassFlux', mdot, 'X', 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);
setMaxJacAge(fl, 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,'h2flame.yaml', 'energy')
solve(fl, 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.
%
enableEnergy(f);
setRefineCriteria(fl, 2, 200.0, 0.05, 0.1);
solve(fl, 1, 1);
saveSoln(fl,'h2flame.yaml','energy',['solution with energy equation']);
%%%%%%%%%% show statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
writeStats(fl);
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');

112
samples/matlab/flame2.m
View File

@ -1,112 +0,0 @@
% FLAME1 - An axisymmetric stagnation-point non-premixed flame
%
% This script simulates a stagnation-point ethane-air flame.
%
% Keywords: combustion, 1D flow, strained flame, diffusion flame, plotting
t0 = cputime; % record the starting time
% 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','mixture-averaged');
% 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');
set(f, 'P', p, 'grid', initial_grid);
set(f, 'tol', tol_ss, 'tol-time', tol_ts);
%%%%%%%%%%%%%%% create the air inlet %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The temperature, mass flux, and composition (relative molar) may be
% specified.
%
inlet_o = Inlet('air_inlet');
set(inlet_o, 'T', tin, 'MassFlux', mdot_o, 'X', comp1);
%%%%%%%%%%%%%% create the fuel inlet %%%%%%%%%%%%%%%%%%%%%%%%%%%%
inlet_f = Inlet('fuel_inlet');
set(inlet_f, 'T', tin, 'MassFlux', mdot_f, 'X', 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,'c2h6-flame.yaml', 'energy')
% solve with fixed temperature profile first
solve(fl, 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.
%
enableEnergy(f);
setRefineCriteria(fl, 2, 200.0, 0.1, 0.1);
solve(fl, loglevel, refine_grid);
saveSoln(fl,'c2h6-flame.yaml','energy',['solution with energy equation']);
%%%%%%%%%% show statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
writeStats(fl);
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]');

102
samples/matlab/ignite.m
View File

@ -1,102 +0,0 @@
function plotdata = ignite(g)
% IGNITE Zero-dimensional kinetics: adiabatic, constant pressure.
%
% This example solves the same problem as 'reactor1,' but does
% it using one of MATLAB's ODE integrators, rather than using the
% Cantera Reactor class.
%
% Keywords: combustion, reactor network, ignition delay, plotting
help ignite
if nargin == 1
gas = g;
else
gas = Solution('gri30.yaml');
end
% 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;
[~,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');

38
samples/matlab/ignite_hp.m
View File

@ -1,38 +0,0 @@
function ignite_hp(gas)
% IGNITE_HP Solves the same ignition problem as 'ignite', but uses function
% conhp instead of reactor.
%
% Keywords: combustion, user-defined model, ignition delay, plotting
help ignite_hp
if nargin == 0
gas = Solution('h2o2.yaml', 'ohmech', 'none');
end
mw = molecularWeights(gas);
set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
y0 = [temperature(gas)
massFractions(gas)];
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 = speciesIndex(gas,'OH');
plot(out.x,out.y(1+ioh,:));
xlabel('time');
ylabel('Mass Fraction');
title('OH Mass Fraction');
end

38
samples/matlab/ignite_uv.m
View File

@ -1,38 +0,0 @@
function ignite_uv(gas)
% IGNITE_UV Solves the same ignition problem as 'ignite2', except
% that function conuv is used instead of reactor.
%
% Keywords: combustion, user-defined model, ignition delay, plotting
help ignite_uv
if nargin == 0
gas = Solution('h2o2.yaml', 'ohmech', 'none');
end
mw = molecularWeights(gas);
set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
y0 = [temperature(gas)
massFractions(gas)];
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 = speciesIndex(gas,'OH');
plot(out.x,out.y(1+ioh,:));
xlabel('time');
ylabel('Mass Fraction');
title('OH Mass Fraction');
end

59
samples/matlab/isentropic.m
View File

@ -1,59 +0,0 @@
function isentropic(g)
% ISENTROPIC isentropic, adiabatic flow example
%
% In this example, the area ratio vs. Mach number curve is
% computed for a hydrogen/nitrogen gas mixture.
%
% Keywords: thermodynamics, compressible flow, plotting
help isentropic
if nargin == 1
gas = g;
else
gas = Solution('h2o2.yaml');
end
% set the stagnation state
set(gas,'T',1200.0,'P',10.0*oneatm,'X','H2:1,N2:0.1');
s0 = entropy_mass(gas);
h0 = enthalpy_mass(gas);
p0 = pressure(gas);
mdot = 1; % arbitrary
mach = [];
a = [];
i = 1;
amin = 1.e14;
% compute values for a range of pressure ratios
for r = 0.005:0.0025:0.995
p = p0*r;
% set the state using (p,s0)
set(gas,'P',p,'S',s0);
h = enthalpy_mass(gas);
rho = density(gas);
v2 = 2.0*(h0 - h); % h + V^2/2 = h0
v = sqrt(v2);
a(i) = mdot/(rho*v); % rho*v*A = constant
if a(i) < amin
amin = a(i);
end
mach(i) = v/soundspeed(gas);
i = i + 1;
end
a = a/amin;
% plot results
clf;
plot(mach,a);
ylabel('Area Ratio');
xlabel('Mach Number');
title('Isentropic Flow: Area Ratio vs. Mach Number');

136
samples/matlab/lithium_ion_battery.m
View File

@ -1,136 +0,0 @@
% This example file calculates the cell voltage of a lithium-ion battery
% at given temperature, pressure, current, and range of state of charge (SOC).
%
% The thermodynamics are based on a graphite anode and a LiCoO2 cathode,
% modeled using the 'BinarySolutionTabulatedThermo' class.
% Further required cell parameters are the electrolyte ionic resistance, the
% stoichiometry ranges of the active materials (electrode balancing), and the
% surface area of the active materials.
%
% The functionality of this example is presented in greater detail in the
% reference (which also describes the derivation of the
% BinarySolutionTabulatedThermo class):
%
% Reference:
% M. Mayur, S. C. DeCaluwe, B. L. Kee, W. G. Bessler, “Modeling and simulation
% of the thermodynamics of lithium-ion battery intercalation materials in the
% open-source software Cantera,” Electrochim. Acta 323, 134797 (2019),
% https://doi.org/10.1016/j.electacta.2019.134797
%
% Keywords: surface chemistry, kinetics, electrochemistry, battery, plotting
% -----------------------------------------------------------------------------
% Input
% -----------------------------------------------------------------------------
% Operation parameters
SOC = 0:0.02:1; % [-] Input state of charge (0...1) (can be a vector)
I_app = -1; % [A] Externally-applied current, negative for discharge
T = 293; % [K] Temperature
P = oneatm; % [Pa] Pressure
% Cell properties
inputFile = 'lithium_ion_battery.yaml'; % Cantera input file name
R_elyt = 0.0384; % [Ohm] Electrolyte resistance
S_ca = 1.1167; % [m^2] Cathode total active material surface area
S_an = 0.7824; % [m^2] Anode total active material surface area
% Electrode balancing: The "balancing" of the electrodes relates the chemical
% composition (lithium mole fraction in the active materials) to the macroscopic
% cell-level state of charge.
X_Li_an_0 = 0.01; % [-] anode Li mole fraction at SOC = 0 %
X_Li_an_1 = 0.75; % [-] anode Li mole fraction at SOC = 100 %
X_Li_ca_0 = 0.99; % [-] cathode Li mole fraction at SOC = 0 %
X_Li_ca_1 = 0.49; % [-] cathode Li mole fraction at SOC = 100 %
% -----------------------------------------------------------------------------
% Calculations
% -----------------------------------------------------------------------------
% Calculate mole fractions from SOC
X_Li_an = (X_Li_an_1-X_Li_an_0)*SOC+X_Li_an_0; % anode balancing
X_Li_ca = (X_Li_ca_0-X_Li_ca_1)*(1-SOC)+X_Li_ca_1; % cathode balancing
% Import all Cantera phases
anode = Solution(inputFile, 'anode');
cathode = Solution(inputFile, 'cathode');
elde = Solution(inputFile, 'electron');
elyt = Solution(inputFile, 'electrolyte');
anode_interface = Interface(inputFile, 'edge_anode_electrolyte', anode, elde, elyt);
cathode_interface = Interface(inputFile, 'edge_cathode_electrolyte', cathode, elde, elyt);
% Set the temperatures and pressures of all phases
set(anode,'T',T,'P',P);
set(cathode,'T',T,'P',P);
set(elde,'T',T,'P',P);
set(elyt,'T',T,'P',P);
set(anode_interface,'T',T,'P',P);
set(cathode_interface,'T',T,'P',P);
% Calculate cell voltage, separately for each entry of the input vectors
V_cell = zeros(length(SOC),1);
phi_l_an = 0;
phi_s_ca = 0;
for i = 1:length(SOC)
% Set anode electrode potential to 0
phi_s_an = 0;
% Calculate anode electrolyte potential
phi_l_an = fzero(@(E) anode_curr(phi_s_an,E,X_Li_an(i),anode,elde,elyt,anode_interface,S_an)-I_app, phi_l_an);
% Calculate cathode electrolyte potential
phi_l_ca = phi_l_an + I_app*R_elyt;
% Calculate cathode electrode potential
phi_s_ca = fzero(@(E) cathode_curr(E,phi_l_ca,X_Li_ca(i),cathode,elde,elyt,cathode_interface,S_ca)-I_app, phi_s_ca);
% Calculate cell voltage
V_cell(i) = phi_s_ca - phi_s_an;
end
% Let's plot the cell voltage, as a function of the state of charge:
figure(1);
plot(SOC*100,V_cell,'linewidth',2.5)
ylim([2.5,4.3])
xlabel('State of charge / %')
ylabel('Cell voltage / V')
set(gca,'fontsize',14)
%--------------------------------------------------------------------------
% Helper functions
% -----------------------------------------------------------------------------
% This function returns the Cantera calculated anode current (in A)
function anCurr = anode_curr(phi_s,phi_l,X_Li_an,anode,elde,elyt,anode_interface,S_an)
% Set the active material mole fraction
set(anode,'X',['Li[anode]:' num2str(X_Li_an) ', V[anode]:' num2str(1-X_Li_an)]);
% Set the electrode and electrolyte potential
setElectricPotential(elde,phi_s);
setElectricPotential(elyt,phi_l);
% Get the net reaction rate at the anode-side interface
% Reaction according to YAML file: Li+[elyt] + V[anode] + electron <=> Li[anode]
r = rop_net(anode_interface); % [kmol/m2/s]
% Calculate the current. Should be negative for cell discharge.
anCurr = r*faradayconstant*S_an; %
end
% This function returns the Cantera calculated cathode current (in A)
function caCurr = cathode_curr(phi_s,phi_l,X_Li_ca,cathode,elde,elyt,cathode_interface,S_ca)
% Set the active material mole fractions
set(cathode,'X',['Li[cathode]:' num2str(X_Li_ca) ', V[cathode]:' num2str(1-X_Li_ca)]);
% Set the electrode and electrolyte potential
setElectricPotential(elde,phi_s);
setElectricPotential(elyt,phi_l);
% Get the net reaction rate at the cathode-side interface
% Reaction according to YAML file: Li+[elyt] + V[cathode] + electron <=> Li[cathode]
r = rop_net(cathode_interface); % [kmol/m2/s]
% Calculate the current. Should be negative for cell discharge.
caCurr = r*faradayconstant*S_ca*(-1); %
end

106
samples/matlab/periodic_cstr.m
View File

@ -1,106 +0,0 @@
function periodic_cstr
%
% A CSTR with steady inputs but periodic interior state.
%
% A stoichiometric hydrogen/oxygen mixture is introduced and reacts to produce
% water. But since water has a large efficiency as a third body in the chain
% termination reaction
%
% H + O2 + M = HO2 + M
%
% as soon as a significant amount of water is produced the reaction
% stops. After enough time has passed that the water is exhausted from
% the reactor, the mixture explodes again and the process
% repeats. This explanation can be verified by decreasing the rate for
% reaction 7 in file 'h2o2.yaml' and re-running the example.
%
% Acknowledgments: The idea for this example and an estimate of the
% conditions needed to see the oscillations came from Bob Kee,
% Colorado School of Mines
%
% Keywords: combustion, reactor network, well-stirred reactor, plotting
help periodic_cstr
% create the gas mixture
gas = Solution('h2o2.yaml');
% pressure = 60 Torr, T = 770 K
p = 60.0*133.3;
t = 770.0;
OneAtm = 1.01325e5;
set(gas,'T', 300.0, 'P', p, 'X', 'H2:2, O2:1');
% create an upstream reservoir that will supply the reactor. The
% temperature, pressure, and composition of the upstream reservoir are
% set to those of the 'gas' object at the time the reservoir is
% created.
upstream = Reservoir(gas);
% Now set the gas to the initial temperature of the reactor, and create
% the reactor object.
set(gas, 'T', t, 'P', p);
cstr = IdealGasReactor(gas);
% Set its volume to 10 cm^3. In this problem, the reactor volume is
% fixed, so the initial volume is the volume at all later times.
setInitialVolume(cstr, 10.0*1.0e-6);
% We need to have heat loss to see the oscillations. Create a
% reservoir to represent the environment, and initialize its
% temperature to the reactor temperature.
env = Reservoir(gas);
% Create a heat-conducting wall between the reactor and the
% environment. Set its area, and its overall heat transfer
% coefficient. Larger U causes the reactor to be closer to isothermal.
% If U is too small, the gas ignites, and the temperature spikes and
% stays high.
w = Wall;
install(w, cstr, env);
setArea(w, 1.0);
setHeatTransferCoeff(w, 0.02);
% Connect the upstream reservoir to the reactor with a mass flow
% controller (constant mdot). Set the mass flow rate to 1.25 sccm.
sccm = 1.25;
vdot = sccm * 1.0e-6/60.0 * ((OneAtm / pressure(gas)) * ( temperature(gas) / 273.15)); % m^3/s
mdot = density(gas) * vdot; % kg/s
mfc = MassFlowController;
install(mfc, upstream, cstr);
setMassFlowRate(mfc, mdot);
% now create a downstream reservoir to exhaust into.
downstream = Reservoir(gas);
% connect the reactor to the downstream reservoir with a valve, and
% set the coefficient sufficiently large to keep the reactor pressure
% close to the downstream pressure of 60 Torr.
v = Valve;
install(v, cstr, downstream);
setValveCoeff(v, 1.0e-9);
% create the network
network = ReactorNet({cstr});
% now integrate in time
tme = 0.0;
dt = 0.1;
n = 0;
while tme < 300.0
n = n + 1;
tme = tme + dt;
advance(network, tme);
tm(n) = tme;
y(1,n) = massFraction(cstr,'H2');
y(2,n) = massFraction(cstr,'O2');
y(3,n) = massFraction(cstr,'H2O');
end
clf
figure(1)
plot(tm,y)
legend('H2','O2','H2O')
title('Mass Fractions')

76
samples/matlab/prandtl1.m
View File

@ -1,76 +0,0 @@
function prandtl1(g)
% PRANDTL1 Prandtl number for an equilibrium H/O gas mixture.
%
% This example computes and plots the Prandtl number for a
% hydrogen / oxygen mixture in chemical equilibrium for P = 1
% atm and a range of temperatures and elemental O/(O+H) ratios.
%
% Keywords: equilibrium, transport, plotting
help prandtl1
if nargin == 1
gas = g;
else
gas = Solution('h2o2.yaml', 'ohmech', 'mixture-averaged');
end
pr = zeros(31,31);
xh2 = zeros(31,31);
visc = zeros(31,31);
lambda = zeros(31,31);
t = [];
xo2 = [];
io2 = speciesIndex(gas,'O2');
ih2 = speciesIndex(gas,'H2');
minT = minTemp(gas);
maxT = maxTemp(gas);
dT = (maxT - minT)/30.0;
t0 = cputime;
for i = 1:31
t(i) = minT + dT*(i-1);
for j = 1:31
xo2(j) = 0.99*(j-1)/30.0;
x = zeros(nSpecies(gas),1);
x(io2) = xo2(j);
x(ih2) = 1.0 - xo2(j);
set(gas,'T',t(i),'P',oneatm,'X',x);
equilibrate(gas,'TP');
visc(i,j) = viscosity(gas);
lambda(i,j) = thermalConductivity(gas);
pr(i,j) = visc(i,j)*cp_mass(gas)/lambda(i,j);
x = moleFractions(gas);
xh2(i,j) = x(ih2);
end
end
disp(['CPU time = ' num2str(cputime - t0)]);
% plot results
clf;
%figure(1);
subplot(2,2,1);
surf(xo2,t,pr);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Prandtl Number');
subplot(2,2,2);
surf(xo2,t,xh2);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('H_2 Mole Fraction');
subplot(2,2,3);
surf(xo2,t,visc);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Viscosity');
subplot(2,2,4);
surf(xo2,t,lambda);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Thermal Conductivity');

74
samples/matlab/prandtl2.m
View File

@ -1,74 +0,0 @@
function prandtl2(g)
% PRANDTL2 Prandtl number for an equilibrium H/O gas mixture.
%
% This example does the same thing as prandtl1, but using
% the multicomponent expression for the thermal conductivity.
%
% Keywords: transport, equilibrium, multicomponent transport, plotting
help prandtl2
if nargin == 1
gas = g;
else
gas = Solution('h2o2.yaml', 'ohmech', 'multicomponent');
end
pr = zeros(31,31);
xh2 = zeros(31,31);
visc = zeros(31,31);
lambda = zeros(31,31);
t = [];
xo2 = [];
io2 = speciesIndex(gas,'O2');
ih2 = speciesIndex(gas,'H2');
minT = minTemp(gas);
maxT = maxTemp(gas);
dT = (maxT - minT)/30.0;
t0 = cputime;
for i = 1:31
t(i) = minT + dT*(i-1);
for j = 1:31
xo2(j) = 0.99*(j-1)/30.0;
x = zeros(nSpecies(gas),1);
x(io2) = xo2(j);
x(ih2) = 1.0 - xo2(j);
set(gas,'T',t(i),'P',oneatm,'X',x);
equilibrate(gas,'TP');
visc(i,j) = viscosity(gas);
lambda(i,j) = thermalConductivity(gas);
pr(i,j) = visc(i,j)*cp_mass(gas)/lambda(i,j);
x = moleFractions(gas);
xh2(i,j) = x(ih2);
end
end
disp(['CPU time = ' num2str(cputime - t0)]);
% plot results
clf;
subplot(2,2,1);
surf(xo2,t,pr);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Prandtl Number');
subplot(2,2,2);
surf(xo2,t,xh2);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('H_2 Mole Fraction');
subplot(2,2,3);
surf(xo2,t,visc);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Viscosity');
subplot(2,2,4);
surf(xo2,t,lambda);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Thermal Conductivity');

61
samples/matlab/rankine.m
View File

@ -1,61 +0,0 @@
function [work, efficiency] = rankine(t1, p2, eta_pump, eta_turbine)
%
% This example computes the efficiency of a simple vapor power cycle.
%
% Keywords: thermodynamics, thermodynamic cycle, non-ideal fluid
help rankine
% create an object representing water
w = Water;
% start with saturated liquid water at t1
set(w,'T',t1,'Liquid',1.0);
p1 = pressure(w);
% pump it to p2
pump_work = pump(w, p2, eta_pump);
h2 = enthalpy_mass(w);
p2 = pressure(w);
% heat to saturated vapor
set(w,'P',p2,'Vapor',1.0);
h3 = enthalpy_mass(w);
heat_added = h3 - h2;
% expand adiabatically back to the initial pressure
work = expand(w, p1, eta_turbine);
% compute the efficiency
efficiency = (work - pump_work)/heat_added;
function w = pump(fluid, pfinal, eta)
% PUMP - Adiabatically pump a fluid to pressure pfinal, using a pump
% with isentropic efficiency eta.
%
h0 = enthalpy_mass(fluid);
s0 = entropy_mass(fluid);
set(fluid, 'S', s0, 'P', pfinal);
h1s = enthalpy_mass(fluid);
isentropic_work = h1s - h0;
actual_work = isentropic_work / eta;
h1 = h0 + actual_work;
set(fluid, 'H',h1, 'P',pfinal);
w = actual_work;
function w = expand(fluid, pfinal, eta)
% EXPAND - Adiabatically expand a fluid to pressure pfinal, using a
% turbine with isentropic efficiency eta.
%
h0 = enthalpy_mass(fluid);
s0 = entropy_mass(fluid);
set(fluid, 'S', s0, 'P', pfinal);
h1s = enthalpy_mass(fluid);
isentropic_work = h0 - h1s;
actual_work = isentropic_work * eta;
h1 = h0 - actual_work;
set(fluid, 'H',h1, 'P',pfinal);
w = actual_work;

83
samples/matlab/reactor1.m
View File

@ -1,83 +0,0 @@
function reactor1(g)
% REACTOR1 Zero-dimensional kinetics: adiabatic, constant pressure.
%
% >>>> For a simpler way to carry out a constant-pressure simulation,
% see example reactor3.m <<<<<
%
% This example illustrates how to use class 'Reactor' for
% zero-dimensional kinetics simulations. Here the parameters are
% set so that the reactor is adiabatic and very close to constant
% pressure.
%
% Keywords: combustion, reactor network, ignition delay, plotting
help reactor1
if nargin == 1
gas = g;
else
gas = Solution('h2o2.yaml', 'ohmech', 'none');
end
P = oneatm;
% set the initial conditions
set(gas,'T',1001.0,'P',P,'X','H2:2,O2:1,N2:4');
% create a reactor, and insert the gas
r = IdealGasReactor(gas);
% create a reservoir to represent the environment
a = Solution('air.yaml','air','none');
set(a,'P',P)
env = Reservoir(a);
% Define a wall between the reactor and the environment and
% make it flexible, so that the pressure in the reactor is held
% at the environment pressure.
w = Wall;
install(w,r,env);
% set expansion parameter. dV/dt = KA(P_1 - P_2)
setExpansionRateCoeff(w, 1.0e6);
% set wall area
setArea(w, 1.0);
% create a reactor network and insert the reactor:
network = ReactorNet({r});
nSteps = 100;
tim(nSteps) = 0;
temp(nSteps) = 0;
x(nSteps,3) = 0;
t = 0.0;
dt = 1.0e-5;
t0 = cputime;
for n = 1:nSteps
t = t + dt;
advance(network, t);
tim(n) = time(network);
temp(n) = temperature(r);
x(n,1:3) = moleFraction(gas,{'OH','H','H2'});
end
disp(['CPU time = ' num2str(cputime - t0)]);
clf;
subplot(2,2,1);
plot(tim,temp);
xlabel('Time (s)');
ylabel('Temperature (K)');
subplot(2,2,2)
plot(tim,x(:,1));
xlabel('Time (s)');
ylabel('OH Mole Fraction (K)');
subplot(2,2,3)
plot(tim,x(:,2));
xlabel('Time (s)');
ylabel('H Mole Fraction (K)');
subplot(2,2,4)
plot(tim,x(:,3));
xlabel('Time (s)');
ylabel('H2 Mole Fraction (K)');
clear all
cleanup

61
samples/matlab/reactor2.m
View File

@ -1,61 +0,0 @@
function reactor2(g)
% REACTOR2 Zero-dimensional kinetics: adiabatic, constant volume.
%
% This example illustrates how to use class 'Reactor' for
% zero-dimensional kinetics simulations. Here the parameters are
% set so that the reactor is adiabatic and constant volume.
%
% Keywords: combustion, reactor network, ignition delay, plotting
help reactor2
if nargin == 1
gas = g;
else
gas = Solution('h2o2.yaml', 'ohmech', 'none');
end
% set the initial conditions
set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
% create a reactor, and insert the gas
r = IdealGasReactor(gas);
% create a reactor network and insert the reactor
network = ReactorNet({r});
nSteps = 100;
tim(nSteps) = 0;
temp(nSteps) = 0;
x(nSteps,3) = 0;
t = 0;
dt = 1.0e-5;
t0 = cputime;
for n = 1:100
t = t + dt;
advance(network, t);
tim(n) = time(network);
temp(n) = temperature(r);
x(n,1:3) = moleFraction(gas,{'OH','H','H2'});
end
disp(['CPU time = ' num2str(cputime - t0)]);
clf;
subplot(2,2,1);
plot(tim,temp);
xlabel('Time (s)');
ylabel('Temperature (K)');
subplot(2,2,2)
plot(tim,x(:,1));
xlabel('Time (s)');
ylabel('OH Mole Fraction (K)');
subplot(2,2,3)
plot(tim,x(:,2));
xlabel('Time (s)');
ylabel('H Mole Fraction (K)');
subplot(2,2,4)
plot(tim,x(:,3));
xlabel('Time (s)');
ylabel('H2 Mole Fraction (K)');
clear all
cleanup

50
samples/matlab/reactor_ode.m
View File

@ -1,50 +0,0 @@
function dydt = reactor_ode(t,y,gas,vdot,area,heatflux)
% REACTOR ODE system for a generic zero-dimensional reactor.
%
% Function REACTOR evaluates the system of ordinary differential
% equations for a zero-dimensional reactor with arbitrary heat
% transfer and volume change.
%
% Solution vector components:
% y(1) Total internal energy U
% y(2) Volume V
% y(3) Mass of species 1
% ....
% y(2+nsp) Mass of last species
%
[m,n] = size(y);
dydt = zeros(m,n);
for j = 1:n
this_y = y(:,j);
int_energy = this_y(1);
vol = this_y(2);
masses = this_y(3:end);
% evaluate the total mass, and the specific internal energy and volume.
total_mass = sum(masses);
u_mass = int_energy/total_mass;
v_mass = vol/total_mass;
% set the state of the gas by specifying (u,v,{Y_k})
setMassFractions(gas,masses);
setState_UV(gas, [u_mass v_mass]);
p = pressure(gas);
% volume equation
vdt = feval(vdot, t, vol, gas);
% energy equation
a = feval(area, t, vol);
q = feval(heatflux, t, gas);
udt = -p * vdt + a * q;
% species equations
ydt = total_mass * ydot(gas);
% set up column vector for dydt
dydt(:,j) = [udt
vdt
ydt ];
end

99
samples/matlab/surfreactor.m
View File

@ -1,99 +0,0 @@
% SURFREACTOR Zero-dimensional reactor with surface chemistry
%
% This example illustrates how to use class 'Reactor' for
% zero-dimensional simulations including both homogeneous and
% heterogeneous chemistry.
%
% Keywords: catalysis, combustion, reactor network, plotting
help surfreactor
t = 870.0;
gas = Solution('ptcombust.yaml','gas');
% set the initial conditions
set(gas,'T',t,'P',oneatm,'X','CH4:0.01, O2:0.21, N2:0.78');
% The surface reaction mechanism describes catalytic combustion of
% methane on platinum, and is from Deutschmann et al., 26th
% Symp. (Intl.) on Combustion,1996, pp. 1747-1754
surf = importInterface('ptcombust.yaml','Pt_surf', gas);
setTemperature(surf, t);
nsp = nSpecies(gas);
nSurfSp = nSpecies(surf);
% create a reactor, and insert the gas
r = IdealGasReactor(gas);
setInitialVolume(r, 1.0e-6)
% create a reservoir to represent the environment
a = Solution('air.yaml','air','none');
set(a,'T',t,'P',oneatm);
env = Reservoir(a);
% Define a wall between the reactor and the environment and
% make it flexible, so that the pressure in the reactor is held
% at the environment pressure.
w = Wall;
install(w,r,env);
A = 1e-4; % Wall area
% Add a reacting surface, with an area matching that of the wall
rsurf = ReactorSurface(surf, r, A);
% set the wall area and heat transfer coefficient.
setArea(w, A);
setHeatTransferCoeff(w,1.0e1); % W/m2/K
% set expansion rate parameter. dV/dt = KA(P_1 - P_2)
setExpansionRateCoeff(w, 1.0);
network = ReactorNet({r});
% setTolerances(network, 1.0e-8, 1.0e-12);
nSteps = 100;
p0 = pressure(r);
names = {'CH4','CO','CO2','H2O'};
x = zeros([nSteps 4]);
tim = zeros(nSteps, 1);
temp = zeros(nSteps, 1);
pres = zeros(nSteps, 1);
cov = zeros([nSteps nSurfSp]);
t = 0;
dt = 0.1;
t0 = cputime;
for n = 1:nSteps
t = t + dt;
advance(network, t);
tim(n) = t;
temp(n) = temperature(r);
pres(n) = pressure(r) - p0;
cov(n,:) = coverages(surf)';
x(n,:) = moleFraction(gas,names);
end
disp(['CPU time = ' num2str(cputime - t0)]);
clf;
subplot(2,2,1);
plot(tim,temp);
xlabel('Time (s)');
ylabel('Temperature (K)');
subplot(2,2,2);
plot(tim,pres);
axis([0 5 -0.1 0.1]);
xlabel('Time (s)');
ylabel('Delta Pressure (Pa)');
subplot(2,2,3);
semilogy(tim,cov);
xlabel('Time (s)');
ylabel('Coverages');
legend(speciesNames(surf));
subplot(2,2,4);
plot(tim,x);
xlabel('Time (s)');
ylabel('Mole Fractions');
legend(names);
clear all
cleanup

29
samples/matlab/test_examples.m
View File

@ -1,29 +0,0 @@
% runs selected examples without pausing
clear all
close all
cleanup
equil();
isentropic();
reactor1();
reactor2();
surfreactor;
periodic_cstr;
Plug_Flow_Reactor;
lithium_ion_battery
rankine(300.0, 2.0*oneatm, 0.8, 0.7);
prandtl1();
prandtl2();
flame1;
flame2;
catcomb;
clear all
close all
diffflame;
ignite_hp;
ignite_uv;
clear all
close all
cleanup

174
samples/matlab/tut1.m
View File

@ -1,174 +0,0 @@
% Tutorial 1: Getting started
%
% Topics:
% - creating a gas mixture
% - setting the state
% - cleaning up
%
% Keywords: tutorial
help tut1
% Start MATLAB, and at the prompt type:
gas1 = GRI30
% If you have successfully installed the Cantera toolbox, you should
%see something like this:
% temperature 300 K
% pressure 101325 Pa
% density 0.081889 kg/m^3
% mean mol. weight 2.01588 amu
%
% 1 kg 1 kmol
% ----------- ------------
% enthalpy 26470.1 5.336e+04 J
% internal energy -1.21088e+06 -2.441e+06 J
% entropy 64914 1.309e+05 J/K
% Gibbs function -1.94477e+07 -3.92e+07 J
% heat capacity c_p 14311.8 2.885e+04 J/K
% heat capacity c_v 10187.3 2.054e+04 J/K
%
% X Y Chem. Pot. / RT
% ------------- ------------ ------------
% H2 1 1 -15.7173
% [ +52 minor] 0 0
%
% What you have just done is to create an object ("gas1") that
% implements GRI-Mech 3.0, the 53-species, 325-reaction natural gas
% combustion mechanism developed by Gregory P. Smith, David M. Golden,
% Michael Frenklach, Nigel W. Moriarty, Boris Eiteneer, Mikhail
% Goldenberg, C. Thomas Bowman, Ronald K. Hanson, Soonho Song, William
% C. Gardiner, Jr., Vitali V. Lissianski, and Zhiwei Qin. (See
% http://www.me.berkeley.edu/gri_mech/ for more information about
% GRI-Mech 3.0.)
%
% The object created by GI30 has properties you would expect for a gas
% mixture - it has a temperature, a pressure, species mole and mass
% fractions, etc. As we'll soon see, it has many other properties too.
%
% The summary of the state of 'gas1' printed above shows that new
% objects created by function GRI30 start out with a temperature of
% 300 K, a pressure of 1 atm, and have a composition that consists of
% only one species, in this case hydrogen. There is nothing special
% about H2 - it just happens to be the first species listed in the
% input file defining GRI-Mech 3.0 that the 'GRI30' function reads. In
% general, the species listed first will initially have a mole
% fraction of 1.0, and all of the others will be zero.
% Setting the state
% -----------------
% The state of the object can be easily changed. For example,
setTemperature(gas1, 1200)
% sets the temperature to 1200 K. (Cantera always uses SI units.)
% Notice in the summary of properties that MATLAB prints after this
% command is executed that the temperature has been changed as
% requested, but the pressure has changed too. The density and
% composition have not.
%
% When setting properties individually, some convention needs to be
% adopted to specify which other properties are held constant. This is
% because thermodynamics requires that *two* properties (not one) in
% addition to composition information be specified to fix the
% intensive state of a substance (or mixture).
%
% Cantera adopts the following convention: only one of the set
% (temperature, density, mass fractions) is altered by setting any
% single property. In particular:
%
% a) Setting the temperature is done holding density and
% composition fixed. (The pressure changes.)
% b) Setting the pressure is done holding temperature and
% composition fixed. (The density changes.)
%
% c) Setting the composition is done holding temperature
% and density fixed. (The pressure changes).
%
% Setting multiple properties: the 'set' method
% ---------------------------------------------
% If you want to set multiple properties at once, use the 'set'
% method. (Note: a 'method' is just the term for a function that acts
% on an object. In MATLAB, methods take the object as the first
% argument.)
set(gas1, 'Temperature', 900.0, 'Pressure', 1.e5);
% This statement sets both temperature and pressure at the same
% time. Any number of property/value pairs can be specified in a
% call to 'set'. For example, the following sets the mole fractions
% too:
set(gas1, 'Temperature', 900.0, 'Pressure', 1.e5, 'MoleFractions', ...
'CH4:1,O2:2,N2:7.52');
% The 'set' method also accepts abbreviated property names:
set(gas1,'T',900.0,'P',1.e5,'X','CH4:1,O2:2,N2:7.52')
% Either version results in
%
% temperature 900 K
% pressure 100000 Pa
% density 0.369279 kg/m^3
% mean mol. weight 27.6332 amu
%
% 1 kg 1 kmol
% ----------- ------------
% enthalpy 455660 1.259e+07 J
% internal energy 184862 5.108e+06 J
% entropy 8529.31 2.357e+05 J/K
% Gibbs function -7.22072e+06 -1.995e+08 J
% heat capacity c_p 1304.4 3.604e+04 J/K
% heat capacity c_v 1003.52 2.773e+04 J/K
%
% X Y Chem. Pot. / RT
% ------------- ------------ ------------
% O2 0.190114 0.220149 -27.9596
% CH4 0.095057 0.0551863 -37.0813
% N2 0.714829 0.724665 -24.935
% [ +50 minor] 0 0
% Other properties may also be set using 'set', including some that
% can't be set individually. The following property pairs may be
% set: (Enthalpy, Pressure), (IntEnergy, Volume), (Entropy,
% Volume), (Entropy, Pressure). In each case, the values of the
% extensive properties must be entered *per unit mass*.
% Setting the enthalpy and pressure:
set(gas1, 'Enthalpy', 2*enthalpy_mass(gas1), 'Pressure', 2*oneatm);
% The composition above was specified using a string. The format is a
% comma-separated list of <species name>:<relative mole numbers>
% pairs. The mole numbers will be normalized to produce the mole
% fractions, and therefore they are 'relative' mole numbers. Mass
% fractions can be set in this way too by changing 'X' to 'Y' in the
% above statement.
% The composition can also be set using an array, which can be
% either a column vector or a row vector but must have the same
% size as the number of species. For example, to set all 53 mole
% fractions to the same value, do this:
x = ones(53,1); % a column vector of 53 ones
set(gas1, 'X', x)
% To set the mass fractions to equal values:
set(gas1, 'Y', x)
% This clears all Matlab objects created
clear all
% and this clears all Cantera objects created
cleanup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of tutorial 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

105
samples/matlab/tut2.m
View File

@ -1,105 +0,0 @@
% Tutorial 2: Working with input files
%
% Topics:
% - using functions 'Solution' and 'importInterface'
% - input files distributed with Cantera
% - the Cantera search path
% - CTML files
% - converting from CK format
%
% Keywords: tutorial, input files
help tut2
t0 = cputime;
% In the last tutorial, we used function GRI30 to create an object
% that models an ideal gas mixture with the species and reactions of
% GRI-Mech 3.0. Another way to do this is shown here, with statements
% added to measure how long this takes:
gas1 = Solution('gri30.yaml', 'gas', 'mixture-averaged');
msg = sprintf('time to create gas1: %f', cputime - t0)
% Function 'Solution' constructs an object representing a phase of
% matter by reading in attributes of the phase from a file, which in
% this case is 'gri30.yaml'. This file contains a complete specification
% of the GRI-Mech 3.0 reaction mechanism, including element data
% (name, atomic weight), species data (name, elemental composition,
% coefficients to compute thermodynamic and transport properties), and
% reaction data (stoichiometry, rate coefficient parameters). The file
% is written in a format understood by Cantera, which is described in
% the document "Defining Phases and Interfaces."
% On some systems, processing long files can be a little slow. Loading
% a file again takes much less time, because Cantera 'remembers' files
% it has already processed and doesn't need to read them in again:
t0 = cputime;
gas1b = Solution('gri30.yaml', 'gas', 'none');
msg = sprintf('time to create gas1b: %f', cputime - t0)
% Since GRI-Mech is a rather small mechanism, there might not be much
% difference in these times.
% YAML files distributed with Cantera
%------------------------------------
% Several reaction mechanism files in this format are included in the
% Cantera distribution, including ones that model high-temperature
% air, a hydrogen/oxygen reaction mechanism, and a few surface
% reaction mechanisms. Under Windows, these files may be located in
% 'C:\Program Files\Common Files\Cantera', or in 'C:\cantera\data',
% depending on how you installed Cantera and the options you
% specified. On a Unix/Linux/macOS machine, they are usually kept
% in the 'data' subdirectory within the Cantera installation
% directory.
% If for some reason Cantera has difficulty finding where these files
% are on your system, set environment variable CANTERA_DATA to the
% directory where they are located. Alternatively, you can call function
% adddir to add a directory to the Cantera search path:
adddir('/usr/local/cantera/my_data_files');
% Cantera input files are plain text files, and can be created with
% any text editor. See the document 'Defining Phases and Interfaces'
% for more information.
% Importing multiple phases or interfaces
% ---------------------------------------
% A Cantera input file may contain more than one phase specification,
% or may contain specifications of interfaces (surfaces). Here we
% import definitions of two bulk phases and the interface between them
% from file diamond.yaml:
gas2 = Solution('diamond.yaml', 'gas'); % a gas
diamond = Solution('diamond.yaml','diamond'); % bulk diamond
diamonnd_surf = importInterface('diamond.yaml','diamond_100',...
gas2, diamond);
% Note that the bulk (that is, 3D) phases that participate in the surface
% reactions must also be passed as arguments to importInterface.
% Converting CK-format files
% --------------------------
% Many existing reaction mechanism files are in "CK format," by which
% we mean the input file format developed for use with the Chemkin-II
% software package. [See R. J. Kee, F. M. Rupley, and J. A. Miller,
% Sandia National Laboratories Report SAND89-8009 (1989).]
% Cantera comes with a converter utility program 'ck2yaml' (or
% 'ck2yaml.exe') that converts CK format into YAML format. This
% program should be run from the command line first to convert any CK
% files you plan to use into YAML format.
%
% Here's an example of how to use it:
%
% ck2yaml --input=mech.inp --thermo=therm.dat --transport=tran.dat --name=mymech
clear all
cleanup

55
samples/matlab/tut3.m
View File

@ -1,55 +0,0 @@
% Tutorial 3: Getting Help
%
% Keywords: tutorial
help tut3
% Suppose you have created a Cantera object and want to know what
% methods are available for it, and get help on using the methods.
g = GRI30
% The first thing you need to know is the MATLAB class object g
% belongs to. Type:
class(g)
% This tells you that g belongs to a class called 'Solution'. To find
% the methods for this class, type
methods Solution
% This command returns only a few method names. These are the ones
% directly defined in this class. But solution inherits many other
% methods from base classes. To see all of its methods, type
methods Solution -full
% Now a long list is printed, along with a specification of the class
% the method is inherited from. For example, 'setPressure' is
% inherited from a class 'ThermoPhase'. Don't be concerned at this
% point about what these base classes are - we'll come back to them
% later.
% Now that you see what methods are available, you can type 'help
% <method_name>' to print help text for any method. For example,
help setTemperature
help setMassFractions
help rop_net
% For help on how to construct objects of a given class, type 'help
% <classname>'
help Solution
% Now that you know how to get help when you need it, you can
% explore using the Cantera Toolbox on your own. But there are a
% few more useful things to know, which are described in the next
% few tutorials.
clear all
cleanup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of tutorial 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

87
samples/matlab/tut4.m
View File

@ -1,87 +0,0 @@
% Tutorial 4: Chemical Equilibrium
%
% Topics:
% - the equilibrate method
% - specifying fixed TP, HP, UV, SV, or SP
% - checking reaction rates of progress
%
% Keywords: tutorial, equilibrium, kinetics
help tut4
% To set a gas mixture to a state of chemical equilibrium, use the
% 'equilibrate' method.
%
g = GRI30('none');
set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52')
equilibrate(g,'TP')
% The above statement sets the state of object 'g' to the state of
% chemical equilibrium holding temperature and pressure
% fixed. Alternatively, the specific enthalpy and pressure can be held
% fixed:
disp('fixed H and P:');
set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2.0,N2:7.52');
equilibrate(g,'HP')
% Other options are
% 'UV' fixed specific internal energy and specific volume
% 'SV' fixed specific entropy and specific volume
% 'SP' fixed specific entropy and pressure
disp('fixed U and V:');
set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52');
equilibrate(g,'UV')
disp('fixed S and V:');
set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52');
equilibrate(g,'SV')
disp('fixed S and P:');
set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52');
equilibrate(g,'SP')
% How can you tell if 'equilibrate' has correctly found the
% chemical equilibrium state? One way is verify that the net rates of
% progress of all reversible reactions are zero.
% Here is the code to do this:
set(g,'T',2000.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52');
equilibrate(g,'TP')
rf = rop_f(g);
rr = rop_r(g);
format short e;
for i = 1:nReactions(g)
if isReversible(g,i)
disp([i, rf(i), rr(i), (rf(i) - rr(i))/rf(i)]);
end
end
% You might be wondering how 'equilibrate' works. (Then again, you might
% not, in which case you can go on to the next tutorial now.) Method
% 'equilibrate' invokes Cantera's chemical equilibrium solver, which
% uses an element potential method. The element potential method is
% one of a class of equivalent 'nonstoichiometric' methods that all
% have the characteristic that the problem reduces to solving a set of
% M nonlinear algebraic equations, where M is the number of elements
% (not species). The so-called 'stoichiometric' methods, on the other
% hand, (including Gibbs minimization), require solving K nonlinear
% equations, where K is the number of species (usually K >> M). See
% Smith and Missen, "Chemical Reaction Equilibrium Analysis" for more
% information on the various algorithms and their characteristics.
%
% Cantera uses a damped Newton method to solve these equations, and
% does a few other things to generate a good starting guess and to
% produce a reasonably robust algorithm. If you want to know more
% about the details, look at the on-line documented source code of
% Cantera C++ class 'ChemEquil' at https://cantera.org.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
cleanup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of tutorial 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

109
samples/matlab/tut5.m
View File

@ -1,109 +0,0 @@
% Tutorial 5: Reaction information and rates
%
% Topics:
% - stoichiometric coefficients
% - reaction rates of progress
% - species production rates
% - reaction equations
% - equilibrium constants
% - rate multipliers
%
% Keywords: tutorial, kinetics
help tut5
g = GRI30('none');
set(g,'T',1500,'P',oneatm,'X',ones(nSpecies(g),1));
% Methods are provided that compute many quantities of interest for
% kinetics. Some of these are:
% 1) Stoichiometric coefficients
nu_r = stoich_r(g) % reactant stoichiometric coefficient mstix
nu_p = stoich_p(g) % product stoichiometric coefficient mstix
nu_net = stoich_net(g) % net (product - reactant) stoichiometric
% coefficient mstix
% For any of these, the (k,i) matrix element is the stoichiometric
% coefficient of species k in reaction i. Since these coefficient
% matrices are very sparse, they are implemented as MATLAB sparse
% matrices.
% 2) Reaction rates of progress
% Methods rop_f, rop_r, and rop_net return column vectors containing
% the forward, reverse, and net (forward - reverse) rates of
% progress, respectively, for all reactions.
qf = rop_f(g);
qr = rop_r(g);
qn = rop_net(g);
rop = [qf, qr, qn]
% This plots the rates of progress
figure(1);
bar(rop);
legend('forward','reverse','net');
% 3) Species production rates
% Methods creationRates, destructionRates, and netProdRates return
% column vectors containing the creation, destruction, and net
% production (creation - destruction) rates, respectively, for all species.
cdot = creationRates(g);
ddot = destructionRates(g);
wdot = netProdRates(g);
rates = [cdot, ddot, wdot]
% This plots the production rates
figure(2);
bar(rates);
legend('creation','destruction','net');
% For comparison, the production rates may also be computed
% directly from the rates of progress and stoichiometric
% coefficients.
cdot2 = nu_p*qf + nu_r*qr;
creation = [cdot, cdot2, cdot - cdot2]
ddot2 = nu_r*qf + nu_p*qr;
destruction = [ddot, ddot2, ddot - ddot2]
wdot2 = nu_net * qn;
net = [wdot, wdot2, wdot - wdot2]
% 4) Reaction equations
e8 = reactionEqn(g,8) % equation for reaction 8
e1_10 = reactionEqn(g,1:10) % equation for rxns 1 - 10
eqs = reactionEqn(g) % all equations
% 5) Equilibrium constants
% The equilibrium constants are computed in concentration units,
% with concentrations in kmol/m^3.
kc = equil_Kc(g);
for i = 1:nReactions(g)
fprintf('%50s %13.5g', eqs{i}, kc(i))
end
% 6) Multipliers
% For each reaction, a multiplier may be specified that is applied
% to the forward rate coefficient. By default, the multiplier is
% 1.0 for all reactions.
for i = 1:nReactions(g)
setMultiplier(g, i, 2*i);
m = multiplier(g, i);
end
clear all
cleanup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

73
samples/matlab/tut6.m
View File

@ -1,73 +0,0 @@
% Tutorial 6: Transport properties
%
% Topics:
% - mixture-averaged and multicomponent models
% - viscosity
% - thermal conductivity
% - binary diffusion coefficients
% - mixture-averaged diffusion coefficients
% - multicomponent diffusion coefficients
% - thermal diffusion coefficients
%
% Keywords: tutorial, transport
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Methods are provided to compute transport properties. By
% default, calculation of transport properties is not enabled. If
% transport properties are required, the transport model must be
% specified when the gas mixture object is constructed.
% Currently, two models are implemented. Both are based on kinetic
% theory expressions, and follow the approach described in Dixon-Lewis
% (1968) and Kee, Coltrin, and Glarborg (2003). The first is a full
% multicomponent formulation, and the second is a simplification that
% uses expressions derived for mixtures with a small number of species
% (1 to 3), using approximate mixture rules to average over
% composition.
% To use the multicomponent model with GRI-Mech 3.0, call function
% GRI30 as follows:
g1 = GRI30('multicomponent')
% To use the mixture-averaged model:
g2 = GRI30('mixture-averaged')
% Both models use a mixture-averaged formulation for the viscosity.
visc = [viscosity(g1), viscosity(g2)]
% The thermal conductivity differs, however.
lambda = [thermalConductivity(g1), thermalConductivity(g2)]
% Binary diffusion coefficients
bdiff1 = binDiffCoeffs(g1)
bdiff2 = binDiffCoeffs(g2)
% Mixture-averaged diffusion coefficients. For convenience, the
% multicomponent model implements mixture-averaged diffusion
% coefficients too.
dmix2 = mixDiffCoeffs(g1)
dmix1 = mixDiffCoeffs(g2)
% Multicomponent diffusion coefficients. These are only implemented
% if the multicomponent model is used.
dmulti = multiDiffCoeffs(g1)
% Thermal diffusion coefficients. These are only implemented with the
% multicomponent model. These will be very close to zero, since
% the composition is pure H2.
dt = thermalDiffCoeffs(g1)
% Now change the composition and re-evaluate
set(g1,'X',ones(nSpecies(g1),1));
dt = thermalDiffCoeffs(g1)
% Note that there are no singularities for pure gases. This is
% because a very small positive value is added to all mole
% fractions for the purpose of computing transport properties.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
cleanup

35
samples/matlab/tut7.m
View File

@ -1,35 +0,0 @@
% Tutorial 7: Thermodynamic Properties
%
% Keywords: tutorial, thermodynamics
help tut7
% A variety of thermodynamic property methods are provided.
gas = Air
set(gas,'T',800,'P',oneatm)
% temperature, pressure, density
T = temperature(gas)
P = pressure(gas)
rho = density(gas)
n = molarDensity(gas)
% species non-dimensional properties
hrt = enthalpies_RT(gas) % vector of h_k/RT
% mixture properties per mole
hmole = enthalpy_mole(gas)
umole = intEnergy_mole(gas)
smole = entropy_mole(gas)
gmole = gibbs_mole(gas)
% mixture properties per unit mass
hmass = enthalpy_mass(gas)
umass = intEnergy_mass(gas)
smass = entropy_mass(gas)
gmass = gibbs_mass(gas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
cleanup

4
samples/matlab_experimental/plotSolution.m
View File

@ -4,7 +4,7 @@ function plotSolution(s, domain, component)
% >> plotSolution(s, domain, component)
%
% :s:
% Instance of class :mat:class:`Sim1D`.
% Instance of class Matlab class `Sim1D`.
% :domain:
% Name of domain from which the component should be retrieved.
% :component:
@ -17,4 +17,4 @@ function plotSolution(s, domain, component)
plot(z, x);
xlabel('z (m)');
ylabel(component);
end
end