1% -------------------------------------------------------------------------
2% EXAMPLE: HP COMPLETE VS INCOMPLETE
3%
4% Compute adiabatic temperature and equilibrium composition at constant
5% pressure (e.g., 1.01325 bar) for lean to rich natural gas-air mixtures at
6% standard conditions, a set of equivalence ratios phi contained in
7% (0.5, 3) [-], and two different sets of species "Complete" or "Incomplete"
8%
9% * Complete:
10% - lean = {'CO2', 'H2O', 'N2', 'Ar', 'O2'}; (equivalence ratio < 1)
11% - rich = {'CO2', 'H2O', 'N2', 'Ar', 'CO', 'H2'}; (equivalence ratio > 1)
12% - soot = {'N2', 'Ar', 'CO', 'H2', 'Cbgrb', 'CO2', 'H2O'}; (equivalence ratio > equivalence ratio soot)
13%
14% * Incomplete:
15% - Soot formation == {'CO2','CO','H2O','H2','O2','N2','Ar','Cbgrb',...
16% 'C2','C2H4','CH','CH3','CH4','CN','H',...
17% 'HCN','HCO','N','NH','NH2','NH3','NO','O','OH'}
18%
19% See wiki or setListspecies method from ChemicalSystem class for more
20% predefined sets of species
21%
22% @author: Alberto Cuadra Lara
23% Postdoctoral researcher - Group Fluid Mechanics
24% Universidad Carlos III de Madrid
25%
26% Last update Jul 25 2024
27% -------------------------------------------------------------------------
28
29% Import packages
30import combustiontoolbox.databases.NasaDatabase
31import combustiontoolbox.core.*
32import combustiontoolbox.equilibrium.*
33import combustiontoolbox.utils.display.*
34
35% Get Nasa database
36DB = NasaDatabase();
37
38% Initialize solver
39solver = EquilibriumSolver('problemType', 'HP');
40
41%% COMPLETE COMBUSTION
42
43% Define chemical system
44system = ChemicalSystem(DB, 'complete');
45
46% Initialize mixture
47mix = Mixture(system);
48
49% Define chemical state
50set(mix, {'CH4', 'C2H6', 'C3H8', 'C4H10_isobutane', 'H2'}, 'fuel', [0.8, 0.05, 0.05, 0.05, 0.05]);
51set(mix, {'N2', 'O2', 'Ar', 'CO2'}, 'oxidizer', [78.084, 20.9476, 0.9365, 0.0319] / 20.9476);
52
53% Define properties
54mixArray1 = setProperties(mix, 'temperature', 300, 'pressure', 1 * 1.01325, 'equivalenceRatio', 0.5:0.02:3);
55
56% Solve problem
57solver.solveArray(mixArray1);
58
59%% INCOMPLETE COMBUSTION
60
61% Define chemical system
62system = ChemicalSystem(DB, 'soot formation');
63
64% Initialize mixture
65mix = Mixture(system);
66
67% Define chemical state
68set(mix, {'CH4', 'C2H6', 'C3H8', 'C4H10_isobutane', 'H2'}, 'fuel', [0.8, 0.05, 0.05, 0.05, 0.05]);
69set(mix, {'N2', 'O2', 'Ar', 'CO2'}, 'oxidizer', [78.084, 20.9476, 0.9365, 0.0319] / 20.9476);
70
71% Define properties
72mixArray2 = setProperties(mix, 'temperature', 300, 'pressure', 1 * 1.01325, 'equivalenceRatio', 0.5:0.02:3);
73
74% Solve problem
75solver.solveArray(mixArray2);
76
77%% COMPARE RESULTS
78ax = solver.report(mixArray1, mixArray2);
79
80% Another possibility is call directly the next functions:
81% ax = plotProperties(repmat({mixArray1(1).rangeName}, 1, 9), mixArray1, {'T', 'rho', 'h', 'e', 'g', 'cp', 's', 'gamma_s', 'sound'}, mixArray1, 'basis', {[], [], 'mi', 'mi', 'mi', 'mi', 'mi', [], []});
82% ax = plotProperties(repmat({mixArray1(1).rangeName}, 1, 9), mixArray1, {'T', 'rho', 'h', 'e', 'g', 'cp', 's', 'gamma_s', 'sound'}, mixArray2, 'basis', {[], [], 'mi', 'mi', 'mi', 'mi', 'mi', [], []}, 'ax', ax);
83% legend(ax.Children(end), {'Complete', 'Incomplete'}, 'Interpreter', 'latex', 'FontSize', ax.Children(end).FontSize);