Skip to content

Commit

Permalink
Merge pull request #293 from AlbertoCuadra/develop
Browse files Browse the repository at this point in the history 
Update: separate the routine dedicated to plotting the solutions of t…
  • Loading branch information
@AlbertoCuadra
AlbertoCuadra committed Mar 23, 2022
2 parents dbb38b8 + 28e94c2 commit a515dd8
Show file tree
Hide file tree
Showing 3 changed files with 109 additions and 64 deletions.
70 changes: 70 additions & 0 deletions Settings/functions/Display/plot_shock_polar.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
function [ax1, ax2, ax3] = plot_shock_polar(self, mix1, mix2)
% Routine to obtain shock polar plots
% * Plot (pressure, deflection)
% * Plot (wave angle, deflection)
% * Plot velocity components

% Abbreviations
config = self.Misc.config;
% Definitions
N = length(mix1);
% Plots
for i = N:-1:1
% Plot (pressure, deflection) - shock polar
ax1 = plot_shock_polar_pressure(mix1{i}, mix2{i}, config);
% Plot (wave angle, deflection) - shock polar
ax2 = plot_shock_polar_wave(mix1{i}, mix2{i}, config);
% Plot velocity components
ax3 = plot_shock_polar_velocities(mix1{i}, mix2{i}, config);
end
end

% SUB-PASS FUNCTIONS
function ax = set_fixed_figure(fixed_number, config)
% Generate a figure with a defined object identifier (fixed_number)
figure(fixed_number); ax = gca;
set(ax, 'LineWidth', config.linewidth, 'FontSize', config.fontsize-2, 'BoxStyle', 'full');
grid(ax, 'off'); box(ax, 'off'); hold(ax, 'on'); ax.Layer = 'Top'; axis(ax, 'tight');
end

function ax = plot_shock_polar_pressure(mix1, mix2, config)
% Plot (pressure, deflection) - shock polar
ax = set_fixed_figure(1, config);

