Area Mach Number Relation
7 次查看(过去 30 天)
显示 更早的评论
The following code shows me how to get converged Mach number solutions for both subsonic and supersonic given Area ratio (ARatio), however, how can i input a range of ARatio instead of just one value so that the solutions are an array of converged mach numbers?
clear;
clc;
%% INPUTS
% Define some paramters
g = 1.4;
gm1 = g-1;
gp1 = g+1;
% Define anonymous function with two inputs (M and ARatio)
% - Will be used in the methods below
% - Pass M and ARatio as arguments to AM_EQN to get function value
% - funVal = AM_EQN(M,ARatio)
AM_EQN = @(M,ARatio) sqrt((1/M^2)*(((2+gm1*M^2)/gp1)^...
(gp1/gm1)))-ARatio;
% Solve for Msub and Msup using this area ratio (A/A*)
ARatio = 1.5;
% Error tolerance
errTol = 1e-4;
% Flags for printing iterations to screen
verboseBisection = 0;
verboseIncremental = 0;
%% SUBSONIC INCREMENTAL SEARCH
% Initial values
dM = 0.1; % Initial M step size
M = 1e-6; % Initial M value
iConvSub = 0; % Initial converge index
if (verboseIncremental == 1)
fprintf('Incremental Search Method: Subsonic\n');
fprintf('-----------------------------------\n');
end
% Iterate to solve for root
iterMax = 100; % Maximum iterations
stepMax = 100; % Maximum step iterations
for i = 1:1:iterMax
for j = 1:1:stepMax
% Evaluate function at j and j+1
fj = AM_EQN(M,ARatio);
fjp1 = AM_EQN(M+dM,ARatio);
% Print iterations to command window
if (verboseIncremental == 1)
fprintf('fj | fjp1: %3.4f\t%3.4f\n',fj,fjp1);
end
% Update M depending on sign change or not
% - If no sign change, then we are not bounding root yet
% - If sign change, then we are bounding the root, update dM
if (fj*fjp1 > 0)
M = M + dM; % Update M
elseif (fj*fjp1 < 0)
dM = dM*0.1; % Refine the M increment
break; % Break out of j loop
end
end % END: j Loop
% Check for convergence
if (abs(fj-fjp1) <= errTol) % If converged
iConvSub = i; % Set converged index
break; % Exit loop
end
end % END: i Loop
% Set subsonic Mach number to final M from iterations
Msub = M;
%% SUPERSONIC INCREMENTAL SEARCH
% Initial values
dM = 1; % Initial M step size
M = 1+1e-6; % Initial M value
iConvSup = 0; % Initial converge index
if (verboseIncremental == 1)
fprintf('\nIncremental Search Method: Supersonic\n');
fprintf('-------------------------------------\n');
end
% Iterate to solve for root
iterMax = 100; % Maximum iterations
stepMax = 100; % Maximum step iterations
for i = 1:1:iterMax
for j = 1:1:stepMax
% Evaluate function at j and j+1
fj = AM_EQN(M,ARatio);
fjp1 = AM_EQN(M+dM,ARatio);
% Print iterations to command window
if (verboseIncremental == 1)
fprintf('fj | fjp1: %3.4f\t%3.4f\n',fj,fjp1);
end
% Update M depending on sign change or not
% - If no sign change, then we are not bounding root yet
% - If sign change, then we are bounding the root, update dM
if (fj*fjp1 > 0)
M = M + dM; % Update M
elseif (fj*fjp1 < 0)
dM = dM*0.1; % Refine the M increment
break; % Break out of j loop
end
end % END: j Loop
% Check for convergence
if (abs(fj-fjp1) <= errTol) % If converged
iConvSup = i; % Set converged index
break; % Exit loop
end
end % END: i Loop
% Set supersonic Mach number to final M from iterations
Msup = M;
% Print solutions to command window
fprintf('==== INCREMENTAL SEARCH SOLVER ====\n');
fprintf('Msub: %3.4f after %i iterations\n',Msub,iConvSub);
fprintf('Msup: %3.4f after %i iterations\n',Msup,iConvSup);
fprintf('===================================\n\n');
0 个评论
回答(2 个)
Chris
2019-9-30
You can format posts to look llike ML text and make them easier to read. Without reading any of your code it sounds like you can write a wraper function if you can not eaisly vectorize the code.
Chris
2019-10-1
Put all that code info a function with inputs/outputs; call it for each ARatio you want. https://en.wikipedia.org/wiki/Wrapper_function
Control statements with vectors are commonly inapropreate.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Dates and Time 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!