Incident detonation polar diagrams for different drive factors#

 1% -------------------------------------------------------------------------
 2% EXAMPLE: DET_POLAR
 3%
 4% Compute detonation polar curves for a stoichiometric H2-air mixture at 
 5% standard conditions (T = 300 K, p = 1 atm), and a set of pre-shock Mach
 6% numbers M1 = [5, 6, 7, 8, 10], or what is the same, drive factors 
 7% M1/M1_cj = [1.0382, 1.2458, 1.4534, 1.6611, 2.0763]
 8%
 9% See wiki or setListspecies method from ChemicalSystem class for predefined
10% sets of species
11%
12% @author: Alberto Cuadra Lara
13%
14% Last update October 06 2025
15% -------------------------------------------------------------------------
16
17% Import packages
18import combustiontoolbox.databases.NasaDatabase
19import combustiontoolbox.core.*
20import combustiontoolbox.shockdetonation.*
21
22% Get Nasa database
23DB = NasaDatabase();
24
25% Define chemical system
26system = ChemicalSystem(DB);
27
28% Initialize mixture
29mix = Mixture(system);
30
31% Define chemical state
32set(mix, {'H2'}, 'fuel', 1);
33set(mix, {'N2', 'O2'}, 'oxidizer', [79, 21] / 21);
34
35% Define properties
36mixArray1 = setProperties(mix, 'temperature', 300, 'pressure', 1.01325, 'equivalenceRatio', 1, 'driveFactor', [1.0382, 1.2458, 1.4534, 1.6611, 2.0763]);
37
38% Initialize solver
39solver = DetonationSolver('problemType', 'DET_POLAR', 'numPointsPolar', 300);
40
41% Solve problem
42[mixArray1, mixArray2] = solver.solveArray(mixArray1);
43
44% Generate report
45report(solver, mixArray1, mixArray2);