1% -------------------------------------------------------------------------
2% EXAMPLE: HP - Humid air
3%
4% Compute adiabatic temperature and equilibrium composition at constant
5% pressure (e.g., 1.01325 bar) for lean to rich CH4-ideal_air mixtures at
6% standard conditions, a set of 10 species considered, and excess of air
7% of 15%, and a set of relative humidity (0, 20, 40, 60, 80, 100) [%]
8%
9% LS == {'CO2', 'CO', 'H2O', 'H2', 'O2', 'N2', 'NO', 'OH', 'O', 'H'}
10%
11% See wiki or setListspecies method from ChemicalSystem class for more
12% predefined sets of species
13%
14% Results compared with [1]
15%
16% [1] Sánchez, Y. A. C., Piedrahita, C. A. R., & Serrano, E. G. F. (2014).
17% Influencia del aire húmedo en la combustión del metano. Scientia et
18% technica, 19(4), 364-370.
19%
20% @author: Alberto Cuadra Lara
21%
22% Last update Aug 05 2024
23% -------------------------------------------------------------------------
24
25% Import packages
26import combustiontoolbox.databases.NasaDatabase
27import combustiontoolbox.core.*
28import combustiontoolbox.equilibrium.*
29import combustiontoolbox.utils.display.*
30import combustiontoolbox.utils.findIndex
31import combustiontoolbox.utils.extensions.brewermap
32
33% User inputs
34airExcess = 15; % [%]
35humidityRelative = 0:20:100; % [%]
36T = [10:2:60] + 273.15; % [K]
37p = 1 * 1.01325; % [bar]
38
39% Get Nasa database
40DB = NasaDatabase();
41
42% Define chemical system
43system = ChemicalSystem(DB, {'CO2', 'CO', 'H2O', 'H2', 'O2', 'N2', 'NO', 'OH', 'O', 'H'});
44
45% Initialize mixture
46mix = Mixture(system);
47
48% Initialize solver
49solver = EquilibriumSolver('problemType', 'HP', 'FLAG_TIME', false);
50
51% Definitions
52stoichiometricMoles = 2; % stoichiometric number of moles of air per mole of fuel
53equivalenceRatio = 1 / (1 + airExcess / 100); % equivalence ratio [-]
54W_air = DB.species.Air.W; % [kg/mol]
55W_H2O = DB.species.H2O.W; % [kg/mol]
56moles_air = 4.76; % total number of moles of stoichiometric ideal dry air (79% N2 and 21% O2 in volume)
57
58% Miscellaneous
59FLAG_FIRST = true;
60config = PlotConfig();
61config.innerposition = [0.2 0.2 0.35 0.5]; % Set figure inner position [normalized]
62config.outerposition = [0.2 0.2 0.35 0.5]; % Set figure outer position [normalized]
63colors = brewermap(length(humidityRelative), 'Greys');
64ax1 = setFigure(config);
65ax2 = setFigure(config);
66ax3 = setFigure(config);
67
68for j = length(humidityRelative):-1:1
69 % Get specific humidity
70 W = humiditySpecific(T, p, humidityRelative(j));
71
72 % Get moles H2O
73 moles_H20 = W * moles_air * (W_air / W_H2O);
74
75 % Solve problem
76 for i = length(T):-1:1
77 % Initialize mixture
78 mix = Mixture(system);
79
80 % Define chemical state
81 set(mix, {'CH4'}, 'fuel', 1);
82 set(mix, {'N2', 'O2', 'H2O'}, 'oxidizer', stoichiometricMoles ./ equivalenceRatio .* [79/21, 1, moles_H20(i)]);
83
84 % Define properties
85 mix1 = setProperties(mix, 'temperature', T(i), 'pressure', p, 'equivalenceRatio', equivalenceRatio);
86 mix2 = setProperties(mix, 'temperature', T(i), 'pressure', p, 'equivalenceRatio', equivalenceRatio);
87
88 % Solve problem
89 solver.solveArray(mix2);
90
91 if FLAG_FIRST
92 index_CO2 = findIndex(system.listSpecies, 'CO2');
93 index_H2O = findIndex(system.listSpecies, 'H2O');
94 index_CH4 = findIndex(system.listSpecies, 'CH4');
95 FLAG_FIRST = false;
96 end
97
98 % Extract data
99 mR(i, j) = mix1.mi;
100 TR(i, j) = mix1.T;
101 Yi_R_CH4(i, j) = mix1.Yi(index_CH4);
102
103 mP(i, j) = mix2.mi;
104 TP(i, j) = mix2.T;
105 Yi_P_H2O(i, j) = mix2.Yi(index_H2O, :);
106 Xi_P_CO2(i, j) = mix2.Xi(index_CO2, :);
107
108 m_CH4(i, j) = mR(i, j) .* Yi_R_CH4(i, j);
109 m_H20(i, j) = mP(i, j) .* Yi_P_H2O(i, j);
110 end
111
112 % Plot results
113
114 % Adiabatic flame temperature vs. reactant's temperature
115 ax1 = plotFigure('T_R\ [\rm K]', TR(:, j), 'T_P\ [\rm K]', TP(:, j), 'ax', ax1, 'color', colors(j, :));
116 % Molar fraction of CO2 vs. reactant's temperature
117 ax2 = plotFigure('T_R\ [\rm K]', TR(:, j), 'X_{\rm CO_2}', Xi_P_CO2(:, j), 'ax', ax2, 'color', colors(j, :));
118 % Water content after the combustion process (kg H2O / kg fuel) vs. reactant's temperature
119 ax3 = plotFigure('T_R\ [\rm K]', TR(:, j), '\rm Water\ content\ [kg_{H_2O}/kg_{CH_4}]', m_H20(:, j)./m_CH4(:, j), 'ax', ax3, 'color', colors(j, :));
120end