Example_SHOCK_POLAR.m#

 1% -------------------------------------------------------------------------
 2% EXAMPLE: SHOCK_POLAR
 3%
 4% Compute shock polar plots at standard conditions, a set of 51 species
 5% considered, and a set of initial shock front Mach numbers = [2, 3, 5, 14]
 6%    
 7% Air_ions == {'eminus', 'Ar', 'Arplus', 'C', 'Cplus', 'Cminus', ...
 8%              'CN', 'CNplus', 'CNminus', 'CNN', 'CO', 'COplus', ...
 9%              'CO2', 'CO2plus', 'C2', 'C2plus', 'C2minus', 'CCN', ...
10%              'CNC', 'OCCN', 'C2N2', 'C2O', 'C3', 'C3O2', 'N', ...
11%              'Nplus', 'Nminus', 'NCO', 'NO', 'NOplus', 'NO2', ...
12%              'NO2minus', 'NO3', 'NO3minus', 'N2', 'N2plus', ...
13%              'N2minus', 'NCN', 'N2O', 'N2Oplus', 'N2O3', 'N2O4', ...
14%              'N2O5', 'N3', 'O', 'Oplus', 'Ominus', 'O2', 'O2plus', ...
15%              'O2minus', 'O3'}
16%   
17% See wiki or setListspecies method from ChemicalSystem class for more
18% predefined sets of species
19%
20% @author: Alberto Cuadra Lara
21%          Postdoctoral researcher - Group Fluid Mechanics
22%          Universidad Carlos III de Madrid
23%                 
24% Last update April 02 2024
25% -------------------------------------------------------------------------
26
27% Import packages
28import combustiontoolbox.databases.NasaDatabase
29import combustiontoolbox.core.*
30import combustiontoolbox.shockdetonation.*
31import combustiontoolbox.utils.display.*
32
33% Get Nasa database
34DB = NasaDatabase();
35
36% Define chemical system
37system = ChemicalSystem(DB, 'air ions');
38
39% Initialize mixture
40mix = Mixture(system);
41
42% Define chemical state
43set(mix, {'N2', 'O2', 'Ar', 'CO2'}, [78.084, 20.9476, 0.9365, 0.0319] / 20.9476);
44
45% Define properties
46mixArray1 = setProperties(mix, 'temperature', 300, 'pressure', 1.01325, 'M1', [2, 3, 5, 14]);
47
48% Initialize solver
49solver = ShockSolver('problemType', 'SHOCK_POLAR');
50
51% Solve problem
52[mixArray1, mixArray2] = solver.solveArray(mixArray1);
53
54% Plot polars
55[ax1, ax2, ax3] = plotPolar(mixArray1, mixArray2);