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% Postdoctoral researcher - Group Fluid Mechanics
31% Universidad Carlos III de Madrid
32%
33% Last update Jun 11 2024
34% -------------------------------------------------------------------------
35
36% Import packages
37import combustiontoolbox.databases.NasaDatabase
38import combustiontoolbox.core.*
39import combustiontoolbox.equilibrium.*
40import combustiontoolbox.shockdetonation.*
41import combustiontoolbox.utils.display.*
42import combustiontoolbox.common.Units
43
44% Definitions
45S_Oxidizer = {'N2', 'O2', 'Ar', 'CO2'};
46N_Oxidizer = [78.084, 20.9476, 0.9365, 0.0319] ./ 20.9476;
47machNumber = 1.75:0.1:20.5;
48z = [0, 15000, 30000];
49numCases = length(machNumber);
50numAltitude = length(z);
51
52% Get Nasa database
53DB = NasaDatabase();
54
55% Define chemical system
56system = ChemicalSystem(DB, 'air_ions');
57
58% Initialize mixture
59mix = Mixture(system);
60
61% Define chemical state
62set(mix, {'N2', 'O2', 'Ar', 'CO2'}, [78.084, 20.9476, 0.9365, 0.0319] / 20.9476);
63
64% Initialization of Mixture objects
65mixArray1 = repmat(mix, [numCases, 3, numAltitude]);
66mixArray2 = repmat(mix, [numCases, 3, numAltitude]);
67mixArray2_1 = repmat(mix, [numCases, 3, numAltitude]);
68mixArray3 = repmat(mix, [numCases, 3, numAltitude]);
69
70% Initialization of EquilibriumSolver
71equilibriumSolver = EquilibriumSolver();
72
73for i = 3:3
74
75 switch i
76 case 1
77 equilibriumSolver.FLAG_FROZEN = false;
78 equilibriumSolver.FLAG_TCHEM_FROZEN = false;
79 case 2
80 equilibriumSolver.FLAG_FROZEN = true;
81 equilibriumSolver.FLAG_TCHEM_FROZEN = false;
82 case 3
83 equilibriumSolver.FLAG_FROZEN = false;
84 equilibriumSolver.FLAG_TCHEM_FROZEN = true;
85 end
86
87 for j = numAltitude
88 % Get temperature and pressure
89 [~, ~, temperature, pressure] = atmos(z(j));
90 pressure = Units.convert(pressure, 'Pa', 'bar');
91
92 % Define properties
93 mixArray1(:, i, j) = setProperties(mix, 'temperature', temperature, 'pressure', pressure, 'M1', machNumber);
94
95 % Initialize solver
96 solver = ShockSolver('problemType', 'SHOCK_POLAR_LIMITRR', 'equilibriumSolver', equilibriumSolver, 'tol0', 1e-7);
97
98 % Solve problem
99 [mixArray1(:, i, j),...
100 mixArray2(:, i, j),...
101 mixArray2_1(:, i, j),...
102 mixArray3(:, i, j)] = solver.solveArray(mixArray1(:, i, j));
103 end
104
105end
106
107% Definitions
108ax = setFigure();
109colors = [175 221 233; 95 188 211; 0 102 128] / 255;
110
111for j = 1:numAltitude
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, :));
113 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, :));
114 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, :));
115end
116
117for j = 1:numAltitude
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, :));
119 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, :));
120 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, :));
121end