Skip to content

Commit

Permalink
Merge pull request #463 from AlbertoCuadra/develop
Browse files Browse the repository at this point in the history
Add: include Finite-Area-Chamber model (rockets)
  • Loading branch information
AlbertoCuadra committed Apr 12, 2022
2 parents 437632f + 6f1115f commit ad2401d
Show file tree
Hide file tree
Showing 15 changed files with 215 additions and 102 deletions.
33 changes: 33 additions & 0 deletions Examples/Example_ROCKET_FAC.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
% -------------------------------------------------------------------------
% EXAMPLE: ROCKET Propellants considering an Infinite-Area-Chamber (IAC)
%
% Compute adiabatic temperature and equilibrium composition at constant
% pressure (e.g., 1.01325 bar) for lean to rich LH2-LOX mixtures at
% standard conditions, a set of 24 species considered and a set of
% equivalence ratios phi contained in (2, 5) [-]
%
% HYDROGEN_L == {'H','H2O','OH','H2','O','O3','O2','HO2','H2O2',...
% 'H2bLb','O2bLb'}
%
% See wiki or ListSpecies() for more predefined sets of species
%
% @author: Alberto Cuadra Lara
% PhD Candidate - Group Fluid Mechanics
% Universidad Carlos III de Madrid
%
% Last update April 12 2022
% -------------------------------------------------------------------------

%% INITIALIZE
self = App('HYDROGEN_L');
%% INITIAL CONDITIONS
self = set_prop(self, 'TR', 90, 'pR', 100 * 1.01325, 'phi', 2);
self.PD.S_Fuel = {'H2bLb'};
self.PD.S_Oxidizer = {'O2bLb'};
self.PD.FLAG_IAC = false;
%% ADDITIONAL INPUTS (DEPENDS OF THE PROBLEM SELECTED)
self = set_prop(self, 'Aratio_c', 2);
%% SOLVE PROBLEM
self = SolveProblem(self, 'ROCKET');
%% DISPLAY RESULTS (PLOTS)
postResults(self);
9 changes: 5 additions & 4 deletions Examples/Example_ROCKET.m → Examples/Example_ROCKET_IAC.m
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
% -------------------------------------------------------------------------
% EXAMPLE: ROCKET Propellants
% EXAMPLE: ROCKET Propellants considering an Infinite-Area-Chamber (IAC)
%
% Compute adiabatic temperature and equilibrium composition at constant
% pressure (e.g., 1.01325 bar) for lean to rich LH2-LOX mixtures at
% standard conditions, a set of 24 species considered and a set of
% equivalence ratios phi contained in (2, 8) [-]
% equivalence ratios phi contained in (2, 5) [-]
%
% HYDROGEN_L == {'H','H2O','OH','H2','O','O3','O2','HO2','H2O2',...
% 'H2bLb','O2bLb'}
Expand All @@ -15,15 +15,16 @@
% PhD Candidate - Group Fluid Mechanics
% Universidad Carlos III de Madrid
%
% Last update March 24 2022
% Last update April 12 2022
% -------------------------------------------------------------------------

%% INITIALIZE
self = App('HYDROGEN_L');
%% INITIAL CONDITIONS
self = set_prop(self, 'TR', 90, 'pR', 1 * 1.01325, 'phi', 2:0.05:8);
self = set_prop(self, 'TR', 90, 'pR', 100 * 1.01325, 'phi', 1:0.05:5);
self.PD.S_Fuel = {'H2bLb'};
self.PD.S_Oxidizer = {'O2bLb'};
self.PD.FLAG_IAC = true;
%% SOLVE PROBLEM
self = SolveProblem(self, 'ROCKET');
%% DISPLAY RESULTS (PLOTS)
Expand Down
24 changes: 13 additions & 11 deletions Settings/display/displayresults.m
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ function displayresults(varargin)
fprintf('max def. [deg]| | %12.4f | %12.4f\n', mix2.theta_max, mix3.theta_max);
fprintf('sonic def.[deg]| | %12.4f | %12.4f\n', mix2.theta_sonic, mix3.theta_sonic);
end
elseif strcmpi(ProblemType, 'ROCKET')
elseif contains(ProblemType, 'ROCKET')
fprintf('------------------------------------------------------------------------\n');
fprintf('PERFORMANCE PARAMETERS\n');
fprintf('CSTAR [m/s] | | | %12.4f \n', mix3.cstar);
Expand Down Expand Up @@ -227,7 +227,7 @@ function displayresults(varargin)
elseif contains(ProblemType, '_R')
fprintf(' | STATE 1 | STATE 2 | STATE 3 | STATE 4\n');
else
fprintf(' | INLET CHAMBER | OUTLET CHAMBER | THROAT | EXIT\n');
fprintf(' | INLET CHAMBER | INJECTOR | OUTLET CHAMBER | THROAT \n');
end

