-
Notifications
You must be signed in to change notification settings - Fork 5
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #323 from AlbertoCuadra/develop
Add: compute Regular Reflections for detonations #322
- Loading branch information
Showing
7 changed files
with
249 additions
and
25 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,36 @@ | ||
| % ------------------------------------------------------------------------- | ||
| % EXAMPLE: DET_OBLIQUE_THETA | ||
| % | ||
| % Compute pre-shock and post-shock state for a oblique overdriven detonation | ||
| % considering Chapman-Jouguet (CJ) theory for a stoichiometric CH4-air | ||
| % mixture at standard conditions, a set of 24 species considered, an | ||
| % overdrive of 4 and a set of deflection angles [10:5:40] [deg]. | ||
| % | ||
| % Soot formation == {'CO2', 'CO', 'H2O', 'H2', 'O2', 'N2', 'He', 'Ar',... | ||
| % 'HCN','H','OH','O','CN','NH3','CH4','C2H4','CH3',... | ||
| % 'NO','HCO','NH2','NH','N','CH','Cbgrb'} | ||
| % | ||
| % 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 March 25 2022 | ||
| % ------------------------------------------------------------------------- | ||
|
|
||
| %% INITIALIZE | ||
| self = App('Soot Formation'); | ||
| %% INITIAL CONDITIONS | ||
| self = set_prop(self, 'TR', 300, 'pR', 1 * 1.01325, 'phi', 1); | ||
| self.PD.S_Fuel = {'CH4'}; | ||
| self.PD.S_Oxidizer = {'O2'}; | ||
| self.PD.S_Inert = {'N2', 'Ar', 'CO2'}; | ||
| self.PD.proportion_inerts_O2 = [78.084, 0.9365, 0.0319] ./ 20.9476; | ||
| %% ADDITIONAL INPUTS (DEPENDS OF THE PROBLEM SELECTED) | ||
| overdriven = 4; | ||
| self = set_prop(self, 'overdriven', overdriven, 'theta', 15); | ||
| %% SOLVE PROBLEM | ||
| self = SolveProblem(self, 'DET_OBLIQUE'); | ||
| %% DISPLAY RESULTS (PLOTS) | ||
| postResults(self); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,60 @@ | ||
| function [mix1, mix2] = det_oblique_beta(varargin) | ||
| % Solve oblique shock for a given shock wave angle | ||
|
|
||
| % Unpack input data | ||
| [self, mix1, mix2] = unpack(varargin); | ||
| % Compute CJ state | ||
| [mix1, ~] = cj_detonation(self, mix1); | ||
| mix1.cj_speed = mix1.u; | ||
| % Definitions | ||
| a1 = soundspeed(mix1); % [m/s] | ||
| u1 = mix1.u; % [m/s] | ||
| M1 = u1 / a1; % [-] | ||
| beta = mix1.beta * pi/180; % [rad] | ||
| beta_min = asin(1/M1); % [rad] | ||
| beta_max = pi/2; % [rad] | ||
| overdriven_n = mix1.overdriven * sin(beta); % [-] | ||
| % 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]'], mix1.beta,... | ||
| beta_min * 180/pi, beta_max * 180/pi); | ||
| end | ||
|
|
||
| % Obtain deflection angle, pre-shock state and post-shock states for | ||
| % an oblique detonation | ||
| [~, mix2] = overdriven_detonation(self, mix1, overdriven_n, mix2); | ||
|
|
||
| u2n = mix2.v_shock; | ||
| theta = beta - atan(u2n / (u1 .* cos(beta))); | ||
| u2 = u2n * csc(beta - theta); | ||
| % Save results | ||
| mix2.beta = beta * 180/pi; % [deg] | ||
| mix2.theta = theta * 180/pi; % [deg] | ||
| mix2.beta_min = beta_min * 180/pi; % [deg] | ||
| mix2.beta_max = beta_max * 180/pi; % [deg] | ||
| mix2.u = u2; % [m/s] | ||
| mix2.u2n = u2n; % [m/s] | ||
| mix2.v_shock = u2; % [m/s] | ||
| end | ||
|
|
||
| % SUB-PASS FUNCTIONS | ||
| function [self, mix1, mix2] = unpack(x) | ||
| % Unpack input data | ||
| self = x{1}; | ||
| mix1 = x{2}; | ||
| overdriven = x{3}; | ||
| beta = x{4}; | ||
| mix1.overdriven = overdriven; | ||
| mix1.beta = beta; % wave angle [deg] | ||
| if length(x) > 4 | ||
| mix2 = x{5}; | ||
| else | ||
| mix2 = []; | ||
| end | ||
| if isempty(mix2) | ||
| mix1.cj_speed = []; | ||
| else | ||
| mix1.cj_speed = mix2.cj_speed; | ||
| end | ||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,115 @@ | ||
| function [mix1, mix2_1, mix2_2] = det_oblique_theta(varargin) | ||
| % Solve oblique shock | ||
|
|
||
| % Unpack input data | ||
| [self, mix1, mix2_1, mix2_2] = unpack(varargin); | ||
| % Abbreviations | ||
| TN = self.TN; | ||
| % Compute CJ state | ||
| [mix1, ~] = cj_detonation(self, mix1); | ||
| mix1.cj_speed = mix1.u; | ||
| % Definitions | ||
| a1 = soundspeed(mix1); % [m/s] | ||
| u1 = mix1.u; % [m/s] | ||
| M1 = u1 / a1; % [-] | ||
| overdriven = mix1.overdriven; % [-] | ||
| theta = mix1.theta * pi/180; % [rad] | ||
| beta_min = asin(1/M1); % [rad] | ||
| beta_max = pi/2; % [rad] | ||
| lambda = 0.7; | ||
| m = 4; | ||
| % Solve first branch (weak shock) | ||
| beta_guess = 0.5 * (beta_min + beta_max); % [rad] | ||
| mix2_1 = solve_det_oblique(mix2_1); | ||
| % Solve second branch (strong shock) | ||
| beta_guess = beta_max * 0.97; % [rad] | ||
| mix2_2 = solve_det_oblique(mix2_2); | ||
| % NESTED FUNCTIONS | ||
| function mix2 = solve_det_oblique(mix2) | ||
| % Obtain wave angle, pre-shock state and post-shock states for an oblique shock | ||
| STOP = 1; it = 0; itMax = TN.it_oblique; | ||
| while STOP > TN.tol_oblique && it < itMax | ||
| it = it + 1; | ||
| % Compute incident normal velocity | ||
| overdriven_n = overdriven * sin(beta_guess); % [-] | ||
| % Obtain post-shock state | ||
| [~, mix2] = overdriven_detonation(self, mix1, overdriven_n, mix2); | ||
| u2n = mix2.v_shock; | ||
| % Compute f0 , df0, and d2f0 | ||
| f0 = f0_beta(beta_guess, theta, u2n, u1); | ||
| df0 = df0_beta(beta_guess, u2n, u1); | ||
|
|
||
| % Apply correction | ||
| % beta = beta_guess - f0 / df0; | ||
|
|
||
| d2f0 = d2f0_beta(beta_guess, u2n, u1); | ||
| beta = beta_guess - (f0 * df0) / (df0^2 - m^-1 * f0 * d2f0); | ||
| % Check value | ||
| if beta < beta_min || beta > beta_max | ||
| % beta = beta_guess - lambda * f0 / df0; | ||
| beta = beta_guess - lambda * (f0 * df0) / (df0^2 - 0.5 * f0 * d2f0); | ||
| end | ||
| % Compute error | ||
| STOP = max(abs(beta - beta_guess) / abs(beta), abs(f0)); | ||
| % Update guess | ||
| beta_guess = beta; | ||
| end | ||
|
|
||
| u2 = u2n * csc(beta - theta); | ||
| M2 = u2 / soundspeed(mix2); | ||
|
|
||
| if beta < beta_min || beta > pi/2 || theta < 0 || M2 >= M1 | ||
| error('There is not solution for the given deflection angle %.2g', mix1.theta); | ||
| end | ||
| if STOP > TN.tol_oblique | ||
| print_error_root(it, itMax, mix2.T, STOP); | ||
| end | ||
|
|
||
| % Save results | ||
| mix2.beta = beta * 180/pi; % [deg] | ||
| mix2.theta = theta * 180/pi; % [deg] | ||
| mix2.beta_min = beta_min * 180/pi; % [deg] | ||
| mix2.beta_max = beta_max * 180/pi; % [deg] | ||
| mix2.u = u2; % [m/s] | ||
| mix2.un = u2n; % [m/s] | ||
| mix2.v_shock = u2; % [m/s] | ||
| end | ||
| end | ||
|
|
||
| % SUB-PASS FUNCTIONS | ||
| function [self, mix1, mix2_1, mix2_2] = unpack(x) | ||
| % Unpack input data | ||
| self = x{1}; | ||
| mix1 = x{2}; | ||
| overdriven = x{3}; | ||
| theta = x{4}; | ||
| mix1.overdriven = overdriven; | ||
| mix1.theta = theta; % deflection angle [deg] | ||
| if length(x) > 4 | ||
| mix2_1 = x{5}; | ||
| mix2_2 = x{6}; | ||
| else | ||
| mix2_1 = []; | ||
| mix2_2 = []; | ||
| end | ||
| if isempty(mix2_2) | ||
| mix1.cj_speed = []; | ||
| else | ||
| mix1.cj_speed = mix2_1.cj_speed; | ||
| end | ||
| end | ||
|
|
||
| function value = f0_beta(beta, theta, u2n, u1) | ||
| % Function to find the roots of beta | ||
| value = theta + atan(u2n / (u1 .* cos(beta))) - beta; | ||
| end | ||
|
|
||
| function value = df0_beta(beta, u2n, u1) | ||
| % Derivative of the function to find the roots of beta | ||
| value = (u2n * u1 * sin(beta)) / (u2n^2 + u1^2 * cos(beta)^2) - 1; | ||
| end | ||
|
|
||
| function value = d2f0_beta(beta, u2n, u1) | ||
| % Derivative of the function to find the roots of beta | ||
| value = 0.5 * (u1 * u2n * cos(beta) * (3*u1^2 + 2*u2n^2 - u1^2 * cos(2*beta))) / (u2n^2 + u1^2 * cos(beta)^2)^2; | ||
| end |