1% -------------------------------------------------------------------------
2% EXAMPLE: SHOCK_I_THERMO
3%
4% Influence of caloric models on jump conditions in normal shocks
5%
6% This script examines the effects of different caloric models on the jump
7% conditions (changes in temperature, density, and adiabatic index)
8% encountered in normal shock waves. The models under consideration are:
9%
10% 1. Thermochemical frozen: assumes a calorically perfect gas, where
11% specific heat values remain constant.
12%
13% 2. Frozen: assumes a thermally perfect gas, where specific heat values
14% vary with temperature.
15%
16% 3. Equilibrium: assumes a calorically imperfect gas, where both
17% thermal and caloric properties vary, accounting for chemical
18% reactions in equilibrium.
19%
20% @author: Alberto Cuadra Lara
21%
22% Last update Oct 19 2025
23% -------------------------------------------------------------------------
24
25% Import packages
26import combustiontoolbox.databases.NasaDatabase
27import combustiontoolbox.core.*
28import combustiontoolbox.shockdetonation.ShockSolver
29import combustiontoolbox.utils.display.*
30
31% Get Nasa database
32DB = NasaDatabase();
33
34% Define chemical system
35system = ChemicalSystem(DB);
36
37% Initialize mixture
38mix = Mixture(system);
39
40% Define chemical state
41set(mix, {'N2', 'O2'}, [79/21, 1]);
42
43% Initialize figure
44plotConfig = PlotConfig();
45plotConfig.innerposition = [0.05, 0.05, 0.45, 0.55];
46plotConfig.outerposition = [0.05, 0.05, 0.45, 0.55];
47ax1 = setFigure(plotConfig);
48ax2 = setFigure(plotConfig);
49ax3 = setFigure(plotConfig);
50
51% Calculations
52for i = 1:3
53
54 % Define caloric model and linestyle
55 switch i
56 case 1
57 caloricGasModel = CaloricGasModel.perfect;
58 linestyle = ':';
59 case 2
60 caloricGasModel = CaloricGasModel.thermallyPerfect;
61 linestyle = '--';
62 case 3
63 caloricGasModel = CaloricGasModel.imperfect;
64 linestyle = '-';
65 end
66
67 % Define properties
68 mixArray1 = setProperties(mix, 'temperature', 300, 'pressure', 1, 'M1', 1:0.1:10);
69
70 % Initialize solver
71 solver = ShockSolver('problemType', 'SHOCK_I', 'caloricGasModel', caloricGasModel);
72
73 % Solve problem
74 [mixArray1, mixArray2] = solver.solveArray(mixArray1);
75
76 % Plots
77 ax1 = plotFigure('M1', [mixArray1.mach], 'T_2/T_1', [mixArray2.T] ./ [mixArray1.T], 'linestyle', linestyle, 'color', [0 0 0], 'ax', ax1);
78 ax2 = plotFigure('M1', [mixArray1.mach], '\rho_2/\rho_1', [mixArray2.rho] ./ [mixArray1.rho], 'linestyle', linestyle, 'color', [0 0 0], 'ax', ax2);
79 ax3 = plotFigure('M1', [mixArray1.mach], '\gamma_s', [mixArray2.gamma_s], 'linestyle', linestyle, 'color', [0 0 0], 'ax', ax3);
80
81 if i ~= 3
82 continue
83 end
84
85 ax1 = plotFigure('M1', [5 5], 'T_2/T_1', [ax1.YLim(1), ax1.YLim(2)], 'linestyle', '--', 'color', [0.5 0.5 0.5], 'ax', ax1);
86 ax2 = plotFigure('M1', [5 5], '\rho_2/\rho_1', [ax2.YLim(1), ax2.YLim(2)], 'linestyle', '--', 'color', [0.5 0.5 0.5], 'ax', ax2);
87 ax3 = plotFigure('M1', [5 5], '\gamma_s', [ax3.YLim(1), ax3.YLim(2)], 'linestyle', '--', 'color', [0.5 0.5 0.5], 'ax', ax3);
88end