Example_DET_POLAR_SONIC_AND_MAX.m#

 1% -------------------------------------------------------------------------
 2% EXAMPLE: DET_POLAR_SONIC_AND_MAX
 3%
 4% Compute detonation polar curves at standard conditions, a set of 30 species
 5% considered, and a set of pre-shock Mach numbers M1 = [1.01*M1cj:15], or
 6% what is the same, drive factors M1/M1_cj = [1.01:2.9069]
 7%    
 8% Hydrogen == {'H2O','H2','O2','N2','He','Ar','H','HNO',...
 9%              'HNO3','NH','NH2OH','NO3','N2H2','N2O3','N3','OH',...
10%              'HNO2','N','NH3','NO2','N2O','N2H4','N2O5','O','O3',...
11%              'HO2','NH2','H2O2','N3H','NH2NO2'}
12%   
13% See wiki or list_species() for more predefined sets of species
14%
15% @author: Alberto Cuadra Lara
16%          PhD Candidate - Group Fluid Mechanics
17%          Universidad Carlos III de Madrid
18%                 
19% Last update Oct 07 2022
20% -------------------------------------------------------------------------
21
22%% INITIALIZE
23self = App('Hydrogen');
24%% INITIAL CONDITIONS
25self = set_prop(self, 'TR', 300, 'pR', 1 * 1.01325, 'phi', 1);
26self.PD.S_Fuel = {'H2'};
27self.PD.S_Oxidizer = {'N2', 'O2'};
28self.PD.ratio_oxidizers_O2 = [79, 21]/21;
29%% ADDITIONAL INPUTS (DEPENDS OF THE PROBLEM SELECTED)
30self = set_prop(self, 'overdriven', 1.01:0.01:2.9069);
31%% TUNING PROPERTIES
32self.TN.N_points_polar = 300;  % Number of points to compute shock polar
33%% SOLVE PROBLEM
34self = solve_problem(self, 'DET_POLAR');
35%% GET POINTS
36theta_sonic = cell2vector(self.PS.strP, 'theta_sonic');
37beta_sonic = cell2vector(self.PS.strP, 'beta_sonic');
38theta_max = cell2vector(self.PS.strP, 'theta_max');
39beta_max = cell2vector(self.PS.strP, 'beta_max');
40%% DISPLAY RESULTS (PLOTS)
41post_results(self);
42%% SMOOTH RESULTS - FOURIER
43start_point_sonic = [0 0 0 0 0 0.0508682282375078];
44start_point_max = [0 0 0 0 0 0.0656724073958934];
45[beta_sonic_smooth, theta_sonic_smooth] = smooth_data([90, beta_sonic], [0, theta_sonic], start_point_sonic);
46[theta_max_smooth, beta_max_smooth] = smooth_data([0, theta_max], [90, beta_max], start_point_max);
47%% PLOT
48f = figure(2); ax = gca;
49plot_figure('$\theta_{\rm cj}$', theta_sonic_smooth, '$\beta_{\rm cj}$', beta_sonic_smooth, 'linestyle', 'k:', 'color', 'auto')
50plot_figure('$\theta_{\rm max}$', theta_max_smooth, '$\beta_{\rm max}$', beta_max_smooth, 'linestyle', 'k:', 'color', 'auto')