From c9a4004ab5845abababa98c32f2ad462d780fcc1 Mon Sep 17 00:00:00 2001 From: Su Sun <56568705+ssun30@users.noreply.github.com> Date: Fri, 21 Jan 2022 16:55:08 -0500 Subject: [PATCH] Added and fixed Examples for the new toolbox Also fix Examples/flame1.m --- samples/matlab_experimental/catcomb.m | 239 ++++++++++++++++++++++ samples/matlab_experimental/conhp.m | 27 +++ samples/matlab_experimental/conuv.m | 31 +++ samples/matlab_experimental/diff_flame.m | 141 +++++++++++++ samples/matlab_experimental/flame.m | 10 +- samples/matlab_experimental/flame1.m | 122 +++++++++++ samples/matlab_experimental/flame2.m | 126 ++++++++++++ samples/matlab_experimental/ignite.m | 23 ++- samples/matlab_experimental/ignite_hp.m | 38 ++++ samples/matlab_experimental/ignite_uv.m | 38 ++++ samples/matlab_experimental/prandtl1.m | 5 + samples/matlab_experimental/prandtl2.m | 5 + samples/matlab_experimental/surfreactor.m | 10 +- 13 files changed, 794 insertions(+), 21 deletions(-) create mode 100644 samples/matlab_experimental/catcomb.m create mode 100644 samples/matlab_experimental/conhp.m create mode 100644 samples/matlab_experimental/conuv.m create mode 100644 samples/matlab_experimental/diff_flame.m create mode 100644 samples/matlab_experimental/flame1.m create mode 100644 samples/matlab_experimental/flame2.m create mode 100644 samples/matlab_experimental/ignite_hp.m create mode 100644 samples/matlab_experimental/ignite_uv.m diff --git a/samples/matlab_experimental/catcomb.m b/samples/matlab_experimental/catcomb.m new file mode 100644 index 000000000..2b786ed58 --- /dev/null +++ b/samples/matlab_experimental/catcomb.m @@ -0,0 +1,239 @@ +% Catalytic combustion of a stagnation flow on a platinum surface +% +% This script solves a catalytic combustion problem. A stagnation flow +% is set up, with a gas inlet 10 cm from a platinum surface at 900 +% K. The lean, premixed methane/air mixture enters at ~ 6 cm/s (0.06 +% kg/m2/s), and burns catalytically on the platinum surface. Gas-phase +% chemistry is included too, and has some effect very near the +% surface. +% +% The catalytic combustion mechanism is from Deutschman et al., 26th +% Symp. (Intl.) on Combustion,1996 pp. 1747-1754 + +%% Initialization + +help catcomb; + +clear all +close all +cleanup +clc + +t0 = cputime; % record the starting time + +%% Set parameter values + +p = oneatm; % pressure +tinlet = 300.0; % inlet temperature +tsurf = 900.0; % surface temperature +mdot = 0.06; % kg/m^2/s +transport = 'mixture-averaged'; % transport model + +% Solve first for a hydrogen/air case for use as the initial estimate for +% the methane/air case. + +% composition of the inlet premixed gas for the hydrogen/air case +comp1 = 'H2:0.05, O2:0.21, N2:0.78, AR:0.01'; + +% composition of the inlet premixed gas for the methane/air case +comp2 = 'CH4:0.095, O2:0.21, N2:0.78, AR:0.01'; + +% the initial grid, in meters. The inlet/surface separation is 10 cm. +initial_grid = [0.0, 0.02, 0.04, 0.06, 0.08, 0.1]; % m + +% numerical parameters +tol_ss = {1.0e-8 1.0e-14}; % {rtol atol} for steady-state problem +tol_ts = {1.0e-4 1.0e-9}; % {rtol atol} for time stepping + +loglevel = 1; % amount of diagnostic output + % (0 to 5) + +refine_grid = 1; % 1 to enable refinement, 0 to + % disable + +%% Create the gas object +% +% This object will be used to evaluate all thermodynamic, kinetic, +% and transport properties +% +% The gas phase will be taken from the definition of phase 'gas' in +% input file 'ptcombust.yaml,' which is a stripped-down version of +% GRI-Mech 3.0. + +gas = Solution('ptcombust.yaml', 'gas', transport); +gas.TPX = {tinlet, p, comp1}; + +%% create the interface object +% +% This object will be used to evaluate all surface chemical production +% rates. It will be created from the interface definition 'Pt_surf' +% in input file 'ptcombust.yaml,' which implements the reaction +% mechanism of Deutschmann et al., 1995 for catalytic combustion on +% platinum. + +surf_phase = importInterface('ptcombust.yaml', 'Pt_surf', gas); +surf_phase.T = tsurf; + +% integrate the coverage equations in time for 1 s, holding the gas +% composition fixed to generate a good starting estimate for the +% coverages. + +surf_phase.advanceCoverages(1.0); + +% The two objects we just created are independent of the problem +% type -- they are useful in zero-D simulations, 1-D simulations, +% etc. Now we turn to creating the objects that are specifically +% for 1-D simulations. These will be 'stacked' together to create +% the complete simulation. + +%% create the flow object +% +% The flow object is responsible for evaluating the 1D governing +% equations for the flow. We will initialize it with the gas +% object, and assign it the name 'flow'. + +flow = AxisymmetricFlow(gas, 'flow'); + +% set some parameters for the flow +flow.setPressure(p); +flow.setupGrid(initial_grid); +flow.setSteadyTolerances('default', tol_ss{:}); +flow.setTransientTolerances('default', tol_ts{:}); + +%% create the inlet +% +% The temperature, mass flux, and composition (relative molar) may be +% specified. This object provides the inlet boundary conditions for +% the flow equations. + +inlt = Inlet('inlet'); + +% set the inlet parameters. Start with comp1 (hydrogen/air) +inlt.T = tinlet; +inlt.setMdot(mdot); +inlt.setMoleFractions(comp1); + +%% create the surface +% +% This object provides the surface boundary conditions for the flow +% equations. By supplying object surface_phase as an argument, the +% coverage equations for its surface species will be added to the +% equation set, and used to compute the surface production rates of +% the gas-phase species. + +surf = Surface('surface', surf_phase); +surf.T = tsurf; + +%% create the stack +% +% Once the component parts have been created, they can be assembled +% to create the 1D simulation. + +sim1D = Stack([inlt, flow, surf]); + +% set the initial profiles. +sim1D.setProfile(2, {'u', 'V', 'T'}, [0.0, 1.0 % z/zmax + 0.06, 0.0 % u + 0.0, 0.0 % V + tinlet, tsurf]); % T +names = gas.speciesNames; + +for k = 1:gas.nSpecies + y = inlt.massFraction(k); + sim1D.setProfile(2, names{k}, [0, 1; y, y]); +end + +sim1D + +%setTimeStep(fl, 1.0e-5, [1, 3, 6, 12]); +%setMaxJacAge(fl, 4, 5); + +%% Solution + +% start with the energy equation on +flow.enableEnergy; + +% disable the surface coverage equations, and turn off all gas and +% surface chemistry + +surf.setCoverageEqs('off'); +surf_phase.setMultiplier(0.0); +gas.setMultiplier(0.0); + +% solve the problem, refining the grid if needed +sim1D.solve(1, refine_grid); + +% now turn on the surface coverage equations, and turn the +% chemistry on slowly + +surf.setCoverageEqs('on'); + +for iter=1:6 + mult = 10.0^(iter - 6); + surf_phase.setMultiplier(mult); + gas.setMultiplier(mult); + sim1D.solve(1, refine_grid); +end + +% At this point, we should have the solution for the hydrogen/air +% problem. Now switch the inlet to the methane/air composition. +inlt.setMoleFractions(comp2); + +% set more stringent grid refinement criteria +sim1D.setRefineCriteria(2, 100.0, 0.15, 0.2); + +% solve the problem for the final time +sim1D.solve(loglevel, refine_grid); + +% show the solution +sim1D + +% save the solution +sim1D.saveSoln('catcomb.xml', 'energy', ['solution with energy equation']); + +%% Show statistics + +sim1D.writeStats; +elapsed = cputime - t0; +e = sprintf('Elapsed CPU time: %10.4g',elapsed); +disp(e); + +%% Make plots + +clf; + +subplot(3, 3, 1); +sim1D.plotSolution('flow', 'T'); +title('Temperature [K]'); + +subplot(3, 3, 2); +sim1D.plotSolution('flow', 'u'); +title('Axial Velocity [m/s]'); + +subplot(3, 3, 3); +sim1D.plotSolution('flow', 'V'); +title('Radial Velocity / Radius [1/s]'); + +subplot(3, 3, 4); +sim1D.plotSolution('flow', 'CH4'); +title('CH4 Mass Fraction'); + +subplot(3, 3, 5); +sim1D.plotSolution('flow', 'O2'); +title('O2 Mass Fraction'); + +subplot(3, 3, 6); +sim1D.plotSolution('flow', 'CO'); +title('CO Mass Fraction'); + +subplot(3, 3, 7); +sim1D.plotSolution('flow', 'CO2'); +title('CO2 Mass Fraction'); + +subplot(3, 3, 8); +sim1D.plotSolution('flow', 'H2O'); +title('H2O Mass Fraction'); + +subplot(3, 3, 9); +sim1D.plotSolution('flow', 'H2'); +title('H2 Mass Fraction'); diff --git a/samples/matlab_experimental/conhp.m b/samples/matlab_experimental/conhp.m new file mode 100644 index 000000000..5e1de06a9 --- /dev/null +++ b/samples/matlab_experimental/conhp.m @@ -0,0 +1,27 @@ +function dydt = conhp(t, y, gas, mw) %#ok + % CONHP - ODE system for a constant-pressure, adiabatic reactor. + % + % Function CONHP evaluates the system of ordinary differential + % equations for an adiabatic, constant-pressure, + % zero-dimensional reactor. It assumes that the 'gas' object + % represents a reacting ideal gas mixture. + + % Set the state of the gas, based on the current solution vector. + gas.Y = y(2: end); + gas.TP = {y(1), gas.P}; + nsp = gas.nSpecies; + + % energy equation + wdot = gas.netProdRates; + gas.basis = 'mass'; + tdot = - gas.T * gasconstant * gas.enthalpies_RT .* wdot / (gas.D * gas.cp); + + % set up column vector for dydt + dydt = [tdot' + zeros(nsp, 1)]; + + % species equations + rrho = 1.0/gas.D; + for i = 1:nsp + dydt(i+1) = rrho * mw(i) * wdot(i); +end diff --git a/samples/matlab_experimental/conuv.m b/samples/matlab_experimental/conuv.m new file mode 100644 index 000000000..1bede4c6e --- /dev/null +++ b/samples/matlab_experimental/conuv.m @@ -0,0 +1,31 @@ +function dydt = conuv(t, y, gas, mw) %#ok + % CONUV ODE system for a constant-volume, adiabatic reactor. + % + % Function CONUV evaluates the system of ordinary differential + % equations for an adiabatic, constant-volume, + % zero-dimensional reactor. It assumes that the 'gas' object + % represents a reacting ideal gas mixture. + + + % Set the state of the gas, based on the current solution vector. + gas.Y = y(2:end); + gas.TD = {y(1), gas.D}; + nsp = gas.nSpecies; + + % energy equation + wdot = gas.netProdRates; + gas.basis = 'mass'; + tdot = - gas.T * gasconstant * (gas.enthalpies_RT - ones(1, nsp)) ... + .* wdot / (gas.D * gas.cv); + + % set up column vector for dydt + dydt = [tdot' + zeros(nsp, 1)]; + + % species equations + rrho = 1.0/gas.D; + for i = 1:nsp + dydt(i+1) = rrho * mw(i) * wdot(i); + end + +end diff --git a/samples/matlab_experimental/diff_flame.m b/samples/matlab_experimental/diff_flame.m new file mode 100644 index 000000000..c2620ec54 --- /dev/null +++ b/samples/matlab_experimental/diff_flame.m @@ -0,0 +1,141 @@ +% DIFF_FLAME - An opposed-flow diffusion flame. +% This example uses the CounterFlowDiffusionFlame function to solve an +% opposed-flow diffusion flame for Ethane in Air. This example is the same +% as the diffusion_flame.py example without radiation. + +%% Initialization + +help diff_flame + +clear all +close all +cleanup +clc + +runtime = cputime; % Record the starting time + +%% Parameter values of inlet streams + +p = oneatm; % Pressure +tin = 300.0; % Inlet temperature +mdot_o = 0.72; % Air mass flux, kg/m^2/s +mdot_f = 0.24; % Fuel mass flux, kg/m^2/s +transport = 'Mix'; % Transport model +% NOTE: Transport model needed if mechanism file does not have transport +% properties. + +%% Set-up initial grid, loglevel, tolerances. Enable/Disable grid refinement + +initial_grid = 0.02*[0.0, 0.2, 0.4, 0.6, 0.8, 1.0]; % Units: m +tol_ss = {1.0e-5, 1.0e-9}; % {rtol atol} for steady-state problem +tol_ts = {1.0e-3, 1.0e-9}; % {rtol atol} for time stepping +loglevel = 1; % Amount of diagnostic output (0 to 5) +refine_grid = 1; % 1 to enable refinement, 0 to disable + +%% Create the gas objects for the fuel and oxidizer streams +% +% These objects will be used to evaluate all thermodynamic, kinetic, and +% transport properties. + +fuel = GRI30(transport); +ox = GRI30(transport); +oxcomp = 'O2:0.21, N2:0.78'; % Air composition +fuelcomp = 'C2H6:1'; % Fuel composition +% Set each gas mixture state with the corresponding composition. +fuel.TPX = {tin, p, fuelcomp}; +ox.TPX = {tin, p, oxcomp}; + +%% Set-up the flow object +% +% For this problem, the AxisymmetricFlow model is needed. Set the state of +% the flow as the fuel gas object. This is arbitrary and is only used to +% initialize the flow object. Set the grid to the initial grid defined +% prior, same for the tolerances. + +f = AxisymmetricFlow(fuel, 'flow'); +f.setPressure(p); +f.setupGrid(initial_grid); +f.setSteadyTolerances('default', tol_ss{:}); +f.setTransientTolerances('default', tol_ts{:}); + +%% Create the fuel and oxidizer inlet steams +% +% Specify the temperature, mass flux, and composition correspondingly. + +% Set the oxidizer inlet. +inlet_o = Inlet('air_inlet'); +inlet_o.T = tin; +inlet_o.setMdot(mdot_o); +inlet_o.setMoleFractions(oxcomp); + +% Set the fuel inlet. +inlet_f = Inlet('fuel_inlet'); +inlet_f.T = tin; +inlet_f.setMdot(mdot_f); +inlet_f.setMoleFractions(fuelcomp); + +%% Create the flame object. +% +% Once the inlets have been created, they can be assembled +% to create the flame object. Function CounterFlorDiffusionFlame +% (in Cantera/1D) sets up the initial guess for the solution using a +% Burke-Schumann flame. The input parameters are: fuel inlet object, flow +% object, oxidizer inlet object, fuel gas object, oxidizer gas object, and +% the name of the oxidizer species as in character format. + +fl = CounterFlowDiffusionFlame(inlet_f, f, inlet_o, fuel, ox, 'O2'); + +%% Solve with fixed temperature profile +% +% Grid refinement is turned off for this process in this example. +% To turn grid refinement on, change 0 to 1 for last input is solve function. + +fl.solve(loglevel, 0); + +%% Enable the energy equation. +% +% The energy equation will now be solved to compute the temperature profile. +% We also tighten the grid refinement criteria to get an accurate final +% solution. The explanation of the setRefineCriteria function is located +% on cantera.org in the Matlab User's Guide and can be accessed by +% help setRefineCriteria + +f.enableEnergy; +fl.setRefineCriteria(2, 200.0, 0.1, 0.2); +fl.solve(loglevel, refine_grid); +fl.saveSoln('c2h6.xml', 'energy', ['solution with energy equation']); + +%% Show statistics of solution and elapsed time. + +fl.writeStats; +elapsed = cputime - runtime; +e = sprintf('Elapsed CPU time: %10.4g',elapsed); +disp(e); + +% Make a single plot showing temperature and mass fraction of select +% species along axial distance from fuel inlet to air inlet. +% + +z = fl.grid('flow'); % Get grid points of flow +spec = fuel.speciesNames; % Get species names in gas +T = fl.solution('flow', 'T'); % Get temperature solution + +for i = 1:length(spec) + % Get mass fraction of all species from solution + y(i, :) = fl.solution('flow', spec{i}); +end + +j = fuel.speciesIndex('O2'); % Get index of O2 in gas object +k = fuel.speciesIndex('H2O'); % Get index of H2O in gas object +l = fuel.speciesIndex('C2H6'); % Get index of C2H6 in gas object +m = fuel.speciesIndex('CO2'); % Get index of CO2 in gas object + +clf; +yyaxis left +plot(z, T) +xlabel('z (m)'); +ylabel('Temperature (K)'); +yyaxis right +plot(z, y(j, :), 'r', z, y(k, :), 'g', z, y(l, :), 'm', z, y(m, :), 'b'); +ylabel('Mass Fraction'); +legend('T', 'O2', 'H2O', 'C2H6', 'CO2'); diff --git a/samples/matlab_experimental/flame.m b/samples/matlab_experimental/flame.m index 8372dde3e..e80c50caf 100644 --- a/samples/matlab_experimental/flame.m +++ b/samples/matlab_experimental/flame.m @@ -6,20 +6,20 @@ function f = flame(gas, left, flow, right) error('wrong number of input arguments.'); end - if ~gas.thermo.isIdealGas + if ~gas.isIdealGas error('gas object must represent an ideal gas mixture.'); end - if ~isInlet(left) + if ~left.isInlet error('burner object of wrong type.'); end - if ~isFlow(flow) + if ~flow.isFlow error('flow object of wrong type.'); end flametype = 0; - if isSurface(right) + if right.isSurface flametype = 1; - elseif isInlet(right) + elseif right.isInlet flametype = 3; end diff --git a/samples/matlab_experimental/flame1.m b/samples/matlab_experimental/flame1.m new file mode 100644 index 000000000..acab9508a --- /dev/null +++ b/samples/matlab_experimental/flame1.m @@ -0,0 +1,122 @@ +% FLAME1 - A burner-stabilized flat flame +% +% This script simulates a burner-stablized lean hydrogen-oxygen flame +% at low pressure. + +%% Initialization + +help flame1 + +clear all +close all +cleanup +clc + +t0 = cputime; % record the starting time + +%% Set parameter values + +p = 0.05*oneatm; % pressure +tburner = 373.0; % burner temperature +mdot = 0.06; % kg/m^2/s + +rxnmech = 'h2o2.yaml'; % reaction mechanism file +comp = 'H2:1.8, O2:1, AR:7'; % premixed gas composition + +initial_grid = [0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.15, 0.2, 0.4,... + 0.49, 0.5]; % m + +tol_ss = {1.0e-5, 1.0e-13}; % {rtol atol} for steady-state + % problem +tol_ts = {1.0e-4, 1.0e-9}; % {rtol atol} for time stepping + +loglevel = 1; % amount of diagnostic output (0 + % to 5) + +refine_grid = 1; % 1 to enable refinement, 0 to + % disable +max_jacobian_age = [5, 10]; + +%% Create the gas object +% +% This object will be used to evaluate all thermodynamic, kinetic, +% and transport properties + +gas = Solution(rxnmech, 'ohmech', 'Mix'); + +% set its state to that of the unburned gas at the burner +gas.TPX = {tburner, p, comp}; + +%% Create the flow object + +f = AxisymmetricFlow(gas, 'flow'); +f.setPressure(p); +f.setupGrid(initial_grid); +f.setSteadyTolerances('default', tol_ss{:}); +f.setTransientTolerances('default', tol_ts{:}); + +%% Create the burner +% +% The burner is an Inlet object. The temperature, mass flux, +% and composition (relative molar) may be specified. +burner = Inlet('burner'); +burner.T = tburner; +burner.setMdot(mdot); +burner.setMoleFractions(comp); + +%% Create the outlet +% +% The type of flame is determined by the object that terminates +% the domain. An Outlet object imposes zero gradient boundary +% conditions for the temperature and mass fractions, and zero +% radial velocity and radial pressure gradient. + +s = Outlet('out'); + +%% Create the flame object +% +% Once the component parts have been created, they can be assembled +% to create the flame object. +% +fl = flame(gas, burner, f, s); +fl.setMaxJacAge(max_jacobian_age(1), max_jacobian_age(2)); + +% if the starting solution is to be read from a previously-saved +% solution, uncomment this line and edit the file name and solution id. +%restore(fl,'h2flame2.xml', 'energy') + +fl.solve(loglevel, refine_grid); + +%% Enable the energy equation +% +% The energy equation will now be solved to compute the +% temperature profile. We also tighten the grid refinement +% criteria to get an accurate final solution. + +f.enableEnergy; +fl.setRefineCriteria(2, 200.0, 0.05, 0.1); +fl.solve(1, 1); +fl.saveSoln('h2fl.xml', 'energy', ['solution with energy equation']); + +%% Show statistics + +fl.writeStats; +elapsed = cputime - t0; +e = sprintf('Elapsed CPU time: %10.4g',elapsed); +disp(e); + +%% Make plots + +clf; +subplot(2, 2, 1); +plotSolution(fl, 'flow', 'T'); +title('Temperature [K]'); +subplot(2, 2, 2); +plotSolution(fl, 'flow', 'u'); +title('Axial Velocity [m/s]'); +subplot(2, 2, 3); +plotSolution(fl, 'flow', 'H2O'); +title('H2O Mass Fraction'); +subplot(2, 2, 4); +plotSolution(fl, 'flow', 'O2'); +title('O2 Mass Fraction'); diff --git a/samples/matlab_experimental/flame2.m b/samples/matlab_experimental/flame2.m new file mode 100644 index 000000000..f9b68cd2f --- /dev/null +++ b/samples/matlab_experimental/flame2.m @@ -0,0 +1,126 @@ +% FLAME2 - An axisymmetric stagnation-point non-premixed flame +% +% This script simulates a stagnation-point ethane-air flame. + +%% Initialization + +help flame2 + +clear all +close all +cleanup +clc + +t0 = cputime; % record the starting time + +%% Set parameter values + +p = oneatm; % pressure +tin = 300.0; % inlet temperature +mdot_o = 0.72; % air, kg/m^2/s +mdot_f = 0.24; % fuel, kg/m^2/s + +rxnmech = 'gri30.yaml'; % reaction mechanism file +comp1 = 'O2:0.21, N2:0.78, AR:0.01'; % air composition +comp2 = 'C2H6:1'; % fuel composition + +initial_grid = 0.02*[0.0, 0.2, 0.4, 0.6, 0.8, 1.0]; % m + +tol_ss = {1.0e-5, 1.0e-13}; % {rtol atol} for steady-state + % problem +tol_ts = {1.0e-4, 1.0e-13}; % {rtol atol} for time stepping + +loglevel = 1; % amount of diagnostic output (0 + % to 5) + +refine_grid = 1; % 1 to enable refinement, 0 to + % disable + +%% Create the gas object +% +% This object will be used to evaluate all thermodynamic, kinetic, +% and transport properties + +gas = Solution(rxnmech, 'gri30', 'Mix'); + +% set its state to that of the fuel (arbitrary) +gas.TPX = {tin, p, comp2}; + +%% Create the flow object + +f = AxisymmetricFlow(gas,'flow'); +f.setPressure(p); +f.setupGrid(initial_grid); +f.setSteadyTolerances('default', tol_ss{:}); +f.setTransientTolerances('default', tol_ts{:}); + +%% Create the air inlet +% +% The temperature, mass flux, and composition (relative molar) may be +% specified. + +inlet_o = Inlet('air_inlet'); +inlet_o.T = tin; +inlet_o.setMdot(mdot_o); +inlet_o.setMoleFractions(comp1); + +%% Create the fuel inlet + +inlet_f = Inlet('fuel_inlet'); +inlet_f.T = tin; +inlet_f.setMdot(mdot_f); +inlet_f.setMoleFractions(comp2); + +%% Create the flame object +% +% Once the component parts have been created, they can be assembled +% to create the flame object. + +fl = flame(gas, inlet_o, f, inlet_f); + +% if the starting solution is to be read from a previously-saved +% solution, uncomment this line and edit the file name and solution id. +%restore(fl,'h2flame2.xml', 'energy') + +% solve with fixed temperature profile first +fl.solve(loglevel, refine_grid); + +%% Enable the energy equation +% +% The energy equation will now be solved to compute the +% temperature profile. We also tighten the grid refinement +% criteria to get an accurate final solution. + +f.enableEnergy; +fl.setRefineCriteria(2, 200.0, 0.1, 0.1); +fl.solve(loglevel, refine_grid); +fl.saveSoln('c2h6.xml', 'energy', ['solution with energy equation']); + +%% Show statistics + +fl.writeStats; +elapsed = cputime - t0; +e = sprintf('Elapsed CPU time: %10.4g',elapsed); +disp(e); + +%% Make plots + +figure(1); +subplot(2, 3, 1); +plotSolution(fl, 'flow', 'T'); +title('Temperature [K]'); +subplot(2, 3, 2); +plotSolution(fl, 'flow', 'C2H6'); +title('C2H6 Mass Fraction'); +subplot(2, 3, 3); +plotSolution(fl, 'flow', 'O2'); +title('O2 Mass Fraction'); +subplot(2, 3, 4); +plotSolution(fl, 'flow', 'CH'); +title('CH Mass Fraction'); +subplot(2, 3, 5); +plotSolution(fl, 'flow', 'V'); +title('Radial Velocity / Radius [s^-1]'); +subplot(2, 3, 6); +plotSolution(fl, 'flow', 'u'); +title('Axial Velocity [m/s]'); diff --git a/samples/matlab_experimental/ignite.m b/samples/matlab_experimental/ignite.m index 296598c88..49a65a89e 100644 --- a/samples/matlab_experimental/ignite.m +++ b/samples/matlab_experimental/ignite.m @@ -16,20 +16,21 @@ function plotdata = ignite(g) % set the initial conditions - gas.TPX = {1001.0, oneatm, 'X','H2:2,O2:1,N2:4'}; + gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'}; gas.basis = 'mass'; y0 = [gas.U 1.0/gas.D gas.Y']; time_interval = [0 0.001]; - options = odeset('RelTol',1.e-5,'AbsTol',1.e-12,'Stats','on'); + options = odeset('RelTol', 1.e-5, 'AbsTol', 1.e-12, 'Stats', 'on'); t0 = cputime; - out = ode15s(@reactor_ode,time_interval,y0,options,gas,@vdot,@area,@heatflux); + out = ode15s(@reactor_ode, time_interval, y0, options, gas, ... + @vdot, @area, @heatflux); disp(['CPU time = ' num2str(cputime - t0)]); - plotdata = output(out,gas); + plotdata = output(out, gas); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % the functions below may be defined arbitrarily to set the reactor @@ -71,7 +72,7 @@ function plotdata = ignite(g) gas.TP = {1001.0, oneatm}; for j = 1:n - ss = soln(:,j); + ss = soln(:, j); y = ss(3:end); mass = sum(y); u_mass = ss(1)/mass; @@ -79,11 +80,11 @@ function plotdata = ignite(g) gas.Y = y; gas.UV = {u_mass, v_mass}; - pv(1,j) = times(j); - pv(2,j) = gas.T; - pv(3,j) = gas.D; - pv(4,j) = gas.P; - pv(5:end,j) = y; + pv(1, j) = times(j); + pv(2, j) = gas.T; + pv(3, j) = gas.D; + pv(4, j) = gas.P; + pv(5:end, j) = y; end % plot the temperature and OH mass fractions. @@ -95,7 +96,7 @@ function plotdata = ignite(g) figure(2); ioh = gas.speciesIndex('OH'); - plot(pv(1,:),pv(4+ioh,:)); + plot(pv(1, :), pv(4+ioh, :)); xlabel('time'); ylabel('Mass Fraction'); title('OH Mass Fraction'); diff --git a/samples/matlab_experimental/ignite_hp.m b/samples/matlab_experimental/ignite_hp.m new file mode 100644 index 000000000..915bef79a --- /dev/null +++ b/samples/matlab_experimental/ignite_hp.m @@ -0,0 +1,38 @@ +function ignite_hp(gas) + % IGNITE_HP Solves the same ignition problem as 'ignite', but uses + % function conhp instead of reactor. + + help ignite_hp + + if nargin == 0 + gas = Solution('gri30.yaml'); + end + + mw = gas.MolecularWeights; + gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'}; + + y0 = [gas.T + gas.X']; + tel = [0, 0.001]; + options = odeset('RelTol', 1.e-5, 'AbsTol', 1.e-12, 'Stats', 'on'); + t0 = cputime; + out = ode15s(@conhp, tel, y0, options, gas, mw); + disp(['CPU time = ' num2str(cputime - t0)]); + + if nargout == 0 + % plot the temperature and OH mole fractions. + figure(1); + plot(out.x, out.y(1,:)); + xlabel('time'); + ylabel('Temperature'); + title(['Final T = ' num2str(out.y(1, end)), ' K']); + + figure(2); + ioh = gas.speciesIndex('OH'); + plot(out.x, out.y(1+ioh, :)); + xlabel('time'); + ylabel('Mass Fraction'); + title('OH Mass Fraction'); + end + +end diff --git a/samples/matlab_experimental/ignite_uv.m b/samples/matlab_experimental/ignite_uv.m new file mode 100644 index 000000000..a475f442b --- /dev/null +++ b/samples/matlab_experimental/ignite_uv.m @@ -0,0 +1,38 @@ +function ignite_uv(gas) + % IGNITE_UV Solves the same ignition problem as 'ignite2', except + % function conuv is used instead of reactor. + % + help ignite_uv + + if nargin == 0 + gas = Solution('gri30.yaml'); + end + + mw = gas.MolecularWeights; + gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'}; + + y0 = [gas.T + gas.X']; + tel = [0, 0.001]; + options = odeset('RelTol', 1.e-5, 'AbsTol', 1.e-12, 'Stats', 'on'); + t0 = cputime; + out = ode15s(@conuv, tel, y0, options, gas, mw); + disp(['CPU time = ' num2str(cputime - t0)]); + + if nargout == 0 + % plot the temperature and OH mole fractions. + figure(1); + plot(out.x, out.y(1, :)); + xlabel('time'); + ylabel('Temperature'); + title(['Final T = ' num2str(out.y(1, end)), ' K']); + + figure(2); + ioh = gas.speciesIndex('OH'); + plot(out.x, out.y(1+ioh, :)); + xlabel('time'); + ylabel('Mass Fraction'); + title('OH Mass Fraction'); + end + +end diff --git a/samples/matlab_experimental/prandtl1.m b/samples/matlab_experimental/prandtl1.m index 8423b081e..3c6169f60 100644 --- a/samples/matlab_experimental/prandtl1.m +++ b/samples/matlab_experimental/prandtl1.m @@ -7,6 +7,11 @@ function prandtl1(g) help prandtl1 + clear all + close all + cleanup + clc + if nargin == 1 gas = g; else diff --git a/samples/matlab_experimental/prandtl2.m b/samples/matlab_experimental/prandtl2.m index 9ca49bbc7..5f2e62dec 100644 --- a/samples/matlab_experimental/prandtl2.m +++ b/samples/matlab_experimental/prandtl2.m @@ -6,6 +6,11 @@ function prandtl2(g) % help prandtl2 + clear all + close all + cleanup + clc + if nargin == 1 gas = g; else diff --git a/samples/matlab_experimental/surfreactor.m b/samples/matlab_experimental/surfreactor.m index 87c71340c..a94675016 100644 --- a/samples/matlab_experimental/surfreactor.m +++ b/samples/matlab_experimental/surfreactor.m @@ -21,10 +21,10 @@ gas.TPX = {t, oneatm, 'CH4:0.01, O2:0.21, N2:0.78'}; % methane on platinum, and is from Deutschman et al., 26th % Symp. (Intl.) on Combustion,1996, pp. 1747-1754 surf = importInterface('ptcombust.yaml','Pt_surf', gas); -surf.th.T = t; +surf.T = t; nsp = gas.nSpecies; -nSurfSp = surf.th.nSpecies; +nSurfSp = surf.nSpecies; % create a reactor, and insert the gas r = IdealGasReactor(gas); @@ -60,9 +60,9 @@ nSteps = 100; p0 = r.P; names = {'CH4','CO','CO2','H2O'}; x = zeros([nSteps 4]); -tim = zeros(nSteps); -temp = zeros(nSteps); -pres = zeros(nSteps); +tim = zeros(nSteps, 1); +temp = zeros(nSteps, 1); +pres = zeros(nSteps, 1); cov = zeros([nSteps nSurfSp]); t = 0; dt = 0.1;