fprintf('T [K] | %12.4f | %12.4f | %12.4f | %12.4f\n', temperature(mix1), temperature(mix2), temperature(mix3), temperature(mix4));
Expand Down Expand Up @@ -255,13 +255,14 @@ function displayresults(varargin)
fprintf('max def. [deg]| | %12.4f | %12.4f | %12.4f\n', mix2.theta_max, mix3.theta_max, mix4.theta_max);
fprintf('sonic def.[deg]| | %12.4f | %12.4f | %12.4f\n', mix2.theta_sonic, mix3.theta_sonic, mix4.theta_sonic);
end
elseif strcmpi(ProblemType, 'ROCKET')
elseif contains(ProblemType, 'ROCKET')
fprintf('--------------------------------------------------------------------------------------\n');
fprintf('PERFORMANCE PARAMETERS\n');
fprintf('CSTAR [m/s] | | | %12.4f | %12.4f |\n', mix3.cstar, mix4.cstar);
fprintf('CF [-] | | | %12.4f | %12.4f |\n', mix3.cf, mix4.cf);
fprintf('Ivac [s] | | | %12.4f | %12.4f |\n', mix3.I_vac, mix4.I_vac);
fprintf('Isp [s] | | | %12.4f | %12.4f |\n', mix3.I_sp, mix4.I_sp);
fprintf('PERFORMANCE PARAMETERS\n');
fprintf('A/At [-] | | | %12.4f | %12.4f \n', mix3.Aratio, mix4.Aratio);
fprintf('CSTAR [m/s] | | | %12.4f | %12.4f \n', mix3.cstar, mix4.cstar);
fprintf('CF [-] | | | %12.4f | %12.4f \n', mix3.cf, mix4.cf);
fprintf('Ivac [s] | | | %12.4f | %12.4f \n', mix3.I_vac, mix4.I_vac);
fprintf('Isp [s] | | | %12.4f | %12.4f \n', mix3.I_sp, mix4.I_sp);
end
fprintf('--------------------------------------------------------------------------------------\n');

Expand All @@ -288,7 +289,7 @@ function displayresults(varargin)
if contains(ProblemType, '_R')
fprintf('STATE 2 Xi [-]\n');
else
fprintf('OUTLET CHAMBER Xi [-]\n');
fprintf('INJECTOR Xi [-]\n');
end
%%%% SORT SPECIES COMPOSITION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[mix2.Xi(:), ind_sort] = sort(mix2.Xi(:),'descend');
Expand All @@ -311,7 +312,7 @@ function displayresults(varargin)
elseif contains(ProblemType, '_R')
fprintf('STATE 3 Xi [-]\n');
else
fprintf('THROAT Xi [-]\n');
fprintf('OUTLET CHAMBER Xi [-]\n');
end
%%%% SORT SPECIES COMPOSITION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[mix3.Xi(:), ind_sort] = sort(mix3.Xi(:), 'descend');
Expand All @@ -334,7 +335,8 @@ function displayresults(varargin)
elseif contains(ProblemType, '_R')
fprintf('STATE 4 Xi [-]\n');
else
fprintf('EXIT Xi [-]\n');
fprintf('THROAT Xi [-]\n');
% fprintf('EXIT Xi [-]\n');
end
%%%% SORT SPECIES COMPOSITION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[mix4.Xi(:), ind_sort] = sort(mix4.Xi(:),'descend');
Expand Down
2 changes: 1 addition & 1 deletion Settings/display/displaysweepresults.m
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
mintol = self.C.mintol_display;

