Merge pull request #310 from AlbertoCuadra/develop

Add: solve Regular Reflection - shock polar for a given theta or beta #309
CombustionToolbox · Mar 24, 2022 · b6ae2e4 · b6ae2e4
2 parents 4407591 + 7d866db
commit b6ae2e4
Show file tree

Hide file tree

Showing 6 changed files with 42 additions and 111 deletions.
diff --git a/Settings/functions/Display/plot_shock_polar.m b/Settings/functions/Display/plot_shock_polar.m
@@ -1,17 +1,27 @@
-function [ax1, ax2, ax3] = plot_shock_polar(self, mix1, mix2)
+function [ax1, ax2, ax3] = plot_shock_polar(varargin)
     % Routine to obtain shock polar plots
     %   * Plot (pressure, deflection) 
     %   * Plot (wave angle, deflection)
     %   * Plot velocity components
-
+
+    self = varargin{1};
+    mix1 = varargin{2};
+    mix2 = varargin{3};
+    if nargin > 3
+        mix2_case = varargin{4};
+        mix0 = varargin{5};
+    else
+        mix2_case = cell(1, length(mix1));
+        mix0 = cell(1, length(mix1));
+    end
     % 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);
+	    ax1 = plot_shock_polar_pressure(mix1{i}, mix2{i}, config, mix2_case{i}, mix0{i});
         % Plot (wave angle, deflection) - shock polar
         ax2 = plot_shock_polar_wave(mix1{i}, mix2{i}, config);
         % Plot velocity components
@@ -29,19 +39,30 @@
     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)
+function ax = plot_shock_polar_pressure(mix1, mix2, config, mix2_case, mix0)
     % 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);
+    p1 = pressure(mix1);
+    Mach_label = '$\mathcal{M}_1 = ';
+    xaxis = 0;
+    if ~isempty(mix2_case)
+        M1 = velocity_relative(mix2_case) / soundspeed(mix2_case);
+        p1 = mix0.p;
+        Mach_label = '$\mathcal{M}_2 = ';
+        xaxis = mix2_case.theta;
+        y_lim = get(ax, 'YLim'); y_lim(2) = max(mix2.polar.p) / p1;
+        plot(ax, [0, 0], y_lim, '-.k', 'LineWidth', config.linewidth)
+    end
+    plot(ax, xaxis + [- flip(mix2.polar.theta), mix2.polar.theta], [flip(mix2.polar.p), mix2.polar.p] / p1, '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, xaxis, max(mix2.polar.p), strcat(Mach_label, 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');
+    ylabel(ax, '$p/p_1$','FontSize', config.fontsize, 'Interpreter', 'latex');
 %     ylim(ax, [1, max(str2.polar.p2)]);
 end
 

diff --git a/Settings/functions/check_inputs.m b/Settings/functions/check_inputs.m
@@ -27,7 +27,7 @@
                 self = set_prop(self, 'pP', self.PD.pR.value); % Guess
             case {'SHOCK_I', 'SHOCK_R'} % * SHOCK_I and SHOCK_R: CALCULATE PLANAR SHOCK WAVE
                 self = check_inputs_prop(self, 'u1');
-            case 'SHOCK_POLAR' % * SHOCK_OBLIQUE: CALCULATE OBLIQUE SHOCK WAVE
+            case {'SHOCK_POLAR', 'SHOCK_POLAR_R'} % * SHOCK_POLAR: CALCULATE OBLIQUE SHOCK POLARS
                 self = check_inputs_prop(self, 'u1');
             case 'SHOCK_OBLIQUE' % * SHOCK_OBLIQUE: CALCULATE OBLIQUE SHOCK WAVE
                 self = check_inputs_prop(self, 'u1');

diff --git a/Settings/functions/postResults.m b/Settings/functions/postResults.m
@@ -78,6 +78,13 @@ function postResults(self)
 %         beta = cell2vector(mix2{i}.polar, 'beta');
 %         displaysweepresults(self, mix2{i}.polar, beta);
 %     end
+elseif strcmp(ProblemType,{'SHOCK_POLAR_R'})
+    mix3 = mix2;
+    mix2 = self.PS.str2;
+    % Shock polars incident
+    plot_shock_polar(self, mix1, mix2);
+    % Shock polars reflected
+    plot_shock_polar(self, mix2, mix3, self.PS.str2_case, mix1);
 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]';

diff --git a/Solver/Shocks and detonations/shock_polar.m b/Solver/Shocks and detonations/shock_polar.m
@@ -53,6 +53,7 @@
     % Range values for the shock polar
     mix2.polar.p = p2;    % [bar]
     mix2.polar.Mach = M2; % [-]
+    mix2.polar.u = u2;    % [m/s]
     mix2.polar.un = u2n;  % [m/s]
     mix2.polar.ux = u2x;  % [m/s]
     mix2.polar.uy = u2y;  % [m/s]

diff --git a/Solver/Shocks and detonations/shock_polar_reflected.m b/Solver/Shocks and detonations/shock_polar_reflected.m
diff --git a/Solver/SolveProblem.m b/Solver/SolveProblem.m
@@ -101,8 +101,9 @@
                     theta = self.PD.theta.value;
                 end
                 self.PD.ProblemType = 'SHOCK_OBLIQUE';
-                [self.PS.strR{i}, self.PS.str2{i}, self.PS.str2_2{i}] = shock_oblique_theta(self, self.PS.strR{i}, u1, theta);
-
+                [self.PS.strR{i}, self.PS.str2{i}, ~] = shock_oblique_theta(self, self.PS.strR{i}, u1, theta);
+                self.PD.ProblemType = 'SHOCK_OBLIQUE_R';
+                [self.PS.strR{i}, self.PS.str2{i}, self.PS.str3{i}, self.PS.str3_2{i}] = shock_oblique_reflected_theta(self, self.PS.strR{i}, self.PS.str2{i}.u, self.PS.str2{i}.theta, self.PS.str2{i});
             else
                 try
                     beta = self.PD.beta.value(i);
@@ -128,7 +129,8 @@
                 u1 = self.PD.u1.value;
             end
                 [self.PS.strR{i}, self.PS.str2{i}] = shock_polar(self, self.PS.strR{i}, u1);
-                [self.PS.strR{i}, self.PS.str2{i}, self.PS.strP{i}] = shock_polar_reflected(self, self.PS.strR{i}, u1, self.PS.str2{i});
+                [~, self.PS.str2_case{i}, ~] = shock_oblique_theta(self, self.PS.strR{i}, u1, self.PD.theta.value);
+                [~, self.PS.strP{i}] = shock_polar(self, self.PS.str2_case{i}, self.PS.str2_case{i}.u);
         case {'DET'}
             if i==self.C.l_phi
                 [self.PS.strR{i}, self.PS.strP{i}] = cj_detonation(self, self.PS.strR{i});