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

Added examples to the new Matlab toolbox

Browse Source
This commit is contained in:
Su Sun 2022-01-11 14:36:53 -05:00 committed by Ray Speth
parent 0189abb705
commit 1e59450ca4
15 changed files with 1277 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

168
samples/matlab_experimental/PFR.m Normal file
View File

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

53
samples/matlab_experimental/PFR_Solver.m Normal file
View File

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

63
samples/matlab_experimental/equil.m Normal file
View File

@ -0,0 +1,63 @@
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.
help equil
if nargin == 1
gas = g;
else
gas = Solution('gri30.yaml');
end
nsp = gas.nSpecies;
% find methane, nitrogen, and oxygen indices
ich4 = gas.speciesIndex('CH4');
io2 = gas.speciesIndex('O2');
in2 = gas.speciesIndex('N2');
nPhis = 50;
phi = linspace(0.2, 2.70, nPhis);
tad(nPhis) = 0;
xeq(nsp, nPhis) = 0;
for i = 1:nPhis
x = zeros(nsp, 1);
x(ich4, 1) = phi(i);
x(io2, 1) = 2.0;
x(in2, 1) = 7.52;
gas.T = 300;
gas.P = 101325;
gas.X = x;
gas.equilibrate('HP');
tad(i) = gas.T;
xeq(:, i) = gas.X;
end
% make plots
clf;
subplot(1, 2, 1);
plot(phi, tad);
xlabel('Equivalence Ratio');
ylabel('Temperature (K)');
title('Adiabatic Flame Temperature');
subplot(1,2,2);
semilogy(phi, xeq);
axis([phi(1), phi(50), 1.0e-14, 1]);
%legend(speciesName(gas,1:nsp),1);
j = 10;
for k = 1:nsp
text(phi(j), 1.5*xeq(k, j), gas.speciesName(k))
j = j + 2;
if j > 46
j = 10;
end
end
xlabel('Equivalence Ratio');
ylabel('Mole Fraction');
title('Equilibrium Composition');
end

73
samples/matlab_experimental/flame.m Normal file
View File

@ -0,0 +1,73 @@
function f = flame(gas, left, flow, right)
% FLAME - create a flame object.
%
% Check input parameters
if nargin ~= 4
error('wrong number of input arguments.');
end
if ~gas.thermo.isIdealGas
error('gas object must represent an ideal gas mixture.');
end
if ~isInlet(left)
error('burner object of wrong type.');
end
if ~isFlow(flow)
error('flow object of wrong type.');
end
flametype = 0;
if isSurface(right)
flametype = 1;
elseif isInlet(right)
flametype = 3;
end
% create the container object
f = Stack([left flow right]);
% set default initial profiles.
rho0 = gas.D;
% find the adiabatic flame temperature and corresponding
% equilibrium composition
gas.equilibrate('HP');
teq = gas.T;
yeq = gas.Y;
z1 = 0.2;
mdot0 = left.massFlux;
mdot1 = right.massFlux;
t0 = left.T;
if flametype == 0
t1 = teq;
mdot1 = -mdot0;
else
t1 = right.T;
end
f.setProfile(2, {'u', 'V'}, [0.0 1.0
mdot0/rho0 -mdot1/rho0
0.0 0.0]);
f.setProfile(2, 'T', [0.0 z1 1.0; t0 2000.0 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, z1, 1
left.massFraction(n), yint, y1]);
end
% set minimal grid refinement criteria
f.setRefineCriteria(2, 10.0, 0.8, 0.8);
end

103
samples/matlab_experimental/ignite.m Normal file
View File