% Function
if length(xvar)>1
if length(xvar) > 1
ax = set_figure(config);
[yvar, all_ind, indlabel] = get_parameters(mix, species, display_species, mintol);
plot_line(self, ax, xvar, yvar, all_ind, indlabel, species)
Expand Down
5 changes: 5 additions & 0 deletions Settings/display/results.m
Original file line number Diff line number Diff line change
@@ -1,8 +1,13 @@
function results(self, i)
% Display results in the command window

if self.Misc.FLAG_RESULTS
if ~isfield(self.PS, 'str2')
displayresults(self.PS.strR{i}, self.PS.strP{i}, self.PD.ProblemType, self.C.mintol_display, self.S.LS);
elseif isfield(self.PS, 'mix2_c') && isempty(self.PS.strP{i}) && self.PD.FLAG_IAC
displayresults(self.PS.strR{i}, self.PS.mix2_c{i}, self.PS.mix3{i}, self.PD.ProblemType, self.C.mintol_display, self.S.LS);
elseif isfield(self.PS, 'mix2_c') && isempty(self.PS.strP{i}) && ~self.PD.FLAG_IAC
displayresults(self.PS.strR{i}, self.PS.str2{i}, self.PS.mix2_c{i}, self.PS.mix3{i}, self.PD.ProblemType, self.C.mintol_display, self.S.LS);
elseif ~isfield(self.PS, 'str3_2')
displayresults(self.PS.strR{i}, self.PS.str2{i}, self.PS.strP{i}, self.PD.ProblemType, self.C.mintol_display, self.S.LS);
elseif ~isfield(self.PS, 'str3')
Expand Down
13 changes: 8 additions & 5 deletions Settings/functions/postResults.m
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,9 @@ function postResults(self)
self.Misc.FLAGS_PROP.TP = false;
end

if isempty(mix2{end})
mix2 =self.PS.mix3;
end
% PLOTS REACTANTS
if self.Misc.FLAGS_PROP.TR && length(phi) > 1
self.Misc.config.labelx = 'Temperature reactants $T$ [K]';
Expand Down Expand Up @@ -186,7 +189,7 @@ function postResults(self)
if strcmp(ProblemType,{'ROCKET'}) && length(phi) > 1
self.Misc.config.labelx = 'Equivalence Ratio $\phi$';
self.Misc.config.labely = 'Characteristic velocity $c^*$ [m/s]';
ax = plot_figure(self.PS.strR, self.PS.strP, 'phi', 'cstar', self.Misc.config, self.PD.CompleteOrIncomplete);
ax = plot_figure(mix1, mix2, 'phi', 'cstar', self.Misc.config, self.PD.CompleteOrIncomplete);
% legend_name = sprintf('$p = %.2f$ [bar]', self.PS.strR.p);
% for i = 2:length(self)
% legend_name{i} = sprintf('$p = %.2f$ [bar]', self{i}.PS.strR.p);
Expand All @@ -203,14 +206,14 @@ function postResults(self)
legend_name = {'$I_{sp}$', '$I_{vac}$'};
self.Misc.config.labelx = 'Equivalence Ratio $\phi$';
self.Misc.config.labely = 'Specific impulse $I$ [s]';
ax = plot_figure(self.PS.strR, self.PS.strP, 'phi', 'I_sp', self.Misc.config, self.PD.CompleteOrIncomplete);
ax = plot_figure(self.PS.strR, self.PS.strP, 'phi', 'I_vac', self.Misc.config, self.PD.CompleteOrIncomplete, ax);
ax = plot_figure(mix1, mix2, 'phi', 'I_sp', self.Misc.config, self.PD.CompleteOrIncomplete);
ax = plot_figure(mix1, mix2, 'phi', 'I_vac', self.Misc.config, self.PD.CompleteOrIncomplete, ax);
set_legends(ax, legend_name, self.Misc.config)

