1% -------------------------------------------------------------------------
2% EXAMPLE: DET_OVERDRIVEN_AND_UNDERDRIVEN
3%
4% Compute pre-shock and post-shock state for a planar underdriven to
5% overdriven detonation considering Chapman-Jouguet (CJ) theory for a
6% stoichiometric CH4-air mixture at standard conditions, a set of 24
7% species considered and a set of overdrives contained in (1,10) [-].
8%
9% Soot formation == {'CO2','CO','H2O','H2','O2','N2','Ar','Cbgrb',...
10% 'C2','C2H4','CH','CH3','CH4','CN','H',...
11% 'HCN','HCO','N','NH','NH2','NH3','NO','O','OH'}
12%
13% See wiki or setListspecies method from ChemicalSystem class for more
14% predefined sets of species
15%
16% @author: Alberto Cuadra Lara
17% Postdoctoral researcher - Group Fluid Mechanics
18% Universidad Carlos III de Madrid
19%
20% Last update Oct 12 2022
21% -------------------------------------------------------------------------
22
23% Import packages
24import combustiontoolbox.databases.NasaDatabase
25import combustiontoolbox.core.*
26import combustiontoolbox.shockdetonation.*
27import combustiontoolbox.utils.display.*
28
29% Get Nasa database
30DB = NasaDatabase();
31
32% Define chemical system
33system = ChemicalSystem(DB, 'soot formation');
34
35% Initialize mixture
36mix = Mixture(system);
37
38% Define chemical state
39set(mix, {'CH4'}, 'fuel', 1);
40set(mix, {'N2', 'O2', 'Ar', 'CO2'}, 'oxidizer', [78.084, 20.9476, 0.9365, 0.0319] / 20.9476);
41
42% Define properties
43mixArray1_over = setProperties(mix, 'temperature', 300, 'pressure', 1.01325, 'equivalenceRatio', 1, 'driveFactor', 1:0.1:10);
44mixArray1_under = setProperties(mix, 'temperature', 300, 'pressure', 1.01325, 'equivalenceRatio', 1, 'driveFactor', 1:0.1:10);
45
46% Initialize solver
47solver1 = DetonationSolver('problemType', 'DET_OVERDRIVEN');
48solver2 = DetonationSolver('problemType', 'DET_UNDERDRIVEN');
49
50% Solve problem
51[mixArray1_over, mixArray2_over] = solver1.solveArray(mixArray1_over);
52[mixArray1_under, mixArray2_under] = solver2.solveArray(mixArray1_under);
53
54% Plot Hugoniot curve
55ax1 = plotFigure('\rho_1 / \rho_2', [mixArray1_over.rho] ./ [mixArray2_over.rho], 'p_2 / p_1', [mixArray2_over.p] ./ [mixArray1_over.p], 'xScale', 'log', 'yScale', 'log', 'linestyle', '-');
56ax1 = plotFigure('\rho_1 / \rho_2', [mixArray1_under.rho] ./ [mixArray2_under.rho], 'p_2 / p_1', [mixArray2_under.p] ./ [mixArray1_under.p], 'xScale', 'log', 'yScale', 'log', 'linestyle', '--', 'ax', ax1);
57legend(ax1, {'Overdriven', 'Underdriven'}, 'Interpreter', 'latex', 'FontSize', ax1.FontSize);
58
59% Plot post-shock temperature
60ax2 = plotFigure('driveFactor', mixArray1_over, 'T', mixArray2_over, 'linestyle', '-');
61ax2 = plotFigure('driveFactor', mixArray1_under, 'T', mixArray2_under, 'linestyle', '--', 'ax', ax2);
62legend(ax2, {'Overdriven', 'Underdriven'}, 'Interpreter', 'latex', 'FontSize', ax1.FontSize);
63
64% Plot molar fractions
65plotComposition(mixArray2_over(1), mixArray1_over, 'driveFactor', 'Xi', 'mintol', 1e-14, 'y_var', mixArray2_over);
66plotComposition(mixArray2_under(1), mixArray1_under, 'driveFactor', 'Xi', 'mintol', 1e-14, 'y_var', mixArray2_under);