From 43bd96b5ba6073bd952dd2231daac3ead40e1ff7 Mon Sep 17 00:00:00 2001 From: KyleLinevitchJr Date: Thu, 4 May 2017 18:44:02 -0400 Subject: [PATCH] [Matlab] Make counterflow diffusion flame simulation more general The CounterFlowDiffusionFlame (CFDF) code is able to perform more general cases of npflame_init for multiple species fuel and oxidizer streams. The stoichiometric mixture fraction in the CFDF code uses the Bilger definition of mixture fraction, using the conservation of elements C, H, and O. This method is used in the python module, but not the MATLAB npflame_init function. Also, the CFDF code uses the fuel stream density to calculate the fuel stream velocity and the oxidizer stream density to calculate the oxidizer stream velocity, where as the npflame_init code uses the fuel density for both velocity calculations. The elementMassFraction code is a MATLAB version of the python function: elemental_mass_fraction, which is needed to run the CFDF code. Update the diffflame.m example to use the more general CFDF function since the input parameters are different than the npflame_init function. This example is the same as the diffusion_flame.py sample in the Python module. --- doc/SConscript | 3 +- .../toolbox/1D/CounterFlowDiffusionFlame.m | 220 ++++++++++++++++++ .../@ThermoPhase/elementalMassFraction.m | 51 ++++ samples/matlab/diffflame.m | 188 ++++++++------- 4 files changed, 372 insertions(+), 90 deletions(-) create mode 100644 interfaces/matlab/toolbox/1D/CounterFlowDiffusionFlame.m create mode 100644 interfaces/matlab/toolbox/@ThermoPhase/elementalMassFraction.m 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');