self.Misc.config.labelx = '$Mixture ratio O/F$';
self.Misc.config.labely = 'Specific impulse $I$ [s]';
ax = plot_figure(self.PS.strR, self.PS.strP, 'FO', 'I_sp', self.Misc.config, self.PD.CompleteOrIncomplete);
ax = plot_figure(self.PS.strR, self.PS.strP, 'FO', 'I_vac', self.Misc.config, self.PD.CompleteOrIncomplete, ax);
ax = plot_figure(mix1, mix2, 'FO', 'I_sp', self.Misc.config, self.PD.CompleteOrIncomplete);
ax = plot_figure(mix1, mix2, 'FO', 'I_vac', self.Misc.config, self.PD.CompleteOrIncomplete, ax);
set_legends(ax, legend_name, self.Misc.config)
% for i = 2:length(self)
% ax = plot_figure(self{i}.PS.strR, self{i}.PS.strP, 'phi', 'I_sp', self.Misc.config, self.PD.CompleteOrIncomplete, ax, legend_name);
Expand Down
9 changes: 6 additions & 3 deletions Settings/self/ProblemDescription/ProblemDescription.m
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,10 @@
self.theta.value = []; % [deg]
self.beta.description = "Wave angle - oblique shocks";
self.beta.value = []; % [deg]
self.Aratio.description = "Area ratio combustion chamber - rocket";
self.Aratio.description = "Area ratio exit/throat - rocket";
self.Aratio.value = []; % [-]
self.Aratio_c.description = "Area ratio combustion chamber/thoat - rocket";
self.Aratio_c.value = []; % [-]
% * Mixture conditions
self.S_Fuel = []; % Cell with the list of fuel species in the mixture
self.N_Fuel = []; % Vector with the number of moles of the fuel species in the mixture
Expand All @@ -55,6 +57,7 @@
self.T_Inert = []; % Vector with the temperature values of the inert species in the mixture
self.proportion_inerts_O2 = []; % Proportion Inerts / O2 [-]
% * Flags
self.FLAG_ION = false; % Flag ionized species in the system
self.FLAG_IAC = true; % Flag use IAC model for rocket computations
self.FLAG_ION = false; % Flag ionized species in the system
self.FLAG_IAC = true; % Flag use IAC model for rocket computations
self.FLAG_SUBSONIC = false; % Flag subsonic Area ratio
end
46 changes: 0 additions & 46 deletions Solver/Rocket/compute_IAC_model.m

This file was deleted.

80 changes: 80 additions & 0 deletions Solver/Rocket/compute_chamber_FAC.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
function [mix2_inj, mix2_c, mix3] = compute_chamber_FAC(self, mix1, varargin)
% Compute chemical equilibria at the injector, outlet of the chamber,
% using the Finite-Area-Chamber (FAC) model
%
% This method is based on Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.
%
% Args:
% self (struct): Data of the mixture, conditions, and databases
% mix1 (struct): Properties of the initial mixture
% mix2_inj (struct): Properties of the mixture at the injector of the chamber (previous calculation)
% mix2_c (struct): Properties of the mixture at the outlet of the chamber (previous calculation)
%
% Returns:
% mix2_inj (struct): Properties of the mixture at the injector of the chamber
% mix2_c (struct): Properties of the mixture at the outlet of the chamber
% mix3 (struct): Properties of the mixture at the throat

% Abbreviations
TN = self.TN;
% Unpack
if nargin > 3, mix2_inj = get_struct(varargin{1}); else, mix2_inj = []; end
if nargin > 4, mix2_c = get_struct(varargin{2}); else, mix2_c = []; end
if nargin > 5, mix3 = get_struct(varargin{3}); else, mix3 = []; end
% Definitions
Aratio_chamber = self.PD.Aratio_c.value;
% Obtain mixture compostion and properties at the injector of the chamber (equivalent to the outlet of the chamber using IAC model)
mix2_inj = compute_chamber_IAC(self, mix1, mix2_inj);
% Get results
pressure_inj = pressure(mix2_inj); % [bar]
% Compute guess pressure_inf
pressure_inf = compute_initial_guess_pressure(pressure_inj, Aratio_chamber); % [bar]
% Initialization
STOP = 1; it = 0;
mix2_inf = [];
pressure_inj = pressure_inj * 1e5; % [Pa]
% Loop
while STOP > TN.tol_rocket && it < TN.it_rocket
it = it + 1;
% Get guess
mix1.p = pressure_inf; % [bar]
% Obtain mixture compostions
[mix2_inf, mix3, mix2_c] = compute_IAC_model(self, mix1, mix2_inf, mix3, mix2_c, Aratio_chamber);
% Get results
pressure_c = pressure(mix2_c) * 1e5; % [Pa]
density_c = density(mix2_c); % [kg/m3]
velocity_c = velocity_relative(mix2_c); % [m/s]
% Compute pressuse_inj_a
pressure_inj_a = (pressure_c + density_c * velocity_c^2); % [Pa]
% Check convergence
STOP = abs(pressure_inj - pressure_inj_a) / pressure_inj;
% Update guess
pressure_inf = pressure_inf * pressure_inj / pressure_inj_a; % [bar]
mix2_inf.p = pressure_inf;
end
% Assign values
mix2_c.Aratio = Aratio_chamber; % [-]
end

