Example_HP_PRESSURE.m#

  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);
 78    legend_name(i+1) = {sprintf('$p_{%d} = %.4g$ bar', i, pressure(mix2_incomplete{1}))};
 79end
 80
 81set_legends(ax, legend_name)
 82ax.Legend.Location = 'northeast';
 83%% Comparison of the adiabatic index
 84yvariable = 'gamma_s';
 85
 86self.Misc.config.labely = 'Adibatic index $\gamma_s$';
 87
 88ax = plot_figure(phi, mix2_complete, 'phi', yvariable, self.Misc.config, self.PD.CompleteOrIncomplete);
 89
 90for i = 1:length(pressure_vector)
 91    mix2_incomplete = results_incomplete{i}.PS.strP; 
 92    ax = plot_figure(phi, mix2_incomplete, 'phi', yvariable, self.Misc.config, self.PD.CompleteOrIncomplete, ax);
 93end
 94
 95set_legends(ax, legend_name)
 96ax.Legend.Location = 'best';
 97%% Comparison of the Gibbs energy
 98yvariable = 'g';
 99
100self.Misc.config.labely = 'Gibbs free energy $g$ [kJ/kg]';
101
102ax = plot_figure(phi, mix2_complete, 'phi', yvariable, self.Misc.config, self.PD.CompleteOrIncomplete);
103
104for i = 1:length(pressure_vector)
105    mix2_incomplete = results_incomplete{i}.PS.strP; 
106    ax = plot_figure(phi, mix2_incomplete, 'phi', yvariable, self.Misc.config, self.PD.CompleteOrIncomplete, ax);
107end
108
109set_legends(ax, legend_name)
110ax.Legend.Location = 'best';