1% -------------------------------------------------------------------------
2% EXAMPLE: SHOCK_POLAR_LIMIT_RR_ALTITUDE
3%
4% Compute limit regular reflections at different altitudes considering a
5% thermochemical frozen gas, a chemically frozen gas, and dissociation,
6% ionization, vibrational excitation and electronic excitation. The
7% calculations are carried out for a set of pre-shock Mach numbers
8% M1 = [1.75:0.1:20]
9%
10% Air_ions == {'eminus', 'Ar', 'Arplus', 'C', 'Cplus', 'Cminus', ...
11% 'CN', 'CNplus', 'CNminus', 'CNN', 'CO', 'COplus', ...
12% 'CO2', 'CO2plus', 'C2', 'C2plus', 'C2minus', 'CCN', ...
13% 'CNC', 'OCCN', 'C2N2', 'C2O', 'C3', 'C3O2', 'N', ...
14% 'Nplus', 'Nminus', 'NCO', 'NO', 'NOplus', 'NO2', ...
15% 'NO2minus', 'NO3', 'NO3minus', 'N2', 'N2plus', ...
16% 'N2minus', 'NCN', 'N2O', 'N2Oplus', 'N2O3', 'N2O4', ...
17% 'N2O5', 'N3', 'O', 'Oplus', 'Ominus', 'O2', 'O2plus', ...
18% 'O2minus', 'O3'}
19%
20% See wiki or setListspecies method from ChemicalSystem class for more
21% predefined sets of species
22%
23% Note: This example uses the Standard Atmosphere Functions package [1]
24%
25% [1] Sky Sartorius (2024). Standard Atmosphere Functions. Available at
26% https://github.com/sky-s/standard-atmosphere, GitHub.
27% Retrieved June 12, 2024.
28%
29% @author: Alberto Cuadra Lara
30%
31% Last update Jun 11 2024
32% -------------------------------------------------------------------------
33
34% Import packages
35import combustiontoolbox.databases.NasaDatabase
36import combustiontoolbox.core.*
37import combustiontoolbox.equilibrium.*
38import combustiontoolbox.shockdetonation.*
39import combustiontoolbox.utils.display.*
40import combustiontoolbox.common.Units
41
42% Definitions
43S_Oxidizer = {'N2', 'O2', 'Ar', 'CO2'};
44N_Oxidizer = [78.084, 20.9476, 0.9365, 0.0319] ./ 20.9476;
45machNumber = 1.75:0.1:20.5;
46z = [0, 15000, 30000];
47numCases = length(machNumber);
48numAltitude = length(z);
49
50% Get Nasa database
51DB = NasaDatabase();
52
53% Define chemical system
54system = ChemicalSystem(DB, 'air_ions');
55
56% Initialize mixture
57mix = Mixture(system);
58
59% Define chemical state
60set(mix, {'N2', 'O2', 'Ar', 'CO2'}, [78.084, 20.9476, 0.9365, 0.0319] / 20.9476);
61
62% Initialization of Mixture objects
63mixArray1 = repmat(mix, [numCases, 3, numAltitude]);
64mixArray2 = repmat(mix, [numCases, 3, numAltitude]);
65mixArray2_1 = repmat(mix, [numCases, 3, numAltitude]);
66mixArray3 = repmat(mix, [numCases, 3, numAltitude]);
67
68% Initialization of EquilibriumSolver
69equilibriumSolver = EquilibriumSolver();
70
71for i = 3:3
72
73 switch i
74 case 1
75 equilibriumSolver.FLAG_FROZEN = false;
76 equilibriumSolver.FLAG_TCHEM_FROZEN = false;
77 case 2
78 equilibriumSolver.FLAG_FROZEN = true;
79 equilibriumSolver.FLAG_TCHEM_FROZEN = false;
80 case 3
81 equilibriumSolver.FLAG_FROZEN = false;
82 equilibriumSolver.FLAG_TCHEM_FROZEN = true;
83 end
84
85 for j = numAltitude
86 % Get temperature and pressure
87 [~, ~, temperature, pressure] = atmos(z(j));
88 pressure = Units.convert(pressure, 'Pa', 'bar');
89
90 % Define properties
91 mixArray1(:, i, j) = setProperties(mix, 'temperature', temperature, 'pressure', pressure, 'M1', machNumber);
92
93 % Initialize solver
94 solver = ShockSolver('problemType', 'SHOCK_POLAR_LIMITRR', 'equilibriumSolver', equilibriumSolver, 'tol0', 1e-7);
95
96 % Solve problem
97 [mixArray1(:, i, j),...
98 mixArray2(:, i, j),...
99 mixArray2_1(:, i, j),...
100 mixArray3(:, i, j)] = solver.solveArray(mixArray1(:, i, j));
101 end
102
103end
104
105% Definitions
106ax = setFigure();
107colors = [175 221 233; 95 188 211; 0 102 128] / 255;
108
109for j = 1:numAltitude
110 ax = plotFigure('Pre-shock Mach number', [mixArray1(:, i, j).mach], 'Wave angle limit RR [deg]', [mixArray2_1(:, i, j).beta], 'linestyle', 'k-', 'ax', ax, 'color', colors(j, :));
111 ax = plotFigure('Pre-shock Mach number', [mixArray1(:, i, j).mach], 'Wave angle limit RR [deg]', [mixArray2_1(:, i, j).beta], 'linestyle', 'k--', 'ax', ax, 'color', colors(j, :));
112 ax = plotFigure('Pre-shock Mach number', [mixArray1(:, i, j).mach], 'Wave angle limit RR [deg]', [mixArray2_1(:, i, j).beta], 'linestyle', 'k:', 'ax', ax, 'color', colors(j, :));
113end
114
115for j = 1:numAltitude
116 ax = plotFigure('Pre-shock Mach number', [mixArray1(:, i, j).mach], 'Deflection angle limit RR [deg]', [mixArray2_1(:, i, j).theta], 'linestyle', 'k-', 'ax', ax, 'color', colors(j, :));
117 ax = plotFigure('Pre-shock Mach number', [mixArray1(:, i, j).mach], 'Deflection angle limit RR [deg]', [mixArray2_1(:, i, j).theta], 'linestyle', 'k--', 'ax', ax, 'color', colors(j, :));
118 ax = plotFigure('Pre-shock Mach number', [mixArray1(:, i, j).mach], 'Deflection angle limit RR [deg]', [mixArray2_1(:, i, j).theta], 'linestyle', 'k:', 'ax', ax, 'color', colors(j, :));
119end