Incident shock waves for different thermodynamic models#

 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%          Postdoctoral researcher - Group Fluid Mechanics
22%          Universidad Carlos III de Madrid
23%                  
24% Last update Sep 05 2024
25% -------------------------------------------------------------------------
26
27% Import packages
28import combustiontoolbox.databases.NasaDatabase
29import combustiontoolbox.core.*
30import combustiontoolbox.shockdetonation.ShockSolver
31import combustiontoolbox.utils.display.*
32
33% Get Nasa database
34DB = NasaDatabase();
35
36% Define chemical system
37system = ChemicalSystem(DB);
38
39% Initialize mixture
40mix = Mixture(system);
41
42% Define chemical state
43set(mix, {'N2', 'O2'}, [79/21, 1]);
44
45% Initialize figure
46plotConfig = PlotConfig();
47plotConfig.innerposition = [0.05, 0.05, 0.45, 0.55];
48plotConfig.outerposition = [0.05, 0.05, 0.45, 0.55];
49ax1 = setFigure(plotConfig);
50ax2 = setFigure(plotConfig);
51ax3 = setFigure(plotConfig);
52
53% Calculations
54for i = 1:3
55
56    switch i
57        case 1
58            FLAG_TCHEM_FROZEN = true;
59            FLAG_FROZEN = false;
60            linestyle = ':';
61        case 2
62            FLAG_TCHEM_FROZEN = false;
63            FLAG_FROZEN = true;
64            linestyle = '--';
65        case 3
66            FLAG_TCHEM_FROZEN = false;
67            FLAG_FROZEN = false;
68            linestyle = '-';
69    end
70
71    % Define properties
72    mixArray1 = setProperties(mix, 'temperature', 300, 'pressure', 1, 'M1', 1:0.1:10);
73
74    % Initialize solver
75    solver = ShockSolver('problemType', 'SHOCK_I', 'FLAG_TCHEM_FROZEN', FLAG_TCHEM_FROZEN, 'FLAG_FROZEN', FLAG_FROZEN);
76
77    % Solve problem
78    [mixArray1, mixArray2] = solver.solveArray(mixArray1);
79    
80    % Plots
81    ax1 = plotFigure('M1', [mixArray1.mach], 'T_2/T_1', [mixArray2.T] ./ [mixArray1.T], 'linestyle', linestyle, 'color', [0 0 0], 'ax', ax1);
82    ax2 = plotFigure('M1', [mixArray1.mach], '\rho_2/\rho_1', [mixArray2.rho] ./ [mixArray1.rho], 'linestyle', linestyle, 'color', [0 0 0], 'ax', ax2);
83    ax3 = plotFigure('M1', [mixArray1.mach], '\gamma_s', [mixArray2.gamma_s], 'linestyle', linestyle, 'color', [0 0 0], 'ax', ax3);
84
85    if i ~= 3
86        continue
87    end
88
89    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);
90    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);
91    ax3 = plotFigure('M1', [5 5], '\gamma_s', [ax3.YLim(1), ax3.YLim(2)], 'linestyle', '--', 'color', [0.5 0.5 0.5], 'ax', ax3);
92end