@ -0,0 +1,103 @@
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.
%
help ignite
if nargin == 1
gas = g;
else
gas = Solution('gri30.yaml');
end
% set the initial conditions
gas.TPX = {1001.0, oneatm, 'X','H2:2,O2:1,N2:4'};
gas.basis = 'mass';
y0 = [gas.U
1.0/gas.D
gas.Y'];
time_interval = [0 0.001];
options = odeset('RelTol',1.e-5,'AbsTol',1.e-12,'Stats','on');
t0 = cputime;
out = ode15s(@reactor_ode,time_interval,y0,options,gas,@vdot,@area,@heatflux);
disp(['CPU time = ' num2str(cputime - t0)]);
plotdata = output(out,gas);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the functions below may be defined arbitrarily to set the reactor
% boundary conditions - the rate of change of volume, the heat
% flux, and the area.
% Rate of change of volume. Any arbirtrary function may be implemented.
% Input arguments:
% t time
% vol volume
% gas ideal gas object
function v = vdot(t, vol, gas)
%v = 0.0; %uncomment for constant volume
v = 1.e11 * (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
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
% not necessarily those desired for output, this function is called
% after the integration is complete to generate the desired
% outputs.
function pv = output(s, gas)
times = s.x;
soln = s.y;
[~,n] = size(times);
pv = zeros(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};
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.
figure(1);
plot(pv(1, :), pv(2, :));
xlabel('time');
ylabel('Temperature');
title(['Final T = ' num2str(pv(2, end)) ' K']);
figure(2);
ioh = gas.speciesIndex('OH');
plot(pv(1,:),pv(4+ioh,:));
xlabel('time');
ylabel('Mass Fraction');
title('OH Mass Fraction');
end
end

59
samples/matlab_experimental/isentropic.m Normal file
View File

@ -0,0 +1,59 @@
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.
%
help isentropic
if nargin == 1
gas = g;
else
gas = Solution('gri30.yaml');
end
% set the stagnation state
gas.TPX = {1200.0, 10.0 * oneatm, 'H2:1,N2:0.1'};
gas.basis = 'mass';
s0 = gas.S;
h0 = gas.H;
p0 = gas.P;
mdot = 1; % arbitrary
mach = [];
a = [];
i = 1;
amin = 1.e14;
% compute values for a range of pressure ratios
for r = 0.005:0.0025:0.995
p = p0*r;
% set the state using (p,s0)
gas.SP = {s0, p};
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
if a(i) < amin
amin = a(i);
end
mach(i) = v/gas.soundspeed;
i = i + 1;
end
a = a/amin;
% plot results
clf;
plot(mach,a);
ylabel('Area Ratio');
xlabel('Mach Number');
title('Isentropic Flow: Area Ratio vs. Mach Number');
end

142
samples/matlab_experimental/lithium_ion_battery.m Normal file
View File

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

105
samples/matlab_experimental/periodic_cstr.m Normal file
View File

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

76
samples/matlab_experimental/prandtl1.m Normal file
View File

@ -0,0 +1,76 @@
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.
help prandtl1
if nargin == 1
gas = g;
else
gas = GRI30('Mix');
end
pr = zeros(31,31);
xh2 = zeros(31,31);
visc = zeros(31,31);
lambda = zeros(31,31);
t = [];
xo2 = [];
io2 = gas.speciesIndex('O2');
ih2 = gas.speciesIndex('H2');
minT = gas.minTemp;
maxT = gas.maxTemp;
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
end
disp(['CPU time = ' num2str(cputime - t0)]);
% plot results
clf;
%figure(1);
subplot(2,2,1);
surf(xo2,t,pr);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Prandtl Number');
subplot(2,2,2);
surf(xo2,t,xh2);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('H_2 Mole Fraction');
subplot(2,2,3);
surf(xo2,t,visc);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Viscosity');
subplot(2,2,4);
surf(xo2,t,lambda);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Thermal Conductivity');
end

74
samples/matlab_experimental/prandtl2.m Normal file
View File

@ -0,0 +1,74 @@
function prandtl2(g)
% PRANDTL2 Prandtl number for an equilibrium H/O gas mixture.
%
% This example does the same thing as prandtl1, but using
% the multicomponent expression for the thermal conductivity.
%
help prandtl2
if nargin == 1
gas = g;
else
gas = GRI30('Multi');
end
pr = zeros(31,31);
xh2 = zeros(31,31);
visc = zeros(31,31);
lambda = zeros(31,31);
t = [];
xo2 = [];
io2 = gas.speciesIndex('O2');
ih2 = gas.speciesIndex('H2');
minT = gas.minTemp;
maxT = gas.maxTemp;
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
end
disp(['CPU time = ' num2str(cputime - t0)]);
% plot results
clf;
subplot(2,2,1);
surf(xo2,t,pr);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Prandtl Number');
subplot(2,2,2);
surf(xo2,t,xh2);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('H_2 Mole Fraction');
subplot(2,2,3);
surf(xo2,t,visc);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Viscosity');
subplot(2,2,4);
surf(xo2,t,lambda);
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Thermal Conductivity');
end

63
samples/matlab_experimental/rankine.m Normal file
View File

@ -0,0 +1,63 @@
function [work, efficiency] = rankine(t1, p2, eta_pump, eta_turbine)
% This example computes the efficiency of a simple vapor power cycle.
help rankine
% create an object representing water
w = Water;
% start with saturated liquid water at t1
w.setState_Tsat(t1, 1.0);
p1 = w.P;
% pump it to p2
basis = 'mass';
pump_work = pump(w, p2, eta_pump);
h2 = w.H;
p2 = w.P;
% heat to saturated vapor
w.setState_Psat(p2, 1.0);
h3 = w.H;
heat_added = h3 - h2;
% expand adiabatically back to the initial pressure
work = expand(w, p1, eta_turbine);
% compute the efficiency
efficiency = (work - pump_work)/heat_added;
end
function w = pump(fluid, pfinal, eta)
% PUMP - Adiabatically pump a fluid to pressure pfinal, using a pump
% with isentropic efficiency eta.
fluid.basis = 'mass';
h0 = fluid.H;
s0 = fluid.S;
fluid.SP = {s0, pfinal};
h1s = fluid.H;
isentropic_work = h1s - h0;
actual_work = isentropic_work / eta;
h1 = h0 + actual_work;
fluid.HP = {h1, pfinal};
w = actual_work;
end
function w = expand(fluid, pfinal, eta)
% EXPAND - Adiabatically expand a fluid to pressure pfinal, using a
% turbine with isentropic efficiency eta.
fluid.basis = 'mass';
h0 = fluid.H;
s0 = fluid.S;
fluid.SP = {s0, pfinal};
h1s = fluid.H;
isentropic_work = h0 - h1s;
actual_work = isentropic_work * eta;
h1 = h0 - actual_work;
fluid.HP = {h1, pfinal};
w = actual_work;
end

85
samples/matlab_experimental/reactor1.m Normal file
View File

@ -0,0 +1,85 @@
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.
help reactor1
if nargin == 1
gas = g;
else
gas = GRI30('None');
end
P = 101325.0;
% set the initial conditions
gas.T = 1001.0;
gas.P = P;
nsp = gas.nSpecies;
xx = zeros(nsp, 1);
xx(1) = 0.285;
xx(4) = 0.142;
xx(48) = 0.573;
gas.X = xx;
% create a reactor, and insert the gas
r = IdealGasReactor(gas);
% create a reservoir to represent the environment
a = Solution('air.yaml','air','None');
a.P = P;
env = Reservoir(a);
% Define a wall between the reactor and the environment and
% make it flexible, so that the pressure in the reactor is held
% at the environment pressure.
w = Wall;
w.install(r,env);
% set expansion parameter. dV/dt = KA(P_1 - P_2)
w.setExpansionRateCoeff(1.0e6);
% set wall area
w.area = 1.0;
% create a reactor network and insert the reactor:
network = ReactorNet({r});
nSteps = 100;
tim(nSteps) = 0;
temp(nSteps) = 0;
x(nSteps,3) = 0;
t = 0.0;
dt = 1.0e-5;
t0 = cputime;
for n = 1:nSteps
t = t + dt;
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);
xlabel('Time (s)');
ylabel('Temperature (K)');
subplot(2,2,2)
plot(tim,x(:,1));
xlabel('Time (s)');
ylabel('OH Mole Fraction (K)');
subplot(2,2,3)
plot(tim,x(:,2));
xlabel('Time (s)');
ylabel('H Mole Fraction (K)');
subplot(2,2,4)
plot(tim,x(:,3));
xlabel('Time (s)');
ylabel('H2 Mole Fraction (K)');
clear all
cleanup
end

60
samples/matlab_experimental/reactor2.m Normal file
View File

@ -0,0 +1,60 @@
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.
help reactor2
if nargin == 1
gas = g;
else
gas = GRI30('None');
end
% set the initial conditions
gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'};
% create a reactor, and insert the gas
r = IdealGasReactor(gas);
% create a reactor network and insert the reactor
network = ReactorNet({r});
nSteps = 100;
tim(nSteps) = 0;
temp(nSteps) = 0;
x(nSteps,3) = 0;
t = 0;
dt = 1.0e-5;
t0 = cputime;
for n = 1:100
t = t + dt;
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);
xlabel('Time (s)');
ylabel('Temperature (K)');
subplot(2,2,2)
plot(tim,x(:,1));
xlabel('Time (s)');
ylabel('OH Mole Fraction (K)');
subplot(2,2,3)
plot(tim,x(:,2));
xlabel('Time (s)');
ylabel('H Mole Fraction (K)');
subplot(2,2,4)
plot(tim,x(:,3));
xlabel('Time (s)');
ylabel('H2 Mole Fraction (K)');
clear all
cleanup
end

