1% -------------------------------------------------------------------------
2% EXAMPLE: SHOCK_POLAR_REFLECTED_THERMO
3%
4% Compute shock polar plots at T1 = 226.65 K and p1 = 0.0117 bar
5% (altitude 30 km), a set of 28 species considered, an initial
6% shock front Mach number = 20, deflection angle theta = 35, and different
7% thermochemical models (chemical equilibrium and frozen chemistry).
8%
9% LS== {'eminus', 'Ar', 'Arplus', 'N', ...
10% 'Nplus', 'Nminus', 'NO', 'NOplus', 'NO2', ...
11% 'NO2minus', 'NO3', 'NO3minus', 'N2', 'N2plus', ...
12% 'N2minus', 'N2O', 'N2Oplus', 'N2O3', 'N2O4', ...
13% 'N2O5', 'N3', 'O', 'Oplus', 'Ominus', 'O2', 'O2plus', ...
14% 'O2minus', 'O3'}
15%
16% See wiki or setListspecies method from ChemicalSystem class for more
17% predefined sets of species
18%
19% @author: Alberto Cuadra Lara
20% Postdoctoral researcher - Group Fluid Mechanics
21% Universidad Carlos III de Madrid
22%
23% Last update April 02 2024
24% -------------------------------------------------------------------------
25
26% Import packages
27import combustiontoolbox.databases.NasaDatabase
28import combustiontoolbox.core.*
29import combustiontoolbox.shockdetonation.*
30import combustiontoolbox.utils.display.*
31
32% Definitions
33listSpecies = {'eminus', 'Ar', 'Arplus', 'N', ...
34 'Nplus', 'Nminus', 'NO', 'NOplus', 'NO2', ...
35 'NO2minus', 'NO3', 'NO3minus', 'N2', 'N2plus', ...
36 'N2minus', 'N2O', 'N2Oplus', 'N2O3', 'N2O4', ...
37 'N2O5', 'N3', 'O', 'Oplus', 'Ominus', 'O2', 'O2plus', ...
38 'O2minus', 'O3'};
39
40% Get Nasa database
41DB = NasaDatabase();
42
43% Define chemical system
44system = ChemicalSystem(DB, listSpecies);
45
46% Initialize mixture
47mix = Mixture(system);
48
49% Define chemical state
50set(mix, {'N2', 'O2', 'Ar'}, [78, 21, 1] / 21);
51
52% Define properties
53mixArray1 = setProperties(mix, 'temperature', 226.65, 'pressure', 0.0117, 'M1', 20, 'theta', 35);
54
55% Initialize solver
56solver = ShockSolver('problemType', 'SHOCK_POLAR_R');
57
58% Solve problem
59[mixArray1, mixArray2, mixArray2_1, mixArray3, mixArray3_1, mixArray3_2] = solver.solveArray(mixArray1);
60
61% Plot polars - incident
62plotPolar(mixArray1, mixArray2);
63
64% Plot polars - reflected
65plotPolar(mixArray2_1, mixArray3, mixArray2_1, mixArray1);