% SUB-PASS FUNCTIONS
function str = get_struct(var)
try
str = var{1,1};
catch
str = var;
end
end

function pressure_inf = compute_initial_guess_pressure(pressure_inj, Aratio_chamber)
% Compute initial guess of the pressure assuming an Infinite-Area-Chamber
% (IAC) for the Finite-Area-Chamber model (FAC)
pressure_inf = pressure_inj * (1.0257 - 1.2318 * Aratio_chamber) / (1 - 1.26505 * Aratio_chamber);
end

function [mix2_inf, mix3, mix2_c] = compute_IAC_model(self, mix1, mix2_inf, mix3, mix2_c, Aratio)
% Solve IAC model
self.PD.FLAG_IAC = true;
self.PD.FLAG_SUBSONIC = true;
mix2_inj = [];
[~, mix2_inf, mix3, mix2_c] = solve_model_rocket(self, mix1, mix2_inj, mix2_inf, mix3, mix2_c, Aratio);
end
6 changes: 3 additions & 3 deletions Solver/Rocket/compute_chamber_IAC.m
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,11 @@
%
% Args:
% self (struct): Data of the mixture, conditions, and databases
% mix2 (struct): Properties of the mixture at the outlet of the chamber
% mix3 (struct): Properties of the mixture at the throat (previous calculation)
% mix1 (struct): Properties of the initial mixture
% mix2 (struct): Properties of the mixture at the outlet of the chamber (previous calculation)
%
% Returns:
% mix3 (struct): Properties of the mixture at the throat
% mix2 (struct): Properties of the mixture at the outlet of the chamber

% Definitions
self.PD.ProblemType = 'HP';
Expand Down
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
function mix4 = compute_exit_IAC(self, mix2, mix3, mix4, Aratio)
% Compute thermochemical composition for the Infinite-Area-Chamber (IAC) model
function mix4 = compute_exit(self, mix2, mix3, mix4, Aratio)
% Compute thermochemical composition for a given Aratio
%
% This method is based on Gordon, S., & McBride, B. J. (1994). NASA reference publication,
% 1311.
Expand All @@ -20,7 +20,7 @@
% Definitions
self.PD.ProblemType = 'SP';
% Compute pressure guess [bar] for Infinite-Area-Chamber (IAC)
logP = guess_pressure_exit_IAC(mix2, mix3, Aratio, false);
logP = guess_pressure_exit_IAC(mix2, mix3, Aratio, self.PD.FLAG_SUBSONIC);
% Initialization
STOP = 1; it = 0; logP0 = logP;
% Loop
Expand All @@ -42,6 +42,7 @@
% Assign values
mix4.p = extract_pressure(logP, mix2.p); % [bar]
mix4.v_shock = mix4.u; % [m/s]
mix4.Aratio = Aratio; % [-]
end

% SUB-PASS FUNCTIONS
Expand Down
1 change: 1 addition & 0 deletions Solver/Rocket/compute_throat_IAC.m
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
% Assign values
mix3.p = pressure; % [bar]
mix3.v_shock = mix3.u; % [m/s]
mix3.Aratio = 1; % [-]
end

function pressure = compute_pressure(mix3)
Expand Down
Loading

0 comments on commit ad2401d

Please sign in to comment.