mirror of
https://github.com/Cantera/cantera.git
synced 2025-02-25 18:55:29 -06:00
Reformatted all matlab interface and sample files to Cantera contribution guide.
Removed most nested IF statements.
This commit is contained in:
ssun30
Ray Speth
parent
0351f3b0f8
commit
9e64da1bcd
@ -10,9 +10,11 @@ function m = AxisymmetricFlow(gas, id)
|
||||
% Domain1D instance representing an axisymmetric flow.
|
||||
%
|
||||
m = Domain1D('StagnationFlow', gas);
|
||||
|
||||
if nargin == 1
|
||||
m.setID('flow');
|
||||
else
|
||||
m.setID(id);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -29,21 +29,27 @@ function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidiz
|
||||
if nargin ~= 6
|
||||
error('CounterFlowDiffusionFlame expects six input arguments.');
|
||||
end
|
||||
|
||||
if ~tp_f.isIdealGas
|
||||
error('Fuel gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
|
||||
if ~tp_o.isIdealGas
|
||||
error('Oxidizer gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
|
||||
if ~left.isInlet
|
||||
error('Left inlet object of wrong type.');
|
||||
end
|
||||
|
||||
if ~flow.isFlow
|
||||
error('Flow object of wrong type.');
|
||||
end
|
||||
|
||||
if ~right.isInlet
|
||||
error('Right inlet object of wrong type.');
|
||||
end
|
||||
|
||||
if ~ischar(oxidizer)
|
||||
error('Oxidizer name must be of format character.');
|
||||
end
|
||||
@ -74,7 +80,7 @@ function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidiz
|
||||
% conditions. The stoichiometric mixture fraction, Zst, is then
|
||||
% calculated.
|
||||
|
||||
sFuel = elMoles(tp_f, 'O')- 2 * elMoles(tp_f, 'C') - 0.5 * elMoles(tp_f, 'H');
|
||||
sFuel = elMoles(tp_f, 'O') - 2 * elMoles(tp_f, 'C') - 0.5 * elMoles(tp_f, 'H');
|
||||
sOx = elMoles(tp_o, 'O') - 2 * elMoles(tp_o, 'C') - 0.5 * elMoles(tp_o, 'H');
|
||||
phi = sFuel / sOx;
|
||||
zst = 1.0 / (1.0 - phi);
|
||||
@ -101,6 +107,7 @@ function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidiz
|
||||
|
||||
% Initialize stoichiometric mass fraction cell with first SP:Y value.
|
||||
ystoich = [y_stoich{1}];
|
||||
|
||||
for i = 2:nsp
|
||||
% Update cell to have format similar to N2:Yst,O2:Yst,...
|
||||
ystoich = [ystoich ',', y_stoich{i}];
|
||||
@ -129,7 +136,7 @@ function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidiz
|
||||
f = sqrt(a / (2.0 * diff(ioxidizer)));
|
||||
x0num = sqrt(uleft * mdotl) * dz;
|
||||
x0den = sqrt(uleft * mdotr) + sqrt(uright * mdotr);
|
||||
x0 = x0num/x0den;
|
||||
x0 = x0num / x0den;
|
||||
|
||||
%% Calculate initial values of temperature and mass fractions.
|
||||
% These values to be used for energy equation solution. Method is based
|
||||
@ -149,18 +156,25 @@ function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidiz
|
||||
zm(j) = zmix;
|
||||
u(j) = a * (x0 - zz(j)); % Axial velocity.
|
||||
v(j) = a; % Radial velocity.
|
||||
|
||||
if zmix > zst
|
||||
|
||||
for n = 1:nsp
|
||||
y(j,n) = yeq(n) + (zmix - zst) * (yf(n) - yeq(n)) / (1.0 - zst);
|
||||
y(j, n) = yeq(n) + (zmix - zst) * (yf(n) - yeq(n)) / (1.0 - zst);
|
||||
end
|
||||
|
||||
t(j) = teq + (tf - teq) * (zmix - zst) / (1.0 - zst);
|
||||
else
|
||||
|
||||
for n = 1:nsp
|
||||
y(j,n) = yox(n) + zmix * (yeq(n) - yox(n)) / zst;
|
||||
y(j, n) = yox(n) + zmix * (yeq(n) - yox(n)) / zst;
|
||||
end
|
||||
|
||||
t(j) = tox + zmix * (teq - tox) / zst;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
zrel = zz / dz;
|
||||
|
||||
%% Create the flame stack.
|
||||
@ -168,11 +182,13 @@ function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidiz
|
||||
% radial velocities, temperature, and mass fractions calculated above.
|
||||
flame = Sim1D([left flow right]);
|
||||
flame.setProfile(2, {'velocity', 'spread_rate'}, [zrel; u; v]);
|
||||
flame.setProfile(2, 'T', [zrel; t] );
|
||||
flame.setProfile(2, 'T', [zrel; t]);
|
||||
|
||||
for n = 1:nsp
|
||||
nm = tp_f.speciesName(n);
|
||||
flame.setProfile(2, nm, [zrel; transpose(y(:,n))])
|
||||
flame.setProfile(2, nm, [zrel; transpose(y(:, n))])
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% Define elMoles function
|
||||
@ -184,9 +200,11 @@ function moles = elMoles(tp, element)
|
||||
if nargin ~= 2
|
||||
error('elMoles expects two input arguments.');
|
||||
end
|
||||
|
||||
if ~tp.isIdealGas
|
||||
error('Gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
|
||||
if ~ischar(element)
|
||||
error('Element name must be of format character.');
|
||||
end
|
||||
|
||||
@ -126,6 +126,7 @@ classdef Domain1D < handle
|
||||
checklib;
|
||||
|
||||
if nargin == 1
|
||||
|
||||
if strcmp(a, 'Inlet1D')
|
||||
d.domainID = callct('inlet_new');
|
||||
elseif strcmp(a, 'Surf1D')
|
||||
@ -139,23 +140,29 @@ classdef Domain1D < handle
|
||||
else
|
||||
error('Not enough arguments for that job number');
|
||||
end
|
||||
|
||||
elseif nargin == 2
|
||||
% a stagnation flow
|
||||
if strcmp(a, 'StagnationFlow')
|
||||
|
||||
if isa(b, 'Solution')
|
||||
d.domainID = callct('stflow_new', ...
|
||||
b.tpID, b.kinID, b.trID, 1);
|
||||
b.tpID, b.kinID, b.trID, 1);
|
||||
else
|
||||
error('Wrong argument type. Expecting instance of class Solution.');
|
||||
end
|
||||
|
||||
elseif strcmp(a, 'AxisymmetricFlow')
|
||||
|
||||
if isa(b, 'Solution')
|
||||
d.domainID = callct('stflow_new', ...
|
||||
b.tpID, b.kinID, b.trID, 2);
|
||||
else
|
||||
error('Wrong argument type. Expecting instance of class Solution.');
|
||||
end
|
||||
|
||||
elseif strcmp(a, 'ReactingSurface')
|
||||
|
||||
if isa(b, 'Interface')
|
||||
d.domainID = callct('reactingsurf_new');
|
||||
callct('reactingsurf_setkineticsmgr', ...
|
||||
@ -163,10 +170,13 @@ classdef Domain1D < handle
|
||||
else
|
||||
error('Wrong argument type. Expecting instance of class Interface.');
|
||||
end
|
||||
|
||||
else
|
||||
error('Wrong object type.');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
d.type = a;
|
||||
end
|
||||
|
||||
@ -269,14 +279,18 @@ classdef Domain1D < handle
|
||||
n = name;
|
||||
else
|
||||
n = callct('domain_componentIndex', ...
|
||||
d.domainID, name);
|
||||
d.domainID, name);
|
||||
|
||||
if n >= 0
|
||||
n = n+1;
|
||||
n = n + 1;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if n <= 0
|
||||
error('Component not found');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function s = componentName(d, index)
|
||||
@ -295,21 +309,26 @@ classdef Domain1D < handle
|
||||
|
||||
n = length(index);
|
||||
s = cell(1, n);
|
||||
|
||||
for i = 1:n
|
||||
id = index(i)-1;
|
||||
id = index(i) - 1;
|
||||
output = callct2('domain_componentName', d.domainID, id);
|
||||
s{i} = output;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function i = get.domainIndex(d)
|
||||
i = callct('domain_index', d.domainID);
|
||||
|
||||
if i >= 0
|
||||
i = i + 1;
|
||||
end
|
||||
|
||||
if i <= 0
|
||||
error('Domain not found');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function i = get.domainType(d)
|
||||
@ -332,16 +351,21 @@ classdef Domain1D < handle
|
||||
if nargin == 1
|
||||
np = d.nPoints;
|
||||
zz = zeros(1, np);
|
||||
|
||||
for i = 1:np
|
||||
zz(i) = callct('domain_grid', d.domainID, i-1);
|
||||
zz(i) = callct('domain_grid', d.domainID, i - 1);
|
||||
end
|
||||
|
||||
else
|
||||
m = length(n);
|
||||
zz = zeros(1, m);
|
||||
|
||||
for i = 1:m
|
||||
zz(i) = callct('domain_grid', d.domainID, n(i)-1);
|
||||
zz(i) = callct('domain_grid', d.domainID, n(i) - 1);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function a = get.isFlow(d)
|
||||
@ -351,6 +375,7 @@ classdef Domain1D < handle
|
||||
a = 1;
|
||||
else a = 0;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function a = get.isInlet(d)
|
||||
@ -360,6 +385,7 @@ classdef Domain1D < handle
|
||||
a = 1;
|
||||
else a = 0;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function a = get.isSurface(d)
|
||||
@ -369,6 +395,7 @@ classdef Domain1D < handle
|
||||
a = 1;
|
||||
else a = 0;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function mdot = get.massFlux(d)
|
||||
@ -397,9 +424,10 @@ classdef Domain1D < handle
|
||||
end
|
||||
|
||||
if d.isInlet
|
||||
y = callct('bdry_massFraction', d.domainID, k-1);
|
||||
y = callct('bdry_massFraction', d.domainID, k - 1);
|
||||
else error('Input domain must be an inlet');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function n = get.nComponents(d)
|
||||
@ -455,7 +483,7 @@ classdef Domain1D < handle
|
||||
%
|
||||
|
||||
n = d.componentIndex(component);
|
||||
callct('domain_setBounds', d.domainID, n-1, lower, upper);
|
||||
callct('domain_setBounds', d.domainID, n - 1, lower, upper);
|
||||
end
|
||||
|
||||
function setCoverageEqs(d, onoff)
|
||||
@ -477,7 +505,9 @@ classdef Domain1D < handle
|
||||
end
|
||||
|
||||
ion = -1;
|
||||
if isa(onoff,'char')
|
||||
|
||||
if isa(onoff, 'char')
|
||||
|
||||
if strcmp(onoff, 'on') || strcmp(onoff, 'yes')
|
||||
ion = 1;
|
||||
elseif strcmp(onoff, 'off') || strcmp(onoff, 'no')
|
||||
@ -485,9 +515,11 @@ classdef Domain1D < handle
|
||||
else
|
||||
error(strcat('unknown option: ', onoff))
|
||||
end
|
||||
|
||||
elseif isa(onoff, 'numeric')
|
||||
ion = onoff;
|
||||
end
|
||||
|
||||
callct('reactingsurf_enableCoverageEqs', d.domainID, ion);
|
||||
end
|
||||
|
||||
@ -508,6 +540,7 @@ classdef Domain1D < handle
|
||||
% is specified.
|
||||
%
|
||||
sz = size(profile);
|
||||
|
||||
if sz(1) == 2
|
||||
l = length(profile(1, :));
|
||||
callct('stflow_setFixedTempProfile', d.domainID, ...
|
||||
@ -518,6 +551,7 @@ classdef Domain1D < handle
|
||||
l, profile(:, 1), l, profile(:, 2));
|
||||
else error('Wrong temperature profile array shape.');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function setID(d, id)
|
||||
@ -578,22 +612,27 @@ classdef Domain1D < handle
|
||||
%
|
||||
if strcmp(component, 'default')
|
||||
nc = d.nComponents;
|
||||
|
||||
for ii = 1:nc
|
||||
callct('domain_setSteadyTolerances', ...
|
||||
d.domainID, ii-1, rtol, atol);
|
||||
d.domainID, ii - 1, rtol, atol);
|
||||
end
|
||||
|
||||
elseif iscell(component)
|
||||
nc = length(component);
|
||||
|
||||
for ii = 1:nc
|
||||
n = d.componentIndex(component{ii});
|
||||
callct('domain_setSteadyTolerances', ...
|
||||
d.domainID, n, rtol, atol);
|
||||
end
|
||||
|
||||
else
|
||||
n = d.componentIndex(component);
|
||||
callct('domain_setSteadyTolerances', ...
|
||||
d.domainID, n, rtol, atol);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function setTransientTolerances(d, component, rtol, atol)
|
||||
@ -614,22 +653,27 @@ classdef Domain1D < handle
|
||||
%
|
||||
if strcmp(component, 'default')
|
||||
nc = d.nComponents;
|
||||
|
||||
for ii = 1:nc
|
||||
callct('domain_setTransientTolerances', ...
|
||||
d.domainID, ii-1, rtol, atol);
|
||||
d.domainID, ii - 1, rtol, atol);
|
||||
end
|
||||
|
||||
elseif iscell(component)
|
||||
nc = length(component);
|
||||
|
||||
for ii = 1:nc
|
||||
n = d.componentIndex(component{ii});
|
||||
callct('domain_setTransientTolerances', ...
|
||||
d.domainID, n, rtol, atol);
|
||||
end
|
||||
|
||||
else
|
||||
n = d.componentIndex(component);
|
||||
callct('domain_setTransientTolerances', ...
|
||||
d.domainID, n, rtol, atol);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function setTransport(d, itr)
|
||||
@ -657,18 +701,23 @@ classdef Domain1D < handle
|
||||
end
|
||||
|
||||
function set.T(d, t)
|
||||
|
||||
if t <= 0
|
||||
error('The temperature must be positive');
|
||||
end
|
||||
|
||||
callct('bdry_setTemperature', d.domainID, t);
|
||||
end
|
||||
|
||||
function set.P(d, p)
|
||||
|
||||
if p <= 0
|
||||
error('The pressure must be positive');
|
||||
end
|
||||
|
||||
callct('stflow_setPressure', d.domainID, p);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -11,9 +11,11 @@ function m = FreeFlame(gas, id)
|
||||
% adiabatic flame
|
||||
%
|
||||
m = Domain1D('StagnationFlow', gas, 2);
|
||||
|
||||
if nargin == 1
|
||||
m.setID('flame');
|
||||
else
|
||||
m.setID(id);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -10,9 +10,11 @@ function m = Inlet(id)
|
||||
% Instance of class :mat:func:`Domain1D` representing an inlet.
|
||||
%
|
||||
m = Domain1D('Inlet1D');
|
||||
|
||||
if nargin == 0
|
||||
m.setID('inlet');
|
||||
else
|
||||
m.setID(id);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -8,9 +8,11 @@ function m = Outlet(id)
|
||||
% Instance of :mat:func:`Domain1D` representing an outlet.
|
||||
%
|
||||
m = Domain1D('Outlet1D');
|
||||
|
||||
if nargin == 0
|
||||
m.setID('outlet');
|
||||
else
|
||||
m.setID(id);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -7,9 +7,11 @@ function m = OutletRes(id)
|
||||
% reservoir.
|
||||
%
|
||||
m = Domain1D('OutletRes');
|
||||
|
||||
if nargin == 0
|
||||
m.setID('outletres');
|
||||
else
|
||||
m.setID(id);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -30,17 +30,21 @@ classdef Sim1D < handle
|
||||
|
||||
s.stID = -1;
|
||||
s.domains = domains;
|
||||
if nargin == 1
|
||||
nd = length(domains);
|
||||
ids = zeros(1, nd);
|
||||
for n=1:nd
|
||||
ids(n) = domains(n).domainID;
|
||||
end
|
||||
s.stID = callct('sim1D_new', nd, ids);
|
||||
else
|
||||
|
||||
if nargin ~= 1
|
||||
help(Sim1D);
|
||||
error('Wrong number of parameters.');
|
||||
end
|
||||
|
||||
nd = length(domains);
|
||||
ids = zeros(1, nd);
|
||||
|
||||
for n = 1:nd
|
||||
ids(n) = domains(n).domainID;
|
||||
end
|
||||
|
||||
s.stID = callct('sim1D_new', nd, ids);
|
||||
|
||||
end
|
||||
|
||||
%% Sim1D Class Destructor
|
||||
@ -59,6 +63,7 @@ classdef Sim1D < handle
|
||||
if nargin == 1
|
||||
fname = '-';
|
||||
end
|
||||
|
||||
callct('sim1D_showSolution', s.stID, fname);
|
||||
end
|
||||
|
||||
@ -132,6 +137,7 @@ classdef Sim1D < handle
|
||||
elseif nargin == 3
|
||||
desc = '--';
|
||||
end
|
||||
|
||||
callct('sim1D_save', s.stID, fname, id, desc);
|
||||
end
|
||||
|
||||
@ -157,23 +163,31 @@ classdef Sim1D < handle
|
||||
idom = s.stackIndex(domain);
|
||||
d = s.domains(idom);
|
||||
np = d.nPoints;
|
||||
|
||||
if nargin == 3
|
||||
icomp = d.componentIndex(component);
|
||||
x = zeros(1, np);
|
||||
|
||||
for n = 1:np
|
||||
x(n) = callct('sim1D_value', s.stID, ...
|
||||
idom - 1, icomp - 1, n - 1);
|
||||
end
|
||||
|
||||
else
|
||||
nc = d.nComponents;
|
||||
x = zeros(nc, np);
|
||||
|
||||
for m = 1:nc
|
||||
|
||||
for n = 1:np
|
||||
x(m, n) = callct('sim1D_value', s.stID, ...
|
||||
idom - 1, m - 1, n - 1);
|
||||
idom - 1, m - 1, n - 1);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function solve(s, loglevel, refineGrid)
|
||||
@ -235,12 +249,15 @@ classdef Sim1D < handle
|
||||
n = name;
|
||||
else
|
||||
n = callct('sim1D_domainIndex', s.stID, name);
|
||||
|
||||
if n >= 0
|
||||
n = n+1;
|
||||
n = n + 1;
|
||||
else
|
||||
error('Domain not found');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function z = grid(s, name)
|
||||
@ -287,12 +304,16 @@ classdef Sim1D < handle
|
||||
|
||||
r = zeros(nc, np);
|
||||
callct('sim1D_eval', s.stID, rdt, count);
|
||||
|
||||
for m = 1:nc
|
||||
|
||||
for n = 1:np
|
||||
r(m, n) = callct('sim1D_workValue', ...
|
||||
s.stID, idom - 1, m - 1, n - 1);
|
||||
s.stID, idom - 1, m - 1, n - 1);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% Sim1D Set Methods
|
||||
@ -310,6 +331,7 @@ classdef Sim1D < handle
|
||||
if T <= 0
|
||||
error('temperature must be positive');
|
||||
end
|
||||
|
||||
callct('sim1D_setFixedTemperature', s.stID, T);
|
||||
end
|
||||
|
||||
@ -343,7 +365,7 @@ classdef Sim1D < handle
|
||||
% Double minimum grid spacing.
|
||||
%
|
||||
|
||||
callct('sim1D_setGridMin', s.stID, domain-1, gridmin);
|
||||
callct('sim1D_setGridMin', s.stID, domain - 1, gridmin);
|
||||
end
|
||||
|
||||
function setMaxJacAge(s, ss_age, ts_age)
|
||||
@ -364,6 +386,7 @@ classdef Sim1D < handle
|
||||
if nargin == 2
|
||||
ts_age = ss_age;
|
||||
end
|
||||
|
||||
callct('sim1D_setMaxJacAge', s.stID, ss_age, ts_age);
|
||||
end
|
||||
|
||||
@ -416,19 +439,23 @@ classdef Sim1D < handle
|
||||
|
||||
np = length(c);
|
||||
sz = size(p);
|
||||
|
||||
if sz(1) == np + 1;
|
||||
|
||||
for j = 1:np
|
||||
ic = d.componentIndex(c{j});
|
||||
callct('sim1D_setProfile', s.stID, ...
|
||||
n - 1, ic - 1, sz(2), p(1, :), sz(2), p(j+1, :));
|
||||
n - 1, ic - 1, sz(2), p(1, :), sz(2), p(j + 1, :));
|
||||
end
|
||||
|
||||
elseif sz(2) == np + 1;
|
||||
ic = d.componentIndex(c{j});
|
||||
callct('sim1D_setProfile', s.stID, ...
|
||||
n - 1, ic - 1, sz(1), p(:, 1), sz(1), p(:, j+1));
|
||||
n - 1, ic - 1, sz(1), p(:, 1), sz(1), p(:, j + 1));
|
||||
else
|
||||
error('Wrong profile shape.');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function setRefineCriteria(s, n, ratio, slope, curve, prune)
|
||||
@ -456,15 +483,19 @@ classdef Sim1D < handle
|
||||
if nargin < 3
|
||||
ratio = 10.0;
|
||||
end
|
||||
|
||||
if nargin < 4
|
||||
slope = 0.8;
|
||||
end
|
||||
|
||||
if nargin < 5
|
||||
curve = 0.8;
|
||||
end
|
||||
|
||||
if nargin < 6
|
||||
prune = -0.1;
|
||||
end
|
||||
|
||||
callct('sim1D_setRefineCriteria', s.stID, ...
|
||||
n - 1, ratio, slope, curve, prune);
|
||||
end
|
||||
@ -516,9 +547,9 @@ classdef Sim1D < handle
|
||||
%
|
||||
|
||||
callct('sim1D_setValue', s.stID, ...
|
||||
n - 1, comp - 1, localPoints - 1, v);
|
||||
n - 1, comp - 1, localPoints - 1, v);
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -11,15 +11,19 @@ function m = Surface(id, surface_mech)
|
||||
% Instance of class :mat:func:`Domain1D` representing a
|
||||
% non-reacting or reacting surface.
|
||||
%
|
||||
if nargin < 2
|
||||
m = Domain1D('Surf1D');
|
||||
if nargin == 0
|
||||
m.setID('surface');
|
||||
elseif nargin == 1
|
||||
m.setID(id);
|
||||
end
|
||||
else
|
||||
|
||||
if nargin > 1
|
||||
m = Domain1D('ReactingSurface', surface_mech);
|
||||
m.setID(id);
|
||||
return
|
||||
end
|
||||
|
||||
m = Domain1D('Surf1D');
|
||||
|
||||
if nargin == 0
|
||||
m.setID('surface');
|
||||
elseif nargin == 1
|
||||
m.setID(id);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -9,9 +9,11 @@ function m = SymmPlane(id)
|
||||
% plane.
|
||||
%
|
||||
m = Domain1D(Symm1D);
|
||||
|
||||
if nargin == 0
|
||||
m.setID('symmetry_plane');
|
||||
else
|
||||
m.setID(id);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -43,6 +43,7 @@ classdef Interface < handle & ThermoPhase & Kinetics
|
||||
|
||||
t = ThermoPhase(src, id);
|
||||
s@ThermoPhase(src, id);
|
||||
|
||||
if nargin == 2
|
||||
args = {};
|
||||
elseif nargin == 3
|
||||
@ -54,6 +55,7 @@ classdef Interface < handle & ThermoPhase & Kinetics
|
||||
elseif nargin == 6
|
||||
args = {p1, p2, p3, p4};
|
||||
end
|
||||
|
||||
s@Kinetics(t, src, id, args{:});
|
||||
s.tpID = t.tpID;
|
||||
end
|
||||
@ -116,17 +118,20 @@ classdef Interface < handle & ThermoPhase & Kinetics
|
||||
|
||||
if isa(cov, 'double')
|
||||
sz = length(cov);
|
||||
if sz == nsp
|
||||
if ((m == nsp && n == 1) || (m == 1 & n == nsp))
|
||||
callct('surf_setCoverages', surfID, cov, norm_flag);
|
||||
else error('wrong size for coverage array');
|
||||
end
|
||||
else
|
||||
|
||||
if sz ~= nsp
|
||||
error('wrong size for coverage array');
|
||||
end
|
||||
|
||||
if ((m == nsp && n == 1) || (m == 1 & n == nsp))
|
||||
callct('surf_setCoverages', surfID, cov, norm_flag);
|
||||
else error('wrong size for coverage array');
|
||||
end
|
||||
|
||||
elseif isa(cov, 'char')
|
||||
callct('surf_setCoveragesByName', surfID, cov);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function set.siteDensity(s, d)
|
||||
@ -145,5 +150,5 @@ classdef Interface < handle & ThermoPhase & Kinetics
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -113,6 +113,7 @@ classdef Kinetics < handle
|
||||
inb2 = -1;
|
||||
inb3 = -1;
|
||||
inb4 = -1;
|
||||
|
||||
if nargin == 2
|
||||
id = '-';
|
||||
end
|
||||
@ -120,18 +121,23 @@ classdef Kinetics < handle
|
||||
% get the integer indices used to find the stored objects
|
||||
% representing the phases participating in the mechanism
|
||||
iph = ph.tpID;
|
||||
|
||||
if nargin > 6
|
||||
inb4 = n4.tpID;
|
||||
end
|
||||
|
||||
if nargin > 5
|
||||
inb3 = n3.tpID;
|
||||
end
|
||||
|
||||
if nargin > 4
|
||||
inb2 = n2.tpID;
|
||||
end
|
||||
|
||||
if nargin > 3
|
||||
inb1 = n1.tpID;
|
||||
if nargin > 4
|
||||
inb2 = n2.tpID;
|
||||
if nargin > 5
|
||||
inb3 = n3.tpID;
|
||||
if nargin > 6
|
||||
inb4 = n4.tpID;
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
kin.kinID = callct('kin_newFromFile', src, id, ...
|
||||
iph, inb1, inb2, inb3, inb4);
|
||||
end
|
||||
@ -168,7 +174,7 @@ classdef Kinetics < handle
|
||||
% :return:
|
||||
% Multiplier of the rate of progress of reaction irxn.
|
||||
|
||||
n = callct('kin_multiplier', kin.kinID, irxn-1);
|
||||
n = callct('kin_multiplier', kin.kinID, irxn - 1);
|
||||
end
|
||||
|
||||
function n = get.nPhases(kin)
|
||||
@ -218,6 +224,7 @@ classdef Kinetics < handle
|
||||
nsp = kin.nTotalSpecies;
|
||||
nr = kin.nReactions;
|
||||
temp = sparse(nsp, nr);
|
||||
|
||||
if nargin == 1
|
||||
krange = 1:nsp;
|
||||
irange = 1:nr;
|
||||
@ -228,13 +235,17 @@ classdef Kinetics < handle
|
||||
end
|
||||
|
||||
for k = krange
|
||||
|
||||
for i = irange
|
||||
t = callct('kin_reactantStoichCoeff', ...
|
||||
kin.kinID, k-1, i-1);
|
||||
kin.kinID, k - 1, i - 1);
|
||||
|
||||
if t ~= 0.0
|
||||
temp(k, i) = t;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
n = temp;
|
||||
@ -258,6 +269,7 @@ classdef Kinetics < handle
|
||||
nsp = kin.nTotalSpecies;
|
||||
nr = kin.nReactions;
|
||||
temp = sparse(nsp, nr);
|
||||
|
||||
if nargin == 1
|
||||
krange = 1:nsp;
|
||||
irange = 1:nr;
|
||||
@ -268,13 +280,17 @@ classdef Kinetics < handle
|
||||
end
|
||||
|
||||
for k = krange
|
||||
|
||||
for i = irange
|
||||
t = callct('kin_productStoichCoeff', ...
|
||||
kin.kinID, k-1, i-1);
|
||||
kin.kinID, k - 1, i - 1);
|
||||
|
||||
if t ~= 0.0
|
||||
temp(k, i) = t;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
n = temp;
|
||||
@ -297,9 +313,11 @@ classdef Kinetics < handle
|
||||
if nargin == 1
|
||||
n = kin.stoichProduct - kin.stoichReactant;
|
||||
elseif nargin == 3
|
||||
n = kin.stoichProduct(species, rxns) - kin.stoichReactant(species, rxns);
|
||||
n = kin.stoichProduct(species, rxns) - ...
|
||||
kin.stoichReactant(species, rxns);
|
||||
else error('stoichNet requires 1 or 3 arguments.');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% Get reaction array attributes
|
||||
@ -375,15 +393,17 @@ classdef Kinetics < handle
|
||||
% :return:
|
||||
% String reaction equation.
|
||||
|
||||
rxn = callct2('kin_getReactionString', kin.kinID, irxn-1);
|
||||
rxn = callct2('kin_getReactionString', kin.kinID, irxn - 1);
|
||||
end
|
||||
|
||||
function rxn = get.reactionEqns(kin)
|
||||
m = kin.nReactions;
|
||||
rxns = cell(1, m);
|
||||
|
||||
for i = 1:m
|
||||
rxn{i} = kin.reactionEqn(i);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function enthalpy = get.dH(kin)
|
||||
@ -492,8 +512,9 @@ classdef Kinetics < handle
|
||||
end
|
||||
|
||||
for i = 1:n
|
||||
callct('kin_setMultiplier', kin.kinID, irxn(i)-1, v);
|
||||
callct('kin_setMultiplier', kin.kinID, irxn(i) - 1, v);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function advanceCoverages(kin, dt)
|
||||
@ -509,4 +530,5 @@ classdef Kinetics < handle
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -256,6 +256,7 @@ classdef ThermoPhase < handle
|
||||
if nargin ~= 2
|
||||
error('ThermoPhase expects 2 input arguments.');
|
||||
end
|
||||
|
||||
tp.tpID = callct('thermo_newFromFile', src, id);
|
||||
tp.basis = 'molar';
|
||||
end
|
||||
@ -277,16 +278,18 @@ classdef ThermoPhase < handle
|
||||
aa = char(ones(1, buflen));
|
||||
ptr = libpointer('cstring', aa);
|
||||
[iok, bb] = calllib(ct, 'thermo_report', tp.tpID, buflen, ptr, 1);
|
||||
|
||||
if iok < 0
|
||||
error(geterr);
|
||||
else
|
||||
disp(bb);
|
||||
end
|
||||
|
||||
clear aa bb ptr
|
||||
end
|
||||
|
||||
function tp = equilibrate(tp, xy, solver, rtol, maxsteps, ...
|
||||
maxiter, loglevel)
|
||||
maxiter, loglevel)
|
||||
% Set the phase to a state of chemical equilibrium.
|
||||
%
|
||||
% tp.equilibrate(xy, solver, rtol, maxsteps, maxiter, loglevel)
|
||||
@ -321,18 +324,23 @@ classdef ThermoPhase < handle
|
||||
if nargin < 3
|
||||
solver = -1;
|
||||
end
|
||||
|
||||
if nargin < 4
|
||||
rtol = 1.0e-9;
|
||||
end
|
||||
|
||||
if nargin < 5
|
||||
maxsteps = 1000;
|
||||
end
|
||||
|
||||
if nargin < 6
|
||||
maxiter = 100;
|
||||
end
|
||||
|
||||
if nargin < 7
|
||||
loglevel = 0;
|
||||
end
|
||||
|
||||
callct('thermo_equilibrate', tp.tpID, xy, solver, ...
|
||||
rtol, maxsteps, maxiter, loglevel);
|
||||
end
|
||||
@ -343,8 +351,7 @@ classdef ThermoPhase < handle
|
||||
nel = tp.nElements;
|
||||
aa = zeros(1, nel);
|
||||
pt = libpointer('doublePtr', aa);
|
||||
callct('thermo_getAtomicWeights', ...
|
||||
tp.tpID, nel, pt);
|
||||
callct('thermo_getAtomicWeights', tp.tpID, nel, pt);
|
||||
amu = pt.Value;
|
||||
end
|
||||
|
||||
@ -352,8 +359,7 @@ classdef ThermoPhase < handle
|
||||
nsp = tp.nSpecies;
|
||||
yy = zeros(1, nsp);
|
||||
pt = libpointer('doublePtr', yy);
|
||||
callct('thermo_getCharges', ...
|
||||
tp.tpID, nsp, pt);
|
||||
callct('thermo_getCharges', tp.tpID, nsp, pt);
|
||||
e = pt.Value;
|
||||
end
|
||||
|
||||
@ -386,28 +392,35 @@ classdef ThermoPhase < handle
|
||||
if iscell(name)
|
||||
[m, n] = size(name);
|
||||
k = zeros(m, n);
|
||||
|
||||
for i = 1:m
|
||||
|
||||
for j = 1:n
|
||||
k(i, j) = callct('thermo_elementIndex', ...
|
||||
tp.tpID, name{i, j}) + 1;
|
||||
tp.tpID, name{i, j}) + 1;
|
||||
|
||||
if k(i, j) > 1e3
|
||||
warning(['Element ', name{i, j}, ...
|
||||
' does not exist in the phase']);
|
||||
' does not exist in the phase']);
|
||||
k(i, j) = -1;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
elseif ischar(name)
|
||||
k = callct('thermo_elementIndex', ...
|
||||
tp.tpID, name) + 1;
|
||||
k = callct('thermo_elementIndex', tp.tpID, name) + 1;
|
||||
|
||||
if k > 1e3
|
||||
warning(['Element ', name, ...
|
||||
' does not exist in the phase']);
|
||||
warning(['Element ', name, ' does not exist in the phase']);
|
||||
k = -1;
|
||||
end
|
||||
|
||||
else
|
||||
error('name must be either a string or cell array of strings')
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function elMassFrac = elementalMassFraction(tp, element)
|
||||
@ -428,9 +441,11 @@ classdef ThermoPhase < handle
|
||||
if nargin ~= 2
|
||||
error('elementalMassFraction expects two input arguments.');
|
||||
end
|
||||
|
||||
if ~tp.isIdealGas
|
||||
error('Gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
|
||||
if ~ischar(element)
|
||||
error('Element name must be of format character.');
|
||||
end
|
||||
@ -453,12 +468,13 @@ classdef ThermoPhase < handle
|
||||
MW = tp.molecularWeights;
|
||||
% Initialize the element mass fraction as zero.
|
||||
elMassFrac = 0.0;
|
||||
% Use loop to perform summation of elemental mass fraction over all species.
|
||||
% Perform summation of elemental mass fraction over all species.
|
||||
for i = 1:n
|
||||
natoms(i) = tp.nAtoms(spec{i}, element);
|
||||
yy(i) = tp.massFraction(spec{i});
|
||||
elMassFrac = elMassFrac + (natoms(i)*Mel*yy(i))/MW(i);
|
||||
elMassFrac = elMassFrac + (natoms(i) * Mel * yy(i)) / MW(i);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function mmw = get.meanMolecularWeight(tp)
|
||||
@ -473,8 +489,7 @@ classdef ThermoPhase < handle
|
||||
nsp = tp.nSpecies;
|
||||
yy = zeros(1, nsp);
|
||||
pt = libpointer('doublePtr', yy);
|
||||
callct('thermo_getMolecularWeights', ...
|
||||
tp.tpID, nsp, pt);
|
||||
callct('thermo_getMolecularWeights', tp.tpID, nsp, pt);
|
||||
mw = pt.Value;
|
||||
end
|
||||
|
||||
@ -493,27 +508,27 @@ classdef ThermoPhase < handle
|
||||
% :return:
|
||||
% Number of atoms of element ``m`` in species ``k``.
|
||||
|
||||
if nargin == 3
|
||||
if ischar(species)
|
||||
k = tp.speciesIndex(species);
|
||||
else k = species;
|
||||
end
|
||||
if k < 0
|
||||
n = -1;
|
||||
return
|
||||
end
|
||||
if ischar(element)
|
||||
m = tp.elementIndex(element);
|
||||
else m = element;
|
||||
end
|
||||
if m < 0
|
||||
n = -1;
|
||||
return
|
||||
end
|
||||
n = callct('thermo_nAtoms', tp.tpID, k-1, m-1);
|
||||
else
|
||||
if nargin ~= 3
|
||||
error('Two input arguments required.')
|
||||
end
|
||||
|
||||
if ischar(species)
|
||||
k = tp.speciesIndex(species);
|
||||
else k = species;
|
||||
end
|
||||
|
||||
if ischar(element)
|
||||
m = tp.elementIndex(element);
|
||||
else m = element;
|
||||
end
|
||||
|
||||
if k < 0 | m < 0
|
||||
n = -1;
|
||||
return
|
||||
end
|
||||
|
||||
n = callct('thermo_nAtoms', tp.tpID, k - 1, m - 1)
|
||||
|
||||
end
|
||||
|
||||
function nel = get.nElements(tp)
|
||||
@ -555,40 +570,53 @@ classdef ThermoPhase < handle
|
||||
if iscell(name)
|
||||
[m, n] = size(name);
|
||||
k = zeros(m, n);
|
||||
|
||||
for i = 1:m
|
||||
|
||||
for j = 1:n
|
||||
k(i, j) = callct('thermo_speciesIndex', ...
|
||||
tp.tpID, name{i, j}) + 1;
|
||||
tp.tpID, name{i, j}) + 1;
|
||||
|
||||
if k(i, j) > 1e6
|
||||
warning(['Species ', name{i, j}, ...
|
||||
' does not exist in the phase']);
|
||||
' does not exist in the phase']);
|
||||
k(i, j) = -1;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
elseif ischar(name)
|
||||
k = callct('thermo_speciesIndex', ...
|
||||
tp.tpID, name) + 1;
|
||||
|
||||
if k > 1e6
|
||||
warning(['Species ', name, ...
|
||||
' does not exist in the phase.']);
|
||||
' does not exist in the phase.']);
|
||||
k = -1;
|
||||
end
|
||||
|
||||
else
|
||||
error('Name must be either a string or cell array of strings.')
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function nm = speciesName(tp, k)
|
||||
[m, n] = size(k);
|
||||
nm = cell(m, n);
|
||||
|
||||
for i = 1:m
|
||||
|
||||
for j = 1:n
|
||||
ksp = k(i, j) - 1;
|
||||
output = callct2('thermo_getSpeciesName', tp.tpID, ksp);
|
||||
nm{i, j} = output;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function n = speciesNames(tp)
|
||||
@ -618,15 +646,14 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
function volume = get.V(tp)
|
||||
volume = 1/tp.D;
|
||||
volume = 1 / tp.D;
|
||||
end
|
||||
|
||||
function moleFractions = get.X(tp)
|
||||
nsp = tp.nSpecies;
|
||||
xx = zeros(1, nsp);
|
||||
pt = libpointer('doublePtr', xx);
|
||||
callct('thermo_getMoleFractions', ...
|
||||
tp.tpID, nsp, pt);
|
||||
callct('thermo_getMoleFractions', tp.tpID, nsp, pt);
|
||||
moleFractions = pt.Value;
|
||||
end
|
||||
|
||||
@ -645,31 +672,38 @@ classdef ThermoPhase < handle
|
||||
% Scalar or vector double mole fractions
|
||||
%
|
||||
xarray = tp.X;
|
||||
|
||||
if isa(species, 'char')
|
||||
k = tp.speciesIndex(species);
|
||||
if k > 0
|
||||
|
||||
if k > 0
|
||||
x = xarray(k);
|
||||
else error("species not found.");
|
||||
end
|
||||
|
||||
elseif isa(species, 'cell')
|
||||
n = length(species);
|
||||
x = zeros(1, n);
|
||||
|
||||
for j = 1:n
|
||||
k = tp.speciesIndex(species{j});
|
||||
|
||||
if k > 0
|
||||
x(j) = xarray(k);
|
||||
else error("species not found.");
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function massFractions = get.Y(tp)
|
||||
nsp = tp.nSpecies;
|
||||
yy = zeros(1, nsp);
|
||||
pt = libpointer('doublePtr', yy);
|
||||
callct('thermo_getMassFractions', ...
|
||||
tp.tpID, nsp, pt);
|
||||
callct('thermo_getMassFractions', tp.tpID, nsp, pt);
|
||||
massFractions = pt.Value;
|
||||
end
|
||||
|
||||
@ -688,23 +722,31 @@ classdef ThermoPhase < handle
|
||||
% Scalar or vector double mass fractions
|
||||
%
|
||||
yy = tp.Y;
|
||||
|
||||
if isa(species, 'char')
|
||||
k = tp.speciesIndex(species);
|
||||
if k > 0
|
||||
|
||||
if k > 0
|
||||
y = yy(k);
|
||||
else error("Error: species not found.")
|
||||
end
|
||||
|
||||
elseif isa(species, 'cell')
|
||||
n = length(species);
|
||||
y = zeros(1, n);
|
||||
|
||||
for j = 1:n
|
||||
k = tp.speciesIndex(species{j});
|
||||
|
||||
if k > 0
|
||||
y(j) = yy(k);
|
||||
else error("Error: species not found.")
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% ThermoGet single methods
|
||||
@ -713,27 +755,30 @@ classdef ThermoPhase < handle
|
||||
nsp = tp.nSpecies;
|
||||
xx = zeros(1, nsp);
|
||||
pt = libpointer('doublePtr', xx);
|
||||
callct('thermo_chemPotentials', ...
|
||||
tp.tpID, nsp, pt);
|
||||
callct('thermo_chemPotentials', tp.tpID, nsp, pt);
|
||||
mu = pt.Value;
|
||||
end
|
||||
|
||||
function c = get.cv(tp)
|
||||
|
||||
if strcmp(tp.basis, 'molar')
|
||||
c = callct('thermo_cv_mole', tp.tpID);
|
||||
elseif strcmp(tp.basis, 'mass')
|
||||
c = callct('thermo_cv_mass', tp.tpID);
|
||||
else error("basis not specified");
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function c = get.cp(tp)
|
||||
|
||||
if strcmp(tp.basis, 'molar')
|
||||
c = callct('thermo_cp_mole', tp.tpID);
|
||||
elseif strcmp(tp.basis, 'mass')
|
||||
c = callct('thermo_cp_mass', tp.tpID);
|
||||
else error("basis not specified");
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function d = get.critDensity(tp)
|
||||
@ -757,11 +802,13 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
function v = get.isIdealGas(tp)
|
||||
|
||||
if strcmp(tp.eosType, 'IdealGas')
|
||||
v = 1;
|
||||
else
|
||||
v = 0;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function b = get.isothermalCompressibility(tp)
|
||||
@ -813,12 +860,13 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
function c = get.soundSpeed(tp)
|
||||
|
||||
if tp.isIdealGas
|
||||
tp.basis = 'mass';
|
||||
gamma = tp.cp/tp.cv;
|
||||
gamma = tp.cp / tp.cv;
|
||||
wtm = tp.meanMolecularWeight;
|
||||
r = 8314.4621/wtm;
|
||||
c = sqrt(gamma * r *tp.T);
|
||||
r = 8314.4621 / wtm;
|
||||
c = sqrt(gamma * r * tp.T);
|
||||
else
|
||||
rho0 = tp.D;
|
||||
p0 = tp.P;
|
||||
@ -829,6 +877,7 @@ classdef ThermoPhase < handle
|
||||
dpdrho_s = (p1 - p0) / (rho1 - rho0);
|
||||
c = sqrt(dpdrho_s);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function a = get.thermalExpansionCoeff(tp)
|
||||
@ -840,12 +889,14 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
function enthalpy = get.H(tp)
|
||||
|
||||
if strcmp(tp.basis, 'molar')
|
||||
enthalpy = callct('thermo_enthalpy_mole', tp.tpID);
|
||||
elseif strcmp(tp.basis, 'mass')
|
||||
enthalpy = callct('thermo_enthalpy_mass', tp.tpID);
|
||||
else error("basis not specified");
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function enthalpy = get.enthalpies_RT(tp)
|
||||
@ -857,30 +908,36 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
function entropy = get.S(tp)
|
||||
|
||||
if strcmp(tp.basis, 'molar')
|
||||
entropy = callct('thermo_entropy_mole', tp.tpID);
|
||||
elseif strcmp(tp.basis, 'mass')
|
||||
entropy = callct('thermo_entropy_mass', tp.tpID);
|
||||
else error("basis not specified");
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function intEnergy = get.U(tp)
|
||||
|
||||
if strcmp(tp.basis, 'molar')
|
||||
intEnergy = callct('thermo_intEnergy_mole', tp.tpID);
|
||||
elseif strcmp(tp.basis, 'mass')
|
||||
intEnergy = callct('thermo_intEnergy_mass', tp.tpID);
|
||||
else error("basis not specified");
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function gibbs = get.G(tp)
|
||||
|
||||
if strcmp(tp.basis, 'molar')
|
||||
gibbs = callct('thermo_gibbs_mole', tp.tpID);
|
||||
elseif strcmp(tp.basis, 'mass')
|
||||
gibbs = callct('thermo_gibbs_mass', tp.tpID);
|
||||
else error("basis not specified");
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% ThermoGet multi methods
|
||||
@ -998,7 +1055,7 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
function output = get.TPX(tp)
|
||||
output ={tp.T, tp.P, tp.X};
|
||||
output = {tp.T, tp.P, tp.X};
|
||||
end
|
||||
|
||||
function output = get.TPY(tp)
|
||||
@ -1010,7 +1067,7 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
function output = get.TVX(tp)
|
||||
output ={tp.T, tp.V, tp.X};
|
||||
output = {tp.T, tp.V, tp.X};
|
||||
end
|
||||
|
||||
function output = get.TVY(tp)
|
||||
@ -1110,13 +1167,15 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
function tp = set.basis(tp, b)
|
||||
|
||||
if strcmp(b, 'mole') || strcmp(b, 'molar') ...
|
||||
|| strcmp(b, 'Mole') || strcmp(b, 'Molar')
|
||||
|| strcmp(b, 'Mole') || strcmp(b, 'Molar')
|
||||
tp.basis = 'molar';
|
||||
elseif strcmp(b, 'mass') || strcmp(b, 'Mass')
|
||||
tp.basis = 'mass';
|
||||
else error("Basis must be mass or molar")
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function set.T(tp, temperature)
|
||||
@ -1133,31 +1192,37 @@ classdef ThermoPhase < handle
|
||||
|
||||
function set.X(tp, xx)
|
||||
tol = 1e-9;
|
||||
|
||||
if isa(xx, 'double')
|
||||
nsp = tp.nSpecies;
|
||||
|
||||
if abs(sum(xx) - 1) <= tol
|
||||
norm = 0;
|
||||
else norm = 1;
|
||||
end
|
||||
callct('thermo_setMoleFractions', tp.tpID, ...
|
||||
nsp, xx, norm);
|
||||
|
||||
callct('thermo_setMoleFractions', tp.tpID, nsp, xx, norm);
|
||||
elseif isa(xx, 'char')
|
||||
callct('thermo_setMoleFractionsByName', tp.tpID, xx);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function set.Y(tp, yy)
|
||||
|
||||
if isa(yy, 'double')
|
||||
nsp = tp.nSpecies;
|
||||
|
||||
if abs(sum(yy) -1) <= 1e-9
|
||||
norm = 0;
|
||||
else norm = 1;
|
||||
end
|
||||
callct('thermo_setMassFractions', tp.tpID, ...
|
||||
nsp, yy, norm);
|
||||
|
||||
callct('thermo_setMassFractions', tp.tpID, nsp, yy, norm);
|
||||
elseif isa(yy, 'char')
|
||||
callct('thermo_setMassFractionsByName', tp.tpID, yy);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% PhaseSet multi methods
|
||||
@ -1389,4 +1454,5 @@ classdef ThermoPhase < handle
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -54,6 +54,7 @@ classdef Transport < handle
|
||||
%
|
||||
checklib;
|
||||
tr.trID = 0;
|
||||
|
||||
if nargin == 2
|
||||
model = 'default'
|
||||
end
|
||||
@ -64,17 +65,17 @@ classdef Transport < handle
|
||||
|
||||
if ~isa(tp, 'ThermoPhase')
|
||||
error(['The first argument must be an ', ...
|
||||
'instance of class ThermoPhase']);
|
||||
else
|
||||
tr.th = tp;
|
||||
if strcmp(model, 'default')
|
||||
tr.trID = callct('trans_newDefault', ...
|
||||
tp.tpID, loglevel);
|
||||
else
|
||||
tr.trID = callct('trans_new', model, ...
|
||||
tp.tpID, loglevel);
|
||||
end
|
||||
'instance of class ThermoPhase']);
|
||||
end
|
||||
|
||||
tr.th = tp;
|
||||
|
||||
if strcmp(model, 'default')
|
||||
tr.trID = callct('trans_newDefault', tp.tpID, loglevel);
|
||||
else
|
||||
tr.trID = callct('trans_new', model, tp.tpID, loglevel);
|
||||
end
|
||||
|
||||
tr.tpID = tp.tpID;
|
||||
end
|
||||
|
||||
@ -159,4 +160,5 @@ classdef Transport < handle
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -89,6 +89,7 @@ classdef Func < handle
|
||||
else
|
||||
nn = callct('func_new', itype, n, 0, p);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if strcmp(typ, 'polynomial')
|
||||
@ -110,6 +111,7 @@ classdef Func < handle
|
||||
x.coeffs = p;
|
||||
x.id = newFunc(itype, n.id, p);
|
||||
else
|
||||
|
||||
if strcmp(typ, 'sum')
|
||||
itype = 20;
|
||||
elseif strcmp(typ, 'diff')
|
||||
@ -121,6 +123,7 @@ classdef Func < handle
|
||||
elseif strcmp(typ, 'composite')
|
||||
itype = 60;
|
||||
end
|
||||
|
||||
x.f1 = n;
|
||||
x.f2 = p;
|
||||
x.id = newFunc(itype, n.id, p.id);
|
||||
@ -143,7 +146,7 @@ classdef Func < handle
|
||||
% Display the equation of the input function on the terminal.
|
||||
|
||||
disp(' ');
|
||||
disp([inputname(1),' = '])
|
||||
disp([inputname(1), ' = '])
|
||||
disp(' ');
|
||||
disp([' ' f.char])
|
||||
disp(' ');
|
||||
@ -162,13 +165,16 @@ classdef Func < handle
|
||||
% Returns the value of the function evaluated at ``s``
|
||||
%
|
||||
if strcmp(s.type, '()')
|
||||
ind= s.subs{:};
|
||||
b = zeros(1, length(ind));
|
||||
for k = 1:length(ind)
|
||||
b(k) = callct('func_value', a.id, ind(k));
|
||||
end
|
||||
ind = s.subs{:};
|
||||
b = zeros(1, length(ind));
|
||||
|
||||
for k = 1:length(ind)
|
||||
b(k) = callct('func_value', a.id, ind(k));
|
||||
end
|
||||
|
||||
else error('Specify value for x as p(x)');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function s = get.char(f)
|
||||
@ -181,94 +187,114 @@ classdef Func < handle
|
||||
% :return:
|
||||
% Formatted string displaying the function
|
||||
%
|
||||
if strcmp(f.typ,'sum')
|
||||
s = ['(' char(f.f1) ') + (' char(f.f2) ')'];
|
||||
elseif strcmp(f.typ,'diff')
|
||||
s = ['(' char(f.f1) ') - (' char(f.f2) ')'];
|
||||
elseif strcmp(f.typ,'prod')
|
||||
s = ['(' char(f.f1) ') * (' char(f.f2) ')'];
|
||||
elseif strcmp(f.typ,'ratio')
|
||||
s = ['(' char(f.f1) ') / (' char(f.f2) ')'];
|
||||
elseif all(p.coeffs == 0)
|
||||
s = '0';
|
||||
else
|
||||
if strcmp(f.typ,'polynomial')
|
||||
d = length(p.coeffs) - 1;
|
||||
s = [];
|
||||
nn = 0;
|
||||
for b = f.coeffs;
|
||||
cc(d + 1 - nn) = b;
|
||||
nn = nn + 1;
|
||||
end
|
||||
for a = cc;
|
||||
if a ~= 0;
|
||||
if ~isempty(s)
|
||||
if a > 0
|
||||
s = [s ' + '];
|
||||
else
|
||||
s = [s ' - '];
|
||||
a = -a;
|
||||
end
|
||||
end
|
||||
if a ~= 1 || d == 0
|
||||
s = [s num2str(a)];
|
||||
if d > 0
|
||||
s = [s '*'];
|
||||
end
|
||||
end
|
||||
if d >= 2
|
||||
s = [s 'x^' int2str(d)];
|
||||
elseif d == 1
|
||||
s = [s 'x'];
|
||||
end
|
||||
end
|
||||
d = d - 1;
|
||||
end
|
||||
elseif strcmp(f.typ, 'gaussian')
|
||||
s = ['Gaussian(' num2str(f.coeffs(1)) ',' ...
|
||||
num2str(f.coeffs(2)) ',' ...
|
||||
num2str(f.coeffs(3)) ')'];
|
||||
elseif strcmp(f.typ, 'fourier')
|
||||
c = reshape(f.coeffs, [], 2);
|
||||
Ao = c(1, 1);
|
||||
w = c(1, 2);
|
||||
A = c(2:end, 1);
|
||||
B = c(2:end, 2);
|
||||
N = size(c, 1) - 1;
|
||||
if Ao ~= 0
|
||||
s = num2str(Ao/2);
|
||||
else
|
||||
s = '';
|
||||
end
|
||||
for n=1:N
|
||||
if A(n) ~= 0
|
||||
if A(n) < 0
|
||||
prefix = ' - ';
|
||||
elseif s
|
||||
prefix = ' + ';
|
||||
else
|
||||
prefix = '';
|
||||
end
|
||||
s = [s prefix num2str(abs(A(n))), ...
|
||||
'*cos(' num2str(n*w) '*x)'];
|
||||
end
|
||||
if B(n) ~= 0
|
||||
if B(n) < 0
|
||||
prefix = ' - ';
|
||||
elseif s
|
||||
prefix = ' + ';
|
||||
else
|
||||
prefix = '';
|
||||
end
|
||||
s = [s prefix num2str(abs(B(n))),...
|
||||
'*sin(' num2str(n*w) '*x)'];
|
||||
end
|
||||
end
|
||||
else
|
||||
s = ['*** char not yet implemented for' f.typ ' ***'];
|
||||
end
|
||||
end
|
||||
if strcmp(f.typ, 'sum')
|
||||
s = ['(' char(f.f1) ') + (' char(f.f2) ')'];
|
||||
elseif strcmp(f.typ, 'diff')
|
||||
s = ['(' char(f.f1) ') - (' char(f.f2) ')'];
|
||||
elseif strcmp(f.typ, 'prod')
|
||||
s = ['(' char(f.f1) ') * (' char(f.f2) ')'];
|
||||
elseif strcmp(f.typ, 'ratio')
|
||||
s = ['(' char(f.f1) ') / (' char(f.f2) ')'];
|
||||
elseif all(p.coeffs == 0)
|
||||
s = '0';
|
||||
elseif strcmp(f.typ, 'polynomial')
|
||||
d = length(p.coeffs) - 1;
|
||||
s = [];
|
||||
nn = 0;
|
||||
|
||||
for b = f.coeffs;
|
||||
cc(d + 1 - nn) = b;
|
||||
nn = nn + 1;
|
||||
end
|
||||
|
||||
for a = cc;
|
||||
|
||||
if a ~= 0;
|
||||
|
||||
if ~isempty(s)
|
||||
|
||||
if a > 0
|
||||
s = [s ' + '];
|
||||
else
|
||||
s = [s ' - '];
|
||||
a = -a;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if a ~= 1 || d == 0
|
||||
s = [s num2str(a)];
|
||||
|
||||
if d > 0
|
||||
s = [s '*'];
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if d >= 2
|
||||
s = [s 'x^' int2str(d)];
|
||||
elseif d == 1
|
||||
s = [s 'x'];
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
d = d - 1;
|
||||
end
|
||||
elseif strcmp(f.typ, 'gaussian')
|
||||
s = ['Gaussian(' num2str(f.coeffs(1)) ',' ...
|
||||
num2str(f.coeffs(2)) ',' num2str(f.coeffs(3)) ')'];
|
||||
elseif strcmp(f.typ, 'fourier')
|
||||
c = reshape(f.coeffs, [], 2);
|
||||
Ao = c(1, 1);
|
||||
w = c(1, 2);
|
||||
A = c(2:end, 1);
|
||||
B = c(2:end, 2);
|
||||
N = size(c, 1) - 1;
|
||||
|
||||
if Ao ~= 0
|
||||
s = num2str(Ao / 2);
|
||||
else
|
||||
s = '';
|
||||
end
|
||||
|
||||
for n = 1:N
|
||||
|
||||
if A(n) ~= 0
|
||||
|
||||
if A(n) < 0
|
||||
prefix = ' - ';
|
||||
elseif s
|
||||
prefix = ' + ';
|
||||
else
|
||||
prefix = '';
|
||||
end
|
||||
|
||||
s = [s prefix num2str(abs(A(n))), ...
|
||||
'*cos(' num2str(n * w) '*x)'];
|
||||
end
|
||||
|
||||
if B(n) ~= 0
|
||||
|
||||
if B(n) < 0
|
||||
prefix = ' - ';
|
||||
elseif s
|
||||
prefix = ' + ';
|
||||
else
|
||||
prefix = '';
|
||||
end
|
||||
|
||||
s = [s prefix num2str(abs(B(n))), ...
|
||||
'*sin(' num2str(n * w) '*x)'];
|
||||
end
|
||||
|
||||
end
|
||||
else
|
||||
s = ['*** char not yet implemented for' f.typ ' ***'];
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -17,6 +17,7 @@ function poly = polynom(coeffs)
|
||||
% Instance of class :mat:func:`Func`
|
||||
|
||||
[n m] = size(coeffs);
|
||||
|
||||
if n == 1
|
||||
poly = Func('polynomial', m - 1, coeffs);
|
||||
elseif m == 1
|
||||
@ -24,4 +25,5 @@ function poly = polynom(coeffs)
|
||||
else
|
||||
error('wrong shape for coefficient array');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -161,23 +161,29 @@ classdef Mixture < handle
|
||||
m.phases = phases;
|
||||
|
||||
% If phases are supplied, add them
|
||||
if nargin == 1
|
||||
if ~isa(phases, 'cell')
|
||||
error('Enter phases as a cell array.');
|
||||
end
|
||||
|
||||
% First column contains the phase objects, and the second
|
||||
% column contains the mole numbers of each phase.
|
||||
[np, nc] = size(phases);
|
||||
if nc ~= 2
|
||||
error('Cell array of phases should have each phase on a new row');
|
||||
end
|
||||
for n = 1:np
|
||||
m.addPhase(phases{n, 1}, phases{n, 2});
|
||||
end
|
||||
m.T = phases{n, 1}.T;
|
||||
m.P = phases{n, 1}.P;
|
||||
if nargin == 0
|
||||
return
|
||||
end
|
||||
|
||||
if ~isa(phases, 'cell')
|
||||
error('Enter phases as a cell array.');
|
||||
end
|
||||
|
||||
% First column contains the phase objects, and the second
|
||||
% column contains the mole numbers of each phase.
|
||||
[np, nc] = size(phases);
|
||||
|
||||
if nc ~= 2
|
||||
error('Cell array of phases should have each phase on a new row');
|
||||
end
|
||||
|
||||
for n = 1:np
|
||||
m.addPhase(phases{n, 1}, phases{n, 2});
|
||||
end
|
||||
|
||||
m.T = phases{n, 1}.T;
|
||||
m.P = phases{n, 1}.P;
|
||||
|
||||
end
|
||||
|
||||
%% Mixture Class Destructor
|
||||
@ -195,13 +201,15 @@ classdef Mixture < handle
|
||||
|
||||
calllib(ct, 'mix_updatePhases', m.mixID);
|
||||
[np, nc] = size(m.phases);
|
||||
|
||||
for n = 1:np
|
||||
s = [sprintf('\n******************* Phase %d', n) ...
|
||||
sprintf(' ******************************\n\n Moles: %12.6g', ...
|
||||
phaseMoles(m, n))];
|
||||
phaseMoles(m, n))];
|
||||
disp(s);
|
||||
display(m.phases{n, 1});
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function addPhase(m, phase, moles)
|
||||
@ -221,18 +229,17 @@ classdef Mixture < handle
|
||||
if ~isa(phase, 'ThermoPhase')
|
||||
error('Phase object of wrong type.');
|
||||
end
|
||||
|
||||
if ~isa(moles, 'numeric')
|
||||
error('Number of moles must be numeric');
|
||||
end
|
||||
|
||||
if moles < 0.0
|
||||
error('Negative moles');
|
||||
end
|
||||
|
||||
iok = calllib(ct, 'mix_addPhase', m.mixID, phase.tp_id, ...
|
||||
moles);
|
||||
if iok < 0
|
||||
error('Error adding phase');
|
||||
end
|
||||
calllib(ct, 'mix_addPhase', m.mixID, phase.tp_id, moles);
|
||||
|
||||
end
|
||||
|
||||
%% Mixture Get methods
|
||||
@ -246,7 +253,7 @@ classdef Mixture < handle
|
||||
end
|
||||
|
||||
function n = get.nAtoms(m, e)
|
||||
n = calllib(ct, 'mix_nPhases', m.mixID, k-1, e-1);
|
||||
n = calllib(ct, 'mix_nPhases', m.mixID, k - 1, e - 1);
|
||||
end
|
||||
|
||||
function n = get.nElements(m)
|
||||
@ -266,48 +273,58 @@ classdef Mixture < handle
|
||||
end
|
||||
|
||||
function n = get.speciesIndex(m, k, p)
|
||||
n = calllib(ct, 'mix_speciesIndex', m.mixID, k-1, p-1) + 1;
|
||||
n = calllib(ct, 'mix_speciesIndex', m.mixID, k - 1, p - 1) + 1;
|
||||
end
|
||||
|
||||
function moles = get.elementMoles(m, e)
|
||||
|
||||
if nargin == 2
|
||||
moles = calllib(ct, 'mix_elementMoles', m.mixID, e)
|
||||
elseif nargin == 1
|
||||
nel = m.nElements;
|
||||
moles = zeros(1, nel);
|
||||
|
||||
for i = 1:nel
|
||||
moles(i) = calllib(ct, 'mix_elementMoles', m.mixID, i);
|
||||
end
|
||||
|
||||
else error('wrong number of arguments');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function moles = get.phaseMoles(m, n)
|
||||
|
||||
if nargin == 2
|
||||
moles = calllib(ct, 'mix_phaseMoles', m.mixID, n);
|
||||
elseif nargin == 1
|
||||
np = m.nPhases;
|
||||
moles = zeros(1, np);
|
||||
|
||||
for i = 1:np
|
||||
moles(i) = calllib(ct, 'mix_phaseMoles', ...
|
||||
m.mixID, i);
|
||||
moles(i) = calllib(ct, 'mix_phaseMoles', m.mixID, i);
|
||||
end
|
||||
|
||||
else error('wrong number of arguments');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function moles = speciesMoles(m, k)
|
||||
|
||||
if nargin == 2
|
||||
moles = calllib(ct, 'mix_speciesMoles', m.mixID, k);
|
||||
elseif nargin == 1
|
||||
nsp = m.nSpecies;
|
||||
moles = zeros(1, nsp);
|
||||
|
||||
for i = 1:nsp
|
||||
moles(i) = calllib(ct, 'mix_speciesMoles', ...
|
||||
m.mixID, i);
|
||||
moles(i) = calllib(ct, 'mix_speciesMoles', m.mixID, i);
|
||||
end
|
||||
|
||||
else error('wrong number of arguments');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function mu = get.chemPotentials(m)
|
||||
@ -340,7 +357,7 @@ classdef Mixture < handle
|
||||
% :param moles:
|
||||
% Number of moles to add. Units: kmol
|
||||
%
|
||||
calllib(ct, 'mix_setPhaseMoles', m.mixID, n-1, moles);
|
||||
calllib(ct, 'mix_setPhaseMoles', m.mixID, n - 1, moles);
|
||||
end
|
||||
|
||||
function setSpeciesMoles(m, moles)
|
||||
@ -369,6 +386,7 @@ classdef Mixture < handle
|
||||
else
|
||||
error('The input must be a vector or string!');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function r = equilibrate(m, XY, err, maxsteps, maxiter, loglevel)
|
||||
@ -420,21 +438,27 @@ classdef Mixture < handle
|
||||
if nargin < 6
|
||||
loglevel = 0;
|
||||
end
|
||||
|
||||
if nargin < 5
|
||||
maxiter = 200;
|
||||
end
|
||||
|
||||
if nargin < 4
|
||||
maxsteps = 1000;
|
||||
end
|
||||
|
||||
if nargin < 3
|
||||
err = 1.0e-9;
|
||||
end
|
||||
|
||||
if nargin < 2
|
||||
XY = 'TP'
|
||||
end
|
||||
|
||||
r = calllib(ct, 'mix_equilibrate', m.mixID, XY, err, ...
|
||||
maxsteps, maxiter, loglevel);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -24,5 +24,6 @@ function r = ConstPressureReactor(contents)
|
||||
if nargin == 0
|
||||
contents = 0;
|
||||
end
|
||||
|
||||
r = Reactor(contents, 'ConstPressureReactor');
|
||||
end
|
||||
|
||||
@ -89,15 +89,18 @@ classdef FlowDevice < handle
|
||||
% Instance of class :mat:func:`FlowDevice`
|
||||
%
|
||||
if nargin == 3
|
||||
|
||||
if ~isa(upstream, 'Reactor') || ~isa(downstream, 'Reactor')
|
||||
error(['Flow devices can only be installed between',...
|
||||
error(['Flow devices can only be installed between', ...
|
||||
'reactors or reservoirs']);
|
||||
end
|
||||
|
||||
i = upstream.id;
|
||||
j = downstream.id;
|
||||
callct('flowdev_install', f.id, i, j);
|
||||
else error('install requires 3 arguments');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% Flowdevice Get Methods
|
||||
@ -135,6 +138,7 @@ classdef FlowDevice < handle
|
||||
else
|
||||
error('Time function can only be set for mass flow controllers.');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function setMassFlowRate(f, mdot)
|
||||
@ -154,6 +158,7 @@ classdef FlowDevice < handle
|
||||
else
|
||||
error('Mass flow rate can only be set for mass flow controllers.');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function setMaster(f, d)
|
||||
@ -172,6 +177,7 @@ classdef FlowDevice < handle
|
||||
else
|
||||
error('Master flow device can only be set for pressure controllers.');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function setValveCoeff(f, k)
|
||||
@ -196,8 +202,10 @@ classdef FlowDevice < handle
|
||||
if ~strcmp(f.type, 'Valve')
|
||||
error('Valve coefficient can only be set for valves.');
|
||||
end
|
||||
|
||||
ok = callct('flowdev_setValveCoeff', f.id, k);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -22,5 +22,6 @@ function r = FlowReactor(contents)
|
||||
if nargin == 0
|
||||
contents = 0;
|
||||
end
|
||||
|
||||
r = Reactor(contents, 'FlowReactor');
|
||||
end
|
||||
|
||||
@ -26,5 +26,6 @@ function r = IdealGasConstPressureReactor(contents)
|
||||
if nargin == 0
|
||||
contents = 0;
|
||||
end
|
||||
|
||||
r = Reactor(contents, 'IdealGasConstPressureReactor');
|
||||
end
|
||||
|
||||
@ -23,5 +23,6 @@ function r = IdealGasReactor(contents)
|
||||
if nargin == 0
|
||||
contents = 0;
|
||||
end
|
||||
|
||||
r = Reactor(contents, 'IdealGasReactor');
|
||||
end
|
||||
|
||||
@ -18,7 +18,9 @@ function m = MassFlowController(upstream, downstream)
|
||||
% Instance of class 'FlowDevice'.
|
||||
|
||||
m = FlowDevice('MassFlowController');
|
||||
|
||||
if nargin == 2
|
||||
m.install(upstream, downstream)
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -130,9 +130,11 @@ classdef Reactor < handle
|
||||
|
||||
r.contents = gas;
|
||||
setThermoMgr(r, gas);
|
||||
|
||||
if ~strcmp(r.type, 'Reservoir')
|
||||
setKineticsMgr(r, gas);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% Reactor Get Methods
|
||||
@ -182,9 +184,11 @@ classdef Reactor < handle
|
||||
function massFractions = get.Y(r)
|
||||
nsp = r.contents.nSpecies;
|
||||
massFractions = zeros(1, nsp);
|
||||
|
||||
for i = 1:nsp
|
||||
massFractions(i) = r.massFraction(i);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% Reactor set methods
|
||||
@ -263,6 +267,7 @@ classdef Reactor < handle
|
||||
%
|
||||
|
||||
iflag = -1;
|
||||
|
||||
if strcmp(flag, 'on')
|
||||
iflag = 1;
|
||||
elseif strcmp(flag, 'off')
|
||||
@ -273,6 +278,7 @@ classdef Reactor < handle
|
||||
callct('reactor_setEnergy', r.id, iflag);
|
||||
else error('Input must be "on" or "off".');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function setThermoMgr(r, t)
|
||||
@ -318,4 +324,5 @@ classdef Reactor < handle
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -60,9 +60,11 @@ classdef ReactorNet < handle
|
||||
|
||||
% add reactors
|
||||
nr = length(reactors);
|
||||
|
||||
for i = 1:nr
|
||||
r.addReactor(reactors{i});
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
%% ReactorNet Class Destructor
|
||||
@ -202,11 +204,12 @@ classdef ReactorNet < handle
|
||||
% Instance of class 'reactor'.
|
||||
|
||||
if isa(component, 'string')
|
||||
callct('reactornet_sensitivity', r.id, component, ...
|
||||
p, rxtr.id);
|
||||
callct('reactornet_sensitivity', r.id, component, p, rxtr.id);
|
||||
end
|
||||
|
||||
% Check back on this one to add cases for component type integer.
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -54,19 +54,23 @@ classdef ReactorSurface < handle
|
||||
end
|
||||
|
||||
if nargin >= 2
|
||||
|
||||
if isa(reactor, 'Reactor')
|
||||
s.install(reactor);
|
||||
else
|
||||
warning('Reactor was not installed due to incorrect type');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if nargin >= 3
|
||||
|
||||
if isnumeric(area)
|
||||
s.area = area;
|
||||
else
|
||||
warning('Area was not a number and was not set');
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
@ -131,12 +135,14 @@ classdef ReactorSurface < handle
|
||||
%
|
||||
|
||||
ikin = 0;
|
||||
|
||||
if isa(kin, 'Kinetics')
|
||||
ikin = kin.kinID;
|
||||
end
|
||||
|
||||
callct('reactorsurface_setkinetics', s.surfID, ikin);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -27,5 +27,6 @@ function r = Reservoir(contents)
|
||||
if nargin == 0
|
||||
contents = 0;
|
||||
end
|
||||
|
||||
r = Reactor(contents, 'Reservoir');
|
||||
end
|
||||
|
||||
@ -26,7 +26,9 @@ function v = Valve(upstream, downstream)
|
||||
% Instance of class 'FlowDevice'.
|
||||
|
||||
v = FlowDevice('Valve');
|
||||
|
||||
if nargin == 2
|
||||
v.install(upstream, downstream)
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -212,5 +212,5 @@ classdef Wall < handle
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -10,15 +10,17 @@ else
|
||||
error('Operating System Not Supported!');
|
||||
return;
|
||||
end
|
||||
|
||||
if ~libisloaded(ct)
|
||||
[~,warnings] = loadlibrary([cantera_root, ctname], ...
|
||||
[~, warnings] = loadlibrary([cantera_root, ctname], ...
|
||||
[cantera_root, '/include/cantera/clib/ctmatlab.h'], ...
|
||||
'includepath', [cantera_root '/include'], ...
|
||||
'addheader','ct','addheader','ctfunc', ...
|
||||
'addheader','ctmultiphase','addheader', ...
|
||||
'ctonedim','addheader','ctreactor', ...
|
||||
'addheader','ctrpath','addheader','ctsurf');
|
||||
'addheader', 'ct', 'addheader', 'ctfunc', ...
|
||||
'addheader', 'ctmultiphase', 'addheader', ...
|
||||
'ctonedim', 'addheader', 'ctreactor', ...
|
||||
'addheader', 'ctrpath', 'addheader', 'ctsurf');
|
||||
end
|
||||
|
||||
ct_ver = canteraVersion;
|
||||
sprintf('%s is ready for use.', ct_ver);
|
||||
clear all
|
||||
|
||||
@ -3,4 +3,5 @@
|
||||
if libisloaded(ct)
|
||||
unloadlibrary(ct);
|
||||
end
|
||||
|
||||
disp('Cantera has been unloaded');
|
||||
|
||||
@ -2,16 +2,13 @@ function output = callct(varargin)
|
||||
% Calls Cantera library functions with single outputs and returns
|
||||
% errors if necessary.
|
||||
|
||||
err1 = -1;
|
||||
err2 = -999.999;
|
||||
err3 = double(intmax('uint64'));
|
||||
errorcode = [-1, -999.999, double(intmax('uint64'))];
|
||||
|
||||
funcName = varargin{1};
|
||||
output = calllib(ct, funcName, varargin{2:end});
|
||||
if output == err2 || output == err3
|
||||
output = err1;
|
||||
end
|
||||
if output == err1
|
||||
|
||||
if ismember(output, errorcode)
|
||||
error(geterr);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -6,6 +6,7 @@ function output = callct2(varargin)
|
||||
|
||||
funcName = varargin{1};
|
||||
buflen = calllib(ct, funcName, varargin{2:end}, 0, '');
|
||||
|
||||
if buflen > 0
|
||||
aa = char(ones(1, buflen));
|
||||
ptr = libpointer('cstring', aa);
|
||||
@ -15,7 +16,9 @@ function output = callct2(varargin)
|
||||
else
|
||||
error(geterr);
|
||||
end
|
||||
|
||||
if iok == -err1
|
||||
error(geterr);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -5,4 +5,5 @@ function str = ct()
|
||||
else
|
||||
str = 'libcantera_shared';
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -2,6 +2,7 @@ function e = geterr()
|
||||
% Get the error message from a Cantera error.
|
||||
%
|
||||
checklib;
|
||||
|
||||
try
|
||||
buflen = calllib(ct, 'ct_getCanteraError', 0, '');
|
||||
aa = char(ones(1, buflen));
|
||||
@ -12,4 +13,5 @@ function e = geterr()
|
||||
catch ME
|
||||
e = getReport(ME);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -1,23 +1,39 @@
|
||||
This experimental Matlab Toolbox for Cantera changes the Matlab interface to the modern Matlab structure and syntaxes for OOP. It replaces the MEX interface with direct function calling from Canter Clib.
|
||||
This experimental Matlab Toolbox for Cantera changes the Matlab interface to the modern
|
||||
Matlab structure and syntaxes for OOP. It replaces the MEX interface with direct
|
||||
function calling from Canter Clib.
|
||||
|
||||
Installation guide:
|
||||
|
||||
1) Install Matlab (any release newer than R2008a).
|
||||
2) Compile Cantera from Source and install in your Conda environment, as directed in this link.
|
||||
https://cantera.org/install/compiling-install.html.
|
||||
2.5) The experimental Matlab Toolbox does not require a SCons option to install at this moment since it's stand-alone. It also does not require the current Matlab Toolbox to be installed.
|
||||
3) For first time users, launch Matlab, then navigate to [/path/to/cantera/source/code] using "Browse for Folder".
|
||||
3.5) Note for Ubuntu users, Matlab must be launched from the terminal using the following command: `LD_PRELOAD='/usr/lib/x86_64-linux-gnu/libstdc++.so.6' matlab`. This is because Matlab does not load the correct GLIBC++ library on start-up and will return an error when loading the Cantera toolbox.
|
||||
4) In the Maltab command window, run `addpath(genpath([pwd, '/interfaces/matlab_experimental']))` to add search path for the experimental toolbox.
|
||||
5) In the Maltab command window, run `addpath(genpath([pwd, '/samples/matlab_experimental']))` to add search path for the sample files.
|
||||
6) In the Matlab command window, run `cd([pwd, '/interfaces/matlab_experimental/Utility'])` to navigate to the Utility folder.
|
||||
7) Open the file named 'cantera_root.m', in the second line, edit `output=` to `output=[/path/to/conda/environment]`, then save the file. This sets the search path for the `LoadCantera` command to find the shared library file for Cantera.
|
||||
8) After steps 3 to 7 are complete, make sure search paths to the current Matlab Toolbox and samples files are removed (if it's already installed). Having both the current and experimental version of the toolbox in the search path will lead to conflicts.
|
||||
9) In the Matlab command window, run `savepath` to save all search paths.
|
||||
10) To start using the experimental toolbox, run `LoadCantera` command.
|
||||
11) To stop using the new Cantera interface, run the following commands:
|
||||
`clear all`
|
||||
`cleanup`
|
||||
`UnloadCantera`
|
||||
11) To switch back to the current matlab toolbox, undo steps 3, 4, 5, 8, and 9. The command to remove search path in Matlab is `rmpath`.
|
||||
12) A future updates will add automated installation.
|
||||
1. Install Matlab (any release newer than R2008a).
|
||||
2. Compile Cantera from Source and install in your Conda environment, as directed in
|
||||
this link. https://cantera.org/install/compiling-install.html. 2.5) The experimental
|
||||
Matlab Toolbox does not require a SCons option to install at this moment since it's
|
||||
stand-alone. It also does not require the current Matlab Toolbox to be installed.
|
||||
3. For first time users, launch Matlab, then navigate to [/path/to/cantera/source/code]
|
||||
using "Browse for Folder". 3.5) Note for Ubuntu users, Matlab must be launched from
|
||||
the terminal using the following command:
|
||||
`LD_PRELOAD='/usr/lib/x86_64-linux-gnu/libstdc++.so.6' matlab`. This is because
|
||||
Matlab does not load the correct GLIBC++ library on start-up and will return an error
|
||||
when loading the Cantera toolbox.
|
||||
4. In the Maltab command window, run
|
||||
`addpath(genpath([pwd, '/interfaces/matlab_experimental']))` to add search path for
|
||||
the experimental toolbox.
|
||||
5. In the Maltab command window, run
|
||||
`addpath(genpath([pwd, '/samples/matlab_experimental']))` to add search path for the
|
||||
sample files.
|
||||
6. In the Matlab command window, run
|
||||
`cd([pwd, '/interfaces/matlab_experimental/Utility'])` to navigate to the Utility
|
||||
folder.
|
||||
7. Open the file named 'cantera_root.m', in the second line, edit `output=` to
|
||||
`output=[/path/to/conda/environment]`, then save the file. This sets the search path
|
||||
for the `LoadCantera` command to find the shared library file for Cantera.
|
||||
8. After steps 3 to 7 are complete, make sure search paths to the current Matlab Toolbox
|
||||
and samples files are removed (if it's already installed). Having both the current
|
||||
and experimental version of the toolbox in the search path will lead to conflicts.
|
||||
9. In the Matlab command window, run `savepath` to save all search paths.
|
||||
10. To start using the experimental toolbox, run `LoadCantera` command.
|
||||
11. To stop using the new Cantera interface, run the following commands: `clear all`
|
||||
`cleanup` `UnloadCantera`
|
||||
12. To switch back to the current matlab toolbox, undo steps 3, 4, 5, 8, and 9. The
|
||||
command to remove search path in Matlab is `rmpath`.
|
||||
13. A future updates will add automated installation.
|
||||
|
||||
@ -8,11 +8,11 @@ function F = PFR_Solver(x, soln_vector, gas, mdot, A_in, dAdx, k)
|
||||
|
||||
rho = soln_vector(1);
|
||||
T = soln_vector(2);
|
||||
Y = soln_vector(3: end);
|
||||
Y = soln_vector(3:end);
|
||||
|
||||
if k==1
|
||||
A = A_in+k * dAdx * x;
|
||||
elseif k==-1
|
||||
if k == 1
|
||||
A = A_in + k * dAdx * x;
|
||||
elseif k == -1
|
||||
A = A_in + k * dAdx * x;
|
||||
dAdx = -dAdx;
|
||||
else
|
||||
@ -31,7 +31,7 @@ function F = PFR_Solver(x, soln_vector, gas, mdot, A_in, dAdx, k)
|
||||
|
||||
gas.basis = 'mass';
|
||||
MW = gas.molecularWeights;
|
||||
h = gas.enthalpies_RT.*R.*T;
|
||||
h = gas.enthalpies_RT .* R .* T;
|
||||
w = gas.netProdRates;
|
||||
Cp = gas.cp;
|
||||
%--------------------------------------------------------------------------
|
||||
@ -39,14 +39,14 @@ function F = PFR_Solver(x, soln_vector, gas, mdot, A_in, dAdx, k)
|
||||
%---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) )...
|
||||
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(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(3:nsp + 2) = w(1:nsp) .* MW(1:nsp) ./ (rho * vx);
|
||||
|
||||
F = F';
|
||||
|
||||
|
||||
@ -20,37 +20,37 @@ clear all
|
||||
close all
|
||||
clc
|
||||
|
||||
t0 = cputime; % record the starting time
|
||||
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
|
||||
p = oneatm; % pressure
|
||||
tinlet = 300.0; % inlet temperature
|
||||
tsurf = 900.0; % surface temperature
|
||||
mdot = 0.06; % kg/m^2/s
|
||||
transport = 'Mix'; % 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';
|
||||
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';
|
||||
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
|
||||
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
|
||||
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)
|
||||
loglevel = 1; % amount of diagnostic output
|
||||
% (0 to 5)
|
||||
|
||||
refine_grid = 1; % 1 to enable refinement, 0 to
|
||||
% disable
|
||||
refine_grid = 1; % 1 to enable refinement, 0 to
|
||||
% disable
|
||||
|
||||
%% Create the gas object
|
||||
%
|
||||
@ -133,15 +133,16 @@ surf.T = tsurf;
|
||||
stack = Sim1D([inlt, flow, surf]);
|
||||
|
||||
% set the initial profiles.
|
||||
stack.setProfile(2, {'velocity', 'spread_rate', 'T'}, [0.0, 1.0 % z/zmax
|
||||
0.06, 0.0 % u
|
||||
0.0, 0.0 % V
|
||||
tinlet, tsurf]); % T
|
||||
stack.setProfile(2, {'velocity', 'spread_rate', '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);
|
||||
stack.setProfile(2, names{k}, [0, 1; y, y]);
|
||||
y = inlt.massFraction(k);
|
||||
stack.setProfile(2, names{k}, [0, 1; y, y]);
|
||||
end
|
||||
|
||||
stack
|
||||
@ -169,11 +170,11 @@ stack.solve(1, refine_grid);
|
||||
|
||||
surf.setCoverageEqs('on');
|
||||
|
||||
for iter=1:6
|
||||
mult = 10.0^(iter - 6);
|
||||
surf_phase.setMultiplier(mult);
|
||||
gas.setMultiplier(mult);
|
||||
stack.solve(1, refine_grid);
|
||||
for iter = 1:6
|
||||
mult = 10.0^(iter - 6);
|
||||
surf_phase.setMultiplier(mult);
|
||||
gas.setMultiplier(mult);
|
||||
stack.solve(1, refine_grid);
|
||||
end
|
||||
|
||||
% At this point, we should have the solution for the hydrogen/air
|
||||
@ -196,7 +197,7 @@ stack.saveSoln('catcomb.xml', 'energy', ['solution with energy equation']);
|
||||
|
||||
stack.writeStats;
|
||||
elapsed = cputime - t0;
|
||||
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
|
||||
e = sprintf('Elapsed CPU time: %10.4g', elapsed);
|
||||
disp(e);
|
||||
|
||||
%% Make plots
|
||||
|
||||
@ -6,7 +6,7 @@ function dydt = conhp(t, y, gas, mw)
|
||||
% 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.Y = y(2:end);
|
||||
gas.TP = {y(1), gas.P};
|
||||
nsp = gas.nSpecies;
|
||||
|
||||
@ -14,15 +14,17 @@ function dydt = conhp(t, y, gas, mw)
|
||||
wdot = gas.netProdRates;
|
||||
H = gas.enthalpies_RT';
|
||||
gas.basis = 'mass';
|
||||
tdot = - gas.T * gasconstant /(gas.D * gas.cp).* wdot * H;
|
||||
tdot =- gas.T * gasconstant / (gas.D * gas.cp) .* wdot * H;
|
||||
|
||||
% set up column vector for dydt
|
||||
dydt = [tdot
|
||||
zeros(nsp, 1)];
|
||||
|
||||
% species equations
|
||||
rrho = 1.0/gas.D;
|
||||
rrho = 1.0 / gas.D;
|
||||
|
||||
for i = 1:nsp
|
||||
dydt(i+1) = rrho * mw(i) * wdot(i);
|
||||
dydt(i + 1) = rrho * mw(i) * wdot(i);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -5,7 +5,6 @@ function dydt = conuv(t, y, gas, mw)
|
||||
% 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};
|
||||
@ -15,16 +14,17 @@ function dydt = conuv(t, y, gas, mw)
|
||||
wdot = gas.netProdRates;
|
||||
H = (gas.enthalpies_RT - ones(1, nsp))';
|
||||
gas.basis = 'mass';
|
||||
tdot = - gas.T * gasconstant / (gas.D * gas.cv).* wdot * H;
|
||||
tdot =- gas.T * gasconstant / (gas.D * gas.cv) .* wdot * H;
|
||||
|
||||
% set up column vector for dydt
|
||||
dydt = [tdot
|
||||
zeros(nsp, 1)];
|
||||
|
||||
% species equations
|
||||
rrho = 1.0/gas.D;
|
||||
rrho = 1.0 / gas.D;
|
||||
|
||||
for i = 1:nsp
|
||||
dydt(i+1) = rrho * mw(i) * wdot(i);
|
||||
dydt(i + 1) = rrho * mw(i) * wdot(i);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -24,8 +24,8 @@ t0 = cputime; % record the starting time
|
||||
|
||||
%% Operation Parameters
|
||||
|
||||
t = 1200.0; % surface temperature
|
||||
p = 20.0 * oneatm / 760.0; % pressure
|
||||
t = 1200.0; % surface temperature
|
||||
p = 20.0 * oneatm / 760.0; % pressure
|
||||
|
||||
%% Create the gas object
|
||||
%
|
||||
|
||||
@ -15,25 +15,25 @@ clc
|
||||
tic % total running time of the script
|
||||
help diff_flame
|
||||
|
||||
runtime = cputime; % Record the starting time
|
||||
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
|
||||
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
|
||||
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
|
||||
%
|
||||
@ -42,8 +42,8 @@ refine_grid = 1; % 1 to enable refinement, 0 to disable
|
||||
|
||||
fuel = Solution('gri30.yaml', 'gri30', transport);
|
||||
ox = Solution('gri30.yaml', 'gri30', transport);
|
||||
oxcomp = 'O2:0.21, N2:0.78'; % Air composition
|
||||
fuelcomp = 'C2H6:1'; % Fuel composition
|
||||
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};
|
||||
@ -112,26 +112,26 @@ fl.saveSoln('c2h6.xml', 'energy', ['solution with energy equation']);
|
||||
|
||||
fl.writeStats;
|
||||
elapsed = cputime - runtime;
|
||||
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
|
||||
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
|
||||
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
|
||||
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
|
||||
|
||||
@ -1,10 +1,10 @@
|
||||
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
|
||||
%% 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
|
||||
|
||||
clear all
|
||||
close all
|
||||
@ -14,17 +14,17 @@ function equil(g)
|
||||
help equil
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
gas = g;
|
||||
else
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
end
|
||||
|
||||
nsp = gas.nSpecies;
|
||||
|
||||
% find methane, nitrogen, and oxygen indices
|
||||
ich4 = gas.speciesIndex('CH4');
|
||||
io2 = gas.speciesIndex('O2');
|
||||
in2 = gas.speciesIndex('N2');
|
||||
io2 = gas.speciesIndex('O2');
|
||||
in2 = gas.speciesIndex('N2');
|
||||
|
||||
nPhis = 50;
|
||||
phi = linspace(0.2, 2.70, nPhis);
|
||||
@ -32,16 +32,16 @@ function equil(g)
|
||||
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;
|
||||
gas.T = 300;
|
||||
gas.P = 101325;
|
||||
gas.X = x;
|
||||
gas.equilibrate('HP');
|
||||
tad(i) = gas.T;
|
||||
xeq(:, i) = gas.X;
|
||||
x = zeros(nsp, 1);
|
||||
x(ich4, 1) = phi(i);
|
||||
x(io2, 1) = 2.0;
|
||||
x(in2, 1) = 7.52;
|
||||
gas.T = 300;
|
||||
gas.P = 101325;
|
||||
gas.X = x;
|
||||
gas.equilibrate('HP');
|
||||
tad(i) = gas.T;
|
||||
xeq(:, i) = gas.X;
|
||||
end
|
||||
|
||||
% make plots
|
||||
@ -52,18 +52,22 @@ function equil(g)
|
||||
ylabel('Temperature (K)');
|
||||
title('Adiabatic Flame Temperature');
|
||||
|
||||
subplot(1,2,2);
|
||||
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), gas.speciesName(k))
|
||||
text(phi(j), 1.5 * xeq(k, j), gas.speciesName(k))
|
||||
j = j + 2;
|
||||
|
||||
if j > 46
|
||||
j = 10;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
xlabel('Equivalence Ratio');
|
||||
ylabel('Mole Fraction');
|
||||
title('Equilibrium Composition');
|
||||
|
||||
@ -1,26 +1,29 @@
|
||||
function f = flame(gas, left, flow, right)
|
||||
% FLAME - create a flame object.
|
||||
%
|
||||
% FLAME - create a flame object.
|
||||
%
|
||||
% Check input parameters
|
||||
if nargin ~= 4
|
||||
error('wrong number of input arguments.');
|
||||
error('wrong number of input arguments.');
|
||||
end
|
||||
|
||||
if ~gas.isIdealGas
|
||||
error('gas object must represent an ideal gas mixture.');
|
||||
error('gas object must represent an ideal gas mixture.');
|
||||
end
|
||||
|
||||
if ~left.isInlet
|
||||
error('burner object of wrong type.');
|
||||
error('burner object of wrong type.');
|
||||
end
|
||||
|
||||
if ~flow.isFlow
|
||||
error('flow object of wrong type.');
|
||||
error('flow object of wrong type.');
|
||||
end
|
||||
|
||||
flametype = 0;
|
||||
|
||||
if right.isSurface
|
||||
flametype = 1;
|
||||
flametype = 1;
|
||||
elseif right.isInlet
|
||||
flametype = 3;
|
||||
flametype = 3;
|
||||
end
|
||||
|
||||
% create the container object
|
||||
@ -39,36 +42,40 @@ function f = flame(gas, left, flow, right)
|
||||
mdot0 = left.massFlux;
|
||||
mdot1 = right.massFlux;
|
||||
t0 = left.T;
|
||||
|
||||
if flametype == 0
|
||||
t1 = teq;
|
||||
mdot1 = -mdot0;
|
||||
t1 = teq;
|
||||
mdot1 = -mdot0;
|
||||
else
|
||||
t1 = right.T;
|
||||
t1 = right.T;
|
||||
end
|
||||
f.setProfile(2, {'velocity', 'spread_rate'}, [0.0 1.0
|
||||
mdot0/rho0 -mdot1/rho0
|
||||
0.0 0.0]);
|
||||
|
||||
f.setProfile(2, {'velocity', 'spread_rate'}, [0.0 1.0
|
||||
mdot0 / rho0 -mdot1 / rho0
|
||||
0.0 0.0]);
|
||||
f.setProfile(2, 'T', [0.0, 1.0
|
||||
t0, t1]);
|
||||
|
||||
for n = 1:gas.nSpecies
|
||||
nm = gas.speciesName(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 = right.massFraction(n);
|
||||
else
|
||||
y1 = yeq(n);
|
||||
end
|
||||
f.setProfile(2, nm, [0, 1.0
|
||||
left.massFraction(n), y1]);
|
||||
nm = gas.speciesName(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 = right.massFraction(n);
|
||||
else
|
||||
y1 = yeq(n);
|
||||
end
|
||||
|
||||
f.setProfile(2, nm, [0, 1.0
|
||||
left.massFraction(n), y1]);
|
||||
end
|
||||
|
||||
% set minimal grid refinement criteria
|
||||
f.setRefineCriteria(2, 10.0, 0.8, 0.8);
|
||||
end
|
||||
|
||||
|
||||
@ -14,29 +14,28 @@ clc
|
||||
tic
|
||||
help flame1
|
||||
|
||||
t0 = cputime; % record the starting time
|
||||
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
|
||||
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
|
||||
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
|
||||
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
|
||||
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)
|
||||
loglevel = 1; % amount of diagnostic output (0
|
||||
% to 5)
|
||||
|
||||
refine_grid = 1; % 1 to enable refinement, 0 to
|
||||
% disable
|
||||
refine_grid = 1; % 1 to enable refinement, 0 to
|
||||
% disable
|
||||
max_jacobian_age = [5, 10];
|
||||
|
||||
%% Create the gas object
|
||||
@ -81,7 +80,7 @@ s = Outlet('out');
|
||||
% to create the flame object.
|
||||
%
|
||||
fl = flame(gas, burner, f, s);
|
||||
fl.setMaxJacAge(max_jacobian_age(1), max_jacobian_age(2));
|
||||
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.
|
||||
@ -104,7 +103,7 @@ fl.saveSoln('h2fl.xml', 'energy', ['solution with energy equation']);
|
||||
|
||||
fl.writeStats;
|
||||
elapsed = cputime - t0;
|
||||
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
|
||||
e = sprintf('Elapsed CPU time: %10.4g', elapsed);
|
||||
disp(e);
|
||||
|
||||
%% Make plots
|
||||
|
||||
@ -13,30 +13,28 @@ clc
|
||||
tic
|
||||
help flame2
|
||||
|
||||
t0 = cputime; % record the starting time
|
||||
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
|
||||
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
|
||||
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
|
||||
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
|
||||
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)
|
||||
loglevel = 1; % amount of diagnostic output (0 to 5)
|
||||
|
||||
refine_grid = 1; % 1 to enable refinement, 0 to
|
||||
% disable
|
||||
refine_grid = 1; % 1 to enable refinement, 0 to disable
|
||||
|
||||
%% Create the gas object
|
||||
%
|
||||
@ -50,7 +48,7 @@ gas.TPX = {tin, p, comp2};
|
||||
|
||||
%% Create the flow object
|
||||
%
|
||||
f = AxisymmetricFlow(gas,'flow');
|
||||
f = AxisymmetricFlow(gas, 'flow');
|
||||
f.P = p;
|
||||
f.setupGrid(initial_grid);
|
||||
f.setSteadyTolerances('default', tol_ss{:});
|
||||
@ -102,7 +100,7 @@ fl.saveSoln('c2h6.xml', 'energy', ['solution with energy equation']);
|
||||
|
||||
fl.writeStats;
|
||||
elapsed = cputime - t0;
|
||||
e = sprintf('Elapsed CPU time: %10.4g',elapsed);
|
||||
e = sprintf('Elapsed CPU time: %10.4g', elapsed);
|
||||
disp(e);
|
||||
|
||||
%% Make plots
|
||||
|
||||
@ -1,11 +1,11 @@
|
||||
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
|
||||
%% 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
|
||||
|
||||
clear all
|
||||
close all
|
||||
@ -15,9 +15,9 @@ function plotdata = ignite(g)
|
||||
help ignite
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
gas = g;
|
||||
else
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
end
|
||||
|
||||
% set the initial conditions
|
||||
@ -25,7 +25,7 @@ function plotdata = ignite(g)
|
||||
gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'};
|
||||
gas.basis = 'mass';
|
||||
y0 = [gas.U
|
||||
1.0/gas.D
|
||||
1.0 / gas.D
|
||||
gas.Y'];
|
||||
|
||||
time_interval = [0 0.001];
|
||||
@ -33,7 +33,7 @@ function plotdata = ignite(g)
|
||||
|
||||
t0 = cputime;
|
||||
out = ode15s(@reactor_ode, time_interval, y0, options, gas, ...
|
||||
@vdot, @area, @heatflux);
|
||||
@vdot, @area, @heatflux);
|
||||
disp(['CPU time = ' num2str(cputime - t0)]);
|
||||
|
||||
plotdata = output(out, gas);
|
||||
@ -50,18 +50,21 @@ function plotdata = ignite(g)
|
||||
% gas ideal gas object
|
||||
|
||||
function v = vdot(t, vol, gas)
|
||||
%v = 0.0; %uncomment for constant volume
|
||||
v = 1.e11 * (gas.P - 101325.0); % holds pressure very
|
||||
% close to 1 atm
|
||||
%v = 0.0; %uncomment for constant volume
|
||||
v = 1.e11 * (gas.P - 101325.0); % holds pressure very
|
||||
% close to 1 atm
|
||||
end
|
||||
|
||||
% heat flux (W/m^2).
|
||||
function q = heatflux(t, gas)
|
||||
q = 0.0; % adiabatic
|
||||
q = 0.0; % adiabatic
|
||||
end
|
||||
|
||||
% surface area (m^2). Used only to compute heat transfer.
|
||||
function a = area(t, vol)
|
||||
a = 1.0;
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
% Since the solution variables used by the 'reactor' function are
|
||||
@ -72,25 +75,25 @@ function plotdata = ignite(g)
|
||||
function pv = output(s, gas)
|
||||
times = s.x;
|
||||
soln = s.y;
|
||||
[~,n] = size(times);
|
||||
[~, n] = size(times);
|
||||
pv = zeros(gas.nSpecies + 4, n);
|
||||
|
||||
gas.TP = {1001.0, 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;
|
||||
gas.Y = y;
|
||||
gas.UV = {u_mass, v_mass};
|
||||
ss = soln(:, j);
|
||||
y = ss(3:end);
|
||||
mass = sum(y);
|
||||
u_mass = ss(1) / mass;
|
||||
v_mass = ss(2) / mass;
|
||||
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.
|
||||
@ -102,7 +105,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');
|
||||
|
||||
@ -12,7 +12,7 @@ function ignite_hp(gas)
|
||||
help ignite_hp
|
||||
|
||||
if nargin == 0
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
end
|
||||
|
||||
mw = gas.molecularWeights;
|
||||
@ -27,19 +27,19 @@ function ignite_hp(gas)
|
||||
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']);
|
||||
% 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');
|
||||
figure(2);
|
||||
ioh = gas.speciesIndex('OH');
|
||||
plot(out.x, out.y(1 + ioh, :));
|
||||
xlabel('time');
|
||||
ylabel('Mass Fraction');
|
||||
title('OH Mass Fraction');
|
||||
end
|
||||
|
||||
toc
|
||||
|
||||
@ -12,7 +12,7 @@ function ignite_uv(gas)
|
||||
help ignite_uv
|
||||
|
||||
if nargin == 0
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
end
|
||||
|
||||
mw = gas.molecularWeights;
|
||||
@ -27,19 +27,19 @@ function ignite_uv(gas)
|
||||
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']);
|
||||
% 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');
|
||||
figure(2);
|
||||
ioh = gas.speciesIndex('OH');
|
||||
plot(out.x, out.y(1 + ioh, :));
|
||||
xlabel('time');
|
||||
ylabel('Mass Fraction');
|
||||
title('OH Mass Fraction');
|
||||
end
|
||||
|
||||
toc
|
||||
|
||||
@ -1,10 +1,10 @@
|
||||
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
|
||||
%% 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
|
||||
|
||||
clear all
|
||||
close all
|
||||
@ -14,9 +14,9 @@ function isentropic(g)
|
||||
help isentropic
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
gas = g;
|
||||
else
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
gas = Solution('gri30.yaml', 'gri30');
|
||||
end
|
||||
|
||||
% set the stagnation state
|
||||
@ -26,7 +26,7 @@ function isentropic(g)
|
||||
h0 = gas.H;
|
||||
p0 = gas.P;
|
||||
|
||||
mdot = 1; % arbitrary
|
||||
mdot = 1; % arbitrary
|
||||
|
||||
mach = [];
|
||||
a = [];
|
||||
@ -35,31 +35,32 @@ function isentropic(g)
|
||||
|
||||
% compute values for a range of pressure ratios
|
||||
for r = 0.005:0.0025:0.995
|
||||
p = p0*r;
|
||||
p = p0 * r;
|
||||
|
||||
% set the state using (p,s0)
|
||||
gas.SP = {s0, p};
|
||||
% set the state using (p,s0)
|
||||
gas.SP = {s0, p};
|
||||
|
||||
h = gas.H;
|
||||
rho = gas.D;
|
||||
h = gas.H;
|
||||
rho = gas.D;
|
||||
|
||||
v2 = 2.0*(h0 - h); % h + V^2/2 = h0
|
||||
v = sqrt(v2);
|
||||
a(i) = mdot/(rho*v); % rho*v*A = constant
|
||||
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/gas.soundSpeed;
|
||||
i = i + 1;
|
||||
if a(i) < amin
|
||||
amin = a(i);
|
||||
end
|
||||
|
||||
mach(i) = v / gas.soundSpeed;
|
||||
i = i + 1;
|
||||
end
|
||||
|
||||
a = a/amin;
|
||||
a = a / amin;
|
||||
|
||||
% plot results
|
||||
|
||||
clf;
|
||||
plot(mach,a);
|
||||
plot(mach, a);
|
||||
ylabel('Area Ratio');
|
||||
xlabel('Mach Number');
|
||||
title('Isentropic Flow: Area Ratio vs. Mach Number');
|
||||
|
||||
@ -55,8 +55,8 @@ 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
|
||||
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');
|
||||
@ -75,24 +75,25 @@ anode_interface.TP = {T, P};
|
||||
cathode_interface.TP = {T, P};
|
||||
|
||||
% Calculate cell voltage, separately for each entry of the input vectors
|
||||
V_cell = zeros(length(SOC),1);
|
||||
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);
|
||||
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;
|
||||
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);
|
||||
cathode, elde, elyt, cathode_interface, S_ca) ...
|
||||
- I_app, phi_s_ca);
|
||||
|
||||
% Calculate cell voltage
|
||||
V_cell(i) = phi_s_ca - phi_s_an;
|
||||
@ -100,7 +101,7 @@ 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)
|
||||
plot(SOC * 100, V_cell, 'linewidth', 2.5)
|
||||
ylim([2.5, 4.3])
|
||||
xlabel('State of charge / %')
|
||||
ylabel('Cell voltage / V')
|
||||
@ -111,9 +112,10 @@ toc
|
||||
%% 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)
|
||||
function anCurr = anode_curr(phi_s, phi_l, X_Li_an, anode, elde, elyt, ...
|
||||
anode_interface, S_an)
|
||||
% Set the active material mole fraction
|
||||
anode.X = ['Li[anode]:' num2str(X_Li_an) ', V[anode]:' num2str(1-X_Li_an)];
|
||||
anode.X = ['Li[anode]:' num2str(X_Li_an) ', V[anode]:' num2str(1 - X_Li_an)];
|
||||
|
||||
% Set the electrode and electrolyte potential
|
||||
elde.setElectricPotential(phi_s);
|
||||
@ -124,13 +126,13 @@ function anCurr = anode_curr(phi_s, phi_l, X_Li_an, anode, elde, elyt, anode_int
|
||||
r = anode_interface.ropNet; % [kmol/m2/s]
|
||||
|
||||
% Calculate the current. Should be negative for cell discharge.
|
||||
anCurr = r*faradayconstant*S_an; %
|
||||
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
|
||||
cathode.X = ['Li[cathode]:' num2str(X_Li_ca) ', V[cathode]:' num2str(1-X_Li_ca)];
|
||||
cathode.X = ['Li[cathode]:' num2str(X_Li_ca) ', V[cathode]:' num2str(1 - X_Li_ca)];
|
||||
|
||||
% Set the electrode and electrolyte potential
|
||||
elde.setElectricPotential(phi_s);
|
||||
@ -141,5 +143,5 @@ function caCurr = cathode_curr(phi_s, phi_l, X_Li_ca, cathode, elde, elyt, catho
|
||||
r = cathode_interface.ropNet; % [kmol/m2/s]
|
||||
|
||||
% Calculate the current. Should be negative for cell discharge.
|
||||
caCurr = r*faradayconstant*S_ca*(-1); %
|
||||
caCurr = r * faradayconstant * S_ca * (-1); %
|
||||
end
|
||||
|
||||
@ -1,23 +1,23 @@
|
||||
function periodic_cstr
|
||||
%% 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
|
||||
%% 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
|
||||
|
||||
clear all
|
||||
close all
|
||||
@ -27,7 +27,7 @@ function periodic_cstr
|
||||
help periodic_cstr
|
||||
|
||||
% create the gas mixture
|
||||
gas = Solution('h2o2.yaml','ohmech');
|
||||
gas = Solution('h2o2.yaml', 'ohmech');
|
||||
|
||||
% pressure = 60 Torr, T = 770 K
|
||||
p = 60.0 * 133.3;
|
||||
@ -48,7 +48,7 @@ function periodic_cstr
|
||||
|
||||
% 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.
|
||||
cstr.setInitialVolume(10.0*1.0e-6);
|
||||
cstr.setInitialVolume(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
|
||||
@ -68,8 +68,8 @@ function periodic_cstr
|
||||
% 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 / gas.P) * ( gas.T / 273.15)); % m^3/s
|
||||
mdot = gas.D * vdot; % kg/s
|
||||
vdot = sccm * 1.0e-6/60.0 * ((oneatm / gas.P) * (gas.T / 273.15)); % m^3/s
|
||||
mdot = gas.D * vdot; % kg/s
|
||||
mfc = MassFlowController;
|
||||
mfc.install(upstream, cstr);
|
||||
mfc.setMassFlowRate(mdot);
|
||||
@ -89,22 +89,24 @@ function periodic_cstr
|
||||
|
||||
% now integrate in time
|
||||
tme = 0.0;
|
||||
dt = 0.1;
|
||||
dt = 0.1;
|
||||
|
||||
n = 0;
|
||||
|
||||
while tme < 300.0
|
||||
n = n + 1;
|
||||
tme = tme + dt;
|
||||
network.advance(tme);
|
||||
tm(n) = tme;
|
||||
y(1,n) = cstr.massFraction('H2');
|
||||
y(2,n) = cstr.massFraction('O2');
|
||||
y(3,n) = cstr.massFraction('H2O');
|
||||
y(1, n) = cstr.massFraction('H2');
|
||||
y(2, n) = cstr.massFraction('O2');
|
||||
y(3, n) = cstr.massFraction('H2O');
|
||||
end
|
||||
|
||||
clf
|
||||
figure(1)
|
||||
plot(tm,y)
|
||||
legend('H2','O2','H2O')
|
||||
plot(tm, y)
|
||||
legend('H2', 'O2', 'H2O')
|
||||
title('Mass Fractions')
|
||||
|
||||
toc
|
||||
|
||||
@ -29,14 +29,14 @@ close all
|
||||
clc
|
||||
|
||||
tic
|
||||
help PFR
|
||||
help plug_flow_reactor
|
||||
|
||||
%% Operation Parameters
|
||||
|
||||
% Temperature of gas, in K
|
||||
T0 = 1473;
|
||||
% Pressure of gas, in Pa
|
||||
P0 = 4.47*101325;
|
||||
P0 = 4.47 * 101325;
|
||||
|
||||
% Equivalence Ratio
|
||||
Phi = 0.2899;
|
||||
@ -44,15 +44,15 @@ Phi = 0.2899;
|
||||
% Import the gas phase, read out key species indices:
|
||||
gas_calc = Solution('gri30.yaml', 'gri30');
|
||||
ich4 = gas_calc.speciesIndex('CH4');
|
||||
io2 = gas_calc.speciesIndex('O2');
|
||||
in2 = gas_calc.speciesIndex('N2');
|
||||
io2 = gas_calc.speciesIndex('O2');
|
||||
in2 = gas_calc.speciesIndex('N2');
|
||||
nsp = gas_calc.nSpecies;
|
||||
x = zeros(nsp,1);
|
||||
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;
|
||||
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:
|
||||
gas_calc.TPX = {T0, P0, x};
|
||||
@ -67,7 +67,7 @@ A_in = 0.018;
|
||||
% Exit Area, in m^2
|
||||
A_out = 0.003;
|
||||
% Length of the reactor, in m
|
||||
L = 1.284*0.0254;
|
||||
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
|
||||
@ -85,9 +85,9 @@ else
|
||||
k = 0;
|
||||
end
|
||||
|
||||
dAdx = abs(A_in-A_out)/L;
|
||||
dAdx = abs(A_in - A_out) / L;
|
||||
% The whole length of the reactor is divided into n small lengths
|
||||
dx = L/n;
|
||||
dx = L / n;
|
||||
|
||||
x_calc = 0:dx:L;
|
||||
nsp = gas_calc.nSpecies;
|
||||
@ -98,7 +98,7 @@ Y_calc = zeros(length(x_calc), nsp);
|
||||
rho_calc = zeros(length(x_calc), 1);
|
||||
|
||||
T_calc(1) = gas_calc.T;
|
||||
Y_calc(1,:) = gas_calc.Y;
|
||||
Y_calc(1, :) = gas_calc.Y;
|
||||
rho_calc(1) = gas_calc.D;
|
||||
|
||||
for i = 2:length(x_calc)
|
||||
@ -106,39 +106,40 @@ 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)];
|
||||
gas_calc.TDY = {T_calc(i-1), rho_calc(i-1), Y_calc(i-1, :)};
|
||||
options = odeset('RelTol', 1.e-10, 'AbsTol', 1e-10,...
|
||||
'InitialStep', 1e-8, 'NonNegative', 1);
|
||||
%--------------------------------------------------------------------------
|
||||
%--------------------------------------------------------------------------
|
||||
%--------------------------------------------------------------------------
|
||||
%------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)];
|
||||
gas_calc.TDY = {T_calc(i - 1), rho_calc(i - 1), 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, ...
|
||||
[~, 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);
|
||||
Y_calc(i, :) = y(end, 3:nsp + 2);
|
||||
end
|
||||
|
||||
A_calc = A_in + k.* x_calc * dAdx;
|
||||
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)
|
||||
|
||||
for i = 1:length(x_calc)
|
||||
% The gas is set to the solved property values at each location
|
||||
gas.TDY = {T_calc(i), rho_calc(i), Y_calc(i, :)};
|
||||
% Velocity is calculated from Mass flow rate, Area and Density
|
||||
vx_calc(i) = mdot_calc./ (A_calc(i) * rho_calc(i));
|
||||
vx_calc(i) = mdot_calc ./ (A_calc(i) * rho_calc(i));
|
||||
% Specific Gas Constant
|
||||
R_calc(i) = gasconstant() / gas_calc.meanMolecularWeight;
|
||||
% Mach No. is calculated from local velocity and local speed of sound
|
||||
|
||||
@ -1,11 +1,11 @@
|
||||
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
|
||||
%% 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
|
||||
|
||||
clear all
|
||||
close all
|
||||
@ -15,15 +15,15 @@ function prandtl1(g)
|
||||
help prandtl1
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
gas = g;
|
||||
else
|
||||
gas = Solution('gri30.yaml', 'gri30', 'Mix');
|
||||
gas = Solution('gri30.yaml', 'gri30', 'Mix');
|
||||
end
|
||||
|
||||
pr = zeros(31,31);
|
||||
xh2 = zeros(31,31);
|
||||
visc = zeros(31,31);
|
||||
lambda = zeros(31,31);
|
||||
pr = zeros(31, 31);
|
||||
xh2 = zeros(31, 31);
|
||||
visc = zeros(31, 31);
|
||||
lambda = zeros(31, 31);
|
||||
t = [];
|
||||
xo2 = [];
|
||||
io2 = gas.speciesIndex('O2');
|
||||
@ -31,52 +31,56 @@ function prandtl1(g)
|
||||
|
||||
minT = gas.minTemp;
|
||||
maxT = gas.maxTemp;
|
||||
dT = (maxT - minT)/30.0;
|
||||
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(gas.nSpecies, 1);
|
||||
x(io2) = xo2(j);
|
||||
x(ih2) = 1.0 - xo2(j);
|
||||
gas.TPX = {t(i), oneatm, x};
|
||||
equilibrate(gas, 'TP');
|
||||
visc(i,j) = gas.viscosity;
|
||||
lambda(i,j) = gas.thermalConductivity;
|
||||
gas.basis = 'mass';
|
||||
pr(i,j) = visc(i,j) * gas.cp / lambda(i,j);
|
||||
x = gas.X;
|
||||
xh2(i,j) = x(ih2);
|
||||
end
|
||||
t(i) = minT + dT * (i - 1);
|
||||
|
||||
for j = 1:31
|
||||
xo2(j) = 0.99 * (j - 1) / 30.0;
|
||||
x = zeros(gas.nSpecies, 1);
|
||||
x(io2) = xo2(j);
|
||||
x(ih2) = 1.0 - xo2(j);
|
||||
gas.TPX = {t(i), oneatm, x};
|
||||
equilibrate(gas, 'TP');
|
||||
visc(i, j) = gas.viscosity;
|
||||
lambda(i, j) = gas.thermalConductivity;
|
||||
gas.basis = 'mass';
|
||||
pr(i, j) = visc(i, j) * gas.cp / lambda(i, j);
|
||||
x = gas.X;
|
||||
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);
|
||||
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);
|
||||
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);
|
||||
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);
|
||||
subplot(2, 2, 4);
|
||||
surf(xo2, t, lambda);
|
||||
xlabel('Elemental O/(O+H)');
|
||||
ylabel('Temperature (K)');
|
||||
zlabel('Thermal Conductivity');
|
||||
|
||||
@ -14,15 +14,15 @@ function prandtl2(g)
|
||||
help prandtl2
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
gas = g;
|
||||
else
|
||||
gas = Solution('gri30.yaml', 'gri30', 'Multi');
|
||||
gas = Solution('gri30.yaml', 'gri30', 'Multi');
|
||||
end
|
||||
|
||||
pr = zeros(31,31);
|
||||
xh2 = zeros(31,31);
|
||||
visc = zeros(31,31);
|
||||
lambda = zeros(31,31);
|
||||
pr = zeros(31, 31);
|
||||
xh2 = zeros(31, 31);
|
||||
visc = zeros(31, 31);
|
||||
lambda = zeros(31, 31);
|
||||
t = [];
|
||||
xo2 = [];
|
||||
io2 = gas.speciesIndex('O2');
|
||||
@ -30,51 +30,55 @@ function prandtl2(g)
|
||||
|
||||
minT = gas.minTemp;
|
||||
maxT = gas.maxTemp;
|
||||
dT = (maxT - minT)/30.0;
|
||||
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(gas.nSpecies,1);
|
||||
x(io2) = xo2(j);
|
||||
x(ih2) = 1.0 - xo2(j);
|
||||
gas.TPX = {t(i), oneatm, x};
|
||||
equilibrate(gas, 'TP');
|
||||
visc(i,j) = gas.viscosity;
|
||||
lambda(i,j) = gas.thermalConductivity;
|
||||
gas.basis = 'mass';
|
||||
pr(i,j) = visc(i,j) * gas.cp / lambda(i,j);
|
||||
x = gas.X;
|
||||
xh2(i,j) = x(ih2);
|
||||
end
|
||||
t(i) = minT + dT * (i - 1);
|
||||
|
||||
for j = 1:31
|
||||
xo2(j) = 0.99 * (j - 1) / 30.0;
|
||||
x = zeros(gas.nSpecies, 1);
|
||||
x(io2) = xo2(j);
|
||||
x(ih2) = 1.0 - xo2(j);
|
||||
gas.TPX = {t(i), oneatm, x};
|
||||
equilibrate(gas, 'TP');
|
||||
visc(i, j) = gas.viscosity;
|
||||
lambda(i, j) = gas.thermalConductivity;
|
||||
gas.basis = 'mass';
|
||||
pr(i, j) = visc(i, j) * gas.cp / lambda(i, j);
|
||||
x = gas.X;
|
||||
xh2(i, j) = x(ih2);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
disp(['CPU time = ' num2str(cputime - t0)]);
|
||||
|
||||
% plot results
|
||||
|
||||
clf;
|
||||
subplot(2,2,1);
|
||||
surf(xo2,t,pr);
|
||||
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);
|
||||
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);
|
||||
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);
|
||||
subplot(2, 2, 4);
|
||||
surf(xo2, t, lambda);
|
||||
xlabel('Elemental O/(O+H)');
|
||||
ylabel('Temperature (K)');
|
||||
zlabel('Thermal Conductivity');
|
||||
|
||||
@ -40,7 +40,7 @@ turbine_work = expand(w, p1, eta_turbine);
|
||||
w
|
||||
|
||||
% compute the efficiency
|
||||
efficiency = (turbine_work - pump_work)/heat_added;
|
||||
efficiency = (turbine_work - pump_work) / heat_added;
|
||||
disp(sprintf('efficiency = %d', efficiency));
|
||||
|
||||
function work = pump(fluid, pfinal, eta)
|
||||
@ -59,7 +59,6 @@ function work = pump(fluid, pfinal, eta)
|
||||
work = actual_work;
|
||||
end
|
||||
|
||||
|
||||
function work = expand(fluid, pfinal, eta)
|
||||
% EXPAND - Adiabatically expand a fluid to pressure pfinal, using a
|
||||
% turbine with isentropic efficiency eta.
|
||||
|
||||
@ -1,11 +1,11 @@
|
||||
function reactor1(g)
|
||||
%% REACTOR1 Zero-dimensional kinetics: adiabatc, constant pressure.
|
||||
%
|
||||
% 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
|
||||
%% REACTOR1 Zero-dimensional kinetics: adiabatc, constant pressure.
|
||||
%
|
||||
% 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
|
||||
|
||||
clear all
|
||||
close all
|
||||
@ -15,9 +15,9 @@ function reactor1(g)
|
||||
help reactor1
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
gas = g;
|
||||
else
|
||||
gas = Solution('gri30.yaml', 'gri30', 'None');
|
||||
gas = Solution('gri30.yaml', 'gri30', 'None');
|
||||
end
|
||||
|
||||
P = 101325.0;
|
||||
@ -35,7 +35,7 @@ function reactor1(g)
|
||||
r = IdealGasReactor(gas);
|
||||
|
||||
% create a reservoir to represent the environment
|
||||
a = Solution('air.yaml','air','None');
|
||||
a = Solution('air.yaml', 'air', 'None');
|
||||
a.P = P;
|
||||
env = Reservoir(a);
|
||||
|
||||
@ -43,7 +43,7 @@ function reactor1(g)
|
||||
% make it flexible, so that the pressure in the reactor is held
|
||||
% at the environment pressure.
|
||||
w = Wall;
|
||||
w.install(r,env);
|
||||
w.install(r, env);
|
||||
|
||||
% set expansion parameter. dV/dt = KA(P_1 - P_2)
|
||||
w.setExpansionRateCoeff(1.0e6);
|
||||
@ -57,34 +57,36 @@ function reactor1(g)
|
||||
nSteps = 100;
|
||||
tim(nSteps) = 0;
|
||||
temp(nSteps) = 0;
|
||||
x(nSteps,3) = 0;
|
||||
x(nSteps, 3) = 0;
|
||||
t = 0.0;
|
||||
dt = 1.0e-5;
|
||||
t0 = cputime;
|
||||
|
||||
for n = 1:nSteps
|
||||
t = t + dt;
|
||||
network.advance(t);
|
||||
tim(n) = network.time;
|
||||
temp(n) = r.T;
|
||||
x(n,1:3) = gas.moleFraction({'OH','H','H2'});
|
||||
t = t + dt;
|
||||
network.advance(t);
|
||||
tim(n) = network.time;
|
||||
temp(n) = r.T;
|
||||
x(n, 1:3) = gas.moleFraction({'OH', 'H', 'H2'});
|
||||
end
|
||||
|
||||
disp(['CPU time = ' num2str(cputime - t0)]);
|
||||
|
||||
clf;
|
||||
subplot(2,2,1);
|
||||
plot(tim,temp);
|
||||
subplot(2, 2, 1);
|
||||
plot(tim, temp);
|
||||
xlabel('Time (s)');
|
||||
ylabel('Temperature (K)');
|
||||
subplot(2,2,2)
|
||||
plot(tim,x(:,1));
|
||||
subplot(2, 2, 2)
|
||||
plot(tim, x(:, 1));
|
||||
xlabel('Time (s)');
|
||||
ylabel('OH Mole Fraction (K)');
|
||||
subplot(2,2,3)
|
||||
plot(tim,x(:,2));
|
||||
subplot(2, 2, 3)
|
||||
plot(tim, x(:, 2));
|
||||
xlabel('Time (s)');
|
||||
ylabel('H Mole Fraction (K)');
|
||||
subplot(2,2,4)
|
||||
plot(tim,x(:,3));
|
||||
subplot(2, 2, 4)
|
||||
plot(tim, x(:, 3));
|
||||
xlabel('Time (s)');
|
||||
ylabel('H2 Mole Fraction (K)');
|
||||
|
||||
|
||||
@ -1,11 +1,11 @@
|
||||
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
|
||||
%% 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
|
||||
|
||||
clear all
|
||||
close all
|
||||
@ -15,9 +15,9 @@ function reactor2(g)
|
||||
help reactor2
|
||||
|
||||
if nargin == 1
|
||||
gas = g;
|
||||
gas = g;
|
||||
else
|
||||
gas = Solution('gri30.yaml', 'gri30', 'None');
|
||||
gas = Solution('gri30.yaml', 'gri30', 'None');
|
||||
end
|
||||
|
||||
% set the initial conditions
|
||||
@ -32,34 +32,36 @@ function reactor2(g)
|
||||
nSteps = 100;
|
||||
tim(nSteps) = 0;
|
||||
temp(nSteps) = 0;
|
||||
x(nSteps,3) = 0;
|
||||
x(nSteps, 3) = 0;
|
||||
t = 0;
|
||||
dt = 1.0e-5;
|
||||
t0 = cputime;
|
||||
|
||||
for n = 1:100
|
||||
t = t + dt;
|
||||
network.advance(t);
|
||||
tim(n) = network.time;
|
||||
temp(n) = r.T;
|
||||
x(n,1:3) = gas.moleFraction({'OH','H','H2'});
|
||||
t = t + dt;
|
||||
network.advance(t);
|
||||
tim(n) = network.time;
|
||||
temp(n) = r.T;
|
||||
x(n, 1:3) = gas.moleFraction({'OH', 'H', 'H2'});
|
||||
end
|
||||
|
||||
disp(['CPU time = ' num2str(cputime - t0)]);
|
||||
|
||||
clf;
|
||||
subplot(2,2,1);
|
||||
plot(tim,temp);
|
||||
subplot(2, 2, 1);
|
||||
plot(tim, temp);
|
||||
xlabel('Time (s)');
|
||||
ylabel('Temperature (K)');
|
||||
subplot(2,2,2)
|
||||
plot(tim,x(:,1));
|
||||
subplot(2, 2, 2)
|
||||
plot(tim, x(:, 1));
|
||||
xlabel('Time (s)');
|
||||
ylabel('OH Mole Fraction (K)');
|
||||
subplot(2,2,3)
|
||||
plot(tim,x(:,2));
|
||||
subplot(2, 2, 3)
|
||||
plot(tim, x(:, 2));
|
||||
xlabel('Time (s)');
|
||||
ylabel('H Mole Fraction (K)');
|
||||
subplot(2,2,4)
|
||||
plot(tim,x(:,3));
|
||||
subplot(2, 2, 4)
|
||||
plot(tim, x(:, 3));
|
||||
xlabel('Time (s)');
|
||||
ylabel('H2 Mole Fraction (K)');
|
||||
|
||||
|
||||
@ -1,31 +1,31 @@
|
||||
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
|
||||
%
|
||||
%% 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);
|
||||
[m, n] = size(y);
|
||||
dydt = zeros(m, n);
|
||||
|
||||
for j = 1:n
|
||||
this_y = y(:,j);
|
||||
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;
|
||||
u_mass = int_energy / total_mass;
|
||||
v_mass = vol / total_mass;
|
||||
|
||||
% set the state of the gas by specifying (u,v,{Y_k})
|
||||
gas.Y = masses;
|
||||
@ -44,8 +44,9 @@ function dydt = reactor_ode(t, y, gas, vdot, area, heatflux)
|
||||
ydt = total_mass * gas.massProdRate;
|
||||
|
||||
% set up column vector for dydt
|
||||
dydt(:,j) = [udt
|
||||
vdt
|
||||
ydt' ];
|
||||
dydt(:, j) = [udt
|
||||
vdt
|
||||
ydt'];
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@ -17,14 +17,14 @@ help surfreactor
|
||||
%% Set the initial conditions
|
||||
|
||||
t = 870.0;
|
||||
gas = Solution('ptcombust.yaml','gas');
|
||||
gas = Solution('ptcombust.yaml', 'gas');
|
||||
|
||||
gas.TPX = {t, oneatm, 'CH4:0.01, O2:0.21, N2:0.78'};
|
||||
|
||||
% The surface reaction mechanism describes catalytic combustion of
|
||||
% methane on platinum, and is from Deutschman et al., 26th
|
||||
% Symp. (Intl.) on Combustion,1996, pp. 1747-1754
|
||||
surf = Interface('ptcombust.yaml','Pt_surf', gas);
|
||||
surf = Interface('ptcombust.yaml', 'Pt_surf', gas);
|
||||
surf.T = t;
|
||||
|
||||
nsp = gas.nSpecies;
|
||||
@ -35,7 +35,7 @@ r = IdealGasReactor(gas);
|
||||
r.setInitialVolume(1.0e-6)
|
||||
|
||||
% create a reservoir to represent the environment
|
||||
a = Solution('air.yaml','air','None');
|
||||
a = Solution('air.yaml', 'air', 'None');
|
||||
a.TP = {t, oneatm};
|
||||
env = Reservoir(a);
|
||||
|
||||
@ -52,7 +52,7 @@ rsurf = ReactorSurface(surf, r, A);
|
||||
|
||||
% set the wall area and heat transfer coefficient.
|
||||
w.area = A;
|
||||
w.setHeatTransferCoeff(1.0e1); % W/m2/K
|
||||
w.setHeatTransferCoeff(1.0e1); % W/m2/K
|
||||
|
||||
% set expansion rate parameter. dV/dt = KA(P_1 - P_2)
|
||||
w.setExpansionRateCoeff(1.0);
|
||||
@ -62,7 +62,7 @@ network = ReactorNet({r});
|
||||
|
||||
nSteps = 100;
|
||||
p0 = r.P;
|
||||
names = {'CH4','CO','CO2','H2O'};
|
||||
names = {'CH4', 'CO', 'CO2', 'H2O'};
|
||||
x = zeros([nSteps 4]);
|
||||
tim = zeros(nSteps, 1);
|
||||
temp = zeros(nSteps, 1);
|
||||
@ -71,36 +71,38 @@ cov = zeros([nSteps nSurfSp]);
|
||||
t = 0;
|
||||
dt = 0.1;
|
||||
t0 = cputime;
|
||||
|
||||
for n = 1:nSteps
|
||||
t = t + dt;
|
||||
network.advance(t);
|
||||
tim(n) = t;
|
||||
temp(n) = r.T;
|
||||
pres(n) = r.P - p0;
|
||||
cov(n,:) = surf.X';
|
||||
x(n,:) = gas.moleFraction(names);
|
||||
t = t + dt;
|
||||
network.advance(t);
|
||||
tim(n) = t;
|
||||
temp(n) = r.T;
|
||||
pres(n) = r.P - p0;
|
||||
cov(n, :) = surf.X';
|
||||
x(n, :) = gas.moleFraction(names);
|
||||
end
|
||||
|
||||
disp(['CPU time = ' num2str(cputime - t0)]);
|
||||
|
||||
%% Plotting
|
||||
|
||||
clf;
|
||||
subplot(2,2,1);
|
||||
plot(tim,temp);
|
||||
subplot(2, 2, 1);
|
||||
plot(tim, temp);
|
||||
xlabel('Time (s)');
|
||||
ylabel('Temperature (K)');
|
||||
subplot(2,2,2);
|
||||
plot(tim,pres);
|
||||
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);
|
||||
subplot(2, 2, 3);
|
||||
semilogy(tim, cov);
|
||||
xlabel('Time (s)');
|
||||
ylabel('Coverages');
|
||||
legend(speciesNames(surf));
|
||||
subplot(2,2,4);
|
||||
plot(tim,x);
|
||||
subplot(2, 2, 4);
|
||||
plot(tim, x);
|
||||
xlabel('Time (s)');
|
||||
ylabel('Mole Fractions');
|
||||
legend(names);
|
||||
@ -12,7 +12,7 @@ surfreactor;
|
||||
periodic_cstr;
|
||||
plug_flow_reactor;
|
||||
lithium_ion_battery
|
||||
rankine(300.0, 2.0*oneatm, 0.8, 0.7);
|
||||
rankine;
|
||||
prandtl1();
|
||||
prandtl2();
|
||||
flame1;
|
||||
@ -24,5 +24,4 @@ ignite_uv;
|
||||
|
||||
clear all
|
||||
close all
|
||||
cleanup
|
||||
UnloadCantera
|
||||
|
||||
Loading…
Reference in New Issue
Block a user