Skip to content

Add: solve SHOCK_PRANDTL_MEYER expansion within the ShockSolver class #1082

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
AlbertoCuadra merged 3 commits into from
Oct 16, 2025
Merged

Add: solve SHOCK_PRANDTL_MEYER expansion within the ShockSolver class #1082

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
662c64c
Update: miscellaneous
AlbertoCuadra Oct 16, 2025
44a3c84
Add: solve `SHOCK_PRANDTL_MEYER` expansion within the `ShockSolver` c…
AlbertoCuadra Oct 16, 2025
8a38476
Update: typo
AlbertoCuadra Oct 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions +combustiontoolbox/+common/@Constants/Constants.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,8 @@
E = 1.602176634e-19 % Electron charge [C]
E0 = 8.854187818e-12 % Vacuum permittivity [F/m]
ME = 9.10938356e-31 % Electron mass [kg]
release = 'v1.2.7' % Release version
date = '14 Oct 2025' % Release date
release = 'v1.2.8' % Release version
date = '16 Oct 2025' % Release date
end

end
4 changes: 4 additions & 0 deletions +combustiontoolbox/+core/@Mixture/print.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,10 @@ function print_properties(ProblemType, numberMixtures, mix)
fprintf(['sonic def.[deg]| |', string_value_2], get_properties('thetaSonic', numberMixtures - 1, mix(2:end)));
end

elseif contains(ProblemType, 'PRANDTL_MEYER')
fprintf('------------------------------------------------------------------------------------------------------------\n');
fprintf('PARAMETERS\n');
fprintf(['deflection[deg]| |', string_value_2], get_properties('theta', numberMixtures - 1, mix(2:end)));
elseif contains(ProblemType, 'ROCKET')
string_value_3 = set_string_value(numberMixtures - 2);

Expand Down
69 changes: 52 additions & 17 deletions +combustiontoolbox/+shockdetonation/@ShockSolver/ShockSolver.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,26 +18,28 @@
% * ``SHOCK_POLAR``: Shock polar diagrams
% * ``SHOCK_POLAR_R``: Shock polar diagrams for incident and reflected states
% * ``SHOCK_POLAR_LIMITRR``: Shock polar within the limit of regular reflection
% * ``SHOCK_PRANDTL_MEYER``: Prandtl-Meyer expansion wave
%
% See also: :mat:func:`Mixture`, :mat:func:`EquilibriumSolver`, :mat:func:`DetonationSolver`, :mat:func:`solve`, :mat:func:`solveArray`, :mat:func:`report`

properties
problemType % Problem type
equilibriumSolver % EquilibriumSolver object
tol0 = 1e-5 % Tolerance of shocks/detonations kernel
itMax = 50 % Max number of iterations - shocks and detonations
machThermo = 2 % Pre-shock Mach number above which T2_guess will be computed considering h2 = h1 + u1^2 / 2
tolOblique = 1e-3 % Tolerance oblique shocks algorithm
itOblique = 20 % Max number of iterations - oblique shocks
numPointsPolar = 100 % Number of points to compute shock/detonation polar curves
tolLimitRR = 1e-4 % Tolerance to calculate the limit of regular reflections
itLimitRR = 10 % Max number of iterations - limit of regular reflections
FLAG_RESULTS = true % Flag to print results
FLAG_TIME = true % Flag to print elapsed time
FLAG_REPORT = false % Flag to print predefined plots
FLAG_CACHE = true % Flag to clear cache after calculations
time % Elapsed time
plotConfig % PlotConfig object
problemType % Problem type
equilibriumSolver % EquilibriumSolver object
tol0 = 1e-5 % Tolerance of shocks/detonations kernel
itMax = 50 % Max number of iterations - shocks and detonations
machThermo = 2 % Pre-shock Mach number above which T2_guess will be computed considering h2 = h1 + u1^2 / 2
tolOblique = 1e-3 % Tolerance oblique shocks algorithm
itOblique = 20 % Max number of iterations - oblique shocks
numPointsPolar = 100 % Number of points to compute shock/detonation polar curves
numPointsPrandtlMeyer = 100 % Number of points to compute Prandtl-Meyer curves
tolLimitRR = 1e-4 % Tolerance to calculate the limit of regular reflections
itLimitRR = 10 % Max number of iterations - limit of regular reflections
FLAG_RESULTS = true % Flag to print results
FLAG_TIME = true % Flag to print elapsed time
FLAG_REPORT = false % Flag to print predefined plots
FLAG_CACHE = true % Flag to clear cache after calculations
time % Elapsed time
plotConfig % PlotConfig object
end

