1% -------------------------------------------------------------------------
2% EXAMPLE: SHOCK_R
3%
4% Compute pre-shock and post-shock state for a planar reflected shock wave
5% at standard conditions, a set of 16 species considered and a set of
6% initial shock front velocities (u1) contained in (360, 9000) [m/s]
7%
8% Air == {'O2','N2','O','O3','N','NO','NO2','NO3','N2O','N2O3','N2O4',...
9% 'N3','C','CO','CO2','Ar'}
10%
11% See wiki or setListspecies method from ChemicalSystem class for more
12% predefined sets of species
13%
14% @author: Alberto Cuadra Lara
15% Postdoctoral researcher - Group Fluid Mechanics
16% Universidad Carlos III de Madrid
17%
18% Last update April 02 2024
19% -------------------------------------------------------------------------
20
21% Import packages
22import combustiontoolbox.databases.NasaDatabase
23import combustiontoolbox.core.*
24import combustiontoolbox.shockdetonation.*
25import combustiontoolbox.utils.display.*
26
27% Get Nasa database
28DB = NasaDatabase();
29
30% Define chemical system
31system = ChemicalSystem(DB, 'air');
32
33% Initialize mixture
34mix = Mixture(system);
35
36% Define chemical state
37set(mix, {'N2', 'O2', 'Ar', 'CO2'}, [78.084, 20.9476, 0.9365, 0.0319] / 20.9476);
38
39% Define properties
40u1 = logspace(2, 5, 500); u1 = u1(u1<9000); u1 = u1(u1>=360);
41mixArray1 = setProperties(mix, 'temperature', 300, 'pressure', 1.01325, 'u1', u1);
42
43% Initialize solver
44solver = ShockSolver('problemType', 'SHOCK_R');
45
46% Solve problem
47[mixArray1, mixArray2, mixArray3] = solver.solveArray(mixArray1);
48
49% Plot Hugoniot curve
50ax1 = plotFigure('\rho_1 / \rho_2', [mixArray1.rho] ./ [mixArray2.rho], 'p_2 / p_1', [mixArray2.p] ./ [mixArray1.p], 'xScale', 'log', 'yScale', 'log', 'linestyle', '-');
51ax1 = plotFigure('\rho_1 / \rho_2', [mixArray1.rho] ./ [mixArray3.rho], 'p_2 / p_1', [mixArray3.p] ./ [mixArray1.p], 'xScale', 'log', 'yScale', 'log', 'linestyle', '--', 'ax', ax1);
52legend(ax1, {'Incident', 'Reflected'}, 'Interpreter', 'latex', 'FontSize', ax1.FontSize);
53
54% Plot post-shock temperature
55ax2 = plotFigure('u1', [mixArray1.u], 'T', [mixArray2.T], 'color', [0 0 0], 'linestyle', '-');
56ax2 = plotFigure('u1', [mixArray1.u], 'T', [mixArray3.T], 'color', [0 0 0], 'linestyle', '--', 'ax', ax2);
57legend(ax2, {'Incident', 'Reflected'}, 'Interpreter', 'latex', 'FontSize', ax2.FontSize);
58
59% Plot molar fractions
60plotComposition(mixArray3(1), mixArray1, 'u', 'Xi', 'mintol', 1e-3, 'y_var', mixArray3);