Merge pull request #268 from AlbertoCuadra/develop

Add: compute shock_oblique for a given deflection angle or for a given wave angle #267
CombustionToolbox · Mar 22, 2022 · 2858fa0 · 2858fa0
2 parents 3da0b29 + 038e34a
commit 2858fa0
Show file tree

Hide file tree

Showing 5 changed files with 80 additions and 12 deletions.
diff --git a/Examples/Example_SHOCK_OBLIQUE.m b/Examples/Example_SHOCK_OBLIQUE.m
@@ -2,8 +2,9 @@
 % EXAMPLE: SHOCK_OBLIQUE
 %
 % Compute pre-shock and post-shock state for a oblique incident shock wave
-% at standard conditions, a set of 20 species considered and a initial 
-% shock front velocities u1 = 1500 [m/s]
+% at standard conditions, a set of 20 species considered, a initial 
+% shock front velocities u1 = a1 * 10 [m/s], and a deflection angle 
+% theta = 20 [deg]
 %    
 % Air == {'O2','N2','O','O3','N','NO','NO2','NO3','N2O','N2O3','N2O4',...
 %         'N3','C','CO','CO2','Ar','H2O','H2','H','He'}
@@ -27,7 +28,7 @@
 self.PD.proportion_inerts_O2 = [78.084, 0.9365, 0.0319] ./ 20.9476;
 %% ADDITIONAL INPUTS (DEPENDS OF THE PROBLEM SELECTED)
 overdriven = 10;
-self = set_prop(self, 'u1', 3.472107491008314e+02 * overdriven, 'theta', 51.2161);
+self = set_prop(self, 'u1', 3.472107491008314e+02 * overdriven, 'theta', 20);
 %% SOLVE PROBLEM
 self = SolveProblem(self, 'SHOCK_OBLIQUE');
 %% DISPLAY RESULTS (PLOTS)

diff --git a/Settings/self/ProblemDescription/ProblemDescription.m b/Settings/self/ProblemDescription/ProblemDescription.m
@@ -29,6 +29,10 @@
     self.u1.value = []; % [m/s]
     self.overdriven.description = "Overdriven shock velocity";
     self.overdriven.value = []; % [-]
+    self.theta.description = "Deflection angle - oblique shocks";
+    self.theta.value = [];
+    self.beta.description = "Wave angle - oblique shocks";
+    self.beta.value = [];
     self.S_Fuel = [];
     self.N_Fuel = [];
     self.T_Fuel = [];

diff --git a/Solver/Shocks and detonations/shock_oblique_beta.m b/Solver/Shocks and detonations/shock_oblique_beta.m
@@ -0,0 +1,55 @@
+function [str1, str2] = shock_oblique_beta(varargin)
+    % Solve oblique shock for a given shock wave angle
+
+    % Unpack input data
+    [self, str1, str2] = unpack(varargin);
+    % Definitions
+    a1 = soundspeed(str1);     % [m/s]
+    u1 = str1.u;               % [m/s]
+    M1 = u1 / a1;              % [-]
+    beta = str1.beta * pi/180; % [rad]
+    beta_min = asin(1/M1);     % [rad]
+    beta_max = pi/2;           % [rad]
+    w1 = u1 * sin(beta);       % [m/s]
+    % Check range beta
+    if beta < beta_min || beta > beta_max
+        error('\nERROR! The given wave angle beta = %.2g is not in the range of possible solutions [%.2g, %2.g]', str1.beta, beta_min * 180/pi, beta_max * 180/pi);
+    end
+
+    % Obtain deflection angle, pre-shock state and post-shock states for an oblique shock
+    if ~contains(self.PD.ProblemType, '_R')
+        [~, str2] = shock_incident(self, str1, w1, str2);
+    else
+        [~, ~, str2] = shock_reflected(self, str1, str2.w2, str2);
+    end
+
+    w2 = str2.v_shock;
+    theta = beta - atan(w2 / (u1 .* cos(beta)));
+    u2 = w2 * csc(beta - theta);
+
+    % Save results
+    str2.beta = beta * 180/pi;   % [deg]
+    str2.theta = theta * 180/pi; % [deg]
+    str2.beta_min = beta_min * 180/pi; % [deg]
+    str2.beta_max = beta_max * 180/pi; % [deg]
+    str2.u = u2; % [m/s]
+    str2.w2 = w2; % [m/s]
+    str2.v_shock = u2; % [m/s]
+end
+
+% SUB-PASS FUNCTIONS
+function [self, str1, str2] = unpack(x)
+    % Unpack input data
+    self  = x{1};
+    str1  = x{2};
+    u1    = x{3};
+    beta = x{4};       
+    str1.u = u1;        % velocity preshock [m/s] - laboratory fixed
+    str1.v_shock = u1;  % velocity preshock [m/s] - shock fixed
+    str1.beta = beta;   % wave angle  [deg]
+    if length(x) > 4
+        str2 = x{5};
+    else
+        str2 = [];
+    end
+end
diff --git a/...er/Shocks and detonations/shock_oblique.m → ...cks and detonations/shock_oblique_theta.m b/...er/Shocks and detonations/shock_oblique.m → ...cks and detonations/shock_oblique_theta.m
@@ -1,4 +1,4 @@
-function [str1, str2_1, str2_2] = shock_oblique(varargin)
+function [str1, str2_1, str2_2] = shock_oblique_theta(varargin)
     % Solve oblique shock
 
     % Unpack input data
@@ -48,7 +48,7 @@
         M2 = u2 / soundspeed(str2);
 
         if STOP > TN.tol_oblique || beta < beta_min || beta > pi/2 || theta < 0 || M2 >= M1
-            error('There is not solution for the given deflection angle');
+            error('There is not solution for the given deflection angle %.2g', str1.theta);
         end
 
         % Save results

diff --git a/Solver/SolveProblem.m b/Solver/SolveProblem.m
@@ -70,16 +70,24 @@
         case {'SHOCK_OBLIQUE'}
             try
                 u1 = self.PD.u1.value(i);
-                theta = self.PD.theta.value;
             catch
                 u1 = self.PD.u1.value;
-                theta = self.PD.theta.value;
             end
-%             if i==self.C.l_phi
-                [self.PS.strR{i}, self.PS.str2{i}, self.PS.strP{i}] = shock_oblique(self, self.PS.strR{i}, u1, theta);
-%             else
-%                 [self.PS.strR{i}, self.PS.strP{i}] = shock_oblique(self, self.PS.strR{i}, u1, self.PS.strP{i+1});
-%             end
+            if ~isempty(self.PD.theta.value)
+                try
+                    theta = self.PD.theta.value(i);
+                catch
+                    theta = self.PD.theta.value;
+                end
+                [self.PS.strR{i}, self.PS.str2{i}, self.PS.strP{i}] = shock_oblique_theta(self, self.PS.strR{i}, u1, theta);
+            else
+                try
+                    beta = self.PD.beta.value(i);
+                catch
+                    beta = self.PD.beta.value;
+                end
+                [self.PS.strR{i}, self.PS.strP{i}] = shock_oblique_beta(self, self.PS.strR{i}, u1, beta);
+            end
         case {'SHOCK_OBLIQUE_R'}
             try
                 u1 = self.PD.u1.value(i);