51
samples/matlab_experimental/reactor_ode.m Normal file
View File

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

102
samples/matlab_experimental/surfreactor.m Normal file
View File

@ -0,0 +1,102 @@
% SURFREACTOR Zero-dimensional reactor with surface chemistry
%
% This example illustrates how to use class 'Reactor' for
% zero-dimensional simulations including both homogeneous and
% heterogeneous chemistry.
help surfreactor
clear all
close all
cleanup
clc
t = 870.0;
gas = Solution('ptcombust.yaml','gas');
% set the initial conditions
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 = importInterface('ptcombust.yaml','Pt_surf', gas);
surf.th.T = t;
nsp = gas.nSpecies;
nSurfSp = surf.th.nSpecies;
% create a reactor, and insert the gas
r = IdealGasReactor(gas);
r.setInitialVolume(1.0e-6)
% create a reservoir to represent the environment
a = Solution('air.yaml','air','None');
a.TP = {t, oneatm};
env = Reservoir(a);
% Define a wall between the reactor and the environment and
% make it flexible, so that the pressure in the reactor is held
% at the environment pressure.
w = Wall;
w.install(r, env);
A = 1e-4; % Wall area
% Add a reacting surface, with an area matching that of the wall
rsurf = ReactorSurface(surf, r, A);
% set the wall area and heat transfer coefficient.
w.area = A;
w.setHeatTransferCoeff(1.0e1); % W/m2/K
% set expansion rate parameter. dV/dt = KA(P_1 - P_2)
w.setExpansionRateCoeff(1.0);
network = ReactorNet({r});
% setTolerances(network, 1.0e-8, 1.0e-12);
nSteps = 100;
p0 = r.P;
names = {'CH4','CO','CO2','H2O'};
x = zeros([nSteps 4]);
tim = zeros(nSteps);
temp = zeros(nSteps);
pres = zeros(nSteps);
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.coverages';
x(n,:) = gas.moleFraction(names);
end
disp(['CPU time = ' num2str(cputime - t0)]);
clf;
subplot(2,2,1);
plot(tim,temp);
xlabel('Time (s)');
ylabel('Temperature (K)');
subplot(2,2,2);
plot(tim,pres);
axis([0 5 -0.1 0.1]);
xlabel('Time (s)');
ylabel('Delta Pressure (Pa)');
subplot(2,2,3);
semilogy(tim,cov);
xlabel('Time (s)');
ylabel('Coverages');
legend(speciesNames(surf));
subplot(2,2,4);
plot(tim,x);
xlabel('Time (s)');
ylabel('Mole Fractions');
legend(names);
clear all
cleanup