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, 'drive_factor', 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
41for i = length(self.PS.strP):-1:1
42 postshock_u2n(:, i) = self.PS.strP{i}.polar.un;
43 postshock_u2(:, i) = self.PS.strP{i}.polar.u;
44 postshock_a(:, i) = self.PS.strP{i}.polar.sound;
45 beta(:, i) = self.PS.strP{i}.polar.beta;
46 theta(:, i) = self.PS.strP{i}.polar.theta;
47end
48
49M2 = postshock_u2 ./ postshock_a;
50M2n = postshock_u2n ./ postshock_a;
51
52for i = length(self.PS.strP):-1:1
53 ind_row_sonic(i) = find(M2(:, i) >= 1, 1);
54 ind_row_sonic_n(i) = find(M2n(:, i) >= 1, 1);
55 beta_sonic_2(i) = beta(ind_row_sonic(i), i);
56 theta_sonic_2(i) = theta(ind_row_sonic(i), i);
57end
58%% DISPLAY RESULTS (PLOTS)
59post_results(self);
60%% SMOOTH RESULTS - FOURIER
61start_point_sonic = [0 0 0 0 0 0.0508682282375078];
62start_point_sonic_2 = [0 0 0 0 0 0.0661385531723257];
63start_point_max = [0 0 0 0 0 0.0656724073958934];
64[beta_sonic_smooth, theta_sonic_smooth] = smooth_data([90, beta_sonic], [0, theta_sonic], start_point_sonic);
65[theta_max_smooth, beta_max_smooth] = smooth_data([0, theta_max], [90, beta_max], start_point_max);
66[theta_sonic_2_smooth, beta_sonic_2_smooth] = smooth_data([0, theta_sonic_2], [90, beta_sonic_2], start_point_sonic_2);
67%% PLOT
68f = figure(2); ax = gca;
69plot_figure('$\theta_{\rm M2n sonic}$', theta_sonic_smooth, '$\beta_{\rm M2n sonic}$', beta_sonic_smooth, 'linestyle', 'k:', 'color', 'auto', 'ax', ax)
70plot_figure('$\theta_{\rm max}$', theta_max_smooth, '$\beta_{\rm max}$', beta_max_smooth, 'linestyle', 'k:', 'color', 'auto', 'ax', ax)
71plot_figure('$\theta_{\rm M2 sonic}$', theta_sonic_2_smooth, '$\beta_{\rm M2 sonic}$', beta_sonic_2_smooth, 'linestyle', 'k:', 'color', 'auto', 'ax', ax)