methods
Expand All @@ -49,6 +51,7 @@
[mix1, mix2, mix5_1, mix5_2] = shockObliqueReflectedTheta(obj, mix1, u2, theta, mix2, varargin)
[mix1, mix2] = shockPolar(obj, mix1, u1)
[mix1, mix2, mix2_1, mix3] = shockPolarLimitRR(obj, mix1, u1)
[mix1, mix2] = shockPrandtlMeyer(obj, mix1, u1, theta2, varargin)
[R, P, T, Gammas, M1] = shockIncidentIdeal(obj, gamma, M1)

function obj = ShockSolver(varargin)
Expand All @@ -61,14 +64,15 @@

% Parse input arguments
p = inputParser;
addOptional(p, 'problemType', defaultProblemType, @(x) ischar(x) && any(strcmpi(x, {'SHOCK_I', 'SHOCK_R', 'SHOCK_OBLIQUE', 'SHOCK_OBLIQUE_R', 'SHOCK_POLAR', 'SHOCK_POLAR_R', 'SHOCK_POLAR_LIMITRR'})));
addOptional(p, 'problemType', defaultProblemType, @(x) ischar(x) && any(strcmpi(x, {'SHOCK_I', 'SHOCK_R', 'SHOCK_OBLIQUE', 'SHOCK_OBLIQUE_R', 'SHOCK_POLAR', 'SHOCK_POLAR_R', 'SHOCK_POLAR_LIMITRR', 'SHOCK_PRANDTL_MEYER'})));
addParameter(p, 'equilibriumSolver', defaultEquilibriumSolver);
addParameter(p, 'tol0', obj.tol0, @(x) isnumeric(x) && x > 0);
addParameter(p, 'itMax', obj.itMax, @(x) isnumeric(x) && x > 0);
addParameter(p, 'machThermo', obj.machThermo, @(x) isnumeric(x) && x >= 1);
addParameter(p, 'tolOblique', obj.tolOblique, @(x) isnumeric(x) && x > 0);
addParameter(p, 'itOblique', obj.itOblique, @(x) isnumeric(x) && x > 0);
addParameter(p, 'numPointsPolar', obj.numPointsPolar, @(x) isnumeric(x) && x > 0);
addParameter(p, 'numPointsPrandtlMeyer', obj.numPointsPrandtlMeyer, @(x) isnumeric(x) && x > 0);
addParameter(p, 'tolLimitRR', obj.tolLimitRR, @(x) isnumeric(x) && x > 0);
addParameter(p, 'itLimitRR', obj.itLimitRR, @(x) isnumeric(x) && x > 0);
addParameter(p, 'FLAG_RESULTS', obj.FLAG_RESULTS, @(x) islogical(x));
Expand All @@ -89,6 +93,7 @@
obj.tolOblique = p.Results.tolOblique;
obj.itOblique = p.Results.itOblique;
obj.numPointsPolar = p.Results.numPointsPolar;
obj.numPointsPrandtlMeyer = p.Results.numPointsPrandtlMeyer;
obj.tolLimitRR = p.Results.tolLimitRR;
obj.itLimitRR = p.Results.itLimitRR;
obj.FLAG_RESULTS = p.Results.FLAG_RESULTS;
Expand Down Expand Up @@ -344,6 +349,26 @@

% Set output
varargout = {mix1, mix2, mix2_1, mix3};

case 'SHOCK_PRANDTL_MEYER'
% Solve problem
if nargin > 2
[mix1, mix2] = obj.shockPrandtlMeyer(mix1, u1, theta, varargin{1});
else
[mix1, mix2] = obj.shockPrandtlMeyer(mix1, u1, theta);
end

% Set problemType
mix1.problemType = obj.problemType;
mix2.problemType = obj.problemType;

% Print results
if obj.FLAG_RESULTS
print(mix1, mix2);
end

% Set output
varargout = {mix1, mix2};
end

end
Expand Down Expand Up @@ -394,6 +419,16 @@
% Set output
varargout = {mixArray1, mixArray2};

case {'SHOCK_PRANDTL_MEYER_THETA'}
[mixArray1(1), mixArray2(1)] = obj.solve(mixArray1(1));

for i = 2:n
[mixArray1(i), mixArray2(i)] = obj.solve(mixArray1(i), mixArray2(i - 1));
end

% Set output
varargout = {mixArray1, mixArray2};

case {'SHOCK_R', 'SHOCK_OBLIQUE_THETA'}
% Initialization
mixArray3 = mixArray1;
Expand Down
Loading