1% -------------------------------------------------------------------------
2% EXAMPLE: ROCKET Propellants considering an Finite-Area-Chamber (FAC)
3%
4% Compute rocket propellant performance considering an Finite-Area-Chamber
5% with an area ratio of the combustion chamber A_chamber/A_throat = 2
6% for lean to rich LH2-LOX mixtures at 101.325 bar, a set of 24 species
7% considered, a set of equivalence ratios phi contained in (2, 5) [-], and
8% an exit area ratio A_exit/A_throat = 3
9%
10% HYDROGEN_L == {'H','H2O','OH','H2','O','O3','O2','HO2','H2O2',...
11% 'H2bLb','O2bLb'}
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 Jul 30 2024
21% -------------------------------------------------------------------------
22
23% Import packages
24import combustiontoolbox.databases.NasaDatabase
25import combustiontoolbox.core.*
26import combustiontoolbox.rocket.*
27import combustiontoolbox.utils.display.*
28
29% Get Nasa database
30DB = NasaDatabase();
31
32% Define chemical system
33system = ChemicalSystem(DB, 'hydrogen_L');
34
35% Initialize mixture
36mix = Mixture(system);
37
38% Define chemical state
39set(mix, {'H2bLb'}, 'fuel', 1);
40set(mix, {'O2bLb'}, 'oxidizer', 1);
41
42% Define properties
43mixArray1 = setProperties(mix, 'temperature', 90, 'pressure', 100 * 1.01325, 'equivalenceRatio', 1:0.05:5, 'areaRatioChamber', 2, 'areaRatio', 3);
44
45% Initialize solver
46solver = RocketSolver('problemType', 'ROCKET_FAC');
47
48% Solve problem
49[mixArray1, mixArray2_inj, mixArray2_c, mixArray3, mixArray4] = solver.solveArray(mixArray1);
50
51% Plot properties
52ax1 = plotProperties(repmat({'equivalenceRatio'}, 1, 12), mixArray4, {'T', 'p', 'h', 'e', 'g', 'cp', 's', 'gamma_s', 'sound', 'u', 'I_sp', 'I_vac'}, mixArray4, 'basis', {[], [], 'mi', 'mi', 'mi', 'mi', 'mi', [], [], [], [], []});
53ax1 = plotProperties(repmat({'equivalenceRatio'}, 1, 12), mixArray3, {'T', 'p', 'h', 'e', 'g', 'cp', 's', 'gamma_s', 'sound', 'u', 'I_sp', 'I_vac'}, mixArray3, 'basis', {[], [], 'mi', 'mi', 'mi', 'mi', 'mi', [], [], [], [], []}, 'ax', ax1);
54ax1 = plotProperties(repmat({'equivalenceRatio'}, 1, 12), mixArray2_c, {'T', 'p', 'h', 'e', 'g', 'cp', 's', 'gamma_s', 'sound', 'u', 'I_sp', 'I_vac'}, mixArray2_c, 'basis', {[], [], 'mi', 'mi', 'mi', 'mi', 'mi', [], [], [], [], []}, 'ax', ax1);
55ax1 = plotProperties(repmat({'equivalenceRatio'}, 1, 10), mixArray2_inj, {'T', 'p', 'h', 'e', 'g', 'cp', 's', 'gamma_s', 'sound'}, mixArray2_inj, 'basis', {[], [], 'mi', 'mi', 'mi', 'mi', 'mi', [], [], []}, 'ax', ax1);
56leg = legend(ax1.Children(end), {'Exit', 'Throat', 'Chamber', 'Injector'}, 'Interpreter', 'latex', 'FontSize', ax1.Children(end).FontSize);
57
58% Plot molar fractions
59plotComposition(mixArray2_inj(1), mixArray1, 'equivalenceRatio', 'Xi', 'mintol', 1e-14, 'y_var', mixArray2_inj, 'title', 'Injector');
60plotComposition(mixArray2_c(1), mixArray1, 'equivalenceRatio', 'Xi', 'mintol', 1e-14, 'y_var', mixArray2_c, 'title', 'Chamber');
61plotComposition(mixArray3(1), mixArray1, 'equivalenceRatio', 'Xi', 'mintol', 1e-14, 'y_var', mixArray3, 'title', 'Throat');
62plotComposition(mixArray4(1), mixArray1, 'equivalenceRatio', 'Xi', 'mintol', 1e-14, 'y_var', mixArray4, 'title', 'Exit');