Area Mach Number Relation

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?

Here my input ARatio is: ARatio = 1.5

and i get two solutions: Msub & Msup

However if i input ARatio as: ARatio = (0.1:0.1:10) i get the following error

Error in AREAMACH2 (line 62)

if (fj*fjp1 > 0)

I need to get Msub & Msup for an array of ARatios from 0.1 to 10, in order to then plot ARatio and the Mach number solutions

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');

[Jim Riggs]

I am curious as to the origin of this function equation.

Can you explain it further?

[Shashwat Shah]

I don't know whether you have got your answer or not but assuming you haven't I'll try to explain. This expression is the relation between Mach number and Area Ratio in a wind tunnel. Area Ratio(AR) is the ratio between cross sectional at any point and area of throat of the nozzle(which is the narrowest part). With this equation, one can determine the Mach variation in WInd tunnel or Rocket nozzle along the change in area. As you can see, there are two Mach numbers soln for every AR. One is Subsonic and one is Supersonic.

Other equation you mentioned of NACA report is just a adiabatic relation in normal shock wave. Actual theories are much more complicated but I tried to explain as easily as possible.

[Jim Riggs]

Thank you for the comment. I spent some more time looking at the NACA 1135 report and I understand.

Answer 1

[Star Strider]

The easiest way is to just enclose the essence of your calculations in a loop:

% Solve for Msub and Msup using this area ratio (A/A*)
ARatio = 0:0.25:2.5;

% 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;

for k = 1:numel(ARatio)
    % 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(k));
            fjp1 = AM_EQN(M+dM,ARatio(k));
            
            % 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(k));
            fjp1 = AM_EQN(M+dM,ARatio(k));
            
            % 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');
end

producing:

==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.5533 after 5 iterations
Msup: 1.5997 after 5 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.4303 after 5 iterations
Msup: 1.8541 after 6 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.3566 after 5 iterations
Msup: 2.0433 after 6 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.3059 after 5 iterations
Msup: 2.1972 after 6 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.2685 after 5 iterations
Msup: 2.3282 after 6 iterations
===================================

==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.2395 after 5 iterations
Msup: 2.4428 after 6 iterations
===================================

It would likely be easier to create a function from your code, then call it with different inputs, producing the desired output.

Make necessary changes to create the result you want.

[Steven Castrillon]

Thank you for your help!

The code you provided succeeded in inputing multiple ARatio values, however although the various iterations are shown, the Msub and Msup values only appear as 1 value each in my workspace. How can the code be modified so that i get an array of values for both Msub and Msup? ie)an Msub and Msup for each ARatio that i input.

[Star Strider]

My pleasure!

Consider:

Msubm(k,:) = [Msub(k), iConvSub(k)];
Msupm(k,:) = [Msup(k), iConvSup(k)];

will collect those values in each matrix to analyse later. Put those assignments in an appropriate place in your code.

With respect to creating a function, it would also be possible to use these as outputs from the function, so it would return them to your workspace.

[Steven Castrillon]

i am not exactly sure where i would place those assignments, i get an error everytime i run.

In the subsonic case, would it be correct to place after

Msub = M; ?

When i do this i get an error

Index exceeds the number of array elements (1).

Error in AMR7 (line 64)

Msubm(k,:) = [Msub(k), iConvSub(k)];

[Star Strider]

I indexed the individual variables in error. (I was thinking one thing — indexing them was one of my solutions — and writing another.)

I put them here, and used this version of the matrices:

    % Set supersonic Mach number to final M from iterations
    Msup = M;
    Msubm(k,:) = [Msub, iConvSub];
    Msupm(k,:) = [Msup, iConvSup];

That worked when I tried it.

It might be more meaningful to include ‘ARatio’:

    % Set supersonic Mach number to final M from iterations
    Msup = M;
    Msubm(k,:) = [ARatio(k), Msub, iConvSub];
    Msupm(k,:) = [ARatio(k), Msup, iConvSup];

Include whatever elements you want to return from the calculations.

[Steven Castrillon]

That did the trick.

Your guidance has been very helpful and greatly appreciated.

Thanks again!

[Star Strider]

As always, my pleasure!

I also was interested in the discussion in: Mach Number Area Relation Several Inputs although aerospace engineering is not an area of my specific (only peripheral) expertise.