M1 = velocity_relative(mix1) / soundspeed(mix1);
plot(ax, [-flip(mix2.polar.theta), mix2.polar.theta], [flip(mix2.polar.p), mix2.polar.p] / pressure(mix1), 'k', 'LineWidth', config.linewidth);
% plot(ax, theta_max, mix2.polar.p(mix2.ind_max), 'ko', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
% plot(ax, mix2.polar.theta(mix2.ind_sonic), mix2.polar.p(mix2.ind_sonic), 'ks', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
text(ax, 0, max(mix2.polar.p), strcat('$\mathcal{M}_1 = ', sprintf('%.2g', M1), '$'), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', 'FontSize', config.fontsize-4, 'Interpreter', 'latex');
% text(ax, theta_max, str2.polar.p2(str2.ind_max), 'detachment point \rightarrow m', 'HorizontalAlignment', 'right', 'FontSize', config.fontsize-6)
% text(ax, theta_sonic, str2.polar.p2(str2.ind_sonic), 'sonic point \rightarrow s', 'HorizontalAlignment', 'right', 'FontSize', config.fontsize-6)
xlabel(ax, 'Deflection angle $[^\circ]$', 'FontSize', config.fontsize, 'Interpreter', 'latex');
ylabel(ax, '$p_2/p_1$','FontSize', config.fontsize, 'Interpreter', 'latex');
% ylim(ax, [1, max(str2.polar.p2)]);
end

function ax = plot_shock_polar_wave(mix1, mix2, config)
% Plot (wave angle, deflection) - shock polar
ax = set_fixed_figure(2, config);

M1 = velocity_relative(mix1) / soundspeed(mix1);
plot(ax, mix2.polar.theta, mix2.polar.beta, 'k', 'LineWidth', config.linewidth);
% plot(ax, mix2.theta_max, mix2.polar.beta(str2.ind_max), 'ko', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
% plot(ax, str2.polar.theta(str2.ind_sonic), mix2.polar.beta(str2.ind_sonic), 'ks', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
% text(ax, mix2.theta_max, mix2.polar.beta(str2.ind_max),
% strcat('$\mathcal{M}_1 = ', sprintf('%.2g', M1), '$'), 'VerticalAlignment', 'middle', 'HorizontalAlignment', 'left', 'FontSize', config.fontsize-6, 'Interpreter', 'latex');
xlabel(ax, 'Deflection angle $[^\circ]$','FontSize', config.fontsize, 'Interpreter', 'latex');
ylabel(ax, 'Wave angle $[^\circ]$','FontSize', config.fontsize, 'Interpreter', 'latex');
% xlim(ax, [0, mix2.theta_max]);
ylim(ax, [0, 90]);
end

function ax = plot_shock_polar_velocities(mix1, mix2, config)
% Plot velocity components
ax = set_fixed_figure(3, config);

plot(ax, mix2.polar.ux, mix2.polar.uy, 'k', 'LineWidth', config.linewidth);
xlabel(ax, '$u_{2x}$ [m/s]', 'FontSize', config.fontsize, 'Interpreter', 'latex');
ylabel(ax, '$u_{2y}$ [m/s]', 'FontSize', config.fontsize, 'Interpreter', 'latex');
% ylim(ax, [0, max(mix2.polar.uy)]);
end
11 changes: 10 additions & 1 deletion Settings/functions/postResults.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,16 @@ function postResults(self)
plot_figure(self.PD.vP_vR.value, mix2,'v_P/v_R','T',self.Misc.config,self.PD.CompleteOrIncomplete);
T = cell2vector(mix2, 'T');
displaysweepresults(self, mix2, T);

elseif strcmp(ProblemType,{'SHOCK_POLAR'})
% Shock polars
plot_shock_polar(self, mix1, mix2)
% % Incident velocity [m/s] againts molar fractions [-]
% self.Misc.config.labelx = 'Incident velocity $u_1$ [m/s]';
% self.Misc.config.labely = 'Molar fraction $X_i$';
% for i = length(mix2):-1:1
% beta = cell2vector(mix2{i}.polar, 'beta');
% displaysweepresults(self, mix2{i}.polar, beta);
% end
elseif strcmp(ProblemType,{'SHOCK_I'}) && length(phi) > 1
self.Misc.config.labelx = 'Incident velocity $u_1$ [m/s]';
self.Misc.config.labely = 'Temperature $T$ [K]';
Expand Down
92 changes: 29 additions & 63 deletions Solver/Shocks and detonations/shock_polar.m
View file
Original file line number Diff line number Diff line change
@@ -1,89 +1,55 @@
function [str1, str2] = shock_polar(varargin)
function [mix1, mix2] = shock_polar(varargin)
% Solve oblique shock

% Unpack input data
[self, str1, str2] = unpack(varargin);
% Abbreviations
config = self.Misc.config;
[self, mix1, mix2] = unpack(varargin);
% Definitions
a1 = soundspeed(str1);
u1 = str1.u;
M1 = u1/a1;
a1 = soundspeed(mix1);
u1 = mix1.u;
step = 5; % [m/s]
u1n = linspace(a1 + 1, u1, (u1 - a1) / step);
N = length(u1n);
beta = asin(u1n ./ u1);
beta_min = asin(a1 / u1);
un = u1 .* cos(beta);

% Solve first case for initialization
[~, str2] = shock_incident(self, str1, u1n(end), str2);
% Loop
for i = N:-1:1
[~, str2] = shock_incident(self, str1, u1n(i), str2);
a2(i) = soundspeed(str2);
u2n(i) = str2.v_shock;
p2(i) = pressure(str2);
[~, mix2] = shock_incident(self, mix1, u1n(i), mix2);
a2(i) = soundspeed(mix2);
u2n(i) = mix2.v_shock;
p2(i) = pressure(mix2);
theta(i) = beta(i) - atan(u2n(i) / un(i));
u2(i) = u2n(i) * csc(beta(i) - theta(i));
Xi(:, i) = mix2.Xi;
end

% Velocity components downstream - laboratory fixed
u2x = u2 .* cos(theta);
u2y = u2 .* sin(theta);
u2x = u2 .* cos(theta); % [m/s]
u2y = u2 .* sin(theta); % [m/s]
% Sonic point
M2 = u2 ./ a2;
ind_sonic = find(M2 > 1, 1, 'last');
theta_sonic = theta(ind_sonic) * 180/pi;
% Max deflection angle
[theta_max, ind_max] = max(theta * 180/pi); % [deg]
% Save results
str2.beta = beta(end) * 180/pi; % [deg]
str2.theta = theta(end) * 180/pi; % [deg]
str2.theta_max = theta_max; % [deg]
str2.beta_min = beta_min * 180/pi; % [deg]
str2.beta_max = beta(ind_max) * 180/pi; % [deg]
str2.theta_sonic = theta_sonic; % [deg]
str2.beta_sonic = beta(ind_sonic) * 180/pi; % [deg]

str2.theta_range = theta * 180/pi; % [deg]
str2.un = u2n;
% Plot (pressure, deflection) - shock polar
figure(1); ax = gca;
set(ax, 'LineWidth', config.linewidth, 'FontSize', config.fontsize-2, 'BoxStyle', 'full')
grid(ax, 'off'); box(ax, 'off'); hold(ax, 'on'); ax.Layer = 'Top'; axis(ax, 'tight')

plot(ax, [-flip(theta), theta] * 180/pi, [flip(p2), p2] / pressure(str1), 'k', 'LineWidth', config.linewidth);
% plot(ax, theta_max, p2(ind_max), 'ko', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
% plot(ax, theta(ind_sonic) * 180/pi, p2(ind_sonic), 'ks', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
text(ax, 0, max(p2), strcat('$\mathcal{M}_1 = ', sprintf('%.2g', M1), '$'), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', 'FontSize', config.fontsize-4, 'Interpreter', 'latex')
% text(ax, theta_max, p2(ind_max), 'detachment point \rightarrow m', 'HorizontalAlignment', 'right', 'FontSize', config.fontsize-6)
% text(ax, theta_sonic * 180/pi, p2(ind_sonic), 'sonic point \rightarrow s', 'HorizontalAlignment', 'right', 'FontSize', config.fontsize-6)
xlabel(ax, 'Deflection angle $[^\circ]$', 'FontSize', config.fontsize, 'Interpreter', 'latex');
ylabel(ax, '$p_2/p_1$','FontSize', config.fontsize, 'Interpreter', 'latex');
% ylim(ax, [1, max(p2)]);
% Plot (wave angle, deflection) - shock polar
figure(2); ax = gca;
set(ax, 'LineWidth', config.linewidth, 'FontSize', config.fontsize-2, 'BoxStyle', 'full')
grid(ax, 'off'); box(ax, 'off'); hold(ax, 'on'); ax.Layer = 'Top'; axis(ax, 'tight');

plot(ax, theta * 180/pi, beta * 180/pi, 'k', 'LineWidth', config.linewidth);
% plot(ax, theta_max, beta(ind_max) * 180/pi, 'ko', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
% plot(ax, theta(ind_sonic) * 180/pi, beta(ind_sonic) * 180/pi, 'ks', 'LineWidth', config.linewidth, 'MarkerFaceColor', 'auto');
% text(ax, theta_max, beta(ind_max) * 180/pi, strcat('$\mathcal{M}_1 = ', sprintf('%.2g', M1), '$'), 'VerticalAlignment', 'middle', 'HorizontalAlignment', 'left', 'FontSize', config.fontsize-6, 'Interpreter', 'latex')
xlabel(ax, 'Deflection angle $[^\circ]$','FontSize', config.fontsize, 'Interpreter', 'latex');
ylabel(ax, 'Wave angle $[^\circ]$','FontSize', config.fontsize, 'Interpreter', 'latex');
% xlim(ax, [0, theta_max]);
ylim(ax, [0, 90]);
% Plot velocity components
figure(3); ax = gca;
set(ax, 'LineWidth', config.linewidth, 'FontSize', config.fontsize-2, 'BoxStyle', 'full')
grid(ax, 'off'); box(ax, 'off'); hold(ax, 'on'); ax.Layer = 'Top'; axis(ax, 'tight');

plot(ax, u2x, u2y, 'k', 'LineWidth', config.linewidth);
xlabel(ax, '$u_{2x}$ [m/s]', 'FontSize', config.fontsize, 'Interpreter', 'latex');
ylabel(ax, '$u_{2y}$ [m/s]', 'FontSize', config.fontsize, 'Interpreter', 'latex');
% ylim(ax, [0, max(u2y)]);
mix2.beta = beta(end) * 180/pi; % [deg]
mix2.theta = theta(end) * 180/pi; % [deg]
mix2.theta_max = theta_max; % [deg]
mix2.beta_min = beta_min * 180/pi; % [deg]
mix2.beta_max = beta(ind_max) * 180/pi; % [deg]
mix2.theta_sonic = theta_sonic; % [deg]
mix2.beta_sonic = beta(ind_sonic) * 180/pi; % [deg]
mix2.ind_max = ind_max;
mix2.ind_sonic = ind_sonic;
% Range values for the shock polar
mix2.polar.p = p2; % [bar]
mix2.polar.Mach = M2; % [-]
mix2.polar.un = u2n; % [m/s]
mix2.polar.ux = u2x; % [m/s]
mix2.polar.uy = u2y; % [m/s]
mix2.polar.beta = beta * 180/pi; % [deg]
mix2.polar.theta = theta * 180/pi; % [deg]
mix2.polar.Xi = Xi; % [-]
end

% SUB-PASS FUNCTIONS
Expand Down

0 comments on commit a515dd8

Please sign in to comment.