1% -------------------------------------------------------------------------
2% EXAMPLE: HP - EFFECT OF PRESSURE
3%
4% Compute adiabatic temperature and equilibrium composition at constant
5% pressure (e.g., 1 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.5) [-], and two different sets of species "Complete" or "Incomplete"
8% The incomplete combustion case is evaluated at different pressures to see
9% Le Chatelier's principle, i.e., an increase of pressure shifts
10% equilibrium to the side of the reaction with fewer number of moles of
11% gas (higher pressures implie less dissociation). A comparison of the
12% adiabatic flame temperature, adiabatic index, and Gibbs free energy is
13% included.
14%
15% * Complete:
16% - lean = {'CO2', 'H2O', 'N2', 'Ar', 'O2'}; (equivalence ratio < 1)
17% - rich = {'CO2', 'H2O', 'N2', 'Ar', 'CO', 'H2'}; (equivalence ratio > 1)
18% - soot = {'N2', 'Ar', 'CO', 'H2', 'Cbgrb', 'CO2', 'H2O'}; (equivalence ratio > equivalence ratio soot)
19%
20% * Incomplete:
21% - Soot formation == {'CO2','CO','H2O','H2','O2','N2','He','Ar','Cbgrb',...
22% 'C2','C2H4','CH','CH','CH3','CH4','CN','H',...
23% 'HCN','HCO','N','NH','NH2','NH3','NO','O','OH'}
24%
25% See wiki or list_species() for more predefined sets of species
26%
27% @author: Alberto Cuadra Lara
28% PhD Candidate - Group Fluid Mechanics
29% Universidad Carlos III de Madrid
30%
31% Last update July 29 2022
32% -------------------------------------------------------------------------
33
34%% COMPLETE COMBUSTION
35% Initialization
36self = App('Complete');
37% Set fuel composition
38self.PD.S_Fuel = {'CH4', 'C2H6', 'C3H8', 'C4H10_isobutane', 'H2'};
39self.PD.N_Fuel = [0.8, 0.05, 0.05, 0.05, 0.05];
40% Set oxidizer composition
41self = set_air(self, false);
42% Set temperature, pressure and equivalence ratio
43self = set_prop(self, 'TR', 300, 'pR', 1, 'phi', 0.5:0.005:3.5);
44% Set constant pressure for products
45self = set_prop(self, 'pP', self.PD.pR.value);
46% Solve Problem
47self = solve_problem(self, 'HP');
48% Save results
49results_complete = self;
50%% INCOMPLETE COMBUSTION
51pressure_vector = self.PD.pR.value * [1, 10, 100];
52for i = length(pressure_vector):-1:1
53 % Set product species at equilibrium
54 self = App('copy', self, 'Soot formation');
55 % Set pressure
56 self = set_prop(self, 'pP', pressure_vector(i), 'pP', pressure_vector(i));
57 % Solve Problem
58 self = solve_problem(self, 'HP');
59 % Save results
60 results_incomplete{i} = self;
61end
62%% Comparison of adiabatic flame temperature
63mix1 = results_complete.PS.strR; % Is the same as the incomplete case (same initial mixture)
64mix2_complete = results_complete.PS.strP;
65
66phi = cell2vector(mix1, 'phi');
67
68self.Misc.config.title = '\rm{Complete\ vs\ Incomplete\ at\ different\ pressures}';
69self.Misc.config.labelx = 'Equivalence ratio $\phi$';
70self.Misc.config.labely = 'Temperature $T$ [K]';
71legend_name = {'Complete'};
72
73ax = plot_figure('phi', phi, 'T', mix2_complete, 'config', self.Misc.config);
74
75for i = 1:length(pressure_vector)
76 mix2_incomplete = results_incomplete{i}.PS.strP;
77 ax = plot_figure('phi', phi, 'T', mix2_incomplete, 'config', self.Misc.config, 'ax', ax, 'color', 'auto');
78 legend_name(i+1) = {sprintf('$p_{%d} = %.4g$ bar', i, pressure(mix2_incomplete{1}))};
79end
80
81set_legends(ax, legend_name, 'FLAG_SPECIES', false)
82ax.Legend.Location = 'northeast';
83%% Comparison of the adiabatic index
84self.Misc.config.labely = 'Adibatic index $\gamma_s$';
85
86ax = plot_figure('phi', phi, 'gamma_s', mix2_complete, 'config', self.Misc.config);
87
88for i = 1:length(pressure_vector)
89 mix2_incomplete = results_incomplete{i}.PS.strP;
90 ax = plot_figure('phi', phi, 'gamma_s', mix2_incomplete, 'config', self.Misc.config, 'ax', ax, 'color', 'auto');
91end
92
93set_legends(ax, legend_name, 'FLAG_SPECIES', false)
94ax.Legend.Location = 'best';
95%% Comparison of the Gibbs energy
96self.Misc.config.labely = 'Gibbs free energy $g$ [kJ/kg]';
97
98ax = plot_figure('phi', phi, 'g', mix2_complete, 'config', self.Misc.config);
99
100for i = 1:length(pressure_vector)
101 mix2_incomplete = results_incomplete{i}.PS.strP;
102 ax = plot_figure('phi', phi, 'g', mix2_incomplete, 'config', self.Misc.config, 'ax', ax, 'color', 'auto');
103end
104
105set_legends(ax, legend_name, 'FLAG_SPECIES', false)
106ax.Legend.Location = 'best';