diff --git a/doc/SConscript b/doc/SConscript index e9aece71b..59a026df1 100644 --- a/doc/SConscript +++ b/doc/SConscript @@ -207,7 +207,8 @@ if localenv['sphinx_docs']: '1D/@Domain1D': ['1D/AxiStagnFlow.m', '1D/AxisymmetricFlow.m', '1D/Inlet.m', '1D/Outlet.m', '1D/OutletRes.m', '1D/Surface.m', '1D/SymmPlane.m'], - '1D/@Stack': ['1D/FreeFlame.m', '1D/npflame_init.m'], + '1D/@Stack': ['1D/FreeFlame.m', '1D/npflame_init.m', + '1D/CounterFlowDiffusionFlame.m'], '@Interface': ['importEdge.m', 'importInterface.m'], '@Data': ['Air.m', 'gasconstant.m', 'GRI30.m', 'Hydrogen.m', 'Methane.m', 'Nitrogen.m', 'oneatm.m', diff --git a/interfaces/matlab/toolbox/1D/CounterFlowDiffusionFlame.m b/interfaces/matlab/toolbox/1D/CounterFlowDiffusionFlame.m new file mode 100644 index 000000000..b485e1019 --- /dev/null +++ b/interfaces/matlab/toolbox/1D/CounterFlowDiffusionFlame.m @@ -0,0 +1,220 @@ +function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidizer) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% COUNTERFLOWDIFFUSIONFLAME Create a counter flow diffusion flame stack. +% flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidizer) +% :param left: +% Object representing the left inlet, which must be +% created using function :mat:func:`Inlet`. +% :param flow: +% Object representing the flow, created with +% function :mat:func:`AxisymmetricFlow`. +% :param right: +% Object representing the right inlet, which must be +% created using function :mat:func:`Inlet`. +% :param tp_f: +% Object representing the fuel inlet gas, instance of class +% :mat:func:`Solution`, and an ideal gas. +% :param tp_o: +% Object representing the oxidizer inlet gas, instance of class +% :mat:func:`Solution`, and an ideal gas. +% :param oxidizer: +% String representing the oxidizer species. Most commonly O2. +% :return: +% Instance of :mat:func:`Stack` object representing the left +% inlet, flow, and right inlet. +% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Check input parameters +% + +if nargin ~= 6 + error('CounterFlowDiffusionFlame expects six input arguments.'); +end +if ~isIdealGas(tp_f) + error('Fuel gas object must represent an ideal gas mixture.'); +end +if ~isIdealGas(tp_o) + error('Oxidizer gas object must represent an ideal gas mixture.'); +end +if ~isInlet(left) + error('Left inlet object of wrong type.'); +end +if ~isFlow(flow) + error('Flow object of wrong type.'); +end +if ~isInlet(right) + error('Right inlet object of wrong type.'); +end +if ~ischar(oxidizer) + error('Oxidizer name must be of format character.'); +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Get the density of both fuel and oxidizer streams. To be used in +% determining velocity of each stream. Also get the temperature of both +% inlet streams. +% + +rhof = density(tp_f); +rho0 = density(tp_o); +tf = temperature(left); +tox = temperature(right); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Find the species index of the oxidizer. To be used in determining initial +% strain rate. +% + +ioxidizer = speciesIndex(tp_o, oxidizer); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Calculate the stoichiometric mixture fraction. Needed for determining +% location of flame edges and composition. elMoles function used to +% calculate the number of moles of C, H, and O atoms in the fuel and +% oxidizer streams: elMoles = elementalMassFraction/element atomic weight. +% From this, the stoichiometric Air/Fuel ratio can be determined. +% 1 Mole of O needs 2 Moles of C and 0.5 Moles of H for stoichiometric +% 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'); +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); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Compute the stoichiometric mass fractions of each species. Use this to +% set the fuel gas object and calculate adiabatic flame temperature and +% equilibrium composition. +% + +spec = speciesNames(tp_f); % Get all of the species names in gas object. +nsp = nSpecies(tp_f); % Get total number of species in gas object. +% Get the current mass fractions of both fuel and inlet streams. +yox = massFractions(tp_o); +yf = massFractions(tp_f); +ystoich_double = zeros(1, nsp); % Create empty vector for stoich mass frac. + +for n = 1:nsp + % Calculate stoichiometric mass fractions. + ystoich_double(n) = zst*yf(n) + (1.0 - zst)*yox(n); + % Convert mass fraction vector to string vector. + ystoich_str{n} = num2str(ystoich_double(n)); + % Convert string vector to cell with SPECIES:MASS FRACTION format. + y_stoich{n} = [spec{n}, ':', ystoich_str{n}]; +end +% 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}]; +end +% Set the fuel gas object as stoichiometric values and use equilibrate +% function to determine stoichiometric equilibrium temperature and mass +% fractions. +set(tp_f, 'T', tf, 'P', pressure(tp_f), 'Y', ystoich); +equilibrate(tp_f, 'HP'); +teq = temperature(tp_f); +yeq = massFractions(tp_f); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Estimate the strain rate based on the inlet stream velocities and +% determine initial "guess" for mixture fraction based on mass flux ratio. +% + +zz = gridPoints(flow); +dz = zz(end) - zz(1); +uleft = massFlux(left)/rhof; +uright = massFlux(right)/rho0; +a = (abs(uleft) + abs(uright))/dz; +diff = mixDiffCoeffs(tp_f); +f = sqrt(a/(2.0*diff(ioxidizer))); +x0num = sqrt(uleft*massFlux(left))*dz; +x0den = sqrt(uleft*massFlux(left)) + sqrt(uright*massFlux(right)); +x0 = x0num/x0den; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Calculate initial values of temperature and mass fraction of species in +% flame at each gridpoint. These values to be used for energy equation +% solution. Method is based on the Burke-Schumann model. +% + +nz = nPoints(flow); +zm = zeros(1, nz); +u = zeros(1, nz); +v = zeros(1, nz); +y = zeros(nz, nsp); +t = zeros(1, nz); +for j = 1:nz + x = zz(j); + zeta = f*(x - x0); + zmix = 0.5*(1.0 - erf(zeta)); % Mixture fraction in flame. + 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); + 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; + end + t(j) = tox + zmix*(teq - tox)/zst; + end +end +zrel = zz/dz; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Create the flame stack with the fuel inlet, flow object, and oxidizer +% inlet. Set the profile of the flame with the estimated axial velocities, +% radial velocities, temperature, and mass fractions calculated above. +% + +flame = Stack([left flow right]); +setProfile(flame, 2, {'u', 'V'}, [zrel; u; v]); +setProfile(flame, 2, 'T', [zrel; t] ); +for n = 1:nsp + nm = speciesName(tp_f, n); + setProfile(flame, 2, nm, [zrel; transpose(y(:,n))]) +end +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Define elMoles function +% + +function moles = elMoles(tp, element) +% ELMOLES Determine the elemental moles in a gas object per unit mass. +% moles = Moles(tp, element) +% :param tp: +% Object representing the gas, instance of class :mat:func:`Solution`, +% and an ideal gas. The state of this object should be set to an +% estimate of the gas state before calling Moles. +% :param element: +% String representing the element name. +% :return: +% Elemental moles within a gas object per unit mass of mixture. +% Units: kmol/kg +% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Check input parameters +% + +if nargin ~= 2 + error('elMoles expects two input arguments.'); +end +if ~isIdealGas(tp) + error('Gas object must represent an ideal gas mixture.'); +end +if ~ischar(element) + error('Element name must be of format character.'); +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Calculate the moles per mass of mixture of an element within a gas +% object. The equation used is: elmoles = elMassFrac/Mel where elMassFrac +% is the elemental mass fraction within the gas object using the +% elementalMassFraction function; Mel is the atomic mass of the element. +% + +elMassFrac = elementalMassFraction(tp, element); +eli = elementIndex(tp, element); +M = atomicMasses(tp); +Mel = M(eli); +moles = elMassFrac/Mel; +end diff --git a/interfaces/matlab/toolbox/@ThermoPhase/elementalMassFraction.m b/interfaces/matlab/toolbox/@ThermoPhase/elementalMassFraction.m new file mode 100644 index 000000000..a2ecbabb2 --- /dev/null +++ b/interfaces/matlab/toolbox/@ThermoPhase/elementalMassFraction.m @@ -0,0 +1,51 @@ +function elMassFrac = elementalMassFraction(tp, element) +% ELEMENTALMASSFRACTION Determine the elemental mass fraction in gas object. +% elMassFrac = elementalMassFraction(tp, element) +% :param tp: +% Object representing the gas, instance of class :mat:func:`Solution`, +% and an ideal gas. The state of this object should be set to an +% estimate of the gas state before calling elementalMassFraction. +% :param element: +% String representing the element name. +% :return: +% Elemental mass fraction within a gas object. +% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Check input parameters +% + +if nargin ~= 2 + error('elementalMassFraction expects two input arguments.'); +end +if ~isIdealGas(tp) + error('Gas object must represent an ideal gas mixture.'); +end +if ~ischar(element) + error('Element name must be of format character.'); +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Calculate the elemental mass fraction in a gas object. Equation used is +% elMassFrac = sum of nAtoms(k,m)*Mel(m)*Y(k)/mw(k) where nAtoms(k,m) is +% the number of atoms of element, m, in species, k; Mel(m) is the atomic +% weight of the element, m; Y(k) is the mass fraction of species,k, in the +% gas object; and mw(k) is the molecular weight of species, k. +% + +n = nSpecies(tp); +massFrac = massFractions(tp); +spec = speciesNames(tp); +eli = elementIndex(tp, element); +M = atomicMasses(tp); +Mel = M(eli); +MW = molecularWeights(tp); +% Initialize the element mass fraction as zero. +elMassFrac = 0.0; +% Use loop to perform summation of elemental mass fraction over all species. +for i = 1:n + natoms(i) = nAtoms(tp,spec{i},element); + mw(i) = MW(i); + Y(i) = massFraction(tp,spec{i}); + elMassFrac = elMassFrac + (natoms(i)*Mel*Y(i))/mw(i); +end +end diff --git a/samples/matlab/diffflame.m b/samples/matlab/diffflame.m index 022b11e43..78bae086e 100644 --- a/samples/matlab/diffflame.m +++ b/samples/matlab/diffflame.m @@ -1,119 +1,129 @@ -% DIFFFLAME - A non-premixed opposed-jet flame. -% +% DIFFFLAME - An opposed-flow diffusion flame. +% This example uses the CounterFlowDiffusionFlame function to solve an +% opposed-flow diffusion flame for Ethane in Air. This example is the same +% as the diffusion_flame.py example without radiation. % -help diffflame -disp('press any key to begin the simulation'); -pause +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% -t0 = cputime; % record the starting time +runtime = cputime; % Record the starting time +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Parameter values of inlet streams +% -% 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 mass flux, kg/m^2/s +mdot_f = 0.24; % Fuel mass flux, kg/m^2/s +rxnmech = 'gri30.xml'; % Reaction mechanism file +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. +% -rxnmech = 'gri30.xml'; % reaction mechanism file -transport = 'Mix'; % transport model -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-9]; % [rtol atol] for steady-state - % problem +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 object %%%%%%%%%%%%%%%%%%%%%%%% +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. % -% This object will be used to evaluate all thermodynamic, kinetic, -% and transport properties + +fuel = GRI30('Mix'); +ox = GRI30('Mix'); +oxcomp = 'O2:0.21, N2:0.78'; % Air composition +fuelcomp = 'C2H6:1'; % Fuel composition +% Set each gas mixture state with the corresponding composition. +set(fuel,'T', tin, 'P', p, 'X', fuelcomp); +set(ox,'T',tin,'P',p,'X', oxcomp); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Set-up the flow object. For this problem, the AxisymmetricFlow model is +% needed. Set the state of the flow as the fuel gas object. This is +% arbitrary and is only used to initialize the flow object. Set the grid to +% the initial grid defined prior, same for the tolerances. % -gas = GRI30('Mix') - -% set its state to that of the fuel (arbitrary) -set(gas,'T', tin, 'P', p, 'X', comp2); - -%%%%%%%%%%%%%%%% create the flow object %%%%%%%%%%%%%%%%%%%%%%% - -f = AxisymmetricFlow(gas,'flow'); +f = AxisymmetricFlow(fuel,'flow'); set(f, 'P', p, 'grid', initial_grid); set(f, 'tol', tol_ss, 'tol-time', tol_ts); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Create the fuel and oxidizer inlet steams. Specify the temperature, mass +% flux, and composition correspondingly. +% -%%%%%%%%%%%%%%% create the air inlet %%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% The temperature, mass flux, and composition (relative molar) may be -% specified. -% +% Set the oxidizer inlet. inlet_o = Inlet('air_inlet'); -set(inlet_o, 'T', tin, 'MassFlux', mdot_o, 'X', comp1); - -%%%%%%%%%%%%%% create the fuel inlet %%%%%%%%%%%%%%%%%%%%%%%%%%%% -% +set(inlet_o, 'T', tin, 'MassFlux', mdot_o, 'X', oxcomp); % +% Set the fuel inlet. inlet_f = Inlet('fuel_inlet'); -set(inlet_f, 'T', tin, 'MassFlux', mdot_f, 'X', comp2); - -%%%%%%%%%%%%% create the flame object %%%%%%%%%%%% +set(inlet_f, 'T', tin, 'MassFlux', mdot_f, 'X', fuelcomp); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Once the inlets have been created, they can be assembled +% to create the flame object. Function CounterFlorDiffusionFlame +% (in Cantera/1D) sets up the initial guess for the solution using a +% Burke-Schumann flame. The input parameters are: fuel inlet object, flow +% object, oxidizer inlet object, fuel gas object, oxidizer gas object, and +% the name of the oxidizer species as in character format. % -% Once the component parts have been created, they can be assembled -% to create the flame object. Function npflame_init (in Cantera/1D) -% sets up the initial guess for the solution using a Burke-Schumann -% flame. + +fl = CounterFlowDiffusionFlame(inlet_f, f, inlet_o, fuel, ox, 'O2'); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Solve with fixed temperature profile first. Grid refinement is turned off +% for this process in this example. To turn grid refinement on, change 0 to +% 1 for last input is solve function. % -fl = npflame_init(gas, inlet_f, f, inlet_o, 'C2H6', 'O2', 3.5); -% 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 -solve(fl, loglevel, 0); %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. +solve(fl, loglevel, 0); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Enable the energy equation. The energy equation will now be solved to +% compute the temperature profile. We also tighten the grid refinement +% criteria to get an accurate final solution. The explanation of the +% setRefineCriteria function is located on cantera.org in the Matlab User's +% Guide and can be accessed by help setRefineCriteria % enableEnergy(f); setRefineCriteria(fl, 2, 200.0, 0.1, 0.2); solve(fl, loglevel, refine_grid); saveSoln(fl,'c2h6.xml','energy',['solution with energy equation']); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Show statistics of solution and elapsed time. +% -%%%%%%%%%% show statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% writeStats(fl); -elapsed = cputime - t0; +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. +% -%%%%%%%%%% make plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - +z = grid(fl, 'flow'); % Get grid points of flow +spec = speciesNames(fuel); % Get species names in gas +T = solution(fl, 'flow', 'T'); % Get temperature solution +for i = 1:length(spec) + % Get mass fraction of all species from solution + y(i,:) = solution(fl, 'flow', spec{i}); +end +j = speciesIndex(fuel, 'O2'); % Get index of O2 in gas object +k = speciesIndex(fuel, 'H2O'); % Get index of H2O in gas object +l = speciesIndex(fuel, 'C2H6'); % Get index of C2H6 in gas object +m = speciesIndex(fuel, 'CO2'); % Get index of CO2 in gas object clf; -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]'); +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');