1% -------------------------------------------------------------------------
2% EXAMPLE: SHOCK_OBLIQUE_THETA
3%
4% Compute pre-shock and post-shock state for a oblique incident shock wave
5% at standard conditions, a set of 51 species considered, a pre-shock Mach
6% number M1 = 10, and a set of deflection angle theta = [5:1:40] [deg]
7%
8% Air_ions == {'eminus', 'Ar', 'Arplus', 'C', 'Cplus', 'Cminus', ...
9% 'CN', 'CNplus', 'CNminus', 'CNN', 'CO', 'COplus', ...
10% 'CO2', 'CO2plus', 'C2', 'C2plus', 'C2minus', 'CCN', ...
11% 'CNC', 'OCCN', 'C2N2', 'C2O', 'C3', 'C3O2', 'N', ...
12% 'Nplus', 'Nminus', 'NCO', 'NO', 'NOplus', 'NO2', ...
13% 'NO2minus', 'NO3', 'NO3minus', 'N2', 'N2plus', ...
14% 'N2minus', 'NCN', 'N2O', 'N2Oplus', 'N2O3', 'N2O4', ...
15% 'N2O5', 'N3', 'O', 'Oplus', 'Ominus', 'O2', 'O2plus', ...
16% 'O2minus', 'O3'}
17%
18% See wiki or setListspecies method from ChemicalSystem class for more
19% predefined sets of species
20%
21% @author: Alberto Cuadra Lara
22% Postdoctoral researcher - Group Fluid Mechanics
23% Universidad Carlos III de Madrid
24%
25% Last update April 02 2024
26% -------------------------------------------------------------------------
27
28% Import packages
29import combustiontoolbox.databases.NasaDatabase
30import combustiontoolbox.core.*
31import combustiontoolbox.shockdetonation.*
32import combustiontoolbox.utils.display.*
33
34% Get Nasa database
35DB = NasaDatabase();
36
37% Define chemical system
38system = ChemicalSystem(DB, 'air ions');
39
40% Initialize mixture
41mix = Mixture(system);
42
43% Define chemical state
44set(mix, {'N2', 'O2', 'Ar', 'CO2'}, [78.084, 20.9476, 0.9365, 0.0319] / 20.9476);
45
46% Define properties
47mixArray1 = setProperties(mix, 'temperature', 300, 'pressure', 1.01325, 'M1', 10, 'theta', 5:1:40);
48
49% Initialize solver
50solver = ShockSolver('problemType', 'SHOCK_OBLIQUE');
51
52% Solve problem
53[mixArray1, mixArray2_1, mixArray2_2] = solver.solveArray(mixArray1);
54
55% Plot Hugoniot curve
56ax1 = plotFigure('\rho_1 / \rho_2', [mixArray1.rho] ./ [mixArray2_1.rho], 'p_2 / p_1', [mixArray2_1.p] ./ [mixArray1.p], 'xScale', 'log', 'yScale', 'log', 'color', 'auto');
57ax1 = plotFigure('\rho_1 / \rho_2', [mixArray1.rho] ./ [mixArray2_2.rho], 'p_2 / p_1', [mixArray2_2.p] ./ [mixArray1.p], 'xScale', 'log', 'yScale', 'log', 'color', 'auto', 'ax', ax1);
58
59% Plot post-shock temperature
60ax2 = plotFigure('theta', [mixArray1.theta], 'T', [mixArray2_1.T], 'color', 'auto');
61ax2 = plotFigure('theta', [mixArray1.theta], 'T', [mixArray2_2.T], 'color', 'auto', 'ax', ax2);
62
63% Plot molar fractions
64plotComposition(mixArray2_1(1), mixArray1, 'theta', 'Xi', 'mintol', 1e-3, 'y_var', mixArray2_1);
65plotComposition(mixArray2_2(1), mixArray1, 'theta', 'Xi', 'mintol', 1e-3, 'y_var', mixArray2_2);
66
67% Plot properties
68properties = {'T', 'p', 'rho', 'h', 'e', 'g', 's', 'gamma_s'};
69propertiesBasis = {[], [], [], 'mi', 'mi', 'mi', 'mi', []};
70ax = plotProperties(repmat({'theta'}, 1, length(properties)), mixArray2_1, properties, mixArray2_1, 'basis', propertiesBasis, 'color', [0 0 0], 'linestyle', '-');
71ax = plotProperties(repmat({'theta'}, 1, length(properties)), mixArray2_2, properties, mixArray2_2, 'basis', propertiesBasis, 'color', [0 0 0], 'linestyle', '--', 'ax', ax);